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:

(1.1)

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

(1.2)

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

(1.3)

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

(1.4)

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

(1.5)

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

(2.1)

It is not difficult to verify that (2.1) implies that

(2.2)

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

(2.3)

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

(2.4)

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

(2.5)

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.6)

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

(2.7)

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.

 

Lemma 2.1
For every compact set |$K \subset{{\mathbb{R}}}^{d} \setminus \{0\}$| there exists |$\sigma =\sigma (K)>0$| such that
(2.8a)
 
(2.8b)
where |$[p,q]$| denotes the segment between |$p$| and |$q$|⁠.

 

Proof.

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

(2.9)

solve (2.6). To this end, for any |$p \in{{\mathbb{R}}}^{d} \setminus \{0\}$| we make the ansatz

(2.10)

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

(2.11)

Direct calculation shows that

(2.12)

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).

 

Lemma 2.2
For |$\mathfrak t \in{{\mathbb{S}}}^{d-1}$|⁠, |$w \in \{\mathfrak t\}^{\perp }$| and |$\alpha \in{{\mathbb{R}}}_{>0}$| let
(2.13)
Then it holds that
(2.14)
Moreover, the matrix |$H$| is positive definite and satisfies
(2.15)

 

Proof.
Direct calculation shows that
which proves that (2.14) is indeed the inverse of (2.13). Moreover, it holds for |$w \neq 0$| that
where we also used that |$w \cdot \mathfrak t=0$|⁠. This proves the desired result (2.15), since it holds trivially for |$w=0$|⁠.

In summary, we have shown the following result.

 

Theorem 2.3
Let |$\varPhi (p) = \tfrac 12 \phi ^{2}(p)$| and define for |$p \in{{\mathbb{R}}}^{d}\setminus \{0\}$|  
(2.16)
with
(2.17)
If |$x: I \times [0,T] \to{{\mathbb{R}}}^{d}$| satisfies (2.9), then |$x$| is a solution to anisotropic curve shortening flow, (2.6). Furthermore, |$p \mapsto H(p)$| belongs to |$C^{3}({{\mathbb{R}}}^{d} \setminus \{0\}, {{\mathbb{R}}}^{d\times d})$| and |$H(p)$| is positive definite with
(2.18)

In addition, we can show that (2.9) is parabolic.

 

Lemma 2.4

The system (2.9) is parabolic in the sense of Petrovsky.

 

Proof.
On inverting the matrix |$H(x_{\rho })$| we may write (2.9) in the form
Hence, by definition we need to show that the eigenvalues of |$H^{-1}(p) \varPhi ^{\prime\prime}(p)$| have positive real parts for every |$p \in{{\mathbb{R}}}^{d} \setminus \{0\}$|⁠; see, e.g. Eidelman et al. (2004, Definition 1.2). Let us fix |$p \in{{\mathbb{R}}}^{d} \setminus \{0\}$|⁠. A straightforward extension of the proof of (2.15) shows that
Let |$(\lambda ,z) \in{{\mathbb{C}}} \times{{\mathbb{C}}}^{d} \setminus \{0\}$| be an eigenpair of |$H^{-1}(p) \varPhi ^{\prime\prime}(p)$|⁠, i.e. |$\varPhi ^{\prime\prime}(p)z=\lambda H(p)z$|⁠. Then we have
where we have used that |$\varPhi ^{\prime\prime}(p)$| is symmetric and positive definite.

 

Remark 2.5
  • (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

 
On the other hand, (2.14), with |$\alpha =|p|^{2}$|⁠, can be re-written as
(2.19)
where we have used that for |$d=2$| it holds that |$w \otimes w + |w|^{2} \mathfrak t \otimes \mathfrak t = |w|^{2} \text{Id}$|⁠. Observing that |$w \cdot \mathfrak t^{\perp } = - \frac{1}{\gamma (p^{\perp })} \left ( \gamma ^{\prime}(p^{\perp }) \cdot p \right )$| and
we finally obtain from (2.19) that
which is precisely (3.21b) from Deckelnick & Nürnberg (2023b).

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]$|⁠,

(3.1)

On choosing |$\eta = x_{t}$| in (3.1) we obtain the natural energy estimate

(3.2)

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

(3.3a)
(3.3b)

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

(3.4)

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).

 

Theorem 3.1
Any solution of (3.4) satisfies the energy estimate
(3.5)
for |$m=0,\ldots ,M-1$|⁠.

 

Proof.
The convexity of |$\varPhi $| implies that |$\varPhi ^{\prime}(p) \cdot (p-q) \geq \varPhi (p) - \varPhi (q)$| for all |$ p,q \in{{\mathbb{R}}}^{d}$| so that
(3.6)
for all |$\eta _{h} \in \underline{V}^{h}$|⁠. Choosing |$\eta _{h} = x^{m+1}_{h} - x^{m}_{h}$| in (3.4) and applying (3.6) yields the bound (3.5).

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})$|⁠.

 

Theorem 3.2
Suppose that (3.1) has a smooth solution |$x:I \times [0,T] \rightarrow{{\mathbb{R}}}^{d}$|⁠, satisfying
(3.7)
and
(3.8)
for some constants |$c_{0}, C_{0} \in{{\mathbb{R}}}_{>0}$|⁠. Then there exist |$\delta>0$| and |$h_{0}> 0$| such that if |$0<h \leq h_{0}$| and |$\varDelta t \leq \delta h$|⁠, then (3.4) has a unique solution |$(x^{m}_{h})_{m=0,\ldots ,M}$|⁠, and the following error bounds hold:
(3.9)

4. Proof of Theorem 3.2

As a crucial ingredient of our error analysis we introduce the following nonlinear Ritz-type projection.

 

Lemma 4.1
Let |$y \in [W^{2,\infty }(I)]^{d}$| with |$\Vert y \Vert _{2,\infty } \leq C_{1}$| and |$c_{1} \leq | y_{\rho } | \leq C_{1}$| in |$I$| for some |$C_{1},c_{1}>0$|⁠. Then there exists a unique function |$Q_{h} y \in \underline{V}^{h}$| such that
(4.1)
Furthermore, there exist |$h_{1}>0$| and |$C> 0$| depending on |$C_{1}, c_{1}$| and |$\varPhi $| such that
(4.2a)
 
(4.2b)
for |$0<h \leq h_{1}$|⁠. Similarly, if |$y \in W^{1,\infty }(0,T; [W^{2,\infty }(I)]^{d})$| with |${ess\,sup} _{0 < t < T} \Vert (y,y_{t})(\cdot ,t) \Vert _{2,\infty } \leq C_{1}$| and |$| y_{\rho } | \geq c_{1}$| a.e. in |$I \times [0,T]$|⁠, then there exist |$h_{1}>0$| and |$C$| depending on |$C_{1}, c_{1}$| and |$\varPhi $| such that for all |$0<h \leq h_{1}$|  
(4.3)

 

Proof.

See Appendix  A.

Let us define

as well as

where we have used Einstein summation convention and where

(4.4)

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

(4.5)

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

(4.6)

Let us abbreviate

(4.7)

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

(4.8)

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

(4.9)

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.

 

Lemma 4.2
There exists |$\delta>0$| such that for |$\varDelta t \leq \delta h$| there exists a unique element |$x^{m+1}_{h} \in \underline{V}^{h}$|⁠, satisfying (3.4), as well as
(4.10)

 

Proof.
First of all, the existence and uniqueness of |$x^{m+1}_{h}$| can be obtained as in Deckelnick & Nürnberg (2023a, Theorem 2.3), taking advantage of the inequality
(4.11)
where |$\widehat c> 0$| is a positive constant. We remark that the proof of (4.11) from Deckelnick & Nürnberg (2023a, Lemma 2.2) can be easily generalized to |$d\geq 2$|⁠.
Next, on choosing |$\eta _{h}=x^{m+1}_{h} - x^{m}_{h}$| in (3.4), and using (4.6), (4.11) together with the fact that |$x^{m}_{h,\rho } \in K$|⁠, we deduce that
(4.12)
where we have used integration by parts in the last step. Since |$0 \notin [x^{m}_{\rho }(\rho ),x^{m}_{h,\rho }(\rho )] \subset K$| for all |$\rho \in I$|⁠, and since |$\varPhi ^{\prime\prime}$| is |$0$|–homogeneous, recall (2.7), we have from (4.2a) and (4.9) that
Clearly,
so that (4.12) implies that

The bound (4.10) then follows with the help of the inverse estimate (3.3a).

 

Remark 4.3

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

(4.13)

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

(4.14)

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

(4.15)

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.

 

Lemma 4.4
It holds that
provided that |$h_{0}$| and |$\delta $| are small enough.

 

Proof.
Throughout this proof we use the shorthand notations |${\mathfrak{X}} $| for |$Q_{h} x$| and |${\mathfrak{X}} ^{m}$| for |$Q_{h} x^{m}$|⁠. Let
Since |$\varPhi ^{\prime}(p) \cdot p = 2 \varPhi (p)$|⁠, we can rewrite this as
(4.16)
where
Using a Taylor expansion, we obtain
for some |$\theta _{1} \in [{\mathfrak{X}} ^{m}_{\rho },{\mathfrak{X}} ^{m+1}_{\rho }] \subset K$|⁠, recall (4.13). Since |$| {\mathfrak{X}} _{t} |_{1,\infty } \leq C$| by (4.3) and (3.7), we infer on recalling (4.13) and (4.7) that
(4.17)
Similarly, we obtain
where, recalling (4.13),
(4.18)
If we combine the expressions for |$D^{m}_{1}$| and |$D^{m}_{2}$| with |$D^{m}_{3}$| and recall (4.7), we obtain
(4.19)
where we used (2.7) in the last step and defined
In order to estimate |$R^{m}_{3}$| we again use a Taylor expansion and obtain
where |$\theta _{2} \in [{\mathfrak{X}} ^{m}_{\rho },{\mathfrak{X}} ^{m+1}_{\rho }] \subset K$|⁠, |$\theta _{3} \in [{\mathfrak{X}} ^{m}_{\rho },x^{m+1}_{h,\rho }]\subset K$|⁠, recall again (4.13). Then, on recalling (2.7) and (4.3), we have that
(4.20)
where we have observed that for a |$\lambda \in [0,1]$| we can write
Combining (4.19), (4.17), (4.18) and (4.20), and recalling (2.8a) and (4.13), then yields
Inserting this estimate into (4.16), integrating the resulting inequality with respect to |$\rho $| and recalling the definitions of |$D^{m}$| and |$E^{m}$| yields
(4.21)
On recalling (4.10), we have
provided that |$h_{0}$| and |$\delta $| are small enough. Inserting the above relation into (4.21) finishes the proof.

Let us next turn to the terms on the right-hand side of (4.15). We have

 

Lemma 4.5
For every |$\varepsilon>0$| there exists a |$C_{\varepsilon }>0$| such that

 

Proof.
We again use the shorthand notations |${\mathfrak{X}} $| for |$Q_{h} x$| and |${\mathfrak{X}} ^{m}$| for |$Q_{h} x^{m}$|⁠. To begin, we write
Clearly, it holds that
(4.22)
In order to treat |$T_{1,2}$| we write, on recalling (4.4),
Then, with the help of (4.2a), (4.3) and (3.7), we have that
(4.23)
Combining (4.22) and (4.23), we have
(4.24)
Next, on recalling (4.3), we obtain
(4.25)
For the third term on the right-hand side of (4.15), we have
(4.26)
We have from (4.4) and (2.7) that |$\zeta = \varPhi ^{\prime}({\mathfrak{X}} _{\rho }) - \varPhi ^{\prime\prime}(x_{\rho }) {\mathfrak{X}} _{\rho }$|⁠, and so
so that
(4.27)
Combining (4.26) and (4.27) yields, on recalling (4.2) and (4.3), that
(4.28)
Lastly, we can bound
(4.29)

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

(4.30)

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

(4.31)

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

(5.1)

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

(5.2)

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$|⁠.

Table 1

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
644.2446e-022.7051e-01
1282.4117e-020.821.5308e-010.82
2561.3110e-020.888.3054e-020.88
5126.8444e-030.944.3326e-020.94
10243.5095e-030.962.2206e-020.96
20481.7757e-030.981.1234e-020.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
644.2446e-022.7051e-01
1282.4117e-020.821.5308e-010.82
2561.3110e-020.888.3054e-020.88
5126.8444e-030.944.3326e-020.94
10243.5095e-030.962.2206e-020.96
20481.7757e-030.981.1234e-020.98
Table 1

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
644.2446e-022.7051e-01
1282.4117e-020.821.5308e-010.82
2561.3110e-020.888.3054e-020.88
5126.8444e-030.944.3326e-020.94
10243.5095e-030.962.2206e-020.96
20481.7757e-030.981.1234e-020.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
644.2446e-022.7051e-01
1282.4117e-020.821.5308e-010.82
2561.3110e-020.888.3054e-020.88
5126.8444e-030.944.3326e-020.94
10243.5095e-030.962.2206e-020.96
20481.7757e-030.981.1234e-020.98
Table 2

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
645.4585e-041.1777e-01
1281.3651e-042.005.8889e-021.00
2563.4129e-052.002.9445e-021.00
5128.5325e-062.001.4723e-021.00
10242.1331e-062.007.3614e-031.00
20485.3328e-072.003.6807e-031.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
645.4585e-041.1777e-01
1281.3651e-042.005.8889e-021.00
2563.4129e-052.002.9445e-021.00
5128.5325e-062.001.4723e-021.00
10242.1331e-062.007.3614e-031.00
20485.3328e-072.003.6807e-031.00
Table 2

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
645.4585e-041.1777e-01
1281.3651e-042.005.8889e-021.00
2563.4129e-052.002.9445e-021.00
5128.5325e-062.001.4723e-021.00
10242.1331e-062.007.3614e-031.00
20485.3328e-072.003.6807e-031.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
645.4585e-041.1777e-01
1281.3651e-042.005.8889e-021.00
2563.4129e-052.002.9445e-021.00
5128.5325e-062.001.4723e-021.00
10242.1331e-062.007.3614e-031.00
20485.3328e-072.003.6807e-031.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).
Fig. 1.

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

(5.3)

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.
Fig. 2.

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

(5.4)

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$.
Fig. 3.

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

(5.5)

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.
Fig. 4.

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

(5.6)

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.
Fig. 5.

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

(5.7)

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$.
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.
Fig. 7.

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.
Fig. 8.

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$.
Fig. 9.

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.
Fig. 10.

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}$|⁠.

Anisotropic curve shortening flow for the helix (5.7) for the anisotropy (5.4). We show $x_{h}^{m}$ at times $t=0,0.1,0.2,0.3,0.25,T=0.39$.
Fig. 11.

Anisotropic curve shortening flow for the helix (5.7) for the anisotropy (5.4). We show |$x_{h}^{m}$| at times |$t=0,0.1,0.2,0.3,0.25,T=0.39$|⁠.

References

Ambrosio
,
L.
&
Soner
,
H. M.
(
1996
)
Level set approach to mean curvature flow in arbitrary codimension
.
J. Differential Geom.
,
43
,
693
737
.

Angenent
,
S.
(
1991
)
Parabolic equations for curves on surfaces. II. Intersections, blow-up and generalized solutions
.
Ann. of Math.
,
133
,
171
215
.

Bänsch
,
E.
,
Deckelnick
,
K.
,
Garcke
,
H.
&
Pozzi
,
P.
(
2023
)
Interfaces: Modeling, Analysis, Numerics
, vol.
51
of
Oberwolfach Seminars
.
Cham
:
Birkhäuser/Springer
.

Barrett
,
J. W.
,
Garcke
,
H.
&
Nürnberg
,
R.
(
2008
)
Numerical approximation of anisotropic geometric evolution equations in the plane
.
IMA J. Numer. Anal.
,
28
,
292
330
.

Barrett
,
J. W.
,
Garcke
,
H.
&
Nürnberg
,
R.
(
2010
)
Numerical approximation of gradient flows for closed curves in
|${\mathbb{R}}^d$|⁠.
IMA J. Numer. Anal.
,
30
,
4
60
.

Barrett
,
J. W.
,
Garcke
,
H.
&
Nürnberg
,
R.
(
2011
)
The approximation of planar curve evolutions by stable fully implicit finite element schemes that equidistribute
.
Numer. Methods Partial Differential Equations
,
27
,
1
30
.

Barrett
,
J. W.
,
Garcke
,
H.
&
Nürnberg
,
R.
(
2020
)
Parametric finite element approximations of curvature driven interface evolutions
.
Handb. Numer. Anal.
, vol. 21. (A. Bonito & R. H. Nochetto eds)
Amsterdam
:
Elsevier
, pp.
275
423
.

Bellettini
,
G.
&
Paolini
,
M.
(
1996
)
Anisotropic motion by mean curvature in the context of Finsler geometry
.
Hokkaido Math. J.
,
25
,
537
566
.

Beneš
,
M.
,
Kolář
,
M.
&
Ševčovič
,
D.
(
2022
)
Qualitative and numerical aspects of a motion of a family of interacting curves in space
.
SIAM J. Appl. Math.
,
82
,
549
575
.

Binz
,
T.
&
Kovács
,
B.
(
2023
)
A convergent finite element algorithm for mean curvature flow in arbitrary codimension
.
Interfaces Free Bound.
,
25
,
373
400
.

Carlini
,
E.
,
Falcone
,
M.
&
Ferretti
,
R.
(
2007
)
A semi-Lagrangian scheme for the curve shortening flow in codimension-2
.
J. Comput. Phys.
,
225
,
1388
1408
.

Davis
,
T. A.
(
2004
)
Algorithm 832: UMFPACK V4.3—an unsymmetric-pattern multifrontal method
.
ACM Trans. Math. Softw.
,
30
,
196
199
.

Deckelnick
,
K.
(
1997
)
Weak solutions of the curve shortening flow
.
Calc. Var. Partial Differential Equations
,
5
,
489
510
.

Deckelnick
,
K.
&
Dziuk
,
G.
(
1995
)
On the approximation of the curve shortening flow
.
Calculus of Variations, Applications and Computations (Pont-à-Mousson, 1994)
(
C. Bandle, J. Bemelmans, M. Chipot, J. S. J. Paulin & I. Shafrir eds
), vol. 326 of
Pitman Res. Notes Math. Ser.
 
Harlow
:
Longman Sci. Tech.
, pp.
100
108
.

Deckelnick
,
K.
&
Nürnberg
,
R.
(
2023a
)
An unconditionally stable finite element scheme for anisotropic curve shortening flow
.
Arch. Math. (Brno)
,
59
,
263
274
.

Deckelnick
,
K.
&
Nürnberg
,
R.
(
2023b
)
A novel finite element approximation of anisotropic curve shortening flow
.
Interfaces Free Bound.
,
25
,
671
708
.

Deckelnick
,
K.
,
Dziuk
,
G.
&
Elliott
,
C. M.
(
2005
)
Computation of geometric partial differential equations and mean curvature flow
.
Acta Numerica
,
14
,
139
232
.

Dziuk
,
G.
(
1994
)
Convergence of a semi-discrete scheme for the curve shortening flow
.
Math. Models Methods Appl. Sci.
,
04
,
589
606
.

Dziuk
,
G.
(
1999
)
Discrete anisotropic curve shortening flow
.
SIAM J. Numer. Anal.
,
36
,
1808
1830
.

Eidelman
,
S. D.
,
Ivasyshen
,
S. D.
&
Kochubei
,
A. N.
(
2004
)
Analytic Methods in the Theory of Differential and Pseudo-Differential Equations of Parabolic Type
, vol.
152
of
Operator Theory: Advances and Applications
.
Basel
:
Birkhäuser
.

Elliott
,
C. M.
&
Fritz
,
H.
(
2017
)
On approximations of the curve shortening flow and of the mean curvature flow based on the DeTurck trick
.
IMA J. Numer. Anal.
,
37
,
543
603
.

Fonseca
,
I.
&
Müller
,
S.
(
1991
)
A uniqueness proof for the Wulff theorem
.
Proc. Roy. Soc. Edinburgh Sect. A
,
119
,
125
136
.

Giga
,
Y.
(
2006
)
Surface Evolution Equations
, vol.
99
of
Monographs in Mathematics
.
Basel
:
Birkhäuser
.

Gurtin
,
M. E.
(
1993
)
Thermomechanics of Evolving Phase Boundaries in the Plane
.
Oxford Mathematical Monographs
.
New York
:
The Clarendon Press; Oxford University Press
.

Li
,
B.
(
2020
)
Convergence of Dziuk’s linearly implicit parametric finite element method for curve shortening flow
.
SIAM J. Numer. Anal.
,
58
,
2315
2333
.

Ma
,
L.
&
Chen
,
D.
(
2007
)
Curve shortening in a Riemannian manifold
.
Ann. Mat. Pura Appl. (4)
,
186
,
663
684
.

Pozzi
,
P.
(
2007
)
Anisotropic curve shortening flow in higher codimension
.
Math. Methods Appl. Sci.
,
30
,
1243
1281
.

Pozzi
,
P.
(
2012
)
On the gradient flow for the anisotropic area functional
.
Math. Nachr.
,
285
,
707
726
.

Schmidt
,
A.
&
Siebert
,
K. G.
(
2005
)
Design of Adaptive Finite Element Software: The Finite Element Toolbox ALBERTA
, vol.
42
of
Lecture Notes in Computational Science and Engineering
.
Berlin
:
Springer
.

Taylor
,
J. E.
,
Cahn
,
J. W.
&
Handwerker
,
C. A.
(
1992
)
Geometric models of crystal growth
.
Acta Metall. Mater.
,
40
,
1443
1474
.

Wheeler
,
M. F.
(
1973
)
A priori |${L_2}$| error estimates for Galerkin approximations to parabolic partial differential equations
.
SIAM J. Numer. Anal.
,
10
,
723
759
.

Ye
,
C.
&
Cui
,
J.
(
2021
)
Convergence of Dziuk’s fully discrete linearly implicit scheme for curve shortening flow
.
SIAM J. Numer. Anal.
,
59
,
2823
2842
.

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

(A.1)

Moreover, similarly to (2.8a) there exists |$\sigma _{1}> 0$| independent of |$y$| such that

(A.2)

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

(A.3)

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

(A.4)

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.

This is an Open Access article distributed under the terms of the Creative Commons Attribution NonCommercial-NoDerivs licence (http://creativecommons.org/licenses/by-nc-nd/4.0/), which permits non-commercial reproduction and distribution of the work, in any medium, provided the original work is not altered or transformed in any way, and that the work properly cited. For commercial re-use, please contact [email protected] for reprints and translation rights for reprints. All other permissions can be obtained through our RightsLink service via the Permissions link on the article page on our site-for further information please contact [email protected].