A novel finite element approximation of anisotropic curve shortening flow

We extend the DeTurck trick from the classical isotropic curve shortening flow to the anisotropic setting. Here the anisotropic energy density is allowed to depend on space, which allows an interpretation in the context of Finsler metrics, giving rise to e.g.\ geodesic curvature flow in Riemannian manifolds. Assuming that the density is strictly convex and smooth, we introduce a novel weak formulation for anisotropic curve shortening flow. We then derive an optimal $H^1$--error bound for a continuous-in-time semidiscrete finite element approximation that uses piecewise linear elements. In addition, we consider some fully practical fully discrete schemes and prove their unconditional stability. Finally, we present several numerical simulations, including some convergence experiments that confirm the derived error bound, as well as applications to crystalline curvature flow and geodesic curvature flow.


Introduction
The aim of this paper is to introduce and analyze a novel approach to approximate the evolution of curves by anisotropic curve shortening flow.The evolution law that we consider arises as a natural gradient flow for the anisotropic, spatially inhomogeneous, energy for a closed curve Γ, with unit normal ν, that is contained in the domain Ω ⊂ R 2 .In the above, γ : Ω×R 2 → R ≥0 denotes the anisotropy function and a : Ω → R >0 is a positive weight function.
In the spatially homogeneous case, i.e.
γ(z, p) = γ 0 (p) and a(z) = 1 the corresponding functional E frequently occurs as an interfacial energy, e.g. in models of crystal growth, [44,31].Our more general setting is motivated by the work [14] of Bellettini and Paolini, who consider the gradient flow for a perimeter functional P φ that is associated with a Finsler metric φ.Equations (2.5), (2.6) in [14] show that E(Γ) = P φ (Γ) if one chooses γ as the dual of φ and a in terms of the 2-dimensional φ-volume.As an important special case we mention that the choices γ(z, p) = G −1 (z)p • p and a(z) = det G(z) (1.3) can be used to describe the length of a curve in a two-dimensional Riemannian manifold (M, g).
In this case G(z) is the first fundamental form arising from a local parameterization of M, cf.Example 2.2(c) below.Apart from being of geometric interest, the functional E also has applications in image processing, [17].
The natural gradient flow for the energy E evolves a family of curves Γ(t) ⊂ Ω according to the law where V γ and κ γ are the anisotropic normal velocity and the anisotropic curvature, respectively.The precise definitions of these quantities are based on a formula for the first variation of E, and they will be given in Section 2, below.In the isotropic case, i.e. when γ(z, p) = |p| and a(z) = 1 we have that (1.4) is just the well-known curve shortening flow, V = κ, with V and κ denoting the normal velocity and the curvature of Γ(t), respectively.For theoretical aspects of the anisotropic evolution law (1.4) we refer to [13,30].Further information on (spatially homogeneous) anisotropic surface energies and the corresponding gradient flow can be found in [21,30] and the references therein.
In this paper we are interested in the numerical solution of (1.4) based on a parametric description of the evolving curves, i.e.Γ(t) = x(I, t) for some mapping x : I × [0, T ] → Ω.Here most of the existing literature has focused on the spatially homogeneous case (1.2).Then the law (1.4) reduces to 1 γ 0 (ν) V = κ γ 0 , see Section 2 for details, which can be viewed as a special case of the weighted anisotropic curvature flow for some mobility function β 0 , see e.g.[21, (8.20)].In [24], a finite element scheme is proposed and analyzed for (1.6) with β 0 = 1.The method uses a variational formulation of the parabolic system and is generalized to higher codimension in [38].A drawback of this approach is that the above system is degenerate in tangential direction and that the numerical analysis requires an additional equation for the length element.A way to circumvent this difficulty consists in replacing (1.7) by a strictly parabolic system with the help of a suitable tangential motion, known as DeTurck's trick in the literature.For the isotropic case, corresponding schemes have been suggested and analyzed in [19,26].We mention that alternative parametric approaches for (1.6) allow for some benign tangential motion, see e.g.[34,36,4,8].Since the choice (1.3) allows to describe the length of a curve in a Riemannian manifold (M, g), it is possible to use (1.4) in order to treat geodesic curvature flow in M within our framework.
Existing parametric approaches for this flow on a hypersurface of R 3 include the work [35] for the flow on a graph, as well as [6] for the case that the hypersurface is given as a level set.
Moreover, numerical schemes for the geodesic curvature flow (and other flows) of curves in a locally flat two-dimensional Riemannian manifold have been proposed in [12].To the best of our knowledge, no error bounds have been derived for a numerical approximation of geodesic curvature flow in Riemannian manifolds in the literature so far.
In this paper we propose and analyze a new method for solving (1.4), which is based on DeTurck's trick, and which applies to general, spatially inhomogeneous anisotropies.Let us outline the contents of the paper and describe the main results of our work.By taking advantage of the fact that the function (z, p) → 1 2 a 2 (z)γ 2 (z, p) is strictly convex in p, we derive in Section 3 a strictly parabolic system whose solution satisfies (1.4).It turns out that this system can be written in a variational form, which makes it accessible to discretization by linear finite elements.In the isotropic case, the resulting numerical scheme is precisely the method proposed and analyzed in [19], while in the anisotropic case we obtain a novel scheme that can be considered as a generalization of the ideas in [19] and [26].As one of the main results of this paper we show in Section 4 an optimal H 1 -error bound in the continuous-in-time semidiscrete case.Unlike in [24] and [38], the corresponding proof does not need an equation for the length element because of the strict parabolicity of the underlying partial differential equation.In order to discretize in time, we use the backward Euler method.In particular, in Section 5, as another important contribution of our work, we introduce unconditionally stable fully discrete finite element approximations for the following scenarios: In particular, the functions in (b) can be used to approximate the case of a crystalline anisotropy, cf.[4].Using these three fully discrete schemes, we present in Section 6 results of test calculations that confirm our error bound and show that the tangential motion that is introduced in our approach has a positive effect on the distribution of grid points along the discrete curve.
As our approach is based on a parameterization of the evolving curves, we only briefly mention numerical methods that employ an implicit description such as the level-set method or the phasefield approach.The interested reader may consult [37,10] for anisotropic curve shortening flow, as well as [16,15,43] for the geodesic curvature flow.These papers also provided additional references.
We end this section with a few comments about notation.Throughout, we let I = R/Z denote the periodic interval [0, 1].We adopt the standard notation for Sobolev spaces, denoting the norm of Furthermore, throughout the paper C will denote a generic positive constant independent of the mesh parameter h, see below.At times ε will play the role of a (small) positive parameter, with C ε > 0 depending on ε, but independent of h.Finally, in this paper we make use of the Einstein summation convention.

Anisotropy and anisotropic curve shortening flow
Let Ω ⊂ R 2 be a domain, and let a ∈ C 2 (Ω, R >0 ).Moreover, we assume that γ which means that γ is positively one-homogeneous with respect to the second variable.It is not difficult to verify that (2.1) implies that Here γ p = (γ p j ) 2 j=1 and γ pp = (γ p i p j ) 2 i,j=1 denote the first and second derivatives of γ with respect to the second argument.Similarly, we let γ z = (γ z j ) 2 j=1 denote the derivatives of γ with respect to the first argument.We note for later use that on differentiating (2.2) with respect to z we immediately obtain that the functions γ z j (z, •) and γ pz j (z, •) are positively one-and zero-homogeneous, respectively, for every z ∈ Ω.In addition, we assume that p → γ(z, p) is strictly convex for every z ∈ Ω in the sense that We are now in a position to define anisotropic curve shortening flow.To this end, with the help of Corollary 4.3 in [22], we first state the first variation of the functional E in (1.1), the proof of which will be given in Appendix A.
Lemma.2.1.Let Γ ⊂ Ω be a smooth curve with unit normal ν, unit tangent τ and scalar curvature κ.Let V be a smooth vector field defined in an open neighbourhood of Γ.Then the first variation of E at Γ in the direction V is given by where denote the anisotropic normal and the anisotropic curvature of Γ, respectively.
We remark that the definitions in (2.5) correspond to (3.5) and (4.1) in [14].Note also that ν γ is a vector that is normal to Γ, but normalized in such a way that γ(z, ν γ (z)) = 1, z ∈ Γ.We remark that although κ γ clearly depends on both γ and a, we prefer to use the simpler notation that drops the dependence on a.
Following [14, (1.1)], we now consider a natural gradient flow induced by (2.4).In particular, given a family of curves (Γ(t)) t∈[0,T ] in Ω, we say that Γ(t) evolves according to anisotropic curve shortening flow, provided that where V, with V denoting the normal velocity of Γ(t), and where κ γ is defined in (2.5).We remark that the name of the flow is inspired by the fact that solutions of (2.6) satisfy the energy relation We note that the higher dimensional analogue of (2.6) is usually called anisotropic mean curvature flow or anisotropic motion by mean curvature.Hence alternative names for the evolution law (2.6) in the planar case treated in this paper are anisotropic curvature flow and anisotropic motion by curvature. Example.2.2.
Let F : Ω → M be a local parameterization of M, {∂ 1 , ∂ 2 } the corresponding basis of the tangent space T F (z) M and g ij (z) = g F (z) (∂ i , ∂ j ), z ∈ Ω, as well as G(z) = (g ij (z)) 2 i,j=1 .We set γ(z, p) = G −1 (z)p • p and a(z) = det G(z).Then we have ) denotes an anti-clockwise rotation of p by π 2 .For a curve Γ ⊂ Ω the vector τ = −ν ⊥ then is a unit tangent and is the Riemannian length of the curve Γ = F (Γ) ⊂ M. We show in Appendix B that the geodesic curvature of Γ at F (z) is equal to κ γ at z ∈ Γ, and also that (Γ(t)) t∈[0,T ] ⊂ Ω is a solution of (2.6), if and only if Γ(t) = F (Γ(t)) evolves according to geodesic curvature flow in M.

DeTurck's trick for anisotropic curve shortening flow
In what follows we shall employ a parametric description of the evolving curves.Let Γ(t) = x(I, t), where x : I × [0, T ] → R and I = R/Z.In order to satisfy (2.6) we require that From now on we fix a normal on Γ(t) induced by the parameterization x, and, as no confusion can arise, we identify ν • x with ν, κ γ • x with κ γ and similarly κ • x with κ.In particular, we define the unit tangent, the unit normal and the curvature of Γ(t) by In place of (3.1) we simply write Clearly, (3.3)only prescribes the normal component of the velocity vector x t , and so there is a certain freedom in the tangential direction.Our aim is to introduce a strictly parabolic system of partial differential equations for the parameterization x, whose solution in normal direction still satisfies (3.3).
By way of motivation, let us briefly review the DeTurck trick in the isotropic setting, recall (1.5) and Example 2.2(a).Then (3.3) collapses to x t • ν = κ, and adjoining a zero tangential velocity leads to the formulation ) ρ , as the isotropic equivalent to (1.7).We recall that optimal error bounds for a semidiscrete continuous-in-time finite element approximation of this formulation have been obtained in the seminal paper [23] by Gerd Dziuk.One difficulty of Dziuk's original approach is that the analysed system is degenerate in the tangential direction.DeTurck's trick addresses this problem by removing the degeneracy through a suitable reparameterization.In fact, it is natural to consider the system recall (3.2), for which a semidiscretization by linear finite elements was analyzed in [19].The appeal of this approach is that the analysis is very elegant and simple.For example, the weak formulation of (3.4) is given by and choosing η = x t immediately gives rise to the estimate which can be mimicked on the discrete level.
Our starting point for extending DeTurck's trick to the anisotropic setting is to define the function Φ : We mention that the square of the anisotropy function plays an important role in the phase field approach to anisotropic mean curvature flow, cf.[27,1,10].
On noting (3.2), (2.1) and (3.7) we compute, similarly to (3.6), that The crucial idea is now to define positive definite matrices H(z, p) ∈ R 2×2 , for (z, p) ∈ Ω × (R 2 \ {0}), such that if a sufficiently smooth x satisfies then x is a solution to anisotropic curve shortening flow, (3.3).For the construction of the matrices H it is important to relate the right hand side in (3.9) to the right hand side in (3.3).We begin by calculating where we use the notation γ ⊥ p (z, p) = [γ p (z, p)] ⊥ .Furthermore, we obtain with the help of (2.1), (2.2) and (3.2) that and therefore and a similar argument for the second component yields we derive Let us now assume that x is a solution of (3.9), where H(z, p) is an invertible matrix of the form In order to determine α(z, p), β(z, p) ∈ R such that x satisfies (3.3), we calculate If we multiply by α 2 (x, x ρ ) + β 2 (x, x ρ ) and insert (3.14) we obtain, on noting , and so where in the last step we have used the one-homogeneity of γ, recall (2.2).Clearly, (3.3) will now hold if we choose In summary, we have shown the following result.
Furthermore, it can be shown that the system (3.9) is strictly parabolic.The proof hinges on the fact that H and Φ pp are positive definite matrices in Ω × (R 2 \ {0}).This property of Φ pp immediately follows from our convexity assumptions on γ, recall (2.3).
(3.18) Here, [p, q] ⊂ R 2 denotes the line segment connecting p and q.
Proof.It is shown in [30,Remark 1.7.5] that (2.3) implies that Φ pp (z, p) is positive definite for all z ∈ Ω and p = 0.The bound (3.17) then follows with the help of a compactness argument, while the elementary identity together with (3.17), implies (3.18).Lemma.3.3.The system (3.9) is parabolic in the sense of Petrovsky.Proof.On inverting the matrix H(x, x ρ ) we may write (3.9) in the form Hence, by definition we need to show that the eigenvalues of since H 11 > 0 and det H > 0, recall (3.15), and since A is symmetric positive definite in view of (3.17).Hence, either both eigenvalues are positive real numbers, or In view of Lemma 3.3 we expect that it is possible to prove the short-time existence of a unique smooth solution to (3.9).Moreover, existence and uniqueness of classical smooth solutions to PDEs of the form x t = b(κ, ν)ν + aτ , arising from closely related curvature driven geometric evolution equations, have been obtained in [34].
The weak formulation of (3.9) now reads as follows.Given x 0 : I → Ω, find x : Clearly (3.20) is the desired anisotropic analogue to (3.6).So together with the fact that x is a solution of the gradient flow (3.3), recall also (2.7), we obtain that both  (12) in [19].
(b) Space-independent anisotropy: We have Φ(z, p) = Φ 0 (p) = 1 2 γ 2 0 (p ⊥ ), so that Hence the weak formulation reads Hence the weak formulation reads Remark.3.5.It is a straightforward matter to extend our approach for (3.3) to the more general flow compare also with (1.6) in the space-independent case.In particular, it can be easily shown that if x is a solution to then it automatically solves (3.24).Extending our analysis in Section 4 to this more general case is straightforward, on making the necessary smoothness assumptions on β.

Finite element approximation
In order to define our finite element approximation, let 0 = q 0 < q 1 < . . .< q J−1 < q J = 1 be a decomposition of [0, 1] into intervals I j = [q j−1 , q j ].Let h j = q j − q j−1 as well as h = max 1≤j≤J h j .We assume that there exists a positive constant c such that so that the resulting family of partitions of [0, 1] is quasiuniform.Within I we identify q J = 1 with q 0 = 0 and define the finite element spaces Let {χ j } J j=1 denote the standard basis of V h .For later use, we let π h : C 0 (I) → V h be the standard interpolation operator at the nodes {q j } J j=1 , and we use the same notation for the interpolation of vector-valued functions.It is well-known that for k ∈ {0, 1}, ∈ {1, 2} and p ∈ [2, ∞] it holds that Our semidiscrete approximation of (3.19) is now given as follows.Find x h (q j , t)χ j we find that (4.2) gives rise to a system of ordinary differential equations (ODEs) in R 2J , which has a unique solution on some interval [0, T h ).By choosing η h = x h,t one also immediately obtains a semidiscrete analogue of (3.20).
In what follows we assume that (3.9) has a smooth solution x : Let S = x(I × [0, T ]).Then there exists δ > 0 such that B δ (S) ⊂ Ω and we define the compact set Theorem.4.1.Suppose that (3.9) has a smooth solution x : Then there exists h 0 > 0 such that for 0 < h ≤ h 0 the semidiscrete problem (4.2) has a unique solution x h : I × [0, T ] → Ω, and the following error bounds hold: Proof.Let us define . Arguing in a similar way for the first derivative, we deduce that Using the identity e = x − π h x + π h e (4.9) and choosing η h = π h e t in (4.8), we obtain Let us begin with the two terms on the left hand side of (4.10).Clearly, (3.16), (4.4), (4.5) and (4.7) imply that where Since p → Φ(z, p) and p → Φ z j (z, p) are positively homogeneous of degree 2, we have and therefore where the last equality follows from the relation Φ p (z, p) = Φ pp (z, p)p, recall (2.2).If we combine (4.12) with (4.11) and use a Taylor expansion together with (4.4) and (4.7), we obtain for the left hand side of (4.10), that Let us next estimate the terms on the right hand side of (4.10).To begin, we obtain from (3.15), (4.4), (4.5) and (4.1b) that The remaining terms involve differences between Φ p , H and Φ z , which will be estimated with the help of (4.4) and (4.7).Using (4.1b) we have as well as Finally, noting once again the identity (4.9) and the estimate (4.1b), we have where in the last inequality we have also used the embedding result (1.8).If we insert (4.13) and (4.14) into (4.10), and choose ε sufficiently small, we obtain where Clearly, we have from (3.18), on noting (4.9) and (4.1b), that and hence where Integrating (4.15) with respect to time, and observing (4.16) as well as µ(0) ≤ Ch 2 , we derive In particular, we have on recalling (1.8) that provided that 0 < h ≤ h 0 .Next, we have from (4.1b), (4.1a), (4.9) and (4.17) that and similarly by (4.3) that By choosing h 0 smaller if necessary we may therefore assume that ( x h,t 0,∞ ds ≤ 2C 0 , contradicting the definition of T h .Thus T h = T and the theorem is proved.

Fully discrete schemes
From now on, let the L 2 -inner product on I be denoted by (•, •).Due to the nonlinearities present in (4.2), for a fully practical scheme we need to introduce numerical quadrature.For our purposes it is sufficient to consider classical mass lumping.Hence for two piecewise continuous functions, with possible jumps at the nodes {q j } J j=1 , we define the mass lumped where u(q ± j ) = lim δ 0 u(q j ± δ).The definition (5.1) naturally extends to vector valued functions.
In particular, we will consider fully discrete approximations of in place of (4.2).Using this quadrature does not affect our derived error estimate (4.6), as can be shown with standard techniques.
In order to discretize (5.2) in time, let t m = m∆t, m = 0, . . ., M , with the uniform time step size ∆t = T M > 0. In the following, we let x 0 h = x h (•, 0) = π h x 0 ∈ V h and for m = 0, . . ., M − 1 let x m+1 h ∈ V h be the solution of a system of algebraic equations, which we will specify.In general, we will attempt to define fully discrete approximations that are unconditionally stable, in the sense that they satisfy the following discrete analogue of (3.20), (5.3) for k = 1, . . ., M .Here, the second inequality is a consequence of (3.16).
Lemma.5.1.A solution (x m h ) M m=0 to (5.4) satisfies the stability bound (5.3).A disadvantage of the scheme (5.4) is that at each time level a nonlinear system of equations needs to be solved.Following the approach in [20], see also [38], one could alternatively consider a linear scheme by introducing a suitable stabilization term and treating the elliptic term in (5.4) fully explicitly.

Proof. Choose η
If we restrict our attention to a special class of anisotropies, then a linear and unconditionally stable approximation can be introduced that does not rely on a stabilization term.This idea goes back to [4], and was extended to the phase field context in [10].In fact, a wide class of anisotropies can either be modelled or at least very well approximated by where Λ ∈ R 2×2 , = 1, . . ., L, are symmetric and positive definite, see [4,5].Hence the assumption (2.3) is satisfied.In order to be able to apply the ideas in [10,11], we define the auxiliary function φ 0 (p) = γ 0 (p ⊥ ), so that Φ 0 (p) = 1 2 γ 2 0 (p ⊥ ) = 1 2 φ 2 0 (p).Observe that φ 0 also falls within the class of densities of the form (5.5), i.e.
Moreover, we recall from [10,11] that Φ 0 (p) = B(p)p, if we introduce the matrices (5.6) On recalling once again the definition of H 0 from (3.21b), we then consider the scheme which is inspired by the treatment of the anisotropy in [10,11] and leads to a system of linear equations.We note that for the case L = 1 the two schemes (5.7) and (5.4) are identical.
Lemma.5.2.Suppose that x m h ∈ V h with x m h,ρ = 0 in I. Then there exists a unique solution x m+1 h ∈ V h to (5.7).Moreover, a solution (x m h ) M m=0 to (5.7) satisfies the stability bound (5.3).
Proof.Existence follows from uniqueness, and so we consider the homogeneous system: Find Choosing η h = X h , and observing that the matrices B(p) are positive definite, we obtain which implies that X h = 0 in view of (3.16).This proves the existence of a unique solution.In order to show the stability bound, we recall from Corollary 2.3 in [10] that with B as defined in (5.6), it holds that Hence choosing η h = x m+1 h − x m h in (5.7) and summing for m = 0 to k − 1 yields (5.3).

Curve shortening flow in Riemannian manifolds
Let us denote by Sym(2, R) the set of symmetric 2×2 matrices over R. For A, B ∈ Sym(2, R) we define A B if and only if A − B is positive semidefinite.We say that a differentiable function We now consider the situation in Example 3.4(c).Let G : Ω → Sym(2, R), so that Φ(z, p) = 1 2 G(z)p • p, where G(z) is positive definite for z ∈ Ω.In order to obtain an unconditionally stable scheme, we adapt an idea from [12] and assume that we can split G into Such a splitting exists if there exists a constant c G ∈ R ≥0 such that In that case one may choose (5.9) We now consider the scheme where H is as defined in (3.23).We note that in general (5.10) is a nonlinear scheme.

Proof. Let us again choose
h in (5.10) and calculate with the help of (5.9) ) is symmetric and positive definite.This yields the desired result similarly to the proof of Lemma 5.2.
In the special case that the Riemannian manifold is conformally equivalent to the Euclidean plane, i.e. when G(z) = g(z)Id for g : Ω → R >0 , several numerical schemes have been proposed in [12, §3.1].In this situation our fully discrete approximation (5.10) collapses to the new scheme where g = g + + g − and ±g ± : Ω → R are convex.We observe that for the special case g(z) = (z 1 ) 2 the scheme (5.11) is in fact very close to the approximation [3, (4.4)], modulo the different time scaling factor that arises in the context of mean curvature flow for axisymmetric hypersurfaces in R 3 considered there.

Numerical results
We implemented our fully discrete schemes within the finite element toolbox Alberta, [41].
Where the systems of equations arising at each time level are nonlinear, they are solved using a Newton method or a Picard-type iteration, while all linear (sub)problems are solved with the help of the sparse factorization package UMFPACK, see [18].For example, for the solution of (5.4) we employ a Newton iteration, while the Picard iteration for the solution of (5.10) is defined through x m+1,0 h = x m h and, for i ≥ 0, by In all our simulations the Newton solver for (5.4) converged in at most one iteration, while the Picard iteration (6.1) always converged in at most three iterations.
For all our numerical simulations we use a uniform partitioning of [0, 1], so that q j = jh, j = 0, . . ., J, with h = 1 J .Unless otherwise stated, we use J = 128 and ∆t = 10 −4 .On recalling (1.1), for χ h ∈ V h we define the discrete energy We also consider the ratio between the longest and shortest element of Γ m h = x m h (I), and are often interested in the evolution of this ratio over time.We stress that no redistribution of vertices was necessary during any of our numerical simulations.In the isotropic case this can be explained by the diffusive character of the tangential motion induced by (3.4), since the flow can be rewritten as x t = κν − ( 1 |xρ| ) ρ τ , as has been pointed out in e.g.[34, p. 1477].Our numerical experiments indicate that while the induced tangential motion from (3.9) may in general not be diffusive, it is sufficiently benevolent to avoid coalescence of vertices in practice.

Space-independent anisotropic curve shortening flow
In this subsection, we consider the situation from Example 2.2(b), see also Example 3.4(b).Anisotropies of the form γ(z, p) = γ 0 (p) can be visualized by their Frank diagram F = {p ∈ R 2 : γ 0 (p) ≤ 1} and their Wulff shape W = {q ∈ R 2 : γ * 0 (q) ≤ 1}, where γ * 0 is the dual to γ 0 defined by γ * 0 (q) = sup We recall from [28] that the boundary of the Wulff shape, ∂W, is the solution of the isoperimetric problem for E(Γ) = Γ γ 0 (ν) dH 1 .Moreover, it was shown in [42] that self-similarly shrinking boundaries of Wulff shapes are a solution to anisotropic curve shortening flow.In particular, solves (2.8).We demonstrate this behaviour in Figure 1 for the "elliptic" anisotropy Observe that (6.4) is a special case of (5.5) with L = 1, so that the scheme (5.4) collapses to the linear scheme (5.7).The evolution in Figure 1 nicely shows how the curve shrinks self-similarly to a point.We can also see that the scheme (5.4) induces a tangential motion that moves the initially equidistributed vertices along the curve, so that eventually a higher density of vertices can be observed in regions of larger curvature.We note that this behaviour is not dissimilar to the behaviour observed in the numerical experiments in [4].
Figure 1: Anisotropic curvature flow for (6.4) using the scheme (5.4).Solution at times t = 0, 0.05, . . ., 0.45, 0.499, as well as the distribution of vertices on x M h (I).Below we show a plot of the discrete energy E h (x m h ) and of the ratio (6.2) over time.1: Errors for the convergence test for (6.5) over the time interval [0, 0.45] for the scheme (5.4) with the additional right hand side (f (t m+1 ), η h ) h from (6.6).We also display the experimental orders of convergence (EOC).
We now use the exact solution (6.3) to perform a convergence experiment for our proposed finite element approximation (4.2).To this end, we choose the particular parameterization and define so that (6.5) is the exact solution of (3.22) with the additional right hand side (f, η).Upon adding the right hand side (f (t m+1 ), η h ) h to (5.4), we can thus use (6.5) as a reference solution for a convergence experiment of our proposed finite element approximation.We report on the observed H 1 -and L 2 -errors for the scheme (5.4) for a sequence of mesh sizes in Table 1.Here we partition the time interval [0, T ], with T = 0.45, into uniform time steps of size ∆t = h 2 , for h = J −1 = 2 −k , k = 5, . . ., 9. The observed numerical results confirm the optimal convergence rate for the H 1 -error from Theorem 4.1.It is not difficult to verify that this anisotropy satisfies (2.3) if and only if δ < 1 k 2 −1 , see also [8, p. 27].In order to visualize this family of anisotropies, we show a Frank diagram and a Wulff shape for the cases k = 3 and k = 6, respectively, in Figure 2. We show the evolutions for anisotropic curve shortening flow induced by these two anisotropies, starting in each case from an equidistributed approximation of a unit circle, in Figure 3.Here we use the fully discrete scheme (5.4).We observe from the evolutions in Figure 3 that the shape of the curve quickly approaches the Wulff shape, while it continuously shrinks towards a point.It is interesting to note that the ratio (6.2) increases only slightly and then appears to remain nearly constant for the remainder of the evolution.
One motivation for choosing a sixfold anisotropy, as in (6.7) with k = 6, is its relevance for modelling ice crystal growth, see e.g.[7,9].Here it is desirable to choose a (nearly) crystalline anisotropy, which means that the Wulff shape exhibits flat sides and sharp corners.With the help of the class of anisotropies (5.5) this is possible.We immediately demonstrate how this can be achieved for a general k-fold symmetry, for even k ∈ N. On choosing L = k/2, we define the rotation matrix Q(θ) = cos θ sin θ − sin θ cos θ and the diagonal matrix D(δ) = diag(1, δ 2 ), and then let We visualize some Wulff shapes of (6.8) for L = 2, 3, 4 in Figure 4 and observe that these Wulff  We show the solution at times t = 0, 0.05, . . ., 0.35, the distribution of vertices on x M h (I), as well as the evolution of the ratio (6.2) over time.
shapes, for δ → 0, will approach a square, a regular hexagon and a regular octagon, respectively.Of course, the associated crystalline anisotropic energy densities, when δ = 0, are no longer differentiable, and so the theory developed in this paper no longer applies.Yet, for a fixed δ > 0 all the assumptions in this paper are satisfied and our scheme (5.7) works extremely well.As an example, we repeat the simulations in Figure 3 for the anisotropy (6.8) with L = 2 and δ = 10 −2 , now using the scheme (5.7).From the evolution shown in Figure 5 it can be seen that the initial curve assumes the shape of a smoothed square that then shrinks to a point.We also observe that after an initial increase the ratio (6.2) decreases and eventually reaches a steady state.The final distribution of mesh points is such that there is a slightly lower density of vertices on the nearly flat parts of the curve.
Inspired by the computations in [37, Fig. 6.1] we now consider evolutions for an initial curve that consists of a 3  2 π-segment of the unit circle merged with parts of a square of side length 2. For our computations we employ the scheme (5.7), with the discretization parameters J = 256 and ∆t = 10 −4 .The evolutions for the three anisotropies visualized in Figure 4 can be seen in Figure 6.We observe that the smooth part of the initial curve transitions into a crystalline shape, while the initial facets of the curve that are aligned with the Wulff shape simply shrink.The other facets disappear, some immediately and some over time, as they are replaced by facets aligned with the Wulff shape.Particularly interesting is the evolution of the left nonconvex corner in the initial curve, which shows three qualitatively very different behaviours for the three chosen anisotropies.
For the final simulations in this subsection we choose as initial data a polygon that is very similar to the initial curve from [2, Fig. 0].In their seminal work, Almgren and Taylor consider motion  by crystalline curvature, which is the natural generalization of anisotropic curve shortening flow to purely crystalline anisotropies, that is when the Wulff shape is a polygon.For motion by crystalline curvature a system of ODEs for the sizes and positions of all the facets of an evolving polygonal curve has to be solved.Here the initial curve needs to be admissible, in the sense that it only exhibits facets that also appear in the Wulff shape, and any two of its neighbouring facets are also neighbours in the Wulff shape.Hence the initial curve for the computations shown in Figure 7, for which we employed the scheme (5.7) with J = 512 and ∆t = 10 −4 , is admissible for an eightfold anisotropy, with a regular octagon as Wulff shape.Our simulation for the smoothed anisotropy (6.8) with L = 4 and δ = 10 −4 agrees remarkably well with the evolution shown in [2, Fig. 0].In fact, it is natural to conjecture that in the limit δ → 0, anisotropic curve shortening flow for the anisotropies (6.8) converges to flow by crystalline curvature with respect to the crystalline surface energies (6.8) with δ = 0. We stress that for the cases L = 2 and L = 3, when the Wulff shape is a square and a regular hexagon, respectively, the initial curve in Figure 7 is no longer admissible in the sense described above.As we only deal with the case δ > 0, our scheme (5.7) has no difficulties in computing the evolutions shown in Figure 7 for L = 2 and L = 3.We observe once again that new facets appear where the initial polygon is not aligned with the Wulff shape, while the admissible facets simply shrink.
Figure 8: Curvature flow in the hyperbolic plane, using the scheme (5.11).Solution at times t = 0, 0.02, . . ., 0.14.We also show plots of the discrete energy E h (x m h ) and of the ratio (6.2) over time.

Curve shortening flow in Riemannian manifolds
In this subsection we consider the setup from Example 2.2(c), see also Example 3.4(c).At first we look at the simpler case of a manifold that is conformally flat, so that we can employ the scheme (5.11).As an example we take G(z) = g(z)Id with g(z) = (z 1 ) −2 and note that with Ω = {z ∈ R 2 : z 1 > 0} we obtain a model for the hyperbolic plane, which is a two-dimensional manifold that cannot be embedded into R 3 , as was proved by Hilbert, [32], see also [40, §11.1].
From [12, Appendix A], and on noting Lemma B.2 in Appendix B, we recall that a true solution for (2.6), i.e. geodesic curvature flow in the hyperbolic plane, is given by a family of translating and shrinking circles in Ω: Γ(t) = a(t) 0 + r(t)S 1 , a(t) = e −t a(0), r(t) = r 2 (0) − a 2 (0) 1 − e −2t 1 2 , (6.9) with a(0) > r(0) > 0 and S 1 = {z ∈ R 2 : |z| = 1}.In Figure 8 we show such an evolution, starting from a unit circle centred at 2 0 , computed with the scheme (5.11), where, since g is convex in Ω, we choose g + = g.We observe that during the evolution the discrete geodesic length is decreasing, while the approximation to the shrinking circle remains nearly equidistributed throughout.At the final time T = 0.14 the maximum difference between r(T ) and |x M h (q j ) − a(T ) 0 |, for 1 ≤ j ≤ J, is less than 6 • 10 −3 , indicating that the polygonal curve Γ M h = x M h (I) is a very good approximation of the true solution Γ(T ) from (6.9).
For the remainder of this subsection we consider general Riemannian manifolds that are not necessarily conformally flat.An example application is the modelling of geodesic curvature flow on a hypersurface in R 3 that is given by a graph.In particular, we assume that

.10)
The induced matrix G is then given by G(z) = Id + ∇ϕ(z) ⊗ ∇ϕ(z), and the splitting (5.8) for the scheme (5.10) can be defined by G + (z) = G(z) + c ϕ |z| 2 Id and G − (z) = −c ϕ |z| 2 Id, with c ϕ ∈ R ≥0 chosen sufficiently large.In all our computations we observed a monotonically decreasing discrete energy when choosing c ϕ = 0, and so we always let G + = G.
We begin with a convergence experiment on the right circular cone defined by ϕ(z) = b|z| and Ω = R 2 \ {0} in (6.10), for some b ∈ R ≥0 .A simple calculation verifies that the family of curves Γ(t) = r(t)(S   flow on M = F (Ω).In fact, it is not difficult to show that the particular parameterization so that Γ(t) = F (x(I, t)), solves (3.9).Similarly to Table 1, we report on the H 1 -and L 2 -errors between (6.11), for b = √ 3 and r(0) = 1, and the discrete solutions for the scheme (5.10) in Table 2.Here for a sequence of mesh sizes we use uniform time steps of size ∆t = h 2 , for h = J −1 = 2 −k , k = 5, . . ., 9. Once again, the observed numerical results confirm the optimal convergence rate from Theorem 4.1.
On the same cone M, we perform two computations for a curve evolving by geodesic curvature flow.For the simulation on the left of Figure 9 it can be observed that as the initial curve F (x 0 h (I)) is homotopic to a point on M, it shrinks to a point away from the apex.On recalling Conjecture 5.1 in [29], due to Charles M. Elliott, on the right of Figure 9 we also show a numerical experiment for a curve that is not homotopic to a point on M. According to the conjecture, any such curve should shrink to a point at the apex in finite time, and this is indeed what we observe.
For the final set of numerical simulations, we model a surface with two mountains.Following Figure 10: Geodesic curvature flow on the graph defined by (6.12) with λ 1 = λ 2 = 1.We show the evolution of F (x m h ) on M at times t = 0, 1, 2, 2.2.
[45], we define and let Ω = R 2 .We show three evolutions for geodesic curvature flow on such surfaces in Figures 10,11 12.In each case we start the evolution from an equidistributed approximation of a circle of radius 2 in Ω, centred at the origin.In the first two simulations the curve manages to continuously decrease its length in R 3 , until it shrinks to a point.To achieve this in the second example, the curve needs to "climb up" the higher mountain.However, in the final example the two mountains are too steep, and so the curve can no longer decrease its length by climbing higher.In fact, the curve approaches a steady state for the flow, that is, a geodesic on M, i.e. a curve with vanishing geodesic curvature.The plot of the discrete energy in Figure 12 confirms that the evolution is approaching a geodesic.

A First variation of the anisotropic energy
Proof of Lemma 2.1.Abbreviating γ(z, p) = a(z)γ(z, p), (z, p) ∈ Ω × R 2 , we temporarily write E in (1.1) as Let us fix a curve Γ ⊂ Ω and a smooth vector field V defined in an open neighbourhood of Γ.We infer from Corollary 4.3 in [22] and (2.2) that the first variation of E(Γ) in the direction V Figure 12: Geodesic curvature flow on the graph defined by (6.12) with λ 1 = λ 2 = 5.We show the evolution of F (x m h ) on M at times t = 0, 1, 2, 4. Below we show a plot of F (x M h ) on M, as well as a plot of the discrete energy E h (x m h ) over time.
is given by Here we note that the differential operators ∂ ν f = f z i ν i and div Γ f = f i,z i − f i,z j ν j ν i on Γ only act on the first variable of functions defined in Ω × R 2 .In addition, we observe that the Weingarten map ∇ Γ ν is given by ∇ Γ ν = −κτ ⊗ τ .We then calculate, on noting (2.2), that If we insert the above relations into (A.1) and recall (2.5), we obtain which is (2.4).

B Geodesic curve shortening flow in Riemannian manifolds
In this appendix we prove the claims formulated at the end of Example 2.2(c).Here we will make use of standard concepts in Riemannian geometry, and we refer the reader to the textbook [33] for further details.Let F : Ω → M be a local parameterization of a two-dimensional Riemannian manifold (M, g) and denote by {∂ 1 , ∂ 2 } the corresponding basis of the tangent space T F (z) M, for z ∈ Ω.We also let g ij (z) = g F (z) (∂ i , ∂ j ), G(z) = (g ij (z)) 2 i,j=1 , (g ij (z)) 2 i,j=1 = G −1 (z), γ(z, p) = G −1 (z)p • p and a(z) = det G(z), for z ∈ Ω and p ∈ R 2 , which induces the energy equivalence (2.10).Let Γ be a smooth curve in M with unit tangent τ g and a unit normal ν g such that {τ g , ν g } is an orthonormal basis of the tangent space T M, i.e. g(τ g , τ g ) = g(ν g , ν g ) = 1 and g(τ g , ν g ) = 0. Then the geodesic curvature κ g of Γ is defined by We infer from (B.4a) that For notational convenience, we drop the dependences on x from now on.It is well-known that g kl,z i = g kr Γ r il + g lr Γ r ik , i, k, l = 1, 2. (B.6) Combining (B.6) with the relation (G −1 ) z i = −G −1 G z i G −1 , we find that [(G −1 ) z i ν] j = −g jk g kl,z i g lm ν m = −g jk g lm g kr Γ r il + g lr Γ r ik ν m = −g lm Γ j il ν m − g jk Γ m ik ν m , as well as (G −1 ) z i ν • ν = [(G −1 ) z i ν] j ν j = − g lm Γ j il ν m + g jk Γ m ik ν m ν j = −2g jk Γ m ik ν m ν j .
As a result, g lk g jr − g jk g lr Γ m jk ν m ν l ν r . Clearly, (a) a spatially homogeneous, smooth anisotropy function γ(z, p) = γ 0 (p); (b) a spatially homogeneous anisotropy function γ(z, p) = L =1 √ Λ p • p, where Λ are symmetric and positive definite matrices; (c) a spatially inhomogeneous anisotropy function of the form (1.3) to model geodesic curvature flow in a two-dimensional Riemannian manifold.
by • ,p and the seminorm by | • | ,p .For p = 2, W ,2 (I) will be denoted by H (I) with the associated norm and seminorm written as, respectively, • and |•| .The above are naturally extended to vector functions, and we will write [W ,p (I)] 2 for a vector function with two components.For later use we recall the well-known Sobolev embedding H 1 (I) → C 0 (I), i.e. there exists C I > 0 such that