1 Introduction

Numerical approximations of curvature flows such as the mean curvature flow and the Gauss curvature flow have been studied intensively during the last 30 years. In many situations the axisymmetry of these geometric flows can be used to reduce the dimension of the governing equations, and so numerical methods have been used frequently in such axisymmetric settings. However, results on the numerical analysis of such schemes so far are rare. In this paper we present parametric finite element approximations for axisymmetric curvature flows, and carefully analyse their properties.

In general, in curvature driven evolution equations the normal velocity of a hypersurface in \({{\mathbb {R}}}^3\) is given by an expression involving the mean and/or the Gauss curvature of the surface. Evolving surfaces are of interest in geometry, and they can appear in application areas such as materials science, for example as grain boundaries. In addition, evolution laws involving the curvature of the surface arise in situations, where surface quantities are coupled to the surrounding volume by additional fields, which for example arises in the evolution of phase boundaries or in two-phase flow. In any case solving the evolution law for the surface with a stable discretization of curvature is a corner stone of a reliable and efficient numerical method.

Approaches to solve surface evolution equations numerically involve different descriptions of the evolving surface. Traditionally level set methods, phase field methods or parametric front tracking methods have been used. For example, parametric finite element approximations of curvature flows have been considered in [7, 20, 23, 35]. We refer to the review paper [18], and the references therein, for further information on numerical methods for general geometric evolution equations.

In this paper we aim to numerically compute a family of hypersurfaces \(({\mathscr {S}}(t))_{t \ge 0}\)\(\subset {{\mathbb {R}}}^3\), which we later assume to be axisymmetric, and which fulfills a geometric evolution law involving its principal curvatures. We will focus on the mean curvature flow, which for \({\mathscr {S}}(t)\) is given by the evolution law

$$\begin{aligned} {\mathscr {V}}_{{\mathscr {S}}} = k_m\quad \text { on } {\mathscr {S}}(t), \end{aligned}$$
(1.1)

and which is the \(L^2\)-gradient flow for \({\mathscr {H}}^2({\mathscr {S}}(t))\), since

$$\begin{aligned} \frac{\mathrm{d}}{{\mathrm{d}}t}\, {\mathscr {H}}^2({\mathscr {S}}(t)) = - \int _{{\mathscr {S}}(t)} k_m\,{\mathscr {V}}_{{\mathscr {S}}} \;{\mathrm{d}}{{\mathscr {H}}}^{2} = - \int _{{\mathscr {S}}(t)} ({\mathscr {V}}_{{\mathscr {S}}})^2 \;{\mathrm{d}}{{\mathscr {H}}}^{2} \end{aligned}$$

for surfaces without boundary. Here \({\mathscr {V}}_{{\mathscr {S}}}\) denotes the normal velocity of \({\mathscr {S}}(t)\) in the direction of the normal \(\vec {\mathrm{n}}_{{\mathscr {S}}}\). Moreover, \(k_m\) is the mean curvature of \({\mathscr {S}}(t)\), i.e. the sum of the principal curvatures of \({\mathscr {S}}(t)\), see [37] for an introduction to the mean curvature flow.

We also consider the nonlinear mean curvature flow

$$\begin{aligned} {\mathscr {V}}_{{\mathscr {S}}} = f(k_m) \quad \text {on }\ {\mathscr {S}}(t), \end{aligned}$$
(1.2)

where \(f{:}\,(a,b)\rightarrow {\mathbb {R}}\) with \(-\,\infty \le a<b\le \infty \), is a strictly monotonically increasing continuous function, as well as the volume preserving variant

$$\begin{aligned} {\mathscr {V}}_{{\mathscr {S}}} = f(k_m) - \frac{\int _{{\mathscr {S}}} f(k_m) \;{\mathrm{d}}{{\mathscr {H}}}^{2}}{\int _{{\mathscr {S}}} 1 \;{\mathrm{d}}{{\mathscr {H}}}^{2}} \quad \text {on }\ {\mathscr {S}}(t). \end{aligned}$$
(1.3)

Possible choices for f are

$$\begin{aligned} f(r) =|r|^{\beta -1}r, \quad \beta \in {\mathbb {R}}_{>0}, \end{aligned}$$
(1.4a)

or

$$\begin{aligned} f(r) = - r^{-1} \end{aligned}$$
(1.4b)

for the inverse mean curvature flow. These two choices have applications for example in image processing or in general relativity, see [33, 42] and the references therein. Of course, (1.2) with

$$\begin{aligned} f(r) = r \end{aligned}$$
(1.4c)

collapses to (1.1).

If \(\varOmega (t)\) denotes the region enclosed by \({\mathscr {S}}(t)\), i.e. \({\mathscr {S}}(t) = \partial \varOmega (t)\), then the flow (1.3) is such that

$$\begin{aligned} \frac{\mathrm{d}}{{\mathrm{d}}t}\,{\mathscr {L}}^3(\varOmega (t)) = \int _{{\mathscr {S}}(t)} {\mathscr {V}}_{{\mathscr {S}}} \;{\mathrm{d}}{{\mathscr {H}}}^{2} = 0, \end{aligned}$$
(1.5)

where here we assume that \(\vec {\mathrm{n}}_{{\mathscr {S}}}\) is the outer normal to \(\varOmega (t)\) on \({\mathscr {S}}(t)\). This justifies the expression volume preserving flow. These flows are of interest in geometry and we refer to [4, 14, 29, 31] for more information.

More generally, we can also consider flows of the form

$$\begin{aligned} {\mathscr {V}}_{{\mathscr {S}}} = F(k_m, k_g) = F(k_1 + k_2, k_1\,k_2) \quad \text {on }\ {\mathscr {S}}(t), \end{aligned}$$
(1.6)

where \(k_g = k_1\,k_2\) denotes the Gauss curvature of \({\mathscr {S}}(t)\), with \(k_1\) and \(k_2\) the two principal curvatures. Of course, (1.6) with \(F(r,s) = f(r)\) reduces to (1.2). On the other hand, the choice \(F(r,s) = -s\), for closed surfaces, leads to the Gauss curvature flow

$$\begin{aligned} {\mathscr {V}}_{{\mathscr {S}}} = - k_g \quad \text {on }\ {\mathscr {S}}(t), \end{aligned}$$
(1.7)

see e.g. [8, (1.14)], where in (1.7) we again assume that \(\vec {\mathrm{n}}_{{\mathscr {S}}}\) is the outer normal to \(\varOmega (t)\) on \({\mathscr {S}}(t)\). Such flows have found considerable interest in geometry recently and we refer to [27, 40, 41, 43] for more information. One reason why the Gauss curvature flow is of particular interest, is because this flow allows to study the fate of the rolling stones, see [2].

In this paper, we consider the case that \({\mathscr {S}}(t)\) is an axisymmetric surface, that is rotationally symmetric with respect to the \(x_2\)-axis. We further assume that \({\mathscr {S}}(t)\) is made up of a single connected component, with or without boundary. Clearly, in the latter case the boundary \(\partial {\mathscr {S}}(t)\) of \({\mathscr {S}}(t)\) consists of either one or two circles that each lie within a hyperplane that is parallel to the \(x_1-x_3\)-plane. For the evolving family of surfaces we allow for the following types of boundary conditions. A boundary circle may assumed to be fixed, it may be allowed to move vertically along the boundary of a fixed infinite cylinder that is aligned with the axis of rotation, or it may be allowed to expand and shrink within a hyperplane that is parallel to the \(x_1-x_3\)-plane. Depending on the postulated free energy, certain angle conditions will arise where \({\mathscr {S}}(t)\) meets the external boundary. If the free energy is just surface area, \({\mathscr {H}}^2({\mathscr {S}}(t))\), then a \(90^\circ \) degree contact angle condition arises. We refer to Sect. 2 below for further details, in particular with regard to more general contact angles.

The dimensionally reduced formulation has several severe advantages both analytically as well as numerically. In analysis it has been used for example to study the onset of singularities, see [22, 32, 36, 38, 40] and other singularity formation mechanisms, see [12, 13]. Numerically it leads to equations which are far easier to solve and at the same time problems with the mesh topology do not occur. Therefore, axisymmetric settings have been frequently used for numerical computations of surface evolutions. For example, graph formulations for axisymmetric geometric evolution laws have been considered in [15, 17, 19], while a finite difference approximation of a parametric description for the evolution of general axisymmetric surfaces has been studied in [39]. Hence the latter is closely related to the presented work, although we stress that it does not contain any numerical analysis. Moreover, also more complex problems such as for example two phase flows or biomembranes, in which also curvature effects play a role, have been treated in an axially symmetric setting. We refer to [16, 26, 30, 44, 45], and we expect that our approach will have an impact on such more complex evolutions as well. In terms of the numerical analysis for the approximation of axisymmetric surface evolutions only very few results have appeared in the literature so far, see e.g. [17, 19] in the context of a graph formulation for the higher order curvature flows surface diffusion and Willmore flow, respectively. To the best of our knowledge, our paper contains the first stability results for fully discrete approximations of axisymmetric mean curvature flow. In addition, we consider the numerical analysis of approximations for axisymmetric higher order flows, such as surface diffusion and Willmore flow, in the recently appeared article [11].

The present authors in the last 10 years introduced parametric finite element methods for geometric evolution equations which have the property that the mesh generically behaves well during the evolution. We also refer to the recent work [23] for a method which also leads to good meshes. This is an advantage compared to earlier front tracking approaches in which often the meshes degenerated during the evolution such that the computations had to be stopped. In a series of papers [5,6,7,8,9], we were able to analyze mesh properties and showed stability results. In particular, in two dimensions a semi-discrete version of the method led to equidistribution of mesh points. In this paper we introduce a parametric finite element method for the axisymmetric formulations of the surface evolution equations discussed above relying on ideas of our earlier work. However, a lot of new techniques have to be introduced stemming partly from the fact that close to the axis of rotation the equations, depending on the formulation, become either singular or degenerate, and partly because one has to decide how to deal with the equidistribution property. We will discuss several ways to handle these issues and will show stability, equidistribution, existence and uniqueness results for the new schemes.

This paper is organised as follows. In Sect. 2 we introduce several weak formulations which will be crucial for the parametric finite element approximations introduced later. In Sect. 3 we derive semidiscrete, i.e. continuous-in-time discrete in space discretizations, and discuss stability and equidistribution properties. Section 4 is devoted to the fully discrete schemes that turn out to be most practical. For the linear schemes we prove existence and uniqueness results, while a stability result can be shown for a mildly nonlinear discretization. Additional fully discrete variants, and some of their properties, are summarized in “Appendix C”. Finally, we present several numerical results demonstrating that the majority of the schemes lead to efficient, reliable results for mean curvature flow as well as for fully nonlinear curvature flows including its mass preserving variants.

Fig. 1
figure 1

Sketch of \(\varGamma \) and \({\mathscr {S}}\), as well as the unit vectors \(\vec {e}_1\), \(\vec {e}_2\) and \(\vec {e}_3\)

2 Weak formulations

Let \({{\mathbb {R}}} \slash {{\mathbb {Z}}}\) be the periodic interval [0, 1], and set

$$\begin{aligned} I = {{\mathbb {R}}} \slash {{\mathbb {Z}}}, \text { with } \partial I = \emptyset ,\quad \text {or}\quad I = (0,1), \text { with } \partial I = \{0,1\}. \end{aligned}$$

We consider the axisymmetric situation, where \(\vec {x}(t){:}\,\overline{I} \rightarrow {{\mathbb {R}}}^2\) is a parameterization of \(\varGamma (t)\). Throughout \(\varGamma (t)\) represents the generating curve of a surface \({\mathscr {S}}(t)\) that is axisymmetric with respect to the \(x_2\)-axis, see Fig. 1. In particular, on defining

$$\begin{aligned} \vec {\varPi }_3^3(r, z, \theta ) = (r \cos \theta , z, r \sin \theta )^T \quad \text {for}\quad r\in {{\mathbb {R}}}_{\ge 0},\ z \in {{\mathbb {R}}},\ \theta \in [0,2\,\pi ] \end{aligned}$$

and

$$\begin{aligned} \varPi _2^3(r, z) = \{\vec {\varPi }_3^3(r, z, \theta ){:}\,\theta \in [0,2\,\pi )\}, \end{aligned}$$

we have that

$$\begin{aligned} {\mathscr {S}}(t) = \bigcup _{(r,z)^T \in \varGamma (t)} \varPi _2^3(r, z) = \bigcup _{\rho \in {{\overline{I}}}} \varPi _2^3(\vec {x}(\rho ,t)). \end{aligned}$$
(2.1)

Here we allow \(\varGamma (t)\) to be either a closed curve, parameterized over \({{\mathbb {R}}} \slash {{\mathbb {Z}}}\), which corresponds to \({\mathscr {S}}(t)\) being a genus-1 surface without boundary. Or \(\varGamma (t)\) may be an open curve, parameterized over [0, 1]. Then \(\varGamma (t)\) has two endpoints, and each endpoint can either correspond to an interior point of \({\mathscr {S}}(t)\), or to a boundary circle of \({\mathscr {S}}(t)\). Endpoints of \(\varGamma (t)\) that correspond to an interior point of the surface \({\mathscr {S}}(t)\) are attached to the \(x_2\)-axis, on which they can freely move up and down. For example, if both endpoints of \(\varGamma (t)\) are attached to the \(x_2\)-axis, then \({\mathscr {S}}(t)\) is a genus-0 surface without boundary. If only one end of \(\varGamma (t)\) is attached to the \(x_2\)-axis, then \({\mathscr {S}}(t)\) is an open surface with boundary, where the boundary consists of a single connected component. If no endpoint of \(\varGamma (t)\) is attached to the \(x_2\)-axis, then \({\mathscr {S}}(t)\) is an open surface with boundary, where the boundary consists of two connected components.

Table 1 The different types of boundary nodes enforced by (2.2b)–(2.2d)

In particular, we always assume that, for all \(t \in [0,T]\),

$$\begin{aligned} \vec {x}(\rho ,t) \cdot \,\vec {e}_1&> 0 \quad \forall \ \rho \in {{\overline{I}}}{\setminus } \partial _0 I, \end{aligned}$$
(2.2a)
$$\begin{aligned} \vec {x}(\rho ,t) \cdot \,\vec {e}_1&= 0 \quad \forall \ \rho \in \partial _0 I, \end{aligned}$$
(2.2b)
$$\begin{aligned} \vec {x}_t(\rho ,t) \cdot \,\vec {e}_i&= 0 \quad \forall \ \rho \in \partial _i I, \ i =1,2, \end{aligned}$$
(2.2c)
$$\begin{aligned} \vec {x}_t(\rho ,t)&= \vec {0} \quad \forall \ \rho \in \partial _D I, \end{aligned}$$
(2.2d)

where \(\partial _D I \cup \bigcup _{i=0}^2 \partial _i I = \partial I\) is a disjoint partitioning of \(\partial I\), with \(\partial _0 I\) denoting the subset of boundary points of I that correspond to endpoints of \(\varGamma (t)\) attached to the \(x_2\)-axis. Moreover, \(\partial _D I \cup \bigcup _{i=1}^2 \partial _i I\) denotes the subset of boundary points of I that model components of the boundary of \({\mathscr {S}}(t)\). Here endpoints in \(\partial _D I\) correspond to fixed boundary circles of \({\mathscr {S}}(t)\), that lie within a hyperplane parallel to the \(x_1-x_3\)-plane \({{\mathbb {R}}}\times \{0\} \times {{\mathbb {R}}}\). Endpoints in \(\partial _1 I\) correspond to boundary circles of \({\mathscr {S}}(t)\) that can move freely along the boundary of an infinite cylinder that is aligned with the axis of rotation. Endpoints in \(\partial _2 I\) correspond to boundary circles of \({\mathscr {S}}(t)\) that can expand/shrink freely within a hyperplane parallel to the \(x_1-x_3\)-plane \({{\mathbb {R}}}\times \{0\} \times {{\mathbb {R}}}\). See Table 1 for a visualization of the different types of boundary nodes.

On assuming that

$$\begin{aligned} |\vec {x}_\rho | \ge c_0 > 0 \quad \forall \ \rho \in {\overline{I}}, \end{aligned}$$
(2.3)

we introduce the arclength s of the curve, i.e. \(\partial _s = |\vec {x}_\rho |^{-1}\,\partial _\rho \), and set

$$\begin{aligned} \vec {\tau }(\rho ,t) = \vec {x}_s(\rho ,t) = \frac{\vec {x}_\rho (\rho ,t)}{|\vec {x}_\rho (\rho ,t)|} \quad \text{ and } \quad \vec {\nu }(\rho ,t) = -[\vec {\tau }(\rho ,t)]^\perp , \end{aligned}$$
(2.4)

where \((\cdot )^\perp \) denotes a clockwise rotation by \(\frac{\pi }{2}\).

On recalling (2.1), we observe that the normal \(\vec {\mathrm{n}}_{{\mathscr {S}}}\) on \({\mathscr {S}}(t)\) is given by

$$\begin{aligned} \vec {\mathrm{n}}_{{\mathscr {S}}}(\vec {\varPi }_3^3(\vec {x}(\rho ,t),\theta )) = \begin{pmatrix} (\vec {\nu }(\rho ,t)\cdot \,\vec {e}_1)\,\cos \theta \\ \vec {\nu }(\rho ,t)\cdot \,\vec {e}_2 \\ (\vec {\nu }(\rho ,t)\cdot \,\vec {e}_1)\,\sin \theta \end{pmatrix} \quad \text {for}\quad \rho \in {\overline{I}},\ t \in [0,T],\ \theta \in [0,2\,\pi ). \end{aligned}$$
(2.5)

Similarly, the normal velocity \({\mathscr {V}}_{{\mathscr {S}}}\) of \({\mathscr {S}}(t)\) in the direction \(\vec {\mathrm{n}}_{{\mathscr {S}}}\) is given by

$$\begin{aligned} {\mathscr {V}}_{{\mathscr {S}}} = \vec {x}_t(\rho ,t)\cdot \,\vec {\nu }(\rho ,t) \quad \text {on } \varPi _2^3(\vec {x}(\rho ,t)) \subset {\mathscr {S}}(t), \quad \forall \ \rho \in {\overline{I}},\ t \in [0,T]. \end{aligned}$$
(2.6)

For the curvature \(\varkappa \) of \(\varGamma (t)\) it holds that

$$\begin{aligned} \varkappa \,\vec {\nu } = \vec {\varkappa } = \vec {\tau }_s = \frac{1}{|\vec {x}_\rho |} \left[ \frac{\vec {x}_\rho }{|\vec {x}_\rho |} \right] _\rho . \end{aligned}$$
(2.7)

An important role in this paper is played by the surface area of the surface \({\mathscr {S}}(t)\), which is equal to

$$\begin{aligned} {\mathscr {H}}^2({\mathscr {S}}(t)) = A(\vec {x}(t)) = 2\,\pi \,\int _I \vec {x}(\rho ,t)\cdot \,\vec {e}_1\,|\vec {x}_\rho (\rho ,t)| \;{\mathrm{d}}\rho . \end{aligned}$$
(2.8)

Often the surface area, \(A(\vec {x}(t))\), will play the role of the free energy in our paper. But for an open surface \({\mathscr {S}}(t)\), with boundary \(\partial {\mathscr {S}}(t)\), we consider contact energy contributions which are discussed in [24], see also [9, (2.21)]. In the axisymmetric setting the relevant energy is given by

$$\begin{aligned} E(\vec {x}(t)) = A(\vec {x}(t)) + 2\,\pi \,\sum _{p\in \partial _1 I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\,(\vec {x}(p,t)\cdot \,\vec {e}_1)\,\vec {x}(p,t)\cdot \,\vec {e}_2 + \pi \,\sum _{p\in \partial _2 I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\,(\vec {x}(p,t)\cdot \,\vec {e}_1)^2, \end{aligned}$$
(2.9)

where we recall from (2.2c) that, for \(i=1,2\), either \(\partial _i I = \emptyset \), \(\{0\}\), \(\{1\}\) or \(\{0,1\}\). In the above \({\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)} \in {{\mathbb {R}}}\), for \(p\in \{0,1\}\), are given constants. Here \({\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\), for \(p \in \partial _1 I\), denotes the change in contact energy density in the direction of \(-\,\vec {e}_2\), that the two phases separated by the interface \({\mathscr {S}}(t)\) have with the infinite cylinder at the boundary circle of \({\mathscr {S}}(t)\) represented by \(\vec {x}(p,t)\). Similarly, \({\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\), for \(p \in \partial _2 I\), denotes the change in contact energy density in the direction of \(-\,\vec {e}_1\), that the two phases separated by the interface \({\mathscr {S}}(t)\) have with the hyperplane \({{\mathbb {R}}}\times \{0\}\times {{\mathbb {R}}}\) at the boundary circle of \({\mathscr {S}}(t)\) represented by \(\vec {x}(p,t)\). These changes in contact energy lead to the contact angle conditions

$$\begin{aligned} (-1)^p\,\vec {\tau }(p,t)\cdot \,\vec {e}_2&= {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)} \qquad p \in \partial _1 I, \end{aligned}$$
(2.10a)
$$\begin{aligned} (-1)^p\,\vec {\tau }(p,t)\cdot \,\vec {e}_1&= {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)} \qquad p \in \partial _2 I, \end{aligned}$$
(2.10b)

for all \(t \in (0,T]\). In most cases, the contact energies are assumed to be the same, so that \({\widehat{\rho }}_{\partial {\mathscr {S}}}^{(0)}={\widehat{\rho }}_{\partial {\mathscr {S}}}^{(1)}=0\), which leads to \(90^\circ \) contact angle conditions in (2.10a, b), and means that (2.9) collapses to (2.8). See [9] for more details on contact angles and contact energies. We note that a necessary condition to admit a solution to (2.10a) or to (2.10b) is that

$$\begin{aligned} |{\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}| \le 1 \quad p \in \{0,1\}. \end{aligned}$$
(2.11)

In addition, we observe that the energy (2.9) is not bounded from below if \({\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)} \not =0 \) for \(p\in \partial _1 I\) or if \({\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}<0\) for \(p\in \partial _2 I\).

For later use we note that

$$\begin{aligned} \frac{\mathrm{d}}{{\mathrm{d}}t}\, E(\vec {x}(t))&= 2\,\pi \,\int _I \left[ \vec {x}_t\cdot \,\vec {e}_1 + \vec {x}\cdot \,\vec {e}_1\,\frac{(\vec {x}_t)_\rho \cdot \,\vec {x}_\rho }{|\vec {x}_\rho |^2} \right] |\vec {x}_\rho | \;{\mathrm{d}}\rho \nonumber \\&\quad +\,2\,\pi \,\sum _{p\in \partial _1 I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\left[ (\vec {x}_t(p,t)\cdot \,\vec {e}_1)\,\vec {x}(p,t)\cdot \,\vec {e}_2 + (\vec {x}(p,t)\cdot \,\vec {e}_1)\,\vec {x}_t(p,t)\cdot \,\vec {e}_2 \right] \nonumber \\&\quad +\,2\,\pi \,\sum _{p\in \partial _2 I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\,(\vec {x}(p,t)\cdot \,\vec {e}_1)\,\vec {x}_t(p,t)\cdot \,\vec {e}_1. \end{aligned}$$
(2.12)

Moreover, we recall that expressions for the mean curvature and the Gauss curvature of \({\mathscr {S}}(t)\) are given by

$$\begin{aligned} \varkappa _{{\mathscr {S}}} = \varkappa - \frac{\vec {\nu }\cdot \,\vec {e}_1}{\vec {x}\cdot \,\vec {e}_1} \quad \text {and}\quad {{\mathscr {K}}}_{{\mathscr {S}}} = - \varkappa \,\frac{\vec {\nu }\cdot \,\vec {e}_1}{\vec {x}\cdot \,\vec {e}_1} \quad \text {on }\ \overline{I}, \end{aligned}$$
(2.13)

respectively; see e.g. [16, (6)]. More precisely, if \(k_m\) and \(k_g\) denote the mean and Gauss curvatures of \({\mathscr {S}}(t)\), then

$$\begin{aligned} k_m = \varkappa _{{\mathscr {S}}}(\rho ,t) \quad \ \text { and }\ k_g = {{\mathscr {K}}}_{{\mathscr {S}}}(\rho ,t) \quad \text {on } \varPi _2^3(\vec {x}(\rho ,t)) \subset {\mathscr {S}}(t), \quad \forall \ \rho \in {\overline{I}},\ t \in [0,T]. \end{aligned}$$

In the literature, the two terms making up \(\varkappa _{{\mathscr {S}}}\) in (2.13) are often referred to as in-plane and azimuthal curvatures, respectively, with their sum being equal to the mean curvature. We note that combining (2.13) and (2.7) yields that

$$\begin{aligned} \varkappa _{{\mathscr {S}}}\,\vec {\nu } = \vec {x}_{ss} - \frac{\vec {\nu }\cdot \,\vec {e}_1}{\vec {x}\cdot \,\vec {e}_1}\,\vec {\nu } = \frac{1}{|\vec {x}_\rho |} \left[ \frac{\vec {x}_\rho }{|\vec {x}_\rho |} \right] _\rho - \frac{\vec {\nu }\cdot \,\vec {e}_1}{\vec {x}\cdot \,\vec {e}_1}\,\vec {\nu }. \end{aligned}$$
(2.14)

It follows from (2.14) and (2.4) that

$$\begin{aligned} \vec {x}\cdot \,\vec {e}_1\,\varkappa _{{\mathscr {S}}}\,\vec {\nu } = (\vec {x}\cdot \,\vec {e}_1)\,\vec {x}_{ss} - (\vec {\nu }\cdot \,\vec {e}_1)\,\vec {\nu } = ((\vec {x}\cdot \,\vec {e}_1)\,\vec {x}_s)_s - \vec {e}_1. \end{aligned}$$
(2.15)

A weak formulation of (2.15) will form the basis of our stable approximations for mean curvature flow and surface diffusion. Clearly, for a smooth surface with bounded mean curvature it follows from (2.14) that

$$\begin{aligned} \vec {\nu }(\rho ,t) \cdot \,\vec {e}_1 = 0 \qquad \forall \ \rho \in \partial _0 I,\quad \forall \ t\in [0,T], \end{aligned}$$
(2.16)

which, on recalling (2.4), is clearly equivalent to

$$\begin{aligned} \vec {x}_\rho (\rho ,t) \cdot \,\vec {e}_2 = 0 \qquad \forall \ \rho \in \partial _0 I,\quad \forall \ t\in [0,T]. \end{aligned}$$
(2.17)

A precise derivation of (2.17) in the context of a weak formulation of (2.14) will be given in the “Appendix A”.

2.1 Mean curvature flow

In terms of the axisymmetric description of \({\mathscr {S}}(t)\), the evolution law (1.1) can be written as

$$\begin{aligned} \vec {x}_t\cdot \,\vec {\nu } = \varkappa _{{\mathscr {S}}} = \varkappa - \frac{\vec {\nu }\cdot \,\vec {e}_1}{\vec {x}\cdot \,\vec {e}_1} \qquad \text {on } I, \end{aligned}$$
(2.18a)

with, on recalling (2.2b–d),

$$\begin{aligned} \vec {x}_t(\rho ,t)\cdot \,\vec {e}_1&= 0 \quad \forall \ \rho \in \partial _0 I, \quad \vec {x}_t(\rho ,t)\cdot \,\vec {e}_i = 0 \quad \forall \ \rho \in \partial _i I,\ i = 1,2, \nonumber \\ \vec {x}_t(\rho ,t)&= \vec {0} \quad \forall \ \rho \in \partial _D I, \qquad \forall \ t\in [0,T], \end{aligned}$$
(2.18b)

as well as (2.17) and (2.10a, b).

Let

$$\begin{aligned} \underline{V}_{\partial _0}&= \left\{ \vec {\eta } \in [H^1(I)]^2{:}\,\vec {\eta }(\rho )\cdot \,\vec {e}_1 = 0 \quad \forall \ \rho \in \partial _0 I\right\} , \nonumber \\ \underline{V}_\partial&= \left\{ \eta \in \underline{V}_{\partial _0}{:}\,\vec {\eta }(\rho )\cdot \,\vec {e}_i = 0 \quad \forall \ \rho \in \partial _i I,\ i=1,2, \quad \vec {\eta }(\rho ) = \vec {0} \quad \forall \ \rho \in \partial _D I\right\} . \end{aligned}$$

Then we consider the following weak formulation of (2.18a, b), on recalling (2.7).

\(({\mathscr {A}})\) Let \(\vec {x}(0) \in \underline{V}_{\partial _0}\). For \(t \in (0,T]\) find \(\vec {x}(t) \in [H^1(I)]^2\), with \(\vec {x}_t(t) \in \underline{V}_\partial \), and \(\varkappa (t)\in L^2(I)\) such that

$$\begin{aligned}&\int _I \vec {x}_t\cdot \,\vec {\nu }\,\chi \,|\vec {x}_\rho |\;{\mathrm{d}}\rho = \int _I \left( \varkappa - \frac{\vec {\nu }\cdot \,\vec {e}_1}{\vec {x}\cdot \,\vec {e}_1} \right) \chi \,|\vec {x}_\rho | \;{\mathrm{d}}\rho \quad \forall \ \chi \in L^2(I), \end{aligned}$$
(2.19a)
$$\begin{aligned}&\int _I \varkappa \,\vec {\nu }\cdot \,\vec {\eta }\, |\vec {x}_\rho | \;{\mathrm{d}}\rho + \int _I (\vec {x}_\rho \cdot \,\vec {\eta }_\rho )\,|\vec {x}_\rho |^{-1} \;{\mathrm{d}}\rho = - \sum _{i=1}^2 \sum _{p \in \partial _i I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\,\vec {\eta }(p)\cdot \,\vec {e}_{3-i} \quad \forall \ \vec {\eta } \in \underline{V}_\partial . \end{aligned}$$
(2.19b)

We note that (2.19b) weakly imposes (2.17) and (2.10a, b). We observe that (2.18a) degenerates for \(\vec {x}\cdot \,\vec {e}_1 = 0\), i.e. when \(\rho \in \partial _0 I\). Hence this degeneracy is balanced by the condition (2.16). In fact, on recalling (2.7) it holds that

$$\begin{aligned} \lim _{\rho \rightarrow \rho _0} \frac{\vec {\nu }(\rho ,t)\cdot \,\vec {e}_1}{\vec {x}(\rho ,t)\cdot \,\vec {e}_1} = \lim _{\rho \rightarrow \rho _0} \frac{\vec {\nu }_\rho (\rho ,t)\cdot \,\vec {e}_1}{\vec {x}_\rho (\rho ,t)\cdot \,\vec {e}_1}= & {} \vec {\nu }_s(\rho _0,t)\cdot \,\vec {\tau }(\rho _0,t) = -\varkappa (\rho _0,t) \nonumber \\&\qquad \qquad \forall \ \rho _0\in \partial _0 I,\ \forall \ t \in [0,T]. \nonumber \\ \end{aligned}$$
(2.20)

We remark that the weak formulation \(({\mathscr {A}})\) is close in spirit to the weak formulations introduced in [5, 7] for mean curvature flow. In particular, the tangential component of \(\vec {x}_t\) is not prescribed, which on the discrete level leads to an equidistribution property.

Choosing \(\vec {\eta } = (\vec {x}\cdot \,\vec {e}_1)\,\vec {x}_t \in \underline{V}_\partial \) in (2.19b) and \(\chi = (\vec {x}\cdot \,\vec {e}_1)\,(\vec {x}_t\cdot \,\vec {\nu })\) in (2.19a), we obtain on recalling (2.12), \(\vec {x}_t \in \underline{V}_\partial \), (2.4) and (2.2a) that

$$\begin{aligned}&\frac{\mathrm{d}}{{\mathrm{d}}t}\, E(\vec {x}(t)) \nonumber \\&\quad = 2\,\pi \,\int _I \left[ \vec {x}_t\cdot \,\vec {e}_1 + \vec {x}\cdot \,\vec {e}_1\,\frac{(\vec {x}_t)_\rho \cdot \,\vec {x}_\rho }{|\vec {x}_\rho |^2} \right] |\vec {x}_\rho | \;{\mathrm{d}}\rho + 2\,\pi \sum _{i=1}^2 \sum _{p\in \partial _i I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\,(\vec {x}(p,t)\cdot \,\vec {e}_1)\,\vec {x}_t(p,t)\cdot \,\vec {e}_{3-i} \nonumber \\&\quad = 2\,\pi \,\int _I \vec {x}_t.\left[ \vec {e}_1 - (\vec {e}_1\cdot \,\vec {\tau })\,\vec {\tau } \right] |\vec {x}_\rho | \;{\mathrm{d}}\rho - 2\,\pi \,\int _I (\vec {x}\cdot \,\vec {e}_1)\,\varkappa \,\vec {\nu }\cdot \, \vec {x}_t\, |\vec {x}_\rho | \;{\mathrm{d}}\rho \nonumber \\&\quad = 2\,\pi \,\int _I (\vec {x}_t\cdot \,\vec {\nu })\,\vec {e}_1\cdot \,\vec {\nu }\,|\vec {x}_\rho | \;{\mathrm{d}}\rho - 2\,\pi \,\int _I (\vec {x}\cdot \,\vec {e}_1)\,\varkappa \, \vec {x}_t\cdot \,\vec {\nu }\, |\vec {x}_\rho | \;{\mathrm{d}}\rho \nonumber \\&\quad = - 2\,\pi \,\int _I \vec {x}\cdot \,\vec {e}_1 \left[ \varkappa - \frac{\vec {\nu }\cdot \,\vec {e}_1}{\vec {x}\cdot \,\vec {e}_1}\right] \vec {x}_t\cdot \,\vec {\nu }\, |\vec {x}_\rho | \;{\mathrm{d}}\rho = - 2\,\pi \,\int _I \vec {x}\cdot \,\vec {e}_1\left( \vec {x}_t\cdot \,\vec {\nu }\right) ^2 |\vec {x}_\rho | \;{\mathrm{d}}\rho \le 0. \end{aligned}$$
(2.21)

An alternative strong formulation of mean curvature flow, in the axisymmetric setting, to (2.18a) is given by

$$\begin{aligned} \vec {x}_t = \vec {\varkappa } - \frac{\vec {\nu }\cdot \,\vec {e}_1}{\vec {x}\cdot \,\vec {e}_1}\,\vec {\nu }, \end{aligned}$$
(2.22)

with (2.18b), where we have recalled (2.7). We consider the following weak formulation of (2.22).

\(({\mathscr {B}})\) Let \(\vec {x}(0) \in \underline{V}_{\partial _0}\). For \(t \in (0,T]\) find \(\vec {x}(t) \in [H^1(I)]^2\), with \(\vec {x}_t(t) \in \underline{V}_\partial \), and \(\vec {\varkappa }(t)\in [L^2(I)]^2\) such that

$$\begin{aligned}&\int _I \vec {x}_t\cdot \,\vec {\chi }\,|\vec {x}_\rho |\;{\mathrm{d}}\rho = \int _I \left( \vec {\varkappa }\cdot \,\vec {\chi } - \frac{\vec {\nu }\cdot \,\vec {e}_1}{\vec {x}\cdot \,\vec {e}_1}\,\vec {\nu }\cdot \, \vec {\chi }\right) |\vec {x}_\rho | \;{\mathrm{d}}\rho \qquad \forall \ \vec {\chi } \in [L^2(I)]^2, \end{aligned}$$
(2.23a)
$$\begin{aligned}&\int _I \vec {\varkappa }\cdot \,\vec {\eta }\, |\vec {x}_\rho | \;{\mathrm{d}}\rho + \int _I (\vec {x}_\rho \cdot \,\vec {\eta }_\rho )\,|\vec {x}_\rho |^{-1} \;{\mathrm{d}}\rho = - \sum _{i=1}^2 \sum _{p \in \partial _i I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\,\vec {\eta }(p)\cdot \,\vec {e}_{3-i} \qquad \forall \ \vec {\eta } \in \underline{V}_\partial . \end{aligned}$$
(2.23b)

Similarly to (2.19b), we observe that (2.23b) weakly imposes (2.17) and (2.10a, b). We remark that the weak formulation \(({\mathscr {B}})\) in some sense is close in spirit to the weak formulations introduced in [20, 21] for mean curvature flow. In particular, the tangential component of \(\vec {x}_t\) is fixed to be zero, as the right hand side of (2.23a) is normal, recall (2.7).

Choosing \(\vec {\chi } = \vec {\eta } = (\vec {x}\cdot \,\vec {e}_1)\,\vec {x}_t \in \underline{V}_\partial \) in (2.23a, b), we obtain, similarly to (2.21), that

$$\begin{aligned}&\frac{\mathrm{d}}{{\mathrm{d}}t}\, E(\vec {x}(t)) \nonumber \\&\quad = 2\,\pi \,\int _I \left[ \vec {x}_t\cdot \,\vec {e}_1 + \vec {x}\cdot \,\vec {e}_1\,\frac{(\vec {x}_t)_\rho \cdot \,\vec {x}_\rho }{|\vec {x}_\rho |^2} \right] |\vec {x}_\rho | \;{\mathrm{d}}\rho + 2\,\pi \sum _{i=1}^2 \sum _{p\in \partial _i I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\,(\vec {x}(p,t)\cdot \,\vec {e}_1)\,\vec {x}_t(p,t)\cdot \,\vec {e}_{3-i} \nonumber \\&\quad = 2\,\pi \,\int _I (\vec {x}_t\cdot \,\vec {\nu })\,\vec {e}_1\cdot \,\vec {\nu }\,|\vec {x}_\rho | \;{\mathrm{d}}\rho - 2\,\pi \,\int _I (\vec {x}\cdot \,\vec {e}_1)\,\vec {\varkappa }\cdot \, \vec {x}_t\, |\vec {x}_\rho | \;{\mathrm{d}}\rho \nonumber \\&\quad = - 2\,\pi \,\int _I \vec {x}\cdot \,\vec {e}_1 \left[ \vec {\varkappa } - \frac{\vec {\nu }\cdot \,\vec {e}_1}{\vec {x}\cdot \,\vec {e}_1}\,\vec {\nu }\right] \cdot \, \vec {x}_t\, |\vec {x}_\rho | \;{\mathrm{d}}\rho = - 2\,\pi \,\int _I \vec {x}\cdot \,\vec {e}_1\,|\vec {x}_t|^2 |\vec {x}_\rho | \;{\mathrm{d}}\rho \le 0. \end{aligned}$$
(2.24)

We remark that it does not appear possible to mimic either (2.21) for \(({\mathscr {A}})\) or (2.24) for \(({\mathscr {B}})\) on the discrete level. Hence, in order to develop stable approximations, we investigate alternative formulations based on (2.15). The first formulation corresponds to the strong formulation \((\vec {x}\cdot \,\vec {e}_1)\,\vec {x}_t\cdot \,\vec {\nu } = \vec {x}\cdot \,\vec {e}_1\,\varkappa _{{\mathscr {S}}}\), together with (2.15).

\(({\mathscr {C}})\) Let \(\vec {x}(0) \in \underline{V}_{\partial _0}\). For \(t \in (0,T]\) find \(\vec {x}(t) \in [H^1(I)]^2\), with \(\vec {x}_t(t) \in \underline{V}_\partial \), and \(\varkappa _{{\mathscr {S}}}(t) \in L^2(I)\) such that

$$\begin{aligned}&\int _I (\vec {x} \cdot \,\vec {e}_1)\,\vec {x}_t\cdot \,\vec {\nu }\, \chi \,|\vec {x}_\rho |\;{\mathrm{d}}\rho = \int _I \vec {x} \cdot \,\vec {e}_1\,\varkappa _{{\mathscr {S}}}\,\chi \, |\vec {x}_\rho |\;{\mathrm{d}}\rho \qquad \forall \ \chi \in L^2(I), \end{aligned}$$
(2.25a)
$$\begin{aligned}&\int _I \vec {x} \cdot \,\vec {e}_1\,\varkappa _{{\mathscr {S}}}\,\vec {\nu }\cdot \, \vec {\eta }\,|\vec {x}_\rho |\;{\mathrm{d}}\rho + \int _I \left[ \vec {\eta } \cdot \,\vec {e}_1 + \vec {x}\cdot \,\vec {e}_1 \,\frac{\vec {x}_\rho \cdot \,\vec {\eta }_\rho }{|\vec {x}_\rho |^2} \right] |\vec {x}_\rho | \;{\mathrm{d}}\rho \nonumber \\&\quad = - \sum _{i=1}^2 \sum _{p \in \partial _i I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\, (\vec {x}(p,t)\cdot \,\vec {e}_1)\,\vec {\eta }(p)\cdot \,\vec {e}_{3-i} \qquad \forall \ \vec {\eta } \in \underline{V}_\partial . \end{aligned}$$
(2.25b)

The second formulation corresponds to the strong formulation \((\vec {x}\cdot \,\vec {e}_1)\,\vec {x}_t = \vec {x}\cdot \,\vec {e}_1\,\vec {\varkappa }_{{\mathscr {S}}}\), where \(\vec {\varkappa }_{{\mathscr {S}}} = \varkappa _{{\mathscr {S}}}\,\vec {\nu }\), together with (2.15).

\(({\mathscr {D}})\) Let \(\vec {x}(0) \in \underline{V}_{\partial _0}\). For \(t \in (0,T]\) find \(\vec {x}(t) \in [H^1(I)]^2\), with \(\vec {x}_t(t) \in \underline{V}_\partial \), and \(\vec {\varkappa }_{{\mathscr {S}}}(t) \in [L^2(I)]^2\) such that

$$\begin{aligned}&\int _I (\vec {x} \cdot \,\vec {e}_1)\,\vec {x}_t\cdot \,\vec {\chi }\,|\vec {x}_\rho |\;{\mathrm{d}}\rho = \int _I (\vec {x} \cdot \,\vec {e}_1)\,\vec {\varkappa }_{{\mathscr {S}}}\cdot \,\vec {\chi }\, |\vec {x}_\rho |\;{\mathrm{d}}\rho \qquad \forall \ \vec {\chi } \in [L^2(I)]^2, \end{aligned}$$
(2.26a)
$$\begin{aligned}&\int _I (\vec {x} \cdot \,\vec {e}_1)\,\vec {\varkappa }_{{\mathscr {S}}}\cdot \, \vec {\eta }\,|\vec {x}_\rho |\;{\mathrm{d}}\rho + \int _I \left[ \vec {\eta } \cdot \,\vec {e}_1 + \vec {x}\cdot \,\vec {e}_1 \,\frac{\vec {x}_\rho \cdot \,\vec {\eta }_\rho }{|\vec {x}_\rho |^2} \right] |\vec {x}_\rho | \;{\mathrm{d}}\rho \nonumber \\&\quad = - \sum _{i=1}^2 \sum _{p \in \partial _i I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\, (\vec {x}(p,t)\cdot \,\vec {e}_1)\,\vec {\eta }(p)\cdot \,\vec {e}_{3-i} \qquad \forall \ \vec {\eta } \in \underline{V}_\partial . \end{aligned}$$
(2.26b)

We note that the variational formulation for \(\vec {\varkappa }_{{\mathscr {S}}}\) in (2.26b) has previously been employed in [26, p. 124].

Choosing \(\vec {\eta } = \vec {x}_t \in \underline{V}_\partial \) in (2.25b) and \(\chi = \varkappa _{{\mathscr {S}}}\) in (2.25a), we obtain for the formulation \(({\mathscr {C}})\), on recalling (2.12), that

$$\begin{aligned} -\frac{1}{2\,\pi }\, \frac{\mathrm{d}}{{\mathrm{d}}t}\, E(\vec {x}(t)) = \int _I \vec {x} \cdot \,\vec {e}_1\, |\varkappa _{{\mathscr {S}}}|^2\,|\vec {x}_\rho |\;{\mathrm{d}}\rho = \int _I \vec {x} \cdot \,\vec {e}_1\, (\vec {x}_t\cdot \,\vec {\nu })^2\,|\vec {x}_\rho |\;{\mathrm{d}}\rho \ge 0. \end{aligned}$$
(2.27)

Similarly, choosing \(\vec {\eta } = \vec {x}_t \in \underline{V}_\partial \) in (2.26b) and \(\vec {\chi } = \vec {\varkappa }_{{\mathscr {S}}}\) in (2.26a), we obtain for the formulation \(({\mathscr {D}})\) that

$$\begin{aligned} -\frac{1}{2\,\pi }\, \frac{\mathrm{d}}{{\mathrm{d}}t}\, E(\vec {x}(t)) = \int _I \vec {x} \cdot \,\vec {e}_1\, |\vec {\varkappa }_{{\mathscr {S}}}|^2\,|\vec {x}_\rho |\;{\mathrm{d}}\rho = \int _I \vec {x} \cdot \,\vec {e}_1\,|\vec {x}_t|^2\,|\vec {x}_\rho |\;{\mathrm{d}}\rho \ge 0. \end{aligned}$$
(2.28)

We observe that (2.25b) and (2.26b) weakly impose (2.10a, b). But, in contrast to (2.19b) and (2.23b), it is not obvious that they also weakly impose (2.17), due to the presence of the degenerate weight \(\vec {x}\cdot \,\vec {e}_1\). However, we show in the “Appendix A” that in fact they also weakly impose (2.17).

We also note that the formulation \(({\mathscr {C}})\) is loosely related to \(({\mathscr {A}})\), in the sense that the tangential component of \(\vec {x}_t\) is not prescribed. But in contrast to \(({\mathscr {A}})\), discretizations of \(({\mathscr {C}})\) cannot be shown to have an equidistribution property. In a similar way, the formulation \(({\mathscr {D}})\) is loosely related to \(({\mathscr {B}})\), in the sense that the velocity \(\vec {x}_t\) is purely in the normal direction, recall (2.15) and (2.7). Finally, we observe that the variable \(\varkappa _{{\mathscr {S}}}\) can be eliminated from \(({\mathscr {C}})\), by choosing \(\chi = \vec {\nu }\cdot \,\vec {\eta }\) in (2.25a) for \(\vec {\eta } \in \underline{V}_\partial \), and then combining (2.25a) and (2.25b). Similarly, \(\vec {\varkappa }_{{\mathscr {S}}}\) can be eliminated from \(({\mathscr {D}})\) by choosing \(\vec {\chi } = \vec {\eta }\) in (2.26a) for \(\vec {\eta } \in \underline{V}_\partial \), and then combining (2.26a) and (2.26b). We remark that the formulation \(({\mathscr {C}})\), with the variable \(\varkappa _{{\mathscr {S}}}\), as well as the formulation \(({\mathscr {A}})\), are useful with a view towards introducing numerical approximations of higher order flows, such as surface diffusion, see [11].

2.2 Nonlinear mean curvature flow

It is a simple matter to extend the formulations \(({\mathscr {A}})\) and \(({\mathscr {C}})\) to the nonlinear flow (1.2). In principle this can also be achieved for \(({\mathscr {B}})\) and \(({\mathscr {D}})\), but as the mean curvature needs to be recovered from the mean curvature vector, the resulting formulations are less natural. Hence we concentrate on \(({\mathscr {A}})\) and \(({\mathscr {C}})\). For the former, replacing the right hand side in (2.19a) with \(\int _I f(\varkappa - \frac{\vec {\nu }\cdot \,\vec {e}_1}{\vec {x}\cdot \,\vec {e}_1} )\, \chi \,|\vec {x}_\rho | \;{\mathrm{d}}\rho \) yields a weak formulation for (1.2), which we call \(({\mathscr {A}}^f)\). Similarly, replacing \(\varkappa _{{\mathscr {S}}}\) with \(f(\varkappa _{{\mathscr {S}}})\) in (2.25a) generalizes \(({\mathscr {C}})\) to \(({\mathscr {C}}^f)\) for (1.2).

Similarly to (2.27), and using the same choices of \(\vec {\eta }\) and \(\chi \), it can be shown that solutions to \(({\mathscr {C}}^f)\) satisfy

$$\begin{aligned} -\frac{1}{2\,\pi }\, \frac{\mathrm{d}}{{\mathrm{d}}t}\, E(\vec {x}(t)) = \int _I \vec {x} \cdot \,\vec {e}_1\, f(\varkappa _{{\mathscr {S}}})\,\varkappa _{{\mathscr {S}}}\,|\vec {x}_\rho |\;{\mathrm{d}}\rho , \end{aligned}$$
(2.29)

which yields stability if f is monotonically increasing with \(f(0)=0\).

Finally, we may also generalize these nonlinear formulations to the volume preserving flow (1.3).

\(({\mathscr {A}}^{f,V})\) Let \(\vec {x}(0) \in \underline{V}_{\partial _0}\). For \(t \in (0,T]\) find \(\vec {x}(t) \in [H^1(I)]^2\), with \(\vec {x}_t(t) \in \underline{V}_\partial \), and \(\varkappa (t)\in L^2(I)\) such that

$$\begin{aligned}&\int _I \vec {x}_t\cdot \,\vec {\nu }\,\chi \,|\vec {x}_\rho |\;{\mathrm{d}}\rho = \int _I f\left( \varkappa - \frac{\vec {\nu }\cdot \,\vec {e}_1}{\vec {x}\cdot \,\vec {e}_1} \right) \chi \,|\vec {x}_\rho | \;{\mathrm{d}}\rho \nonumber \\&\qquad - \frac{\int _I \vec {x}\cdot \,\vec {e}_1\, f(\varkappa - \frac{\vec {\nu }\cdot \,\vec {e}_1}{\vec {x}\cdot \,\vec {e}_1})\, |\vec {x}_\rho | \;{\mathrm{d}}\rho }{\int _I \vec {x}\cdot \,\vec {e}_1\,|\vec {x}_\rho | \;{\mathrm{d}}\rho }\, \int _I \chi \,|\vec {x}_\rho | \;{\mathrm{d}}\rho \quad \forall \ \chi \in L^2(I), \end{aligned}$$
(2.30a)
$$\begin{aligned}&\int _I \varkappa \,\vec {\nu }\cdot \,\vec {\eta }\, |\vec {x}_\rho | \;{\mathrm{d}}\rho + \int _I (\vec {x}_\rho \cdot \,\vec {\eta }_\rho )\,|\vec {x}_\rho |^{-1} \;{\mathrm{d}}\rho = - \sum _{i=1}^2 \sum _{p \in \partial _i I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\,\vec {\eta }(p)\cdot \,\vec {e}_{3-i} \quad \forall \ \vec {\eta } \in \underline{V}_\partial . \end{aligned}$$
(2.30b)

Choosing \(\chi = 2\,\pi \,\vec {x}\cdot \,\vec {e}_1\) in (2.30a) yields, on recalling (1.5), that

$$\begin{aligned} \pm \frac{\mathrm{d}}{{\mathrm{d}}t}\,{\mathscr {L}}^3(\varOmega (t)) = \int _{{\mathscr {S}}(t)} {\mathscr {V}}_{{\mathscr {S}}} \;{\mathrm{d}}{{\mathscr {H}}}^{2} = 2\,\pi \,\int _I (\vec {x}\cdot \,\vec {e}_1) \,\vec {x}_t\cdot \,\vec {\nu }\,|\vec {x}_\rho |\;{\mathrm{d}}\rho = 0, \end{aligned}$$
(2.31)

where \({\mathscr {S}}(t) = \partial \varOmega (t)\), and where the sign in (2.31) depends on whether \(\vec {\mathrm{n}}_{{\mathscr {S}}}\) is the outer or inner normal to \(\varOmega (t)\) on \({\mathscr {S}}(t)\), recall (2.5) and (2.6).

\(({\mathscr {C}}^{f,V})\) Let \(\vec {x}(0) \in \underline{V}_{\partial _0}\). For \(t \in (0,T]\) find \(\vec {x}(t) \in [H^1(I)]^2\), with \(\vec {x}_t(t) \in \underline{V}_\partial \), and \(\varkappa _{{\mathscr {S}}}(t) \in L^2(I)\) such that

$$\begin{aligned}&\int _I (\vec {x} \cdot \,\vec {e}_1)\,\vec {x}_t\cdot \,\vec {\nu }\, \chi \,|\vec {x}_\rho |\;{\mathrm{d}}\rho = \int _I \vec {x} \cdot \,\vec {e}_1\,f(\varkappa _{{\mathscr {S}}})\,\chi \, |\vec {x}_\rho |\;{\mathrm{d}}\rho \nonumber \\&\qquad - \frac{\int _I \vec {x}\cdot \,\vec {e}_1\,f(\varkappa _{{\mathscr {S}}})\, |\vec {x}_\rho | \;{\mathrm{d}}\rho }{\int _I \vec {x}\cdot \,\vec {e}_1\,|\vec {x}_\rho | \;{\mathrm{d}}\rho }\, \int _I \vec {x} \cdot \,\vec {e}_1\,\chi \,|\vec {x}_\rho | \;{\mathrm{d}}\rho \qquad \forall \ \chi \in L^2(I), \end{aligned}$$
(2.32a)
$$\begin{aligned}&\int _I \vec {x} \cdot \,\vec {e}_1\,\varkappa _{{\mathscr {S}}}\,\vec {\nu }\cdot \, \vec {\eta }\,|\vec {x}_\rho |\;{\mathrm{d}}\rho + \int _I \left[ \vec {\eta } \cdot \,\vec {e}_1 + \vec {x}\cdot \,\vec {e}_1 \,\frac{\vec {x}_\rho \cdot \,\vec {\eta }_\rho }{|\vec {x}_\rho |^2} \right] |\vec {x}_\rho | \;{\mathrm{d}}\rho \nonumber \\&\quad = - \sum _{i=1}^2 \sum _{p \in \partial _i I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\, (\vec {x}(p,t)\cdot \,\vec {e}_1)\,\vec {\eta }(p)\cdot \,\vec {e}_{3-i} \qquad \forall \ \vec {\eta } \in \underline{V}_\partial . \end{aligned}$$
(2.32b)

Choosing \(\chi = 2\,\pi \) in (2.32a) yields (2.31), as before. Moreover, and similarly to (2.29), it can be shown for solutions of \(({\mathscr {C}}^{f,V})\) in the case (1.4c) that

$$\begin{aligned} -\frac{1}{2\,\pi }\, \frac{\mathrm{d}}{{\mathrm{d}}t}\, E(\vec {x}(t))&= \int _I \vec {x} \cdot \,\vec {e}_1\, |\varkappa _{{\mathscr {S}}}|^2\,|\vec {x}_\rho |\;{\mathrm{d}}\rho - \left[ \int _I \vec {x}\cdot \,\vec {e}_1\,|\vec {x}_\rho | \;{\mathrm{d}}\rho \right] ^{-1} \left| \int _I \vec {x}\cdot \,\vec {e}_1\,\varkappa _{{\mathscr {S}}}\, |\vec {x}_\rho | \;{\mathrm{d}}\rho \right| ^2 \nonumber \\&\ge 0, \end{aligned}$$
(2.33)

where we have noted the Cauchy–Schwarz inequality. It does not appear possible to extend the stability result (2.33) to the case of more general f.

3 Semidiscrete schemes

Let \([0,1]=\bigcup _{j=1}^J I_j\), \(J\ge 3\), be a decomposition of [0, 1] into intervals given by the nodes \(q_j\), \(I_j=[q_{j-1},q_j]\). For simplicity, and without loss of generality, we assume that the subintervals form an equipartitioning of [0, 1], i.e. that

$$\begin{aligned} q_j = j\,h,\quad \text{ with }\quad h = J^{-1},\quad j=0,\ldots , J. \end{aligned}$$
(3.1)

Clearly, if \(I={{\mathbb {R}}} \slash {{\mathbb {Z}}}\) we identify \(0=q_0 = q_J=1\).

The necessary finite element spaces are given by

$$\begin{aligned} V^h = \left\{ \chi \in C({\overline{I}}){:}\,\chi \!\mid _{I_j} \text{ is } \text{ linear }\ \forall \ j=1\rightarrow J\right\} , \end{aligned}$$

\(\underline{V}^h= [V^h]^2\), \(\underline{V}^h_{\partial _0}= \underline{V}^h\cap \underline{V}_{\partial _0}\) and \(\underline{V}^h_\partial = \underline{V}^h\cap \underline{V}_\partial \). Let \(\{\chi _j\}_{j=j_0}^J\) denote the standard basis of \(V^h\), where \(j_0 = 0\) if \(I = (0,1)\) and \(j_0 = 1\) if \(I={{\mathbb {R}}} \slash {{\mathbb {Z}}}\). For later use, we let \(\pi ^h{:}\,C({\overline{I}})\rightarrow V^h\) be the standard interpolation operator at the nodes \(\{q_j\}_{j=0}^J\). We extend \(\pi ^h\) to \(\vec \pi ^h\) for vector valued functions. Let \((\cdot ,\cdot )\) denote the \(L^2\)-inner product on I, and define the mass lumped \(L^2\)-inner product \((f,g)^h\), for two piecewise continuous functions, with possible jumps at the nodes \(\{q_j\}_{j=1}^J\), via

$$\begin{aligned} ( f, g )^h = \tfrac{1}{2}\,h\,\sum _{j=1}^J \left[ (f\,g)(q_j^-) + (f\,g)(q_{j-1}^+)\right] , \end{aligned}$$
(3.2)

where we define \(f(q_j^\pm )=\lim \nolimits _{\delta \searrow 0}f(q_j\pm \delta )\). The definition (3.2) naturally extends to vector valued functions. It is easily shown that

$$\begin{aligned} (\eta , \eta ) \le (\eta ,\eta )^h \le 3\,(\eta ,\eta )\qquad \forall \ \eta \in V^h. \end{aligned}$$
(3.3)

Let \((\vec {X}^h(t))_{t\in [0,T]}\), with \(\vec {X}^h(t)\in \underline{V}^h_{\partial _0}\), be an approximation to \((\vec {x}(t))_{t\in [0,T]}\) and define \(\varGamma ^h(t) = \vec {X}^h(t)({\overline{I}})\). Throughout this section we assume that

$$\begin{aligned} \vec {X}^h(\rho ,t) \cdot \,\vec {e}_1 > 0 \quad \forall \ \rho \in {\overline{I}}{\setminus } \partial _0 I, \quad \forall \ t \in [0,T]. \end{aligned}$$
(3.4)

Assuming that \(|\vec {X}^h_\rho | > 0\) almost everywhere on I, and similarly to (2.4), we set

$$\begin{aligned} \vec {\tau }^h = \vec {X}^h_s = \frac{\vec {X}^h_\rho }{|\vec {X}^h_\rho |} \qquad \text{ and } \qquad \vec {\nu }^h = -(\vec {\tau }^h)^\perp . \end{aligned}$$

We note that

$$\begin{aligned} (\vec {\tau }^h, (\vec {\pi }^h[\chi \,\vec {\eta }])_\rho ) = (\vec {\tau }^h, (\chi \,\vec {\eta })_\rho ) \qquad \forall \ \chi \in C({\overline{I}}),\ \vec {\eta } \in [C({\overline{I}})]^2. \end{aligned}$$
(3.5)

For later use, we let \(\vec {\omega }^h \in {\underline{V}}^h\) be the mass-lumped \(L^2\)-projection of \(\vec {\nu }^h\) onto \({\underline{V}}^h\), i.e.

$$\begin{aligned} \left( \vec {\omega }^h, \vec {\varphi } \, |\vec {X}^h_\rho | \right) ^h = \left( \vec {\nu }^h, \vec {\varphi } \, |\vec {X}^h_\rho | \right) = \left( \vec {\nu }^h, \vec {\varphi } \, |\vec {X}^h_\rho | \right) ^h \qquad \forall \ \vec {\varphi }\in {\underline{V}}^h. \end{aligned}$$
(3.6)

Recall from (2.8) and (2.9) that

$$\begin{aligned} E(\vec {X}^h(t))&= 2\,\pi \left( \vec {X}^h(t)\cdot \,\vec {e}_1, |\vec {X}^h_\rho (t)|\right) \nonumber \\&\quad + 2\,\pi \sum _{p\in \partial _1 I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\,(\vec {X}^h(p,t)\cdot \,\vec {e}_1)\,\vec {X}^h(p,t)\cdot \,\vec {e}_2 + \pi \sum _{p\in \partial _2 I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\,(\vec {X}^h(p,t)\cdot \,\vec {e}_1)^2. \end{aligned}$$
(3.7)

We have, similarly to (2.12), that

$$\begin{aligned} \frac{\mathrm{d}}{{\mathrm{d}}t}\, E(\vec {X}^h(t))&= 2\,\pi \left( \left[ \vec {X}^h_t\cdot \,\vec {e}_1 + \vec {X}^h\cdot \,\vec {e}_1 \,\frac{(\vec {X}^h_t)_\rho \cdot \,\vec {X}^h_\rho }{|\vec {X}^h_\rho |^2} \right] , |\vec {X}^h_\rho | \right) \nonumber \\&\quad +\,2\,\pi \,\sum _{p\in \partial _1 I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\left[ (\vec {X}^h_t(p,t)\cdot \,\vec {e}_1)\, \vec {X}^h(p,t)\cdot \,\vec {e}_2 + (\vec {X}^h(p,t)\cdot \,\vec {e}_1)\, \vec {X}^h_t(p,t)\cdot \,\vec {e}_2 \right] \nonumber \\&\quad +\, 2\,\pi \,\sum _{p\in \partial _2 I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\,(\vec {X}^h(p,t)\cdot \,\vec {e}_1)\, \vec {X}^h_t(p,t)\cdot \,\vec {e}_1. \end{aligned}$$
(3.8)

3.1 Mean curvature flow

In view of the degeneracy on the right hand side of (2.18a), and on recalling (2.20) and (3.6), we introduce, given a \(\kappa ^h(t) \in V^h\), the function \({\mathfrak {K}}^h(\kappa ^h(t),t) \in V^h\) such that

$$\begin{aligned}{}[{\mathfrak {K}}^h (\kappa ^h(t),t)](q_j) = {\left\{ \begin{array}{ll} \dfrac{\vec {\omega }^h(q_j,t)\cdot \,\vec {e}_1}{\vec {X}^h(q_j,t)\cdot \,\vec {e}_1} &{}\quad q_j \in {\overline{I}} {\setminus } \partial _0 I, \\ - \kappa ^h(q_j,t) &{}\quad q_j \in \partial _0 I. \end{array}\right. } \end{aligned}$$
(3.9)

Our semidiscrete finite element approximation of \(({\mathscr {A}})\), (2.19a, b), is given as follows.

\(({\mathscr {A}}_h)^h\) Let \(\vec {X}^h(0) \in \underline{V}^h_{\partial _0}\). For \(t \in (0,T]\) find \(\vec {X}^h(t) \in \underline{V}^h\), with \(\vec {X}^h_t(t) \in \underline{V}^h_\partial \), and \(\kappa ^h(t) \in V^h\) such that

$$\begin{aligned}&\left( \vec {X}^h_t, \chi \,\vec {\nu }^h\,|\vec {X}^h_\rho |\right) ^h = \left( \kappa ^h - {\mathfrak {K}}^h (\kappa ^h), \chi \,|\vec {X}^h_\rho | \right) ^h \qquad \forall \ \chi \in V^h, \end{aligned}$$
(3.10a)
$$\begin{aligned}&\left( \kappa ^h\,\vec {\nu }^h, \vec {\eta }\,|\vec {X}^h_\rho |\right) ^h + \left( \vec {X}^h_\rho , \vec {\eta }_\rho \,|\vec {X}^h_\rho |^{-1}\right) = - \sum _{i=1}^2 \sum _{p \in \partial _i I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\,\vec {\eta }(p)\cdot \,\vec {e}_{3-i} \qquad \forall \ \vec {\eta } \in \underline{V}^h_\partial . \end{aligned}$$
(3.10b)

Remark 3.1

Let \(\vec {h}_j(t) = \vec {X}^h(q_j,t) - \vec {X}^h(q_{j-1},t)\) for \(j = 1,\ldots , J\), and set \(\vec {h}_0 = \vec {h}_J\) if \(\partial I = \emptyset \). Then, if \((\vec {X}^h(t), \kappa ^h(t)) \in \underline{V}^h\times V^h\) satisfies (3.10b), it holds that

$$\begin{aligned} |\vec {h}_j(t)| = |\vec {h}_{j - 1}(t)| \quad \text{ if } \quad \vec {h}_j(t) \nparallel \vec {h}_{j - 1}(t) \quad {\left\{ \begin{array}{ll} j = 1,\ldots , J &{}\quad \partial I = \emptyset , \\ j = 2,\ldots , J &{}\quad \partial I \not = \emptyset . \end{array}\right. } \end{aligned}$$
(3.11)

The equidistribution property (3.11) can be shown by choosing \(\vec {\eta }=\chi _{j-1}\,[\vec {\omega }^h(q_{j-1},t)]^\perp \in \underline{V}^h_\partial \) in (3.10b), recall (3.6). See also [6, Remark 2.4] and [5, Remark 2.5] for more details.

We also remark that it follows from (3.6) that

$$\begin{aligned} |\vec {\omega }^h(q_{j-1},t)| < 1 \quad \text{ if } \quad \vec {h}_j(t) \nparallel \vec {h}_{j - 1}(t) \quad {\left\{ \begin{array}{ll} j = 1,\ldots , J &{}\quad \partial I = \emptyset , \\ j = 2,\ldots , J &{}\quad \partial I \not = \emptyset . \end{array}\right. } \end{aligned}$$
(3.12)

We note that mass lumping in (3.10b) is crucial for the proof of the equidistribution property (3.11). Hence we only consider the variant \(({\mathscr {A}}_h)^h\) with mass lumping. Of course, in the case \(\partial _0 I= \emptyset \), an alternative scheme to (3.10a, b) is

$$\begin{aligned} \left( \vec {X}^h_t, \chi \,\vec {\nu }^h\,|\vec {X}^h_\rho |\right) ^h = \left( \kappa ^h - \frac{\vec {\nu }^h\cdot \,\vec {e}_1}{\vec {X}^h\cdot \,\vec {e}_1}, \chi \,|\vec {X}^h_\rho | \right) ^h \qquad \forall \ \chi \in V^h, \end{aligned}$$
(3.13)

together with (3.10b). Note that if \(\partial _0 I= \emptyset \) then (3.10a) collapses to (3.13) with \(\vec {\nu }^h\) replaced by \(\vec {\omega }^h\). Unfortunately, neither choice appears to lead to a stability proof.

In an attempt to prove stability, we choose \(\vec {\eta } = \vec {\pi }^h[(\vec {X}^h\cdot \,\vec {e}_1)\,\vec {X}^h_t]\) in (3.10b). Then it follows from (3.8), \(\vec {X}^h_t \in \underline{V}^h_\partial \), (3.5) and (3.6) that

$$\begin{aligned} \frac{\mathrm{d}}{{\mathrm{d}}t}\,E(\vec {X}^h(t))&= 2\,\pi \left( \left[ \vec {X}^h_t\cdot \,\vec {e}_1 + \vec {X}^h\cdot \,\vec {e}_1\,\frac{(\vec {X}^h_t)_\rho \cdot \,\vec {X}^h_\rho }{|\vec {X}^h_\rho |^2} \right] , |\vec {X}^h_\rho | \right) \nonumber \\&\qquad + 2\,\pi \sum _{i=1}^2 \sum _{p\in \partial _i I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\,(\vec {X}^h(p,t)\cdot \,\vec {e}_1)\, \vec {X}^h_t(p,t)\cdot \,\vec {e}_{3-i}. \nonumber \\&= 2\,\pi \left( \vec {X}^h_t,\left[ \vec {e}_1 - (\vec {e}_1\cdot \,\vec {\tau }^h)\,\vec {\tau }^h \right] |\vec {X}^h_\rho |\right) - 2\,\pi \left( \vec {X}^h\cdot \,\vec {e}_1\,\kappa ^h\,\vec {\nu }^h, \vec {X}^h_t\, |\vec {X}^h_\rho |\right) ^h \nonumber \\&= 2\,\pi \left( \vec {X}^h_t\cdot \,\vec {\nu }^h,\vec {e}_1\cdot \,\vec {\nu }^h\, |\vec {X}^h_\rho |\right) - 2\,\pi \left( \vec {X}^h\cdot \,\vec {e}_1\,\kappa ^h, \vec {X}^h_t\cdot \,\vec {\nu }^h\, |\vec {X}^h_\rho |\right) ^h \nonumber \\&= 2\,\pi \left( \vec {X}^h_t\cdot \,\vec {\nu }^h,\vec {e}_1\cdot \,\vec {\nu }^h\,|\vec {X}^h_\rho |\right) ^h - 2\,\pi \left( \vec {X}^h\cdot \,\vec {e}_1\,\kappa ^h, \vec {X}^h_t\cdot \,\vec {\nu }^h\, |\vec {X}^h_\rho |\right) ^h \nonumber \\&= - 2\,\pi \left( \vec {X}^h\cdot \,\vec {e}_1\left[ \kappa ^h - \frac{\vec {\nu }^h\cdot \,\vec {e}_1}{\vec {X}^h\cdot \,\vec {e}_1}\right] , \vec {X}^h_t\cdot \,\vec {\nu }^h\, |\vec {X}^h_\rho |\right) ^h \nonumber \\&= - 2\,\pi \left( \vec {X}^h\cdot \,\vec {e}_1, \kappa ^h\, \vec {X}^h_t\cdot \,\vec {\omega }^h\, |\vec {X}^h_\rho |\right) ^h + 2\,\pi \left( \vec {\nu }^h\cdot \,\vec {e}_1, \vec {X}^h_t\cdot \,\vec {\nu }^h\, |\vec {X}^h_\rho |\right) ^h. \end{aligned}$$
(3.14)

Moreover, considering for simplicity the case \(\partial _0 I= \emptyset \), and choosing \(\chi = -\pi ^h[2\,\pi \,(\vec {X}^h\cdot \,\vec {e}_1)\,(\vec {X}^h_t\cdot \,\vec {\omega }^h)]\) in (3.13) yields, on noting (3.6), that

$$\begin{aligned} 0&\ge - 2\,\pi \left( \vec {X}^h\cdot \,\vec {e}_1,(\vec {X}^h_t\cdot \,\vec {\omega }^h)^2\, |\vec {X}^h_\rho | \right) ^h =-2\,\pi \left( \vec {X}^h\cdot \,\vec {e}_1\left[ \kappa ^h - \frac{\vec {\omega }^h\cdot \,\vec {e}_1}{\vec {X}^h\cdot \,\vec {e}_1}\right] , \vec {X}^h_t\cdot \,\vec {\omega }^h\, |\vec {X}^h_\rho |\right) ^h \nonumber \\&= -2\,\pi \left( \vec {X}^h\cdot \,\vec {e}_1, \kappa ^h\, \vec {X}^h_t\cdot \,\vec {\omega }^h\, |\vec {X}^h_\rho |\right) ^h + 2\,\pi \left( \vec {\omega }^h\cdot \,\vec {e}_1, \vec {X}^h_t\cdot \,\vec {\omega }^h\, |\vec {X}^h_\rho |\right) ^h \nonumber \\&= -2\,\pi \left( \vec {X}^h\cdot \,\vec {e}_1, \kappa ^h\, \vec {X}^h_t\cdot \,\vec {\omega }^h\, |\vec {X}^h_\rho |\right) ^h + 2\,\pi \left( \vec {\nu }^h\cdot \,\vec {e}_1, \vec {X}^h_t\cdot \,\vec {\omega }^h\, |\vec {X}^h_\rho |\right) ^h. \end{aligned}$$
(3.15)

Unfortunately, the right hand sides in (3.14) and (3.15) are not equal, recall (3.6), and so combining (3.14) and (3.15) does not yield a stability result. On the other hand, the function \((\vec {X}^h\cdot \,\vec {e}_1)\,(\vec {X}^h_t\cdot \,\vec {\nu }^h)\) is discontinuous, and so \(\pi ^h[(\vec {X}^h\cdot \,\vec {e}_1)\,(\vec {X}^h_t\cdot \,\vec {\nu }^h)]\) is not well-defined, and cannot be chosen as a test function in (3.10a) or (3.13).

However, the fully discrete variant of \(({\mathscr {A}}_h)^h\), (3.10a, b), in our numerical simulations did not exhibit any stability issues.

A semidiscrete approximation of \(({\mathscr {B}})\), (2.23a, b), is given as follows.

\(({\mathscr {B}}_h)^h\) Let \(\vec {X}^h(0) \in \underline{V}^h_{\partial _0}\). For \(t \in (0,T]\) find \(\vec {X}^h(t) \in \underline{V}^h\), with \(\vec {X}^h_t(t) \in \underline{V}^h_\partial \), and \(\vec {\kappa }^h(t) \in \underline{V}^h\), such that

$$\begin{aligned}&\left( \vec {X}^h_t, \vec {\chi }\,|\vec {X}^h_\rho | \right) ^h = \left( \vec {\kappa }^h - \vec {{\mathfrak {K}}}^h(\vec {\kappa }^h), \vec {\chi }\,|\vec {X}^h_\rho |\right) ^h \qquad \forall \ \vec {\chi } \in \underline{V}^h, \end{aligned}$$
(3.16a)
$$\begin{aligned}&\left( \vec {\kappa }^h, \vec {\eta }\,|\vec {X}^h_\rho |\right) ^h + \left( \vec {X}^h_\rho , \vec {\eta }_\rho \,|\vec {X}^h_\rho |^{-1}\right) = - \sum _{i=1}^2 \sum _{p \in \partial _i I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\,\vec {\eta }(p)\cdot \,\vec {e}_{3-i} \qquad \forall \ \vec {\eta } \in \underline{V}^h_\partial , \end{aligned}$$
(3.16b)

where \(\vec {{\mathfrak {K}}}^{h}(\vec {\kappa }^h) \in \underline{V}^h\) is such that

$$\begin{aligned}{}[\vec {{\mathfrak {K}}}^{h}(\vec {\kappa }^h(t),t))](q_j) = {\left\{ \begin{array}{ll} \dfrac{\vec {\omega }^h(q_j,t)\cdot \,\vec {e}_1}{\vec {X}^h(q_j,t)\cdot \,\vec {e}_1}\, \dfrac{\vec {\omega }^h(q_j,t)}{|\vec {\omega }^h(q_j,t)|^2} &{}\quad q_j \in {\overline{I}} {\setminus } \partial _0 I, \\ - \vec {\kappa }^h(q_j,t) &{}\quad q_j \in \partial _0 I. \end{array}\right. } \end{aligned}$$
(3.17)

The rescaling factor \(|\vec {\omega }^h(q_j,t)|^2\) in (3.17) normalizes the discrete vertex normals \(\vec {\omega }^h(q_j,t)\), recall (3.12), which is the most natural approach. Similarly to \(({\mathscr {A}}_h)^h\), it does not appear possible to prove a stability result for \(({\mathscr {B}}_h)^h\).

However, it turns out that approximations of the formulations \(({\mathscr {C}})\) and \(({\mathscr {D}})\) can be shown to be stable. In particular, our semidiscrete approximations of \(({\mathscr {C}})\), (2.25a, b), and \(({\mathscr {D}})\), (2.26a, b), are given as follows, where we first define

$$\begin{aligned} W^h = V^h,\quad W^h_{\partial _0} = \left\{ \chi \in V^h{:}\,\chi (\rho ) = 0 \quad \forall \ \rho \in \partial _0 I\right\} , \quad {\underline{W}}^h = {\underline{V}}^h,\quad {\underline{W}}^h_{\partial _0} = [W^h_{\partial _0}]^2. \end{aligned}$$

\(({\mathscr {C}}_h)^{(h)}\) Let \(\vec {X}^h(0) \in \underline{V}^h_{\partial _0}\). For \(t \in (0,T]\) find \(\vec {X}^h(t) \in \underline{V}^h\), with \(\vec {X}^h_t(t) \in \underline{V}^h_\partial \), and \(\kappa _{{\mathscr {S}}}^h(t) \in W^h_{(\partial _0)}\) such that

$$\begin{aligned}&\left( (\vec {X}^h\cdot \,\vec {e}_1)\, \vec {X}^h_t, \chi \,\vec {\nu }^h\,|\vec {X}^h_\rho |\right) ^{(h)} = \left( \vec {X}^h\cdot \,\vec {e}_1\,\kappa _{{\mathscr {S}}}^h, \chi \,|\vec {X}^h_\rho |\right) ^{(h)} \qquad \forall \ \chi \in W^h_{(\partial _0)}, \end{aligned}$$
(3.18a)
$$\begin{aligned}&\left( \vec {X}^h\cdot \,\vec {e}_1\,\kappa _{{\mathscr {S}}}^h\,\vec {\nu }^h, \vec {\eta }\,|\vec {X}^h_\rho |\right) ^{(h)} + \left( \vec {\eta } \cdot \,\vec {e}_1, |\vec {X}^h_\rho |\right) + \left( (\vec {X}^h\cdot \,\vec {e}_1)\, \vec {X}^h_\rho ,\vec {\eta }_\rho \, |\vec {X}^h_\rho |^{-1} \right) \nonumber \\&\quad = - \sum _{i=1}^2 \sum _{p \in \partial _i I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\, (\vec {X}^h(p,t)\cdot \,\vec {e}_1)\,\vec {\eta }(p)\cdot \,\vec {e}_{3-i} \qquad \forall \ \vec {\eta } \in \underline{V}^h_\partial . \end{aligned}$$
(3.18b)

Here and throughout we use the notation \(\cdot ^{(h)}\) to denote an expression with or without the superscript h, and similarly for the subscripts \(\cdot _{(\partial _0)}\). I.e. the scheme \(({\mathscr {C}}_h)^h\) employs mass lumping, recall (3.2), and seeks \(\kappa _S(t) \in W^h_{\partial _0}\), while the scheme \(({\mathscr {C}}_h)\) employs true integration throughout and seeks \(\kappa _S(t) \in W^h = V^h\).

\(({\mathscr {D}}_h)^{(h)}\) Let \(\vec {X}^h(0) \in \underline{V}^h_{\partial _0}\). For \(t \in (0,T]\) find \(\vec {X}^h(t) \in \underline{V}^h\), with \(\vec {X}^h_t(t) \in \underline{V}^h_\partial \), and \(\vec {\kappa }_{{\mathscr {S}}}^h(t) \in \underline{W}_{(\partial _0)}^h\) such that

$$\begin{aligned}&\left( (\vec {X}^h\cdot \,\vec {e}_1)\, \vec {X}^h_t,\vec {\chi }\,|\vec {X}^h_\rho |\right) ^{(h)} = \left( (\vec {X}^h\cdot \,\vec {e}_1)\,\vec {\kappa }_{{\mathscr {S}}}^h, \vec {\chi }\,|\vec {X}^h_\rho |\right) ^{(h)} \qquad \forall \ \vec {\chi } \in \underline{W}_{(\partial _0)}^h, \end{aligned}$$
(3.19a)
$$\begin{aligned}&\left( (\vec {X}^h\cdot \,\vec {e}_1)\,\vec {\kappa }_{{\mathscr {S}}}^h, \vec {\eta }\,|\vec {X}^h_\rho |\right) ^{(h)} + \left( \vec {\eta } \cdot \,\vec {e}_1, |\vec {X}^h_\rho |\right) + \left( (\vec {X}^h\cdot \,\vec {e}_1)\, \vec {X}^h_\rho ,\vec {\eta }_\rho \, |\vec {X}^h_\rho |^{-1} \right) \nonumber \\&\quad = - \sum _{i=1}^2 \sum _{p \in \partial _i I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\, (\vec {X}^h(p,t)\cdot \,\vec {e}_1)\,\vec {\eta }(p)\cdot \,\vec {e}_{3-i} \quad \forall \ \vec {\eta } \in \underline{V}^h_\partial . \end{aligned}$$
(3.19b)

We observe that \(({\mathscr {C}}_h)^{h}\) and \(({\mathscr {D}}_h)^{h}\) do not depend on the values of \(\kappa _{{\mathscr {S}}}^h\) and \(\vec {\kappa }_{{\mathscr {S}}}^h\), respectively, on \(\partial _0 I\). Hence we fix these values to be zero by requiring that \(\kappa _{{\mathscr {S}}}^h \in W^h_{\partial _0}\) and \(\vec {\kappa }_{{\mathscr {S}}}^h \in {\underline{W}}^h_{\partial _0}\), and by using a reduced set of test functions in (3.18a) and (3.19a). As a consequence, it seems at first that \(\vec {X}^h_t\) is not defined on \(\partial _0 I\). However, \(\vec {X}^h\) on \(\partial _0 I\) is determined through (3.18b) and (3.19b), respectively.

We have on choosing \(\chi = \kappa _{{\mathscr {S}}}^h\) in (3.18a), \(\vec {\chi } = \vec {\kappa }_{{\mathscr {S}}}^h\) in (3.19a) and \(\vec {\eta } = \vec {X}^h_t\) in (3.18b) and (3.19b), on recalling (3.8), that

$$\begin{aligned} -\frac{1}{2\,\pi }\, \frac{\mathrm{d}}{{\mathrm{d}}t}\, E(\vec {X}^h(t)) = {\left\{ \begin{array}{ll} \left( \vec {X}^h\cdot \,\vec {e}_1\, |\kappa _{{\mathscr {S}}}^h|^2, |\vec {X}^h_\rho |\right) ^{(h)},\\ \left( \vec {X}^h\cdot \,\vec {e}_1\, |\vec {\kappa }_{{\mathscr {S}}}^h|^2, |\vec {X}^h_\rho |\right) ^{(h)},\\ \end{array}\right. } \end{aligned}$$
(3.20)

respectively. This shows that both methods are stable, where we recall (3.4). Similarly to (2.28) and (2.27), we observe that (3.20) implies for \(({\mathscr {D}}_h)^{(h)}\) and \(({\mathscr {C}}_h)^h\) that

$$\begin{aligned} -\frac{1}{2\,\pi }\, \frac{\mathrm{d}}{{\mathrm{d}}t}\, E(\vec {X}^h(t)) = {\left\{ \begin{array}{ll} \left( \vec {X}^h\cdot \,\vec {e}_1\,|\vec {X}^h_t|^2, |\vec {X}^h_\rho |\right) ^{(h)},\\ \left( \vec {X}^h\cdot \,\vec {e}_1\, (\vec {X}^h_t\cdot \,\vec {\omega }^h)^2, |\vec {X}^h_\rho |\right) ^h, \end{array}\right. } \end{aligned}$$

respectively, where we have recalled (3.6). This shows that they can be interpreted as natural \(L^2\)-gradient flows of (3.7).

Table 2 Properties of the different semidiscrete schemes for mean curvature flow

We observe that it is possible to eliminate \(\vec {\kappa }^h_{{\mathscr {S}}}\) from the schemes \(({\mathscr {D}}_h)^{(h)}\), which yields (3.19b) with \(\vec {\kappa }^h_{{\mathscr {S}}}\) replaced by \(\vec {X}^h_t\). Similarly, \(\kappa ^h_{{\mathscr {S}}}\) can be removed from the scheme \(({\mathscr {C}}_h)^{h}\) to yield (3.18b) with \(\kappa ^h_{{\mathscr {S}}}\,\vec {\nu }^h\) replaced by \((\vec {X}^h_t\cdot \,\vec {\omega }^h)\,\vec {\omega }^h\), on recalling (3.6). For the scheme \(({\mathscr {C}}_h)\) this elimination procedure is not possible.

For the reader’s convenience, Table 2 summarises the main properties of all the schemes introduced in this section.

3.2 Nonlinear mean curvature flow

Replacing \(\kappa ^h - {\mathfrak {K}}^h(\kappa ^h)\) with \(f(\kappa ^h - {\mathfrak {K}}^h(\kappa ^h))\) in (3.10a) yields the scheme \(({\mathscr {A}}_h^f)^h\). Similarly, the scheme \(({\mathscr {A}}_h^{f,V})^h\) is given by (3.10a, b) with the right hand side in (3.10a) replaced by

$$\begin{aligned} \left( f(\kappa ^h - {\mathfrak {K}}^h(\kappa ^h)), \chi \,|\vec {X}^h_\rho |\right) ^h - \frac{\left( \vec {X}^h\cdot \,\vec {e}_1, f(\kappa ^h - {\mathfrak {K}}^h(\kappa ^h))\, |\vec {X}^h_\rho |\right) ^{h}}{\left( \vec {X}^h\cdot \,\vec {e}_1, |\vec {X}^h_\rho |\right) } \left( \chi ,|\vec {X}^h_\rho |\right) ^{h}. \end{aligned}$$
(3.21)

These two schemes inherit the equidistribution property, recall (3.11). Replacing \(\kappa ^h_{{\mathscr {S}}}\) with \(\pi ^h[f(\kappa ^h_{{\mathscr {S}}})]\) in (3.18a) yields the schemes \(({\mathscr {C}}_h^f)^{(h)}\) and similarly we can define \(({\mathscr {C}}_h^{f,V})^{(h)}\) by replacing the right hand side in (3.18a) by

$$\begin{aligned} \left( \vec {X}^h\cdot \,\vec {e}_1\,\pi ^h[f(\kappa ^h_{{\mathscr {S}}})], \chi \,|\vec {X}^h_\rho |\right) ^{(h)} - \frac{\left( \vec {X}^h\cdot \,\vec {e}_1, \pi ^h[f(\kappa ^h_{{\mathscr {S}}})]\, |\vec {X}^h_\rho |\right) ^{(h)}}{\left( \vec {X}^h\cdot \,\vec {e}_1, |\vec {X}^h_\rho |\right) } \left( \vec {X}^h\cdot \,\vec {e}_1,\chi \,|\vec {X}^h_\rho |\right) ^{(h)}. \end{aligned}$$
(3.22)

Similarly to (3.20), and using the same choices of \(\vec {\eta }\) and \(\chi \), it can be shown that solutions to the scheme \(({\mathscr {C}}_h^f)^{(h)}\) satisfy \(-\frac{1}{2\,\pi }\,\frac{\mathrm{d}}{{\mathrm{d}}t}\, E(\vec {X}^h(t)) = \left( (\vec {X}^h\cdot \,\vec {e}_1)\,f(\kappa _{{\mathscr {S}}}^h), \kappa _{{\mathscr {S}}}^h\, |\vec {X}^h_\rho |\right) ^{(h)}\), which yields a stability bound for \(({\mathscr {C}}_h^f)^{h}\) if f is monotonically increasing with \(f(0) = 0\). Of course, (3.22) is a discrete analogue of (2.29). Moreover, solutions to \(({\mathscr {C}}_h^{f,V})^{(h)}\), in the case (1.4c), satisfy

$$\begin{aligned}&-\frac{1}{2\,\pi }\,\frac{\mathrm{d}}{{\mathrm{d}}t}\, E(\vec {X}^h(t)) \nonumber \\&\qquad = \left( \vec {X}^h \cdot \,\vec {e}_1\, |\kappa _{{\mathscr {S}}}^h|^2,|\vec {X}^h_\rho | \right) ^{(h)} - \left[ \left( \vec {X}^h\cdot \,\vec {e}_1,|\vec {X}^h_\rho | \right) \right] ^{-1} \left| \left( \vec {X}^h\cdot \,\vec {e}_1,\kappa _{{\mathscr {S}}}^h\, |\vec {X}^h_\rho | \right) ^{(h)}\right| ^2 \ge 0, \end{aligned}$$

similarly to (2.33), where here we have also used a Cauchy–Schwarz inequality for the mass lumped inner product (3.2). Finally, solutions to the scheme \(({\mathscr {C}}_h^{f,V})\) conserve the volume of the domain \(\varOmega ^h(t) \subset {{\mathbb {R}}}^3\) that is enclosed by the three-dimensional axisymmetric surface \({\mathscr {S}}^h(t)\) that is generated by the curve \(\varGamma ^h(t)\). To see this, choose \(\chi = 2\,\pi \) in (3.18a), with the modified right hand side (3.22), to obtain

$$\begin{aligned} 0 = 2\,\pi \left( \vec {X}^h\cdot \,\vec {e}_1, \vec {X}^h_t\cdot \, \vec {\nu }^h\,|\vec {X}^h_\rho |\right) = \int _{{\mathscr {S}}^h(t)} {\mathscr {V}}^h_{{\mathscr {S}}^h} \;{\mathrm{d}}{{\mathscr {H}}}^{2} = \frac{\mathrm{d}}{{\mathrm{d}}t}\,{\mathscr {L}}^3(\varOmega ^h(t)), \end{aligned}$$
(3.23)

recall (1.5). Here \({\mathscr {V}}^h_{{\mathscr {S}}^h}(t)\) denotes the normal velocity of \({\mathscr {S}}^h(t)\) in the direction of \(\vec {\nu }^h_{{\mathscr {S}}^h}(t)\), the outer normal to \(\varOmega ^h(t)\) on \({\mathscr {S}}^h(t)\), where \(\vec {\nu }^h_{{\mathscr {S}}^h}(t)\) is induced by \(\vec {\nu }^h\) through a discrete analogue of (2.5). Using the same testing procedure for the scheme \(({\mathscr {C}}_h^{f,V})^h\) yields that

$$\begin{aligned} 0 = 2\,\pi \left( \vec {X}^h\cdot \,\vec {e}_1, \vec {X}^h_t\cdot \, \vec {\nu }^h\,|\vec {X}^h_\rho |\right) ^h, \end{aligned}$$
(3.24)

and so the enclosed volume is only approximately preserved, compare with (3.23). Finally, choosing \(\chi = \vec {X}^h\cdot \,\vec {e}_1\) in \(({\mathscr {A}}_h^{f,V})^h\), recall (3.21), also yields (3.24), and so an approximate volume preservation property.

4 Fully discrete schemes

In this section we present and analyse fully discrete variants of the approximations discussed in Sect. 3 that turn out to be most useful in practice. Other fully discrete schemes, which we will use for comparison purposes in the numerical results section, are summarized in “Appendix C”.

Let \(0= t_0< t_1< \cdots< t_{M-1} < t_M = T\) be a partitioning of [0, T] into possibly variable time steps \(\varDelta t_m = t_{m+1} - t_{m}\), \(m=0\rightarrow M-1\). We set \(\varDelta t= \max _{m=0\rightarrow M-1}\varDelta t_m\). For a given \(\vec {X}^m\in \underline{V}^h_{\partial _0}\), assuming that \(|\vec {X}^m_\rho | > 0\) almost everywhere on I, we set \(\vec {\nu }^m = - \frac{[\vec {X}^m_\rho ]^\perp }{|\vec {X}^m_\rho |}\). Let \(\vec {\omega }^m \in \underline{V}^h\) be the natural fully discrete analogue of \(\vec {\omega }^h \in \underline{V}^h\), recall (3.6), i.e.

$$\begin{aligned} \left( \vec {\omega }^m, \vec {\varphi } \, |\vec {X}^m_\rho | \right) ^h = \left( \vec {\nu }^m, \vec {\varphi } \, |\vec {X}^m_\rho | \right) = \left( \vec {\nu }^m, \vec {\varphi } \, |\vec {X}^m_\rho | \right) ^h \qquad \forall \ \vec {\varphi }\in {\underline{V}}^h. \end{aligned}$$
(4.1)

4.1 Mean curvature flow

Similarly to (3.9), and given a \(\kappa ^{m+1} \in V^h\), we introduce \({\mathfrak {K}}^{m}(\kappa ^{m+1}) \in V^h\) such that

$$\begin{aligned}{}[{\mathfrak {K}}^{m}(\kappa ^{m+1})](q_j) = {\left\{ \begin{array}{ll} \dfrac{\vec {\omega }^m(q_j)\cdot \,\vec {e}_1}{\vec {X}^m(q_j)\cdot \,\vec {e}_1} &{}\quad q_j \in {\overline{I}} {\setminus } \partial _0 I, \\ - \kappa ^{m+1}(q_j) &{}\quad q_j \in \partial _0 I. \end{array}\right. } \end{aligned}$$
(4.2)

Then our fully discrete analogue of \(({\mathscr {A}}_h)^h\), (3.10a, b), is given as follows.

\(({\mathscr {A}}_m)^h\) Let \(\vec {X}^0 \in \underline{V}^h_{\partial _0}\). For \(m=0,\ldots ,M-1\), find \((\delta \vec {X}^{m+1}, \kappa ^{m+1}) \in \underline{V}^h_\partial \times V^h\), where \(\vec {X}^{m+1} = \vec {X}^m + \delta \vec {X}^{m+1}\), such that

$$\begin{aligned} \left( \frac{\vec {X}^{m+1} - \vec {X}^m}{\varDelta t_m}, \chi \,\vec {\nu }^m\,|\vec {X}^m_\rho |\right) ^h =&\left( \kappa ^{m+1} - {\mathfrak {K}}^{m}(\kappa ^{m+1}), \chi \,|\vec {X}^m_\rho |\right) ^h \quad \forall \ \chi \in V^h, \end{aligned}$$
(4.3a)
$$\begin{aligned} \left( \kappa ^{m+1}\,\vec {\nu }^m, \vec {\eta }\,|\vec {X}^m_\rho |\right) ^h + \left( \vec {X}^{m+1}_\rho , \vec {\eta }_\rho \,|\vec {X}^m_\rho |^{-1}\right) =&- \sum _{i=1}^2 \sum _{p \in \partial _i I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\,\vec {\eta }(p)\cdot \,\vec {e}_{3-i} \nonumber \\&\qquad \forall \ \vec {\eta } \in \underline{V}^h_\partial . \end{aligned}$$
(4.3b)

We make the following assumptions.

figure b

Note that the assumption \(({\mathfrak {B}})^h\), on recalling (3.6), is equivalent to assuming that \(\dim \text {span}\{\vec {\omega }^m(q_j)\}_{j=0}^{J}= 2\). Clearly, it can only be violated if all the \(J+1\) vectors lie on a single line. It was shown in [6, Remark 2.2] that the assumption \(({\mathfrak {B}})^h\) always holds for closed curves without self-intersections.

Lemma 4.1

Let \(\vec {X}^m \in \underline{V}^h_{\partial _0}\) satisfy the assumptions \(({\mathfrak {A}})\) and \(({\mathfrak {B}})^h\). Then there exists a unique solution \((\delta \vec {X}^{m+1}, \kappa ^{m+1}) \in \underline{V}^h_\partial \times V^h\) to \(({\mathscr {A}}_m)^h\).

Proof

We note that since \(\vec {X}^m \in \underline{V}^h_{\partial _0}\) satisfies the assumption \(({\mathfrak {A}})\), the right hand side of (4.3a) is well-defined. As (4.3a, b) is linear, existence follows from uniqueness. To investigate the latter, we consider the system: find \((\delta \vec {X},\kappa ) \in \underline{V}^h_\partial \times V^h\) such that

$$\begin{aligned} \left( \frac{\delta \vec {X}}{\varDelta t_m}, \chi \,\vec {\nu }^m\,|\vec {X}^m_\rho |\right) ^h = \left( \lambda \,\kappa , \chi \,|\vec {X}^m_\rho |\right) ^h \qquad \forall \ \chi \in V^h, \end{aligned}$$
(4.4a)
$$\begin{aligned} \left( \kappa \,\vec {\nu }^m, \vec {\eta }\,|\vec {X}^m_\rho |\right) ^h + \left( (\delta \vec {X})_\rho , \vec {\eta }_\rho \,|\vec {X}^m_\rho |^{-1}\right) = 0 \qquad \forall \ \vec {\eta } \in \underline{V}^h_\partial , \end{aligned}$$
(4.4b)

where we recall from (4.2) that \(\lambda \in V^h\) with

$$\begin{aligned} \lambda (q_j) = {\left\{ \begin{array}{ll} 1 &{}\quad q_j \in {\overline{I}} {\setminus } \partial _0 I, \\ 2 &{}\quad q_j \in \partial _0 I. \end{array}\right. } \end{aligned}$$
(4.5)

Choosing \(\chi =\kappa \in V^h\) in (4.4a) and \(\vec {\eta }= \delta \vec {X} \in \underline{V}^h_\partial \) in (4.4b) yields that

$$\begin{aligned} \left( |(\delta \vec {X})_\rho |^2, |\vec {X}^m_\rho |^{-1}\right) +\varDelta t_m \left( \lambda \,|\kappa |^2, |\vec {X}^m_\rho |\right) ^h = 0. \end{aligned}$$
(4.6)

It follows from (4.6) that \(\kappa = 0\) and that \(\delta \vec {X} \equiv \vec {X}^c\in {{\mathbb {R}}}^2\); and hence that

$$\begin{aligned} 0 = \left( \vec {X}^c, \chi \,\vec {\nu }^m\,|\vec {X}^m_\rho |\right) ^h = \vec {X}^c \cdot \left( \vec {\nu }^m, \chi \,|\vec {X}^m_\rho |\right) ^h \quad \forall \ \chi \in V^h. \end{aligned}$$
(4.7)

It follows from (4.7) and assumption \(({\mathfrak {B}})^h\) that \(\vec {X}^c=\vec {0}\). Hence we have shown that (4.3a, b) has a unique solution \((\delta \vec {X}^{m+1},\kappa ^{m+1}) \in \underline{V}^h_\partial \times V^h\).

A fully discrete analogue of the scheme \(({\mathscr {C}}_h)^{(h)}\), (3.18a, b), is given as follows. Here we introduce the notation \([ r ]_\pm = \pm \max \{ \pm r, 0 \}\) for \(r \in {{\mathbb {R}}}\).

\(({\mathscr {C}}_{m,\star })^{(h)}\) Let \(\vec {X}^0 \in \underline{V}^h_{\partial _0}\). For \(m=0,\ldots ,M-1\), find \((\delta \vec {X}^{m+1}, \kappa _{{\mathscr {S}}}^{m+1}) \in \underline{V}^h_\partial \times W^h_{(\partial _0)}\), where \(\vec {X}^{m+1} = \vec {X}^m + \delta \vec {X}^{m+1}\), such that

$$\begin{aligned}&\left( \vec {X}^m\cdot \,\vec {e}_1\, \frac{\vec {X}^{m+1} - \vec {X}^m}{\varDelta t_m}, \chi \,\vec {\nu }^m\, |\vec {X}^m_\rho |\right) ^{(h)} = \left( \vec {X}^m\cdot \,\vec {e}_1\,\kappa _{{\mathscr {S}}}^{m+1}, \chi \,|\vec {X}^m_\rho |\right) ^{(h)} \quad \forall \ \chi \in W^h_{(\partial _0)}, \end{aligned}$$
(4.8a)
$$\begin{aligned}&\left( \vec {X}^m\cdot \,\vec {e}_1\,\kappa _{{\mathscr {S}}}^{m+1}\,\vec {\nu }^m, \vec {\eta }\,|\vec {X}^m_\rho |\right) ^{(h)} + \left( \vec {\eta } \cdot \,\vec {e}_1, |\vec {X}^{m+1}_\rho |\right) + \left( (\vec {X}^m\cdot \,\vec {e}_1)\, \vec {X}^{m+1}_\rho ,\vec {\eta }_\rho \, |\vec {X}^m_\rho |^{-1} \right) \nonumber \\&\quad = - \sum _{p \in \partial _1 I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\, (\vec {X}^m(p)\cdot \,\vec {e}_1)\,\vec {\eta }(p)\cdot \,\vec {e}_2 \nonumber \\&\qquad - \sum _{p \in \partial _2 I} ( ( [{\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}]_+ \, \,\vec {X}^{m+1}(p) + [{\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}]_- \, \,\vec {X}^{m}(p)) \cdot \,\vec {e}_1)\,\vec {\eta }(p)\cdot \,\vec {e}_1 \quad \forall \ \vec {\eta } \in \underline{V}^h_\partial . \end{aligned}$$
(4.8b)

We note that in practice we solve the nonlinear scheme \(({\mathscr {C}}_{m,\star })^{(h)}\) with a Newton iteration, which in all our experiments always converged with at most three iterations.

Remark 4.1

Similarly to the semidiscrete variants, we observe that it is possible to eliminate the discrete curvature, \(\kappa _{{\mathscr {S}}}^{m+1}\) from \(({\mathscr {C}}_{m,\star })^h\). On recalling (3.6) and on choosing \(\chi = \pi ^h[\vec {\eta }\cdot \,\vec {\omega }^m] \in W^h_{(\partial _0)}\) in (4.8a) for \(\vec {\eta } \in \underline{V}^h_\partial \), the scheme reduces to: Find \(\delta \vec {X}^{m+1}\in \underline{V}^h_\partial \) such that, for all \(\vec {\eta } \in \underline{V}^h_\partial \),

$$\begin{aligned}&\left( \vec {X}^m\cdot \,\vec {e}_1\, \frac{\vec {X}^{m+1} - \vec {X}^m}{\varDelta t_m}\cdot \,\vec {\omega }^m, \vec {\eta }\cdot \, \vec {\omega }^m\,|\vec {X}^m_\rho |\right) ^h + \left( \vec {\eta } \cdot \,\vec {e}_1, |\vec {X}^{m+1}_\rho |\right) \nonumber \\&\qquad +\,\left( (\vec {X}^m\cdot \,\vec {e}_1)\, \vec {X}^{m+1}_\rho ,\vec {\eta }_\rho \, |\vec {X}^m_\rho |^{-1} \right) = - \sum _{i=1}^2 \sum _{p \in \partial _i I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\, (\vec {X}^{m+1}(p)\cdot \,\vec {e}_1)\, \vec {\eta }(p)\cdot \,\vec {e}_{3-i}. \end{aligned}$$
(4.9)

For the scheme \(({\mathscr {C}}_{m,\star })\) this elimination procedure is not possible. A related variant to (4.9) is given by: find \(\delta \vec {X}^{m+1}\in \underline{V}^h_\partial \) such that, for all \(\vec {\eta } \in \underline{V}^h_\partial \),

$$\begin{aligned}&\left( \vec {X}^m\cdot \,\vec {e}_1\, \frac{\vec {X}^{m+1} - \vec {X}^m}{\varDelta t_m}\cdot \,\vec {\nu }^m, \vec {\eta }\cdot \, \vec {\nu }^m\,|\vec {X}^m_\rho |\right) ^h + \left( \vec {\eta }\cdot \,\vec {e}_1, |\vec {X}^{m+1}_\rho |\right) \nonumber \\&\quad +\,\left( (\vec {X}^m\cdot \,\vec {e}_1)\, \vec {X}^{m+1}_\rho ,\vec {\eta }_\rho \, |\vec {X}^m_\rho |^{-1} \right) = - \sum _{i=1}^2 \sum _{p \in \partial _i I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\, (\vec {X}^{m+1}(p)\cdot \,\vec {e}_1)\,\vec {\eta }(p)\cdot \,\vec {e}_{3-i}. \end{aligned}$$
(4.10)

Theorem 4.1

Let \(\vec {X}^m \in \underline{V}^h_{\partial _0}\) satisfy the assumption \((\mathfrak A)\). Let \((\vec {X}^{m+1},\kappa _{{\mathscr {S}}}^{m+1})\) be a solution to \(({\mathscr {C}}_{m,\star })^{(h)}\). Then it holds that

$$\begin{aligned} E(\vec {X}^{m+1}) + 2\,\pi \,\varDelta t_m \left( \vec {X}^m\cdot \,\vec {e}_1\,|\kappa _{{\mathscr {S}}}^{m+1}|^2, |\vec {X}^m_\rho |\right) ^{(h)} \le E(\vec {X}^m), \end{aligned}$$
(4.11)

where we recall the definition (3.7).

Proof

Choosing \(\chi = \varDelta t_m\,\kappa _{{\mathscr {S}}}^{m+1}\) in (4.8a) and \(\vec {\eta } = \vec {X}^{m+1} - \vec {X}^m \in \underline{V}^h_\partial \) in (4.8b) yields, on noting that \(\vec {X}^{m}(p)\cdot \,\vec {e}_1 = \vec {X}^{m+1}(p)\cdot \,\vec {e}_1\) for \(p \in \partial _1 I\), that

$$\begin{aligned}&- \varDelta t_m \left( \vec {X}^m\cdot \,\vec {e}_1\,|\kappa _{{\mathscr {S}}}^{m+1}|^2, |\vec {X}^m_\rho |\right) ^{(h)} \nonumber \\&\qquad = \left( \vec {X}^{m+1} - \vec {X}^m,\vec {e}_1\, |\vec {X}^{m+1}_\rho |\right) + \left( (\vec {X}^m\cdot \,\vec {e}_1)\, (\vec {X}^{m+1} - \vec {X}^m)_\rho ,\vec {X}^{m+1}_\rho \, |\vec {X}^m_\rho |^{-1} \right) \nonumber \\&\qquad \qquad + \sum _{p \in \partial _1 I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\, (\vec {X}^{m}(p)\cdot \,\vec {e}_1)\,(\vec {X}^{m+1}(p) - \vec {X}^m(p))\cdot \,\vec {e}_2 \nonumber \\&\qquad \qquad + \sum _{p \in \partial _2 I} ([{\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}]_+\,\vec {X}^{m+1}(p) + [{\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}]_-\,\vec {X}^{m}(p)]\cdot \,\vec {e}_1)\, (\vec {X}^{m+1}(p) - \vec {X}^m(p))\cdot \,\vec {e}_1 \nonumber \\&\qquad \ge \left( \vec {X}^{m+1} - \vec {X}^m,\vec {e}_1\, |\vec {X}^{m+1}_\rho |\right) + \left( \vec {X}^m\cdot \,\vec {e}_1, |\vec {X}^{m+1}_\rho | - |\vec {X}^m_\rho | \right) \nonumber \\&\qquad \qquad + \sum _{p \in \partial _1 I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\, (\vec {X}^{m}(p)\cdot \,\vec {e}_1)\,\vec {X}^{m+1}(p) \cdot \,\vec {e}_2 - \sum _{p \in \partial _1 I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\, (\vec {X}^{m}(p)\cdot \,\vec {e}_1)\,\vec {X}^{m}(p) \cdot \,\vec {e}_2 \nonumber \\&\qquad \qquad + \tfrac{1}{2} \sum _{p \in \partial _2 I} [{\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}]_+\, (\vec {X}^{m+1}(p)\cdot \,\vec {e}_1)^2\, - \tfrac{1}{2} \sum _{p \in \partial _2 I} [{\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}]_+\, (\vec {X}^{m}(p)\cdot \,\vec {e}_1)^2\, \nonumber \\&\qquad \qquad + \tfrac{1}{2} \sum _{p \in \partial _2 I} [{\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}]_-\, (\vec {X}^{m+1}(p)\cdot \,\vec {e}_1)^2\, - \tfrac{1}{2} \sum _{p \in \partial _2 I} [{\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}]_-\, (\vec {X}^{m}(p)\cdot \,\vec {e}_1)^2\, \nonumber \\&\qquad = \left( \vec {X}^{m+1} \cdot \,\vec {e}_1, |\vec {X}^{m+1}_\rho |\right) - \left( \vec {X}^m\cdot \,\vec {e}_1, |\vec {X}^m_\rho | \right) \nonumber \\&\qquad \qquad + \sum _{p \in \partial _1 I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\, (\vec {X}^{m+1}(p)\cdot \,\vec {e}_1)\,\vec {X}^{m+1}(p) \cdot \,\vec {e}_2 - \sum _{p \in \partial _1 I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\, (\vec {X}^{m}(p)\cdot \,\vec {e}_1)\,\vec {X}^{m}(p) \cdot \,\vec {e}_2 \nonumber \\&\qquad \qquad + \tfrac{1}{2} \sum _{p \in \partial _2 I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\, (\vec {X}^{m+1}(p)\cdot \,\vec {e}_1)^2\, - \tfrac{1}{2} \sum _{p \in \partial _2 I} {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(p)}\, (\vec {X}^{m}(p)\cdot \,\vec {e}_1)^2\, \nonumber \\&\qquad = \frac{1}{2\,\pi }\,E(\vec {X}^{m+1}) - \frac{1}{2\,\pi }\,E(\vec {X}^m), \end{aligned}$$
(4.12)

where we have used the two inequalities \(\vec {a}\cdot \,(\vec {a} - \vec {b}) \ge |\vec {b}|\,(|\vec {a}| - |\vec {b}|)\) for \(\vec {a}\), \(\vec {b} \in {{\mathbb {R}}}^2\), and \(2\,\gamma \,(\gamma - \alpha ) \ge \gamma ^2 - \alpha ^2\) for \(\alpha ,\gamma \in {{\mathbb {R}}}\). This proves the desired result.

We note that the scheme (4.10) can also be shown to be unconditionally stable, i.e. a solution to (4.10) satisfies

$$\begin{aligned} E(\vec {X}^{m+1}) + 2\,\pi \,\varDelta t_m \left( \vec {X}^m \cdot \,\vec {e}_1 \left| \dfrac{\vec {X}^{m+1} - \vec {X}^m}{\varDelta t_m}\cdot \,\vec {\nu }^m\right| ^2, |\vec {X}^m_\rho | \right) ^h \le E(\vec {X}^m). \end{aligned}$$

4.2 Nonlinear mean curvature flow

It is a simple matter to extend the presented fully discrete approximations to the nonlinear mean curvature flow (1.2) and the volume preserving variant (1.3). We recall the fully 3d parametric finite element schemes (2.4a, b) and (2.5) in [7] for the approximation of (1.2) and (1.3), respectively. We can now define their natural axisymmetric analogues. For example, the natural adaptation \(({\mathscr {A}}^f_{m})^{h}\) of the scheme \(({\mathscr {A}}_{m})^{h}\) to (1.2) is given by (4.3a, b) with (4.3a) replaced by

$$\begin{aligned} \left( \frac{\vec {X}^{m+1} - \vec {X}^m}{\varDelta t_m}, \chi \,\vec {\nu }^m\,|\vec {X}^m_\rho |\right) ^h = \left( f(\kappa ^{m+1} - {\mathfrak {K}}^{m}(\kappa ^{m+1})), \chi \,|\vec {X}^m_\rho |\right) ^h \qquad \forall \ \chi \in V^h. \end{aligned}$$
(4.13)

Similarly, the natural adaptation \(({\mathscr {C}}^f_{m,\star })^{(h)}\) of the scheme \(({\mathscr {C}}_{m,\star })^{(h)}\) to (1.2) is given by (4.8a, b) with (4.8a) replaced by

$$\begin{aligned} \left( \vec {X}^m\cdot \,\vec {e}_1\, \frac{\vec {X}^{m+1} - \vec {X}^m}{\varDelta t_m}, \chi \,\vec {\nu }^m\, |\vec {X}^m_\rho |\right) ^{(h)} =&\left( \vec {X}^m\cdot \,\vec {e}_1\,\pi ^h[f(\kappa _{{\mathscr {S}}}^{m+1})], \chi \,|\vec {X}^m_\rho |\right) ^{(h)} \nonumber \\&\qquad \forall \ \chi \in W^h_{(\partial _0)}. \end{aligned}$$
(4.14)

Similarly to (4.11), with the same choices of \(\chi \) and \(\vec {\eta }\), it is then possible to prove that solutions to \(({\mathscr {C}}^f_{m,\star })^{h}\) satisfy \(E(\vec {X}^{m+1}) + 2\,\pi \,\varDelta t_m\, \left( \vec {X}^m\cdot \,\vec {e}_1\, f( \kappa _{{\mathscr {S}}}^{m+1}), \kappa _{{\mathscr {S}}}^{m+1} \,|\vec {X}^m_\rho | \right) ^{h} \le E(\vec {X}^m)\), which provides a stability bound if f is monotonically increasing with \(f(0)=0\).

Finally, replacing (4.13) with

$$\begin{aligned} \left( \frac{\vec {X}^{m+1} - \vec {X}^m}{\varDelta t_m}, \chi \,\vec {\nu }^m\,|\vec {X}^m_\rho |\right) ^h&= \left( f(\kappa ^{m+1} - {\mathfrak {K}}^{m}(\kappa ^{m+1})), \chi \,|\vec {X}^m_\rho |\right) ^h \nonumber \\&\quad -\,\frac{\left( \vec {X}^m\cdot \,\vec {e}_1, f(\kappa ^m - {\mathfrak {K}}^{m}(\kappa ^m))\, |\vec {X}^m_\rho |\right) ^{h}}{\left( \vec {X}^m\cdot \,\vec {e}_1, |\vec {X}^m_\rho |\right) } \left( \chi ,|\vec {X}^m_\rho |\right) ^{h} \quad \forall \ \chi \in V^h \end{aligned}$$
(4.15)

and replacing (4.14) with

$$\begin{aligned}&\left( \vec {X}^m\cdot \,\vec {e}_1\, \frac{\vec {X}^{m+1} - \vec {X}^m}{\varDelta t_m}, \chi \,\vec {\nu }^m\, |\vec {X}^m_\rho |\right) ^{(h)} \nonumber \\&\quad = \left( \vec {X}^m\cdot \,\vec {e}_1\,\pi ^h[f(\kappa _{{\mathscr {S}}}^{m+1})], \chi \,|\vec {X}^m_\rho |\right) ^{(h)} \nonumber \\&\qquad - \frac{\left( \vec {X}^m\cdot \,\vec {e}_1, \pi ^h[f(\kappa _{{\mathscr {S}}}^{m+1})]\, |\vec {X}^m_\rho |\right) ^{(h)}}{\left( \vec {X}^m\cdot \,\vec {e}_1, |\vec {X}^m_\rho |\right) } \left( \vec {X}^m\cdot \,\vec {e}_1, \chi \,|\vec {X}^m_\rho |\right) ^{(h)} \quad \forall \ \chi \in W^h_{(\partial _0)}\end{aligned}$$
(4.16)

gives the fully discrete approximations \(({\mathscr {A}}^{f,V}_{m})^{h}\) and \(({\mathscr {C}}^{f,V}_{m,\star })^{(h)}\), respectively, of (1.3). For the case (1.4c) it is possible to prove a stability bound for \(({\mathscr {C}}^{f,V}_{m,\star })^{(h)}\). In particular, solutions to \(({\mathscr {C}}^{f,V}_{m,\star })^{(h)}\) satisfy, similarly to (4.12), that

$$\begin{aligned}&\frac{1}{2\,\pi }\,E(\vec {X}^{m}) - \frac{1}{2\,\pi }\,E(\vec {X}^{m+1}) \nonumber \\&\quad \ge \varDelta t_m \left( \vec {X}^m\cdot \,\vec {e}_1\,|\kappa _{{\mathscr {S}}}^{m+1}|^2, |\vec {X}^m_\rho |\right) ^{(h)} - \left[ \left( \vec {X}^m\cdot \,\vec {e}_1,|\vec {X}^m_\rho | \right) \right] ^{-1} \left| \left( \vec {X}^m\cdot \,\vec {e}_1,\kappa _{{\mathscr {S}}}^{m+1}\, |\vec {X}^m_\rho | \right) ^{(h)}\right| ^2 \nonumber \\&\quad \ge 0, \end{aligned}$$

where we have noted the Cauchy–Schwarz inequality. For the fully discrete approximations \(({\mathscr {A}}^{f,V}_{m})^{h}\) and \(({\mathscr {C}}^{f,V}_{m,\star })^{(h)}\) it is not possible to prove a volume conservation property. However, in practice all three schemes preserve the enclosed volume well, with the relative volume loss decreasing as the discretization parameters become smaller.

Finally, we note that for the schemes \(({\mathscr {A}}_m^f)^h\) and \(({\mathscr {A}}_m^{f,V})^h\), depending on the choice of f, existence and uniqueness results can be shown, see “Appendix B”.

Remark 4.2

In order to be able to compute evolutions for the general flow (1.6), we propose the scheme \(({\mathscr {A}}^F_{m})^{h}\), which can be obtained from the scheme \(({\mathscr {A}}_{m})^{h}\) (4.3a, 4.3b), by replacing (4.3a) with

$$\begin{aligned} \left( \frac{\vec {X}^{m+1} - \vec {X}^m}{\varDelta t_m}, \chi \,\vec {\nu }^m\,|\vec {X}^m_\rho |\right) ^h= & {} \left( F(\kappa ^{m} - {\mathfrak {K}}^{m}(\kappa ^{m}), -\kappa ^{m}\,{\mathfrak {K}}^{m}(\kappa ^m)), \chi \,|\vec {X}^m_\rho |\right) ^h \nonumber \\&\qquad \qquad \qquad \qquad \qquad \quad \qquad \qquad \forall \ \chi \in V^h. \end{aligned}$$
(4.17)

This is linear scheme for which existence of a unique solution, provided that the assumption \(({\mathfrak {B}})^h\) holds, can easily be shown. Moreover, solutions to the semidiscrete variant of \(({\mathscr {A}}^F_{m})^{h}\) satisfy the equidistribution property (3.11).

5 Numerical results

As the fully discrete energy, we consider \(E(\vec {X}^m)\), recall (3.7). We always employ uniform time steps, \(\varDelta t_m = \varDelta t\), \(m=0,\ldots ,M-1\).

We also consider the ratio

$$\begin{aligned} {{\mathfrak {r}}}^m = \dfrac{\max _{j=1\rightarrow J} |\vec {X}^m(q_j) - \vec {X}^m(q_{j-1})|}{\min _{j=1\rightarrow J} |\vec {X}^m(q_j) - \vec {X}^m(q_{j-1})|} \end{aligned}$$
(5.1)

between the longest and shortest element of \(\varGamma ^m\), and are often interested in the evolution of this ratio over time.

On recalling (2.7), and given \(\varGamma ^0 = \vec {X}^0({\overline{I}})\), for the scheme \(({\mathscr {A}}_m^{f,V})^h\) we define \(\kappa ^0 \in V^h\) via

$$\begin{aligned} \kappa ^0 = \pi ^h\left[ \frac{\vec {\kappa }^0\cdot \,\vec {\omega }^0}{|\vec {\omega }^0|} \right] , \end{aligned}$$

recall (4.1), where \(\vec {\kappa }^0\in \underline{V}^h\) is such that

$$\begin{aligned} \left( \vec {\kappa }^{0},\vec {\eta }\, |\vec {X}^0_\rho | \right) ^h + \left( \vec {X}^{0}_\rho , \vec {\eta }_\rho \,|\vec {X}^0_\rho |^{-1} \right) = 0 \quad \forall \ \vec {\eta } \in \underline{V}^h. \end{aligned}$$
Table 3 Errors for the convergence test for (5.2) with \(r_0 = 1\) over the time interval [0, 0.125]

5.1 Numerical results for mean curvature flow

5.1.1 Sphere

It is easy to show that a sphere of radius r(t), with

$$\begin{aligned} r(t) = [r^2(0) - 4\,t]^\frac{1}{2}, \end{aligned}$$
(5.2)

is a solution to (1.1). We use this true solution for a convergence test for the various schemes for mean curvature flow, similarly to Table 1 in [7]. Here we start with a nonuniform partitioning of a semicircle of radius \(r(0)=r_0=1\) and compute the flow until time \(T = 0.125\). In particular, we have \(\partial _0 I = \partial I = \{0,1\}\) and we choose \(\vec {X}^0 \in \underline{V}^h_{\partial _0}\) with

$$\begin{aligned} \vec {X}^0(q_j) = r_0 \begin{pmatrix} \cos [(q_j-\tfrac{1}{2})\,\pi + 0.1\,\cos ((q_j-\tfrac{1}{2})\,\pi )] \\ \sin [(q_j-\tfrac{1}{2})\,\pi + 0.1\,\cos ((q_j-\tfrac{1}{2})\,\pi )] \end{pmatrix}, \quad j = 0,\ldots ,J, \end{aligned}$$
(5.3)

recall (3.1). We compute the error

$$\begin{aligned} \Vert \varGamma - \varGamma ^h\Vert _{L^\infty }= \max _{m=1,\ldots ,M} \max _{j=0,\ldots ,J} | |\vec {X}^m(q_j)| - r(t_m)| \end{aligned}$$
(5.4)

over the time interval [0, T] between the true solution (5.2) and the discrete solutions for the schemes \(({\mathscr {A}}_m)^h\), \(({\mathscr {B}}_m)^h\), \(({\mathscr {C}}_m)^{(h)}\), \(({\mathscr {D}}_m)^{(h)}\), \(({\mathscr {C}}_{m,\star })^{(h)}\) and \(({\mathscr {D}}_{m,\star })^{(h)}\). Here we used the time step size \(\varDelta t=0.1\,h^2_{\varGamma ^0}\), where \(h_{\varGamma ^0}\) is the maximal edge length of \(\varGamma ^0\). The computed errors are reported in Tables 3, 4 and 5. Comparing the reported numbers with the values in Table 1 in [7], we see that the errors for the axisymmetric schemes are significantly smaller than the values in Table 1 in [7] for similar discretization parameters. It is clear from Table 3 that the schemes \(({\mathscr {A}}_m)^h\) and \(({\mathscr {B}}_m)^h\) appear to converge with the optimal convergence rate of \({\mathscr {O}}(h^2_{\varGamma ^0})\). Similarly, Tables 4 and 5 suggest that the schemes \(({\mathscr {C}}_{m,\star })^{(h)}\), \(({\mathscr {D}}_m)^{(h)}\) and \(({\mathscr {D}}_{m,\star })^{(h)}\) converge with an order slightly less than quadratic. We note that the linear schemes \(({\mathscr {C}}_m)^{(h)}\) lead to solutions \(\vec {X}^m\) with \(\min _{\rho \in {\overline{I}}} \vec {X}^m(\rho )\cdot \,\vec {e}_1 < 0\) and so we cannot complete the evolutions. In particular, in practice the two boundary elements shrink in size due to the scheme’s tangential motion. Once the element has shrunk to a length almost zero, the freely moving vertex can become negative. The linear schemes \(({\mathscr {D}}_m)^{(h)}\) behave well, on the other hand. But as they are very close to the nonlinear schemes \(({\mathscr {D}}_{m,\star })^{(h)}\), from now on we concentrate on the schemes \(({\mathscr {D}}_{m,\star })^{(h)}\), \(({\mathscr {C}}_{m,\star })^{(h)}\), \(({\mathscr {A}}_{m})^{h}\) and \(({\mathscr {B}}_{m})^{h}\).

Table 4 Errors for the convergence test for (5.2) with \(r_0 = 1\) over the time interval [0, 0.125]
Table 5 Errors for the convergence test for (5.2) with \(r_0 = 1\) over the time interval [0, 0.125]
Fig. 2
figure 2

\(({\mathscr {A}}_m)^h\) Evolution for a torus with radii \(R=1\), \(r=0.7\). Plots are at times \(t=0,0.01,\ldots ,0.08\). We also show a plot at time \(t=0.082\), together with a plot of the discrete energy. Below we visualize the axisymmetric surface \({\mathscr {S}}^m\) generated by \(\varGamma ^m\) at times \(t=0\) and \(t=0.082\)

5.1.2 Torus

We repeat the two torus experiments in Figures 5 and 6 in [7]. To this end, we let \(\partial I = \emptyset \). For an initial torus with radii \(R=1\), \(r=0.7\), we obtain a surface that closes up towards a genus-0 surface, as in [7, Fig. 5]. See Fig. 2 for the simulation results for the scheme \(({\mathscr {A}}_m)^h\), for the discretization parameters \(J=256\) and \(\varDelta t= 10^{-4}\). On the other hand, for an initial torus with radii \(R=1\), \(r=0.5\), we obtain a shrinking evolution towards a circle, as in [7, Fig. 6]. See Fig. 3 for the evolution for the scheme \(({\mathscr {A}}_m)^h\), again for \(J=256\) and \(\varDelta t= 10^{-4}\). On repeating the numerical experiment for the \(({\mathscr {B}}_m)^h\), we observe strong oscillations, as shown on the right of Fig. 3. These oscillations become smaller in magnitude as \(\varDelta t\) is decreased. The remaining schemes can integrate the evolution shown in Fig. 3 in a stable way, and their numerical results are very close to the ones displayed in Fig. 3 for the scheme \(({\mathscr {A}}_m)^h\). However, the schemes differ in the exhibited tangential motions, which leads to diverse evolutions of the ratio \({{\mathfrak {r}}}^m\), see Fig. 4. The best distribution of mesh points is shown by the schemes \(({\mathscr {A}}_m)^h\) and \(({\mathscr {C}}_{m,\star })^h\), followed by \(({\mathscr {C}}_{m,\star })\). The most nonuniform distribution of mesh points can be observed for the two schemes \(({\mathscr {D}}_{m,\star })^{(h)}\).

Fig. 3
figure 3

\(({\mathscr {A}}_m)^h\) Evolution for a torus with radii \(R=1\), \(r=0.5\). Plots are at times \(t=0,0.01,\ldots ,0.13\). We also show a plot of the discrete energy and, below, we visualize the axisymmetric surface \({\mathscr {S}}^m\) generated by \(\varGamma ^m\) at times \(t=0\) and \(t=0.135\). On the top far right, the evolution for the scheme \(({\mathscr {B}}_m)^h\), with plots at times \(t=0,0.01,\ldots ,0.13\)

Fig. 4
figure 4

Plots of the ratio \({{\mathfrak {r}}}^m\) for the schemes \(({\mathscr {A}}_m)^h\), \(({\mathscr {C}}_{m,\star })^h\), \(({\mathscr {C}}_{m,\star })\), \(({\mathscr {D}}_{m,\star })^h\), \(({\mathscr {D}}_{m,\star })\)

5.1.3 Cylinder

For the scheme \(({\mathscr {A}}_m)^h\) we repeat the singular evolution from [22, Fig. 1]. To this end, we set \(\partial _D I = \partial I = \{0,1\}\). In particular, starting with a cylinder, mean curvature flow leads to a pinch-off. We show the results for the scheme \(({\mathscr {A}}_m)^h\), with the discretization parameters \(J=128\) and \(\varDelta t= 10^{-4}\), in Fig. 5.

Fig. 5
figure 5

\(({\mathscr {A}}_m)^h\)\([\partial _D I = \partial I = \{0,1\}]\) Evolution for a cylinder with fixed boundary. Plots are at times \(t=0,0.1,\ldots ,0.5\). We also show a plot at time \(t=0.51\), as well as a plot of the discrete energy. On the right we visualize the axisymmetric surface \({\mathscr {S}}^m\) generated by \(\varGamma ^m\) at times \(t=0\) and \(t=0.51\)

For the next two experiments, we consider a cylinder attached to two parallel hyperplanes, with prescribed contact angle conditions, recall (2.10b). To this end, we set \(\partial _2 I = \partial I = \{0,1\}\) and use the discretization parameters \(J=128\) and \(\varDelta t= 10^{-3}\). Letting \({\widehat{\rho }}_{\partial {\mathscr {S}}}^{(0)} = {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(1)} = -\frac{1}{2}\) and starting with a cylinder, the evolution yields a growing catenoid-like surface, see Fig. 6. We observe convergence to a travelling wave type solution, with the associated energy unbounded from below. We conjecture that the profile of the curve approaches in the limit the so-called grim reaper solution, see [37, p. 15] and [28],

$$\begin{aligned} \vec {g} (\rho ,t) = (z_0 + \tfrac{\pi }{3}\,t)\,\vec {e}_1 + \left( -\tfrac{3}{\pi }\,\ln \cos \left( \tfrac{\pi }{3}\,(\rho -\tfrac{1}{2})\right) , \rho \right) ^T, \end{aligned}$$
(5.5)

where \(z_0 \in {{\mathbb {R}}}\) specifies the position of the travelling wave solution at time \(t=0\). In fact, plotting \(\vec {g} (\rho ,t) - (z_0 + \tfrac{\pi }{3}\,t)\,\vec {e}_1\) at time \(t=4\) versus \(\vec {X}^m - (\min _{\rho \in {\overline{I}}} \vec {X}^m\cdot \,\vec {e}_1)\,\vec {e}_1\) for our final solution in Fig. 6, yields perfect agreement between the two graphs. We conjecture that the speed of the travelling wave type solution will approach \(\frac{\pi }{3}\), the speed of (5.5). To test this conjecture, we continue the evolution until \(t=100\) and plot the evolution of \(\vec {X}^m(0)\cdot \,\vec {e}_1\) over time, comparing the graph with a suitably chosen line with slope \(\frac{\pi }{3}\), see Fig. 7. As we can see, the speed of the curve does indeed approach \(\frac{\pi }{3}\).

Fig. 6
figure 6

\(({\mathscr {A}}_m)^h \left[ \partial _2 I = \partial I = \{0,1\}, {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(0)} = {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(1)} = -\frac{1}{2}\right] \) Evolution for an open cylinder attached to \({{\mathbb {R}}}\times \{ 0 \} \times {{\mathbb {R}}}\) and \({{\mathbb {R}}}\times \{ 1 \} \times {{\mathbb {R}}}\). Solution at times \(t=0,0.5,\ldots ,4\), as well as a plot of the discrete energy over time. We also visualize the axisymmetric surface \({\mathscr {S}}^m\) generated by \(\varGamma ^m\) at time \(t=4\)

Fig. 7
figure 7

We plot \(\vec {X}^m(0)\cdot \,\vec {e}_1\) over time, compared with the linear function \(t \mapsto \frac{\pi }{3}\,t - \frac{13}{2}\), for the evolution in Fig. 6 over the larger time interval [0, 100]

If we let \({\widehat{\rho }}_{\partial {\mathscr {S}}}^{(0)} = {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(1)} = \frac{1}{2}\), on the other hand, we observe a shrinking surface, with the radius of the contact circles eventually converging to zero. On reaching two single contact points with the external substrates, we allow the discrete surface to detach from the two hyperplanes and to continue the evolution as a closed genus 0 surface, see Fig. 8 for the evolution. To allow for an accurate resolution of the detaching, we employ the smaller time step size \(\varDelta t=10^{-6}\) for this simulation.

Fig. 8
figure 8

\(({\mathscr {A}}_m)^h \left[ \partial _2 I = \partial I = \{0,1\}, {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(0)} = {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(1)} = \frac{1}{2}\right] \) Evolution for an open cylinder attached to \({{\mathbb {R}}}\times \{ 0 \} \times {{\mathbb {R}}}\) and \({{\mathbb {R}}}\times \{ 1 \} \times {{\mathbb {R}}}\). Solution at times \(t=0,0.1,\ldots , 1.5,1.53,1.54,1.55,1.56\), as well as a plot of the discrete energy over time. We also visualize the axisymmetric surface \({\mathscr {S}}^m\) generated by \(\varGamma ^m\) at times \(t=1.5\), \(t=1.55\) and \(t=1.56\)

5.1.4 Surface patch within a cylinder

For the next experiment, we consider a disk attached to an infinite cylinder of radius 1, with prescribed contact angle conditions, recall (2.10a). To this end, we set \(\partial _0 I = \{0\}\) and \(\partial _1 I = \{1\}\), and use the discretization parameters \(J=128\) and \(\varDelta t= 10^{-3}\). Letting \({\widehat{\rho }}_{\partial {\mathscr {S}}}^{(1)} = -\frac{1}{2}\) and starting with a disk, the evolution seems to converge to a translating surface patch, see Fig. 9. Taking the angle condition (2.10a), there is a unique convex scaled surface grim reaper profile moving with constant speed by translation. We conjecture that a general class of initial data will converge to this shape for large times. We refer to [1] for more information on the grim reaper analogues in higher dimensions.

5.2 Numerical results for conserved mean curvature flow

5.2.1 Sphere

Clearly, a sphere is a stationary solution for conserved mean curvature flow, (1.3) with (1.4c). Hence, setting \(\partial _0 I = \partial I = \{0,1\}\) and choosing as initial data the nonuniform approximation of a semicircle (5.3) with \(J=64\), we now investigate the different tangential motions exhibited by our proposed schemes. The initial data \(\vec {X}^0\) has a ratio \({{\mathfrak {r}}}^0=1.22\), recall (5.1). We set \(\varDelta t= 10^{-4}\) and integrate the evolution until time \(T=1\). For the three schemes \(({\mathscr {A}}_m)^h\), \(({\mathscr {C}}_{m,\star })^h\) and \(({\mathscr {C}}_{m,\star })\) the element ratios \({{\mathfrak {r}}}^m\) at time \(T=1\) are 1.01, 73.13, 2.94, and the enclosed volume is preserved almost exactly by all the schemes. We show the final distributions of vertices, and plots of \({{\mathfrak {r}}}^m\) over time in Fig. 10. An insight that we gain from this set of experiments is that the tangential motion displayed by the scheme \(({\mathscr {C}}_{m,\star })^h\) can lead to very nonuniform meshes. Hence, for the remainder of this paper, we will only present numerical results for the two schemes \(({\mathscr {A}}_m)^h\) and \(({\mathscr {C}}_{m,\star })\) and their nonlinear variants. Note that the former is a linear fully discrete approximation of \(({\mathscr {A}}_h)^h\), for which the equidistribution property (3.11) holds. The latter, on the other hand, is a nonlinear scheme that is unconditionally stable, recall Theorem 4.1. As the results for \(({\mathscr {A}}_m)^h\) and \(({\mathscr {C}}_{m,\star })\) are often indistinguishable, we only visualize the numerical results for the former from now on.

Fig. 9
figure 9

\(({\mathscr {A}}_m)^h \left[ \partial _0 I = \{0\}, \partial _1 I = \{1\}, {\widehat{\rho }}_{\partial {\mathscr {S}}}^{(1)} = -\frac{1}{2}\right] \) Evolution for a disk attached to an infinite cylinder of radius 1. Solution at times \(t=0,0.5,\ldots ,2\), as well as a plot of the discrete energy over time. We also visualize the axisymmetric surface \({\mathscr {S}}^m\) generated by \(\varGamma ^m\) at time \(t=2\)

Fig. 10
figure 10

Comparison of the different schemes for conserved mean curvature flow, (1.3) with (1.4c), of the unit sphere. Left to right: \(({\mathscr {A}}_m^{f,V})^h\), \(({\mathscr {C}}_{m,\star }^{f,V})^h\) and \(({\mathscr {C}}_{m,\star }^{f,V})\). Plots are for \(\vec {X}^m\) at time \(t=1\) and for the ratio \({{\mathfrak {r}}}^m\) over time. The element ratios \({{\mathfrak {r}}}^m\) at time \(t=1\) are 1.01, 73.13 and 2.94, respectively

5.2.2 Genus 0 surface

An experiment for a cigar shape can be seen in Fig. 11. Here we have once again that \(\partial _0 I = \partial I = \{0,1\}\). The discretization parameters are \(J=128\) and \(\varDelta t= 10^{-4}\). The relative volume loss for this experiment for \(({\mathscr {A}}_m^{f,V})^h\) is \(0.09\%\), while for \(({\mathscr {C}}_{m,\star }^{f,V})\) it is \(-\,0.01\%\). The same experiment for the scheme \(({\mathscr {C}}_{m,\star }^{f,V})^h\) yields a very nonuniform mesh, with the final ratio \({{\mathfrak {r}}}^m > 430\).

Fig. 11
figure 11

\(({\mathscr {A}}_m^{f,V})^h\) for (1.4c). Conserved mean curvature flow for a cigar. Plots are at times \(t=0,0.1,\ldots ,1\). We also visualize the axisymmetric surface \({\mathscr {S}}^m\) generated by \(\varGamma ^m\) at time \(t=0.3\). On the right are plots of the discrete energy and the ratio \({{\mathfrak {r}}}^m\) and, as a comparison, a plot of the ratio \({{\mathfrak {r}}}^m\) for the scheme \(({\mathscr {C}}_{m,\star }^{f,V})\)

An experiment for a disc shape is shown in Fig. 12. The discretization parameters are \(J=128\) and \(\varDelta t= 10^{-4}\). The relative volume loss for this experiment for \(({\mathscr {A}}_m^{f,V})^h\) is \(-\,0.02\%\), while for \(({\mathscr {C}}_{m,\star }^{f,V})\) it is \(-\,0.01\%\). Once again, the scheme \(({\mathscr {C}}_{m,\star }^{f,V})^h\) yields a very nonuniform mesh for this simulation, with the final ratio \({{\mathfrak {r}}}^m > 145\).

Fig. 12
figure 12

\(({\mathscr {A}}_m^{f,V})^h\) for (1.4c). Conserved mean curvature flow for a disc. Plots are at times \(t=0,0.1,\ldots ,4\). We also visualize the axisymmetric surface \({\mathscr {S}}^m\) generated by \(\varGamma ^m\) at time \(t=0.5\). On the right are plots of the discrete energy and the ratio \({{\mathfrak {r}}}^m\) and, as a comparison, a plot of the ratio \({{\mathfrak {r}}}^m\) for the scheme \(({\mathscr {C}}_{m,\star }^{f,V})\)

5.2.3 Genus 1 surface

We repeat the simulation in Fig. 3 for conserved mean curvature flow, i.e. (1.3) with (1.4c), using the scheme \(({\mathscr {A}}_m^{f,V})^h\). Conservation of the enclosed volume means that the torus can no longer shrink to a circle. Hence the torus now attempts to close up and change topology, as can be seen from the numerical results in Fig. 13. As for the original experiment, we use the discretization parameters \(J=256\) and \(\varDelta t= 10^{-4}\). The relative enclose volume loss for this experiment is \(-\,0.00\%\). The evolutions for the schemes \(({\mathscr {C}}_{m,\star }^{f,V})^{h}\) and \(({\mathscr {C}}_{m,\star }^{f,V})\) are nearly identical to what is shown in Fig. 13, with a relative volume loss of \(0.01\%\) in both cases.

Fig. 13
figure 13

\(({\mathscr {A}}_m^{f,V})^h\) for (1.4c). Conserved mean curvature flow for a torus with radii \(R=1\), \(r=0.5\). Plots are at times \(t=0,0.01,\ldots ,0.14\). We also show a plot at time \(t=0.145\), together with a plot of the discrete energy over time. Below we visualize the axisymmetric surface \({\mathscr {S}}^m\) generated by \(\varGamma ^m\) at times \(t=0\) and \(t=0.145\)

Finally, we present an example for conserved mean curvature flow, (1.3) with (1.4c), for the scheme \(({\mathscr {A}}_m^{f,V})^h\) with the initial data \(\vec {X}^0\) parameterizing a closed spiral, so that the approximated surface has genus 1. As can be seen from Fig. 14, the spiral slowly untangles, until the surface becomes a torus. For this experiment we use the discretization parameters \(J=1024\) and \(\varDelta t= 10^{-6}\). The relative enclosed volume loss for this experiment is \(0.01\%\). The evolutions for the schemes \(({\mathscr {C}}_{m,\star }^{f,V})^{h}\) and \(({\mathscr {C}}_{m,\star }^{f,V})\) are nearly identical to what is shown in Fig. 14, with a relative volume loss of \(0.01\%\) in both cases.

Fig. 14
figure 14

\(({\mathscr {A}}_m^{f,V})^h\) for (1.4c). Conserved mean curvature flow. Plots are at times \(t=0,0.01,\ldots ,0.05,0.1\). We also show a plot of the discrete energy over time. Below we visualize the part of the axisymmetric surface \({\mathscr {S}}^m\) generated by \(\varGamma ^m \cap {{\mathbb {R}}}\times [-0.2,\infty )\) at times \(t=0\), \(t=0.01\), \(t=0.02\) and \(t=0.05\)

5.3 Numerical results for nonlinear mean curvature flow

Similarly to (5.2), it is easy to show that a sphere of radius r(t), with

$$\begin{aligned} r(t) = [1 - 2^\beta \,(\beta +1)\,t]^\frac{1}{\beta +1},\quad r(0) = 1, \end{aligned}$$
(5.6)

is a solution to (1.2) with (1.4a). We use this true solution for a convergence test for \(\beta = \tfrac{1}{2}\), similarly to Table 2 in [7]. Here we start with the nonuniform partitioning (5.3) of a semicircle of radius \(r(0)=r_0=1\) and compute the flow until time \(T = \tfrac{1}{2}\,{\overline{T}}\), where \({\overline{T}} = \tfrac{2}{3}\,2^{-\frac{1}{2}}\) denotes that extinction time of the shrinking sphere. We compute the error \(\Vert \varGamma - \varGamma ^h\Vert _{L^\infty }\), recall (5.4), over the time interval [0, T] between the true solution (5.6) and the discrete solutions for the schemes \(({\mathscr {A}}_m^f)^h\) and \(({\mathscr {C}}_{m,\star }^f)\). Here we used the time step size \(\varDelta t=0.1\,h^2_{\varGamma ^0}\), where \(h_{\varGamma ^0}\) is the maximal edge length of \(\varGamma ^0\). The computed errors are reported in Table 6.

Table 6 Errors for the convergence test for (5.6) over the time interval \([0,\tfrac{1}{2}\,\overline{T}]\)

We repeat the same convergence experiment for the inverse mean curvature flow, where we note that a sphere of radius r(t), with

$$\begin{aligned} r(t) = \exp (\tfrac{1}{2}\,t),\quad r(0) = 1, \end{aligned}$$
(5.7)

is a solution to (1.2) with (1.4b). The errors are reported in Table 7. We recall that these numbers can be compared to the fully 3d results in Table 3 in [7]. It is clear from Tables 6 and 7 that the solutions to the scheme \(({\mathscr {A}}_m^f)^h\) appear to converge with the optimal convergence rate of \({\mathscr {O}}(h^2_{\varGamma ^0})\). For the scheme \(({\mathscr {C}}_{m,\star }^f)\), on the other hand, the solutions appear to converge with an order less than quadratic, and closer to \(\frac{3}{2}\). We believe that these lower convergence rates are caused by the nonuniform meshes induced by the scheme \(({\mathscr {C}}_{m,\star }^f)\), recall Fig. 10, and also by the degeneracy of the coefficients \(\vec {x}\cdot \,\vec {e}_1\) in \(({\mathscr {C}}^f)\).

Table 7 Errors for the convergence test for (5.7) over the time interval [0, 1]

In the next experiment we repeat the simulation in [7, Fig. 8] for the inverse mean curvature of a torus with radii \(R=1\), \(r=0.25\). We recall that for this nonconvex initial data, with \(\varkappa _{{\mathscr {S}}}(\cdot ,0)<0\), the classical inverse mean curvature develops a singularity in finite time, see also [27, 43]. For the axisymmetric setting we use \(I = {{\mathbb {R}}}/{\mathbb {Z}}\), so that \(\partial I = \emptyset \). As the discretization parameters for the scheme \(({\mathscr {A}}_m^f)^h\) we use \(J=256\) and \(\varDelta t= 10^{-4}\). See Fig. 15 for the simulation results. Similarly to the results in [7, Fig. 8], the discrete solution becomes unphysical after around time 0.52, where we conjecture that the singularity for the continuous flow occurs.

Fig. 15
figure 15

\(({\mathscr {A}}_m^{f})^h\) for (1.4b). Inverse mean curvature flow for a torus with radii \(R=1\), \(r=0.25\). Plots are at times \(t=0,0.05,\ldots ,0.55\). We also visualize the axisymmetric surface \({\mathscr {S}}^m\) generated by \(\varGamma ^m\) at times \(t=0\) and \(t=0.5\)

5.4 Numerical results for Gauss curvature flow

An experiment for Gauss curvature flow, (1.7), for the same initial data as in Fig. 11, can be seen in Fig. 16. Here we have once again that \(\partial _0 I = \partial I = \{0,1\}\). The discretization parameters for the scheme \(({\mathscr {A}}^F_{m})^{h}\) from Remark 4.2 are \(J=128\) and \(\varDelta t= 10^{-5}\). As a comparison, we also show the evolution for standard mean curvature flow, computed with the scheme \(({\mathscr {A}}_m)^{h}\), in Fig. 16. It was suggested by Firey [25], that surfaces of stones, which are pounded by waves and other stones, move according to Gauss curvature flow. It is more likely that parts of the surface, where both principal curvature directions are highly curved, will be hit by waves and other stones. He hence proposed the Gauss curvature flow as the governing equation for the evolution of the stone’s surface. In Fig. 16 it is clearly seen that the upper and lower part, which have two highly curved principal curvature directions, move faster within Gauss curvature flow when compared to mean curvature flow. The parts closer to the origin have a nearly flat principal curvature direction. Hence they move far slower under Gauss curvature flow than under mean curvature flow, as can be clearly seen in Fig. 16.

Fig. 16
figure 16

\(({\mathscr {A}}_m^{F})^h\) for (1.7). Gauss curvature flow for a cigar, on the left. Plots are at times \(t=0,0.01,\ldots ,0.12\). We also visualize the axisymmetric surface \({\mathscr {S}}^m\) generated by \(\varGamma ^m\) at time \(t=0.12\). As a comparison, we show the evolution of \(({\mathscr {A}}_m)^h\) for mean curvature flow (1.1), on the right. Here the plots are at times \(t=0,0.01,\ldots ,0.06\), and we visualize the axisymmetric surface \({\mathscr {S}}^m\) generated by \(\varGamma ^m\) at time \(t=0.06\)

The Gauss curvature flow is not well-defined for general hypersurfaces as for non-convex hypersurfaces the resulting equation is not parabolic, see [3, 34]. In the axisymmetric situation the degrees of freedom are reduced and the resulting equation is

$$\begin{aligned} \vec {x}_t\cdot \,\vec {\nu } = \varkappa \,\frac{\vec {\nu }\cdot \,\vec {e}_1}{\vec {x}\cdot \,\vec {e}_1} \quad \text {on }\ I, \end{aligned}$$

which is parabolic as long as \({\vec {\nu }\cdot \,\vec {e}_1}\) is positive. In conclusion, even in the axisymmetric case the evolution is not well-defined if the initial surface has the topology of a torus, and so we do not present results for genus 1 surfaces. However, although this would not be possible in the general formulation we can start the Gauss curvature flow in the axisymmetric case with some nonconvex initial data, and we do so in the simulation in Fig. 17, where we used the discretization parameters \(J=128\) and \(\varDelta t= 10^{-5}\). For a mathematical analysis for Gauss curvature flow in the axisymmetric case we refer to [34].

Fig. 17
figure 17

\(({\mathscr {A}}_m^{F})^h\) for (1.7). Gauss curvature flow for nonconvex initial data. Plots are at times \(t=0,0.05,\ldots ,0.6\). We also visualize the axisymmetric surface \({\mathscr {S}}^m\) generated by \(\varGamma ^m\) at times \(t=0\), \(t=0.2\) and \(t=0.6\)

6 Conclusions

We have derived and analysed various numerical schemes for the parametric approximation of axisymmetric mean curvature flow, its nonlinear and volume conserving variants, as well as more general curvature flows. The main fully discrete schemes to consider for standard mean curvature flow are \(({\mathscr {A}}_m)^h\) and \(({\mathscr {C}}_{m,\star })\). Here we have dismissed the scheme \(({\mathscr {B}}_m)^h\), as it can show oscillations in practice, recall Fig. 3, as well as the scheme \(({\mathscr {C}}_{m,\star })^h\), as it can display very nonuniform meshes in simulations where the discrete curves are attached to the \(x_2\)-axis, recall Fig. 10 for its variant \(({\mathscr {C}}_{m,\star }^{f,V})^h\). We also do not consider the schemes \(({\mathscr {D}}_{m,\star })^{(h)}\), as they have no advantage over \(({\mathscr {C}}_{m,\star })\) and as they can also exhibit very nonuniform meshes. Of the two schemes we consider, the scheme \(({\mathscr {A}}_m)^h\) is a linear scheme that asymptotically leads to an equidistribution of mesh points, recall Remark 3.1. In addition, even though there is no stability proof for \(({\mathscr {A}}_m)^h\), in practice the discrete energy is always monotonically decreasing. The scheme \(({\mathscr {C}}_{m,\star })\), on the other hand, is a nonlinear scheme that is unconditionally stable. The nonlinearity is only very mild, and so a Newton solver never takes more than 3 iterations in practice. Moreover, the distribution of vertices for \(({\mathscr {C}}_{m,\star })\) may be worse than for \(({\mathscr {A}}_m)^h\), but coalescence of vertices is not observed in practice. Similar statements hold for the nonlinear variants \(({\mathscr {A}}_m^f)^h\), \(({\mathscr {A}}_m^{f,V})^h\), \(({\mathscr {C}}_{m,\star }^f)\) and \(({\mathscr {C}}_{m,\star }^{f,V})\), where the two conserving schemes show excellent volume conservation properties in practice, with the observed relative volume changes always being less than \(0.1\%\). Finally, for general curvature flows of the form (1.6), we propose the linear scheme \(({\mathscr {A}}_m^F)^h\), which asymptotically exhibits equidistributed mesh points.