1 Introduction↩︎

Surface evolution under geometric flows, including mean curvature flow, Willmore flow, and surface diffusion flow, has been intensively studied in the last decades [1]–[5] both from the viewpoint of applications and mathematical analysis. However, approximating surface evolution under geometric flows by numerical methods remains a challenging task. Numerical approximation of the mean curvature flow was first addressed by Dziuk [6] in 1990 through a finite element method (FEM) with moving nodes that determine the approximate surface. In the literature, such FEMs appear under various synonymous names: Evolving surface FEM, evolving FEM, and parametric FEM. Here we will often just speak of FEM for brevity. Many researchers have adopted and extended FEMs to various geometric flows of curves and surfaces, including:

• The development of evolving surface FEMs by Dziuk and Elliott [7].

• The development of parametric FEMs with re-meshing techniques by Bänsch, Morin, and Nochetto [8].

• The introduction of tangential motions to improve the mesh quality (i.e., mesh points distribution and shape of triangles) of computed surfaces by Barrett, Garcke, and Nürnberg [9], [10].

• The development of mesh quality improving tangential motions using harmonic map heat flow and the DeTurck trick [11], [12]. Semi-discrete error estimates for a DeTurck-based FEM discretisation were recently shown in [13].

• The development of tangential motions which make the evolving surfaces have minimal deformation rate (MDR) [14].

• The development of tangential motions ensuring minimal deformation from the initial surface [15].

• The development of tangential motions via a dual formulation of the minimal deformation rate (dual-MDR) scheme [16].

The first rigorous proof of convergence for evolving FEMs applied to mean curvature flow on closed surfaces was established in [17]. This approach involves reformulating the mean curvature flow into an equivalent system of solution-driven surface evolution using Huisken’s evolution equations for the normal vector and mean curvature [18], as well as analyzing the error of the evolving FEM using the matrix-vector framework and its properties established in [19], [20]. This approach was later extended to forced mean-curvature flow [21], [22] and further to Willmore flow and surface diffusion [23], reformulating Willmore flow and surface diffusion as a solution-driven surface evolution and proving the convergence of the FEM for this equivalent formulation. These advances of numerical analysis also allowed proving the convergence of Dziuk’s original semi-implicit FEM (with finite elements of degree \(k\ge 3\)) for mean curvature flow [24], as well as the convergence of a stabilized FEM incorporating artificial tangential motions of the Barrett–Garcke–Nürnberg type for curve shortening flow [25].

Despite these advances in the numerical analysis of mean curvature flow, existing algorithms with proven convergence so far lack at least one of the following desirable features: (i) energy stability and (ii) an explicit tangential motion mechanism for improving the quality of the evolving mesh. Both properties have been achieved by several alternative approaches, for example:

• Energy-stable FEMs augmented with artificial tangential motion of the Barrett–Garcke–Nürnberg type to improve mesh quality (e.g., node distribution and triangle shape) on computed surfaces; see [9], [10], [26].

• Energy-stable FEMs based on Lagrange-multiplier formulations [27], [28], including tangential motions of the Barrett–Garcke–Nürnberg type and of minimum-deformation-rate type.

• High-order energy-stable algorithm with artificial tangential motion based on a dual-MDR formulation [16] for both mean curvature flow and surface diffusion.

• Structure-preserving high-order continuous-in-time Petrov–Galerkin discretisation of geometric surface flows [29]. This approach uses auxiliary variables and is provably energy-stable and volume-conserving.

However, a rigorous convergence and error analysis for these methods remains open and challenging.

In this paper, we develop a rigorous error analysis for energy-stable FEMs based on a Lagrange multiplier approach, building upon the existing error analysis of evolving FEMs in [17]. This approach can be combined with incorporating tangential motions to improve the mesh quality, building additionally on [14]. We introduce new frameworks for the numerical methods and geometric flows and integrate them into the existing error analysis. In this way our results provide a unified theoretical foundation for understanding and advancing numerical methods in this field.

While we present the Langrange multiplier approach to preserving energy decay only for mean curvature flow, the approach could be similarly used to obtain a fully discrete energy-decaying method for Willmore flow and an energy-decaying and volume-preserving method for surface diffusion, allowing for a rigorous convergence analysis by combining the techniques of this paper with the error analysis in [23], [30], and presumably this can be extended to further geometric gradient flows.

The paper is organized as follows. In Section 2 we introduce continuous Lagrange-multiplier formulations of mean curvature flow, first without tangential motion and then with an artificial tangential velocity of minimal-deformation-rate type. Section 3 recalls the evolving surface finite element discretization of the solution-driven formulation of mean curvature flow used in [17]. In Section 4 we formulate the corresponding semi-discrete Lagrange-multiplier methods and show how the multiplier enforces the discrete area decay. Section 5 presents the fully discrete linearly-implicit BDF and implicit Adams schemes, including the scalar nonlinear equation for the Lagrange multiplier, and states the main error estimates. The proof for the method without tangential motion is given in Section 6, where the main additional ingredient is a stability estimate for the multiplier. Section 7 extends the argument to the method with minimal-deformation-rate tangential motion. Finally, Section 8 reports implementation details and numerical experiments illustrating the convergence rates, the area decay, and the computational overhead of the Lagrange-multiplier correction.

2 Continuous formulations with Lagrange multipliers↩︎

In this paper, we shall develop a Lagrange multiplier approach for constructing energy-stable FEMs for mean curvature flow and related problems. Consider the evolution of a surface under mean curvature flow with initial surface \(\Gamma^0\), and let the surface \(\Gamma[X(\cdot,t)]\) be the image of a flow map \(X(\cdot,t): \Gamma^0 \rightarrow\mathbb{R}^3\). The original continuous reformulation of mean curvature flow that is discretized by Kovács, Li, and Lubich in [17] incorporates the evolution equations for the normal vector \(n\) and mean curvature \(H\) derived by Huisken [18] and is written as follows: \[\begin{align} \label{KLL-form} \begin{aligned} \partial_t X &= v \circ X, \\ v &= -Hn, \\ \partial^{\tiny\bullet}n - \Delta_{\Gamma[X]} n &= |\nabla_{\Gamma[X]} n|^2 n, \\ \partial^{\tiny\bullet}H - \Delta_{\Gamma[X]} H &= |\nabla_{\Gamma[X]} n|^2 H, \end{aligned} \end{align}\tag{1}\] where \(\partial^{\tiny\bullet}\) is the material derivative, i.e. \(\partial^{\tiny\bullet}u(x,t)= \frac{\textrm{d}}{\textrm{d}t}u(X(q,t),t)\) at \(x=X(q,t)\) with \(q\in\Gamma^0\), and where \(v, n\), and \(H\) denote the velocity, normal vector, and mean curvature of the evolving surface \(\Gamma[X] = \Gamma[X(\cdot,t)]\). However, direct discretization of 1 by FEMs does not preserve the energy (area) decrease in the continuous mean curvature flow.

We introduce an additional Lagrange multiplier \(\lambda\) to this formulation, rewriting it into the following form: \[\tag{2} \begin{align} \partial_t X &= (1+\lambda) \, v \circ X, \\ v &= -Hn , \notag \\ \partial^{\tiny\bullet}n - \Delta_{\Gamma[X]} n &= |\nabla_{\Gamma[X]} n|^2 n, \\ \partial^{\tiny\bullet}H - \Delta_{\Gamma[X]} H &= |\nabla_{\Gamma[X]} n|^2 H, \\ \frac{\textrm{d}}{\textrm{d}t}\int_{\Gamma[X]} 1&= -\int_{\Gamma[X]} H^2 . \tag{3} \end{align}\] The last equation is a constraint corresponding to the Lagrange multiplier \(\lambda\), as the exact solution of mean curvature flows satisfies the following relation: \[\label{area-decay} \frac{\textrm{d}}{\textrm{d}t} \int_{\Gamma[X]} 1 = \int_{\Gamma[X]} \nabla_{\Gamma[X]} \cdot v = \int_{\Gamma[X]} Hn \cdot v = -\int_{\Gamma[X]} H^2 ,\tag{4}\] where the first equality is the Leibniz formula and the second equality is obtained from \(\nabla_{\varGamma}\cdot v = \nabla_{\varGamma}\cdot (- H n) = - \nabla_{\varGamma}H \cdot n - H \nabla_{\varGamma}\cdot n = - H^2\) ; cf. [18]. Therefore, the exact solution of mean curvature flow, together with \(\lambda=0\), satisfies 2 . In other words, the solution of the continuous formulation in 1 automatically satisfies 2 with \(\lambda = 0\).

The numerical solution of 1 may not satisfy the energy-decay property \(\frac{\textrm{d}}{\textrm{d}t} \int_{\Gamma[X]} 1 \leq 0,\) while the numerical solution of 2 still retains this energy-decay property due to the constraint in 3 . One challenge to be overcome in this paper is the convergence of FEM together with an appropriate time discretization for 2 , which is not obvious due to the introduction of the Lagrange multiplier and constraint.

Moreover, we shall further consider a Lagrange multiplier formulation that incorporates an artificial tangential motion of MDR-type: \[\tag{5} \begin{align} \partial_t X &= (1+\lambda)\, v \circ X, \\ v \cdot n &= -H , \tag{6}\\ -\Delta_{\Gamma[X]} v &= -\mu n, \tag{7}\\ \partial^{\tiny\bullet}n - v\cdot\nabla_{\Gamma[X]}n - \Delta_{\Gamma[X]} n &= |\nabla_{\Gamma[X]} n|^2 n, \\ \partial^{\tiny\bullet}H - v\cdot\nabla_{\Gamma[X]}H - \Delta_{\Gamma[X]} H &= |\nabla_{\Gamma[X]} n|^2 H, \\ \frac{\textrm{d}}{\textrm{d}t}\int_{\Gamma[X]} 1&= -\int_{\Gamma[X]} H^2 . \tag{8} \end{align}\] The second and third equations generate a tangential motion that minimizes the deformation rate \[\int_{\Gamma[X]} |\nabla_{\Gamma[X]} v|^2\] of the evolving surface under the constraint 6 . The unknown function \(\mu\in H^{-1}(\Gamma[X])\) is a Lagrange multiplier of this constrained optimization problem, which is to be determined. Discretization of 5 maintains mesh quality as the surface evolves; see the motivation and discussions in [14]. Here we can naturally incorporate this MDR tangential motion into our Lagrange multiplier approach. Again, the convergence of FEMs together with an appropriate time discretization for this continuous formulation is a challenge to be addressed in this paper.

3 Evolving FEM for mean curvature flow: method of [17]↩︎

In this section, we recall the evolving FEM for 1 discussed in [17], using the notion of evolving surface FEM described in [7], [31]. In this method, the given smooth initial surface \(\Gamma^0\) is approximated by a piecewise curved triangulated surface that interpolates \(\Gamma^0\), with each triangular piece of \(\Gamma_h^0\) being the image of the reference plane triangle under a polynomial map of degree \(k\); see [7], [32], [33]. The nodal vector \({\mathbf{x}}^0=(x_1^0, \cdots, x_N^0) \in \mathbb{R}^{3N}\), which collects all the nodes \(x_j^0\) \((j=1,\dots,N)\) of \(\Gamma_h^0\), is evolved in time, and its value at time \(t\) is denoted by \({\mathbf{x}}(t)\), determining a triangulated surface \(\Gamma_h[{\mathbf{x}}(t)]\) that approximates \(\Gamma(t)\).

The finite element spaces on \(\varGamma_h[{\mathbf{x}}(t)]\) and \(\varGamma_h[{\mathbf{x}}^*(t)]\) are denoted by \(S_h[{\mathbf{x}}(t)]\) and \(S_h[{\mathbf{x}}^*(t)]\), respectively. The evolving FEM for mean curvature flow discussed in [17] essentially seeks nodal vectors \({\mathbf{x}}(t) \in \mathbb{R}^{3N}\), \({\mathbf{v}}(t)\in \mathbb{R}^{3N}\) and \({\mathbf{u}}(t)=({\mathbf{n}}(t);{\mathbf{H}}(t))\in \mathbb{R}^{4N}\) which approximate the interpolated nodal vectors \({\mathbf{x}}^*(t)\), \({\mathbf{v}}^*(t)\) and \({\mathbf{u}}^*(t)\), respectively.

The semi-discrete evolving finite element scheme is to find nodal vectors \({\mathbf{x}}(t) \in \mathbb{R}^{3N}\) determining the surface \(\varGamma_h[{\mathbf{x}}(t)]\) and their velocities \({\mathbf{v}}(t)\in \mathbb{R}^{3N}\) and the nodal vector \({\mathbf{u}}(t)=({\mathbf{n}}(t);{\mathbf{H}}(t))\in \mathbb{R}^{4N}\) of a finite element function \(u_h(\cdot,t)=(n_h(\cdot,t);H_h(\cdot,t))\in S_h[{\mathbf{x}}(t)]^4\) satisfying the following weak formulation for all finite element test functions \({\varphi_h}\in S_h[{\mathbf{x}}(t)]^4\) and all \(t\in[0,T]\): omitting the argument \(t\), \[\label{eq:MCF32discrete} \begin{align} \int_{\varGamma_h[{\mathbf{x}}]} \!\!\! \partial^{\tiny\bullet}_h u_h\cdot \varphi_h + \int_{\varGamma_h[{\mathbf{x}}]} \!\!\!\! \nabla_{\varGamma_h[{\mathbf{x}}]} u_h \cdot \nabla_{\varGamma_h[{\mathbf{x}}]} \varphi_h = &\;\int_{\varGamma_h[{\mathbf{x}}]} \!\!\! |\nabla_{\Gamma_h[{\mathbf{x}}]} n_h|^2 u_h \cdot \varphi_h \\[2mm] {\mathbf{v}}= &\;- {\mathbf{H}}\bullet {\mathbf{n}}\\ {\boldsymbol{\dot{\mathbf{x}}}} = &\;{\mathbf{v}}, \end{align}\tag{9}\] where \(\partial^{\tiny\bullet}_h u_h\) is the material derivative of \(u_h\) on \(\varGamma_h[{\mathbf{x}}]\), which equals the finite element function with nodal vector \({\boldsymbol{\dot{\mathbf{u}}}}=d{\mathbf{u}}/dt\). Furthermore, \(\bullet\) denotes the componentwise product of vectors, i.e. \(({\mathbf{H}}\bullet {\mathbf{n}})_i=H_i n_i = H_h(x_i) n_h(x_i)\). We mention that in [17] the discrete velocity was determined via a Ritz projection, but it was found later in [23] that the simpler construction of \({\mathbf{v}}\) as given above also yields a stable discretization of the same order of convergence.

For \(t \in [0,T]\), the globally continuous finite element nodal basis functions which span \(S_h[{\mathbf{x}}(t)]\) are denoted by \(\phi_j[{\mathbf{x}}(t)]\), \(j=1,\dotsc,N\). They are the unique functions whose piecewise pullback to the flat reference triangle is a polynomial of degree \(k\) and which have the nodal values \(\phi_j[{\mathbf{x}}(t)](x_i(t)) = \delta_{ij}\) for all \(i,j = 1, \dotsc, N\).

Matrix–vector formulation. The surface mass matrix \({\mathbf{M}}({\mathbf{x}})\in \mathbb{R}^{N\times N}\) and the stiffness matrix \({\mathbf{A}}({\mathbf{x}})\in \mathbb{R}^{N\times N}\)are defined in terms of these finite element basis functions via their \((i,j)\)th entries, for \(i, j = 1, \ldots, N\), \[\begin{align} \label{def-matrix-M-A} \begin{aligned} {\mathbf{M}}({\mathbf{x}})|_{ij} = &\;\int_{\varGamma_h[{\mathbf{x}}]} \!\! \phi_i[{\mathbf{x}}] \, \phi_j[{\mathbf{x}}] ,\\ {\mathbf{A}}({\mathbf{x}})|_{ij} = &\;\int_{\varGamma_h[{\mathbf{x}}]} \!\!\!\! \nabla_{\varGamma_h[{\mathbf{x}}]} \phi_i[{\mathbf{x}}] \cdot \nabla_{\varGamma_h[{\mathbf{x}}]} \phi_j[{\mathbf{x}}] . \end{aligned} \end{align}\tag{10}\] We introduce the vector \({\mathbf{f}}({\mathbf{x}},{\mathbf{u}})\in \mathbb{R}^{4N}\) with entries \[\begin{align} {\mathbf{f}}({\mathbf{x}},{\mathbf{u}})|_{j+(\ell-1)N} = &\int_{\varGamma_h[{\mathbf{x}}]} \!\!\! |\nabla_{\varGamma_h[{\mathbf{x}}]}n_h|^2 u_h \cdot \phi_j[{\mathbf{x}}] , \end{align}\] for \(j = 1,\dotsc,N\) and \(\ell = 1,2,3,4\).

The formulation of the semi-discretized coupled system for mean curvature flow 9 then leads to the matrix–vector formulation of [17]: writing \({\mathbf{M}}({\mathbf{x}})\) and \({\mathbf{A}}({\mathbf{x}})\) in short for \({\mathbf{M}}({\mathbf{x}})\otimes I_4\) and \({\mathbf{A}}({\mathbf{x}})\otimes I_4\) with the 4-dimensional identity matrix \(I_4\), we have with \({\mathbf{u}}=({\mathbf{n}};{\mathbf{H}})\) \[\label{eq:MCF32matrix-vector} \begin{align} {\mathbf{M}}({\mathbf{x}}) {\boldsymbol{\dot{\mathbf{u}}}} + {\mathbf{A}}({\mathbf{x}}) {\mathbf{u}}= &\; {\mathbf{f}}({\mathbf{x}},{\mathbf{u}}) , \\ {\mathbf{v}}= &\;- {\mathbf{H}}\bullet {\mathbf{n}}, \\ {\boldsymbol{\dot{\mathbf{x}}}} = &\;{\mathbf{v}}. \end{align}\tag{11}\]

4 Semi-discrete methods with area decay↩︎

4.1 Area-decreasing method without tangential velocity↩︎

The area-decreasing evolving FEM for the coupled system 3 reads as follows. Find nodal vectors \({\mathbf{x}}(t) \in \mathbb{R}^{3N}\) determining the surface \(\varGamma_h[{\mathbf{x}}(t)]\) and their velocities \({\mathbf{v}}(t)\in \mathbb{R}^{3N}\) and the nodal vector \({\mathbf{u}}(t)=({\mathbf{n}}(t);{\mathbf{H}}(t))\in \mathbb{R}^{4N}\) of a finite element function \(u_h(t)=u_h(\cdot,t)=(n_h(t);H_h(t))\in S_h[{\mathbf{x}}(t)]^4\) satisfying the following weak formulation: for all \(t\in[0,T]\) and for all finite element test functions \({\varphi_h}\in S_h[{\mathbf{x}}(t)]^4\) we require (omitting the argument \(t\)) \[\label{eq:MCF32discrete32area-decreasing} \begin{align} \int_{\varGamma_h[{\mathbf{x}}]} \!\!\! \partial^{\tiny\bullet}_h u_h\cdot \varphi_h + \int_{\varGamma_h[{\mathbf{x}}]} \!\!\!\! \nabla_{\varGamma_h[{\mathbf{x}}]} u_h \cdot \nabla_{\varGamma_h[{\mathbf{x}}]} \varphi_h = &\;\int_{\varGamma_h[{\mathbf{x}}]} \!\!\! |\nabla_{\Gamma_h[{\mathbf{x}}]} n_h|^2 u_h \cdot \varphi_h \\[2mm] {\mathbf{v}}= &\;- {\mathbf{H}}\bullet {\mathbf{n}}\\ {\boldsymbol{\dot{\mathbf{x}}}} = &\;\big( 1 + \lambda_h \big){\mathbf{v}}, \end{align}\tag{12}\] where \(\lambda_h(t) \in \mathbb{R}\) is a Lagrange multiplier which is to ensure that the area (energy) of the computed surface, i.e., \[\mathcal{E}_h[{\mathbf{x}}(t)] := \mathcal{E}(\varGamma_h[{\mathbf{x}}(t)]) = \int_{\varGamma_h[{\mathbf{x}}(t)]} 1 ,\] decreases along the semi-discrete flow analogously to the continuous flow, see 3 : \[\label{eq:discrete32MCF32energy32decay32-32a} \frac{\textrm{d}}{\textrm{d}t}\mathcal{E}_h[{\mathbf{x}}] = - \| H_h \|_{L^2(\varGamma_h[{\mathbf{x}}])}^2 \leq 0.\tag{13}\] With the finite element functions \(x_h\) and \(v_h\) on \(\Gamma_h[{\mathbf{x}}]\) with nodal vectors \({\mathbf{x}}\) and \({\mathbf{v}}\), respectively, we obtain \[\label{eq:discrete32MCF32energy32decay32-32b} \begin{align} \frac{\textrm{d}}{\textrm{d}t}\int_{\varGamma_h[{\mathbf{x}}]} 1 & = \int_{\varGamma_h[{\mathbf{x}}]} \nabla_{\varGamma_h[{\mathbf{x}}]} \cdot {\partial^{\tiny\bullet}_h x_h} = (1+\lambda_h) \int_{\varGamma_h[{\mathbf{x}}]}\nabla_{\varGamma_h[{\mathbf{x}}]} \cdot v_h , \end{align}\tag{14}\] because of the Leibniz formula in the first equality and because the third relation in 12 and the transport property \(\partial^{\tiny\bullet}_h \phi_i[{\mathbf{x}}]=0\) of the basis functions imply \({\partial^{\tiny\bullet}_h x_h} = (1+\lambda_h) v_h,\) which yields the second equality. Therefore, the combination of 13 and 14 yields the following explicit formula for the semi-discrete Lagrange multiplier: \[\label{eq:def32lambda95h} \lambda_h(t) = -1 - \frac{\| H_h \|_{L^2(\varGamma_h[{\mathbf{x}}])}^2}{\int_{\varGamma_h[{\mathbf{x}}]}\nabla_{\varGamma_h[{\mathbf{x}}]} \cdot v_h} ,\tag{15}\] provided that the denominator is nonzero. Here, \(\| H_h \|_{L^2(\varGamma_h[{\mathbf{x}}])}^2 = {\mathbf{H}}^{\rm T} {\mathbf{M}}({\mathbf{x}}) {\mathbf{H}}\) with the mass matrix \({\mathbf{M}}({\mathbf{x}})\) of 10 .

The matrix–vector formulation of the area-decreasing semi-discretization 12 reads: with \({\mathbf{u}}=({\mathbf{n}};{\mathbf{H}})\), \[\tag{16} \begin{align} {\mathbf{M}}({\mathbf{x}}) {\boldsymbol{\dot{\mathbf{u}}}} + {\mathbf{A}}({\mathbf{x}}) {\mathbf{u}}= &\; {\mathbf{f}}({\mathbf{x}},{\mathbf{u}}) , \\ {\mathbf{v}}= &\;- {\mathbf{H}}\bullet {\mathbf{n}}, \\ {\boldsymbol{\dot{\mathbf{x}}}} = &\;(1 + \lambda_h) {\mathbf{v}}, \\ \tag{17} \frac{\textrm{d}}{\textrm{d}t} \mathcal{E}_h[{\mathbf{x}}] = &\;- {\mathbf{H}}^{\rm T} {\mathbf{M}}({\mathbf{x}}) {\mathbf{H}}. \end{align}\] Note that 17 is equivalent to the explicit expression of \(\lambda_h\) in 15 . It is formulated to keep the semi-discrete scheme analogous to the fully discrete scheme that will be presented in the next section.

4.2 Area-decreasing method with artificial tangential motion↩︎

The area-decreasing evolving FEM with the MDR tangential motion introduced in [14] for the coupled system 5 reads as follows. Find a nodal vector \({\mathbf{x}}(t) \in \mathbb{R}^{3N}\) which determines the approximate surface \(\Gamma_h[{\mathbf{x}}(t)]\), as well as nodal vectors \({\mathbf{v}}(t)\in \mathbb{R}^{3N}\), \({\mathbf{u}}(t)=({\mathbf{n}}(t);{\mathbf{H}}(t))\in \mathbb{R}^{4N}\) and \(\boldsymbol{\mu}(t)\in\mathbb{R}^{N}\) such that the corresponding finite element functions \(v_h(t)=v_h(\cdot,t)\in S_h[{\mathbf{x}}(t)]^3\), \(u_h(t)=(n_h(t);H_h(t))\in S_h[{\mathbf{x}}(t)]^4\) and \(\mu_h(t)\in S_h[{\mathbf{x}}(t)]\) satisfy the following weak formulation: for all \(t\in [0,T]\) and for all \(\varphi_h\in S_h[{\mathbf{x}}(t)]^4\) and \(\psi_h\in S_h[{\mathbf{x}}(t)]^3\), \(\chi_h\in S_h[{\mathbf{x}}(t)]\), \[\tag{18} \begin{align} \int_{\varGamma_h[{\mathbf{x}}]} {\partial^{\tiny\bullet}_h u_h}\cdot {\varphi_h} + \int_{\varGamma_h[{\mathbf{x}}]} \!\!\!\! \nabla_{\Gamma_h[{\mathbf{x}}]} u_h \cdot \nabla_{\Gamma_h[{\mathbf{x}}]} {\varphi_h} = &\;\int_{\varGamma_h[{\mathbf{x}}]} \!\!\! |\nabla_{\Gamma_h[{\mathbf{x}}]} n_h|^2 u_h \cdot {\varphi_h} \notag\\ &\; + \int_{\varGamma_h[{\mathbf{x}}]} \!\!\! [(v_h\cdot \nabla_{\Gamma_h[{\mathbf{x}}]}) u_h] \cdot {\varphi_h} , \tag{19}\\[2mm] \int_{\varGamma_h[{\mathbf{x}}]} \nabla_{\Gamma_h[{\mathbf{x}}]} v_h\cdot\nabla_{\Gamma_h[{\mathbf{x}}]}\psi_h = &\;-\int_{\varGamma_h[{\mathbf{x}}]} \mu_h n_h\cdot \psi_h , \\ \int_{\varGamma_h[{\mathbf{x}}]} v_h\cdot n_h\, \chi_h = &\;-\int_{\varGamma_h[{\mathbf{x}}]} H_h {\chi_h} , \tag{20}\\[5mm] {\boldsymbol{\dot{\mathbf{x}}}} = &\;\big( 1 + \lambda_h \big) {\mathbf{v}}, \\ \frac{\textrm{d}}{\textrm{d}t}\mathcal{E}_h[{\mathbf{x}}(t)] = &\;- \| H_h \|_{L^2(\varGamma_h[{\mathbf{x}}])}^2 , \end{align}\] where again \(\lambda_h(t) \in \mathbb{R}\) is a Lagrange multiplier which ensures that the area of the computed surface decreases as time grows. According to the discussions in 13 –15 , \(\lambda_h\) is given by the explicit formula 15 .

The matrix–vector formulation of this area-decreasing evolving FEM with MDR tangential motion reads: with \({\mathbf{u}}=({\mathbf{n}};{\mathbf{H}})\), \[\label{eq:MCF-matrix-vector-MDR} \begin{align} {\mathbf{M}}({\mathbf{x}}) {\boldsymbol{\dot{\mathbf{u}}}} + {\mathbf{A}}({\mathbf{x}}) {\mathbf{u}}= &\; {\mathbf{f}}({\mathbf{x}},{\mathbf{u}},{\mathbf{v}}) , \\[5pt] \begin{pmatrix} {\boldsymbol{A}}({\mathbf{x}}) & {\boldsymbol{B}}({\mathbf{x}},{\mathbf{n}})^{\rm T} \\ {\boldsymbol{B}}({\mathbf{x}},{\mathbf{n}}) & {\mathbf{0}} \end{pmatrix} \begin{pmatrix} {\mathbf{v}}\\ \boldsymbol{\mu} \end{pmatrix} =&\; \begin{pmatrix} {\mathbf{0}}\\ - {\boldsymbol{M}}({\mathbf{x}}){\mathbf{H}} \end{pmatrix}, \\[5pt] {\boldsymbol{\dot{\mathbf{x}}}} = &\;(1 + \lambda_h) {\mathbf{v}}, \\[-2mm] \frac{\textrm{d}}{\textrm{d}t} \mathcal{E}_h[{\mathbf{x}}] = &\;- {\mathbf{H}}^{\rm T} {\mathbf{M}}({\mathbf{x}}) {\mathbf{H}}, \end{align}\tag{21}\] where \({\boldsymbol{f}}({\boldsymbol{x}},{\boldsymbol{u}},{\boldsymbol{v}})\) is the unique nodal vector in \(\mathbb{R}^{4N}\) satisfying \[{\boldsymbol{f}}({\boldsymbol{x}},{\boldsymbol{u}},{\boldsymbol{v}})^{\rm T} \boldsymbol{\varphi} = \int_{\varGamma_h[{\mathbf{x}}]} \!\!\! |\nabla_{\Gamma_h[{\mathbf{x}}]} n_h|^2 u_h \cdot {\varphi_h} + \int_{\varGamma_h[{\mathbf{x}}]} \!\!\! [(v_h\cdot \nabla_{\Gamma_h[{\mathbf{x}}]}) u_h] \cdot {\varphi_h}\] for all finite element functions \(\varphi_h\in S_h[{\boldsymbol{x}}]^4\) with nodal vector \(\boldsymbol{\varphi}\), and \({\boldsymbol{B}}({\mathbf{x}},{\mathbf{n}})\) is the unique \(N\times3N\) matrix satisfying the following condition: \[\boldsymbol{\chi}^{\rm T} {\boldsymbol{B}}({\mathbf{x}},{\mathbf{n}})\boldsymbol{\psi} = \int_{\varGamma_h[{\mathbf{x}}]} \psi_h\cdot n_h \chi_h \quadfor all\,\, \boldsymbol{\psi} \in\mathbb{R}^{3N}\,\,and\,\,\, \boldsymbol{\chi} \in \mathbb{R}^{N},\] for all finite element functions \(\psi_h\in S_h[{\mathbf{x}}]^3\) and \(\chi_h\in S_h[{\mathbf{x}}]\) with nodal vectors \(\boldsymbol{\psi}\) and \(\boldsymbol{\chi}\), respectively.

In the next section, we introduce the corresponding fully discrete FEMs with Lagrange multipliers. Error bounds for the semi-discrete FEMs in 16 and 21 (without and with tangential motion, respectively) follow directly from the fully discrete error bounds proved in Theorems 1 and 2 by taking the limit \(\tau \rightarrow 0\). Therefore, we do not state the semi-discrete error bounds separately.

5 Fully discrete methods with area decay↩︎

For the time discretization of the system of ordinary differential equations 16 we use a multistep approach: a \(q\)-step linearly implicit backward difference formula (BDF method), with \(2 \le q \leq 5\), for the (stiff) geometric equation, and a \(q\)-step implicit Adams method for the (nonstiff) differential equation for the positions.

5.1 Area-decreasing method without tangential velocity↩︎

The \(q\)-step backward difference formula (BDF) method is determined by its coefficients \(\delta_j\), given by \(\delta(\zeta)=\sum_{j=0}^q \delta_j \zeta^j=\sum_{\ell=1}^q \frac{1}{\ell}(1-\zeta)^\ell\). In its linearly implicit version, the BDF method uses extrapolation coefficients \(\gamma_j\) defined by \(\gamma(\zeta) = \sum_{j=0}^{q-1} \gamma_j \zeta^j = (1 - (1-\zeta)^q)/\zeta\). The BDF method is known to be zero-stable and \(A(\alpha)\)-stable with \(\alpha>0\) for \(q\leq 6\), see [34], and for \(q\le 5\), to converge with order \(q\) for linear parabolic equations on moving surfaces [35]. This order is retained by the linearly implicit variant for quasilinear parabolic equations using the above coefficients \(\gamma_j\); cf. [36], [37].

For a fixed step size \(\tau>0\), we define \(t_n = n \tau\) for \(n=0,1,2,\dots\), and denote by \(\widetilde{\mathbf{x}}^n\) the \(q\)-step extrapolated values defined by the following formula: \[\label{eq:extrapolation32def} \widetilde{\mathbf{x}}^n := \left\{\begin{align} &\sum_{j=0}^{q-1} \gamma_j {\mathbf{x}}^{n-1-j} &&for\,\,\, n \geq q , \\ &{\mathbf{x}}^n &&for\,\,\, 0\le n \le q-1 . \end{align}\right.\tag{22}\] We use the analogous definition for \(\widetilde{\mathbf{u}}^n=(\widetilde{\mathbf{n}}^n;\widetilde{\mathbf{H}}^n)\). Moreover, we denote by \({\boldsymbol{\dot{\mathbf{u}}}}^n\) the \(q\)-step backward finite difference defined by \[\begin{align} \label{eq:BDF32def} {\boldsymbol{\dot{\mathbf{u}}}}^n := \frac{1}{\tau} \sum_{j=0}^q \delta_j {\mathbf{u}}^{n-j} . \end{align}\tag{23}\] We further define the weighted average of the velocity values that appears in the \(q\)-step implicit Adams method, \[\begin{align} \label{def-v-beta} {\bar{\mathbf{v}}}^n = \sum_{j=0}^{q}\beta_j\,\mathbf{v}^{\,n-j}, \end{align}\tag{24}\] with the weights \(\beta_j=\frac{1}{\tau}\int_{t_{n-1}}^{t_n} L_j^{n}(t)\,\mathrm{d}t,\) where \(L_j^{n}(t)\) are the Lagrange basis polynomials of degree \(q\) associated with the nodes \(\{t_n,t_{n-1},\ldots,t_{n-q}\}\), i.e., \(L_j^{n}(t_{n-i})=\delta_{ij}\) for \(i,j=0,\ldots,q\). The coefficients \((\beta_j)_{j=0}^{q}\) define the \(q\)-step implicit Adams method and can be precomputed and stored for use in the implementation; see Table 1. It is well known that the \(q\)-step implicit Adams method has convergence order \(q+1\) when applied to nonstiff ordinary differential equations; see [38].

For \(n\ge q\), we determine the approximations \({\mathbf{x}}^n\) to \({\mathbf{x}}^*(t_n)\), \({\mathbf{v}}^n\) to \({\mathbf{v}}^*(t_n)\), \({\mathbf{u}}^n\) to \({\mathbf{u}}^*(t_n)\) by the fully discrete linearly implicit method \[\tag{25} \begin{align} \tag{26} {\mathbf{M}}(\widetilde{{\mathbf{x}}}^n) {\boldsymbol{\dot{\mathbf{u}}}}^n + {\mathbf{A}}(\widetilde{{\mathbf{x}}}^n) {\mathbf{u}}^n = &\; {\mathbf{f}}(\widetilde{{\mathbf{x}}}^n,\widetilde{{\mathbf{u}}}^n) , \\ \tag{27} {\mathbf{v}}^n = &\;- {\mathbf{H}}^n \bullet {{\mathbf{n}}}^n , \\ \tag{28} \frac{1}{\tau} \big( {\mathbf{x}}^n - {\mathbf{x}}^{n-1} \big) = &\;(1 + \lambda^n) {\bar{\mathbf{v}}}^n, \\ \tag{29} \mathcal{E}_h[{\mathbf{x}}^n] - \mathcal{E}_h[{\mathbf{x}}^{n-1}] = &\;- \int_{t_{n-1}}^{t_n} p^n(t)^2 \,\textrm{d}t , \end{align}\] where \(\lambda^n \in \mathbb{R}\) is a Lagrange multiplier which enforces 29 , \(\bar{\mathbf{v}}^n\) is defined by 24 , and \(p^n(t)\) is the interpolation polynomial through the points \((t_{n-j},\|H_h^{n-j}\|_{L^2(\varGamma_h[\widetilde{{\mathbf{x}}}^{n-j}])})\) for \(0\le j \le q\), where \(\|H_h^{n-j}\|_{L^2(\varGamma_h[\widetilde{{\mathbf{x}}}^{n-j}])} = \linebreak \bigl(({\mathbf{H}}^{n-j})^T{\mathbf{M}}(\widetilde{{\mathbf{x}}}^{n-j}){\mathbf{H}}^{n-j}\bigr)^{1/2}\). The surface area thus decreases with a rate consistent with the approximate mean curvature given by the numerical scheme; cf. 4 .

The starting values \({\mathbf{x}}^i\) and \({\mathbf{u}}^i\) (\(i=0,\dotsc,q-1\)) are assumed to be given. They can be precomputed using either a lower order method with smaller step sizes or an implicit Runge–Kutta method.

Due to the linearly implicit time discretisation, each time step requires the solution of a decoupled system, solved sequentially: first solving a linear system for the geometry \({\mathbf{u}}^n=({\mathbf{n}}^n;{\mathbf{H}}^n)\), then determining the velocity \({\mathbf{v}}^n\), then a simplified Newton iteration for \(\boldsymbol{\lambda}^n\), which finally determines \({\mathbf{x}}^n\). For practical details, we refer to Section 8.1.

The following result shows the existence of a Lagrange multiplier \(\lambda^n\) and its uniqueness in a small interval around \(0\), together with the convergence of the simplified Newton method under conditions that are substantially weaker than the strong regularity conditions required for the optimal-order error bounds in Theorem 1 below. We need some notation. Let \(\bar{\mathbf{x}}^n = {\mathbf{x}}^{n-1} + \tau \bar {\mathbf{v}}^{n}\) and denote the defect in 29 as \[\label{delta-n} \delta^n := \left| \frac{1}{\tau} \bigl(\mathcal{E}_h[\bar{\mathbf{x}}^n] - \mathcal{E}_h[{\mathbf{x}}^{n-1}]\bigr) + \frac{1}{\tau} \int_{t_{n-1}}^{t_n} p^n(t)^2 \,\textrm{d}t \right|.\tag{30}\] We assume that \[\label{sigma-n} \sigma^n := \Big| \int_{\Gamma_h[\bar {\mathbf{x}}^n]} \nabla_{\Gamma_h[\bar {\mathbf{x}}^n]}\cdot \bar v_h^n \Big| >0,\tag{31}\] where \(\bar v_h^n \in S_h[\bar {\mathbf{x}}^n]^3\) is the finite element function on \(\Gamma_h[\bar {\mathbf{x}}^n]\) with the nodal vector \(\bar {\mathbf{v}}^n\) defined in 24 . In view of 4 , this assumption is fulfilled as long as the numerical positions and velocities remain \(H^1\)-close to the exact positions and velocities.

With an absolute constant \(c_2\) (with \(c_2\le 4\)), we further define \[\label{eta-n} \eta^n := 3\biggl(\frac{c_2\, \| \nabla_{\Gamma_h[\bar{\mathbf{x}}^n]}\bar v_h^n \|_{L^2(\Gamma_h[\bar {\mathbf{x}}^n])^{3\times 3}}}{\sigma^n} \biggr)^2.\tag{32}\]

The following result states that the above area-decreasing fully discrete method satisfies optimal-order error bounds in the \(H^1\)-norm of the same type as the original full discretization in [17], which is not guaranteed to decrease the area under weak regularity assumptions.

5.2 Area-decreasing fully discrete FEM with tangential velocity↩︎

Analogous to 25 , while using the matrix-vector formulation in 21 , the fully discrete method with tangential velocity can be described as follows. For \(n\ge q\), we determine the approximations \({\mathbf{x}}^n\) to \({\mathbf{x}}^*(t_n)\), \({\mathbf{v}}^n\) to \({\mathbf{v}}^*(t_n)\), \(\boldsymbol{\mu}^n\) to \(\boldsymbol{\mu}^*(t_n)\), \({\mathbf{u}}^n\) to \({\mathbf{u}}^*(t_n)\) by the fully discrete system of linear equations: \[\tag{35} \begin{align} \tag{36} {\mathbf{M}}(\widetilde{{\mathbf{x}}}^n) {\boldsymbol{\dot{\mathbf{u}}}}^n + {\mathbf{A}}(\widetilde{{\mathbf{x}}}^n) {\mathbf{u}}^n = &\; {\mathbf{f}}(\widetilde{{\mathbf{x}}}^n,\widetilde{{\mathbf{u}}}^n,\widetilde{{\mathbf{v}}}^n) , \\[5pt] \tag{37} \begin{pmatrix} {\boldsymbol{A}}(\widetilde{\mathbf{x}}^n) & {\boldsymbol{B}}(\widetilde{\mathbf{x}}^n,{\mathbf{n}}^n)^{\rm T} \\ {\boldsymbol{B}}(\widetilde{\mathbf{x}}^n,{\mathbf{n}}^n) & {\mathbf{0}} \end{pmatrix} \begin{pmatrix} {\mathbf{v}}^n \\ \boldsymbol{\mu}^n \end{pmatrix} =&\; \begin{pmatrix} {\mathbf{0}}\\ - {\boldsymbol{M}}(\widetilde{\mathbf{x}}^n){\mathbf{H}}^n \end{pmatrix}, \\[5pt] \tag{38} \frac{1}{\tau} \big( {\mathbf{x}}^n - {\mathbf{x}}^{n-1} \big) = &\;(1 + \lambda^n) {\bar{\mathbf{v}}}^n , \\ \tag{39} \mathcal{E}_h[{\mathbf{x}}^n] - \mathcal{E}_h[{\mathbf{x}}^{n-1}] = &\;- \int_{t_{n-1}}^{t_n} p^n(t)^2 \, \textrm{d}t , \end{align}\] where \(\widetilde{\boldsymbol{v}}^n=\sum_{j=0}^{q-1}\gamma_j {\boldsymbol{v}}^{n-1-j}\) and \({\bar{\mathbf{v}}}^n:=\sum_{j=0}^{q} \beta_j {\mathbf{v}}^{n-j}\) and \(p^n(t)\) is again the interpolation polynomial through the points \((t_{n-j},\|{\mathbf{H}}^{n-j}\|_{{\mathbf{M}}(\widetilde{{\mathbf{x}}}^{n-j})})\) for \(0\le j \le q-1\). The solvability of the saddle point system 37 can be shown through error analysis for sufficiently small \(h\), as shown in [14].

6 Proof of Theorem 1↩︎

The proof is a modification of the proof of the fully discrete error bounds in [17]. The only new aspect is the effect on the stability of the discrete Lagrange multipliers \(\lambda^n\) of the numerical method, which is shown to be mild. We focus on this new aspect and recall results and techniques from the stability and consistency analysis of [17] as far as needed.

6.1 Preparations↩︎

6.1.1 Notation of norms of finite element functions↩︎

For the simplicity of notation, we identify a nodal vector \({\boldsymbol{w}}\) with the corresponding finite element function on a specified triangulated surface, such as \(\Gamma_h[{\mathbf{x}}]\) or \(\Gamma_h[{\mathbf{x}}^*]\), whenever we consider the norm of the finite element function on the triangulated surface. For example, for \(s\in\{0,1\}\) and \(1\le p\le\infty\), we denote by \(\| {\boldsymbol{w}} \|_{H^s(\Gamma_h[{\mathbf{x}}])}\) and \(\| {\boldsymbol{w}} \|_{W^{s,p}(\Gamma_h[{\mathbf{x}}])}\) the Sobolev norms of the finite element function on \(\Gamma_h[{\mathbf{x}}]\) having nodal vector \({\boldsymbol{w}}\). In particular, \[\label{norm-w} \begin{align} \| {\boldsymbol{w}} \|_{L^2(\Gamma_h[{\mathbf{x}}])} &= \sqrt{{\boldsymbol{w}}^{\rm T}{\mathbf{M}}({\mathbf{x}}){\boldsymbol{w}}} ,\\ \| \nabla_{\Gamma_h[{\mathbf{x}}]}{\boldsymbol{w}} \|_{L^2(\Gamma_h[{\mathbf{x}}])} &= \sqrt{{\boldsymbol{w}}^{\rm T}{\mathbf{A}}({\mathbf{x}}){\boldsymbol{w}}} ,\\ \| {\boldsymbol{w}} \|_{H^1(\Gamma_h[{\mathbf{x}}])} &= \sqrt{{\boldsymbol{w}}^{\rm T}\bigl({\mathbf{M}}({\mathbf{x}})+{\mathbf{A}}({\mathbf{x}})\bigr){\boldsymbol{w}}}, \end{align}\tag{40}\] for the matrices \({\mathbf{M}}({\mathbf{x}})\) and \({\mathbf{A}}({\mathbf{x}})\) defined in 10 .

6.1.2 Finite element interpolation of position and velocity↩︎

Let \({\mathbf{x}}^*(t)=(x_1^*(t), \cdots, x_N^*(t))\) be the image of \({\mathbf{x}}^0\) under the exact flow map \(X(\cdot,t):\Gamma^0\rightarrow\Gamma(t)\), i.e., \(x_j^*(t)=X(x_j^0,t)\) for \(j=1,\dots,N\). Thus the nodal vector \({\mathbf{x}}^*(t)\) determines a triangulated surface \(\Gamma_h[{\mathbf{x}}^*(t)]\) which interpolates \(\Gamma(t)\). The nodal values of \(v(\cdot,t)\) at the interpolated nodes \(x_j^*(t)\), \(j=1,\dots,N\), are collected into the nodal vector \({\mathbf{v}}^*(t)\). We note that \(\dot{\mathbf{x}}^*={\mathbf{v}}^*\). We denote \({\mathbf{x}}_*^n={{\mathbf{x}}}^*(t_n)\) and \({\mathbf{v}}_*^n={{\mathbf{v}}}^*(t_n)\).

6.1.3 Ritz projection of normal vector and mean curvature↩︎

For the solution \(u=(n;H) \in H^1(\varGamma[X])^4\), we denote by \(u_h^* \in S_h[{\mathbf{x}}^\ast]^4\) the Ritz projection of \(u\), i.e., the unique finite element solution of the following weak formulation: \[\label{eq:Ritzmap} \begin{align} & \int_{\varGamma_h[{\mathbf{x}}^\ast]} \!\! u_h^* \, \varphi_h + \int_{\varGamma_h[{\mathbf{x}}^\ast]} \!\! \nabla_{\varGamma_h[{\mathbf{x}}^\ast]} u^*_h \cdot \nabla_{\varGamma_h[{\mathbf{x}}^\ast]} \varphi_h \\ &= \int_{\varGamma[X]} \!\! u \varphi_h^\ell + \int_{\varGamma[X]} \!\! \nabla_{\varGamma[X]} u \cdot \nabla_{\varGamma[X]} \varphi_h^\ell \qquad \forall\, \varphi_h \in S_h[{\mathbf{x}}^\ast]^4. \end{align}\tag{41}\] Such Ritz projections were used in [40], [41], [42] and played an important role in the error analysis of [17]. We recall that (cf. [43]) \[\begin{align} \tag{42} \Vert (u^*_h)^\ell - u \Vert_{L^2(\varGamma[X])^4} + h \Vert (u^*_h)^\ell - u \Vert_{H^1(\varGamma[X])^4} \leq C h^{k+1}, \\ \Vert (\partial^\bullet_h u^*_h)^\ell - \partial^\bullet u \Vert_{L^2(\varGamma[X])^4} + h \Vert (\partial^\bullet_h u^*_h)^\ell - \partial^\bullet u \Vert_{H^1(\varGamma[X])^4} \leq C h^{k+1}. \tag{43} \end{align}\] Here \(\ell\) is the lift from the interpolated surface to the exact surface, as introduced for linear and higher-order surface approximations in [31] and [32], respectively. Let \({\mathbf{u}}_\ast^n={{\mathbf{u}}}^*(t_n)\) be the nodal vector associated with the Ritz projection \(u_h^*(\cdot,t_n)\).

6.1.4 Equivalences of neighbouring norms↩︎

We have the norm-equivalence relation with constants independent of \(h\) and \(\tau\) (when \(\tau\) is sufficiently small): \[\|{\mathbf{w}}\|_{H^1(\Gamma_h[{\mathbf{x}}_*(t)])} \simeq \|{\mathbf{w}}\|_{H^1(\Gamma_h[{\mathbf{x}}_*^n])} \quadfor\,\,\, t\in[t_{n-1},t_n] ,\] which holds as long as \(\|{\mathbf{x}}_*(t)-{\mathbf{x}}_*^n\|_{W^{1,\infty}(\Gamma_h[{\mathbf{x}}_*^n])} \le c\tau\le \frac{1}{2}\) [20]. For the same reason, we also have the following norm equivalence: \[\label{norm-equiv-n-j} \|{\mathbf{w}}\|_{H^1(\Gamma_h[{\mathbf{x}}_*^n])} \simeq \|{\mathbf{w}}\|_{H^1(\Gamma_h[{\mathbf{x}}_*^{n-j}])} \quadfor\,\,\, j=1,\dots,q-1.\tag{44}\] Analogous norm equivalences hold for the corresponding \(L^2\)-norms.

6.2 Error equations↩︎

We can write down the following matrix-vector equations satisfied by the interpolated values \({\mathbf{x}}_\ast^n={{\mathbf{x}}}^*(t_n)\), \({\mathbf{v}}_\ast^n={{\mathbf{v}}}^*(t_n)\) and Ritz-projected values \({\mathbf{u}}_\ast^n= ({{\mathbf{n}}}_*^n;{\mathbf{H}}_*^n)={{\mathbf{u}}}^*(t_n)\) of the exact solution: \[\begin{align} {\mathbf{M}}(\widetilde{{\mathbf{x}}}_*^n) {\boldsymbol{\dot{\mathbf{u}}}}_*^n + {\mathbf{A}}(\widetilde{{\mathbf{x}}}_*^n) {\mathbf{u}}_*^n = &\; {\mathbf{f}}(\widetilde{{\mathbf{x}}}_*^n,\widetilde{{\mathbf{u}}}_*^n) + {\mathbf{M}}(\widetilde{{\mathbf{x}}}_*^n){\mathbf{d}}_{\mathbf{u}}^n , \\[3pt] {\mathbf{v}}_*^n = &\;- {\mathbf{H}}_*^n \bullet {{\mathbf{n}}}_*^n + {\mathbf{d}}_{\mathbf{v}}^n , \\ \frac{1}{\tau} \big( {\mathbf{x}}_*^n - {\mathbf{x}}_*^{n-1} \big) = &\; \bar{\mathbf{v}}_{*}^n + {{\mathbf{d}}}_{{\mathbf{x}}}^n , \end{align}\] where \({\boldsymbol{\dot{\mathbf{u}}}}_*^n = \tau^{-1} \sum_{j=0}^q \delta_j {\mathbf{u}}_*^{n-j}\) and \(\widetilde{{\mathbf{x}}}_{*}^n=\sum_{j=0}^{q-1} \gamma_j {\mathbf{x}}_*^{n-1-j}\) and \(\bar{\mathbf{v}}_{*}^n=\sum_{j=0}^q \beta_j {\mathbf{v}}_*^{n-j}\), and where \({\mathbf{d}}_{\mathbf{u}}^n\) and \({\mathbf{d}}_{\mathbf{x}}^n\) correspond to the defects defined in [17], and \({\mathbf{d}}_{\mathbf{v}}^n\) is the defect as in [23]. The following estimates have been shown in [17] for \({\mathbf{d}}_{\mathbf{u}}^n\), in [23] for \({\mathbf{d}}_{\mathbf{v}}^n\), and \({\mathbf{d}}_{\mathbf{x}}^n\) is the consistency error of the Adams multistep method: \[\begin{align} \label{defects-du-dx-uv} \| {\mathbf{d}}_{\mathbf{u}}^n \|_{L^2(\Gamma_h[{\mathbf{x}}_*^n])} + \| {\mathbf{d}}_{\mathbf{x}}^n \|_{H^1(\Gamma_h[{\mathbf{x}}_*^n])} + \| {\mathbf{d}}_{\mathbf{v}}^n \|_{H^1(\Gamma_h[{\mathbf{x}}_*^n])} \le c \, (\tau^q+h^k) . \end{align}\tag{45}\]

In [17] and [23], the defect estimates were shown in the \(L^2(\Gamma_h[\widetilde{{\mathbf{x}}}_*^n])\) and \(H^1(\Gamma_h[\widetilde{{\mathbf{x}}}_*^n])\) norms. The change to the \(L^2(\Gamma_h[{\mathbf{x}}_*^n])\) and \(H^1(\Gamma_h[{\mathbf{x}}_*^n])\) norms made here is possible because of the equivalence relation in 44 .

We introduce the fully discrete errors at time \(t_n = n \tau\): \[{\mathbf{e}}_{\mathbf{x}}^n = {\mathbf{x}}^n - {\mathbf{x}}_\ast^n, \quad {\mathbf{e}}_{\mathbf{v}}^n = {\mathbf{v}}^n - {\mathbf{v}}_\ast^n , \quad {\mathbf{e}}_{\mathbf{u}}^n = {\mathbf{u}}^n - {\mathbf{u}}_\ast^n .\] The error equations can then be written as \[\tag{46} \begin{align} \tag{47} {\mathbf{M}}(\widetilde{\mathbf{x}}^n){\boldsymbol{\dot{\mathbf{e}}}}_{{\mathbf{u}}}^n+ {\mathbf{A}}(\widetilde{\mathbf{x}}^n){\mathbf{e}}_{\mathbf{u}}^n = &\;- \big( {\mathbf{M}}(\widetilde{\mathbf{x}}^n)-{\mathbf{M}}(\widetilde{\mathbf{x}}_*^n) \big) {\boldsymbol{\dot{\mathbf{u}}}}_*^n \nonumber \\[3pt] &\;- \big( {\mathbf{A}}(\widetilde{\mathbf{x}}^n)-{\mathbf{A}}(\widetilde{\mathbf{x}}_*^n) \big) {\mathbf{u}}_*^n \nonumber \\[3pt] &\;+ \big({\mathbf{f}}(\widetilde{\mathbf{x}}^n,\widetilde{\mathbf{u}}^n) - {\mathbf{f}}(\widetilde{\mathbf{x}}_*^n,\widetilde{\mathbf{u}}_*^n)\big) - {\mathbf{M}}(\widetilde{{\mathbf{x}}}_\ast^n) {\mathbf{d}}_{\mathbf{u}}^n \\[2pt] \tag{48} {\mathbf{e}}_{\mathbf{v}}^n = &\;-\Big( {\mathbf{H}}^n \bullet {\mathbf{n}}^n - {\mathbf{H}}_\ast^n \bullet {\mathbf{n}}_\ast^n \Big) - {\mathbf{d}}_{\mathbf{v}}^n , \\ \tag{49} \frac{1}{\tau} \big( {\mathbf{e}}_{\mathbf{x}}^n - {\mathbf{e}}_{\mathbf{x}}^{n-1} \big) = &\; \bar{\mathbf{e}}_{{\mathbf{v}}}^n + \lambda^n \bar{\mathbf{v}}^n - {{\mathbf{d}}}_{{\mathbf{x}}}^n , \end{align}\] where \(\bar{\mathbf{e}}_{{\mathbf{v}}}^n=\sum_{j=0}^q \beta_j {\mathbf{e}}_{\mathbf{v}}^{n-j}\).

6.3 Stability estimates for \({\mathbf{e}}_{\mathbf{u}}^n\) and \({\mathbf{e}}_{\mathbf{v}}^n\)↩︎

The proof of the error bounds uses an induction argument. We assume that for some positive integer \(m\le T/\tau\) the following estimate holds for all integers \(n\) with \(1\le n\le m\): \[\begin{align} \label{math-ind} & \| {\mathbf{e}}_{\mathbf{u}}^{n-1} \|_{H^1(\varGamma_h[{\mathbf{x}}_\ast^{n-1}])} + \| {\mathbf{e}}_{\mathbf{x}}^{n-1} \|_{H^1(\varGamma_h[{\mathbf{x}}_\ast^{n-1}])} + \| {\mathbf{e}}_{\mathbf{v}}^{n-1} \|_{H^1(\varGamma_h[{\mathbf{x}}_\ast^{n-1}])} + |\lambda^{n-1}| \notag\\ & \le h^{3/2} , \end{align}\tag{50}\] where \(\lambda^i = 0\) for \(i=0,\dotsc,q-1\). Under the assumptions of Theorem 1, for \(\tau^q\le ch^2\) and sufficiently small \(h\), 50 holds at least for \(m=q\). We shall prove that if this result holds for some \(q\le m<T/\tau\) then it also holds for \(m+1\). This will prove that 50 holds for all \(q\le m\le T/\tau\).

Note that by applying the finite element inverse estimates (see, e.g., [44]), the following estimates for \(n\le m\) follows from 50 : \[\begin{align} \label{math-ind-W1inf} & \| {\mathbf{e}}_{\mathbf{u}}^{n-1} \|_{W^{1,\infty}(\varGamma_h[{\mathbf{x}}_\ast^{n-1}])} + \| {\mathbf{e}}_{\mathbf{x}}^{n-1} \|_{W^{1,\infty}(\varGamma_h[{\mathbf{x}}_\ast^{n-1}])} + \| {\mathbf{e}}_{\mathbf{v}}^{n-1} \|_{W^{1,\infty}(\varGamma_h[{\mathbf{x}}_\ast^{n-1}])} \notag\\ & \le c_{1,\infty} h^{1/2} . \end{align}\tag{51}\] Under this condition, stability estimates for 47 and 48 for \(q\le n\le m\) have been established in [17]: \[\begin{align} \label{stability-eu-2} \|{\mathbf{e}}_{\mathbf{u}}^n\|_{H^1(\Gamma_h[{\mathbf{x}}_*^n])}^2 \le &\;c\sum_{j=0}^{q-1} \| {\mathbf{e}}_{\mathbf{x}}^{n-1-j} \|_{H^1(\Gamma_h[{\mathbf{x}}_*^{n-1-j}])}^2 \notag\\ &\;+ c\tau\sum_{j=q}^{n-1} \Bigl( \|{\mathbf{e}}_{\mathbf{u}}^j\|_{H^1(\Gamma_h[{\mathbf{x}}_*^{j}])}^2 + \|{\mathbf{e}}_{\mathbf{x}}^j\|_{H^1(\Gamma_h[{\mathbf{x}}_*^{j}])}^2 + \|{\mathbf{e}}_{\mathbf{v}}^j\|_{H^1(\Gamma_h[{\mathbf{x}}_*^{j}])}^2 \Bigr) \notag\\ &\;+ c \sum_{j=0}^{q-1} \Bigl( \|{\mathbf{e}}_{\mathbf{u}}^{j}\|_{H^1(\Gamma_h[{\mathbf{x}}_*^{j}])}^2 + \|{\mathbf{e}}_{\mathbf{x}}^{j}\|_{H^1(\Gamma_h[{\mathbf{x}}_*^{j}])}^2 \Bigr) \notag\\ &\;+ c \tau \sum_{j=1}^{n-1} \|({\mathbf{e}}_{\mathbf{x}}^j -{\mathbf{e}}_{\mathbf{x}}^{j-1})/\tau \|_{H^1(\Gamma_h[{\mathbf{x}}_*^{j}])}^2 \notag\\ &\;+ c \tau \sum_{j=q}^{n}\|{\mathbf{d}}_{\mathbf{u}}^j\|_{L^2(\Gamma_h[{\mathbf{x}}_*^j])}^2 \qquadfor\, q\le n\le m. \end{align}\tag{52}\] We note that in [17] this error bound is formulated in the \(H^1(\Gamma_h[\widetilde{{\mathbf{x}}}^n])\) norm. In view of 51 and the norm equivalence discussed in Section 6.1, the error bound (with different constants) is valid also for the \(H^1(\Gamma_h[{\mathbf{x}}_*^n])\) norm as stated here.

The stability estimates for 48 are obtained by taking the \(H^1(\Gamma_h[{\mathbf{x}}_*^n])\) norm of both sides of 48 . Using [23] together with 51 , we obtain \[\begin{align} \label{stability-ev} \|{\mathbf{e}}_{\mathbf{v}}^n \|_{H^1(\Gamma_h[{\mathbf{x}}_*^n])} &\le c\|{\mathbf{e}}_{\mathbf{u}}^n\|_{H^1(\Gamma_h[{\mathbf{x}}_*^n])} + \| {\mathbf{d}}_{\mathbf{v}}^n \|_{H^1(\Gamma_h[{\mathbf{x}}_*^n])} \quadfor\, q\le n\le m. \end{align}\tag{53}\]

6.4 Bound for \(\lambda^n\)↩︎

We recall the definition of \(\delta^n\) in 30 , viz., \[\delta^n = \left| \frac{1}{\tau} \bigl(\mathcal{E}_h[\bar{\mathbf{x}}^n] - \mathcal{E}_h[{\mathbf{x}}^{n-1}]\bigr) + \frac{1}{\tau} \int_{t_{n-1}}^{t_n} p^n(t)^2 \,\textrm{d}t \right|,\] and we define the analogous quantity along the exact solution, \[\label{delta-star} \delta^n_* = \left| \frac{1}{\tau} \bigl(\mathcal{E}_h[\bar{\mathbf{x}}^n_*] - \mathcal{E}_h[{\mathbf{x}}^{n-1}_*]\bigr) + \frac{1}{\tau} \int_{t_{n-1}}^{t_n} p^n_*(t)^2 \,\textrm{d}t \right|,\tag{54}\] where \(\bar{\mathbf{x}}^n_*={\mathbf{x}}^{n-1}_*+\tau \bar{\mathbf{v}}^n_*\) and \(p^n_*(t)\) is the interpolation polynomial of the values \(\|H(.,t_{n-j})\|_{L^2(\Gamma(t_{n-j}))}\) for \(j=0,\dots,q-1\). Under the smoothness assumptions of Theorem 1, we have \[\label{eq:delta95star32bound} \delta^n_* \le c \, (h^k+\tau^q),\tag{55}\] as follows from 4 and standard interpolation error bounds.

We prove the following lemma, which is based on Proposition [prop:lambda].

6.5 Stability estimate for \({\mathbf{e}}_{\mathbf{x}}^n\)↩︎

Inserting the bound of Lemma 1 for \(\lambda^n\) into 49 and using the equivalence of shifted norms yields \[\begin{align} \label{stability-ex} \| {\mathbf{e}}_{\mathbf{x}}^n \|_{H^1(\varGamma_h[{\mathbf{x}}_*^n])} &\le (1+c\tau) \| {\mathbf{e}}_{\mathbf{x}}^{n-1} \|_{H^1(\varGamma_h[{\mathbf{x}}_*^{n-1}])} \\ \nonumber &\quad + c\tau \sum_{j=0}^q \Bigl( \| {\mathbf{e}}_{\mathbf{v}}^{n-j} \|_{H^1(\Gamma_h[{\mathbf{x}}_*^{n-j}])} + \| {\mathbf{e}}_{\mathbf{u}}^{n-j} \|_{L^2(\Gamma_h[{\mathbf{x}}_*^{n-j}])} \Bigr) + c\tau \delta^n_*. \end{align}\tag{57}\]

6.6 Combined stability estimate↩︎

Combining the stability estimates 52 , 53 , ?? and 57 , and using a discrete Gronwall inequality, norm equivalences and the induction argument in the same way as in [17], we arrive at the following stability result, which bounds the errors in terms of defects and initial errors. Up to different constants and the extra defects \(\delta^n_*\) of 54 , this is the same as in [17].

6.7 Error bound of Theorem 1↩︎

Together with the known \(O(h^k+\tau^q)\) bounds for the defects [17], the bound 55 , and the interpolation and Ritz projection errors in the case of smooth solutions, Proposition [proposition:stability32-32coupled32problem32-32full] yields the error bounds of Theorem 1 for \(q\leq n\leq m\) in the same way as is shown in [17]. For sufficiently small \(h\) and \(\tau\leq c_0h\), applying the inverse estimates of finite element functions recovers 50 for \(n=m+1\). This proves the mathematical induction on 50 and therefore completes the proof of Theorem 1.0◻

7 Proof of Theorem 2↩︎

The proof of Theorem 2 follows the same line as that of Theorem 1, except for the stability estimate for the velocity. We therefore highlight this difference and omit the remaining, essentially identical arguments.

7.1 Reformulation of the MDR equations↩︎

We define a discrete Laplace–Beltrami operator \(\Delta_{h,\Gamma_h[{\mathbf{x}}]}: S_h[ {\boldsymbol{x}}]\rightarrow S_h[ {\boldsymbol{x}}]\) via duality by \[\int_{\Gamma_h[{\mathbf{x}}]} \Delta_{h,\Gamma_h[ {\boldsymbol{x}}]} v_h \, w_h = -\int_{\Gamma_h[{\mathbf{x}}]} \nabla_{ \Gamma_h[ {\boldsymbol{x}}]} v_h \cdot \nabla_{ \Gamma_h[ {\boldsymbol{x}}]} w_h \quad \forall w_h\in S_h[{\mathbf{x}}] .\] Let \(P_{\Gamma_h[{\mathbf{x}}]} : L^2(\Gamma_h[{\mathbf{x}}])\rightarrow S_h[{\mathbf{x}}]\) be the \(L^2\)-orthogonal projection onto the finite element space. Then, replacing the test function \(\chi_h\) in 20 by \(-\Delta_{h,\Gamma_h[{\mathbf{x}}]}\chi_h\), we obtain \[\begin{align} \label{grad-v46n46eq} & \int_{\Gamma_h[{\mathbf{x}}]} \nabla_{\Gamma_h[{\mathbf{x}}]} (v_h \cdot n_h) \cdot \nabla_{\Gamma_h[{\mathbf{x}}]} \chi_h \notag\\ =& - \int_{\Gamma_h[{\mathbf{x}}]} \nabla_{\Gamma_h[{\mathbf{x}}]} H_h \cdot \nabla_{\Gamma_h[{\mathbf{x}}]} \chi_h \notag \\ & + \int_{\Gamma_h[{\mathbf{x}}]} \nabla_{\Gamma_h[{\mathbf{x}}]} ( v_h \cdot n_h- P_{\Gamma_h[{\mathbf{x}}]} (v_h \cdot n_h)) \cdot \nabla_{\Gamma_h[{\mathbf{x}}]} \chi_h . \end{align}\tag{62}\] This weak formulation can be written as the following matrix-vector form: \[\begin{align} \label{Mform-Dv} {\boldsymbol{D}}({\mathbf{x}},{\mathbf{n}}){\mathbf{v}}&= -{\boldsymbol{A}}({\mathbf{x}}){\boldsymbol{H}} +{\boldsymbol{\rho}}({\mathbf{x}},{\mathbf{n}},{\mathbf{v}}) . \end{align}\tag{63}\] where \({\boldsymbol{\rho}}({\mathbf{x}},{\mathbf{n}},{\mathbf{v}})\in \mathbb{R}^N\) is the nodal vector satisfying the following relation: \[\label{IntroL2} {\boldsymbol{\rho}}({\mathbf{x}},{\mathbf{n}},{\mathbf{v}})\cdot {\boldsymbol{\chi}} = \int_{\Gamma_h[{\mathbf{x}}]} \nabla_{\Gamma_h[{\mathbf{x}}]} [ v_h \cdot n_h- P_{\Gamma_h[{\mathbf{x}}]} (v_h \cdot n_h)] \cdot \nabla_{\Gamma_h[{\mathbf{x}}]} \chi_h\tag{64}\] for all \(\chi_h\in S_h[{\mathbf{x}}]\) with nodal vector \(\boldsymbol{\chi}\in\mathbb{R}^N\), and \[\boldsymbol{\chi}^T D(\mathbf{x},\mathbf{n})\boldsymbol{\psi} = \int_{\Gamma_h[\mathbf{x}]} \nabla_{\Gamma_h[\mathbf{x}]}(\psi_h\cdot n_h) \cdot \nabla_{\Gamma_h[\mathbf{x}]}\chi_h \qquad \forall \boldsymbol{\chi}\in\mathbb{R}^N,\; \boldsymbol{\psi}\in\mathbb{R}^{3N}\] for all \(\chi_h\in S_h[{\mathbf{x}}]\) and \(\psi_h\in S_h[{\mathbf{x}}]^3\) with nodal vectors \(\boldsymbol{\chi}\in\mathbb{R}^N\) and \(\boldsymbol{\psi}\in\mathbb{R}^{3N}\), respectively. The system in 21 is used in the implementation, whereas 63 is a direct consequence of 20 and will be employed in the error analysis.

7.2 Defects and error equations↩︎

Let \(u_h^*(t) \in S_h[{\mathbf{x}}^\ast(t)]^4\) be the Ritz projection of \(u\), and let \(v_h^*(t) \in S_h[{\mathbf{x}}^\ast(t)]^3\) be the Ritz projection of \(v\). Their nodal vectors are denoted by \({\mathbf{u}}^*(t)\) and \({\mathbf{v}}^*(t)\), respectively. For the fully discrete FEM with tangential motion defined in 35 , the interpolated values \({\mathbf{x}}_\ast^n={{\mathbf{x}}}^*(t_n)\) and \(\boldsymbol{\mu}_*^n=\boldsymbol{\mu}^*(t_n)\), and the Ritz projections \({\mathbf{v}}_\ast^n={{\mathbf{v}}}^*(t_n)\) and \({\mathbf{u}}_\ast^n={{\mathbf{u}}}^*(t_n)\), satisfy the following matrix-vector equations: \[\label{eq:MDR-defects} \begin{align} {\mathbf{M}}(\widetilde{{\mathbf{x}}}_*^n) {\boldsymbol{\dot{\mathbf{u}}}}_*^n + {\mathbf{A}}(\widetilde{{\mathbf{x}}}_*^n) {\mathbf{u}}_*^n = &\; {\mathbf{f}}(\widetilde{{\mathbf{x}}}_*^n,\widetilde{{\mathbf{u}}}_*^n,\widetilde{{\mathbf{v}}}_*^n) + {\mathbf{M}}(\widetilde{{\mathbf{x}}}_*^n) {\mathbf{d}}_{\mathbf{u}}^n , \\[3pt] {\boldsymbol{A}}(\widetilde{{\mathbf{x}}}_*^n){\mathbf{v}}_*^n + {\boldsymbol{B}}(\widetilde{{\mathbf{x}}}_*^n,{{\mathbf{n}}}_*^n)^{\rm T}\boldsymbol{\mu}_*^n = &\;{\mathbf{M}}(\widetilde{{\mathbf{x}}}_*^n){\boldsymbol{d}}_{\boldsymbol{\mu}}^n , \\ {\boldsymbol{B}}(\widetilde{{\mathbf{x}}}_*^n,{{\mathbf{n}}}_*^n) {\mathbf{v}}_*^n = &\;- {\boldsymbol{M}}(\widetilde{{\mathbf{x}}}_*^n){\mathbf{H}}_*^n + {\mathbf{M}}(\widetilde{{\mathbf{x}}}_*^n){\boldsymbol{d}}_{\boldsymbol{v}}^n , \\ \frac{1}{\tau} \big( {\mathbf{x}}_*^n - {\mathbf{x}}_*^{n-1} \big) = &\;\bar{\mathbf{v}}_{*}^n + {{\mathbf{d}}}_{{\mathbf{x}}}^n , \end{align}\tag{65}\] where defect terms \({{\mathbf{d}}}_{{\mathbf{x}}}^n\) and \({\mathbf{d}}_{\mathbf{u}}^n\) have the same estimates as that in the method without tangential motion (though it is caused by the Ritz projection now): \[\begin{align} \| {\mathbf{d}}_{\mathbf{u}}^n \|_{L^2(\Gamma_h[{\mathbf{x}}_*^n])} + \| {{\mathbf{d}}}_{{\mathbf{x}}}^n \|_{H^1(\Gamma_h[{\mathbf{x}}_*^n])} \le c \, (\tau^q+h^k) , \end{align}\] while the defect terms \({\mathbf{d}}_{\mathbf{v}}^n\) and \({\boldsymbol{d}}_{\boldsymbol{\mu}}^n\) are defined in [14], satisfying the following estimates: \[\begin{align} \| {\mathbf{d}}_{\mathbf{v}}^n \|_{L^2(\Gamma_h[{\mathbf{x}}_*^n])} + \| {\boldsymbol{d}}_{\boldsymbol{\mu}}^n\|_{L^2(\Gamma_h[{\mathbf{x}}_*^n])} \le c \, (\tau^q+h^k) . \end{align}\] This estimate was established in [14] for the semi-discrete (space-discretized) scheme. In our fully discrete setting, we include an additional factor \(\tau^q\) to account for the truncation error introduced by the time discretization.

Similarly, \(({\mathbf{x}}_*^n,{\mathbf{u}}_*^n,{\mathbf{v}}_*^n)\) satisfies the following relation corresponding to 63 : \[\begin{align} \label{Mform-dv} {\boldsymbol{D}}(\widetilde{{\mathbf{x}}}_*^n,{\mathbf{n}}_*^n){\mathbf{v}}_*^n &= -{\boldsymbol{A}}(\widetilde{{\mathbf{x}}}_*^n){\boldsymbol{H}}_*^n +{\boldsymbol{\rho}}(\widetilde{{\mathbf{x}}}_*^n,{\mathbf{n}}_*^n,{\mathbf{v}}_*^n) + {\mathbf{M}}(\widetilde{{\mathbf{x}}}_*^n) \widetilde{\mathbf{d}}_{{\mathbf{v}}}^n. \end{align}\tag{66}\] where the defect term, \(\widetilde{\mathbf{d}}_{\mathbf{v}}^n\), is a nodal vector satisfying the following relation, obtained by substituting the nodal vectors \({\boldsymbol{v}}_*^n\), \({\boldsymbol{n}}_*^n\), \({\boldsymbol{H}}_*^n\) and \(\widetilde{{\mathbf{x}}}_*^n\) (viewed as finite element functions on \(\Gamma_h[\widetilde{{\mathbf{x}}}_*^n]\)) into 62 : \[\begin{align} \label{tdv-weak} ({\mathbf{M}}(\widetilde{{\mathbf{x}}}_*^n) \widetilde{\mathbf{d}}_{\mathbf{v}}^n)^{\rm T} \boldsymbol{\chi} := &\;\int_{\Gamma_h[\widetilde{{\mathbf{x}}}_*^n]} \nabla_{\Gamma_h[\widetilde{{\mathbf{x}}}_*^n]} P_{\Gamma_h[\widetilde{{\mathbf{x}}}_*^n]} ({\boldsymbol{v}}_*^n \cdot {{\mathbf{n}}}_*^n) \cdot \nabla_{\Gamma_h[\widetilde{{\mathbf{x}}}_*^n]} \boldsymbol{\chi} \notag \\ & + \int_{\Gamma_h[\widetilde{{\mathbf{x}}}_*^n]} \nabla_{\Gamma_h[\widetilde{{\mathbf{x}}}_*^n]} {\boldsymbol{H}}_*^n \cdot \nabla_{\Gamma_h[\widetilde{{\mathbf{x}}}_*^n]} \boldsymbol{\chi} \quad\forall\, \boldsymbol{\chi}\in\mathbb{R}^N. \end{align}\tag{67}\] The following estimate holds: \[\begin{align} \| \widetilde{{\mathbf{d}}}_{{\mathbf{v}}}^n \|_{H^{-1}(\Gamma_h[\widetilde{{\mathbf{x}}}_*^n])} \le c \, (\tau^q+h^k) . \end{align}\] Again, this result was established in [14] for the semi-discrete (space-discretized) scheme. In our fully discrete setting, we include an additional factor \(\tau^q\) to account for the truncation error introduced by the time discretization.

By considering the difference between 65 and 35 , we obtain the following equations for the error functions \({\mathbf{e}}_{\mathbf{x}}^n = {\mathbf{x}}^n - {\mathbf{x}}_\ast^n\), \({\mathbf{e}}_{\mathbf{v}}^n = {\mathbf{v}}^n - {\mathbf{v}}_\ast^n\), \({\boldsymbol{e}}_{\boldsymbol{\mu}}^n = \boldsymbol{\mu}^n - \boldsymbol{\mu}_*^n\) and \({\mathbf{e}}_{\mathbf{u}}^n = {\mathbf{u}}^n - {\mathbf{u}}_\ast^n\): \[\tag{68} \begin{align} \tag{69} {\mathbf{M}}(\widetilde{\mathbf{x}}^n){\boldsymbol{\dot{\mathbf{e}}}}_{{\mathbf{u}}}^n+ {\mathbf{A}}(\widetilde{\mathbf{x}}^n){\mathbf{e}}_{\mathbf{u}}^n = &\;- \big( {\mathbf{M}}(\widetilde{\mathbf{x}}^n)-{\mathbf{M}}(\widetilde{\mathbf{x}}_*^n) \big) {\boldsymbol{\dot{\mathbf{u}}}}_*^n \notag \\[3pt] &\;- \big( {\mathbf{A}}(\widetilde{\mathbf{x}}^n)-{\mathbf{A}}(\widetilde{\mathbf{x}}_*^n) \big) {\mathbf{u}}_*^n \notag \\[3pt] &\;+ \big({\mathbf{f}}(\widetilde{\mathbf{x}}^n,\widetilde{\mathbf{u}}^n,\widetilde{\mathbf{v}}^n) - {\mathbf{f}}(\widetilde{\mathbf{x}}_*^n,\widetilde{\mathbf{u}}_*^n,\widetilde{\mathbf{v}}_*^n)\big) \notag\\ &\;- {\mathbf{M}}(\widetilde{{\mathbf{x}}}_\ast^n) {\mathbf{d}}_{\mathbf{u}}^n \\[2pt] \tag{70} {\boldsymbol{B}}(\widetilde{{\mathbf{x}}}_*^n,{\mathbf{n}}_*^n)^{\rm T} {\mathbf{e}}_{\boldsymbol{\mu}}^n + {\boldsymbol{A}}(\widetilde{{\mathbf{x}}}_*^n){\mathbf{e}}_{\mathbf{v}}^n =& \;({\boldsymbol{B}}(\widetilde{{\mathbf{x}}}_*^n,{\mathbf{n}}_*^n)^{\rm T}-{\boldsymbol{B}}(\widetilde{{\mathbf{x}}}^n,{\mathbf{n}}^n)^{\rm T}){\boldsymbol{\mu}}^n \notag\\ & + ({\boldsymbol{A}}(\widetilde{{\mathbf{x}}}_*^n)-{\boldsymbol{A}}(\widetilde{{\mathbf{x}}}^n)){\mathbf{v}}^n - {\boldsymbol{M}}(\widetilde{{\mathbf{x}}}_*^n) {\mathbf{d}}_{\boldsymbol{\mu}}^n ,\\[5pt] \tag{71} {\boldsymbol{B}}(\widetilde{{\mathbf{x}}}_*^n,{\mathbf{n}}_*^n){\mathbf{e}}_{\mathbf{v}}^n =&\;({\boldsymbol{B}}(\widetilde{{\mathbf{x}}}_*^n,{\mathbf{n}}_*^n)-{\boldsymbol{B}}(\widetilde{{\mathbf{x}}}^n,{\mathbf{n}}^n)){\mathbf{v}}^n \notag \\ &\;- {\boldsymbol{M}}(\widetilde{{\mathbf{x}}}_*^n){\mathbf{e}}_{\mathbf{H}}^n + ({\boldsymbol{M}}(\widetilde{{\mathbf{x}}}_*^n)-{\boldsymbol{M}}(\widetilde{{\mathbf{x}}}^n)){\mathbf{H}}^n \notag \\ &\;- {\boldsymbol{M}}(\widetilde{{\mathbf{x}}}_*^n) {\mathbf{d}}_{{\mathbf{v}}}^n, \\[3pt] \tag{72} \frac{1}{\tau} \big( {\mathbf{e}}_{\mathbf{x}}^n - {\mathbf{e}}_{\mathbf{x}}^{n-1} \big) = &\; \bar{\mathbf{e}}_{{\mathbf{v}}}^n + \lambda^n \bar{\mathbf{v}}^n - {{\mathbf{d}}}_{{\mathbf{x}}}^n . \end{align}\] Similarly, the error equation corresponding to 66 is given by \[\begin{align} \label{error-ev-grad} {\boldsymbol{D}}(\widetilde{{\mathbf{x}}}_*^n,{\mathbf{n}}_*^n){\mathbf{e}}_{\mathbf{v}}^n =& \;({\boldsymbol{D}}(\widetilde{{\mathbf{x}}}_*^n,{\mathbf{n}}_*^n)-{\boldsymbol{D}}(\widetilde{{\mathbf{x}}}^n,{\mathbf{n}}^n)){\mathbf{v}}^n \notag\\ &\; - {\boldsymbol{A}}(\widetilde{{\mathbf{x}}}_*^n){\mathbf{e}}_{\mathbf{H}}^n + ({\boldsymbol{A}}(\widetilde{{\mathbf{x}}}_*^n)-{\boldsymbol{A}}(\widetilde{{\mathbf{x}}}^n)){\mathbf{H}}^n \notag \\ &\;- \big(\boldsymbol{\rho}(\widetilde{{\mathbf{x}}}_*^n,{\mathbf{n}}_*^n,{\mathbf{v}}_*^n) - \boldsymbol{\rho}(\widetilde{{\mathbf{x}}}^n,{\mathbf{n}}^n,{\mathbf{v}}^n)\big) - {\boldsymbol{M}}(\widetilde{{\mathbf{x}}}_*^n) \widetilde{\mathbf{d}}_{{\mathbf{v}}}^n . \end{align}\tag{73}\]

7.3 Stability estimates↩︎

The stability estimate for 69 is the same as 47 as they have the same structure. Therefore, 52 still holds, i.e., \[\begin{align} \label{stability-MDR-eu-2} \|{\mathbf{e}}_{\mathbf{u}}^n\|_{H^1(\Gamma_h[{\mathbf{x}}_*^n])}^2 \le &\;c\sum_{j=0}^{q-1} \| {\mathbf{e}}_{\mathbf{x}}^{n-1-j} \|_{H^1(\Gamma_h[{\mathbf{x}}_*^{n-1-j}])}^2 \notag\\ &\;+ c\tau\sum_{j=q}^{n-1} \Bigl( \|{\mathbf{e}}_{\mathbf{u}}^j\|_{H^1(\Gamma_h[{\mathbf{x}}_*^{j}])}^2 + \|{\mathbf{e}}_{\mathbf{x}}^j\|_{H^1(\Gamma_h[{\mathbf{x}}_*^{j}])}^2 + \|{\mathbf{e}}_{\mathbf{v}}^j\|_{H^1(\Gamma_h[{\mathbf{x}}_*^{j}])}^2 \Bigr) \notag\\ &\;+ c \sum_{j=0}^{q-1} \Bigl( \|{\mathbf{e}}_{\mathbf{u}}^{j}\|_{H^1(\Gamma_h[{\mathbf{x}}_*^{j}])}^2 + \|{\mathbf{e}}_{\mathbf{x}}^{j}\|_{H^1(\Gamma_h[{\mathbf{x}}_*^{j}])}^2 \Bigr) \notag\\ &\;+ c \tau \sum_{j=1}^{n-1} \|({\mathbf{e}}_{\mathbf{x}}^j -{\mathbf{e}}_{\mathbf{x}}^{j-1})/\tau \|_{H^1(\Gamma_h[{\mathbf{x}}_*^{j}])}^2 \notag\\ &\;+ c \tau \sum_{j=q}^{n} \|{\mathbf{d}}_{\mathbf{u}}^j\|_{L^2(\Gamma_h[{\mathbf{x}}_*^j])}^2 \qquadfor\, q\le n\le m. \end{align}\tag{74}\]

The following stability estimate for 70 71 has been established in [14]: \[\begin{align} \|{\mathbf{e}}_{\mathbf{v}}^n\|_{H^1(\Gamma_h[\widetilde{{\mathbf{x}}}_*^n])} & \le c \, (\|{\mathbf{e}}_{\mathbf{u}}^n\|_{H^1(\Gamma_h[\widetilde{{\mathbf{x}}}_*^n])} + \|\widetilde{{\mathbf{e}}}_{{\mathbf{x}}}^n\|_{H^1(\Gamma_h[\widetilde{{\mathbf{x}}}_*^n])} + \| {\mathbf{d}}_{\mathbf{v}}^n\|_{L^2(\Gamma_h[\widetilde{{\mathbf{x}}}_*^n])}) \notag\\ &\quad\; + c \, (\| {\widetilde{\mathbf{d}}}_{\mathbf{v}}^n\|_{H^{-1}(\Gamma_h[\widetilde{{\mathbf{x}}}_*^n])} + \| {\mathbf{d}}_{\boldsymbol{\mu}}^n\|_{L^2(\Gamma_h[\widetilde{{\mathbf{x}}}_*^n])} ) . \end{align}\] In view of 51 and the norm equivalence discussed in Section 6.1, this error bound (with different constants) is valid also on the surface \(\Gamma_h[{{\mathbf{x}}}_*^n]\), i.e., \[\begin{align} \label{stability-MDR-ev} \|{\mathbf{e}}_{\mathbf{v}}^n\|_{H^1(\Gamma_h[{{\mathbf{x}}}_*^n])} & \le c \, (\|{\mathbf{e}}_{\mathbf{u}}^n\|_{H^1(\Gamma_h[{{\mathbf{x}}}_*^n])} + \sum_{j=1}^q\|{\mathbf{e}}_{{\mathbf{x}}}^{n-j}\|_{H^1(\Gamma_h[{{\mathbf{x}}}_*^n])} + \| {\mathbf{d}}_{\mathbf{v}}^n\|_{L^2(\Gamma_h[{{\mathbf{x}}}_*^n])}) \notag\\ &\quad\; + c \, (\| {\widetilde{\mathbf{d}}}_{\mathbf{v}}^n\|_{H^{-1}(\Gamma_h[\widetilde{{\mathbf{x}}}_*^n])} + \| {\mathbf{d}}_{\boldsymbol{\mu}}^n\|_{L^2(\Gamma_h[{{\mathbf{x}}}_*^n])} ) . \end{align}\tag{75}\]

The stability estimate for 72 is identical to that for 49 , since the two relations have the same structure. Moreover, the stability bound for \(e_\lambda^n\) is unchanged as in Lemma 1. Consequently, the estimate 57 remains valid for the scheme with tangential motion, i.e., \[\begin{align} \label{stability-ex-tan} \| {\mathbf{e}}_{\mathbf{x}}^n \|_{H^1(\varGamma_h[{\mathbf{x}}_*^n])} &\le (1+c\tau) \| {\mathbf{e}}_{\mathbf{x}}^{n-1} \|_{H^1(\varGamma_h[{\mathbf{x}}_*^{n-1}])} \nonumber\\ &\quad + c\tau \sum_{j=0}^q \Bigl( \| {\mathbf{e}}_{\mathbf{v}}^{n-j} \|_{H^1(\Gamma_h[{\mathbf{x}}_*^{n-j}])} + \| {\mathbf{e}}_{\mathbf{u}}^{n-j} \|_{L^2(\Gamma_h[{\mathbf{x}}_*^{n-j}])} \Bigr) \nonumber\\ &\quad + c\tau (h^k+\tau^q). \end{align}\tag{76}\] The remainder of the argument proceeds exactly as in Section 6.6: combining 74 , 75 and 76 yields the desired error bound. This completes the proof of Theorem 2. 0◻

8 Numerical experiments↩︎

We performed numerical simulations and experiments for mean curvature flow using various algorithms of this paper. We report on:

• some important details on implementation;

• convergence experiments with and without a Lagrange-multiplier to illustrate the theoretical results of Theorem 1 and 2;

• efficiency plots comparing the computational overhead for the Lagrange-multiplier algorithm.

The numerical experiments use quadratic evolving surface finite elements, implemented in the Matlab package \(\ell\)FEM [45]. For the computation of all finite element matrices and vectors it uses high-order quadratures so that the resulting quadrature error does not feature in the discussion of the accuracies of the schemes. The temporal discretization uses linearly-implicit BDF–Adams methods of order \(1,\dotsc,5\). The initial meshes were all generated using DistMesh [46], without taking advantage of any symmetry of the surfaces.

8.1 Implementation↩︎

Let us recall the fully-discrete algorithm 25 determining \({\mathbf{x}}^n\), \({\mathbf{v}}^n\), and \({\mathbf{u}}^n = ({\mathbf{n}}^n, {\mathbf{H}}^n)\).

(1) We determine the (normalised) extrapolations: \[\widehat{{\mathbf{n}}}^n = {\mathbf{n}}^n / |{\mathbf{n}}^n| \quad and \quad \widetilde{{\mathbf{u}}}^n = \big( \widehat{ {\widetilde{\mathbf{n}}}}^n , \widetilde{\mathbf{H}}^n \big) .\]

(2) We solve the linear system for the geometric variable \({\mathbf{u}}^n\): \[\Big( \delta_0 {\mathbf{M}}(\widetilde{{\mathbf{x}}}^n) + \tau {\mathbf{A}}(\widetilde{{\mathbf{x}}}^n) \Big) {\mathbf{u}}^n = \tau {\mathbf{f}}(\widetilde{{\mathbf{x}}}^n,\widetilde{{\mathbf{u}}}^n) - {\mathbf{M}}(\widetilde{{\mathbf{x}}}^n) \sum_{j=1}^q \delta_j {\mathbf{u}}^{n-j}.\]

(3) We compute the velocity \({\mathbf{v}}^n\) (using a normalisation once again): \[{\mathbf{v}}^n = - {\mathbf{H}}^n \bullet \widehat{{\mathbf{n}}}^n .\]

(4) We determine \(\lambda^n \in \mathbb{R}\) via a (simplified) Newton iteration, which in turn gives \({\mathbf{x}}^n\): To this end, we start by rewriting 28 , separating Lagrange-parameter-dependent parts, for arbitrary \(\lambda \in \mathbb{R}\): \[\begin{align} {\mathbf{x}}^n(\lambda) := &\;{\mathbf{x}}^n(0) + \lambda \tau {\bar{\mathbf{v}}}^n , \\ \text{with} \qquad &\;{\mathbf{x}}^n(0) := {\mathbf{x}}^{n-1} + \tau {\bar{\mathbf{v}}}^n \quad and \quad{\bar{\mathbf{v}}}^n := \sum_{j=0}^{q} \beta_j {\mathbf{v}}^{n-j} , \end{align}\] in particular \[\begin{align} &\;{\mathbf{x}}^n := {\mathbf{x}}^n(\lambda^n), \quad \text{for a \lambda^n \in \mathbb{R} satisfying} ~ \mathcal{E}^n(\lambda^n) = \mathcal{E}^{n-1} - b^n , \\ &\;\text{with} ~ \mathcal{E}^n(\lambda) := \int_{\varGamma_h[{\mathbf{x}}^n(\lambda)]} 1 \quad and \quad b^n := \int_{t_{n-1}}^{t_n} p^n(t)^2 \textrm{d}t \, . \end{align}\] The integral of the Lagrange interpolation polynomial \(p^n\) is computed exactly, via a high-order Gauss quadrature.

The simplified Newton method for \(\lambda^n\) reads: \[\begin{align} \lambda^{(k+1)} = &\;\lambda^{(k)} + \Delta \lambda^{(k)} ,\\ \partial_\lambda \mathcal{E}^n( 0 ) \, \Delta \lambda^{(k)} = &\;- \big( \mathcal{E}^n( \lambda^{(k)} ) - ( \mathcal{E}^{n-1} - b^n ) \big) , \end{align}\] where the \(\lambda\)-derivative of \(\mathcal{E}^n\colon \mathbb{R}\rightarrow\mathbb{R}\) is given, using the Leibniz formula, by \[\begin{align} \partial_\lambda \mathcal{E}^n( \lambda ) = &\;\partial_\lambda \int_{\varGamma_h[{\mathbf{x}}^n(\lambda)]} 1 = \int_{\varGamma_h[{\mathbf{x}}^n(\lambda)]} \nabla_{\varGamma_h[{\mathbf{x}}^n(\lambda)]} \cdot \big( \tau v_{\beta,h}^n \big) , \\ \partial_\lambda \mathcal{E}^n( 0 ) = &\;\partial_\lambda \mathcal{E}^n( \lambda )|_{\lambda = 0} = \tau \int_{\varGamma_h[{\mathbf{x}}^n(0)]} \nabla_{\varGamma_h[{\mathbf{x}}^n(0)]} \cdot v_{\beta,h}^n , \end{align}\] where \(v_{\beta,h}^n\) denotes the appropriate finite element function with nodal values \({\bar{\mathbf{v}}}^n\).

We use the stopping criteria \(\big| \mathcal{E}^n( \lambda^{(k)} ) - ( \mathcal{E}^{n-1} - b^n ) \big|\leq \delta b^n\) with a small \(\delta \in (0,1)\).

8.2 Comparing the algorithms without tangential motion with and without Lagrange-multiplier↩︎

We perform numerical experiments for mean curvature flow of a sphere, following [17].

We report on convergence experiments for the unit sphere as initial surface \(\varGamma^0\). We start the algorithms from the nodal interpolations of the exact initial values \(\varGamma^0\), \(\nu^0\), and \(H^0 = 2\), subsequent initial values \(i=1,\dotsc,q-1\) were chosen to be the exact values. In order to illustrate the convergence results of Theorem 1, we have computed the errors between the numerical solution 25 and the exact solution of 1 .

In Figure 1 and 2 we report the errors between the numerical solution and the interpolation of the exact solution until the final time \(T=0.24\), for a sequence of meshes (see plots) and for a sequence of time steps \(\tau_{k+1} = \tau_k / 2\). The gray lines correspond to the algorithm without Lagrange multiplier 11 , the black lines correspond to the algorithm with Lagrange multiplier 25 . The lines marked with different symbols correspond to different time step sizes and to different mesh refinements in Figure 1 and 2, respectively. The double-logarithmic plots in the first four panels report on the \(L^\infty(H^1)\) norm of the errors against the mesh width \(h\) in Figure 1, and against the time step size \(\tau\) in Figure 2. The semi-logarithmic plot in the rightmost panel reports on the \(L^\infty\)-norm of the Lagrange-multiplier.

In Figure 1 we can observe two regions: a region where the spatial discretization error dominates, matching the \({\mathcal{O}}(h^2)\) order of convergence of Theorem 1 and [17] (see the reference lines), and a region, with small mesh size, where the temporal discretization error dominates (the error curves flatten out). For Figure 2, the same description applies, but with reversed roles, the errors matching the \({\mathcal{O}}(\tau^3)\) order of convergence of Theorem 1 and [17] (see the reference lines).

In Figure 3 we report on computational efficiency, that is, we plot the errors in the surface and area against the CPU times required by both algorithms. As before, the gray and black lines respectively correspond to the algorithm without and with Lagrange multiplier. The lines correspond to various spatial refinements (see plots) while the nodes on each line correspond to time step sizes \(\tau_i\) (from the convergence test). In Figure 3 it can be observed that the new structure-preserving algorithm 25 only requires a slight overhead when compared to algorithm 11 without Lagrange-multiplier. Furthermore, it is seen that this computational overhead is independent of mesh size and time step size.

Acknowledgements↩︎

The foundations of this work have been laid when the three authors were participating in the Oberwolfach Research Fellows Programme at the Mathematisches Forschungsinstitut Oberwolfach (MFO) in March 2025.

The exchange between Germany and Hong Kong was funded in part by the Germany/Hong Kong Joint Research Scheme sponsored by the Research Grants Council of Hong Kong and the German Academic Exchange Service of Germany (Ref. No. G-PolyU501/24 and DAAD PPP Projekt-ID: 57750652).

The work of B. Kovács was funded by the DFG Heisenberg Programme – Project-ID 446431602, and by the DFG Research Unit FOR 3013 Vector- and tensor-valued surface PDEs (BA2268/6–1).

The work of B. Li was funded in part by the National Natural Science Foundation of China (NSFC Project No. 12525111) and the Research Grants Council of Hong Kong (PolyU/RFS2324-5S03).

References↩︎