-
PDF
- Split View
-
Views
-
Cite
Cite
Klaus Deckelnick, Robert Nürnberg, Discrete anisotropic curve shortening flow in higher codimension, IMA Journal of Numerical Analysis, Volume 45, Issue 1, January 2025, Pages 36–67, https://doi.org/10.1093/imanum/drae015
- Share Icon Share
Abstract
We introduce a novel formulation for the evolution of parametric curves by anisotropic curve shortening flow in |${{\mathbb{R}}}^{d}$|, |$d\geq 2$|. The reformulation hinges on a suitable manipulation of the parameterization’s tangential velocity, leading to a strictly parabolic differential equation. Moreover, the derived equation is in divergence form, giving rise to a natural variational numerical method. For a fully discrete finite element approximation based on piecewise linear elements we prove optimal error estimates. Numerical simulations confirm the theoretical results and demonstrate the practicality of the method.
1. Introduction
The aim of this paper is to define and analyse a finite element approximation for the evolution of a curve in |${{\mathbb{R}}}^{d}$|, for an arbitrary |$d\geq 2$|, by anisotropic curve shortening flow. The evolution law we consider is a natural gradient flow for the following anisotropic energy:
where |$\tau $| denotes the unit tangent of |$\varGamma $| and |$\phi $| is a given, |$1$|–homogeneous energy density; cf. Pozzi (2007). For the case |$d=2$|, i.e. curves in the plane, the energy (1.1) can play the role of an interfacial energy, e.g. in materials science (Taylor et al., 1992; Gurtin, 1993). A possible application in differential geometry is curve shortening flow in a Riemannian manifold. In this case the energy density in (1.1) is required to have a spatial dependence; see Bellettini & Paolini (1996); Ma & Chen (2007). For more details on anisotropic surface energies, we refer to Deckelnick et al. (2005) and Giga (2006), and the references therein. We shall see in Section 2 that a family of curves |$(\varGamma (t))_{t \in (0,T)}$| evolves according to anisotropic curve shortening flow provided that
where |$\mathcal V_{nor}$| denotes the vector of normal velocities, i.e. |$\mathcal V_{nor} \cdot \tau = 0$|, and |$\varkappa _{\phi }$| denotes the anisotropic curvature vector. In the isotropic case, |$\phi (p)=|p|$|, we have that |$\varkappa _{\phi }=\varkappa =\tau _{s}$|, with |$\cdot _{s}$| denoting differentiation with respect to arclength, and (1.2) is just the usual curve shortening flow. If |$\varGamma (t)= x(I,t)$| is described by a parameterization |$x:I \times [0,T] \to{{\mathbb{R}}}^{d}$|, where |$I={{\mathbb{R}}}/{{\mathbb{Z}}}$| is the periodic interval |$[0,1]$|, then (1.2) can be equivalently formulated as
where here and throughout we use a slight abuse of notation, in that we do not distinguish between geometric quantities of the curve being defined on |$\varGamma (t)$| or in |$I \times [0,T]$|. We observe that (1.3) does not prescribe the tangential velocity. From a geometric point of view it is natural to consider the system
see, e.g. Dziuk (1994) and Ambrosio & Soner (1996), for the isotropic setting. Generalizing the earlier works (Dziuk, 1994, 1999), Pozzi (2007) introduced and analysed a semidiscrete finite element method and obtained an |$\mathcal O(h)$| error bound in |$L^{\infty }(H^{1})$| for the position vector. Here, the fact that the anisotropic curvature vector is invariant with respect to reparameterization leads to a degeneracy of the elliptic part on the right-hand side of (1.4), thus complicating the error analysis. Barrett et al. (2010) suggested a discretization of (1.3) that introduces a tangential velocity at the discrete level, which leads to a nice distribution of vertices in practice. In Deckelnick & Nürnberg (2023b), the present authors considered anisotropic curve shortening flow in the planar case |$d=2$|. The crucial idea from Deckelnick & Nürnberg (2023b) was to define positive definite matrices |$H(x_{\rho }) \in{{\mathbb{R}}}^{2 \times 2}$| such that any solution to the PDE
parameterizes anisotropic curve shortening flow. We would like to stress three appealing aspects of (1.5). First, it can be shown that this PDE is strictly parabolic in the sense of Petrovsky; see Deckelnick & Nürnberg (2023b, Lemma 3.3), so that (1.5) can be interpreted as a kind of DeTurck trick, Elliott & Fritz (2017). Secondly, the PDE is in divergence form, making its variational numerical approximation straightforward, and for a semidiscrete finite element approximation an |$\mathcal O(h)$| error bound in |$L^{\infty }(H^{1})$| for the position vector is obtained. And thirdly, the tangential velocities induced by (1.5) lead to well distributed vertices in practice. It is the aim of this paper to generalise these results, in that we
propose an analogue of (1.5) for curves in higher codimension;
provide an error analysis for a fully discrete finite element scheme; and
prove an optimal |$L^{\infty }(L^{2})$| error bound.
Let us emphasize that all three of the above results are new in the literature. We note that in Deckelnick & Nürnberg (2023b) the anisotropic energy was allowed to depend on space. We expect that it is possible to treat such a spatial dependence also in higher codimensions, allowing to study, for example, curve shortening flow in higher dimensional Riemannian manifolds. However, in this paper we restrict our attention to the simpler energy (1.1) in order not to overburden the presentation.
Let us mention related work on the numerical approximation of (anisotropic) curve shortening flow in higher codimension. Dziuk (1994) and Deckelnick & Dziuk (1995) considered the isotropic case and proved optimal |$H^{1}$|-error bounds for semidiscrete approximations based on piecewise linear elements. In the planar case, error estimates for a fully discrete variant of Dziuk (1994) have recently been obtained in Li (2020) and Ye & Cui (2021). The method and numerical analysis in Dziuk (1994) were generalized to the anisotropic setting in Pozzi (2007). Both Dziuk (1994) and Pozzi (2007) may suffer from coalescence of mesh points in practice, since the discrete curves are only updated in normal direction. The so-called BGN schemes, on the other hand, are characterized by an implicit tangential motion that leads to a nice distribution of vertices (Barrett et al., 2020). Their application to (anisotropic) curve shortening flow in higher codimensions has been considered in Barrett et al. (2010), with an error analysis for these schemes still lacking. A semi-Lagrangian scheme in the context of level-set methods was considered in Carlini et al. (2007). A finite volume scheme for possibly interacting curves driven by curvature forces in |${{\mathbb{R}}}^{3}$| was introduced in Beneš et al. (2022). Finally, optimal error estimates for semi-discrete and fully discrete approximations of a system of PDEs for |$\varkappa $| and |$\tau $| describing isotropic curve shortening in higher codimension have been obtained in Binz & Kovács (2023).
The remainder of the paper is organized as follows. In Section 2, we give a rigorous statement of the partial differential equation we wish to study, together with a derivation of the required matrices |$H$|. In Section 3, we state a natural weak formulation and introduce our fully discrete finite element approximation. We also prove an unconditional stability result for the scheme. Section 4 is devoted to the proof of our main error estimates, which include an |$\mathcal O(h+\varDelta t)$| bound for a discrete |$H^{1}$|–norm, and an |$\mathcal O(h^{2}+\varDelta t)$| |$L^{2}$|–error bound. Finally, in Section 5, we present the results of some numerical simulations.
We end this section with a few comments about notation. We adopt the standard notation for Sobolev spaces, denoting the norm of |$W^{\ell ,p}(I)$| (|$\ell \in{{\mathbb{N}}}_{0}$|, |$p \in [1, \infty ]$|) by |$\|\cdot \|_{\ell ,p}$| and the semi-norm by |$|\cdot |_{\ell ,p}$|. For |$p=2$|, |$W^{\ell ,2}(I)$| will be denoted by |$H^{\ell }(I)$| with the associated norm and semi-norm written as, respectively, |$\|\cdot \|_{\ell }$| and |$|\cdot |_{\ell }$|. The above are naturally extended to vector functions, and we will write |$[W^{\ell ,p}(I)]^{d}$| for a vector function with |$d$| components. In addition, we adopt the standard notation |$W^{\ell ,p}(a,b;B)$| (|$\ell \in{{\mathbb{N}}}$|, |$p \in [1, \infty ]$|, |$(a,b)$| an interval in |${{\mathbb{R}}}$|, |$B$| a Banach space) for time dependent spaces with norm |$\|\cdot \|_{W^{\ell ,p}(a,b;B)}$|. Once again, we write |$H^{\ell }(a,b;B)$| if |$p=2$|. In addition, throughout |$C$| denotes a generic positive constant independent of the mesh parameter |$h$| and the time step size |$\varDelta t$|. At times |$\varepsilon $| will play the role of a (small) positive parameter, with |$C_{\varepsilon }>0$| depending on |$\varepsilon $|, but independent of |$h$| and |$\varDelta t$|. Finally, in this paper we make use of the Einstein summation convention.
2. Mathematical formulation
2.1 Anisotropic curve shortening flow
Let |$\phi \in C^{0}({{\mathbb{R}}}^{d}, {{\mathbb{R}}}_{\geq 0}) \cap C^{4}({{\mathbb{R}}}^{d} \setminus \{0\}, {{\mathbb{R}}}_{>0})$|, as well as
It is not difficult to verify that (2.1) implies that
where |$\phi ^{\prime}$| and |$\phi ^{\prime\prime}$| denote the gradient and the Hessian of |$\phi $|, respectively. In addition, we assume that |$\phi $| is strictly convex in the sense that
Let us consider the anisotropic energy |$E_{\phi }$| defined in terms of |$\phi $| via (1.1) and assume that |$\varGamma =\{ x(\rho ): \rho \in I \}$|. In view of (2.1), we have
Using (2.2) and integration by parts, the first variation of |$E_{\phi }$| in the direction of a vectorfield |$V$| defined on |$\varGamma $| can be derived as
where |$\varkappa =\tau _{s}$| denotes the (principal) curvature vector of |$\varGamma $|. The quantity
can then be viewed as an anisotropic curvature vector. Note that |$\varkappa _{\phi } \cdot \tau = \varkappa \cdot \phi ^{\prime\prime}(\tau ) \tau =0$|, so that |$\varkappa _{\phi } \in \{\tau \}^{\perp }$|. In order to define a gradient flow for |$E_{\phi }$| we observe that only the normal part of a vectorfield |$V$| will contribute to a change in the shape of |$\varGamma $|. This motivates to consider the gradient flow with respect to the inner product
see, e.g. Barrett et al. (2010), where |$P= \text{Id} - \tau \otimes \tau $| denotes the projection onto the normal part of |$\varGamma $|. A family of curves |$(\varGamma (t))_{t \in (0,T)}$| in |${{\mathbb{R}}}^{d}$| then evolves according to anisotropic curve shortening flow provided that
where |${{\mathfrak{m}}}: {{\mathbb{S}}}^{d-1} \to{{\mathbb{R}}}_{>0}$| is a given mobility function in order to allow for a more general setting. For simplicity we extend |${{\mathfrak{m}}}$| 0-homogeneously to |${{\mathbb{R}}}^{d} \setminus \{0\}$| and require for our analysis that |${{\mathfrak{m}}}\in C^{3}({{\mathbb{R}}}^{d} \setminus \{0\})$|. The case |${{\mathfrak{m}}}(\tau )=1$| is most frequently treated in the literature, while the choice |${{\mathfrak{m}}}(\tau ) = \frac 1{\phi (\tau )}$| also has some nice properties and was considered, for example, in Deckelnick et al. (2005), Pozzi (2012) and Deckelnick & Nürnberg (2023b). In what follows we assume a parametric description of |$\varGamma (t)$|, i.e. |$\varGamma (t)=x(I,t)$|, with |$I={{\mathbb{R}}} / {\mathbb{Z}}$|. Then the unit tangent is given by |$\tau = \frac{x_{\rho }}{|x_{\rho }|}$|, while the anisotropic curvature vector is calculated as
Thus, the family |$(\varGamma (t))_{t \in (0,T)}$| evolves according to (2.5) provided that
2.2 DeTurck’s trick for anisotropic curve shortening flow
In what follows we shall make frequent use of the function
It is not difficult to verify that |$\varPhi \in C^{1}({{\mathbb{R}}}^{d}) \cap C^{4}({{\mathbb{R}}}^{d} \setminus \{0\})$| is convex and that
where we think of |$\varPhi ^{\prime\prime\prime}(p)$| as a symmetric trilinear form on |${{\mathbb{R}}}^{d} \times{{\mathbb{R}}}^{d} \times{{\mathbb{R}}}^{d}$|. Moreover, we have the following result.
It is shown in Giga (2006, Remark 1.7.5) that (2.3) implies that |$\varPhi ^{\prime\prime}(p)$| is positive definite for all |$p \in{{\mathbb{R}}}^{d} \setminus \{0\}$|. Hence, the estimates (2.8a) and (2.8b) follow.
As described in the introduction, our aim is to construct positive definite matrices |$H(p) \in{{\mathbb{R}}}^{d \times d}$| such that solutions of the system
solve (2.6). To this end, for any |$p \in{{\mathbb{R}}}^{d} \setminus \{0\}$| we make the ansatz
and where |$\alpha (p) \in{{\mathbb{R}}}_{>0}$| and |$w(p) \in{{\mathbb{R}}}^{d}$| have to be chosen appropriately. Here the two terms |$\text{Id} + \mathfrak t \otimes w$| in the definition of |$H^{-1}$| are guided by the insight that tangential changes to |$x_{t}$| do not change the parameterized flow. Adding the third term |$- w\otimes \mathfrak t$| in (2.10) then ensures positive definiteness. Remarkably, the simple scaling factor |$1/\alpha $| allows enough freedom to guarantee that in the normal directions the correct flow is obtained, i.e. that (2.6) holds.
We will now derive values for |$\alpha (p)$| and |$w(p)$| so that (2.6) holds. If we assume that (2.9) holds, then the ansatz (2.10) yields
Direct calculation shows that
Combining (2.11) and (2.12), on noting |$[\phi ^{\prime}(x_{\rho })]_{\rho } \cdot \tau = 0$|, recall (2.4), and (2.2), yields that
provided that
which can be achieved by setting |$w = \frac 1{\phi (\tau )} P \phi ^{\prime}(x_{\rho }) = \frac{|x_{\rho }|}{\phi (x_{\rho })} \Big(\text{Id} - \frac{x_{\rho }\otimes x_{\rho }}{|x_{\rho }|^{2}}\Big) \phi ^{\prime}(x_{\rho })$|. Hence, with this choice (2.6) will be satisfied if we let
Before we summarize our results, we state some properties of |$H$| that immediately follow from the ansatz (2.10).
In summary, we have shown the following result.
In addition, we can show that (2.9) is parabolic.
The system (2.9) is parabolic in the sense of Petrovsky.
(a) In the isotropic case |$\phi (p) = |p|$| with |${{\mathfrak{m}}} = 1$|, we have |$w(p) = 0$| and |$\alpha (p) = |p|^{2}$|, so that |$H(p)=|p|^{2} \text{Id}$|. Therefore, (2.9) becomes |$| x_{\rho } |^{2} x_{t} = x_{\rho \rho }$|, which is precisely the equation considered in Deckelnick & Dziuk (1995).
(b) We would like to compare (2.16) with the two-dimensional analogue from Deckelnick & Nürnberg (2023b). In fact, for the setting in Deckelnick & Nürnberg (2023b), we have |${{\mathfrak{m}}}(p) = \frac{|p|}{\phi (p)}$| and |$\phi (p)=\gamma (p^{\perp })$|, where |$p^{\perp }= { {p_{1}}\choose{p_{2}}}^{\perp } = {{-p_{2}}\choose{p_{1}}}$| and |$\gamma :{{\mathbb{R}}}^{2} \to{{\mathbb{R}}}_{\geq 0}$| is a normal-dependent anisotropic density function. Hence, we have from (2.17) that |$\alpha (p) = |p|^{2}$| and
3. Fully discrete finite element approximation
The weak formulation corresponding to (2.9) reads as follows. Given |$x_{0}: I \to{{\mathbb{R}}}^{d}$|, find |$x:I \times [0,T] \rightarrow{{\mathbb{R}}}^{d}$| such that |$x(\cdot ,0)= x_{0}$| and, for |$t \in (0,T]$|,
On choosing |$\eta = x_{t}$| in (3.1) we obtain the natural energy estimate
in view of (2.18).
In order to define our finite element approximation, let |$0=q_{0} < q_{1} < \ldots < q_{J-1}< q_{J}=1$| be a decomposition of |$[0,1]$| into intervals |$I_{j}=[q_{j-1},q_{j}]$|. Let |$h_{j}=q_{j} - q_{j-1}$|, as well as |$h=\max _{1 \leq j \leq J} h_{j}$|. We assume that there exists a positive constant |$c$| such that
so that the resulting family of partitions of |$[0,1]$| is quasi-uniform. Within |$I$| we identify |$q_{J}=1$| with |$q_{0}=0$| and define the finite element spaces
Let |$\{\chi _{j}\}_{j=1}^{J}$| denote the standard basis of |$V^{h}$|. For later use, we let |$\pi ^{h}:C^{0}(I)\to V^{h}$| be the standard interpolation operator at the nodes |$\{q_{j}\}_{j=1}^{J}$|, and we use the same notation for the interpolation of vector-valued functions. It is well known that for |$k \in \{ 0,1 \}$|, |$\ell \in \{ 1,2 \}$| and |$p \in [2,\infty ]$| it holds that
In order to discretize in time, let |$t_{m}=m \varDelta t$|, |$m=0,\ldots ,M$|, with the uniform time step |$\varDelta t = \frac TM>0$|. Then our finite element scheme is defined as follows. Let |$x^{0}_{h} = Q_{h} x_{0}$|, where the nonlinear projection |$Q_{h}$| is defined in Lemma 4.1 below. Then, for |$m=0,\ldots ,M-1$|, find |$x^{m+1}_{h} \in \underline{V}^{h}$| such that
We note that in the case |$d=2$|, and for the mobility |${{\mathfrak{m}}}(\tau )=\frac 1{\phi (\tau )}$|, the scheme (3.4), with mass lumping used in the first integral, is identical to the fully discrete approximation (5.4) from Deckelnick & Nürnberg (2023b), recall also Remark 2.5.
We begin by stating an unconditional stability result for the scheme (3.4).
Observe that (3.5) is a fully discrete analogue of (3.2). We note that a discrete analogue of |$\frac{\text{d}}{\text{d}t} \int _{I} \phi (x_{\rho }) \;\text{d}\rho \leq 0$|, in analogy to this property holding for solutions of the continuous problem (3.1), is much harder to prove. For the isotropic case such a discrete analogue can be found in Bänsch et al. (2023, Lemma 4.1.3) for the scheme proposed in Deckelnick & Dziuk (1995) with mass lumping. However, extending these techniques to the anisotropic problem studied here appears to be highly nontrivial. Nevertheless, we remark that in all our numerical experiments, both |$\tfrac 12 \int _{I} \varPhi \big (x^{m}_{h,\rho }\big ) \;\text{d}\rho $| and |$\int _{I} \phi \big (x^{m}_{h,\rho }\big ) \;\text{d}\rho $| are monotonically decreasing.
Our main result is stated in the following theorem. Here, and from now on, for a function |$f \in C([0,T];B)$|, with some Banach space |$B$|, we let |$f^{m}= f(t_{m})$|.
4. Proof of Theorem 3.2
As a crucial ingredient of our error analysis we introduce the following nonlinear Ritz-type projection.
See Appendix A.
Let us define
as well as
where we have used Einstein summation convention and where
We note in view of (2.8b) that |$E^{m}$| behaves like |$|Q_{h} x^{m} - x^{m}_{h} |_{1}^{2}$|, while |$F^{m}$| does not have a sign, but will be controlled by |$E^{m}$| and with the help of Lemma 4.1. Our aim is to obtain the superconvergence bound |$\mathcal O(h^{4}+(\varDelta t)^{2})$| for |$E^{m}$|, which in turn will lead to the optimal |$L^{2}$|-error estimate in (3.9). This procedure can be seen as a nonlinear variant of a technique introduced by Wheeler (1973) for the heat equation.
We shall prove Theorem 3.2 with the help of an induction argument. In particular, we will prove that there exist |$h_{0}>0$|, |$0< \delta \leq 1$| and |$\mu>0$| such that if |$0< h \leq h_{0}$| and |$\varDelta t \leq \delta h$|, then for |$m \in \{ 0,\ldots ,M \}$| the discrete solution |$x^{m}_{h}$| exists and satisfies
Since |$x^{0}_{h} = Q_{h} x^{0}$|, the assertion (4.5) clearly holds for |$m=0$| in view of (3.8) and (4.2), for |$h_{0}$| chosen sufficiently small and for arbitrary |$0 < \delta \leq 1$| and |$\mu>0$|. On assuming that (4.5) holds for a fixed |$m \in \{ 0,\ldots ,M-1 \}$|, we will now show that it also holds for |$m+1$|. Let us define
We infer from (2.8), (2.18) and (2.17) that there exists |$\sigma>0$| such that (2.8) holds for this |$K$|, as well as
Let us abbreviate
We have for any |$z \in [x^{m}_{h,\rho }(\rho ),(Q_{h} x^{m})_{\rho }(\rho )]$|, say |$z=\lambda x^{m}_{h,\rho }(\rho )+(1-\lambda ) (Q_{h} x^{m})_{\rho }(\rho )$|, that
provided that |$h_{0}$| is small enough, where we have used (3.8), (4.5) and (4.2b). Deriving an upper bound for |$z$| in a similar way one obtains that |$[x^{m}_{h,\rho }(\rho ),(Q_{h} x^{m})_{\rho }(\rho )] \subset K$| for all |$\rho \in I$|. Thus, we deduce with the help of (2.8b), (4.2) and (4.5) that
and therefore
provided that |$ h_{0}^{2}\, e^{\mu T} \leq 1$| and |$\delta ^{2}\, e^{\mu T} \leq 1$|. Let us now begin with the induction step.
The bound (4.10) then follows with the help of the inverse estimate (3.3a).
We remark that absolute 1-homogeneity assumption in (2.1) is only used for the proof of Lemma 4.2 via the estimate (4.11). We expect that the proof, and hence the results in this paper, can be extended to positively homogeneous anisotropies, i.e. |$\phi (\lambda p) = \lambda \phi (p)$| for |$p \in{{\mathbb{R}}}^{d}$|, |$\lambda \in{{\mathbb{R}}}_{>0}$|, under a more restrictive time step condition.
For later use we remark that in |$I$| it holds that
provided that |$h_{0}$| and |$\delta $| are sufficiently small. For example, if |$z=(1-\lambda ) x^{m}_{h,\rho } + \lambda x^{m+1}_{h,\rho } \in [x^{m}_{h,\rho },x^{m+1}_{h,\rho }]$|, then we have similarly to (4.8) that
where we have used (3.8), (4.5) and (4.10). The other inclusions can be shown in a similar way, on also making use of (3.7). In particular, we obtain in a similar way as in (4.9) that
Evaluating (3.1) for |$t=t_{m+1}$|, we have
Combining this relation with (3.4), and recalling (4.7), (4.4) and (4.1), we obtain the error equation
Choosing |$\eta _{h}=e^{m+1}_{h}-e^{m}_{h}$| and recalling (4.6), we obtain
The treatment of the second term on the left-hand side of (4.15) is quite involved. That is why we deal with it in the following lemma.
Let us next turn to the terms on the right-hand side of (4.15). We have
Collecting (4.24), (4.25), (4.28) and (4.29) and recalling the definition of |$F^{m}$| yields the result.
Let us now insert the estimates obtained in Lemma 4.4 and Lemma 4.5 into (4.15). After choosing |$\varepsilon $| sufficiently small, we obtain
in view of (4.9) and (4.14). If we choose |$h_{0}$| so small that |$C_{3} \varDelta t \leq \frac 12$| for |$\varDelta t \leq \delta h_{0}$|, then |$0 < (1 - C_{3} \varDelta t)^{-1} \leq 1 + 2 C_{3} \varDelta t$|, and so it follows from (4.5) that
if we choose |$\mu = 6C_{3}$|. This proves the second estimate in (4.5). In order to show the first estimate in (4.5), we observe from (4.7), (3.3a), (4.2b) and (4.9) with |$m$| replaced by |$m+1$|, that
provided that |$h_{0}$| is chosen sufficiently small. Since |$\mu $| and |$\delta $| were chosen independently of |$h$| and |$\varDelta t$|, we have shown (4.5) by induction.
It remains to show that (4.5) implies the desired result (3.9). The second bound in (3.9) follows from (4.2a) and (4.9), since |$x^{m}- x^{m}_{h}=x^{m} - Q_{h} x^{m}+ e^{m}_{h}$|. In order to prove the |$L^{2}$|-error bound in (3.9), we first remark that (4.30) together with (4.5) implies that
Since |$e^{0}_{h}=0$|, we obtain with the help of (4.31) for |$1 \leq \ell \leq M$|
The discrete Gronwall lemma yields that |$\max _{0 \leq m \leq M} \Vert e^{m}_{h} \Vert _{0}^{2} \leq C(h^{4} + (\varDelta t)^{2})$|, so that the |$L^{2}$|–bound in (3.9) now follows again from (4.2a).
5. Numerical results
We implemented the scheme (3.4) within the finite element toolbox Alberta (Schmidt & Siebert, 2005). The systems of nonlinear equations arising at each time level are solved using a Newton iteration, where the solutions to the linear subsystems are found with the help of the sparse factorization package UMFPACK; see Davis (2004). For all our numerical simulations we use a uniform partitioning of |$[0,1]$|, so that |$q_{j} = jh$|, |$j=0,\ldots ,J$|, with |$h = \frac 1J$|. Unless otherwise stated, we use |$J=512$| and |$\varDelta t = 10^{-4}$|. For our numerical simulations we will often be interested in the evolution of the discrete energy |$E_{\phi }(x^{m}_{h}) = \int _{I} \phi \big (x^{m}_{h,\rho }\big ) \;\text{d}\rho $|, recall (1.1). We also consider the ratio
between the longest and shortest element of |$\varGamma ^{m}_{h} = x_{h}^{m}(I)$|, and are often interested in the evolution of this ratio over time. We stress that no redistribution of vertices was necessary during any of our numerical simulations.
Moreover, we will at times be interested in a possible blow-up in curvature. To this end, given |$x^{m}_{h} \in \underline{V}^{h}$|, we introduce the discrete curvature vector |$\kappa ^{m}_{h} \in \underline{V}^{h}$| such that
In practice, we will then monitor the quantity
as an approximation to the maximal value of |$|\varkappa | = \frac{|\tau _{\rho }|}{|x_{\rho }|}$|.
5.1 Convergence experiment
We begin with a convergence experiment in order to confirm our theoretical results from Theorem 3.2. To this end we fix |$\phi (p) = \sqrt{{\widehat{\delta}} ^{2} p_{1}^{2} + p_{2}^{2} + p_{3}^{2}}$| and let |${{\mathfrak{m}}}(\tau ) = \frac 1{\phi (\tau )}$|, and then construct a suitable right-hand side for the related flow
in such a way, that the exact solution is given by the family of self-similarly shrinking ellipses parameterized by
Upon adding the correction term
to the right-hand side of (3.4), we can perform a convergence experiment for our scheme, comparing the obtained discrete solutions with (5.2). The results are displayed in Tables 1 and 2, where we observe optimal convergence error estimates, in line with the results proven in Theorem 3.2. Here we partition the time interval |$[0,T]$|, with |$T=0.45$|, into uniform time steps of size either |$\varDelta t = h$| or |$\varDelta t = h^{2}$|, for |$h = J^{-1} = 2^{-k}$|, |$k=6,\ldots ,11$|.
Errors for the convergence test for (5.2) over the time interval |$[0,0.45]$|, using |$\varDelta t = h$|. We also display the experimental orders of convergence (EOC)
|$J$| . | |$\max _{m=0,\ldots ,M} \left \| y(\cdot ,t_{m}) - x^{m}_{h}\right \|_{0}$| . | EOC . | |$\max _{m=0,\ldots ,M} \left \| y(\cdot ,t_{m}) - x^{m}_{h}\right \|_{1}$| . | EOC . |
---|---|---|---|---|
64 | 4.2446e-02 | — | 2.7051e-01 | — |
128 | 2.4117e-02 | 0.82 | 1.5308e-01 | 0.82 |
256 | 1.3110e-02 | 0.88 | 8.3054e-02 | 0.88 |
512 | 6.8444e-03 | 0.94 | 4.3326e-02 | 0.94 |
1024 | 3.5095e-03 | 0.96 | 2.2206e-02 | 0.96 |
2048 | 1.7757e-03 | 0.98 | 1.1234e-02 | 0.98 |
|$J$| . | |$\max _{m=0,\ldots ,M} \left \| y(\cdot ,t_{m}) - x^{m}_{h}\right \|_{0}$| . | EOC . | |$\max _{m=0,\ldots ,M} \left \| y(\cdot ,t_{m}) - x^{m}_{h}\right \|_{1}$| . | EOC . |
---|---|---|---|---|
64 | 4.2446e-02 | — | 2.7051e-01 | — |
128 | 2.4117e-02 | 0.82 | 1.5308e-01 | 0.82 |
256 | 1.3110e-02 | 0.88 | 8.3054e-02 | 0.88 |
512 | 6.8444e-03 | 0.94 | 4.3326e-02 | 0.94 |
1024 | 3.5095e-03 | 0.96 | 2.2206e-02 | 0.96 |
2048 | 1.7757e-03 | 0.98 | 1.1234e-02 | 0.98 |
Errors for the convergence test for (5.2) over the time interval |$[0,0.45]$|, using |$\varDelta t = h$|. We also display the experimental orders of convergence (EOC)
|$J$| . | |$\max _{m=0,\ldots ,M} \left \| y(\cdot ,t_{m}) - x^{m}_{h}\right \|_{0}$| . | EOC . | |$\max _{m=0,\ldots ,M} \left \| y(\cdot ,t_{m}) - x^{m}_{h}\right \|_{1}$| . | EOC . |
---|---|---|---|---|
64 | 4.2446e-02 | — | 2.7051e-01 | — |
128 | 2.4117e-02 | 0.82 | 1.5308e-01 | 0.82 |
256 | 1.3110e-02 | 0.88 | 8.3054e-02 | 0.88 |
512 | 6.8444e-03 | 0.94 | 4.3326e-02 | 0.94 |
1024 | 3.5095e-03 | 0.96 | 2.2206e-02 | 0.96 |
2048 | 1.7757e-03 | 0.98 | 1.1234e-02 | 0.98 |
|$J$| . | |$\max _{m=0,\ldots ,M} \left \| y(\cdot ,t_{m}) - x^{m}_{h}\right \|_{0}$| . | EOC . | |$\max _{m=0,\ldots ,M} \left \| y(\cdot ,t_{m}) - x^{m}_{h}\right \|_{1}$| . | EOC . |
---|---|---|---|---|
64 | 4.2446e-02 | — | 2.7051e-01 | — |
128 | 2.4117e-02 | 0.82 | 1.5308e-01 | 0.82 |
256 | 1.3110e-02 | 0.88 | 8.3054e-02 | 0.88 |
512 | 6.8444e-03 | 0.94 | 4.3326e-02 | 0.94 |
1024 | 3.5095e-03 | 0.96 | 2.2206e-02 | 0.96 |
2048 | 1.7757e-03 | 0.98 | 1.1234e-02 | 0.98 |
Errors for the convergence test for (5.2) over the time interval |$[0,0.45]$|, using |$\varDelta t = h^{2}$|. We also display the experimental orders of convergence (EOC)
|$J$| . | |$\max _{m=0,\ldots ,M} \left \| y(\cdot ,t_{m}) - x^{m}_{h}\right \|_{0}$| . | EOC . | |$\max _{m=0,\ldots ,M} \left \| y(\cdot ,t_{m}) - x^{m}_{h}\right \|_{1}$| . | EOC . |
---|---|---|---|---|
64 | 5.4585e-04 | — | 1.1777e-01 | — |
128 | 1.3651e-04 | 2.00 | 5.8889e-02 | 1.00 |
256 | 3.4129e-05 | 2.00 | 2.9445e-02 | 1.00 |
512 | 8.5325e-06 | 2.00 | 1.4723e-02 | 1.00 |
1024 | 2.1331e-06 | 2.00 | 7.3614e-03 | 1.00 |
2048 | 5.3328e-07 | 2.00 | 3.6807e-03 | 1.00 |
|$J$| . | |$\max _{m=0,\ldots ,M} \left \| y(\cdot ,t_{m}) - x^{m}_{h}\right \|_{0}$| . | EOC . | |$\max _{m=0,\ldots ,M} \left \| y(\cdot ,t_{m}) - x^{m}_{h}\right \|_{1}$| . | EOC . |
---|---|---|---|---|
64 | 5.4585e-04 | — | 1.1777e-01 | — |
128 | 1.3651e-04 | 2.00 | 5.8889e-02 | 1.00 |
256 | 3.4129e-05 | 2.00 | 2.9445e-02 | 1.00 |
512 | 8.5325e-06 | 2.00 | 1.4723e-02 | 1.00 |
1024 | 2.1331e-06 | 2.00 | 7.3614e-03 | 1.00 |
2048 | 5.3328e-07 | 2.00 | 3.6807e-03 | 1.00 |
Errors for the convergence test for (5.2) over the time interval |$[0,0.45]$|, using |$\varDelta t = h^{2}$|. We also display the experimental orders of convergence (EOC)
|$J$| . | |$\max _{m=0,\ldots ,M} \left \| y(\cdot ,t_{m}) - x^{m}_{h}\right \|_{0}$| . | EOC . | |$\max _{m=0,\ldots ,M} \left \| y(\cdot ,t_{m}) - x^{m}_{h}\right \|_{1}$| . | EOC . |
---|---|---|---|---|
64 | 5.4585e-04 | — | 1.1777e-01 | — |
128 | 1.3651e-04 | 2.00 | 5.8889e-02 | 1.00 |
256 | 3.4129e-05 | 2.00 | 2.9445e-02 | 1.00 |
512 | 8.5325e-06 | 2.00 | 1.4723e-02 | 1.00 |
1024 | 2.1331e-06 | 2.00 | 7.3614e-03 | 1.00 |
2048 | 5.3328e-07 | 2.00 | 3.6807e-03 | 1.00 |
|$J$| . | |$\max _{m=0,\ldots ,M} \left \| y(\cdot ,t_{m}) - x^{m}_{h}\right \|_{0}$| . | EOC . | |$\max _{m=0,\ldots ,M} \left \| y(\cdot ,t_{m}) - x^{m}_{h}\right \|_{1}$| . | EOC . |
---|---|---|---|---|
64 | 5.4585e-04 | — | 1.1777e-01 | — |
128 | 1.3651e-04 | 2.00 | 5.8889e-02 | 1.00 |
256 | 3.4129e-05 | 2.00 | 2.9445e-02 | 1.00 |
512 | 8.5325e-06 | 2.00 | 1.4723e-02 | 1.00 |
1024 | 2.1331e-06 | 2.00 | 7.3614e-03 | 1.00 |
2048 | 5.3328e-07 | 2.00 | 3.6807e-03 | 1.00 |
5.2 Anisotropic curve shortening flow in the plane
We begin this subsection with a short investigation into the tangential motion induced by our scheme (3.4), in comparison to other schemes in the literature. To this end, we repeat the exact same experiment from Barrett et al. (2008, Figs 14 and 15), which compared the numerical methods from Dziuk (1999) and Barrett et al. (2008) for the anisotropic curve shortening flow starting from a unit circle, for the anisotropy |$\phi (p) = \sqrt{ \frac 14 p_{1}^{2} + p_{2}^{2}}$|. In particular, we let |${{\mathfrak{m}}}(\tau )=1$| and use the discretization parameters |$J=128$| and |$\varDelta t = 10^{-3}$|. The evolution obtained from our scheme (3.4) is visually indistinguishable from the results shown in Barrett et al. (2008, Fig. 14), and so we only compare plots of the ratio |${{\mathfrak{r}}}^{m}$| for the three schemes; see Figure 1. Note that the plots shown in Barrett et al. (2008, Fig. 15) compared a |$\phi $|-weighted ratio, for which an equidistribution result can be shown in this particular case for the scheme from Barrett et al. (2008). In contrast, the isotropic ratio |${{\mathfrak{r}}}^{m}$| shown in Fig. 1 increases for all three schemes: the most for the scheme from Dziuk (1999), the least for the scheme from Barrett et al. (2008). For completeness, we remark that repeating the experiment with the finer discretization parameters |$J=512$| and |$\varDelta t = 10^{-4}$| leads to almost unchanged plots for the scheme (3.4) and the scheme from Barrett et al. (2008). However, due to coalescence of mesh points, the scheme from Dziuk (1999) is not able to integrate the evolution until the final time |$T=0.6$|.

Plots of the ratio |${{\mathfrak{r}}}^{m}$| for (3.4) over time. Left: (3.4); middle: Dziuk (1999); right: Barrett et al. (2008).
For the remainder of this subsection we again consider the mobility function |${{\mathfrak{m}}}(\tau ) = \frac 1{\phi (\tau )}$|. We begin with an anisotropy that is not absolutely 1-homogeneous, to underline that our numerical method can also be applied in these situations, recall Remark 4.3. In particular, we consider an anisotropy as in Dziuk (1999, (7.1)) and Barrett et al. (2011, (4.4a)). To this end, let
It is not difficult to verify that this anisotropy satisfies (2.3) if and only if |$|{\widehat{\delta}} | < \frac 1{k^{2} - 1}$|. An example computation for |$k=3$| and |${\widehat{\delta}} =0.124$| can be seen in Fig. 2, where we use as initial data an equidistributed 2:1 ellipse with unit semi major axis. During the evolution the curve approaches as the limiting shape the so-called Wulff shape of the chosen anisotropy. Here we recall from Fonseca & Müller (1991) that the Wulff shape is the solution of the isoperimetric problem for the energy (1.1) in the plane.

Anisotropic curvature flow for an ellipse for the anisotropy (5.3) with |$(k,{\widehat{\delta}} )=(3, 0.124)$|. Solution at times |$t=0,0.05,\ldots ,0.25$| on the left. We also show plots of the discrete energy |$E_{\phi }(x_{h}^{m})$| (middle) and of the ratio |${{\mathfrak{r}}}^{m}$| (right) over time.
With our next numerical experiment we would like to demonstrate that our scheme can also deal with nearly crystalline anisotropies. To this end, we choose as anisotropy the regularized |$\ell ^{1}$|–norm
A simulation starting from a spiral shape is shown in Fig. 3, where we notice the developing facets and corners due to the crystalline nature of the chosen anisotropy. For this experiment we use the finer discretization parameters |$J=1024$| and |$\varDelta t = 5 \times 10^{-8}$|.

Anisotropic curvature flow for a spiral for the anisotropy (5.4). Solution at times |$t=0,0.01,0.015,0.019$|.
5.3 Isotropic curve shortening flow in |${{\mathbb{R}}}^{3}$|
The remainder of our numerical experiments are for the two-codimensional case, and from now on we always choose the constant mobility |${{\mathfrak{m}}}(\tau ) = 1$|. In this particular subsection, we in addition consider the isotropic case, |$\phi (p)=|p|$|. The first experiment is for a trefoil knot in |${{\mathbb{R}}}^{3}$|, and in particular the initial curve is given by
See Figure 4 for the numerical results, which agree very well with the results from Barrett et al. (2010, Fig. 1). Observe that the knot approaches a double covering of a circle within a hyperplane of |${{\mathbb{R}}}^{3}$|.

Isotropic curve shortening flow for the trefoil knot (5.5). On the left, |$x_{h}^{m}$| at times |$t=0,0.5,\ldots ,2,T=2.45$|, with |$x_{h}^{M}$| on the right. Below we show plots of |$E_{\phi }(x_{h}^{m})$|, |${{\mathfrak{r}}}^{m}$| and |$1/K^{m}_{\infty }$| over time.
The next experiment is for two interlocked rings in |${{\mathbb{R}}}^{3}$|, and in particular the initial curve is given by
See Figure 5 for the numerical results, where we can observe a singularity in the flow. The curve forms two loops that shrink, developing into two dove-tails. Here the continuous problem develops a singularity, with the curvature blowing up. From the plot of the inverse of the magnitude of the discrete curvature, (5.1), we can clearly see the effect of the singularity also on the discrete level. Let us emphasize that existence and error estimates for the discrete solution only hold as long as the required smoothness assumptions are satisfied, i.e. before the formation of the singularity. Our numerical scheme, however, simply integrates through the singularity and eventually approaches a circle that shrinks to a point within a hyperplane. Notions of weak solutions that allow an extension of the solution beyond the singularity have been proposed in Angenent (1991) and Deckelnick (1997) for the isotropic case, but it is very difficult to obtain a convergence result after the onset of singularities.

Isotropic curve shortening flow for the two interlocked rings (5.6). On the left, |$x_{h}^{m}$| at times |$t=0,0.25,\ldots ,2,T=2.1$|, in the middle |$x^{M}_{h}$|, with |$x_{h}^{m}$| at time |$t=0.5$| on the right. Below we also show plots of |$E_{\phi }(x_{h}^{m})$|, |${{\mathfrak{r}}}^{m}$| and |$1/K^{m}_{\infty }$| over time.
In our final experiment for the isotropic setting we consider a closed helix in |${{\mathbb{R}}}^{3}$|, as in Barrett et al. (2010, Fig. 2). Here the open helix is defined by
and the initial curve is constructed from (5.7) by connecting |$x_{0}(0)$| and |$x_{0}(1)$| with a polygon that visits the origin and |$(0,0,1)^{T}$|. The evolution of the helix under curve shortening flow can be seen in Fig. 6.

Isotropic curve shortening flow for the helix (5.7). We show |$x_{h}^{m}$| at times |$t=0,0.1,\ldots ,T=0.5$|.
5.4 Anisotropic curve shortening flow in |${{\mathbb{R}}}^{3}$|
Still for the constant mobility |${{\mathfrak{m}}}(\tau ) = 1$|, we now repeat the experiments from the previous subsection, but now for the anisotropy
This means that for the evolving curves it will be energetically favourable to have tangents in the |$y-z$| plane, which should result in the curve itself trying to migrate into that hyperplane. See the results in Figs 7, 8 and 9.

Anisotropic curve shortening flow for the trefoil knot (5.5) for the anisotropy |$\phi (p) = \sqrt{p_{1}^{2} + \tfrac 14 \big(p_{2}^{2} + p_{3}^{2}\big)}$|. On the left, |$x_{h}^{m}$| at times |$t=0,0.5,\ldots ,T=3.5$|, with |$x_{h}^{M}$| on the right.

Anisotropic curve shortening flow for the two interlocked rings (5.6) for the anisotropy |$\phi (p) = \sqrt{p_{1}^{2} + \tfrac 14 \big(p_{2}^{2} + p_{3}^{2}\big)}$|. On the left, |$x_{h}^{m}$| at times |$t=0,0.25,\ldots ,2.5,T=2.75$|, in the middle |$x_{h}^{M}$|, with |$x_{h}^{m}$| at time |$t=0.5$| on the right. Below we also show plots of |$E_{\phi }(x_{h}^{m})$|, |${{\mathfrak{r}}}^{m}$| and |$1/K^{m}_{\infty }$| over time.

Anisotropic curve shortening flow for the helix (5.7) for the anisotropy |$\phi (p) = \sqrt{p_{1}^{2} + \tfrac 14 \big(p_{2}^{2} + p_{3}^{2}\big)}$|. We show |$x_{h}^{m}$| at times |$t=0,0.2,\ldots ,T=1$|.
We observe that in Fig. 8 it is no longer clear that a singularity occurs for the continuous problem. In particular, the evolution of the inverse of the maximal discrete curvature appears to indicate that the curvature for the continuous problem remains bounded.
For the evolution of the initial helix there seems to be little qualitative difference compared to the isotropic results.
However, on repeating the simulation for the stronger anisotropy |$\phi (p) = \sqrt{p_{1}^{2} + 0.01 \big(p_{2}^{2} + p_{3}^{2}\big)}$| yields the evolution in Fig. 10, which clearly shows the effect of the anisotropy. In particular, it can be seen that the curve attempts to avoid any tangent vectors that have a non-zero |$x$|-component, as is to be expected from the chosen anisotropy. We also note the large values of |$K^{m}_{\infty }$| around time |$t=1$|, which may indicate the development of a possible singularity in the flow.

Anisotropic curve shortening flow for the helix (5.7) for the anisotropy |$\phi (p) = \sqrt{p_{1}^{2} + \tfrac 14 \big(p_{2}^{2} + p_{3}^{2}\big)}$|. We show |$x_{h}^{m}$| at times |$t=0,1,\ldots ,T=5$|. Below we also show plots of |$E_{\phi }(x_{h}^{m})$|, |${{\mathfrak{r}}}^{m}$| and |$1/K^{m}_{\infty }$| over time.
We end this section with another numerical simulation for the regularized |$\ell ^{1}$|–norm (5.4). This nearly crystalline anisotropy forces the curve to have tangents aligned with the three coordinate axes. Starting from an initial helix as before, we can observe that effect in Fig. 11. For this experiment we used the smaller time step size |$\varDelta t = 10^{-6}$|.
References
Appendix. A.
Proof of Lemma 4.1: Let |$y$| be given as in Lemma 4.1. We define |$A \in [L^{\infty }(I)]^{d \times d}$| by |$A_{| I_{j}} = \frac{1}{| I_{j}|} \int _{I_{j}} \varPhi ^{\prime\prime}(y_{\rho }) \;\text{d}\rho $|, |$j=1,\ldots J$|. Clearly, it holds that
Moreover, similarly to (2.8a) there exists |$\sigma _{1}> 0$| independent of |$y$| such that
Observing that
we may write (4.1) in the form
On testing the above relation with |$\eta _{h} = Q_{h} y - \pi ^{h} y$|, and recalling (A.2), (A.1) and (3.3b), we infer that
provided that |$0<h \leq h_{1}$| and |$h_{1}>0$| is sufficiently small. Hence, we have that
The estimate (4.2a) now follows with the help of (3.3b), while (3.3b), (3.3a) and (A.3) imply that
In particular, |$|(Q_{h} y)_{\rho }| \geq | y_{\rho } | - | y-Q_{h} y|_{1,\infty } \geq c_{1} -Ch \geq \tfrac 12 c_{1}$| provided that |$h_{1}$| is small enough, and similarly |$| (Q_{h} y)_{\rho } | \leq 2 C_{1}$|. Hence, we have shown (4.2b).
Let us now consider the case where |$y$| depends in addition on time. Differentiating (4.1) with respect to |$t$| yields that
Let |$A \in L^{\infty }(0,T;[L^{\infty }(I)]^{d \times d})$| be defined analogously to above, and |$B \in L^{\infty }(0,T;[L^{\infty }(I)]^{d \times d})$| via |$B_{| I_{j}} =\frac{1}{|I_{j}|} \int _{I_{j}} \varPhi ^{\prime\prime\prime}(y_{\rho }) (y_{t,\rho },\cdot ,\cdot )\;\text{d}\rho $| on |$[0,T]$|. Then we can write (A.4) in the form
where we have used the relation |$\int _{I} B y_{\rho } \cdot \eta _{h,\rho } \;\text{d}\rho = \int _{I} B (\pi ^{h} y)_{\rho } \cdot \eta _{h,\rho } \;\text{d}\rho $|. If we choose |$\eta _{h}= (Q_{h} y - \pi ^{h} y)_{t} \in \underline{V}^{h}$|, use (A.3), the bound |$\| B - \varPhi ^{\prime\prime\prime}(y_{\rho }) (y_{t,\rho },\cdot ,\cdot ) \|_{0,\infty } \leq CC_{1} h$| and argue similarly as above we obtain
from which we infer the first estimate in (4.3). Combining this bound with (3.3a) and (3.3b), we finally obtain
which completes the proof of Lemma 4.1.