\emails

[email protected] (E.Miyazaki), chan@email (A. Chan), zhao@email (A. Zhao)

A structure-preserving numerical method
for the fourth-order geometric evolution equations for planar curves

Eiji Miyazaki 1 Tomoya Kemmochi 1 Tomohiro Sogabeand Shao-Liang Zhang 1 1 ¹¹affiliationmark: Graduate School of Engineering, Nagoya University, Furo-cho, Chikusa-ku, Nagoya, Aichi 464-8601 Japan

Abstract

For fourth-order geometric evolution equations for planar curves with the dissipation of the bending energy, including the Willmore and the Helfrich flows, we consider a numerical approach. In this study, we construct a structure-preserving method based on a discrete variational derivative method. Furthermore, to prevent the vertex concentration that may lead to numerical instability, we discretely introduce Deckelnick’s tangential velocity. Here, a modification term is introduced in the process of adding tangential velocity. This modified term enables the method to reproduce the equations’ properties while preventing vertex concentration. Numerical experiments demonstrate that the proposed approach captures the equations’ properties with high accuracy and avoids the concentration of vertices.

keywords:

Geometric evolution equation, Willmore flow, Helfrich flow, Numerical calculation, Structure-preserving, Discrete variational derivative method, Tangential velocity

\ams

37M15, 65M06

1 Introduction

In this paper, we design a numerical method for the Willmore and the Helfrich flows for planar curves. The Willmore flow is a geometric evolution equation that models the behavior of elastic bodies [1]. It is a gradient flow with regard to the bending energy

B(t)=\frac{1}{2}\int_{C(t)}(k-c_{0})^{2}ds,

(1)

where $C(t)$ is a closed curve, $s$ is the arc-length parameter of $C(t)$ , $k$ is the curvature of $C(t)$ , and $c_{0}$ is a given constant [2, 3]. The Willmore flow is given by

\frac{\partial\bm{X}}{\partial{t}}=-\bm{\delta}B,\quad\bm{\delta}B=-\left\{\frac{\partial^{2}k}{\partial s^{2}}+\frac{1}{2}(k-c_{0})^{3}\right\}\bm{N},

(2)

where $\bm{\delta}B$ is the gradient of $B$ , $\bm{X}$ is a point on the closed curve $C(t)$ , and $\bm{N}$ is the unit outward normal vector of $C(t)$ .

The Helfrich flow is a gradient flow for the bending energy under the constraint that the length and the enclosed area of the curve are conserved. The Helfrich flow is given by

\begin{split}&\frac{\partial\bm{X}}{\partial{t}}=\left\{\frac{\partial^{2}k}{\partial s^{2}}+\frac{1}{2}(k-c_{0})^{3}+k\lambda+\mu\right\}\bm{N},\\ &\begin{pmatrix}\lambda\\ \mu\end{pmatrix}=\frac{1}{\langle k\rangle^{2}-\langle k^{2}\rangle}\begin{pmatrix}1&-\langle k\rangle\\ -\langle k\rangle&-\langle k^{2}\rangle\end{pmatrix}\begin{pmatrix}\langle k\bm{\delta}B\rangle+c_{0}^{2}\langle k^{2}\rangle/2\\ \langle\bm{\delta}B\rangle+c_{0}^{2}\langle k\rangle/2\end{pmatrix},\end{split}

(3)

where

\langle F\rangle=\frac{1}{L}\int_{C(t)}Fds,\quad L=\int_{C(t)}ds.

(4)

Helfrich proposed an optimization problem that mathematically models the shape of red blood cells [4, 5]. The Helfrich flow was proposed to solve this optimization problem [6, 7]. Note that the Willmore and the Helfrich flows are fourth-order nonlinear evolution equations.

Some numerical approaches have been investigated for (2) and (LABEL:helf1). There are two primary problems in the numerical computation of these flows.

•

Numerical computation using general-purpose approaches (e.g. the Runge–Kutta method) can become unstable if the time step size is too large.
•

When one approximates curves by polygonal curves, the vertices may be concentrated as the time step proceeds. The concentration of vertices may make numerical computations unstable.

Some numerical methods regarding the above problems for (2) and (LABEL:helf1) are presented. In [8], a linear numerical method for (2) and (LABEL:helf1) are proposed. It is based on the finite element method and approximates a closed curve by a closed polygonal curve. Additionally, the method implicitly involves a tangential velocity through a mass lumped inner product. In [9], a semi-implicit numerical method for (2) is proposed. The method also uses a polygonal curve and introduces tangential velocity that makes the distribution of vertices asymptotically uniform, which was derived in [10, 11]. These tangential velocities’ impact in [8, 9] can prevent the polygonal curve’s vertices from concentrating. In [12], a structure-preserving numerical method for (2) is proposed. The method approximates the curve by B-spline curves [13]. It is based on the discrete partial derivative method [14] and discretely reproduces the dissipation of the bending energy. The method is a nonlinear scheme. A method that discretely reproduces a property of an equation is called the structure-preserving numerical method. The numerical solutions obtained using these methods are not only physically correct but also known to be less likely to diverge even when relatively large time step sizes are used [15].

This study aims to construct a structure-preserving numerical method with tangential velocity for (2) and (LABEL:helf1). That is, for (2), we construct a numerical method that discretely reproduces the bending energy’s dissipation, and for (LABEL:helf1), we construct a numerical method that discretely reproduces the bending energy’s dissipation, conservation of the length, and the enclosed area.

Throughout this paper, the polygonal curve is used to discretize the curve. Additionally, the variational derivative of the bending energy, the length, and the area are discretized based on the discrete variational derivative method (DVDM) [16]. Then, these variational derivatives and the Lagrange multipliers are used to build a structure-preserving method. Additionally, to prevent vertex concentration, we discretely introduce Deckelnick’s tangential velocity [17]. Here, a modification term is introduced to make our method structure-preserving.

The proposed method enables us to conduct numerical computations with relatively large time step sizes while avoiding the polygonal curve’s vertex concentration. Furthermore, since the energy properties of equation (2) and (LABEL:helf1) are discretely reproduced, it is anticipated that the solutions obtained by using our method are physically correct. The method is a nonlinear scheme. Therefore, the cost of numerical computations tends to be high.

The rest of this paper is organized as follows. In Section 2, we collect the bending energy’s dissipation of (2), the bending energy’s dissipation, and the conservation of the length and the area of (LABEL:helf1). In Section 3, we construct structure-preserving numerical methods based on the DVDM. Additionally, discrete tangential velocity is introduced. We present numerical examples in Section 4, and finally, give concluding remarks in Section 5.

2 Properties of the target equations

In this study, we construct a structure-preserving numerical method by reproducing the process of deriving the properties of the gradient flows in a continuous system. Specifically, we discretize variational derivatives based on DVDM, and discretely reproduce the approach of determining the Lagrange multipliers in the continuous system. In this section, we will collect the variational structure and Lagrange multipliers of target equations, and derive the properties of these equations.

2.1 Willmore flow

The Willmore flow is given by

\frac{\partial\bm{X}}{\partial{t}}=-\bm{\delta}B,\quad\bm{\delta}B=-\left\{\frac{\partial^{2}k}{\partial s^{2}}+\frac{1}{2}(k-c_{0})^{3}\right\}\bm{N},

(5)

where $\bm{\delta}B$ represents the gradient of $B$ , $\bm{X}$ represents a point on a closed curve $C(t)$ , $k$ represents the curvature, $c_{0}$ represents a given constant and $\bm{N}$ represents the unit outward normal vector. The Willmore flow has the bending energy’s dissipation:

\begin{split}\frac{dB}{dt}&=\int_{C(t)}\bm{\delta}B\cdot\frac{\partial\bm{X}}{\partial t}ds\\ &=-\int_{C(t)}\left|\bm{\delta}B\right|^{2}ds\leq 0.\end{split}

(6)

2.2 Helfrich flow

The Helfrich flow is given by

\begin{split}&\frac{\partial\bm{X}}{\partial{t}}=\left\{\frac{\partial^{2}k}{\partial s^{2}}+\frac{1}{2}(k-c_{0})^{3}+k\lambda+\mu\right\}\bm{N},\\ &\begin{pmatrix}\lambda\\ \mu\end{pmatrix}=\frac{1}{\langle k\rangle^{2}-\langle k^{2}\rangle}\begin{pmatrix}1&-\langle k\rangle\\ -\langle k\rangle&-\langle k^{2}\rangle\end{pmatrix}\begin{pmatrix}\langle k\bm{\delta}B\rangle+c_{0}^{2}\langle k^{2}\rangle/2\\ \langle\bm{\delta}B\rangle+c_{0}^{2}\langle k\rangle/2\end{pmatrix},\end{split}

(7)

where

\langle F\rangle=\frac{1}{L}\int_{C(t)}Fds,\quad L=\int_{C(t)}ds.

(8)

The Helfrich flow can be rewritten in the following form using the variational derivative of the bending energy $\bm{\delta}B$ , the variational derivative of the length of $C(t)$ (denoted by $\bm{\delta}L$ ) and the enclosed erea of $C(t)$ (denoted by $\bm{\delta}A$ ):

\frac{\partial\bm{X}}{\partial{t}}=-\bm{\delta}B+\lambda\bm{\delta}L+\mu\bm{\delta}A,\quad\bm{\delta}L=k\bm{N},\quad\bm{\delta}A=\bm{N}.

(9)

Hereafter, for simplicity, we denote the inner product over a curve $C(t)$ by

(\bm{X},\bm{Y})=\int_{C(t)}\bm{X}\cdot\bm{Y}ds.

(10)

We obtain the expression of time derivatives of the length and the area as follows:

	$\displaystyle\frac{dL}{dt}$	$\displaystyle=\int_{C(t)}\bm{\delta}L\cdot\frac{\partial\bm{X}}{\partial t}ds=\left(\bm{\delta}L,\frac{\partial\bm{X}}{\partial t}\right),$		(11)
	$\displaystyle\frac{dA}{dt}$	$\displaystyle=\int_{C(t)}\bm{\delta}A\cdot\frac{\partial\bm{X}}{\partial t}ds=\left(\bm{\delta}A,\frac{\partial\bm{X}}{\partial t}\right).$		(12)

Since the length and the area of $C(t)$ are conserved, we have $\frac{dL}{dt}=\frac{dA}{dt}=0$ . By substituting (9) into (11) and (12), we obtain

		$\displaystyle\left(\bm{\delta}L,-\bm{\delta}B+\lambda\bm{\delta}L+\mu\bm{\delta}A\right)=0,$		(13)
		$\displaystyle\left(\bm{\delta}A,-\bm{\delta}B+\lambda\bm{\delta}L+\mu\bm{\delta}A\right)=0.$		(14)

Solving (13) and (14) as equations for $\lambda$ and $\mu$ , we can obtain $\lambda$ and $\mu$ that conserve the curve length and the area of $C(t)$ .

Additionally, The Helfrich flow also has the dissipation of the bending energy:

\begin{split}\frac{dB}{dt}&=\int_{C(t)}\bm{\delta}B\cdot\frac{\partial\bm{X}}{\partial t}ds\\ &=\int_{C(t)}\bm{\delta}B\cdot\left(-\bm{\delta}B+\lambda\bm{\delta}L+\mu\bm{\delta}A\right)ds\\ &=-\mathrm{det}\left(\begin{matrix}(\bm{\delta}B,\bm{\delta}B)&(\bm{\delta}B,\bm{\delta}L)&(\bm{\delta}B,\bm{\delta}A)\\ (\bm{\delta}L,\bm{\delta}B)&(\bm{\delta}L,\bm{\delta}L)&(\bm{\delta}L,\bm{\delta}A)\\ (\bm{\delta}A,\bm{\delta}B)&(\bm{\delta}A,\bm{\delta}L)&(\bm{\delta}A,\bm{\delta}A)\end{matrix}\right)/\mathrm{det}\left(\begin{matrix}(\bm{\delta}L,\bm{\delta}L)&(\bm{\delta}L,\bm{\delta}A)\\ (\bm{\delta}A,\bm{\delta}L)&(\bm{\delta}A,\bm{\delta}A)\end{matrix}\right)\\ &\leq 0.\end{split}

(15)

In the process of above transformation, we solved (13) and (14) as a linear system for $\lambda$ and $\mu$ , and then substituted them.

3 Construction of structure-preserving
numerical method for each equation

In this section, various quantities related to closed curves are discretized based on the approximation of the closed curve by a closed polygonal curve. Let us discretize each variational derivative according to the concept of DVDM and determine the Lagrange multipliers to preserve the equation’s structure. The tangential velocity is also introduced for stabilization. Note that in the process of introducing the tangential velocity, we add a modification term to preserve the discrete bending energy’s dissipation even after the tangential velocity is added.

3.1 Discretization of the curve

Let us discretize a curve and geometric quantities according to [18]. We represent a time-evolving polygonal curve in the plane by $\Gamma(t)$ and the domain of inside $\Gamma(t)$ by $\Omega(t)$ . Note that in this subsection, $t$ is excluded because we fix the time $t$ and consider various quantities defined on the polygonal curve $\Gamma(t)$ . Let $N$ be the number of edges of $\Gamma$ and $\bm{X}_{i},\;i=1,\dots,N$ be the $i$ th vertex. The index $i$ shall increase counterclockwise.

Next, let $\Gamma_{i}$ be the line segment connecting $\bm{X}_{i-1}$ and $\bm{X}_{i}$ . We represent the length of $\Gamma_{i}$ by $r_{i}=|\bm{X}_{i}-\bm{X}_{i-1}|$ . Then, the unit tangent vector $\bm{t}_{i}$ on $\Gamma_{i}$ is defined by

\bm{t}_{i}=\left(\bm{X}_{i}-\bm{X}_{i-1}\right)/{r_{i}}.

(16)

Additionally, the unit tangent vector $\bm{T}_{i}$ on the vertex $\bm{X}_{i}$ is defined by

\bm{T}_{i}=\frac{\bm{t}_{i}+\bm{t}_{i+1}}{|\bm{t}_{i}+\bm{t}_{i+1}|}.

(17)

Then, the curvature at the $i$ th vertex $\bm{X}_{i}$ is defined by

k_{i}=\frac{\mathrm{det}\left[\delta_{i}^{\left<1\right>}\bm{X}_{i},\delta_{i}^{\left<2\right>}\bm{X}_{i}\right]}{|\delta_{i}^{\left<1\right>}\bm{X}_{i}|^{3}},

(18)

where

\delta_{i}^{\left<1\right>}\bm{X}_{i}=\frac{\bm{X}_{i+1}-\bm{X}_{i-1}}{2\Delta u},\quad\delta_{i}^{\left<2\right>}\bm{X}_{i}=\frac{\bm{X}_{i+1}-2\bm{X}_{i}+\bm{X}_{i-1}}{\Delta u},\quad\Delta u=1/N.

(19)

Here, $\mathrm{det}[\cdot,\cdot]$ means the following operations for two vectors $\bm{A}$ and $\bm{B}$ :

\bm{A}=\left(\begin{matrix}a_{1}\\ a_{2}\end{matrix}\right),\quad\bm{B}=\left(\begin{matrix}b_{1}\\ b_{2}\end{matrix}\right)\implies\mathrm{det}[\bm{A},\bm{B}]=\mathrm{det}\left[\begin{matrix}a_{1}&b_{1}\\ a_{2}&b_{2}\end{matrix}\right].

(20)

As above, the curvature is discretized by the central differencing of the derivatives. Then, using discretized curvature (18), the discrete bending energy is defined by

B_{d}(\Gamma)=\frac{1}{2}\sum_{i=1}^{N}(k_{i}-c_{0})^{2}\hat{r}_{i},

(21)

where $\hat{r}_{i}=(r_{i}+r_{i+1})/2$ . The length $L_{d}$ and the enclosed area $A_{d}$ of $\Gamma$ are defined by

\displaystyle L_{d}(\Gamma)=\sum_{i=1}^{N}r_{i},\quad A_{d}(\Gamma)=\frac{1}{2}\sum_{i=1}^{N}\mathrm{det}\left[\bm{X}_{i-1},\bm{X}_{i}\right].

(22)

3.2 Discrete tangential velocity

In this section, discrete tangential velocity is defined. It is known that tangential velocity in the continuous system does not influence the shape of the curve [19]. Thus, those are used as a means of examining the equations [17]. In a discrete system, if the number of vertices $N$ is sufficiently large, the impact on the shape of the polygonal curve is considered small. Numerical computation of time-evolving curves usually leads to concentration of vertices, which induces numerical instability if the tangent velocity is not taken properly. Various numerical methods have demonstrated the efficiency of introducing discrete tangential velocities[10, 11, 18, 20].

In this study, the tangent velocity of Deckelnick[17] is used in the discrete system to prevent vertex concentration. Deckelnick’s tangential velocity is given by

w(u,t)=-\frac{\partial}{\partial u}\left(\frac{1}{|\partial_{u}\bm{X}|}\right).

(23)

We discretize (23) as follows

w_{i}(t)=-\alpha\delta_{i}^{-}\left(\frac{1}{|\delta_{i}^{+}\bm{X}_{i}(t)|}\right),

(24)

where

\quad\delta_{i}^{+}\bm{X}_{i}=\frac{\bm{X}_{i+1}-\bm{X}_{i}}{\Delta u},\quad\delta_{i}^{-}\bm{X}_{i}=\frac{\bm{X}_{i}-\bm{X}_{i-1}}{\Delta u},

and $w_{i}(t)$ represents the discretized tangential velocity at the vertex $\bm{X}_{i}(t)$ , and $\alpha$ is a positive constant. The tangential velocity for a vertex $\bm{X}_{i}(t)$ reduces the difference in distance to its two neighboring vertices $\bm{X}_{i-1}(t)$ and $\bm{X}_{i+1}(t)$ . The larger the value of $\alpha$ , the greater the impact of tangential velocity. We must modify the value of $\alpha$ according to the equations and conditions of the numerical computation, including the curves’ initial arrangement.

3.3 Discretization of time and construction of structure-preserving schemes

We present the $n$ -step time by $t^{(n)}=n\Delta t$ , where $\Delta t$ represents the time step size. The quantity at time $t^{(n)}$ is described by the superscript $n$ . In this subsection, (5) and (9) are discretized temporally to preserve the structure. According to the concept of DVDM, the variational derivative in the discrete system is defined as follows:

Definition 3.1.

If the time-difference of the energy $E^{(n)}$ in the discrete system can be transformed into the following form,

\frac{E^{(n+1)}-E^{(n)}}{\Delta t}=\sum_{i=1}^{N}\bm{\delta}_{d}E_{i}^{(n+1)}\cdot\frac{\bm{X}_{i}^{(n+1)}-\bm{X}_{i}^{(n)}}{\Delta t}{\hat{r}_{i}}^{(n+1/2)},\quad\hat{r}_{i}^{\left(n+1/2\right)}=\frac{\hat{r}_{i}^{(n+1)}+\hat{r}_{i}^{(n)}}{2},

(25)

then, we call $\bm{\delta}_{d}E_{i}^{(n+1)}$ discrete variational derivative of the energy $E^{(n)}$ .

3.3.1 The Willmore flow

We discretize the Wilmore flow (5) with tangential velocity as follows:

\begin{split}&\frac{\bm{X}_{i}^{(n+1)}-\bm{X}^{(n)}_{i}}{\Delta t}=-\bm{\delta}_{d}B_{i}^{(n+1)}+w_{i}^{(n)}\bm{T}_{i}^{(n)}+\gamma^{(n+1)}\bm{\delta}_{d}B_{i}^{(n+1)},\\ &\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\>\>\>\>\>\>i=1,\dots,N,\quad n=0,1,\dots,\end{split}

(26)

where $\bm{\delta}_{d}B_{i}^{(n+1)}$ represents the bending energy’s discrete variational derivative, which is obtained according to definition 3.1. Detailed derivation and the concrete form of $\bm{\delta}_{d}B_{i}^{(n+1)}$ is presented in the Appendix. In (LABEL:sc), a modification term $\gamma^{(n+1)}\bm{\delta}_{d}B_{i}^{(n+1)}$ is added to preserve the bending energy’s discrete dissipation. This is because, in a continuous system, the variations are parallel to the normal vector $\bm{N}(u,t)$ and orthogonal to the tangent vector $\bm{T}(u,t)$ , but the discrete variation $\bm{\delta}_{d}B_{i}^{(n+1)}$ and the tangent vector $\bm{T}_{i}$ obtained from the formula (17) are not necessarily orthogonal. Thus, the modification term $\gamma^{(n+1)}\bm{\delta}_{d}B_{i}^{(n+1)}$ is added to make the scheme structure-preserving.

The Lagrange multiplier $\gamma^{(n+1)}$ is then determined so that the scheme conserves the equation’s energy structure.

For simplicity, we represent the discrete inner product as

(\bm{X},\bm{Y})_{d,n}=\sum_{i=1}^{N}\bm{X}_{i}\cdot\bm{Y}_{i}\hat{r}_{i}^{\left(n+1/2\right)},

(27)

where

\bm{X}=(\bm{X}_{1},\bm{X}_{2},\dots,\bm{X}_{N}),\quad\bm{Y}=(\bm{Y}_{1},\bm{Y}_{2},\dots,\bm{Y}_{N})\in\mathbb{R}^{2\times N}

(28)

represent polygonal curves with vertices $\bm{X}_{i}$ and $\bm{Y}_{i}$ . By the definition of the discrete variational derivative, we obtain

\begin{split}&\frac{B_{d}^{(n+1)}-B_{d}^{(n)}}{\Delta t}=\left(\bm{\delta}_{d}B^{(n+1)},\frac{\bm{X}^{(n+1)}-\bm{X}^{(n)}}{{\Delta t}}\right)_{d,n}\\ &=\left(\bm{\delta}_{d}B^{(n+1)},-\bm{\delta}_{d}B^{(n+1)}+w^{(n)}\bm{T}^{(n)}+\gamma^{(n+1)}\bm{\delta}_{d}B^{(n+1)}\right)_{d,n}\\ &=-\left(\bm{\delta}_{d}B^{(n+1)},\bm{\delta}_{d}B^{(n+1)}\right)_{d,n}+\left(\bm{\delta}_{d}B^{(n+1)},w^{(n)}\bm{T}^{(n)}+\gamma^{(n+1)}\bm{\delta}_{d}B^{(n+1)}\right)_{d,n}.\end{split}

(29)

Therefore, determining $\gamma^{(n+1)}$ by

\left(\bm{\delta}_{d}B^{(n+1)},w^{(n)}\bm{T}^{(n)}+\gamma^{(n+1)}\bm{\delta}_{d}B^{(n+1)}\right)_{d,n}=0,

(30)

we obtain

\frac{B^{(n+1)}-B^{(n)}}{\Delta t}=-\left(\bm{\delta}_{d}B^{(n+1)},\bm{\delta}_{d}B^{(n+1)}\right)_{d,n}\leq 0.

(31)

The dissipative nature of the discrete bending energy is verified from (31). The resulting $\gamma^{(n+1)}$ is anticipated to be quite small. This is because the directions of the discrete variational derivative and the tangent vector are considered almost orthogonal. Thus, it is reasonable to introduce a new variable $\gamma$ in the above way.

From the above, we propose a structure-preserving scheme for the Wilmore flow as follows: for $n=0,1,\dots,$

\begin{split}&\frac{\bm{X}_{i}^{(n+1)}-\bm{X}^{(n)}_{i}}{\Delta t}=-\bm{\delta}_{d}B_{i}^{(n+1)}+w_{i}^{(n)}\bm{T}_{i}^{(n)}-\gamma^{(n+1)}\bm{\delta}_{d}B_{i}^{(n+1)},\quad i=1,2,...N,\\ &\left(\bm{\delta}_{d}B^{(n+1)},w^{(n)}\bm{T}^{(n)}+\gamma^{(n+1)}\bm{\delta}_{d}B^{(n+1)}\right)_{d,n}=0.\end{split}

(32)

Note that the above scheme are nonlinear equations for the unknowns $\bm{X}_{i}^{(n+1)}(i=1,\dots,N)$ and $\gamma^{(n+1)}$ .

3.3.2 The Helfrich flow

The Helfrich flow (9) is discretized into the following form, and we introduce tangential velocities:

\begin{split}&\frac{\bm{X}_{i}^{(n+1)}-\bm{X}_{i}^{(n)}}{\Delta t}\\ &=-\bm{\delta}_{d}B_{i}^{(n+1)}+\lambda^{(n+1)}\bm{\delta}_{d}L_{i}^{(n+1)}+\mu^{(n+1)}\bm{\delta}_{d}A_{i}^{(n+1)}+w_{i}^{(n)}\bm{T}_{i}^{(n)}+\gamma^{(n+1)}\>\bm{\delta}_{d}B_{i}^{(n+1)},\\ &\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\>\>\>\>\>\>i=1,\dots,N,\quad n=0,1,\dots,\end{split}

(33)

where $\bm{\delta}_{d}B_{i}^{(n+1)}$ , $\bm{\delta}_{d}L_{i}^{(n+1)}$ , and $\bm{\delta}_{d}A_{i}^{(n+1)}$ are, respectively, the discrete variational derivatives of the bending energy $B_{d}$ , the length $L_{d}$ , and the enclosed area $A_{d}$ . The detailed deformations and their respective concrete forms are presented in the Appendix. Tangential velocities and their modification terms are also added, as in the Willmore flow’s scheme.

The Lagrange multipliers $\lambda^{(n+1)},\mu^{(n+1)}$ , and $\gamma^{(n+1)}$ are determined so that the scheme reproduces the equation’s energy structure. We first determine the Lagrange multipliers $\lambda^{(n+1)}$ and $\mu^{(n+1)}$ for the conservation of the length $L_{d}$ and the enclosed area $A_{d}$ . By the definition of the discrete variational derivatives, the time-differences of the length and the area are expressed as

		$\displaystyle\frac{L_{d}^{(n+1)}-L_{d}^{(n)}}{\Delta t}=\left(\bm{\delta}_{d}L^{(n+1)},\frac{\bm{X}^{(n+1)}-\bm{X}^{(n)}}{\Delta t}\right)_{d,n},$		(34)
		$\displaystyle\frac{A_{d}^{(n+1)}-A_{d}^{(n)}}{\Delta t}=\left(\bm{\delta}_{d}A^{(n+1)},\frac{\bm{X}^{(n+1)}-\bm{X}^{(n)}}{\Delta t}\right)_{d,n}.$		(35)

By substituting (LABEL:hel_d) for time-difference of vertices $\bm{X}$ in (34) and (35), we have

\frac{L_{d}^{(n+1)}-L_{d}^{(n)}}{\Delta t}\\ =(\bm{\delta}_{d}L^{(n+1)},-\bm{\delta}_{d}B_{i}^{(n+1)}+\lambda^{(n+1)}\bm{\delta}_{d}L^{(n+1)}+\mu^{(n+1)}\bm{\delta}_{d}A^{(n+1)}+w^{(n)}\bm{T}^{(n)}+\gamma^{(n+1)}\>\bm{\delta}_{d}B^{(n+1)})_{d,n},

(36)

\frac{A_{d}^{(n+1)}-A_{d}^{(n)}}{\Delta t}\\ =(\bm{\delta}_{d}A^{(n+1)},-\bm{\delta}_{d}B_{i}^{(n+1)}+\lambda^{(n+1)}\bm{\delta}_{d}L^{(n+1)}+\mu^{(n+1)}\bm{\delta}_{d}A^{(n+1)}+w^{(n)}\bm{T}^{(n)}+\gamma^{(n+1)}\>\bm{\delta}_{d}B^{(n+1)})_{d,n}.

(37)

We can reproduce the conservation of the length and the area by letting (36) and (37) equal to 0.

Next, the Lagrange multiplier $\gamma^{(n+1)}$ for the bending energy’s dissipative nature is determined. By substituting (LABEL:hel_d) for the time-difference of vertices $\bm{X}$ and transforming, we obtain

\begin{split}&\frac{B_{d}^{(n+1)}-B_{d}^{(n)}}{\Delta t}=\left(\bm{\delta}_{d}B^{(n+1)},\frac{\bm{X}^{(n+1)}-\bm{X}^{(n)}}{{\Delta t}}\right)_{d,n}\\ =&\left(\bm{\delta}_{d}B^{(n+1)},-\bm{\delta}_{d}B_{i}^{(n+1)}+\lambda^{(n+1)}\bm{\delta}_{d}L_{i}^{(n+1)}+\mu^{(n+1)}\bm{\delta}_{d}A_{i}^{(n+1)}+w_{i}^{(n)}\bm{T}_{i}^{(n)}+\gamma^{(n+1)}\>\bm{\delta}_{d}B_{i}^{(n+1)}\right)_{d,n}\\ =&\left(\bm{\delta}_{d}B^{(n+1)},-\bm{\delta}_{d}B_{i}^{(n+1)}+\lambda^{(n+1)}\bm{\delta}_{d}L_{i}^{(n+1)}+\mu^{(n+1)}\bm{\delta}_{d}A_{i}^{(n+1)}+w_{i}^{(n)}\bm{T}_{i}^{(n)}+\gamma^{(n+1)}\>\bm{\delta}_{d}B_{i}^{(n+1)}\right)_{d,n}\\ &+\>\mathrm{det}\left(\begin{matrix}(\bm{\delta}B_{d},\bm{\delta}B_{d})&(\bm{\delta}_{d}B,\bm{\delta}_{d}L)&(\bm{\delta}_{d}B,\bm{\delta}_{d}A)\\ (\bm{\delta}_{d}L,\bm{\delta}_{d}B)&(\bm{\delta}_{d}L,\bm{\delta}_{d}L)&(\bm{\delta}_{d}L,\bm{\delta}_{d}A)\\ (\bm{\delta}_{d}A,\bm{\delta}_{d}B)&(\bm{\delta}_{d}A,\bm{\delta}_{d}L)&(\bm{\delta}_{d}A,\bm{\delta}_{d}A)\end{matrix}\right)/\mathrm{det}\left(\begin{matrix}(\bm{\delta}_{d}L,\bm{\delta}_{d}L)&(\bm{\delta}_{d}L,\bm{\delta}_{d}A)\\ (\bm{\delta}_{d}A,\bm{\delta}_{d}L)&(\bm{\delta}_{d}A,\bm{\delta}_{d}A)\end{matrix}\right)\\ &-\>\mathrm{det}\left(\begin{matrix}(\bm{\delta}B_{d},\bm{\delta}B_{d})&(\bm{\delta}_{d}B,\bm{\delta}_{d}L)&(\bm{\delta}_{d}B,\bm{\delta}_{d}A)\\ (\bm{\delta}_{d}L,\bm{\delta}_{d}B)&(\bm{\delta}_{d}L,\bm{\delta}_{d}L)&(\bm{\delta}_{d}L,\bm{\delta}_{d}A)\\ (\bm{\delta}_{d}A,\bm{\delta}_{d}B)&(\bm{\delta}_{d}A,\bm{\delta}_{d}L)&(\bm{\delta}_{d}A,\bm{\delta}_{d}A)\end{matrix}\right)/\mathrm{det}\left(\begin{matrix}(\bm{\delta}_{d}L,\bm{\delta}_{d}L)&(\bm{\delta}_{d}L,\bm{\delta}_{d}A)\\ (\bm{\delta}_{d}A,\bm{\delta}_{d}L)&(\bm{\delta}_{d}A,\bm{\delta}_{d}A)\end{matrix}\right).\end{split}

(38)

Here, for simplicity, the following substitutions are made.

D(\bm{\delta}_{d}B,\bm{\delta}_{d}L,\bm{\delta}_{d}A)=\mathrm{det}\left(\begin{matrix}(\bm{\delta}B_{d},\bm{\delta}B_{d})&(\bm{\delta}_{d}B,\bm{\delta}_{d}L)&(\bm{\delta}_{d}B,\bm{\delta}_{d}A)\\ (\bm{\delta}_{d}L,\bm{\delta}_{d}B)&(\bm{\delta}_{d}L,\bm{\delta}_{d}L)&(\bm{\delta}_{d}L,\bm{\delta}_{d}A)\\ (\bm{\delta}_{d}A,\bm{\delta}_{d}B)&(\bm{\delta}_{d}A,\bm{\delta}_{d}L)&(\bm{\delta}_{d}A,\bm{\delta}_{d}A)\end{matrix}\right)/\mathrm{det}\left(\begin{matrix}(\bm{\delta}_{d}L,\bm{\delta}_{d}L)&(\bm{\delta}_{d}L,\bm{\delta}_{d}A)\\ (\bm{\delta}_{d}A,\bm{\delta}_{d}L)&(\bm{\delta}_{d}A,\bm{\delta}_{d}A)\end{matrix}\right).

(39)

Determining the Lagrange multiplier $\gamma^{(n+1)}$ from

\begin{split}&(\bm{\delta}_{d}B^{(n+1)},-\bm{\delta}_{d}B^{(n+1)}+\lambda^{(n+1)}\bm{\delta}_{d}L^{(n+1)}+\mu^{(n+1)}\bm{\delta}_{d}A^{(n+1)}+w^{(n)}\bm{T}^{(n)}+\gamma^{(n+1)}\>\bm{\delta}_{d}B^{(n+1)})_{d,n}\\ &+D(\bm{\delta}_{d}B,\bm{\delta}_{d}L,\bm{\delta}_{d}A)=0,\end{split}

(40)

we obtain the discrete bending energy dissipation. From the above, we propose a structure-preserving scheme for the Helfrich flow as follows: for $n=0,1,\dots,$

\begin{split}&\frac{\bm{X}_{i}^{(n+1)}-\bm{X}_{i}^{(n)}}{\Delta t}=-\bm{\delta}_{d}B_{i}^{(n+1)}+\lambda^{(n+1)}\bm{\delta}_{d}L_{i}^{(n+1)}+\mu^{(n+1)}\bm{\delta}_{d}A_{i}^{(n+1)}+w_{i}^{(n)}\bm{T}_{i}^{(n)}+\gamma^{(n+1)}\>\bm{\delta}_{d}B_{i}^{(n+1)},\\ &i=1,2,...N,\\ &(\bm{\delta}_{d}L^{(n+1)},-\bm{\delta}_{d}B^{(n+1)}+\lambda^{(n+1)}\bm{\delta}_{d}L^{(n+1)}+\mu^{(n+1)}\bm{\delta}_{d}A^{(n+1)}+w^{(n)}\bm{T}^{(n)}+\gamma^{(n+1)}\>\bm{\delta}_{d}B^{(n+1)})_{d,n}=0,\\ &(\bm{\delta}_{d}A^{(n+1)},-\bm{\delta}_{d}B^{(n+1)}+\lambda^{(n+1)}\bm{\delta}_{d}L^{(n+1)}+\mu^{(n+1)}\bm{\delta}_{d}A^{(n+1)}+w^{(n)}\bm{T}^{(n)}+\gamma^{(n+1)}\>\bm{\delta}_{d}B^{(n+1)})_{d,n}=0,\\ &(\bm{\delta}_{d}B^{(n+1)},-\bm{\delta}_{d}B^{(n+1)}+\lambda^{(n+1)}\bm{\delta}_{d}L^{(n+1)}+\mu^{(n+1)}\bm{\delta}_{d}A^{(n+1)}+w^{(n)}\bm{T}^{(n)}+\gamma^{(n+1)}\>\bm{\delta}_{d}B^{(n+1)})_{d,n}\\ &+D(\bm{\delta}_{d}B,\bm{\delta}_{d}L,\bm{\delta}_{d}A)=0.\end{split}

(41)

Note that the above scheme are nonlinear equations for the unknowns $\bm{X}_{i}^{(n+1)}(i=1,\dots,N)$ , $\lambda^{(n+1)},\mu^{(n+1)}$ , and $\gamma^{(n+1)}$ .

Remark 3.2.

By determining $\gamma^{(n+1)}$ in (LABEL:helf_gamma), the difference quotient of the bending energy can be expressed as

\begin{split}&\frac{B_{d}^{(n+1)}-B_{d}^{(n)}}{\Delta t}\\ &=-\>\mathrm{det}\left(\begin{matrix}(\bm{\delta}B_{d},\bm{\delta}B_{d})_{n,d}&(\bm{\delta}_{d}B,\bm{\delta}_{d}L)_{n,d}&(\bm{\delta}_{d}B,\bm{\delta}_{d}A)_{n,d}\\ (\bm{\delta}_{d}L,\bm{\delta}_{d}B)_{n,d}&(\bm{\delta}_{d}L,\bm{\delta}_{d}L)_{n,d}&(\bm{\delta}_{d}L,\bm{\delta}_{d}A)_{n,d}\\ (\bm{\delta}_{d}A,\bm{\delta}_{d}B)_{n,d}&(\bm{\delta}_{d}A,\bm{\delta}_{d}L)_{n,d}&(\bm{\delta}_{d}A,\bm{\delta}_{d}A)_{n,d}\end{matrix}\right)/\mathrm{det}\left(\begin{matrix}(\bm{\delta}_{d}L,\bm{\delta}_{d}L)_{n,d}&(\bm{\delta}_{d}L,\bm{\delta}_{d}A)_{n,d}\\ (\bm{\delta}_{d}A,\bm{\delta}_{d}L)_{n,d}&(\bm{\delta}_{d}A,\bm{\delta}_{d}A)_{n,d}\end{matrix}\right),\end{split}

(42)

which corresponds to (15). Here, superscripts $(n+1)$ for $\bm{\delta}B_{d}$ , $\bm{\delta}L_{d}$ , and $\bm{\delta}A_{d}$ are omitted. Therefore, the derivation of dissipativity in the continuous system can be reproduced in the discrete system by determining the $\gamma$ by (LABEL:helf_gamma). Moreover, explicitly obtaining $\gamma^{(n+1)}$ from the expression(LABEL:hell_ski), we have

\gamma^{(n+1)}\mathrm{det}\left(\begin{matrix}(\bm{\delta}_{d}B,\bm{\delta}_{d}B)_{n,d}&(\bm{\delta}_{d}B,\bm{\delta}_{d}L)_{n,d}&(\bm{\delta}_{d}B,\bm{\delta}_{d}A)_{n,d}\\ (\bm{\delta}_{d}L,\bm{\delta}_{d}B)_{n,d}&(\bm{\delta}_{d}L,\bm{\delta}_{d}L)_{n,d}&(\bm{\delta}_{d}L,\bm{\delta}_{d}A)_{n,d}\\ (\bm{\delta}_{d}A,\bm{\delta}_{d}B)_{n,d}&(\bm{\delta}_{d}A,\bm{\delta}_{d}L)_{n,d}&(\bm{\delta}_{d}A,\bm{\delta}_{d}A)_{n,d}\end{matrix}\right)=-w^{(n)}\mathrm{det}\left(\begin{matrix}(\bm{\delta}_{d}B,\bm{T})_{n,d}&(\bm{\delta}_{d}B,\bm{T})_{n,d}&(\bm{\delta}_{d}B,\bm{T})_{n,d}\\ (\bm{\delta}_{d}L,\bm{T})_{n,d}&(\bm{\delta}_{d}L,\bm{T})_{n,d}&(\bm{\delta}_{d}L,\bm{T})_{n,d}\\ (\bm{\delta}_{d}A,\bm{T})_{n,d}&(\bm{\delta}_{d}A,\bm{T})_{n,d}&(\bm{\delta}_{d}A,\bm{T})_{n,d}\end{matrix}\right).

(43)

By (43), the value of $\gamma$ is considered to approach zero in the continuous limit.

4 Numerical experiments

In this section, we provide numerical experiments of the schemes proposed in Section 3 for the Willmore and the Helfrich flows. Nonlinear equations (LABEL:will_sk) and (LABEL:hell_ski) are solved by Scipy.optimize.root, which is a common nonlinear solver in Python with a relative residual tolerance of $10^{-6}$ .

4.1 Willmore flow

In this subsection, the findings of numerical computations by the scheme (LABEL:will_sk) are presented.

4.1.1 Examination. 1

Numerical experiments are carried out for different values of the parameter $\alpha$ on tangential velocity to verify the effect of the value of $\alpha$ . The initial curve is set to

\begin{split}x_{1}(t)&=0.5a_{1}(t),\qquad x_{2}(t)=0.54a_{3}(t),\qquad t\in[0,1],\\ a_{1}(t)&=1.8\cos(2\pi t),\qquad a_{2}(t)=0.2+\sin(\pi t)\sin(6\pi t)\sin(2a_{1}(t)),\\ a_{3}(t)&=0.5\sin(2\pi t)+\sin(a_{1}(t))+a_{2}(t)\sin(2\pi t),\end{split}

(44)

which is taken from [20]. We set 30 vertices $\bm{X}_{i}^{0}=(X_{i}^{0},Y_{i}^{0})^{\mathrm{T}}$ to

X_{i}^{0}=x_{1}(i/N),\quad Y_{i}^{0}=x_{2}(i/N),\quad i=1,2,\dots,30.

(45)

Refer to caption — Figure 1: Initial polygonal curve by (45).

The vertices’ arrangement in Fig. 2 is not uniform, and numerical computations are prone to instability. Thus, we adjust them by conducting a numerical computation with $\bm{\delta}_{d}B_{i}\;(i=1,2,...,N)$ equal to $\bm{0}$ in (LABEL:will_sk). The vertices were rearranged under the following conditions

\alpha=5,\quad\Delta t=10^{-4},\quad\text{up to}\>\;T=0.1.

(46)

Fig. 2 shows a initial polygonal curve after relocation. Numerical experiments are carried out with the initial curve illustrated in Fig. 2, time step size $\Delta t=10^{-4}$ and the constant $c_{0}=2$ . We set $\alpha=$ 0,10,50, and 200.

When tangential velocity is absent, i.e., $\alpha=0$ , the concentration of vertices occurs at the step 22. Fig. 3 demonstrates that the vertices are concentrated at $(x,y)=(-0.25,-0.5),(0.8,0.6)$ .

When $\alpha=10$ , the concentration of vertices occurs the step 30. Fig. 4 demonstrates that the vertices are concentrated at $(x,y)=(-0.75,-0.75),(0.8,0.6)$ . This result suggests that the value of $\alpha$ is too small to get stable numerical results.

When $\alpha=50$ , the numerical computation demonstrates that the nonlinear solver converges up to the 3000th step. Fig. 6 demonstrates the initial polygonal curve and the numerical solutions at $n=10$ , 50, 100, 300, 500, 1000, and 3000 steps, and we can see that the solution is approaching a regular circle with a radius $1/c_{0}=0.5$ . Fig. 6 demonstrates the time evolution of the bending energy, and we can see the discrete version of energy dissipative nature. Furthermore, the value of $\gamma$ at each step approached 0 as the time step progressed (Fig. 7).

When $\alpha=200$ , the vertices started to oscillate more tangentially from the fifth step, and the nonlinear solver stopped converging in the seventh step (Fig. 8).

The above results suggest that the value of $\alpha$ is an important parameter in numerical calculations. If the value of $\alpha$ is too small, the concentration of vertices may occur. On the other hand, if the value of $\alpha$ is too big, vertices may oscillate, which causes the concentration of vertices again. Thus, we must choose the value of $\alpha$ appropriately.

4.1.2 Examination. 2

Numerical experiments are conducted with $N=40$ vertices and the rectangular initial curve shown in Fig. 9 . The time step size $\Delta t=10^{-4}$ , the constant $c_{0}=2$ , and the tangential velocity parameter $\alpha=30$ are used.

Fig. 11 demonstrates the initial polygonal curve and the numerical solutions at $n=50$ , 500, 1000, 2500, 5000, and 10000 steps, and we can see that the solution is approaching a regular circle with radius $1/c_{0}=0.5$ . Fig. 11 demonstrates the time evolution of the bending energy, and we can observe the discrete version of the energy dissipative nature.

4.2 Helfrich flow

In this subsection, the findings of numerical computations using the scheme (LABEL:hell_ski) are presented.

4.2.1 Examination. 3

Numerical experiments are conducted for various values of the parameter $\alpha$ on tangential velocity to comfirm the impact of $\alpha$ with the initial curve illustrated in Fig. 2, time step size $\Delta t=10^{-4}$ and the constant $c_{0}=2$ . We set $\alpha=$ 0,10,100, and 200.

When tangential velocity is absent, i.e., $\alpha=0$ , the nonlinear solver stopped converging at the step 15 because of the concentration of the vertices. Fig. 12 demonstrates that the vertices are concentrated at $(x,y)=(-0.4,-0.45)$ .

When $\alpha=10$ , the nonlinear solver stopped converging at the step 16 because of the concentration of vertices. Fig. 13 demonstrates that the vertices are concentrated at $(x,y)=(-0.4,-0.5)$ .

When $\alpha=100$ , the numerical computations demonstrates that the nonlinear solver converges up to the 100th step. Fig. 14 demonstrates the initial polygonal curve and the numerical solutions at $n=10$ , 20, and 100 steps. Fig. 15–16 respectively demonstrates the time evolution of the bending energy, the relative error of the length and the enclosed area for each step. The discrete bending energy’s dissipative nature and the conservation of the length and the enclosed area can be observed. Note that the relative error of the energy $E$ ( $E=A_{d},L_{d}$ ) at the $n$ th step is defined by

\frac{E^{(n)}-E^{(0)}}{E^{(0)}}.

(47)

When $\alpha=200$ , the vertices started to oscillate more tangentially from the 8th step, and the nonlinear solver stopped converging at the 10th step (Fig. 17).

Even in the case of the Helfrich flow, the numerical results suggest that the value of $\alpha$ is a crucial parameter in numerical computations and that it must be selected appropriately.

4.2.2 Examination. 4

Numerical experiments are conducted with $N=40$ vertices and the rectangular initial curve illustrated in Fig. 9. For the time step size $\Delta t=10^{-4}$ , the constant $c_{0}=2$ and the tangential velocity parameter $\alpha=100$ are used. Numerical computations are stopped when the condition

\frac{B^{(n)}-B^{(n+1)}}{B^{(n)}}<\epsilon

(48)

is satisfied. In this examination, we set $\epsilon=10^{-5}$ . Numerical computation stopped at the 2567th step by condition (48). Fig. 19 demonstrates the initial polygonal curve and the numerical solutions at $n=50$ , 100, 300, 500, and 1000 steps. Fig. 19–20 respectively show the time evolution of the bending energy and the relative error of the length and the enclosed area for each step. The discrete bending energy’s dissipative nature and the conservation of the length and the enclosed area can be seen.

4.2.3 Examination. 5

By (43), the value of $\gamma^{(n)}$ obtained by the scheme (LABEL:hell_ski) is considered to approach 0 in the continuous limit. To confirm this numerically, two experiments are conducted: one with a fixed number of vertices and different time step sizes, and the other with a fixed time step size and different number of vertices. For each experiment, we check the value of $\gamma$ when the numerical experiment stops by condition (48). The conditions changing time step size are as follows: the number of vertices $N=50$ , the constant $c_{0}=2$ , the tangential velocity parameter $\alpha=100$ , and the stopping condition $\epsilon=10^{-5}$ . We use five time step sizes $\Delta t=10^{-6}$ , $2.5\times 10^{-6}$ , $5\times 10^{-6}$ , $7.5\times 10^{-6}$ , and $10^{-5}$ . The condition changing the number of vertices is conducted with the time step size $\Delta t=10^{-5}$ , the constant $c_{0}=2$ , the tangential velocity parameter $\alpha=100$ , and the stopping condition $\epsilon=10^{-4}$ . We use five initial polygonals with vertices $N=30$ , 40, 50, 60, and 70. In these examination, we set initial polygonal curves by (44).

In the experiments with various time step sizes, the value of $\gamma$ at the final step became smaller as the time step size decreased (Fig. 22). In the experiments with a different number of vertices, the value of $\gamma$ at the final step also became smaller as the number of vertices increased (Fig. 22). Therefore, in both experiments, the value of $\gamma$ may approach 0 in the continuous limit.

5 Conclusion

In this study, structure-preserving numerical methods for the Willmore and the Helfrich flows are developed. The energy structure was reproduced by discretizing the variational derivatives of the bending energy, the length, and the enclosed area based on DVDM and by determining Lagrange multipliers properly. Moreover, the tangential velocity by Deckelnick was appropriately introduced to prevent the concentration of vertices. Through numerical experiments, we can verify that the Willmore and the Helfrich flows can be computed correctly and our method reproduces the energy structures.

One of the problems of this method is that it needs the solution of nonlinear equations. Therefore, the well-posedness of the scheme is nontrivial and the computational cost is expensive. Thus, theoretical treatment and linearization of the scheme are future tasks. Additionally, the appropriate choice of $\alpha$ for tangential velocity is also a future task.

Acknowledgment

This work was supported by JSPS KAKENHI Grant Numbers 19K14590 and 21K18301, Japan. The authors would like to thank Enago (www.enago.jp) for the English language review.

Appendix A Appendix

A.1 Discrete variational derivative of the bending energy

Let us calculate the discrete variational derivative of the bending energy. For simplicity, we denote $x=x^{(n+1)}$ and $\tilde{x}=x^{(n)}$ for $x=\bm{X}_{i},\;k_{i},\;r_{i},$ etc. The bending energy of the $(n+1)$ -th step can be transformed as

\begin{split}B&=\frac{1}{2}\sum_{i=1}^{N}(k_{i}-c_{0})^{2}\hat{r_{i}}=\frac{1}{2}\sum_{i=1}^{N}(k_{i}-c_{0})^{2}\frac{r_{i}+r_{i+1}}{2}\\ &=\sum_{i=1}^{N}\frac{(k_{i}-c_{0})^{2}+(k_{i-1}-c_{0})^{2}}{4}r_{i}\\ &=\sum_{i=1}^{N}C_{i}r_{i},\end{split}

(49)

where

C_{i}=\frac{(k_{i}-c_{0})^{2}+(k_{i-1}-c_{0})^{2}}{4}.

(50)

From the above form, the energy difference is

\begin{split}B-\tilde{B}=\sum_{i=1}^{N}\left(C_{i}r_{i}-\tilde{C}_{i}\tilde{r}_{i}\right)=\sum_{i=1}^{N}\frac{C_{i}+\tilde{C}_{i}}{2}(r_{i}-\tilde{r}_{i})+\sum_{i=1}^{N}\frac{r_{i}+\tilde{r}_{i}}{2}(C_{i}-\tilde{C}_{i}),\end{split}

(51)

where the formula

ab-cd=\frac{a+c}{2}(b-d)+(a-c)\frac{b+d}{2}

is used. Now let

A_{1}=\sum_{i=1}^{N}\frac{C_{i}+\tilde{C}_{i}}{2}(r_{i}-\tilde{r}_{i}),\quad A_{2}=\sum_{i=1}^{N}\frac{r_{i}+\tilde{r}_{i}}{2}(C_{i}-\tilde{C}_{i}).

Then, we obtain

B-\tilde{B}=A_{1}+A_{2}.

The first term $A_{1}$ can be transformed as

\begin{split}A_{1}&=\sum_{i=1}^{N}\frac{C_{i}+\tilde{C}_{i}}{2}(r_{i}-\tilde{r}_{i})\\ &=\sum_{i=1}^{N}\frac{C_{i}+\tilde{C}_{i}}{2}\frac{|\bm{X}_{i}-\bm{X}_{i-1}|^{2}-|\tilde{\bm{X}}_{i}-\tilde{\bm{X}}_{i-1}|^{2}}{|\bm{X}_{i}-\bm{X}_{i-1}|+|\tilde{\bm{X}}_{i}-\tilde{\bm{X}}_{i-1}|}\\ &=\sum_{i=1}^{N}\frac{C_{i}+\tilde{C}_{i}}{2}\frac{r_{i}\bm{t}_{i}+\tilde{r}_{i}\tilde{\bm{t}}_{i}}{r_{i}+\tilde{r}_{i}}\cdot(\bm{X}_{i}-\tilde{\bm{X}}_{i}-\bm{X}_{i-1}+\tilde{\bm{X}}_{i-1})\\ &=\sum_{i=1}^{N}\left(\frac{C_{i}+\tilde{C}_{i}}{2}\frac{r_{i}\bm{t}_{i}+\tilde{r}_{i}\tilde{\bm{t}}_{i}}{r_{i}+\tilde{r}_{i}}-\frac{C_{i+1}+\tilde{C}_{i+1}}{2}\frac{r_{i+1}\bm{t}_{i+1}+\tilde{r}_{i+1}\tilde{\bm{t}}_{i+1}}{r_{i+1}+\tilde{r}_{i+1}}\right)\cdot(\bm{X}_{i}-\tilde{\bm{X}_{i}})\\ &=\sum_{i=1}^{N}\bm{D}_{i}\cdot(\bm{X}_{i}-\tilde{\bm{X}_{i}}),\end{split}

(52)

where

\bm{D}_{i}=\frac{C_{i}+\tilde{C}_{i}}{2}\frac{r_{i}\bm{t}_{i}+\tilde{r}_{i}\tilde{\bm{t}}_{i}}{r_{i}+\tilde{r}_{i}}-\frac{C_{i+1}+\tilde{C}_{i+1}}{2}\frac{r_{i+1}\bm{t}_{i+1}+\tilde{r}_{i+1}\tilde{\bm{t}}_{i+1}}{r_{i+1}+\tilde{r}_{i+1}}.

Subsequently, $A_{2}$ can be transformed as

\begin{split}A_{2}&=\sum_{i=1}^{N}\frac{r_{i}+\tilde{r}_{i}}{2}(C_{i}-\tilde{C}_{i})\\ &=\sum_{i=1}^{N}\frac{r_{i}+\tilde{r}_{i}}{2}\left\{\frac{(k_{i}-c_{0})^{2}+(k_{i-1}-c_{0})^{2}}{4}-\frac{(\tilde{k}_{i}-c_{0})^{2}+(\tilde{k}_{i-1}-c_{0})^{2}}{4}\right\}\\ &=\sum_{i=1}^{N}\left(\frac{r_{i}+\tilde{r}_{i}}{8}+\frac{r_{i+1}+\tilde{r}_{i+1}}{8}\right)(k_{i}-c_{0}+\tilde{k}_{i}-c_{0})(k_{i}-\tilde{k}_{i})\\ &=\sum_{i=1}^{N}E_{i}(k_{i}-\tilde{k}_{i}),\end{split}

(53)

where

E_{i}=\left(\frac{r_{i}+\tilde{r}_{i}}{8}+\frac{r_{i+1}+\tilde{r}_{i+1}}{8}\right)(k_{i}-c_{0}+\tilde{k}_{i}-c_{0}).

Furthermore, $A_{2}$ can be transformed as

\begin{split}A_{2}&=\sum_{i=1}^{N}E_{i}(k_{i}-\tilde{k}_{i})\\ &=\sum_{i=1}^{N}E_{i}\left\{\frac{\mathrm{det}\left[\delta_{i}^{\left<1\right>}\bm{X}_{i},\delta_{i}^{\left<2\right>}\bm{X}_{i}\right]}{g_{i}^{3}}-\frac{\mathrm{det}\left[\delta_{i}^{\left<1\right>}\tilde{\bm{X}}_{i},\delta_{i}^{\left<2\right>}\tilde{\bm{X}}_{i}\right]}{\tilde{g}_{i}^{3}}\right\}\\ &=\sum_{i=1}^{N}E_{i}\left(\frac{1}{g_{i}^{3}}+\frac{1}{\tilde{g}_{i}^{3}}\right)\frac{1}{2}\left(\mathrm{det}\left[\delta_{i}^{\left<1\right>}\bm{X}_{i},\delta_{i}^{\left<2\right>}\bm{X}_{i}\right]-\mathrm{det}\left[\delta_{i}^{\left<1\right>}\tilde{\bm{X}}_{i},\delta_{i}^{\left<2\right>}\tilde{\bm{X}}_{i}\right]\right)\\ &+\sum_{i=1}^{N}E_{i}\frac{1}{2}\left(\mathrm{det}\left[\delta_{i}^{\left<1\right>}\bm{X}_{i},\delta_{i}^{\left<2\right>}\bm{X}_{i}\right]+\mathrm{det}\left[\delta_{i}^{\left<1\right>}\tilde{\bm{X}}_{i},\delta_{i}^{\left<2\right>}\tilde{\bm{X}}_{i}\right]\right)\left(\frac{1}{g_{i}^{3}}-\frac{1}{\tilde{g}_{i}^{3}}\right)\\ &=\sum_{i=1}^{N}E_{i}F_{i}+\sum_{i=1}^{N}E_{i}G_{i},\end{split}

(54)

where

F_{i}=\left(\frac{1}{g_{i}^{3}}+\frac{1}{\tilde{g}_{i}^{3}}\right)\frac{1}{2}\left(\mathrm{det}\left[\delta_{i}^{\left<1\right>}\bm{X}_{i},\delta_{i}^{\left<2\right>}\bm{X}_{i}\right]-\mathrm{det}\left[\delta_{i}^{\left<1\right>}\tilde{\bm{X}}_{i},\delta_{i}^{\left<2\right>}\tilde{\bm{X}}_{i}\right]\right),

G_{i}=\frac{1}{2}\left(\mathrm{det}\left[\delta_{i}^{\left<1\right>}\bm{X}_{i},\delta_{i}^{\left<2\right>}\bm{X}_{i}\right]+\mathrm{det}\left[\delta_{i}^{\left<1\right>}\tilde{\bm{X}}_{i},\delta_{i}^{\left<2\right>}\tilde{\bm{X}}_{i}\right]\right)\left(\frac{1}{g_{i}^{3}}-\frac{1}{\tilde{g}_{i}^{3}}\right).

Transformations are then performed for $F_{i}$ and $G_{i}$ . Note that the component of $\bm{X}_{i}$ is

\bm{X}_{i}=\left(\begin{matrix}x_{i}\\ y_{i}\end{matrix}\right).

(55)

Firstly, $F_{i}$ can be transformed as

\begin{split}F_{i}&=\frac{1}{2}\left(\frac{1}{g_{i}^{3}}+\frac{1}{\tilde{g}_{i}^{3}}\right)\left(\mathrm{det}\left[\delta_{i}^{\left<1\right>}\bm{X}_{i},\delta_{i}^{\left<2\right>}\bm{X}_{i}\right]-\mathrm{det}\left[\delta_{i}^{\left<1\right>}\tilde{\bm{X}}_{i},\delta_{i}^{\left<2\right>}\tilde{\bm{X}}_{i}\right]\right)\\ &=\frac{1}{2}\left(\frac{1}{g_{i}^{3}}+\frac{1}{\tilde{g}_{i}^{3}}\right)(\delta_{i}^{\left<1\right>}x_{i}\delta_{i}^{\left<2\right>}y_{i}-\delta_{i}^{\left<2\right>}x_{i}\delta_{i}^{\left<1\right>}y_{i}-\delta_{i}^{\left<1\right>}\tilde{x}_{i}\delta_{i}^{\left<2\right>}\tilde{y}_{i}+\delta_{i}^{\left<2\right>}\tilde{x}_{i}\delta_{i}^{\left<1\right>}\tilde{y}_{i})\\ &=\frac{1}{2}\left(\frac{1}{g_{i}^{3}}+\frac{1}{\tilde{g}_{i}^{3}}\right)\left(\begin{matrix}\frac{\delta_{i}^{\left<2\right>}{y}_{i}+\delta_{i}^{\left<2\right>}\tilde{y}_{i}}{2}\\ -\frac{\delta_{i}^{\left<2\right>}{x}_{i}+\delta_{i}^{\left<2\right>}\tilde{x}_{i}}{2}\end{matrix}\right)\cdot\left(\begin{matrix}\delta_{i}^{\left<1\right>}{x}_{i}-\delta_{i}^{\left<1\right>}\tilde{x}_{i}\\ \delta_{i}^{\left<1\right>}{y}_{i}-\delta_{i}^{\left<1\right>}\tilde{y}_{i}\end{matrix}\right)\\ &+\frac{1}{2}\left(\frac{1}{g_{i}^{3}}+\frac{1}{\tilde{g}_{i}^{3}}\right)\left(\begin{matrix}-\frac{\delta_{i}^{\left<1\right>}{y}_{i}+\delta_{i}^{\left<1\right>}\tilde{y}_{i}}{2}\\ \frac{\delta_{i}^{\left<1\right>}{x}_{i}+\delta_{i}^{\left<1\right>}\tilde{x}_{i}}{2}\end{matrix}\right)\cdot\left(\begin{matrix}\delta_{i}^{\left<2\right>}{x}_{i}-\delta_{i}^{\left<2\right>}\tilde{x}_{i}\\ \delta_{i}^{\left<2\right>}{y}_{i}-\delta_{i}^{\left<2\right>}\tilde{y}_{i}\end{matrix}\right)\\ &=\bm{H}_{i}\cdot\delta_{i}^{\left<1\right>}(\bm{X}_{i}-\tilde{\bm{X}}_{i})+\bm{I}_{i}\cdot\delta_{i}^{\left<2\right>}(\bm{X}_{i}-\tilde{\bm{X}}_{i}),\end{split}

(56)

where

\bm{H}_{i}=\frac{1}{2}\left(\frac{1}{g_{i}^{3}}+\frac{1}{\tilde{g}_{i}^{3}}\right)\left(\begin{matrix}\frac{\delta_{i}^{\left<2\right>}{y}_{i}+\delta_{i}^{\left<2\right>}\tilde{y}_{i}}{2}\\ -\frac{\delta_{i}^{\left<2\right>}{x}_{i}+\delta_{i}^{\left<2\right>}\tilde{x}_{i}}{2}\end{matrix}\right),

\bm{I}_{i}=\frac{1}{2}\left(\frac{1}{g_{i}^{3}}+\frac{1}{\tilde{g}_{i}^{3}}\right)\left(\begin{matrix}-\frac{\delta_{i}^{\left<1\right>}{y}_{i}+\delta_{i}^{\left<1\right>}\tilde{y}_{i}}{2}\\ \frac{\delta_{i}^{\left<1\right>}{x}_{i}+\delta_{i}^{\left<1\right>}\tilde{x}_{i}}{2}\end{matrix}\right).

Subsequently, $G_{i}$ can be transformed as

\begin{split}&G_{i}=\frac{1}{2}\left(\mathrm{det}\left[\delta_{i}^{\left<1\right>}\bm{X}_{i},\delta_{i}^{\left<2\right>}\bm{X}_{i}\right]+\mathrm{det}\left[\delta_{i}^{\left<1\right>}\tilde{\bm{X}}_{i},\delta_{i}^{\left<2\right>}\tilde{\bm{X}}_{i}\right]\right)\left(\frac{1}{g_{i}^{3}}-\frac{1}{\tilde{g}_{i}^{3}}\right)\\ &=-\frac{(g_{i}^{2}+g_{i}\tilde{g}_{i}+\tilde{g}_{i}^{2})}{2g_{i}^{3}\tilde{g}_{i}^{3}}\left(\mathrm{det}\left[\delta_{i}^{\left<1\right>}\bm{X}_{i},\delta_{i}^{\left<2\right>}\bm{X}_{i}\right]+\mathrm{det}\left[\delta_{i}^{\left<1\right>}\tilde{\bm{X}}_{i},\delta_{i}^{\left<2\right>}\tilde{\bm{X}}_{i}\right]\right)({g_{i}}-\tilde{g}_{i})\\ &=G_{2,i}({g_{i}}-\tilde{g}_{i}),\end{split}

(57)

where

G_{2,i}=-\frac{(g_{i}^{2}+g_{i}\tilde{g}_{i}+\tilde{g}_{i}^{2})}{2g_{i}^{3}\tilde{g}_{i}^{3}}\left(\mathrm{det}\left[\delta_{i}^{\left<1\right>}\bm{X}_{i},\delta_{i}^{\left<2\right>}\bm{X}_{i}\right]+\mathrm{det}\left[\delta_{i}^{\left<1\right>}\tilde{\bm{X}}_{i},\delta_{i}^{\left<2\right>}\tilde{\bm{X}}_{i}\right]\right).

Furthermore,

\begin{split}G_{i}&=G_{2,i}({g_{i}}-\tilde{g}_{i})\\ &=\frac{G_{2,i}}{g_{i}+\tilde{g}_{i}}(|\delta_{i}^{\left<1\right>}{\bm{X}}_{i}|^{2}-|\delta_{i}^{\left<1\right>}\tilde{\bm{X}}_{i}|^{2})\\ &=\frac{G_{2,i}}{g_{i}+\tilde{g}_{i}}(\delta_{i}^{\left<1\right>}{\bm{X}}_{i}+\delta_{i}^{\left<1\right>}\tilde{\bm{X}}_{i})\cdot(\delta_{i}^{\left<1\right>}{\bm{X}}_{i}-\delta_{i}^{\left<1\right>}\tilde{\bm{X}}_{i})\\ &=\frac{G_{2,i}}{g_{i}+\tilde{g}_{i}}(\delta_{i}^{\left<1\right>}{\bm{X}}_{i}+\delta_{i}^{\left<1\right>}\tilde{\bm{X}}_{i})\cdot\delta_{i}^{\left<1\right>}({\bm{X}}_{i}-\tilde{\bm{X}}_{i})\\ &=\bm{L}_{i}\cdot\delta_{i}^{\left<1\right>}({\bm{X}}_{i}-\tilde{\bm{X}}_{i}),\end{split}

(58)

where

\bm{L}_{i}=\frac{G_{2,i}}{g_{i}+\tilde{g}_{i}}(\delta_{i}^{\left<1\right>}{\bm{X}}_{i}+\delta_{i}^{\left<1\right>}\tilde{\bm{X}}_{i}).

Substituting $F_{i}$ and $G_{i}$ into (54) and using the summation-by-parts formulas, we obtain

\begin{split}A_{2}&=\sum_{i=1}^{N}\left\{E_{i}\bm{H}_{i}\cdot\delta_{i}^{\left<1\right>}(\bm{X}_{i}-\tilde{\bm{X}}_{i})+E_{i}\bm{I}_{i}\cdot\delta_{i}^{\left<2\right>}(\bm{X}_{i}-\tilde{\bm{X}}_{i})+E_{i}\bm{L}_{i}\cdot\delta_{i}^{\left<1\right>}({\bm{X}}_{i}-\tilde{\bm{X}}_{i})\right\}\\ &=\sum_{i=1}^{N}\left\{-\delta_{i}^{\left<1\right>}(E_{i}\bm{H}_{i})+\delta_{i}^{\left<2\right>}(E_{i}\bm{I}_{i})-\delta_{i}^{\left<1\right>}(E_{i}\bm{L}_{i})\right\}\cdot(\bm{X}_{i}-\tilde{\bm{X}}_{i}).\end{split}

(59)

By transforming the equation so far, the difference of the bending energy can be transformed as

\begin{split}&\frac{B^{(n+1)}-\tilde{B}^{(n)}}{\Delta t}\\ &=\sum_{i=1}^{N}\frac{\bm{D}^{(n+1)}_{i}-\delta_{i}^{\left<1\right>}(E^{(n+1)}_{i}\bm{H}^{(n+1)}_{i})+\delta_{i}^{\left<2\right>}(E^{(n+1)}_{i}\bm{I}^{(n+1)}_{i})-\delta_{i}^{\left<1\right>}(E^{(n+1)}_{i}\bm{L}^{(n+1)}_{i})}{\hat{r}^{(n+1/2)}_{i}}\cdot\frac{\bm{X}^{(n+1)}_{i}-\tilde{\bm{X}^{(n)}}_{i}}{\Delta t}\hat{r}^{(n+1/2)}_{i}.\end{split}

(60)

We can define $\bm{\delta}_{d}B$ by

\bm{\delta}_{d}B_{i}^{(n+1)}=\frac{\bm{D}^{(n+1)}_{i}-\delta_{i}^{\left<1\right>}(E^{(n+1)}_{i}\bm{H}^{(n+1)}_{i})+\delta_{i}^{\left<2\right>}(E^{(n+1)}_{i}\bm{I}^{(n+1)}_{i})-\delta_{i}^{\left<1\right>}(E^{(n+1)}_{i}\bm{L}^{(n+1)}_{i})}{\hat{r}^{(n+1/2)}_{i}}.

(61)

A.2 Discrete variational derivative of the length

By calculating the time difference of the discrete length as well as the bending energy, we obtain

\begin{split}&\frac{L^{(n+1)}-L^{(n)}}{\Delta t}=\sum_{i=1}^{N}\left(r_{i}^{(n+1)}-r_{i}^{(n)}\right)\\ &=\sum_{i=1}^{N}\frac{|\bm{X}^{(n+1)}_{i}-\bm{X}^{(n+1)}_{i-1}|^{2}-|\bm{X}_{i}^{(n)}-\bm{X}^{(n)}_{i-1}|^{2}}{|\bm{X}^{(n+1)}_{i}-\bm{X}^{(n+1)}_{i-1}|-|\bm{X}^{(n)}_{i}-\bm{X}^{(n)}_{i-1}|}\\ &=\sum_{i=1}^{N}\frac{r_{i}^{(n+1)}\bm{t}^{(n+1)}_{i}+r^{(n)}_{i}\bm{t}_{i}^{(n)}}{r_{i}^{(n+1)}+r_{i}^{(n)}}\cdot(\bm{X}^{(n+1)}_{i}-\bm{X}^{(n)}_{i}-\bm{X}^{(n+1)}_{i-1}+\bm{X}^{(n)}_{i-1})\\ &=\sum_{i=1}^{N}\left(\frac{r_{i}^{(n+1)}\bm{t}^{(n+1)}_{i}+r^{(n)}_{i}\bm{t}_{i}^{(n)}}{r_{i}^{(n+1)}+r_{i}^{(n)}}-\frac{r_{i+1}^{(n+1)}\bm{t}^{(n+1)}_{i+1}+r^{(n)}_{i+1}\bm{t}_{i+1}^{(n)}}{r_{i+1}^{(n+1)}+r_{i+1}^{(n)}}\right)\cdot(\bm{X}^{(n+1)}_{i}-\bm{X}_{i}^{(n)})\\ &=\sum_{i=1}^{N}\frac{1}{\hat{r}_{i}^{(n+1/2)}}\left(\frac{r_{i}^{(n+1)}\bm{t}^{(n+1)}_{i}+r^{(n)}_{i}\bm{t}_{i}^{(n)}}{r_{i}^{(n+1)}+r_{i}^{(n)}}-\frac{r_{i+1}^{(n+1)}\bm{t}^{(n+1)}_{i+1}+r^{(n)}_{i+1}\bm{t}_{i+1}^{(n)}}{r_{i+1}^{(n+1)}+r_{i+1}^{(n)}}\right)\cdot(\bm{X}^{(n+1)}_{i}-\bm{X}_{i}^{(n)})\hat{r}_{i}^{(n+1/2)}.\end{split}

(62)

Thus, the discrete variational derivative of the length is defined by

\bm{\delta}_{d}L_{i}^{(n+1)}=\frac{1}{\hat{r}_{i}^{(n+1/2)}}\left(\frac{r_{i}^{(n+1)}\bm{t}^{(n+1)}_{i}+r^{(n)}_{i}\bm{t}_{i}^{(n)}}{r_{i}^{(n+1)}+r_{i}^{(n)}}-\frac{r_{i+1}^{(n+1)}\bm{t}^{(n+1)}_{i+1}+r^{(n)}_{i+1}\bm{t}_{i+1}^{(n)}}{r_{i+1}^{(n+1)}+r_{i+1}^{(n)}}\right).

(63)

A.3 Discrete variational derivative of the area

Finally, we calculate the discrete variation of the area. The energy change can be transformed as

\begin{split}&\frac{A^{(n+1)}-A^{(n)}}{\Delta t}=\frac{1}{2}\sum_{i=1}^{N}\left({\bm{X}^{(n+1)}_{i-1}}^{\perp}\cdot\bm{X}^{(n+1)}_{i}-{\bm{X}^{(n)}_{i-1}}^{\perp}\cdot\bm{X}^{(n)}_{i}\right)\\ &=\frac{1}{2}\sum_{i=1}^{N}\left(x^{(n+1)}_{i-1}y^{(n+1)}_{i}-x^{(n+1)}_{i}y^{(n+1)}_{i-1}-x^{(n)}_{i-1}y^{(n)}_{i}+x^{(n)}_{i}y^{(n)}_{i-1}\right)\\ &=\sum_{i=1}^{N}\left(\begin{matrix}\left(-y^{(n+1)}_{i-1}-y^{(n)}_{i-1}+y^{(n+1)}_{i+1}+y^{(n)}_{i+1}\right)/4\\ \left(x^{(n+1)}_{i-1}+x^{(n)}_{i-1}-x^{(n+1)}_{i+1}-x^{(n)}_{i+1}\right)/4\end{matrix}\right)\cdot(\bm{X}^{(n+1)}_{i}-\bm{X}_{i}^{(n)})\\ &=\sum_{i=1}^{N}\frac{1}{\hat{r}_{i}^{(n+1/2)}}\left(\begin{matrix}\left(-y^{(n+1)}_{i-1}-y^{(n)}_{i-1}+y^{(n+1)}_{i+1}+y^{(n)}_{i+1}\right)/4\\ \left(x^{(n+1)}_{i-1}+x^{(n)}_{i-1}-x^{(n+1)}_{i+1}-x^{(n)}_{i+1}\right)/4\end{matrix}\right)\cdot(\bm{X}^{(n+1)}_{i}-\bm{X}_{i}^{(n)})\hat{r}_{i}^{(n+1/2)},\end{split}

(64)

where $A^{\perp}$ denotes $\begin{pmatrix}-a_{2}\\ a_{1}\end{pmatrix}$ for $\bm{A}=\begin{pmatrix}a_{1}\\ a_{2}\end{pmatrix}$ . Thus, discrete variational derivative of the area is defined by

\bm{\delta}_{d}A_{i}^{(n+1)}=\frac{1}{\hat{r}_{i}^{(n+1/2)}}\left(\begin{matrix}\left(-y^{(n+1)}_{i-1}-y^{(n)}_{i-1}+y^{(n+1)}_{i+1}+y^{(n)}_{i+1}\right)/4\\ \left(x^{(n+1)}_{i-1}+x^{(n)}_{i-1}-x^{(n+1)}_{i+1}-x^{(n)}_{i+1}\right)/4\end{matrix}\right).

(65)

References

[1] C. Truesdell. The influence of elasticity on analysis: The classic heritage. Bull. Amer. Math. Soc. (N.S.), 1983, 12, 293–310.
[2] K. Deckelnick and H. Ch. Grunau. Boundary value problems for the one-dimensional Willmore equation-almost explicit solutions. Calc. Var., 2007, 30, 293–314.
[3] G. Dziuk, E. Kuwert, and R. Schatzle. Evolution of elastic curves in $\mathbb{R}^{n}$ : existence and computation. SIAM J. Math. Anal., 2002, 33, 1228–1245.
[4] W. Helfrich. Elastic properties of lipid bilayers: theory and possible experiments. Z. Naturforschung C., 1973, 28, 693–703.
[5] P. B. Canham. The minimum energy of bending as a possible explanation of the biconcave shape of the human red blood cell. J. Theor. Biol., 1970, 26, 61–81.
[6] T. Kurihara and T. Nagasawa. On the gradient flow for a shape optimization problem of plane curves as a singular limit. Saitama Math. J., 2007, 24, 43–75.
[7] Y. Kohsaka and T. Nagasawa. On the existence for the Helfrich flow and its center manifold near spheres. Differ. Integr. Equations, 2006, 19, 121–142.
[8] J. W. Barrett, H. Garcke, and R. Nürnberg. Parametric approximation of Willmore flow and related geometric evolution equations. SIAM. J. Sci. Comput., 2009, 31, 225–253.
[9] M. Beneš, K. Mikula, T. Oberhuber, and D. Ševčovič. Comparison study for level set and direct lagrangian methods for computing Willmore flow of closed planar curves. Comput. Visual Sci., 2009, 12, 307–317.
[10] K. Mikula and D. Ševčovič. A direct method for solving an anisotropic mean curvature flow of planar curve with an external force. Math. Methods Appl. Sci., 2004, 27, 1545–1565.
[11] K. Mikula and D. Ševčovič. Computational and qualitative aspects of evolution of curves driven by curvature and external force. Comput. Visual. Sci., 2004, 6, 211–225.
[12] T. Kemmochi. Energy dissipative numerical schemes for gradient flows of planar curves. BIT Numer. Math., 2017, 57, 991–1017.
[13] L. L. Schumaker. Spline functions: basic theory. Cambridge University Press., third edition, 2007.
[14] T. Matsuo. Dissipative/conservative galerkin method using discrete partial derivative for nonlinear evolution equations. J. Comput. Appl. Math., 2008, 218, 506–521.
[15] D. Furihata and T. Matsuo. Discrete variational derivative method. A structure-preserving numerical method for partial differential equations. New York: CRC Press, Boca Raton, FL, 2011.
[16] D. Furihata. Finite difference schemes for $\frac{\partial u}{\partial t}=\left(\frac{\partial u}{\partial t}\right)^{\alpha}\frac{\delta{{G}}}{\delta u}$ that inherit energy conservation or dissipation property. CJ. Comput. Phys., 1999, 156, 181–205.
[17] K. Deckelnick. Weak solutions of the curve shortening flow. Calc. Var. Partial Differential Equations., 1997, 5, 489–510.
[18] D. Ševčovič. and S. Yazaki. On a gradient flow of plane curves minimizing the anisoperimetric ratio. IAENG Int. J. Appl. Math., 2013, 43, 160–171.
[19] C. L. Epstein and M. Gage. The curve shortening flow, in: Wave motion: Theory, modelling, and computation, Berkeley, Calif. Math. Sci. Res. Inst., 1987, 7, 15–59.
[20] K. Sakakibara and Y. Miyatake. A fully discrete curve-shortening polygonal evolution law for moving boundary problems. J. Comput. Phys., 2021, 424, 109857.

A structure-preserving numerical method for the fourth-order geometric evolution equations for planar curves

Abstract

keywords:

1 Introduction

2 Properties of the target equations

2.1 Willmore flow

2.2 Helfrich flow

3 Construction of structure-preserving numerical method for each equation

3.1 Discretization of the curve

3.2 Discrete tangential velocity

3.3 Discretization of time and construction of structure-preserving schemes

Definition 3.1.

3.3.1 The Willmore flow

3.3.2 The Helfrich flow

Remark 3.2.

4 Numerical experiments

4.1 Willmore flow

4.1.1 Examination. 1

4.1.2 Examination. 2

4.2 Helfrich flow

4.2.1 Examination. 3

4.2.2 Examination. 4

4.2.3 Examination. 5

5 Conclusion

Acknowledgment

Appendix A Appendix

A.1 Discrete variational derivative of the bending energy

A.2 Discrete variational derivative of the length

A.3 Discrete variational derivative of the area

References

A structure-preserving numerical method
for the fourth-order geometric evolution equations for planar curves

3 Construction of structure-preserving
numerical method for each equation