Numerical approximation of boundary value problems for curvature flow and elastic flow in Riemannian manifolds

We present variational approximations of boundary value problems for curvature flow (curve shortening flow) and elastic flow (curve straightening flow) in two-dimensional Riemannian manifolds that are conformally flat. For the evolving open curves we propose natural boundary conditions that respect the appropriate gradient flow structure. Based on suitable weak formulations we introduce finite element approximations using piecewise linear elements. For some of the schemes a stability result can be shown. The derived schemes can be employed in very different contexts. For example, we apply the schemes to the Angenent metric in order to numerically compute rotationally symmetric self-shrinkers for the mean curvature flow. Furthermore, we utilise the schemes to compute geodesics that are relevant for optimal interface profiles in multi-component phase field models.


Introduction
In this paper we consider numerical approximations for gradient flows of curves evolving in Riemannian manifolds that are conformally flat. Here we allow both closed and open curves, where in the latter case appropriate boundary conditions need be considered in order to respect the required gradient flow structure.
We define the Riemannian manifold with the help of its metric tensor as follows. On a domain H ⊂ R 2 we let the metric tensor be given by [( v, w) for z ∈ H , (1.1) where v . w = v T w is the standard Euclidean inner product, and where g : H → R >0 is a smooth positive weight function. This is the setting one obtains for a two-dimensional Riemannian manifold that is conformally equivalent to the Euclidean plane. Of course, if g is constant we recover the case of a Euclidean ambient space. In local coordinates the metric is precisely given by (1.1), see e.g. [46,58,50]. Examples of such situations are the hyperbolic plane, the hyperbolic disc and the elliptic plane. Other examples are given by two-dimensional manifolds in R d , d ≥ 3, that can be conformally parameterised, such as spheres without pole(s), catenoids and tori. Coordinates (x 1 , x 2 ) ∈ H together with a metric g as in (1.1) are called isothermal coordinates, i.e. in all situations considered in this paper we assume that we have isothermal coordinates. We refer to Section 2 and [50, 3.29 in Section 3D] for more information.
The metric tensor (1.1) induces a notion of length in H. In particular, the length of a vector v ∈ R 2 at the location z ∈ H is defined by whereas the length of a curve γ ∈ C 1 ([0, 1], H) is given by We note that L g , which is also called geodesic length, naturally induces a distance function in H, with the distance between two points z 0 , z 1 ∈ H defined as dist g ( z 0 , z 1 ) = inf L g ( γ) : γ ∈ C 1 ([0, 1], H) , γ(0) = z 0 , γ(1) = z 1 .
It can be shown that (H, dist g ) is a metric space, see [46,Section 1.4].
We will present the mathematical details of curvature flow and elastic flow in the next section, together with the derivation of suitable boundary conditions. For now we mention that curvature flow, for a family of curves (Γ(t)) t∈[0,T ] , can be defined as the L 2gradient flow of geodesic length, L g (Γ) = Γ g 1 2 ds, with respect to the L 2 -inner product v, w = Γ v w g 1 2 ds. It is often called curve shortening flow. In particular, on letting the geodesic curvature κ g be the first variation of L g (Γ), with respect to ·, · , then curvature flow is given by where V g = g 1 2 V, and V is the Euclidean normal velocity of Γ in R 2 . Moreover, the geodesic elastic energy of Γ is defined by W g (Γ) = 1 2 Γ (κ g ) 2 g 1 2 ds, and elastic flow arises as the L 2 -gradient flow of W g (Γ). The elastic energy is often called bending energy of Γ, and elastic flow is also known by the name curve straightening flow. Critical points of W g (Γ) are called elastica, and a lot of interest in elastic flow is driven by the fact that stationary solutions to the flow are elastica. Moreover, elastic flow can be a viable strategy to obtain unstable geodesics, i.e. critical points of the length L g (Γ) that are unstable under curvature flow. In fact, all geodesics satisfy κ g = 0, and so they represent global minimisers for the bending energy.
There is a tremendous amount of work in the literature on curvature flow and elastic flow in the Euclidean plane, both from an analytical and a numerical point of view. Curvature flow in more complex ambient spaces has been treated analytically in e.g. [37,20,1], while numerical approximations have been considered in [21,56,60,55,6,12], for the case of closed curves, and in [16] for the case of open curves. Using the flow along the negative gradient of the total squared geodesic curvature functional to obtain geodesics as long-time limits has been first suggested in a seminal paper by Langer and Singer, [51]. Later the same authors used Morse theory to investigate stable critical points of the functional, [52]. The curve straightening flow has been used in [53] to compute periodic geodesics, and in [48] a variant of the flow was used to study the evolution of an elastic wire in a Riemannian manifold. Short-time existence, and in certain cases long-time existence, for elastic flow for closed curves in Riemannian manifolds has been studied analytically in [24], for the case that the manifold is a sphere, and in [29,30,57] for the hyperbolic plane. We are not aware of existing work on boundary value problems for elastic curves in Riemannian manifolds, but note that the case of a Euclidean ambient space has been considered in e.g. [43,25,28,44,26,27,39,40]. As far as the numerical approximation of elastic flow is concerned, we remark that the case of a Euclidean ambient space have been treated in [32,8,15]. In [32] error estimates are shown, while [15] contains a partial convergence result under a regularity assumption on the velocity. The case of a Riemannian manifold has been studied in [6,12,10]. All these approaches use finite element discretisations and are of variational structure.
The present authors, together with John W. Barrett, have considered the evolution of closed curves in Riemannian manifolds that are conformally equivalent to the Euclidean plane in the recent works [12,10]. Building on these works, in this paper we derive boundary value problems for curvature flow and elastic flow in such manifolds, and we believe that for elastic flow the obtained formulations are new in the literature. Using appropriate variations, different boundary value problems are derived in Section 2.1 for curvature flow, see (2.22), and in Section 2.2 for elastic flow, see (2.29), (2.30) and (2.31). In the case of elastic flow, the obtained conditions generalise classical Navier conditions as well as so called clamped conditions and semi-free type conditions. We will also introduce variational formulations, which lead to natural spatially discrete and fully discrete approximations for the highly nonlinear problems. In particular, the variational treatment allows for a natural discretisation of boundary conditions, which in the case of elastic flow are highly non-trivial. We introduce finite element schemes with good mesh properties as well as schemes which allow for a stability result. We also present several numerical results, which include computations that are the first for boundary value problems for elastic flow in Riemannian manifolds. We believe that the presented numerical simulations make a strong case for the usefulness of curvature flow and elastic flow in computing geodesics in Riemannian manifolds, both in the case of closed curves, and in the case of curves with boundary conditions. We end this introduction with the presentation of some example metrics for (1.1). To this end, we define the half-plane with closure H 2 = { z ∈ R 2 : z . e 1 ≥ 0}. Metrics that the authors have considered in their recent works on closed curves, see [12,10], are Recall that (1.4a) with µ = 1 models the hyperbolic plane, while µ = 0 corresponds to the Euclidean plane. Additional metrics that we consider in this paper are g( z) = ( z . e 1 ) 2 (n−1) e − 1 2 | z| 2 , n ∈ Z ≥2 , and H = H 2 , (1.5a) (1.5c) The metric (1.5a) is also called the Angenent metric, see [2, (5)], with a small mistake in the exponent, and [47, (1.3)], and is of interest in differential geometry. Here we mention the fact that for a rotationally symmetric hypersurface S ⊂ R n+1 , with generating curve Γ ⊂ H 2 , the geodesic length of Γ, with respect to the metric (1.5a), collapses, up to a constant factor, to Huisken's F -functional see [45,23,17]. It can be shown, [45], that critical points of (1.6) are self-shrinkers for mean curvature flow, and so geodesics for the metric (1.5a) generate axisymmetric self-shrinkers, such as the Angenent torus, see [2].
Finally, metrics of the type (1.5c), together with the choices play a role in determining optimal interface profiles in multi-component Ginzburg-Landau phase field models, see e.g. [41]. For example, the choice where σ 12 , σ 13 , σ 23 ∈ R >0 and σ 123 ∈ R ≥0 , corresponds to [41, (24), (25)], where u = (u 1 , u 2 , u 3 ) T represents a three-phase order parameter and u i stands for the fraction of phase i. We recall that physically meaningful values for the order parameter have to lie within the Gibbs simplex (1.10) In order to rigorously relate phase field parameters to their sharp interface limits, it is necessary to establish if the only geodesics connecting the three pure phases, e 1 = (1, 0, 0) T , e 2 = (0, 1, 0) T and e 3 = (0, 0, 1) T , are given by straight line segments. Of course, generalisations to other types of potentials are also possible. We refer to [41,38,19] for more details.
The remainder of this paper is organised as follows. In Section 2 we present strong and weak formulations of curvature flow and elastic flow. The semidiscrete continuous-intime finite element approximations introduced in Section 3 will be based on these weak formulations. Stability of some of the schemes is also shown in Section 3. Fully discrete approximations are presented in Section 4, together with results on their well-posedness and stability, where applicable. Finally, in Section 5 we present several numerical simulations for the derived schemes and the various metrics considered in this paper.

Mathematical formulations
We let R/Z be the periodic interval [0, 1], and set Consider a family of curves (Γ(t)) t∈[0,T ] , T > 0, that can be either open, I = (0, 1), or closed, I = R/Z. Given some smooth parameterisation x : where · ⊥ denotes a clockwise rotation by π 2 . We let V = x t . ν denote the normal velocity, and let the Euclidean curvature κ be defined by see [33]. We also let We introduce so that τ g . ν g = 0 and | τ g | 2 g = | ν g | 2 g = ( ν g , ν g ) g = g( x) ν g . ν g = 1, and let At this stage we would like to draw the reader's attention to the different usages of subscripts in this paper. The subscripts · g above, and throughout the paper, denote quantities associated with the metric g. The subscripts · t and · ρ , on the other hand, denote partial derivatives with respect to t and ρ, respectively. Finally, · s and · sg denote weighted partial derivatives, and are defined in (2.1) and in (2.3), respectively.
The geodesic curvature can be defined as see [12]. We note that V g and κ g , similarly to V and κ, only depend on Γ, but not on the chosen parameterisation x. For the case of an evolving closed curve, ∂I = ∅, we recall from [12] that curvature flow, is the L 2 -gradient flow of the geodesic length,  Elastic flow, on the other hand, is the L 2 -gradient flow of the elastic energy It was shown in [51], see also [12], that for closed curves this flow is given by where the Gaussian curvature S 0 is defined by see, e.g., [49,Definition 2.4]. In particular, closed curves evolving by (2.11) satisfy We state the value of the Gaussian curvature S 0 for the metrics (1.4) and (1.5) in Table 1.
Here we note that for the Euclidean case, (1.4a) with µ = 0, the geodesic elastic flow (2.11) collapses to classical elastic flow V = −κ ss − 1 2 κ 3 . In the remainder of this section, we would like to derive suitable boundary conditions for curvature flow and elastic flow that respect the gradient flow structures (2.9) and (2.12), and then introduce weak formulations for the obtained boundary value problems. In general, in the case of an open curve, I = (0, 1), we would like to consider the following types of boundary conditions on ∂I:  (1.4a) with µ ≤ −1 and for (1.5a). Having boundary points move freely on the x 2 -axis in that case is of particular interest, most notably when the evolving curve plays the role of the generating curve for an axisymmetric surface. Altogether, and for later use, we consider the disjoint partition (2.14c) In the above ∂ 0 I denotes the subset of boundary points of I that correspond to endpoints of Γ(t) where g is set to vanish, and only in that case does it make sense to consider (2.14a) separately from (2.14c). I.e. from now on we will assume that g( x(0)) = 0 on ∂ 0 I, so that (2.14a) implies g( x(t)) = 0 on ∂ 0 I for all t. In our paper, we will consider (2.14a) only for (1.4a), with µ ≤ −1, and (1.5a). The subscripts D, C, N relate to Dirichlet, clamped and Navier boundary conditions, respectively, with the former relevant for curvature flow, and the latter two having applications for elastic flow. In Table 2 we visualise the different types of boundary nodes that we consider in this paper.
For some of the weak formulations, it will be useful to have an analogue of (2.2) for the geodesic curvature κ g , recall (2.6). To this end, we note that it can be easily shown from (2.2) that Let (·, ·) denote the L 2 -inner product on I, and let
It holds that where, we have recalled (2.6). Clearly, the curvature κ g is the first variation of the length (2.8).
For the case that ∂ 0 I = ∅, we note that in order for the right hand side of (2.19) to remain bounded, it is appropriate to require that the term ν . ∇ g 1 2 ( x) in the second line remains bounded as we approach ∂ 0 I. In view of with ∇ g( x) . e 2 = 0 on ∂ 0 I, since g = 0 on ∂ 0 I, we see that is a natural assumption to make. Moreover, in situations where the curve Γ(t) = x(I, t) models the generating curve of an axisymmetric surface that is rotationally symmetric with respect to the x 2 -axis, the condition (2.20) ensures that the modelled surface is smooth; see also [11,9] for more details.
Similarly to the closed curve case, as discussed in [12], we note from (2.19) that (2.7) is the natural L 2 -gradient flow of L g with respect to the metric induced by g, i.e. it satisfies (2.9), if the boundary conditions (2.14) hold, together with This condition holds automatically on ∂ 0 I and ∂ D I, while on the remainder of ∂I we require x ρ . e 2 = 0 on ∂ 1 I and x ρ . e 1 = 0 on ∂ 2 I .
In terms of the Euclidean properties of the curve Γ(t), geodesic curvature flow, i.e. the evolution equation (2.7), can be written as This formulation gives another interpretation for the boundary condition (2.20). In fact, imposing (2.20) is necessary to allow the right hand side of (2.21) to remain bounded as we approach ∂ 0 I. In this way, we restrict ourselves to the class of solutions to the evolution equation (2.21), where the normal velocity and curvature remain bounded. Hence in this paper, geodesic curvature flow for an open or closed curve is given by: We remark that the condition x ρ . e 3−i = 0 in (2.22d) corresponds to a 90 • contact angle condition, which is the same as for classical Euclidean curvature flow. In particular, it does not depend on the chosen metric g. That is because in this paper we consider conformal metrics, and so compared to the Euclidean case only the measurement of length changes, but the measurement of angles remains the same.
A weak formulation of curvature flow, (2.22), based on the strong formulation (2.21) in place of (2.22a), is given as follows, where we also recall (2.2) and (2.17).
We note that the boundary conditions for x t in (2.22) are enforced through the trial space X, recall (2.17), while the boundary conditions on x ρ in (2.22) follow from (2.23b).
An alternative weak formulation, based directly on (2.22), together with (2.16), is given as follows.
Once again, the boundary conditions for x t in (2.22) are enforced through the trial space X, while the conditions on x ρ in (2.22d) follow directly from (2.24b). In addition, it can be shown that (2.24b), for the metrics (1.4a), with µ ≤ −1, and (1.5a), also enforces the condition on x ρ in (2.22b). In fact, this follows by using the techniques in [11, Appendix A], and noting that for both metrics it holds that ∇ g 1 2 ( x) . e 2 = 0 on ∂ 0 I, see Table 3 below.

Elastic flow
For elastic flow we assume that ∂I = ∂ 0 I ∪ ∂ C I ∪ ∂ N I, where, as before, ∂ 0 I will only be nonempty for the metrics (1.4a), with µ ≤ −1, and (1.5a).
In order to discuss appropriate boundary conditions consistent with the gradient flow structure (2.12), we need to re-visit the derivation of recall (2.11), as presented in [12, §2.3]. Summarising the authors' procedure there, they inferred by careful calculation that for closed curves, cf. [12, (2.55)], which concurs with the result of [51, p. 3]. In order to generalise their work to the case of open curves, we observe that boundary terms on the right hand side of (2.25) would only be created through applications of integration by parts. In the derivation in [12, §2.3] this occurs only in the third line of (2.47) and in the second line of (2.52). Regarding the former, we note that for open curves we obtain for the term in question that In addition, the three applications of integration by parts in [12, (2.52)] give rise to the following boundary terms: Hence, overall, we obtain the boundary terms An alternative way to derive the boundary terms uses the approach of Langer and Singer, see [51, p. 3]. We now derive natural conditions that make (2.26) vanish for the boundary conditions considered in (2.14). On ∂ C I ∪ ∂ N I, the first two terms vanish. On ∂ N I we require κ g = 0 to make the third term zero, while the clamped boundary conditions Lemma 37(ii) in [13]. For the compact notation in (2.27) we have used the fact that the identity function id always maps elements of ∂I to either 0 or 1. On ∂ 0 I the first term in (2.26) is zero, and on requiring we can make the second term vanish. The third term vanishes if V s = 0 on ∂ 0 I, which similarly to the clamped boundary conditions follows from ensuring the natural assumption (2.20).
Overall, we obtain the following strong formulation for elastic flow for open or closed curves, consistent with the gradient flow structure (2.12).
Other types of boundary conditions, corresponding to so-called free and semi-free boundary nodes, see e.g. [14], are also possible. For example, in the semi-free cases one could require These conditions will also lead to vanishing boundary terms in (2.26). Observe that (2.30) involves a 90 • contact angle condition, while (2.31) does not fix the angle but requires curvature to be zero, similarly to the Navier condition (2.29d). The boundary condition [κ g ] ρ − 1 2 κ g [ln g( x)] ρ = 0 seems to be completely new in the literature, as well as the various combinations of it with other boundary conditions. In order to simplify the presentation of what follows, we will concentrate on the conditions in (2.29) from now on. However, using the techniques from [8,14] it is straightforward to extend our weak formulations and finite element approximations to these other types of boundary conditions as well.
Combining the techniques in [10,14], it is not difficult to derive the following two weak formulations of elastic flow. Here the equations in the interior of I are the same as for (P) and (Q) in [10], while the treatment of the boundary conditions, achieved through the spaces X and Y from (2.17), is very similar to the approach taken in [14] in the context of axisymmetric Willmore flow.
The consistency of (2.32), in the case I = R/Z, was shown in [10, Appendix A.1]. As far as the boundary conditions are concerned, we note that the conditions on x t in (2.29) are enforced through the trial space X, recall (2.17). The second conditions in (2.29b) and (2.29c), respectively, are both enforced through (2.32c). In addition, we note from (2.32b) and (2.6) that y . ν = κ g , and so the trial space Y yields the second condition in (2.29d). It remains to validate that (2.32a) weakly enforces the third condition in (2.29b). This can be achieved on closely following the argument in [10, (A.3)-(A.9)], noting in particular that the integration by parts in [10, (A.5)] gives rise to the boundary term (κ g ) s − 1 2 κ g (ln g( x)) s on ∂ 0 I, which enforces (2.28), as required.
The consistency of (2.33), in the case I = R/Z, was shown in [ τ on ∂ 0 I, which thanks to y g ∈ Y and x ρ . e 1 = 0 collapses to enforcing (2.28).

Semidiscrete finite element approximations
, be a decomposition of [0, 1] into intervals given by the nodes q j , I j = [q j−1 , q j ]. For simplicity, and without loss of generality, we assume that the subintervals form an equipartitioning of [0, 1], i.e. that Clearly, if I = R/Z we identify 0 = q 0 = q J = 1. In addition, we let q J+1 = q 1 .
The necessary finite element spaces are defined as follows: In addition, we define (2.17). We define the mass lumped L 2 -inner product (u, v) h , for two piecewise continuous functions, with possible jumps at the nodes where we define u(q ± j ) = lim δ 0 u(q j ± δ). The definition (3.1) naturally extends to vector valued functions. Moreover, let (·, ·) denote a discrete L 2 -inner product based on some numerical quadrature rule. In particular, for two piecewise continuous functions, with possible jumps at the nodes {q j } J j=1 , we let (u, v) = I (u v), where 2) with K ≥ 2, K k=1 w k = 1, and with distinct α k , k = 1, . . . , K. A special case is (·, ·) = (·, ·) h , recall (3.1), but we also allow for more accurate quadrature rules.
Let ( X h (t)) t∈[0,T ] , with X h (t) ∈ V h , be an approximation to ( x(t)) t∈[0,T ] . Then, similarly to (2.1), we set For later use, we let ω h ∈ V h be the mass-lumped

Curvature flow
We consider the following finite element approximation of (A), recall (2.23). It is closely related to the approximation [12, (3.3), (3.10)] for closed curve evolutions.
To motivate the choice (3.6), we observe that it follows from (2.20), Table 3 and L'Hospital's rule that as ρ → ρ 0 ∈ ∂ 0 I . On recalling (3.4), we note that replacing ω h with ν h in (3.5b) does not change the scheme. On the other hand, for nonconstant g we must use ω h in the first term in (3.5a) to allow for an existence and uniqueness proof on the fully discrete level, see Lemma 4.1 below. In addition, as (3.6) is defined vertex-wise, we use the vertex normals ω h there.
It does not appear possible to prove a stability result for the scheme (3.5). However, thanks to (3.5b) the scheme (A h ) h satisfies a weak equidistribution property, i.e. it can be shown that neighbouring elements of Γ h (t) = X h (I, t) have the same length if they are not parallel. In effect, the side condition (3.5b) induces some tangential motion, which ensures that the equidistribution property holds at all times. The authors first introduced and discussed schemes with such an implicit tangential motion in [5] and have used it in a series of works since. We refer to the recent review article [13] for more details on that aspect of the scheme.
As an alternative approximation, we propose the following finite element approximation of (C), recall (2.24). It is the natural extension to the open curve case of the semidiscrete analogue of [12, (3.5), (3.18)].
Observe that we fix κ h g (t) to be zero on ∂ 0 I, since these values can be arbitrary. Setting them to zero allows for a uniqueness result on the fully discrete level. Note furthermore that on replacing ν h with ω h we obtain a related, but different, approximation of (C), that will have the analogous theoretical properties. We prefer the more natural variant as stated, in agreement with the closed curve scheme from [12]. Similarly to (A h ) h , the scheme (C h ) h also exhibits some implicit tangential motion, but, in contrast to (A h ) h , it does not appear possible to derive rigorous results on it. However, the scheme does admit a stability bound. To formulate this result, and on recalling (2.8), for Z ∈ V h we let Then we can prove the following discrete analogue of (2.9) for the scheme (3.7).
Proof. Choosing χ = κ h g ∈ W h ∂ 0 in (3.7a) and η = X h t ∈ X h in (3.7b) yields and so we obtain the desired result (3.9).
The identity (3.9), after integration in time, yields a stability bound for the discrete length, and it is a direct discrete analogue of (2.9).

Elastic flow
Following the approach by the authors in [10], it is straightforward to derive a semidiscrete approximation of (P), in the case that ∂ 0 I = ∅. The derived scheme, which will correspond to [10, (P h ) h ] with the natural changes to the test and trial spaces, can be shown to be stable. Moreover, and similarly to (A h ) h , the derived scheme will satisfy an equidistribution property. However, as it appears to be highly nontrivial to extend the approximation to the case ∂ 0 I = ∅, we do not pursue this variant any further in this paper.
On the other hand, the following discretisation based on the formulation (Q) can naturally deal with all the considered boundary conditions. It is the natural extension to the open curve case of the semidiscrete scheme [10, (Q h ) ].
If we replace ω h with ν h in (3.10a) we obtain a related, but different, approximation of (Q) with the analogous theoretical properties. We prefer the scheme as stated because it simplifies the presentation of the existence and uniqueness proof on the fully discrete level, see Lemma 4.3 below. In addition, as stated the scheme is consistent with the closed curve variant from [10].
We have the following discrete analogue of (2.12).
, be a solution to (Q h ) . Then the solution satisfies Proof. Similarly to the proof of [10, Theorem 4.4], on choosing χ = X h t ∈ X h in (3.10a) we can show that (3.12) Moreover, differentiating (3.10c) with respect to time, noting that X h t = 0 on ∂ C I, and then choosing η = Y h g , yields that Finally, choosing χ = (κ h g ) t in (3.10b), and combining with (3.12) and (3.13), yields the desired result (3.11).
The identity (3.11), after integration in time, yields a stability bound for the discrete elastic energy. For the implementation of the presented schemes some metric-dependent quantities need to be calculated. For the metrics in (1.4) and (1.5a), (1.5b) we list these expressions for the convenience of the reader in Table 3. For the metric (1.5c) all the necessary quantities can be calculated with the help of the chain rule, using the fact that e.g. (∇ g)( z) = U T (∇ u Ψ)(u 0 + U z).

Curvature flow
We consider the following linear fully discrete analogue of (A h ) h .  Observe that we have chosen the explicit and implicit terms in (4.1) such that we obtain a linear scheme for which existence and uniqueness can be shown. The basic structure of the scheme goes back to [4, (2.3)]. In fact, for a closed curve in the Euclidean plane the scheme (A m ) h collapses to that approximation of curvature flow.
We make the following mild assumption.

Lemma. 4.1. Let the assumption (A) h hold. Then there exists a unique solution
Proof. As (4.1a), (4.1b) is linear, recall (3.6), existence follows from uniqueness. To investigate the latter, we consider the system: where we recall from (3.6) that λ ∈ V h with λ > 0 in I. It immediately follows from (4.2a) that κ = 0 on ∂ 0 I, and so κ ∈ W h ∂ 0 . Choosing η = δ X ∈ X h in (4.2b) and It follows from (4.3) that κ = κ = 0 and that δ X is constant. Hence (4.2a) and (3.1) imply that It follows from (4.4) and assumption (A) h that δ X = 0. Hence we have shown that (A m ) h has a unique solution ( In order to present an unconditionally stable fully discrete approximation of (C h ) h , we assume that we can split g It follows from the splitting in (4.5) that Then we introduce the following nonlinear scheme.
Note that the scheme (C m, ) h is a natural generalisation of the scheme [12, (C m, ) h ] to the case of open curves.
We can prove the following fully discrete analogue of Theorem 3.1.
where we have used (4.6) and the inequality a .
Splittings of the form (4.5) for the metrics (1.4) have been derived in [12], and we repeat them for the benefit of the reader in Table 3. In the same table we also list, where possible, such splittings for the metrics (1.5). In particular, for the metric (1.5b) we note that D 2 g 1 2 e b z . e 1 e 1 ⊗ e 1 is clearly positive semidefinite, and so we can choose g 1 2 + = g 1 2 and g 1 2 − = 0. Moreover, we now demonstrate how to obtain a splitting of the form (4.5) for the metric (1.5a) in the case n = 2. We leave the case n ≥ 3 to the reader. If n = 2, then we note that

Elastic flow
We consider the following linear fully discrete analogue of the scheme (Q h ) , (3.10).
Note that the scheme (Q m ) is a natural generalisation of the scheme [10, (Q m ) ] to the case of open curves. Observe that the second term in (3.10a) is approximated by the last two terms on the left hand side of (4.9a). This is done in order to allow for an existence and uniqueness proof, see Lemma 4.3 below. In particular, the spatial differential operators acting on Y m+1 g and X m+1 in (4.9a) and (4.9c), respectively, are now the same. This technique is in line with the authors' earlier work in e.g. [8,7,10,14].
We make the following mild assumptions.
(B) Let | X m ρ | > 0 for almost all ρ ∈ I, and let dim span In the case (·, ·) = (·, ·) h the above assumption collapses to (A) h . When dealing with clamped boundary conditions, we also need the following assumption, which is similar to [14,Assumption 5.9].
(C) If Z ∈ Y h with ( Z s , χ s | X m ρ | g ) = 0 for all χ ∈ X h and (g , if the quadrature rule (3.2) has at least one interior sampling point, α k ∈ (0, 1). If (·, ·) = (·, ·) h , on the other hand, then there exists a solution that can be made unique by requiring that κ m+1 Proof. We first consider the case that (3.2) is such that α k ∈ (0, 1) for some 1 ≤ k ≤ K. As (4.9) is linear, existence follows from uniqueness. To investigate the latter, we consider the system: First of all it follows from (4.11), our assumption on (3.2) and the positivities of g( X m ) and | X m ρ | g on I \ ∂ 0 I, that κ g = 0 ∈ V h . As a consequence, we obtain by choosing η = δ X ∈ X h ⊂ Y h in (4.10c) that δ X is a constant vector. Now (4.11) implies that this constant is such that δ X . ω m (q j ) = 0 for all q j ∈ I \ ∂ 0 I, and so the assumption (A) h yields that δ X = 0.
It remains to show that Y g = 0. If ∂ C I = ∅, then we can choose χ = Y g ∈ Y h ⊂ X h in (4.10a) to obtain that Y g is constant in I. Combining (4.10b) with assumption (B) then gives that Y g = 0. If ∂ C I = ∅, on the other hand, then assumption (C) directly gives that Y g = 0, in view of (4.10a) and (4.10b). Hence there exists a unique solution ( X m+1 , κ m+1 For the case (·, ·) = (·, ·) h we note that V h in (4.9b) can be equivalently replaced by to this new system can then be shown similarly to the above proof, which gives all the remaining desired results.
Remark. 4.4. We note that in the examples (1.4d), (1.4e) and (1.5b), any closed curve x(I) in H will correspond to a curve Φ( x(I)) on the hypersurface M that is homotopic to a point. In order to model other curves, the domain H needs to be embedded in an algebraic structure different to R 2 . In particular, H = R × R/(2 π Z) for (1.4d) and (1.5b), and H = R/(2 π s Z) × R/(2 π Z) for (1.4e), respectively.
For the implementation of the presented schemes, this only affects the calculation of differences of vectors in H. For example, for each interval X m (I j ) some care needs to be taken when selecting representatives of the endpoints for the computation of X m ρ , which then naturally yields | X m ρ | and ν m . We will present some numerical simulations for closed curves that are not homotopic to a point in Section 5.

Numerical results
We used the finite element toolbox Alberta, [59], to implement our schemes. The arising linear systems are solved with the sparse factorisation package UMFPACK, see [31]. Solutions to the nonlinear equations for the scheme (C m, ) h are computed with a Newton iteration. The two schemes (A m ) h and (C m, ) h for curvature flow produce very similar results, and can be used interchangeably. We include numerical results for both, in order to demonstrate that they work well in practice. However, for evolutions where numerical stability is crucial, we often prefer to employ the unconditionally stable scheme (C m, ) h .
We note from (3.8) and (2.8) that L h g ( X m ) acts as a discrete energy for (A m ) h and (C m, ) h , while on recalling Theorem 3.2 we define W m+1 g = 1 2 ((κ m+1 g ) 2 , | X m ρ | g ) as a discrete analogue of (2.10) for the scheme (Q m ) . As the quadrature rule for the scheme (Q m ) we either consider (3.1), leading to (Q m ) h , or a quadrature that is exact for polynomials of degree up to five, denoted by (·, ·) , and so leading to (Q m ) . In order to avoid the non-uniqueness issue in Lemma 4.3, we always use the latter in the case of open curves.
The initial data for the scheme (Q m ) , given Γ 0 = X 0 (I), is defined as follows. First Then let κ 0 g ∈ W h ∂ 0 with κ 0 g (q j ) = K(κ 0 , ω 0 , X 0 )(q j ) for q j ∈ I \ ∂ 0 I, recall (3.6). In addition, let Y 0 In most of the presented simulations we use uniform time steps, ∆t m = ∆t, m = 0, . . . , M − 1. For some simulations, however, we use an adaptive time step strategy satisfying ∆t min ≤ ∆t m ≤ ∆t max , m = 0, . . . , M − 1, with smaller time steps at the beginning of the evolution. Unless otherwise stated, in all the simulations we use the discretisation parameters J = 256 and uniform time steps ∆t = 10 −4 .

The metric (1.4a)
For the scheme (A m ) h we show the evolution of two cigar shapes in Figure 1 for the metric (1.4a) with µ = 1. We note that in both cases the curve shrinks to a point. Repeating the same evolutions for the metric (1.4a) with µ = −1, now using the scheme (C m, ) h , leads to the results shown in Figure 2. While the horizontally aligned curve again shrink to a point, the vertically aligned curve approaches the x 2 -axis in order to minimise its geodesic length. The degeneracy of g on the axis leads to a breakdown of the evolution. In practice this means that the Newton iteration to find a solution for (C m, ) h no longer converges. Here we note that we used the smaller uniform time step size ∆t = 10 −5 for this experiment.
We stress that the evolution is well defined, however, if we assign boundary points to lie on the x 2 -axis and to move freely on it. This is not dissimilar to the modelling of mean curvature flow for axisymmetric surfaces of genus zero, see [11] for details. As an example, we show the evolution of a semicircle with radius 1 and ∂ 0 I = ∂I in Figure 3. As a comparison, we also show the same evolution for the case ∂ 1 I = ∂I. In both cases, the semicircle shrinks to extinction, but the shape and time scale of the two evolutions differ.
For completeness, we also show some evolutions for the cases ∂ D I = ∂I and ∂ 2 I = ∂I in Figure 4. The first evolution for the Dirichlet, or no-slip, boundary conditions leads to the curve trying to reach the x 2 -axis in order to reduce its length. Similarly to Figure 2 this leads to a breakdown of the scheme. The second evolution for the Dirichlet conditions yields a straight line segment as geodesic, while for the free-slip condition the initial semicircle shrinks to a point on the x 1 -axis.
Evolutions for elastic flow with Navier and clamped boundary conditions, respectively, are shown in Figure 5. Here, for the clamped boundary conditions, recall (2.27), we choose ζ(p) = (sin ϑ(p), cos ϑ(p)) T , with ϑ(0) = 210 • and ϑ(1) = 150 • . While in the Navier case the curve appears to grow unboundedly, in the clamped case the curve seems to approach an optimal shape aligned with the chosen metric.
We remind the reader that many more numerical simulations for closed curves moving under curvature flow or elastic flow in the Riemannian manifold defined by (1.4a),     The solutions X m at times t = 0, 2, . . . , 6. Below we visualise Φ( X m ) at times t = 0 (blue), t = 2 (red) and t = 6 (black), for (1.4e) with s = 1, and also show a plot of the discrete energy L h g ( X m ).
including for the case case µ = 1 for the hyperbolic plane, can be found in [12,10].

The torus metric (1.4e)
A geodesic between two fixed points on the Clifford torus is computed in Figure 6. To this end, we employ the metric induced by (1.4e) with s = 1, so that the torus has radii r = 1 and R = √ 2. We observe that the evolution eventually settles on a geodesic, that is clearly not the shortest path connected the two points on the torus. That is because of a topological restriction stemming from the fact that the curve must stay within the equivalence class that is prescribed by the initial data.
On recalling Remark 4.4, we also present an evolution for geodesic curvature flow of a closed curve that is not homotopic to a point. See Figure 7 for a presentation of the numerical results for the scheme (C m, ) h .

The Angenent metric (1.5a)
Unless otherwise stated, we choose n = 2 in (1.5a). First we show the evolution under curvature flow of an eongated cigar shape that shrinks to a point, see in Figure 8.
In a second experiment, we show the evolution under elastic flow of a circle towards   the generating curve of the Angenent torus in an axisymmetric setting. We recall that the Angenent torus is a critical point of Huisken's F -functional (1.6), and hence a self-shrinker for mean curvature flow in R 3 , with extinction time 1. As a consequence, the generating curve of the Angenent torus, which from now on we will also simply call Angenent torus, is a critical point of the geodesic length L g , and hence a geodesic. For the evolution shown in Figure 9, we observe that the discrete curvature energy W m+1 g reduces from about 3.5 to about 10 −5 , giving a strong indication that we have indeed found a geodesic. Note also that the final shape in Figure 9 agrees well with the numerical results in [22,17,3]. We have also performed simulations for elastic flow of initial curves with a winding number larger than one, with respect to the point 2 e 1 , and they always settle as a stationary solution on a multiple covering of the Angenent torus.
It is known that the Angenent torus is an unstable critical point of the geodesic length L g , see [23,54,18], and this is confirmed by our numerical experiments. Hence it is practically impossible to obtain an approximation to it as a long-time limit of curvature flow. We demonstrate this phenomenon by starting two simulations for the stable scheme (C m, ) h from slightly shifted Angenent tori. Our numerical results in Figure 10 confirm that the stationary solution is unstable, and we see the curve either moving monotonically towards the x 2 -axis, or towards infinity, with a significant decrease in the geodesic length of the curve in each case. For these experiments we used the finer discretisation parameters J = 2048 and ∆t = 10 −5 .
We highlight the capabilities of our numerical method by computing the "Angenent tori" in dimensions four and five, that is hypersurfaces in R n+1 that are topologically equivalent to S 1 × S n−1 , n = 3, 4, and that are self-shrinkers for mean curvature flow with extinction time 1. In particular, in Figure 11 we show the numerical steady states for approximations of elastic flow for the metric (1.5a), with n = 2, 3, 4. In each case the final discrete energy satisfies | W M g | < 10 −9 , confirming that we are indeed approximating geodesics.
Hence in order to deduce the value of Huisken's F -functional for the self-shrinkers that we have found computationally, it is sufficient to report on the length of the corresponding geodesics, which we will do from now on within the captions of the relevant figures. Here we note that the pre-factor in (5.1) for the cases n = 2, 3, 4 is given by 2 1−n [Γ( n 2 )] −1 = 1 2 , 1 2 √ π and 1 8 , respectively. For the three geodesics displayed in Figure 11, our computed values for the final discrete length L h g ( X M ) are 3.70, 6.39 and 14.27, giving approximate values for the entropy of the associated surfaces of revolution of 1.85, 1.80 and 1.78, respectively. We remark that the value for n = 2 agrees very well with the results reported in [17,3].
Denoting the entropy of the n-dimensional "Angenent torus" with E n , we have so far established that E 2 ≈ 1.85, E 3 ≈ 1.80 and E 4 ≈ 1.78. Continuing the above procedure for increasing values of n, we numerically obtain that E 6 ≈ 1.77, E 8 ≈ 1.76 and E 12 ≈ 1.75. This leads us to conjecture that E n is a strictly monotonically decreasing sequence with E n √ 3 as n → ∞. Note that this conjecture is in the spirit of Lemma A.4 in [61].
Of course, the most famous self-shrinker for mean curvature flow in R n+1 , with extinction time 1, is the sphere of radius √ 2 n, see e.g. [23]. In the context of the metric (1.5a), these correspond to geodesics in the shape of semicircles with radius √ 2 n. For n = 2 and n = 3 we show an evolution each for elastic flow towards these geodesics, see Figure 12 for details, where in each case as initial data we choose a semicircle of radius n − 1. For the two geodesics displayed in Figure 12 In order to investigate the numerical accuracy of our proposed method, in Table 4 we compare the numerical steady state solutions of elastic flow in the case n = 2 with a semicircle of radius 2, as well as their respective induced entropy values. That is, we list the errors E M = max j=0,...,J |2−| X M (q j )|| and e M = | 1 2 L h g ( X M )− 4 e | for increasing values of J. Here we fix T = 10 and ∆t = 10 −5 , and always start from the approximation of a unit semicircle. The results presented in Table 4 appear to show a convergence rate of 1.5 for the discrete L ∞ error E M , while the entropy error e M appears to converge quadratically.
The final simulations in this subsection are devoted to finding self-shrinkers for mean curvature flow that are non-embedded, inspired by the work [34]. We begin with an experiment for a closed curve with seven self-intersections, see Figure 13. Under elastic flow the curve evolves towards the generating curve of a non-embedded self-shrinker for mean curvature flow. In fact, the steady state corresponds to the shape in [34, Figure 6]. Due to the large energy decrease at the beginning of the evolution, we use an adaptive   Table 4: Errors and experimental orders of convergence for the convergence test corresponding to the final shape of the evolution on the left of Figure 12.
time stepping strategy with ∆t min = 10 −7 and ∆t max = 10 −6 . The spatial discretisation uses J = 512. The discrete energy of the final solution satisfies | W M g | < 10 −9 , confirming that we are indeed approximating a geodesic. We note that the discrete geodesic length of the final curve is 10.53, giving an entropy of 5.26.
We also investigate, what happens to the geodesic from Figure 13 if we change the metric to (1.5a) with n = 3. See Figure 14 for a plot of the obtained numerical result, which compared to the geodesic for n = 2 has shifted further to the right. For this experiment we once again used an adaptive time stepping strategy. We note that the discrete energy of the final solution satisfies | W M g | < 10 −8 . The discrete geodesic length of the final curve is 18.24, giving an entropy of 5.15.
Inspired by [34, Figure 3], we now perform a numerical simulation to find a nonembedded self-shrinker of genus zero for mean curvature flow. Starting from an initial curve with three self-intersections, we observe the evolution for elastic flow shown in Figure 15, where the final discrete energy satisfies | W M g | < 10 −9 . We note the excellent agreement with [34, Figure 3]. Here we again make use of an adaptive time stepping strategy with ∆t min = 10 −7 and ∆t max = 10 −4 . The spatial discretisation uses J = 512.  We note that the discrete geodesic length of the final curve is 7.33, giving an entropy of 3.66.

The cone metric (1.5b)
In a first experiment for the metric (1.5b), we look at (geodesic) curvature flow for a curve on a cone with b = 0.5, and so β = 3 − 1 2 in (1.7). For the simulation in Figure 16 it can be observed that in H the initial circle of radius 2 deforms and shrinks to a point. On the hypersurface M = Φ(H), the initial curve is homotopic to a point, and so shrinks to a point away from the apex.
The following conjecture on geodesic curvature flow on a cone was formulated by Charles M. Elliott, [36].
Conjecture. 5.1. A closed curve on a cone M, that is not homotopic to a point on M, will under geodesic curvature flow converge to the apex in finite time.  (A m ) h Geodesic curvature flow on a cone. The solutions X m at times t = 0, 0.5, 1. On the right we visualise Φ( X m ) at times t = 0, 0.5, 1, for (1.5b) with b = 0.5. A plot of the discrete energy L h g ( X m ) is shown on the right.
The conjecture says, in particular, that the singularity of the flow will happen at the apex. Indeed we expect that all the points of the curve converge to the apex at the singular time. Moreover, we expect that a similar conjecture holds on more general surfaces on which a curve encloses a singularity.
On recalling Remark 4.4, we now numerically test the conjecture by starting an evolution for geodesic curvature flow with a closed curve that is very close to the apex, but not uniformly so. That is, we vary the x 2 -coordinate of the initial curve in H between ±2. During the evolution, the parts of the curve closest to the apex first start to rise, making the curve becoming more circle-like, before the whole curve sinks towards to apex. See Figure 17, where we also show a plot of the lowest point of the curve on the cone over time, highlighting the rise and fall of the curve on the cone. The observed behaviour is consistent with Conjecture 5.1  An experiment for (geodesic) elastic flow on the same cone is shown in Figure 18. Here the closed curve first approaches a circle, which then rises along the cone. By computing the energy one observes that a circle with increasing radius reduces the elastic energy.

The metric (1.5c)
We end the section on the numerical results for our presented schemes with some simulations for the metric (1.5c) with (1.8) and (1.9). Recall that now geodesics in H correspond to optimal interface profiles in multi-component phase field models. Of particular interest are geodesics, or shortest paths, that connect the vertices e 1 , e 2 , e 3 of the Gibbs simplex G, recall (1.10). To this end, we note that with the choice (1.8), it holds that the map z → f ( z) = u 0 + U z satisfies f (0, 0) = e 1 , f (−2 the first experiment we set (σ 12 , σ 13 , σ 23 , σ 123 ) = (4, 6, 9, 0), and numerically compute a geodesic connecting e 1 and e 2 with the help of geodesic curvature flow. Here we always use the scheme (A m ) h with the uniform time step size ∆t = 10 −5 . The results are shown in Figure 19, where we see that the flow quickly settles on a curved geodesic. We repeat the same simulation also for the paths connecting the pure phases e 1 and e 3 , as well as e 2 and e 3 , and plot all three solutions within the Gibbs simplex G, recall (1.10), at the bottom of Figure 19. In [41] numerical computations indicated that on choosing σ 123 > 0 in (1.9) larger and larger, the minimising profiles can be forced to approach the edges of the Gibbs simplex. To confirm this effect with our numerical method, we now choose σ 123 ∈ {10, 100, 1000} and plot the obtained geodesics in Figure 20. It can be seen that for an increasing value of σ 123 , the geodesics are pushed further and further towards the edges of the Gibbs simplex.
We remark that in [19] a novel approach for multi-component phase field models has been considered, where the Ginzburg-Landau energy can be defined such that the minimising paths connecting the pure phases are always given by the edges of the Gibbs simplex. The metric that would arise in the form of (1.5c) in order to model this situation is in general no longer conformal, and so would be outside the context of this paper. However, following the approach in [42,38], we can consider the following replacement of  (1.9) to achieve the same effect: Ψ(u 1 , u 2 , u 3 ) = σ 12 u 2 1 u 2 2 + σ 13 u 2 1 u 2 3 + σ 23 u 2 2 u 2 3 + σ 123 u 1 u 2 u 2 3 + σ 231 u 2 u 3 u 2 1 + σ 312 u 3 u 1 u 2 2 , (5.2) where σ 123 , σ 231 , σ 312 ∈ R ≥0 . We perform a computation for (5.2) with σ 12 = σ 13 = σ 23 = σ 123 = σ 231 = σ 312 = 1 and show the obtained results in Figure 21. It can be seen that now the geodesics lie on the edges of the Gibbs simplex, confirming the analysis in [42,38].