{'fraction_non_alphanumeric': 0.04734258721697367, 'fraction_numerical': 0.014751675727028025, 'mean_word_length': 3.778688524590164, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Convergence of a crystalline algorithm for the heat equation in one dimension and for the motion of a graph by weighted curvature

Motion by (weighted) mean curvature is a geometric evolution law for surfaces, representing steepest descent with respect to (an)isotropic surface energy. It has been proposed that this motion could be computed by solving the analogous evolution law using a"crystalline"approximation to the surface energy. We present the first convergence analysis for a numerical scheme of this type. Our treatment is restricted to one dimensional surfaces (curves in the plane) which are graphs. In this context, the scheme amounts to a new algorithm for solving quasilinear parabolic equations in one space dimension.


Introduction
In the modeling of phase transformations it is often of interest to consider surfaceenergy-driven motion of interfaces (see the work of Gurtin [14] [15] [16], the recent review by Taylor, Cahn, and Handwerker [27], and the references therein). If the surface energy is isotropic this leads to motion by mean curvature, i.e. the normal velocity of the interface equals its mean curvature. If the surface energy is anisotropic the associated evolution equation has been called "motion by weighted mean curvature" (see the review Taylor [25]).
When the surface energy is "strictly convex" the evolution law is a quasilinear parabolic equation. The isotropic case falls in this class. When the surface energy is "crystalline" the surface must be faceted and its evolution law reduces to a family of ordinary differential equations for the lengths of the faces. The mathematical theory of surface-energy-driven motion of interfaces is discussed at length in Angenent and Gurtin [2], developing upon an extensive materials science literature of which Herring [17] is representative. The essential aspects are summarized in an appendix, for the reader's convenience. We emphasize that familiarity with this theory is not assumed in the present work, though it provides the motivation and the context for what we do.
There has recently been intense mathematical activity concerning the analysis of motion by (weighted) mean curvature (see for example Angenent and Gurtin [2] [3], Brakke [4], Chen, Giga, and Goto [5], de Mottoni and Schatzman [6], Evans and Spruck [7] [8], Gage and Hamilton [11], Grayson [13], Gurtin [16], Huisken [18], and Sethian [22]). The analysis of crystalline surface energies has also received considerable attention both with regard to statics (energy minimization) (see Sullivan [23] and Taylor [24]) and with regard to dynamics (motion by weighted mean curvature) (see Almgren, Taylor and Wang [1], Angenent and Gurtin [2], Ohnuma and Sato [20], and Roosen and Taylor [21]). We provide these references for the interested reader; the only one we actually need in this paper is [2]. It is a geometrically natural idea to approximate a strictly convex energy by a crystalline one. This idea has recently been analyzed by Sullivan in [23], but only for problems of energy minimization. "Crystalline approximation" has been used in applications without rigorous justification as a method for computing motion by weighted mean curvature in the context of closed curves and surfaces, because it replaces a parabolic differential equation with a system of ordinary differential equations. This paper represents a first attempt to establish its convergence. We study the simplest possible case: the interface is the graph of a function of one variable. Our main result is the convergence of a "crystalline" scheme for solving a quasilinear parabolic equation in one space dimension. We prove convergence in H 1 , with a specified rate. Our method is somewhat similar to the convergence analysis for a Galerkin approximation (see e.g. Thomée [28]).
We emphasize that the "crystalline" discretization considered here is quite different from any of the more standard schemes from numerical analysis. The approximate solution is a piecewise linear function of the spatial variable, with "pieces" that have fixed slopes and variable lengths. This is in sharp contrast to finite difference and finite element methods (where it is the spatial grid that is fixed), to front tracking (which places no constraint on the slopes of the approximate solution), and to spectral methods.
The recent paper of Fukui and Giga [10] is closely related to the present work. They prove a general existence and uniqueness theorem for motion of graphs by weighted curvature (by adapting the theory of nonlinear semigroups), which applies to a wider class of interfacial energies and even to "incompatible" initial data. When the surface energy is "crystalline" and the initial data is piecewise linear and compatible (i.e. has the "right" slopes) their solution is the same as ours. They prove continuous dependence on the initial data and on the form of the interfacial energy. This fact would suffice to prove convergence of our approximation scheme. The result we prove here is sharper, however, because we get a specific convergence rate in the H 1 norm.
In another paper [12] we prove convergence of the crystalline algorithm for the motion of a simple closed convex curve by weighted curvature. This is done by parametrizing the weighted curvature by the angle between the interface normal and a fixed coordinate axis, and comparing it with the weighted curvature of the approximate solution. It turns out that the crystalline approximation scheme corresponds to a standard finite difference scheme for the (nonlinear) evolution equation for the weighted curvature.
The organization of this paper is as follows: in Section 2 we set up the scheme and make some preliminary observations. In Section 3 we study the case of the ordinary heat equation. In Section 4 we prove convergence for general convex energies and constant Dirichlet boundary conditions. In Section 5 we prove convergence for constant Neumann boundary conditions. In Section 6 we set up the general Dirichlet problem. However, for general Dirichlet boundary conditions we do not prove convergence of the scheme because we have been unable to prove that the approximate solution exists up to a fixed time as the discretization gets finer. It seems plausible that, under suitable growth conditions on the energy, one should be able to bound the slope of the approximate solutions. Then convergence would follow as in Section 4. Finally, Section 7 is an appendix which outlines the physical and mathematical context of our work.
This paper is part of the first author's Ph.D. thesis.

Setup and preliminary observations
This paper is concerned with the convergence of an approximation scheme for the equation . (1.1) For now we focus on the case of homogeneous Dirichlet boundary conditions, postponing consideration of other boundary conditions to Sections 5 and 6. If u xx is the curvature of the graph of u, and in general it is the negative of the gradient of or the weighted curvature of the graph of u. So Eq. (1) says that the normal velocity of the graph of u equals its weighted curvature.
We assume throughout that W is strictly convex and C 3 , and that Eq. (1) has a C 3 solution in [0, 1] × [0, T ]. The existence of a C 3 solution is guaranteed by Theorem VI.5.2 and the regularity results of Ladyženskaja, Solonnikov, and Ural'ceva [19] when W and u 0 are sufficiently smooth, u 0 satisfies compatibility conditions at zero and one, and W satisfies certain growth conditions. Appropriate growth conditions are discussed in Section 6 (see especially conditions (32)).
The numerical scheme considered in this paper arises by deriving an analogue of Eq. (1) in the setting where W is substituted by a piecewise linear function, W , which coincides with W at its corners (see Figure 1). We now carry out this derivation. As explained in the Appendix, this has been done in a slightly different way by Gibbs (see the review Taylor [25]). The equations we will arrive at have been obtained in a physical context by Angenent and Gurtin [2]. The approximation scheme should deal with continuous piecewise linear functions, u m , where each segment, which from now on will be called face, has an admissible slope, i.e. has slope equal to a value of u x corresponding to one of the corners ofW . We will refer to the set of all such functions as the set of admissible functions. For now, the parameter m has no meaning by itself, but later on it will be the maximum distance between two such adjacent values of u x , and it will determine the rate of convergence of the scheme, thus playing a role similar to the one that the mesh size does in a finite difference scheme.
We order the faces from left to right. At selected times faces will disappear and each time this happens we reorder them. So the number of faces is a piecewise constant function of time denoted simply by N . We usually keep the time dependence of such functions implicit. The exception will be u's time dependence and that of its first argument, when that argument is itself a function, which we often make explicit. Let be the abscissa of the ith corner of u m , l i be the length of the projection of the ith face on the x-axis, and (u m x ) i and (u m t ) i be the values of u m x and u m t , respectively, for the ith face of u m .
First of all, we focus on the initial data. Define u m 0 (x) △ = u m (x, 0). We assume that u m 0 belongs to the set of admissible functions, satisfies the boundary conditions, and jumps in ((u m 0 ) x ) i correspond to adjacent corners ofW (see Figure 2).
This is a reasonable assumption, and we make it because our scheme will move but not create faces. Figure 2. A function u m 0 ( · ). Second, we impose the boundary conditions. Since u(0, t) = u(1, t) = 0 we have to choose (u m t ) 1 = 0 and (u m t ) N = 0 for them to be satisfied. Our next goal is to derive the evolution equation of u m . Before doing this, however, we use the fact that u m is continuous to derive some important relations. The functions u m x and u m t are discontinuous, but the vertical velocity of the ith corner of u m is the same whether computed using the ith or (i + 1)th face. So, differentiating u m (x i (t), t) with respect to t, It is sometimes more convenient to work with the lengths of the faces rather than with the positions of the corners of u m . Eqs. (4) implẏ but these can be simplified using the boundary conditions. We geṫ and also (see Eqs. (10.18) of [2]) Finally, we turn to determining the evolution equation of u m . We compute the normal velocities of the faces of u m in terms of the l i 's. They should equal the negative of the gradient ofĒ on the space of admissible functions, at u m . By the definition ofW , We now compute the negative of the gradient ofĒ. Eqs. (4) imply that the rate of change of energy "at" the ith corner is The boundary terms cancel because of the boundary conditions. The preceding expression is the directional derivative ofĒ, in the direction of the piecewise constant function which has value (u m t ) i for x in (x i−1 , x i ), at u m . In other words, it is equal to where φ h is admissible, φ 0 = u m , and Our space of admissible functions is somewhat unfamiliar, so perhaps some elaboration is appropriate. Suppose φ h is a curve in the space of admissible functions such that φ 0 = u m . Since φ h is admissible and φ 0 = u m , we see that and for 2 ≤ i ≤ N − 1, and has value zero for x in (x 0 , x 1 ) and (x N −1 , x N ). We conclude that, for 2 ≤ i ≤ N − 1, the ith face should have In summary, u m is determined by solving Eqs. (5) with the c j i 's given by Eqs. (6), the (u m t ) i 's by Eqs. (3) and (8), and the ∆ i 's by Eqs. (7); the initial data should satisfy condition (2). This amounts to solving a nonlinear system of ordinary differential equations for the l i 's.
We make here a few comments.
• For reference we note the following elementary "summation by parts" formulas: 1 + ( · ) 2 then ∆i li may be interpreted as the "curvature of u m at its ith segment." We now make some preliminary observations about this scheme.
First, note that one can think of u m as being defined on the whole real line, the extension being odd with period two. Then (u m This is important because it shows that for certain qualitative features (such as our next observation) there is no need to consider the end intervals separately. It will also be useful in our computations when we want to extend the range of some summations by adding and subtracting boundary terms. For example, in the expression for the total rate of change of energy we can let the sum run from 1 to N . Then the new boundary term at Remark: Eqs. (7) and (8) now hold for i equal to 1 and N .
Second, we show that even as time evolves and certain faces disappear, jumps in (u m x ) i correspond to adjacent corners ofW , and the maximum principle holds for u m t , i.e. u m t is bounded above by its maximum at time zero. Statement (9) has also been noted by Taylor [26]. Notice the simple but important consequence of statement (9) that if (u m x ) i is smaller or greater than both (u m Preparing for the proof, we note that if statement (9) is true at time t then, at that time, due to the convexity of W . We proceed in four steps: • Statement (9) holds at time zero, by condition (2). It clearly continues to hold until the time t 1 when one or more faces disappear. • The maximum principle holds for u m t for t < t 1 . This is a consequence of the fact that if max 1≤i≤N (u m t ) i is positive then it is non-increasing. For the proof, suppose the maximum occurs at the ith face.
; the first of these inequalities gives c −1 i < 0 and the second c 1 i < 0. It follows that the right hand side of Eq.
At time t 1 the faces i − 1 and i + 1 join to from a single face. (Note that in this case, for , or else the opposite of both of these inequalities.) -Case (ii) The faces i − 1 and i disappear without their adjacent neighbors doing so.
• If the maximum of u m t is discontinuous at time t 1 then it can only jump down. In fact, consider what happens in each of the cases above: -Case (i) If the (i − 1)th and (i + 1)th faces join to form an inflection face then for this new face u m t = 0. Otherwise the new face has a value of ∆ different from zero. If ∆ i−1 ∆ i+1 < 0 for t < t 1 then the new face has a ∆ equal to ∆ i−1 or ∆ i+1 (and the length of the new face at time t 1 is greater than the limit of the lengths of both l i−1 and l i+1 as t approaches t 1 from below). If ∆ i−1 ∆ i+1 = 0 for t < t 1 then the new face might have a ∆ different from ∆ i−1 and ∆ i+1 . But then its ∆ must be equal to the value of ∆ j for some j such that the jth face joins the (i − 1)th and (i + 1)th ones at time t 1 to form the new face. (In this case the length of the new face at time t 1 is greater than the limit of the length of l j as t approaches t 1 .) -Case (ii) The faces i − 2 and i + 1 do not join to form a new face. Therefore at time t 1 the maximum of u m t can only decrease.
By applying the above reasoning repeatedly, we conclude that statement (9) and the maximum principle for u m t hold for all time. The above shows in particular that u m exists for all time.
Third, we note the discrete analogue of d dt . This is an identity valid for any function u. The discrete version applies to any admissible function u m , i.e. one which is piecewise linear in space with slopes (u m x ) i , such that the jumps in (u m x ) i correspond to adjacent corners ofW . The proof uses only Eqs. (5.1 ′ ) and (5.2), not the discretized differential equation (8) or the Dirichlet boundary condition (5.1): For the solution of Eq. (1) we have d Similarly, for the solution of the discretized problem with homogeneous Dirichlet boundary condition we see from Eq.
Fourth,Ē(u m )( · ) is a continuous function of t but is only piecewise C 1 because faces disappear at selected times. We will not mention this restriction again, although it will come up repeatedly, since it is of no consequence for our analysis.
Fifth and last, the maximum principle holds for u m because an interior maximum will occur at the ith corner ( and in this case both (u m t ) i and (u m t ) i+1 are nonpositive.

Convergence for the heat equation
Before proving convergence of the approximation scheme for Eq. (1), let us treat the simpler case of the corresponding approximation scheme for the heat equation, In this case We shall estimate the growth of the H 1 norm of u − u m . The calculation in this section would be simplest if we took full advantage of the boundary conditions, dropping terms such as u t (0, t) or (u m t ) 1 which vanish due to Eqs. (1.2) and (3). This would be inefficient because it would force us to repeat many calculations in Sections 5 and 6, where we consider other boundary conditions. We therefore proceed as follows: • The Dirichlet boundary conditions u(1, t) = 0 and (u m t ) N = 0 will be used at x = 1, and we shall drop terms which vanish as a result.
• The Dirichlet boundary conditions will also be used at x = 0, but we shall label ((i) through (vii)) the associated boundary terms the first time they appear for future reference. Then these terms will be dropped since they vanish due to u(0, t) = 0 and (u m t ) 1 = 0. • In Sections 5 and 6 we will deal with boundary conditions other than the homogeneous one u(0, t) = (u m t ) 1 = 0. Specialized to the case of the heat equation, they lead to relations of the form where r 1 is a function of time, 0 < r 1 ≤ 1, and ∆ 1 is defined by Eq. (7) for corresponding to the odd, periodic extension of u m . This yields ∆ 1 = 0 and (u m t ) 1 = 0, as noted above. In this case we may as well set r 1 = 1 since both sides of Eq. (14) vanish.
One further remark. The right hand side of Eq. (11) can be written For the case of the heat equation, discretized as above, this becomes Now we begin the convergence analysis. We have d dt We have dropped the terms which vanish due to the homogeneous boundary condition, as we announced we would do. We rewrite A and B: Note that we have used Eqs. (4). So, On the other hand, .
Hence, A + B + C equals for some α i and β i belonging to (x i−1 , x i ). The last term in this sum can be written as for some γ 1 belonging to (0, x 1 ). For zero boundary conditions each of (i) through (vii) vanishes. We stop to note that our calculation was exact up to here. We will now bound A + B + C. The final bound for the L 2 norm of u x − u m x will be sharp inasmuch as the next bound is optimum. The terms corresponding to α 1 and γ 1 combine and we obtain here || || ∞ is the supremum norm in x on the interval [0, 1]. Thus the accuracy of the discretization is governed by The parameter m controls the convergence of the scheme in much the way the mesh size controls the behavior of a finite difference approximation. Notice that our estimates do not require any uniformity on the quantities | (u m To recapitulate, we have shown that if u solves the heat equation (12) and u m solves its discretized version, both with homogeneous Dirichlet boundary conditions, then d dt Suppose that the discrete initial data u m 0 have been chosen so that Then the error at time t is By Poincaré's inequality the L 2 norm of u − u m is controlled by the L 2 norm of its spatial derivative. We conclude that with c depending on u and T but not on m.
Since we are assuming u 0 is C 3 , one can choose u m 0 such that the hypothesis (19) is satisfied. For example, one can proceed as follows. First, identify all the points where (u 0 ) x is admissible and draw the tangents to u 0 through those points. Next, draw a line through each inflection pointx of u m with the admissible slope (u m x ) i closest to the value of (u 0 ) x (x) and such that | (u m Then, draw a line through (0, 0) with the admissible slope closest to the value of (u 0 ) x (0) and greater than (respectively, less than) Remark. We collect the boundary terms which entered in the calculation above because we will need them later. They were: If (u m t ) 1 = u t (0, t), as in Section 6, then the sum (i) through (vii) is zero. Alternatively, if (u m x ) 1 = u x (0, t), as in Section 5, then Before proceeding to the next section we make a final remark. If the distance between any two adjacent admissible slopes is a fixed constant, i.e. the mesh in the u x -axis of Figure 1 is uniform, then it might be possible to improve the bound for A + B + C. In fact, if the ith is not an inflection face of u m then either (u m x ) i−1 < (u m x ) i < (u m x ) i+1 or both these inequalities are reversed. In either case the absolute value of the ith term in the expression for the sum of A + B + C is equal to for some µ i belonging to (x i−1 , x i ), so less than If the ith is an inflection face of u m then the ith term of A + B + C is bounded by This suggests that to get a sharper estimate than the one obtained above one should bound max i l i on an interval [0, T ] in terms of m. Note that for the present case Eqs. (5.2) reduce to between inflection faces.

Convergence for general convex energies
We turn now to proving convergence of the approximation scheme described in Section 2, for Eq. (1). The argument is in many respects similar to that for the heat equation. Of course, the functions u and u m are defined differently than they were in Section 3, because we are solving different differential equations. However, we shall keep the notation and labels used in the previous section to highlight the parallels between the two calculations. We shall also keep the conventions concerning the handling of boundary terms.
The main idea is, once again, to control the evolution of the H 1 norm of u − u m . Arguing as before, we get d dt We have used expression (15) and Eqs. (16) and (17). It is still convenient to rewrite B: It is easy to check that expression (21) can be written as the sum to be zero. For such an i the ith term in I + II is The term I is analogous to the term − 1 0 |u t − u m t | 2 dx in Section 3. It is negative. The analogue of II is identically zero for the heat equation. So II is an error term due to the nonlinearity of the equation. The terms C and D will combine much as A + B + C did above.
To estimate C + D, it is convenient to introduceW , Notice that Eq. (1.1) can be written in divergence form We proceed essentially as we did for the heat equation. On the one hand, where of course .
On the other hand, D equals Therefore C + D equals for some ξ i and ζ i belonging to (x i−1 , x i ). But, for some ϕ 1 belonging to (0, x 1 ). For zero boundary conditions each of (i) through (vii) vanishes. The terms corresponding to ξ 1 and ϕ 1 combine and we obtain It remains to estimate II, which we write with δ an arbitrary positive number. We shall control (II 1 ) i later, by showing that it is dominated by the corresponding term of I, if δ is sufficiently small. In order to estimate the term (II 2 ) i , we rewrite the two terms on the right hand side of Eqs. (7) using Taylor's expansion, where θ i , ϑ i belong to the interval with endpoints (u m x ) i and (u m x ) i+1 . We have not yet used the discretized differential equation; recall that (u m where We use here the fact that | (u m x ) i and u x (x, t), and thus These terms are clearly controllable. We reorganize, I + II + C + D = ( I + II 1 ) + ( II 2 − E ) + C + D + E (there is no risk of confusing this E with the energy). We have shown that d dt We can now prove the following Suppose that the initial data for u m satisfies and satisfies a differential inequality of the form when m is sufficiently small. The coefficients Λ and Γ are independent of m, though they depend upon the underlying solution u. In particular, for any T > 0, Proof. Before we start, note that it is always possible to choose initial data satisfying conditions (28) and (29). This was discussed in Section 3.
Recall from Section 2 that the discrete evolution never introduces new faces after the initial time. This fact together with hypothesis (28) implies that max i {(u m x ) i } is bounded with the bound uniform in m. Furthermore, because of (9) on page 6 and hypothesis (29), for all time, not just for time zero. Almost all the work for the theorem has already been done. The principal task that remains is to choose the parameter δ in such a way that Ξ i ≥ 0 whenever (u m x ) i+1 = (u m x ) i−1 . From Eq. (26) and inequality (27) we see that there exist constants c 1 and c 2 such that with both these quantities strictly positive. Their values depend only on the maximum and minimum slopes of the initial data, and on the properties of the "surface energy" function W . We restrict m to m ≤ c 1 /(2c 2 ) and take δ = c 1 . Then Ξ i ≥ 0 and III is nonpositive.
We may now estimate C, D, E, and IV . By inequality (25), by inequality (27), and from the expression for IV i , The maximum in the last formula is not taken over all values of η, but rather over the range of possible values of u x and u m x ; these values lie in a bounded interval which depends only on the initial data for u. Clearly, with the function Γ independent of m. The estimate for IV can be expressed as This yields inequality (30), and a standard application of Gronwall's inequality gives inequality (31).
Eq. (1) has the form with ψ(ν) = (1 + ν 2 ) 1 2 . We have made no particular use of the form of ψ. A similar convergence theorem holds for Eq. (32) with any ψ, continuously differentiable, positive, and bounded away from zero.
One final remark. For the case of the heat equation it was clear that m = max i | (u m x ) i − (u m x ) i−1 | was the parameter controlling the accuracy of the approximation. In the general case one might have expected something different, for example that one should take a coarser discretization where W ′′ is smaller. Our estimates do not support such an idea. In fact, the argument at the end of Section 3 suggests that if one chooses the admissible slopes so that the distance between any two adjacent ones is a fixed constant, then one might be able to prove quadratic convergence in m.

Convergence for the Neumann problem
Here we study the Neumann problem with a and b constants. We consider only discretizations for which a and b are admissible slopes and we choose (u m x ) 1 = a and (u m x ) N = b for the boundary conditions to be satisfied.
We want to insure that faces 1 and N do not disappear and that statement (9) holds. So at each time we extend u m in such a way that (u m x ) 0 and (u m x ) 1 correspond to adjacent corners ofW , (u m x ) 1 lies between (u m x ) 0 and (u m x ) 2 , and the length of the first face of the extension of u m is l 1 / r 1 where (note that 0 < r 1 < 1 and if the value of (u m x ) 2 changes then r 1 might jump). The reason for this choice of r 1 will become clear below, when we prove convergence of the approximation scheme. Similar conditions apply to the extension at the right endpoint with Furthermore, we impose that for the extension of u m Eqs. (u m it is natural to do so since we want our scheme to amount to motion by weighted curvature. The idea is that the law governing the evolution of u m should be the same in both the interior and boundary of the interval [0, 1] (recall also, from the Remark on page 6, that for homogeneous Dirichlet boundary conditions one could think of u m as being defined on the whole real line, the extension being odd with period two; in that case Eqs. (7) and (8) (8), the ∆ i 's by Eqs. (7), and r 1 and r N by Eqs. (34); the initial data should satisfy condition (2).
Arguing as in Section 2, we see that faces 1 and N of the extension of u m and of u m do not disappear and that statement (9) holds, namely, even as time evolves and certain faces disappear, jumps in (u m x ) i correspond to adjacent corners ofW . Suppose now that we want u m to approximate the motion of u. Then we assume, in addition, that the number of faces of u m 0 is with c a constant. One can, for example, construct u m 0 using the method described on page 13 (basically, by taking the union of segments on lines tangent to u 0 and with admissible slopes) to get N ≤ [Total Variation of u ′ 0 ] /m. We want to estimate the H 1 norm of u − u m . The values of r 1 and r N have been chosen so that the sum of the boundary terms (i) through (vii) of the previous section vanish (see the Remark on page 13). So the estimate of the L 2 norm of u x − u m x given in the previous section remains valid. However, to control the H 1 norm of u − u m we need some additional information, since we cannot use Poincaré's inequality. For any g : [0, 1] → R we have with g = 1 0 g dx, so it suffices to control u − u m . The rate of change of the average of u is d dt withW as in Eq. (22). On the other hand, the rate of change of the average of u m is d dt since u m is continuous andẋ 0 =ẋ N = 0. Hence, by Eqs. (26) and (27), By assumption (35), (Observe that in this calculation we only need that 0 < r 1 , r N ≤ 1; if the distance between any two adjacent admissible slopes is a fixed constant then, with the choice of r 1 and r N made in Eqs. (34), one can improve this result to d The constant C depends only on max i | (u m x ) i | and on the constant appearing in hypothesis (35). If the initial data satisfies inequality (28) then max i | (u m x ) i | is uniformly bounded. If, in addition, the initial data satisfies inequality (29) then, by combining Eqs. (31), (36), and (37), we easily deduce that

Setup of the general Dirichlet problem
In this section we show how one can set up a crystalline algorithm for the general Dirichlet problem There is a new feature in this case: the time dependent boundary condition can lead to the creation of new faces at x = 0 and x = 1. This causes our convergence argument to break down, and, in fact, we do not prove that the algorithm converges. It seems natural to choose instead of Eqs. (3), so the boundary conditions are satisfied.
We also want to insure statement (9), i.e. that even as time evolves and certain faces disappear, jumps in (u m x ) i correspond to adjacent corners ofW . To do this we impose the following condition: there should be an extension of u m for which Eqs. (7) and (8) if (u m t ) 1 = 0 (otherwise there is no restriction on l 1 ). These conditions might break down when a ′ changes sign or when equality holds in inequality (39). More specifically: • There might not exist a (u m x ) 0 adjacent to (u m x ) 1 such that ∆ 1 and a ′ have the same sign. Under these circumstances, we have to allow a new face to appear at the boundary. The introduction of the new face makes it possible to meet the requirements above; it is necessary in two cases: Case (i): If a ′ is zero and about to become negative and (u m x ) 2 > (u m x ) 1 , or a ′ is zero and about to become positive and (u m x ) 2 < (u m x ) 1 , we introduce a new face with slope (u m x ) 2 , now (u m x ) 3 (see Figure 3). With the new ordering (u m x ) 1 > (u m x ) 2 and (u m x ) 1 < (u m x ) 2 , respectively, enabling us to pick the new (u m x ) 0 uniquely. Case (ii): If equality holds in inequality (39) andl 1 is positive we introduce a new face with slope (u m x ) 0 , now (u m x ) 1 (see Figure 4). Note that when a new face appears it has zero length so thereafter Eq. (39) is satisfied during some nonzero time interval.
is not zero, otherwise it is not determined uniquely). We take We can take l1 , as ∆ 1 = 0. Note that for the extension of u m we considered in Section 4 (odd and periodic) the ratio of the length of the first face of u m to the length of the first face of the extension of u m is one half, but we could have extended u m so that this ratio was one. Eq. (39) says that r 1 ≤ 1. The length of the first face of the extension of u m is l 1 / r 1 .
With the r i 's as in the previous paragraph we can estimate the growth of the H 1 norm of u − u m as was done in Section 4. Now the boundary terms do not vanish, but using Eqs. (3 ′′ ) we see that they add up to zero (see the Remark on page 13). However, our argument also used that max 1≤i≤N | (u m x ) i | stays bounded as m → 0, for 0 ≤ t ≤ T . It seems plausible that this should be the case if the solution of the differential equation (38) is such that sup 0≤t≤T || u x || ∞ is finite. This quantity is finite if, for example, we assume that for some c 1 , c 2 , and c 3 > 0.  [19] just mentioned are satisfied with their m equal to two instead of zero, since in this case W ′′ (u x ) = 1 / 1 + u 2 x .] Conditions (40) hold for motion by weighted curvature, as we verify in the Appendix.
The proof of Section 4 would go through if one could bound max 1≤i≤N | (u m x ) i | for 0 ≤ t ≤ T by a constant independent independent of m.

Appendix: Physical and mathematical context
We summarize here the relation between this work and the literature on surface energy driven motion of phase boundaries, especially the papers by Angenent and Gurtin [2] and Taylor [25]. Consider an interface between two phases moving isothermally according to the balance of capillary forces and constitutive equations compatible with thermodynamics. An evolution equation for the interface is derived in [2]. When the phases have the same energy and the kinetic coefficient (which measures the drag opposing interfacial motion) is one, it has the form (see Eq. (4.11) of [2]) where V is the normal velocity of the smooth interface and K is its curvature, θ is the angle from a fixed coordinate axis to the normal to the interface, and f (assumed smooth) is the interfacial energy per unit length. So fix a coordinate system. Let u( · , t) be the interface at time t and We are led to consider Eq. (1) by taking because The function f is the energy per unit length (of the interface) whereas W is energy per unit length of the projection of the interface on the x-axis. Note that the right hand side of Eq. (41) is the negative of the gradient of Taylor calls the negative of the gradient of E the weighted curvature of the interface (see Sections 2.2 and 2.3 of [25]). The interfacial energy per unit length, f , is usually represented in a polar diagram for 1/f , called Frank diagram. The function f is said to be strictly stable when f + f ′′ > 0. This condition corresponds to a strictly convex Frank diagram, and to W ′′ > 0. On the other hand, f is usually said to be crystalline if its convexified Frank diagram is a polygon, and if the vertices of this polygon form the complete set of globally convex sections of the Frank diagram (i.e. the original diagram and the convexified one meet only at the vertices of the latter) (see Section 10.3 of [2]). We prefer a slightly broader definition. We shall call such an energy strictly crystalline and do not require that a crystalline energy satisfy the second condition. Hence, if the Frank diagram is a polygon then f is crystalline, but not strictly crystalline. In general, we denote byf the function whose Frank diagram is the convexification of the Frank diagram of f .
There is an equivalent characterization of crystalline energies. It is obtained as follows. Let n(θ) △ = (cos θ, sin θ). For any surface energy f , we may extend f to R 2 as a homogeneous function of degree one, Recall that the Fenchel transform of f 0 is the function f * 0 : R 2 −→ [−∞, +∞] given by and that the Wulff set of f is The function f is crystalline if W f is polyhedral, and it is strictly crystalline if in addition f (θ) > f * * 0 (n(θ)) (= sup y∈W f {y · n(θ)}) unless n(θ) is normal to ∂W f (see, for example, Fonseca [9]). One can check thatf (θ) = f * * 0 (n(θ)). An evolution equation for an interface with a strictly crystalline energy is derived in [2] using the same physical laws which gave Eq. (41). For such an energy Eq. (41) is backward-parabolic. Therefore, one restricts crystalline interfaces to a space M, consisting of continuous piecewise linear functions such that the normal to each face makes an angle with the x-axis corresponding to one of the corners of the polygon 1/f × n. The normal velocity of a face is (see Eq. (10.12) of [2]) where∆ and σ are constants for each face, and L is the length of the face. For a face with normal n(θ)∆ = [f ′ (θ + 0) −f ′ (θ − 0) ], and σ = +1 if θ increases across the face, σ = −1 if θ decreases across the face, and σ = 0 if the face is nontransitional. Geometrically,∆ is the length of the segment in the Wulff set of f with normal n(θ). Using this formula, one can compute the velocity V i in terms of f , the θ i 's, and L i (the subindex i refers to the ith face): and l i = L i sin θ i (we recall that l i is the length of the projection of the ith face on the x-axis), one easily checks that this is Eq. (8), i.e. ∆ l = σ∆ L . Notice that wheref is as above andW is as in Section 2. Hence the right hand side of Eq. (44) is the negative of the gradient of on M, in other words the weighted curvature of the interface (see Sections 4.2 and 4.3 of [25]).
In summary, we have checked that if W is given by Eq. (42), Eq. (1) is Eq. (41) (with initial data and boundary conditions), and Eq. (8) is Eq. (44). Our approximation scheme is crystalline in the sense that it approximates the motion of an interface with a strictly convex energy by the motion of an interface with a strictly crystalline energy and in the sense thatf is crystalline.
In Sections 4 and 5, we proved convergence of the crystalline approximation scheme for homogeneous Dirichlet and Neumann boundary conditions, respectively. In such cases it follows from the maximum principle that || u x || ∞ is bounded. In Section 6 we discussed the general Dirichlet problem, and noted that growth conditions of the form (40) are sufficient to prove L ∞ bounds on u x . Let us verify that these conditions hold when W is determined by Eq. (43), with f strictly stable and smooth. Then condition (40.1) is obvious and condition (40.2) follows from the identity .