Abstract
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.
Similar content being viewed by others
1 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\subset {\mathbb {R}}^2\) we let the metric tensor be given by
where \({\vec {v}}\,.\,{\vec {w}} = {\vec {v}}^T\,{\vec {w}}\) is the standard Euclidean inner product, and where \(g:H \rightarrow {{\mathbb {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, 50, 58]. Examples of such situations are the hyperbolic plane, the hyperbolic disc and the elliptic plane. Other examples are given by two-dimensional manifolds in \({\mathbb {R}}^d\), \(d\ge 3\), that can be conformally parameterised, such as spheres without pole(s), catenoids and tori. Coordinates \((x_1,x_2)\in 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 Sect. 2 and [50, 3.29 in Sect. 3D] for more information.
The metric tensor (1.1) induces a notion of length in H. In particular, the length of a vector \({\vec {v}} \in {{\mathbb {R}}}^2\) at the location \({\vec {z}} \in H\) is defined by
whereas the length of a curve \(\vec {\gamma } \in 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 \({\vec {z}}_0, {\vec {z}}_1 \in H\) defined as
It can be shown that \((H,{\text {dist}}_g)\) is a metric space, see [46, Sect. 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 \((\Gamma (t))_{t\in [0,T]}\), can be defined as the \(L^2\)–gradient flow of geodesic length, \(L_g(\Gamma ) = \int _\Gamma g^\frac{1}{2} \;\mathrm{d}s\), with respect to the \(L^2\)–inner product \(\langle v,w\rangle = \int _\Gamma v\,w\,g^\frac{1}{2} \;\mathrm{d}s\). It is often called curve shortening flow. In particular, on letting the geodesic curvature \(\varkappa _g\) be the first variation of \(L_g(\Gamma )\), with respect to \(\langle \cdot ,\cdot \rangle \), then curvature flow is given by
where \({\mathcal {V}}_g = g^{\frac{1}{2}}\,{\mathcal {V}}\), and \({\mathcal {V}}\) is the Euclidean normal velocity of \(\Gamma \) in \({{\mathbb {R}}}^2\). Moreover, the geodesic elastic energy of \(\Gamma \) is defined by \(W_g(\Gamma ) = \tfrac{1}{2}\,\int _\Gamma (\varkappa _g)^2\,g^\frac{1}{2} \;\mathrm{d}s\), and elastic flow arises as the \(L^2\)–gradient flow of \(W_g(\Gamma )\). The elastic energy is often called bending energy of \(\Gamma \), and elastic flow is also known by the name curve straightening flow. Critical points of \(W_g(\Gamma )\) 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(\Gamma )\) that are unstable under curvature flow. In fact, all geodesics satisfy \(\kappa _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. [1, 20, 37], while numerical approximations have been considered in [6, 12, 21, 55, 56, 60], 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. [25,26,27,28, 39, 40, 43, 44]. As far as the numerical approximation of elastic flow is concerned, we remark that the case of a Euclidean ambient space has been treated in [7, 15, 32]. 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, 11, 12]. 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 [11, 12]. 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 Sect. 2.1 for curvature flow, see (2.22), and in Sect. 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 \(\overline{{{\mathbb {H}}}^2} = \{ {\vec {z}} \in {{\mathbb {R}}}^2 : {\vec {z}} \,.\,\vec {e}_1 \ge 0 \}\). Metrics that the authors have considered in their recent works on closed curves, see [11, 12], are
Recall that (1.4a) with \(\mu =1\) models the hyperbolic plane, while \(\mu =0\) corresponds to the Euclidean plane. The metric (1.4b) for \(\alpha =1\) models the hyperbolic disk, while \(\alpha =-1\) yields the elliptic plane. Moreover, (1.4c), (1.4d) and (1.4e) arise from conformal parameterisations of spheres, catenoids and tori, respectively, where in the latter case the torus has large radius \([1 + {\mathfrak {s}}^2]^\frac{1}{2}\) and small radius 1.
Additional metrics that we consider in this paper are
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 \({\mathcal {S}} \subset {{\mathbb {R}}}^{n+1}\), with generating curve \(\Gamma \subset {{\mathbb {H}}}^2\), the geodesic length of \(\Gamma \), with respect to the metric (1.5a), collapses, up to a constant factor, to Huisken’s F-functional
see [17, 23, 45]. 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].
The metric (1.5b), on the other hand, arises from a conformal parameterisation of a right circular cone without the apex as follows. Let
We recall from [12, §2.4] that \(\vec \Phi : H \rightarrow {\mathcal {M}}\) is a conformal parameterisation of \({\mathcal {M}}\), if \({\mathcal {M}} = \vec \Phi (H)\), \(|\partial _{\vec {e}_1} \vec \Phi ({\vec {z}})|^2 = |\partial _{\vec {e}_2} \vec \Phi ({\vec {z}})|^2\) and \(\partial _{\vec {e}_1} \vec \Phi ({\vec {z}}) \,.\, \partial _{\vec {e}_2} \vec \Phi ({\vec {z}}) = 0\) for all \({\vec {z}} \in H\). Using the ansatz
for some function \(r\in C^\infty ({{\mathbb {R}}},{{\mathbb {R}}}_{>0})\), it is easy to see that \(\partial _{\vec {e}_1}\vec \Phi \,.\, \partial _{\vec {e}_2}\vec \Phi = 0\), as well as
which shows that \(r(u) = e^{{\mathfrak {b}}\,u}\), where \({\mathfrak {b}} = [\frac{\beta ^2}{1+\beta ^2}]^{\frac{1}{2}}\), leads to a conformal parameterisation of \({\mathcal {M}}\). In particular, setting \(g = |\partial _{\vec {e}_1} \vec \Phi |^2 = |\partial _{\vec {e}_2} \vec \Phi |^2\) leads to (1.5b).
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 \(\sigma _{12}\,,\sigma _{13}\,,\sigma _{23}\in {{\mathbb {R}}}_{>0}\) and \(\sigma _{123} \in {{\mathbb {R}}}_{\ge 0}\), corresponds to [41, (24), (25)], where \(\mathbf{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
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, \(\mathbf{e}_1= (1,0,0)^T\), \(\mathbf{e}_2=(0,1,0)^T\) and \(\mathbf{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 [19, 38, 41] for more details.
The remainder of this paper is organised as follows. In Sect. 2 we present strong and weak formulations of curvature flow and elastic flow. The semidiscrete continuous-in-time finite element approximations introduced in Sect. 3 will be based on these weak formulations. Stability of some of the schemes is also shown in Sect. 3. Fully discrete approximations are presented in Sect. 4, together with results on their well-posedness and stability, where applicable. Finally, in Sect. 5 we present several numerical simulations for the derived schemes and the various metrics considered in this paper.
2 Mathematical formulations
We let \({{\mathbb {R}}} / {{\mathbb {Z}}}\) be the periodic interval [0, 1], and set
Consider a family of curves \((\Gamma (t))_{t\in [0,T]}\), \(T > 0\), that can be either open, \(I=(0,1)\), or closed, \(I={{\mathbb {R}}} / {{\mathbb {Z}}}\). Given some smooth parameterisation \({\vec {x}}: I\times [0,T] \ni (\rho ,t) \mapsto {\vec {x}}(\rho ,t)\in {{\mathbb {R}}}^2\), with \(|{\vec {x}}_\rho | > 0\) in \({\overline{I}} \times [0,T]\), we introduce the arclength s of the curve, i.e. \(\partial _s = |\vec {x}_\rho |^{-1}\,\partial _\rho \), and set
where \(\cdot ^\perp \) denotes a clockwise rotation by \(\frac{\pi }{2}\). We let \({\mathcal {V}} = {\vec {x}}_t\,.\,\vec \nu \) denote the normal velocity, and let the Euclidean curvature \(\varkappa \) be defined by
see [33]. We also let
We introduce
so that \(\vec \tau _g\,.\,\vec \nu _g = 0\) and \(|\vec \tau _g|_g^2 = |\vec \nu _g|_g^2 = (\vec \nu _g, \vec \nu _g)_g = g({\vec {x}})\,\vec \nu _g\,.\,\vec \nu _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 \(\cdot _g\) above, and throughout the paper, denote quantities associated with the metric g. The subscripts \(\cdot _t\) and \(\cdot _\rho \), on the other hand, denote partial derivatives with respect to t and \(\rho \), respectively. Finally, \(\cdot _s\) and \(\cdot _{s_g}\) 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 \({\mathcal {V}}_g\) and \(\varkappa _g\), similarly to \({\mathcal {V}}\) and \(\varkappa \), only depend on \(\Gamma \), but not on the chosen parameterisation \({\vec {x}}\). For the case of an evolving closed curve, \(\partial I = \emptyset \), we recall from [12] that curvature flow,
is the \(L^2\)–gradient flow of the geodesic length,
recall (1.3). In particular, for closed curves evolving by (2.7) it holds that
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 \(\mu =0\), the geodesic elastic flow (2.11) collapses to classical elastic flow \({\mathcal {V}} = -\varkappa _{ss} - \tfrac{1}{2}\,\varkappa ^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 \(\partial I\):
Clearly, (2.13)(i) means that we consider the endpoint fixed in time, while in (2.13)(ii) and (2.13)(iii) we allow the boundary point to move freely parallel to the \(x_2\)– and \(x_1\)–axis, respectively. For some of the metrics in (1.4) and (1.5) it is possible to \(C^1\)–continuously extend the metric g to \(\overline{{{\mathbb {H}}}^2}\) such that \(g=0\) on the \(x_2\)–axis. In fact, this holds precisely for (1.4a) with \(\mu \le -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 \(\partial I = \partial _0 I \cup \partial _{1} I \cup \partial _{2} I\cup \partial _C I \cup \partial _D I \cup \partial _N I\) with the conditions
In the above \(\partial _0 I\) denotes the subset of boundary points of I that correspond to endpoints of \(\Gamma (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({\vec {x}}(0)) = 0\) on \(\partial _0 I\), so that (2.14a) implies \(g({\vec {x}}(t)) = 0\) on \(\partial _0 I\) for all t. In our paper, we will consider (2.14a) only for (1.4a), with \(\mu \le -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 \(\varkappa _g\), recall (2.6). To this end, we note that it can be easily shown from (2.2) that
see [12, (2.16)]. Combining (2.6) and (2.15) yields that
Let \((\cdot ,\cdot )\) denote the \(L^2\)–inner product on I, and let
We also define
2.1 Curvature flow
For curvature flow we assume that \(\partial I = \partial _0 I \cup \partial _1 I \cup \partial _2 I\cup \partial _C I\), where \(\partial _0 I\) will only be nonempty for the metrics (1.4a), with \(\mu \le -1\), and (1.5a).
It holds that
Combining (2.18), (2.15), (2.4), (2.5) and (1.2) yields that
where we have recalled (2.6). Clearly, the curvature \(\varkappa _g\) is the first variation of the length (2.8).
For the case that \(\partial _0 I \not = \emptyset \), we note that in order for the right hand side of (2.19) to remain bounded, it is appropriate to require that the term \(\vec \nu \,.\,\nabla \,g^\frac{1}{2}({\vec {x}})\) in the second line remains bounded as we approach \(\partial _0 I\). In view of
with \(\nabla \,g({\vec {x}})\,.\,\vec {e}_2 = 0\) on \(\partial _0 I\), since \(g=0\) on \(\partial _0 I\), we see that
is a natural assumption to make. Moreover, in situations where the curve \(\Gamma (t) = {\vec {x}}({\overline{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 [9, 10] 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 \(\partial _0 I\) and \(\partial _D I\), while on the remainder of \(\partial I\) we require
In terms of the Euclidean properties of the curve \(\Gamma (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 \(\partial _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 \({\vec {x}}_\rho \,.\,\vec {e}_{3-i}=0\) in (2.22d) corresponds to a \(90^\circ \) 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).
\(({\mathcal {A}})\): Let \({\vec {x}}(0) \in \underline{V}_{\partial _0}\). For \(t \in (0,T]\) find \({\vec {x}}(t) \in [H^1(I)]^2\), with \({\vec {x}}_t(t) \in {\mathbb {X}}\), and \(\varkappa (t)\in L^2(I)\) such that
We note that the boundary conditions for \({\vec {x}}_t\) in (2.22) are enforced through the trial space \({\mathbb {X}}\), recall (2.17), while the boundary conditions on \({\vec {x}}_\rho \) in (2.22) follow from (2.23b).
An alternative weak formulation, based directly on (2.22), together with (2.16), is given as follows.
\(({\mathcal {C}})\): Let \({\vec {x}}(0) \in \underline{V}_{\partial _0}\). For \(t \in (0,T]\) find \({\vec {x}}(t) \in [H^1(I)]^2\), with \({\vec {x}}_t(t) \in {\mathbb {X}}\), and \(\varkappa _g(t) \in L^2(I)\) such that
Once again, the boundary conditions for \({\vec {x}}_t\) in (2.22) are enforced through the trial space \({\mathbb {X}}\), while the conditions on \({\vec {x}}_\rho \) in (2.22d) follow directly from (2.24b). In addition, it can be shown that (2.24b), for the metrics (1.4a), with \(\mu \le -1\), and (1.5a), also enforces the condition on \({\vec {x}}_\rho \) in (2.22b). In fact, this follows by using the techniques in [10, Appendix A], and noting that for both metrics it holds that \(\nabla \,g^\frac{1}{2}({\vec {x}})\,.\,\vec {e}_2 = 0\) on \(\partial _0 I\), see Table 3 below.
2.2 Elastic flow
For elastic flow we assume that \(\partial I = \partial _0 I \cup \partial _C I \cup \partial _N I\), where, as before, \(\partial _0 I\) will only be nonempty for the metrics (1.4a), with \(\mu \le -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 \(\partial _C I \cup \partial _N I\), the first two terms vanish. On \(\partial _N I\) we require \(\varkappa _g=0\) to make the third term zero, while the clamped boundary conditions
with \(\vec \zeta : \partial _C I \rightarrow {{\mathbb {S}}}^1 = \{{\vec {z}} : {{\mathbb {R}}}^2 : |{\vec {z}}| = 1\}\), ensure as usual that \({\mathcal {V}}_s = 0\) on \(\partial _C I\); see e.g. Lemma 37(ii) in [13]. For the compact notation in (2.27) we have used the fact that the identity function \(\mathrm{id}\) always maps elements of \(\partial I\) to either 0 or 1. On \(\partial _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 \({\mathcal {V}}_s = 0\) on \(\partial _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
or
These conditions will also lead to vanishing boundary terms in (2.26). Observe that (2.30) involves a \(90^\circ \) 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 \([\varkappa _g]_\rho -\tfrac{1}{2} \varkappa _g[\ln g({\vec {x}})]_\rho = 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 [7, 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 [11, 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 \(({\mathcal {P}})\) and \(({\mathcal {Q}})\) in [11], while the treatment of the boundary conditions, achieved through the spaces \({\mathbb {X}}\) and \({\mathbb {Y}}\) from (2.17), is very similar to the approach taken in [14] in the context of axisymmetric Willmore flow.
\(({\mathcal {P}})\): Let \({\vec {x}}(0) \in \underline{V}_{\partial _0}\). For \(t \in (0,T]\) find \({\vec {x}}(t) \in [H^1(I)]^2\), with \({\vec {x}}_t(t) \in {\mathbb {X}}\), \({\vec {y}}(t) \in {\mathbb {Y}}\) and \(\varkappa \in L^2(I)\) such that
The consistency of (2.32), in the case \(I = {{\mathbb {R}}} / {{\mathbb {Z}}}\), was shown in [11, Appendix A.1]. As far as the boundary conditions are concerned, we note that the conditions on \({\vec {x}}_t\) in (2.29) are enforced through the trial space \({\mathbb {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 \({\vec {y}}\,.\,\vec \nu = \varkappa _g\), and so the trial space \({\mathbb {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 [11, (A.3)–(A.9)], noting in particular that the integration by parts in [11, (A.5)] gives rise to the boundary term \((\varkappa _g)_s - \tfrac{1}{2}\,\varkappa _g\,(\ln g({\vec {x}}))_s\) on \(\partial _0 I\), which enforces (2.28), as required.
\(({\mathcal {Q}})\): Let \({\vec {x}}(0) \in \underline{V}_{\partial _0}\). For \(t \in (0,T]\) find \({\vec {x}}(t) \in [H^1(I)]^2\), with \({\vec {x}}_t(t) \in {\mathbb {X}}\), \({\vec {y}}_g(t) \in {\mathbb {Y}}\) and \(\varkappa _g \in L^2(I)\) such that
The consistency of (2.33), in the case \(I = {{\mathbb {R}}} / {{\mathbb {Z}}}\), was shown in [11, Appendix A.2]. We note that the boundary conditions on \({\vec {x}}_t\) in (2.29) are once again enforced through the trial space \({\mathbb {X}}\). The second condition in (2.29c) is enforced through (2.33c), and the same can be shown for the second condition in (2.29b), similarly to (2.24). In addition, (2.33b) together with the trial space \({\mathbb {Y}}\) yields the second condition in (2.29d). It remains to validate that (2.33a) weakly enforces the third condition in (2.29b). As before, this can be done on closely following the argument in [11, (A.10)–(A.19)], and collecting the terms that appear on \(\partial _0 I\) due to integration by parts. In fact, [11, (A.13), (A.14)] yield the boundary term \((\varkappa _g)_s - \tfrac{1}{2}\,\varkappa _g\,(\ln g({\vec {x}}))_s + (g^\frac{1}{2}({\vec {x}})\,\varkappa - g({\vec {x}})\,\varkappa _g) \,{\vec {y}}_g\,.\,\vec \tau \) on \(\partial _0 I\), which thanks to \({\vec {y}}_g \in {\mathbb {Y}}\) and \({\vec {x}}_\rho \,.\,\vec {e}_1 = 0\) collapses to enforcing (2.28).
3 Semidiscrete finite element approximations
Let \([0,1]=\cup _{j=1}^J I_j\), \(J\ge 3\), 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={{\mathbb {R}}} / {{\mathbb {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 \(\underline{V}^h_{\partial _0}= \underline{V}^h\cap \underline{V}_{\partial _0}\) and \(W^h_{\partial _0} = \{ \chi \in V^h : \chi (\rho ) = 0 \quad \forall \ \rho \in \partial _0 I\}\), as well as \({\mathbb {X}}^h = {\mathbb {X}}\cap \underline{V}^h\) and \({\mathbb {Y}}^h = {\mathbb {Y}}\cap \underline{V}^h\), recall (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 \(\{q_j\}_{j=1}^J\), via
where we define \(u(q_j^\pm )=\underset{\delta \searrow 0}{\lim }\ u(q_j\pm \delta )\). The definition (3.1) naturally extends to vector valued functions. Moreover, let \((\cdot ,\cdot )^\diamond \) 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=1}^J\), we let \((u,v)^\diamond = I^\diamond (u\,v)\), where
with \(K\ge 2\), \(\sum _{k=1}^K w_k = 1\), and with distinct \(\alpha _k\), \(k=1,\ldots ,K\). A special case is \((\cdot ,\cdot )^\diamond = (\cdot ,\cdot )^h\), recall (3.1), but we also allow for more accurate quadrature rules.
Let \(({\vec {X}}^h(t))_{t\in [0,T]}\), with \({\vec {X}}^h(t)\in \underline{V}^h\), be an approximation to \(({\vec {x}}(t))_{t\in [0,T]}\). Then, similarly to (2.1), we set
For later use, we let \(\vec \omega ^h \in \underline{V}^h\) be the mass-lumped \(L^2\)–projection of \(\vec \nu ^h\) onto \(\underline{V}^h\), i.e.
3.1 Curvature flow
We consider the following finite element approximation of \(({\mathcal {A}})\), recall (2.23). It is closely related to the approximation [12, (3.3), (3.10)] for closed curve evolutions.
\(({\mathcal {A}}_h)^h\): Let \({\vec {X}}^h(0) \in \underline{V}^h_{\partial _0}\). For \(t \in (0,T]\), find \(({\vec {X}}^h(t), \kappa ^h(t)) \in \underline{V}^h\times V^h\), such that \({\vec {X}}^h_t(t) \in {\mathbb {X}}^h\), and such that
where we have defined \({\mathcal {K}}(\kappa ^h, \vec \omega ^h, {\vec {X}}^h) \in V^h\) by
To motivate the choice (3.6), we observe that it follows from (2.20), Table 3 and L’Hospital’s rule that
as \(\rho \rightarrow \rho _0 \in \partial _0 I\). On recalling (3.4), we note that replacing \(\vec \omega ^h\) with \(\vec \nu ^h\) in (3.5b) does not change the scheme. On the other hand, for nonconstant g we must use \(\vec \omega ^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 \(\vec \omega ^h\) there.
It does not appear possible to prove a stability result for the scheme (3.5). However, thanks to (3.5b) the scheme \(({\mathcal {A}}_h)^h\) satisfies a weak equidistribution property, i.e. it can be shown that neighbouring elements of \(\Gamma ^h(t) = {\vec {X}}^h({\overline{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 [4] 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 \(({\mathcal {C}})\), recall (2.24). It is the natural extension to the open curve case of the semidiscrete analogue of [12, (3.5), (3.18)].
\(({\mathcal {C}}_{h})^{h}\): Let \({\vec {X}}^h(0) \in \underline{V}^h_{\partial _0}\). For \(t = (0,T]\), find \(({\vec {X}}^h(t), \kappa _g^h(t)) \in \underline{V}^h\times W^h_{\partial _0}\), such that \({\vec {X}}^h_t(t) \in {\mathbb {X}}^h\), and such that
Observe that we fix \(\kappa _g^h(t)\) to be zero on \(\partial _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 \(\vec \nu ^h\) with \(\vec \omega ^h\) we obtain a related, but different, approximation of \(({\mathcal {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 \(({\mathcal {A}}_h)^h\), the scheme \(({\mathcal {C}}_{h})^{h}\) also exhibits some implicit tangential motion, but, in contrast to \(({\mathcal {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 \({\vec {Z}} \in \underline{V}^h\) we let
Then we can prove the following discrete analogue of (2.9) for the scheme (3.7).
Theorem 3.1
Let \(({\vec {X}}^h(t),\kappa ^h_g(t)) \in \underline{V}^h\times W^h_{\partial _0}\), for \(t\in (0,T]\), be a solution to \(({\mathcal {C}}_h)^h\). Then the solution satisfies the stability bound
Proof
Choosing \(\chi = \kappa _g^h \in W^h_{\partial _0}\) in (3.7a) and \(\vec \eta = {\vec {X}}^h_t \in {\mathbb {X}}^h\) in (3.7b) yields
and so we obtain the desired result (3.9). \(\square \)
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).
3.2 Elastic flow
Following the approach by the authors in [11], it is straightforward to derive a semidiscrete approximation of \(({\mathcal {P}})\), in the case that \(\partial _0 I = \emptyset \). The derived scheme, which will correspond to [11, \(({\mathcal {P}}_h)^h\)] with the natural changes to the test and trial spaces, can be shown to be stable. Moreover, and similarly to \(({\mathcal {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 \(\partial _0 I \not = \emptyset \), we do not pursue this variant any further in this paper.
On the other hand, the following discretisation based on the formulation \(({\mathcal {Q}})\) can naturally deal with all the considered boundary conditions. It is the natural extension to the open curve case of the semidiscrete scheme [11, \(({\mathcal {Q}}_h)^\diamond \)].
\(({\mathcal {Q}}_h)^\diamond \): Let \({\vec {X}}^h(0) \in \underline{V}^h_{\partial _0}\). For \(t \in (0,T]\), find \(({\vec {X}}^h(t), \kappa ^h_g(t), {\vec {Y}}^h_g(t)) \in \underline{V}^h\times V^h \times {\mathbb {Y}}^h\), with \({\vec {X}}^h_t(t) \in {\mathbb {X}}^h\), such that
If we replace \(\vec \omega ^h\) with \(\vec \nu ^h\) in (3.10a) we obtain a related, but different, approximation of \(({\mathcal {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 [11].
We have the following discrete analogue of (2.12).
Theorem 3.2
Let \(({\vec {X}}^h(t),\kappa ^h_g(t),{\vec {Y}}_g^h(t)) \in \underline{V}^h\times V^h \times {\mathbb {Y}}^h\), for \(t\in (0,T]\), be a solution to \(({\mathcal {Q}}_h)^\diamond \). Then the solution satisfies
Proof
Similarly to the proof of [11, Theorem 4.4], on choosing \(\vec \chi = {\vec {X}}^h_t \in {\mathbb {X}}^h\) in (3.10a) we can show that
Moreover, differentiating (3.10c) with respect to time, noting that \({\vec {X}}^h_t = 0\) on \(\partial _C I\), and then choosing \(\vec \eta = {\vec {Y}}^h_g\), yields that
Finally, choosing \(\chi = (\kappa _g^h)_t\) in (3.10b), and combining with (3.12) and (3.13), yields the desired result (3.11). \(\square \)
The identity (3.11), after integration in time, yields a stability bound for the discrete elastic energy.
4 Fully discrete finite element approximations
Let \(0= t_0< t_1< \ldots< t_{M-1} < t_M = T\) be a partitioning of [0, T] into possibly variable time steps \(\Delta t_m = t_{m+1} - t_{m}\), \(m=0,\ldots , M-1\). For a given \(\vec {X}^m\in \underline{V}^h\) we let \(\vec \nu ^m\) and \(\vec \omega ^m\) be the fully discrete analogues to (3.3) and (3.4), respectively.
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. \((\nabla \,g)({\vec {z}}) = U^T\, (\nabla _\mathbf{u}\,\Psi )(\mathbf{u}_0 + U\,{\vec {z}})\).
4.1 Curvature flow
We consider the following linear fully discrete analogue of \(({\mathcal {A}}_h)^h\).
\(({\mathcal {A}}_m)^h\): Let \({\vec {X}}^0 \in \underline{V}^h_{\partial _0}\). For \(m=0,\ldots ,M-1\), find \(({\vec {X}}^{m+1}, \kappa ^{m+1}) \in \underline{V}^h\times V^h\), with \({\vec {X}}^{m+1} - {\vec {X}}^m \in {\mathbb {X}}^h\), such that
Note that the scheme \(({\mathcal {A}}_m)^h\) is a natural generalisation of the scheme [12, \(({\mathcal {A}}_m)^h\)] to the case of open curves. The scheme \(({\mathcal {A}}_m)^h\) has the advantage that it is linear, recall (3.6), and that it asymptotically inherits the equidistribution property from \(({\mathcal {A}}_h)^h\), (3.5). 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 [5, (2.3)]. In fact, for a closed curve in the Euclidean plane the scheme \(({\mathcal {A}}_m)^h\) collapses to that approximation of curvature flow.
We make the following mild assumption.
Lemma 4.1
Let the assumption \(({\mathfrak {A}})^h\) hold. Then there exists a unique solution \(({\vec {X}}^{m+1}, \kappa ^{m+1}) \in \underline{V}^h\times V^h\) to \(({\mathcal {A}}_m)^h\).
Proof
As (4.1a), (4.1b) is linear, recall (3.6), existence follows from uniqueness. To investigate the latter, we consider the system: Find \((\delta {\vec {X}},\kappa ) \in {\mathbb {X}}^h \times V^h\) such that
where we recall from (3.6) that \(\lambda \in V^h\) with \(\lambda > 0\) in \({\overline{I}}\). It immediately follows from (4.2a) that \(\kappa = 0\) on \(\partial _0 I\), and so \(\kappa \in W^h_{\partial _0}\). Choosing \(\vec \eta = \delta {\vec {X}} \in {\mathbb {X}}^h\) in (4.2b) and \(\chi ={{\hat{\kappa }}} \in W^h_{\partial _0}\) in (4.2a), with \({{\hat{\kappa }}}(q_j) = g^{-1}({\vec {X}}^m(q_j))\,\kappa (q_j)\) for \(q_j \in {\overline{I}} \setminus \partial _0 I\), yields that
It follows from (4.3) that \(\kappa = {{\hat{\kappa }}} = 0\) and that \(\delta {\vec {X}}\) is constant. Hence (4.2a) and (3.1) imply that
It follows from (4.4) and assumption \(({\mathfrak {A}})^h\) that \(\delta {\vec {X}}=\vec {0}\). Hence we have shown that \(({\mathcal {A}}_m)^h\) has a unique solution \(({\vec {X}}^{m+1},\kappa ^{m+1}) \in \underline{V}^h\times V^h\). \(\square \)
In order to present an unconditionally stable fully discrete approximation of \(({\mathcal {C}}_{h})^{h}\), we assume that we can split \(g^\frac{1}{2}\) into
It follows from the splitting in (4.5) that
Then we introduce the following nonlinear scheme.
\(({\mathcal {C}}_{m,\star })^{h}\): Let \({\vec {X}}^0 \in \underline{V}^h_{\partial _0}\). For \(m=0,\ldots ,M-1\), find \(({\vec {X}}^{m+1}, \kappa _g^{m+1}) \in \underline{V}^h\times W^h_{\partial _0}\), with \({\vec {X}}^{m+1} - {\vec {X}}^m \in {\mathbb {X}}^h\), such that
Note that the scheme \(({\mathcal {C}}_{m,\star })^h\) is a natural generalisation of the scheme [12, \(({\mathcal {C}}_{m,\star })^h\)] to the case of open curves.
We can prove the following fully discrete analogue of Theorem 3.1.
Theorem 4.2
Let \(({\vec {X}}^{m+1},\kappa _g^{m+1})\) be a solution to \(({\mathcal {C}}_{m,\star })^{h}\). Then it holds that
Proof
Choosing \(\chi = \Delta t_m\,\kappa _g^{m+1}\in W^h_{\partial _0}\) in (4.7a) and \(\vec \eta = {\vec {X}}^{m+1} - {\vec {X}}^m \in {\mathbb {X}}^h\) in (4.7b) yields that
where we have used (4.6) and the inequality \({\vec {a}}\,.\,({\vec {a}} - {\vec {b}}) \ge |{\vec {b}}|\,(|{\vec {a}}| - |{\vec {b}}|)\) for \({\vec {a}}\), \({\vec {b}} \in {{\mathbb {R}}}^2\).
\(\square \)
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^\frac{1}{2}({\vec {z}}) = \frac{{\mathfrak {b}}^3}{[1 - {\mathfrak {b}}^2]^\frac{1}{2}}\, e^{{\mathfrak {b}}\,{\vec {z}}\,.\,\vec {e}_1}\,\vec {e}_1\otimes \vec {e}_1\) is clearly positive semidefinite, and so we can choose \(g^\frac{1}{2}_+ = g^\frac{1}{2}\) and \(g^\frac{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\ge 3\) to the reader. If \(n=2\), then we note that
and we observe that the eigenvalues of \(\begin{pmatrix} 3\,{\vec {z}}\,.\,\vec {e}_1 &{} {\vec {z}}\,.\,\vec {e}_2 \\ {\vec {z}}\,.\,\vec {e}_2 &{} {\vec {z}}\,.\,\vec {e}_1 \end{pmatrix}\) are \(2\,{\vec {z}}\,.\,\vec {e}_1 \pm |{\vec {z}}|\). Moreover, the function \({\mathcal {F}}({\vec {z}}) = \tfrac{1}{2}\,e^{-\frac{1}{4}\, |{\vec {z}}|^2} \, (2\,{\vec {z}}\,.\,\vec {e}_1 + |{\vec {z}}|)\) attains its maximum at \({\vec {z}} = \sqrt{2}\,\vec {e}_1\) with \(\max _{{\vec {z}} \in {{\mathbb {R}}}^2} {\mathcal {F}}({\vec {z}}) = \frac{3}{\sqrt{2\,e}} \approx 1.2866\). Hence the matrix
is positive definite for all \({\vec {z}} \in H\), and so we can choose
4.2 Elastic flow
We consider the following linear fully discrete analogue of the scheme \(({\mathcal {Q}}_h)^\diamond \), (3.10).
\(({\mathcal {Q}}_m)^\diamond \): Let \(({\vec {X}}^0,\kappa ^0_g,{\vec {Y}}^0_g) \in \underline{V}^h_{\partial _0}\times V^h\times \underline{V}^h\). For \(m=0,\ldots ,M-1\), find \(({\vec {X}}^{m+1}, \kappa ^{m+1}_g, {\vec {Y}}^{m+1}_g) \in \underline{V}^h\times V^h \times {\mathbb {Y}}^h\), with \({\vec {X}}^{m+1} - {\vec {X}}^m \in {\mathbb {X}}^h\), such that
Note that the scheme \(({\mathcal {Q}}_{m})^\diamond \) is a natural generalisation of the scheme [11, \(({\mathcal {Q}}_{m})^\diamond \)] 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 \({\vec {Y}}^{m+1}_g\) and \({\vec {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. [7, 8, 11, 14].
We make the following mild assumptions.
In the case \((\cdot ,\cdot )^\diamond = (\cdot ,\cdot )^h\) the above assumption collapses to \(({\mathfrak {A}})^h\). When dealing with clamped boundary conditions, we also need the following assumption, which is similar to [14, Assumption 5.9].
Lemma 4.3
Let the assumptions \(({\mathfrak {A}})^h\) and \(({\mathfrak {B}})^\diamond \) hold. Moreover, if \(\partial _C I \not = \emptyset \) then let assumptions \(({\mathfrak {C}})^\diamond \) hold. Then there exists a unique solution \(({\vec {X}}^{m+1}, \kappa ^{m+1}_g, {\vec {Y}}^{m+1}_g) \in \underline{V}^h\times V^h \times {\mathbb {Y}}^h\) to \(({\mathcal {Q}}_m)^\diamond \), if the quadrature rule (3.2) has at least one interior sampling point, \(\alpha _k\in (0,1)\). If \((\cdot ,\cdot )^\diamond = (\cdot ,\cdot )^h\), on the other hand, then there exists a solution that can be made unique by requiring that \(\kappa ^{m+1}_g \in W^h_{\partial _0}\).
Proof
We first consider the case that (3.2) is such that \(\alpha _k \in (0,1)\) for some \(1\le k\le K\). As (4.9) is linear, existence follows from uniqueness. To investigate the latter, we consider the system: Find \((\delta {\vec {X}},\kappa _g, {\vec {Y}}_g) \in {\mathbb {X}}^h \times V^h \times {\mathbb {Y}}^h\) such that
Choosing \(\vec \chi = \delta {\vec {X}}\in {\mathbb {X}}^h\) in (4.10a), \(\chi = \kappa _g\) in (4.10b) and \(\vec \eta = {\vec {Y}}_g \in {\mathbb {Y}}^h\) in (4.10c) yields that
First of all it follows from (4.11), our assumption on (3.2) and the positivities of \(g({\vec {X}}^m)\) and \(|{\vec {X}}^m_\rho |_g\) on \({\overline{I}}\setminus \partial _0 I\), that \(\kappa _g = 0 \in V^h\). As a consequence, we obtain by choosing \(\vec \eta = \delta {\vec {X}} \in {\mathbb {X}}^h \subset {\mathbb {Y}}^h\) in (4.10c) that \(\delta {\vec {X}}\) is a constant vector. Now (4.11) implies that this constant is such that \(\delta {\vec {X}}\,.\,\vec \omega ^m(q_j) = 0\) for all \(q_j \in {\overline{I}}\setminus \partial _0 I\), and so the assumption \(({\mathfrak {A}})^h\) yields that \(\delta {\vec {X}} = {\vec {0}}\).
It remains to show that \({\vec {Y}}_g = {\vec {0}}\). If \(\partial _C I = \emptyset \), then we can choose \(\vec \chi = {\vec {Y}}_g \in {\mathbb {Y}}^h \subset {\mathbb {X}}^h\) in (4.10a) to obtain that \({\vec {Y}}_g\) is constant in \({\overline{I}}\). Combining (4.10b) with assumption \(({\mathfrak {B}})^\diamond \) then gives that \({\vec {Y}}_g = {\vec {0}}\). If \(\partial _C I \not = \emptyset \), on the other hand, then assumption \(({\mathfrak {C}})^\diamond \) directly gives that \({\vec {Y}}_g={\vec {0}}\), in view of (4.10a) and (4.10b). Hence there exists a unique solution \(({\vec {X}}^{m+1}, \kappa ^{m+1}_g, {\vec {Y}}^{m+1}_g) \in \underline{V}^h\times V^h \times {\mathbb {Y}}^h\) to \(({\mathcal {Q}}_m)^\diamond \).
For the case \((\cdot ,\cdot )^\diamond = (\cdot ,\cdot )^h\) we note that \(V^h\) in (4.9b) can be equivalently replaced by \(W^h_{\partial _0}\). Existence of a unique \(({\vec {X}}^{m+1}, \kappa ^{m+1}_g, {\vec {Y}}^{m+1}_g) \in \underline{V}^h\times W^h_{\partial _0} \times {\mathbb {Y}}^h\) to this new system can then be shown similarly to the above proof, which gives all the remaining desired results. \(\square \)
Remark 4.4
We note that in the examples (1.4d), (1.4e) and (1.5b), any closed curve \({\vec {x}}(I)\) in H will correspond to a curve \(\vec \Phi ({\vec {x}}(I))\) on the hypersurface \({\mathcal {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 \({{\mathbb {R}}}^2\). In particular, \(H = {{\mathbb {R}}}\times {{\mathbb {R}}} / (2\,\pi \,{{\mathbb {Z}}})\) for (1.4d) and (1.5b), and \(H = {{\mathbb {R}}} / (2\,\pi \,{\mathfrak {s}}\,{{\mathbb {Z}}})\times {{\mathbb {R}}} / (2\,\pi \,{{\mathbb {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 \({\vec {X}}^m(I_j)\) some care needs to be taken when selecting representatives of the endpoints for the computation of \({\vec {X}}^m_\rho \), which then naturally yields \(|{\vec {X}}^m_\rho |\) and \(\vec \nu ^m\). We will present some numerical simulations for closed curves that are not homotopic to a point in Sect. 5.
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 \(({\mathcal {C}}_{m,\star })^{h}\) are computed with a Newton iteration. The two schemes \(({\mathcal {A}}_{m})^{h}\) and \(({\mathcal {C}}_{m,\star })^{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 \(({\mathcal {C}}_{m,\star })^{h}\).
We note from (3.8) and (2.8) that \(L^h_g({\vec {X}}^m)\) acts as a discrete energy for \(({\mathcal {A}}_{m})^{h}\) and \(({\mathcal {C}}_{m,\star })^{h}\), while on recalling Theorem 3.2 we define \({\widetilde{W}}_{g}^{m+1} = \tfrac{1}{2}\, ( (\kappa _g^{m+1})^2,|{\vec {X}}^m_\rho |_g )^\diamond \) as a discrete analogue of (2.10) for the scheme \(({\mathcal {Q}}_m)^\diamond \). As the quadrature rule for the scheme \(({\mathcal {Q}}_m)^\diamond \) we either consider (3.1), leading to \(({\mathcal {Q}}_m)^h\), or a quadrature that is exact for polynomials of degree up to five, denoted by \((\cdot ,\cdot )^\star \), and so leading to \(({\mathcal {Q}}_m)^\star \). 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 \(({\mathcal {Q}}_m)^\diamond \), given \(\Gamma ^0 = {\vec {X}}^0({\overline{I}})\), is defined as follows. First we define \(\kappa ^0 \in V^h\) via \(\kappa ^0(q_j) = [|\vec \omega ^0|^{-2}\,\vec \kappa ^0\,.\,\vec \omega ^0](q_j)\) for \(j=0,\ldots ,J\), where \(\vec \kappa ^0\in \underline{V}^h\) is such that
Then let \(\kappa _g^0 \in W^h_{\partial _0}\) with \(\kappa _g^0(q_j) = {\mathcal {K}}(\kappa ^0, \vec \omega ^0, {\vec {X}}^0)(q_j)\) for \(q_j \in {\overline{I}} \setminus \partial _0 I\), recall (3.6). In addition, let \({\vec {Y}}_g^0 \in [W^h_{\partial _0}]^2\) with \({\vec {Y}}_g^0 = [g^{-\frac{1}{2}}({\vec {X}}^0)\,|\vec \omega ^0|^{-2}\, \kappa _g^0\,\vec \omega ^0](q_j)\) for \(q_j \in {\overline{I}} \setminus \partial _0 I\).
In most of the presented simulations we use uniform time steps, \(\Delta t_m = \Delta t\), \(m=0,\ldots ,M-1\). For some simulations, however, we use an adaptive time step strategy satisfying \(\Delta t_{\min } \le \Delta t_m \le \Delta t_{\max }\), \(m=0,\ldots ,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 \(\Delta t= 10^{-4}\).
5.1 The metric (1.4a)
For the scheme \(({\mathcal {A}}_m)^h\) we show the evolution of two cigar shapes in Fig. 1 for the metric (1.4a) with \(\mu =1\). We note that in both cases the curve shrinks to a point.
\(({\mathcal {A}}_m)^h\) Curvature flow towards extinction for (1.4a) with \(\mu =1\). Solution at times \(t=0,0.05,\ldots ,0.2\) (left), and at times \(t=0,0.1,\ldots ,0.5,0.55\) (right)
Repeating the same evolutions for the metric (1.4a) with \(\mu =-1\), now using the scheme \(({\mathcal {C}}_{m,\star })^{h}\), leads to the results shown in Fig. 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 \(({\mathcal {C}}_{m,\star })^{h}\) no longer converges. Here we note that we used the smaller uniform time step size \(\Delta t=10^{-5}\) for this experiment.
\(({\mathcal {C}}_{m,\star })^{h}\) Curvature flow towards extinction for (1.4a) with \(\mu =-1\). Solution at times \(t=0,1,\ldots ,4,4.5\) (left) and at times \(t=0,0.01,0.015,0.0156\) (right)
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 [10] for details. As an example, we show the evolution of a semicircle with radius 1 and \(\partial _0 I = \partial I\) in Fig. 3. As a comparison, we also show the same evolution for the case \(\partial _1 I = \partial I\). In both cases, the semicircle shrinks to extinction, but the shape and time scale of the two evolutions differ.
\(({\mathcal {C}}_{m,\star })^{h}\) Curvature flow for (1.4a) with \(\mu =-1\). Solution for \(\partial I = \partial _0 I = \{0,1\}\) at times \(t=0,0.02,\ldots ,0.08,0.085\) (left) and for \(\partial I = \partial _1 I = \{0,1\}\) at times \(t=0,0.1,\ldots ,0.3,0.34\) (right)
For completeness, we also show some evolutions for the cases \(\partial _D I = \partial I\) and \(\partial _2 I = \partial I\) in Fig. 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 Fig. 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.
\(({\mathcal {C}}_{m,\star })^{h}\) Curvature flow for (1.4a) with \(\mu =-1\) with \(\partial I = \partial _D I = \{0,1\}\) (left and middle) and \(\partial I = \partial _2 I = \{0,1\}\) (right). Solution at times \(t=0,0.1,\ldots ,0.4,0.447\) (left), \(t=0,2,\ldots ,10\) (middle) and \(t=0,1,\ldots ,3,3.3\) (right). Below we show plots of the discrete energies \(L_g^h({\vec {X}}^m)\)
Evolutions for elastic flow with Navier and clamped boundary conditions, respectively, are shown in Fig. 5. Here, for the clamped boundary conditions, recall (2.27), we choose \(\vec \zeta (p) = (\sin \vartheta (p), \cos \vartheta (p))^T\), with \(\vartheta (0) = 210^\circ \) and \(\vartheta (1) = 150^\circ \). 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.
\(({\mathcal {Q}}_m)^{\star }\) Elastic flow for (1.4a) with \(\mu =-1\) and \(\partial I = \partial _N I = \{0,1\}\) (top) and \(\partial I = \partial _D I = \{0,1\}\) (bottom). Solution at times \(t=0,1,\ldots ,5\) (above) and at times \(t=0,1,\ldots ,5\) (below). We also show plots of the discrete energy \({\widetilde{W}}_g^{m+1}\) over time
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), including for the case case \(\mu =1\) for the hyperbolic plane, can be found in [11, 12].
5.2 The torus metric (1.4e)
A geodesic between two fixed points on the Clifford torus is computed in Fig. 6. To this end, we employ the metric induced by (1.4e) with \({\mathfrak {s}} = 1\), so that the torus has radii \(r=1\) and \(R = \sqrt{2}\). We observe that the evolution eventually settles on a geodesic, that is clearly not the shortest path connecting 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.
\(({\mathcal {A}}_{m})^h\) Geodesic curvature flow on a Clifford torus, with \(\partial _D I = \partial I = \{0,1\}\). The solutions \({\vec {X}}^m\) at times \(t = 0, 2,\ldots , 6\). Below we visualise \(\vec \Phi ({\vec {X}}^m)\) at times \(t=0\) (blue), \(t=2\) (red) and \(t=6\) (black), for (1.4e) with \({\mathfrak {s}}=1\), and also show a plot of the discrete energy \(L_g^h({\vec {X}}^m)\)
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 Fig. 7 for a presentation of the numerical results for the scheme \(({\mathcal {C}}_{m,\star })^h\).
\(({\mathcal {C}}_{m,\star })^h\) Geodesic curvature flow on a Clifford torus. We visualise \(\vec \Phi ({\vec {X}}^m)\) at times \(t=0,1,3\), for (1.4e) with \({\mathfrak {s}}=1\). A plot of the discrete energy \(L_g^h({\vec {X}}^m)\) below
5.3 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 elongated cigar shape that shrinks to a point, see in Fig. 8.
\(({\mathcal {C}}_{m,\star })^{h}\) Curvature flow towards extinction for (1.5a). Solution at times \(t=0,0.1,0.2,0.25\)
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 \({{\mathbb {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 Fig. 9, we observe that the discrete curvature energy \(W_g^{m+1}\) 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 Fig. 9 agrees well with the numerical results in [3, 17, 22].
\(({\mathcal {Q}}_{m})^\star \) Elastic flow for (1.5a) towards the Angenent torus. Plots are at times \(t=0,0.1,0.5\). We also show a plot of the discrete energy \({\widetilde{W}}_g^{m+1}\) over time
We have also performed simulations for elastic flow of initial curves with a winding number larger than one, with respect to the point \(2\,\vec {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 [18, 23, 54], 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 \(({\mathcal {C}}_{m,\star })^h\) from slightly shifted Angenent tori. Our numerical results in Fig. 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 \(\Delta t=10^{-5}\).
\(({\mathcal {C}}_{m,\star })^h\) Curvature flow for (1.5a), starting from horizontally shifted Angenent tori. Above shifted by 0.05 to the right, below shifted by 0.05 to the left. Plots are at times \(t=0,0.18,0.2\) (above) and at times \(t=0,0.02,0.03,0.038\) (below). We also show plots of the discrete energy \(L_g^h({\vec {X}}^m)\) over time
We highlight the capabilities of our numerical method by computing the “Angenent tori” in dimensions four and five, that is hypersurfaces in \({{\mathbb {R}}}^{n+1}\) that are topologically equivalent to \({{\mathbb {S}}}^1 \times {{\mathbb {S}}}^{n-1}\), \(n=3,4\), and that are self-shrinkers for mean curvature flow with extinction time 1. In particular, in Fig. 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 \(|{\widetilde{W}}_g^M| < 10^{-9}\), confirming that we are indeed approximating geodesics.
\(({\mathcal {Q}}_{m})^\star \) Steady states for elastic flow for (1.5a) with \(n=2,3,4\). Their discrete geodesic lengths are 3.70, 6.39 and 14.27
We note that in the study of mean curvature flow self-shrinkers the value of Huisken’s F-functional itself is also a relevant quantity of interest, see e.g. [17, 23, 35, 61]. For example, for a self-shrinker its value of F is equal to its entropy. On recalling (1.6), (1.5a) and (2.8), we note that if \({\vec {y}}({\overline{I}}) \subset H\) is the generating curve for a rotationally symmetric hypersurface \({\mathcal {S}} \subset {{\mathbb {R}}}^{n+1}\), then
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}\,[{\Gamma }(\tfrac{n}{2})]^{-1} = \frac{1}{2}\), \(\frac{1}{2\,\sqrt{\pi }}\) and \(\tfrac{1}{8}\), respectively. For the three geodesics displayed in Fig. 11, our computed values for the final discrete length \(L_g^h({\vec {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 [3, 17].
Denoting the entropy of the n-dimensional “Angenent torus” with \(E_n\), we have so far established that \(E_2 \approx 1.85\), \(E_3 \approx 1.80\) and \(E_4 \approx 1.78\). Continuing the above procedure for increasing values of n, we numerically obtain that \(E_6 \approx 1.77\), \(E_8 \approx 1.76\) and \(E_{12} \approx 1.75\). This leads us to conjecture that \(E_n\) is a strictly monotonically decreasing sequence with \(E_n \searrow \sqrt{3}\) as \(n\rightarrow \infty \). 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 \({{\mathbb {R}}}^{n+1}\), with extinction time 1, is the sphere of radius \(\sqrt{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 \(\sqrt{2\,n}\). For \(n=2\) and \(n=3\) we show an evolution each for elastic flow towards these geodesics, see Fig. 12 for details, where in each case as initial data we choose a semicircle of radius \(n-1\). For the two geodesics displayed in Fig. 12, our computed values for the final discrete length \(L_g^h({\vec {X}}^M)\) are 2.94 and 5.15, giving approximate values for the entropy of the associated spheres of 1.47 and 1.45. We remark that these values agree very well with the known values from [61, Examples A.3], which are known to converge to \(\sqrt{2}\) as \(n\rightarrow \infty \), recall [61, Lemma A.4].
\(({\mathcal {Q}}_m)^{\star }\) Elastic flow for (1.5a), \(n=2\) (left) and \(n=3\) (right), and \(\partial I = \partial _0 I = \{0,1\}\). Solutions at times \(t=0,0.1,1\) (left) and at times \(t=0,0.1,1,3\) (right). We also show a plot of the discrete energy \({\widetilde{W}}_g^{m+1}\) over time. The discrete lengths of the obtained geodesics are 2.94 (\(n=2\)) and 5.15 (\(n=3\))
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 \({\mathcal {E}}^M = \max _{j=0,\ldots ,J} | 2 - |{\vec {X}}^M(q_j)| |\) and \(e^M = |\frac{1}{2} L_g^h({\vec {X}}^M) - \frac{4}{e}|\) for increasing values of J. Here we fix \(T=10\) and \(\Delta 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^\infty \) error \({\mathcal {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 Fig. 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, Fig. 6]. Due to the large energy decrease at the beginning of the evolution, we use an adaptive time stepping strategy with \(\Delta t_{\min } = 10^{-7}\) and \(\Delta t_{\max } = 10^{-6}\). The spatial discretisation uses \(J=512\). The discrete energy of the final solution satisfies \(|{\widetilde{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.
\(({\mathcal {Q}}_{m})^\star \) Elastic flow for (1.5a). Solution at times \(t=0,0.01\), \(t=1\) and \(t=20\). The discrete geodesic length of the final curve is 10.53
We also investigate what happens to the geodesic from Fig. 13 if we change the metric to (1.5a) with \(n=3\). See Fig. 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 \(|{\widetilde{W}}^M_g| < 10^{-8}\). The discrete geodesic length of the final curve is 18.24, giving an entropy of 5.15.
\(({\mathcal {Q}}_{m})^\star \) Elastic flow for (1.5a) with \(n=3\). Solution at times \(t=0,1,100\). The discrete geodesic length of the final curve is 18.24
Inspired by [34, Fig. 3], we now perform a numerical simulation to find a non-embedded 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 Fig. 15, where the final discrete energy satisfies \(|{\widetilde{W}}^M_g| < 10^{-9}\). We note the excellent agreement with [34, Fig. 3]. Here we again make use of an adaptive time stepping strategy with \(\Delta t_{\min } = 10^{-7}\) and \(\Delta 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.
\(({\mathcal {Q}}_{m})^\star \) Elastic flow for (1.5a) and \(\partial I = \partial _0 I = \{0,1\}\). Solution at times \(t=0,0.01\), \(t=0.1\) and \(t=10\). The discrete geodesic length of the final curve is 7.33
5.4 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 \({\mathfrak {b}} = 0.5\), and so \(\beta = 3^{-\frac{1}{2}}\) in (1.7). For the simulation in Fig. 16 it can be observed that in H the initial circle of radius 2 deforms and shrinks to a point. On the hypersurface \({\mathcal {M}} = \vec \Phi (H)\), the initial curve is homotopic to a point, and so shrinks to a point away from the apex.
\(({\mathcal {A}}_m)^h\) Geodesic curvature flow on a cone. The solutions \({\vec {X}}^m\) at times \(t = 0, 0.5, 1\). On the right we visualise \(\vec \Phi ({\vec {X}}^m)\) at times \(t=0, 0.5, 1\), for (1.5b) with \({\mathfrak {b}}=0.5\). A plot of the discrete energy \(L_g^h({\vec {X}}^m)\) is shown on the right
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 \({\mathcal {M}}\), that is not homotopic to a point on \({\mathcal {M}}\), will under geodesic curvature flow converge to the apex in finite time.
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 \(\pm 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 Fig. 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
\(({\mathcal {A}}_m)^{h}\) Geodesic curvature flow on a cone. We visualise \(\vec \Phi ({\vec {X}}^m)\) at times \(t=0, 0.5, 1.5\), for (1.5b) with \({\mathfrak {b}}=0.5\). A plot of the discrete energy \(L_g^h({\vec {X}}^m)\) over the time interval [0, 1.52] in the middle. On the right a plot of the lowest point of the curve on the cone, \(\exp ({\mathfrak {b}}\,\min _{{\overline{I}}} {\vec {X}}^m\,.\,\vec {e}_1)\)
An experiment for (geodesic) elastic flow on the same cone is shown in Fig. 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.
\(({\mathcal {Q}}_{m,\star })^{h}\) Geodesic elastic flow on a cone. We visualise \(\vec \Phi ({\vec {X}}^m)\) at times \(t=0, 10, \ldots , 50\), for (1.5b) with \({\mathfrak {b}}=0.5\). A plot of the discrete energy \({\widetilde{W}}_g^{m+1}\) on the right
5.5 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 \(\mathbf{e}_1\), \(\mathbf{e}_2\), \(\mathbf{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 \({\vec {z}} \mapsto f({\vec {z}}) = \mathbf{u}_0 + U\,{\vec {z}}\) satisfies \(f(0,0) = \mathbf{e}_1\), \(f(-2^\frac{1}{2},0) = \mathbf{e}_2\) and \(f(-2^\frac{1}{2},-(\frac{3}{2})^\frac{1}{2}) = \mathbf{e}_3\). For the first experiment we set \((\sigma _{12},\sigma _{13},\sigma _{23},\sigma _{123}) = (4,6,9,0)\), and numerically compute a geodesic connecting \(\mathbf{e}_1\) and \(\mathbf{e}_2\) with the help of geodesic curvature flow. Here we always use the scheme \(({\mathcal {A}}_{m})^{h}\) with the uniform time step size \(\Delta t= 10^{-5}\). The results are shown in Fig. 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 \(\mathbf{e}_1\) and \(\mathbf{e}_3\), as well as \(\mathbf{e}_2\) and \(\mathbf{e}_3\), and plot all three solutions within the Gibbs simplex G, recall (1.10), at the bottom of Fig. 19.
\(({\mathcal {A}}_{m})^{h}\) Geodesic curvature flow for the metric (1.5c) with (1.9) and \((\sigma _{12},\sigma _{13},\sigma _{23},\sigma _{123}) = (4,6,9,0)\). The solutions \({\vec {X}}^m\) at times \(t = 0, 0.01, 0.1\). A plot of the discrete energy \(L_g^h({\vec {X}}^m)\) on the right. Below a plot of the three minimisers connecting the vertices of the Gibbs simplex
In [41] numerical computations indicated that on choosing \(\sigma _{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 \(\sigma _{123} \in \{10,100,1000\}\) and plot the obtained geodesics in Fig. 20. It can be seen that for an increasing value of \(\sigma _{123}\), the geodesics are pushed further and further towards the edges of the Gibbs simplex.
\(({\mathcal {A}}_{m})^{h}\) The minimisers obtained within the Gibbs simplex, for (1.9) with \((\sigma _{12},\sigma _{13},\sigma _{23}) = (4,6,9)\) and \(\sigma _{123} = 10,100,1000\) (from left to right)
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 [38, 42], we can consider the following replacement of (1.9) to achieve the same effect:
where \({{\hat{\sigma }}}_{123},{{\hat{\sigma }}}_{231},{{\hat{\sigma }}}_{312} \in {{\mathbb {R}}}_{\ge 0}\). We perform a computation for (5.2) with \(\sigma _{12}=\sigma _{13}=\sigma _{23}={{\hat{\sigma }}}_{123}={{\hat{\sigma }}}_{231}= {{\hat{\sigma }}}_{312}=1\) and show the obtained results in Fig. 21. It can be seen that now the geodesics lie on the edges of the Gibbs simplex, confirming the analysis in [38, 42].
\(({\mathcal {A}}_{m})^{h}\) The minimisers obtained for (5.2) with \(\sigma _{12}=\sigma _{13}=\sigma _{23}={{\hat{\sigma }}}_{123}={{\hat{\sigma }}}_{231}= {{\hat{\sigma }}}_{312}=1\) are precisely the edges of the Gibbs simplex
References
Andrews, B., Chen, X.: Curvature flow in hyperbolic spaces. J. Reine Angew. Math. 729, 29–49 (2017)
Angenent, S.B.: Shrinking doughnuts. In: Nonlinear Diffusion Equations and their Equilibrium States, 3 (Gregynog, 1989), Progr. Nonlinear Differential Equations Appl., vol. 7, pp. 21–38. Birkhäuser Boston, Boston, MA (1992)
Barrett, J.W., Deckelnick, K., Nürnberg, R.: A finite element error analysis for axisymmetric mean curvature flow. IMA J. Numer. Anal. 41, 1641–1667 (2021)
Barrett, J.W., Garcke, H., Nürnberg, R.: A parametric finite element method for fourth order geometric evolution equations. J. Comput. Phys. 222, 441–462 (2007)
Barrett, J.W., Garcke, H., Nürnberg, R.: On the variational approximation of combined second and fourth order geometric evolution equations. SIAM J. Sci. Comput. 29, 1006–1041 (2007)
Barrett, J.W., Garcke, H., Nürnberg, R.: Numerical approximation of gradient flows for closed curves in \({{\mathbb{R}}}^d\). IMA J. Numer. Anal. 30, 4–60 (2010)
Barrett, J.W., Garcke, H., Nürnberg, R.: Parametric approximation of isotropic and anisotropic elastic flow for closed and open curves. Numer. Math. 120, 489–542 (2012)
Barrett, J.W., Garcke, H., Nürnberg, R.: Elastic flow with junctions: variational approximation and applications to nonlinear splines. Math. Models Methods Appl. Sci. 22, 1250,037 (2012)
Barrett, J.W., Garcke, H., Nürnberg, R.: Finite element methods for fourth order axisymmetric geometric evolution equations. J. Comput. Phys. 376, 733–766 (2019)
Barrett, J.W., Garcke, H., Nürnberg, R.: Variational discretization of axisymmetric curvature flows. Numer. Math. 141, 791–837 (2019)
Barrett, J.W., Garcke, H., Nürnberg, R.: Stable discretizations of elastic flow in Riemannian manifolds. SIAM J. Numer. Anal. 57, 1987–2018 (2019)
Barrett, J.W., Garcke, H., Nürnberg, R.: Numerical approximation of curve evolutions in Riemannian manifolds. IMA J. Numer. Anal. 40, 1601–1651 (2020)
Barrett, J.W., Garcke, H., Nürnberg, R.: Parametric finite element approximations of curvature driven interface evolutions. In: Bonito, A., Nochetto, R.H. (eds.) Handbook of Numerical Analysis, vol. 21, pp. 275–423. Elsevier, Amsterdam (2020)
Barrett, J.W., Garcke, H., Nürnberg, R.: Stable approximations for axisymmetric Willmore flow for closed and open surfaces. Math. Model. Numer. Anal. 55, 833–885 (2021)
Bartels, S.: A simple scheme for the approximation of the elastic flow of inextensible curves. IMA J. Numer. Anal. 33, 1115–1125 (2013)
Benninghoff, H., Garcke, H.: Segmentation and restoration of images on surfaces by parametric active contours with topology changes. J. Math. Imaging Vis. 55, 105–124 (2016)
Berchenko-Kogan, Y.: The entropy of the Angenent torus is approximately 1.85122. Experiment. Math. (2019). https://doi.org/10.1080/10586458.2019.1583616
Berchenko-Kogan, Y.: Numerically computing the index of mean curvature flow self-shrinkers. arXiv:2007.06094 (2020)
Bretin, E., Masnou, S.: A new phase field model for inhomogeneous minimal partitions, and applications to droplets dynamics. Interfaces Free Bound. 19, 141–182 (2017)
Cabezas-Rivas, E., Miquel, V.: Volume preserving mean curvature flow in the hyperbolic space. Indiana Univ. Math. J. 56, 2061–2086 (2007)
Cheng, L.T., Burchard, P., Merriman, B., Osher, S.: Motion of curves constrained on surfaces using a level-set approach. J. Comput. Phys. 175, 604–644 (2002)
Chopp, D.L.: Computation of self-similar solutions for mean curvature flow. Exp. Math. 3, 1–15 (1994)
Colding, T.H., Minicozzi, W.P., II.: Generic mean curvature flow I: generic singularities. Ann. Math. 2(175), 755–833 (2012)
Dall’Acqua, A., Laux, T., Lin, C.C., Pozzi, P., Spener, A.: The elastic flow of curves on the sphere. Geom. Flows 3, 1–13 (2018)
Dall’Acqua, A., Lin, C.C., Pozzi, P.: Evolution of open elastic curves in \(\mathbb{R}^n\) subject to fixed length and natural boundary conditions. Analysis (Berlin) 34, 209–222 (2014)
Dall’Acqua, A., Lin, C.C., Pozzi, P.: A gradient flow for open elastic curves with fixed length and clamped ends. Ann. Sc. Norm. Super. Pisa Cl. Sci. 5(17), 1031–1066 (2017)
Dall’Acqua, A., Lin, C.C., Pozzi, P.: Elastic flow of networks: long-time existence result. Geom. Flows 4, 83–136 (2019)
Dall’Acqua, A., Pozzi, P.: A Willmore-Helfrich \(L^2\)-flow of curves with natural boundary conditions. Commun. Anal. Geom. 22, 617–669 (2014)
Dall’Acqua, A., Spener, A.: The elastic flow of curves in the hyperbolic plane. arXiv:1710.09600 (2017)
Dall’Acqua, A., Spener, A.: Circular solutions to the elastic flow in hyperbolic space. In: Proceedings of analysis on shapes of solutions to partial differential equations, (2017), RIMS Kôkyûroku, vol. 2082. Kyoto, Japan (2018)
Davis, T.A.: Algorithm 832: UMFPACK V4.3—an unsymmetric-pattern multifrontal method. ACM Trans. Math. Software 30, 196–199 (2004)
Deckelnick, K., Dziuk, G.: Error analysis for the elastic flow of parametrized curves. Math. Comput. 78, 645–671 (2009)
Deckelnick, K., Dziuk, G., Elliott, C.M.: Computation of geometric partial differential equations and mean curvature flow. Acta Numer. 14, 139–232 (2005)
Drugan, G., Kleene, S.J.: Immersed self-shrinkers. Trans. Am. Math. Soc. 369, 7213–7250 (2017)
Drugan, G., Nguyen, X.H.: Shrinking doughnuts via variational methods. J. Geom. Anal. 28, 3725–3746 (2018)
Elliott, C.M.: Private communication (2009)
Epstein, C.L., Gage, M.: The curve shortening flow. In: Wave Motion: Theory, Modelling, and Computation (Berkeley, Calif., 1986), Math. Sci. Res. Inst. Publ., vol. 7, pp. 15–59. Springer, New York (1987)
Garcke, H., Haas, R.: Modeling of Nonisothermal Multi-Component, Multi-Phase Systems with Convection, vol. chap. 20, pp. 325–338. Wiley, London (2008)
Garcke, H., Menzel, J., Pluda, A.: Willmore flow of planar networks. J. Differ. Equ. 266, 2019–2051 (2019)
Garcke, H., Menzel, J., Pluda, A.: Long time existence of solutions to an elastic flow of networks. Commun. Partial Differ. Equ. 45, 1253–1305 (2020)
Garcke, H., Stoth, B., Nestler, B.: A multiphase field concept: numerical simulations of moving phase boundaries and multiple junctions. SIAM J. Appl. Math. 60, 295–315 (1999)
Haas, R.: Modeling and analysis for general non-isothermal convective phase field systems. University Regensburg, Regensburg (2007) (Ph.D. thesis)
Helmers, M.: Kinks in two-phase lipid bilayer membranes. Calc. Var. Partial Differ. Equ. 48, 211–242 (2013)
Helmers, M.: Convergence of an approximation for rotationally symmetric two-phase lipid bilayer membranes. Q. J. Math. 66, 143–170 (2015)
Huisken, G.: Asymptotic behavior for singularities of the mean curvature flow. J. Differ. Geom. 31, 285–299 (1990)
Jost, J.: Riemannian Geometry and Geometric Analysis. Springer, Berlin (2005)
Kleene, S., Møller, N.M.: Self-shrinkers with a rotational symmetry. Trans. Am. Math. Soc. 366, 3943–3963 (2014)
Koiso, N.: On motion of an elastic wire in a Riemannian manifold and singular perturbation. Osaka J. Math. 52, 453–473 (2015)
Kraus, D., Roth, O.: Conformal metrics. In: Topics in Modern Function Theory, Ramanujan Math. Soc. Lect. Notes Ser., vol. 19, pp. 41–83. Ramanujan Math. Soc., Mysore, India (2013), (see also https://arxiv.org/abs/0805.2235)
Kühnel, W.: Differential Geometry: Curves - Surfaces - Manifolds, Student Mathematical Library, vol. 77. American Mathematical Society, Providence, RI (2015)
Langer, J., Singer, D.A.: The total squared curvature of closed curves. J. Differ. Geom. 20, 1–22 (1984)
Langer, J., Singer, D.A.: Curve-straightening in Riemannian manifolds. Ann. Global Anal. Geom. 5, 133–150 (1987)
Linnér, A.: Curve-straightening and the Palais–Smale condition. Trans. Am. Math. Soc. 350, 3743–3765 (1998)
Liu, Z.H.: The Morse index of mean curvature flow self-shrinkers. Massachusetts Institute of Technology, Ann Arbor, MI (2016) (Ph.D. thesis)
Maitre, E., Santosa, F.: Level set methods for optimization problems involving geometry and constraints. II. Optimization over a fixed surface. J. Comput. Phys. 227, 9596–9611 (2008)
Mikula, K., Ševčovič, D.: Evolution of curves on a surface driven by the geodesic curvature and external force. Appl. Anal. 85, 345–362 (2006)
Müller, M., Spener, A.: On the convergence of the elastic flow in the hyperbolic plane. Geom. Flows 5, 40–77 (2020)
Schippers, E.: The calculus of conformal metrics. Ann. Acad. Sci. Fenn. Math. 32, 497–521 (2007)
Schmidt, A., Siebert, K.G.: Design of Adaptive Finite Element Software: The Finite Element Toolbox ALBERTA. Lecture Notes in Computational Science and Engineering, vol. 42. Springer, Berlin (2005)
Spira, A., Kimmel, R.: Geometric curve flows on parametric manifolds. J. Comput. Phys. 223, 235–249 (2007)
Stone, A.: A density function and the structure of singularities of the mean curvature flow. Calc. Var. Partial Differ. Equ. 2, 443–480 (1994)
Acknowledgements
The authors gratefully acknowledge the support of the Regensburger Universitätsstiftung Hans Vielberth.
Funding
Open access funding provided by Università degli Studi di Trento within the CRUI-CARE Agreement.
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Garcke, H., Nürnberg, R. Numerical approximation of boundary value problems for curvature flow and elastic flow in Riemannian manifolds. Numer. Math. 149, 375–415 (2021). https://doi.org/10.1007/s00211-021-01231-6
Received:
Revised:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s00211-021-01231-6