Stable discretizations of elastic flow in Riemannian manifolds

# Stable discretizations of elastic flow in Riemannian manifolds

## Abstract

The elastic flow, which is the -gradient flow of the elastic energy, has several applications in geometry and elasticity theory. We present stable discretizations for the elastic flow in two-dimensional Riemannian manifolds that are conformally flat, i.e. conformally equivalent to the Euclidean space. Examples include the hyperbolic plane, the hyperbolic disk, the elliptic plane as well as any conformal parameterization of a two-dimensional manifold in , . Numerical results show the robustness of the method, as well as quadratic convergence with respect to the space discretization.

E

lastic flow, hyperbolic plane, hyperbolic disk, elliptic plane, Riemannian manifolds, geodesic elastic flow, finite element approximation, stability, equidistribution

{AMS}

65M60, 53C44, 53A30, 35K55

## 1 Introduction

Elastic flow of curves in a two-dimensional Riemannian manifold is given as the -gradient flow of the elastic energy , where is the geodesic curvature. It has been shown, see [6] for the general case and [11] for the hyperbolic plane, that the gradient flow of the elastic energy is given as

 Vg=−(ϰg)sgsg−12ϰ3g−S0ϰg, (1)

where is the normal velocity of the curve with respect to the metric , , denoting arclength, and is the sectional curvature of . The evolution law (1) decreases the curvature energy , and long term limits are expected to be critical points of this energy. These critical points are called free elasticae, see [16], and are of interest in geometry and mechanics. In particular, let us mention that a curve is an absolute minimizer if and only if it is a geodesic. Recently the flow (1) was studied in [11, 12], for the case of the hyperbolic plane, relying on earlier results in [14] for a flat background metric. The hyperbolic plane is a particular case of a manifold with non-positive sectional curvature, which is of particular interest as the set of free elasticae is much richer, see [16].

In this paper, we allow for a general conformally flat metric. Examples include the hyperbolic plane, the hyperbolic disk, the elliptic plane, as well as any conformal parameterization of a two-dimensional manifold in , . For parameterized hypersurfaces in , earlier authors, see e.g. [9, 17, 18, 2, 4], used the surrounding space in their numerical approximations, which leads to errors in directions normal to the hypersurface. This will be avoided by the intrinsic approach used in this paper. In particular, our numerical method leads to approximate solutions which remain on the hypersurface after application of the parameterization map. In addition, in this paper we will present a first numerical analysis for elastic flow in manifolds not embedded in . This in particular makes it possible to compute elastic flow of curves in the hyperbolic plane in a stable way.

For finite element approximations of (1) introduced in [6] it does not appear possible to prove a stability result. It is the aim of this paper to introduce novel approximations for (1) that can be shown to be stable. In particular, we will show that the semidiscrete continuous-in-time approximations admit a gradient flow structure. For relevant literature on conformal metrics we refer to [19, 15]. Curvature driven flows in hyperbolic spaces have been studied by [10, 1, 11, 12, 6], and related numerical approximations of elastic flow of curves can be found in [14, 13, 5, 8] for the Euclidean case, and in [9, 17, 18, 2, 4] for the case of curves on hypersurfaces in .

The outline of this paper is as follows. After formulating the problem in detail in the next section, we will derive in Section 3 weak formulations which will be the basis for our finite element approximation. In Section 4 we introduce continuous-in-time, discrete-in-space discretizations which are based on the weak formulations. For these semidiscrete formulations a stability result will be shown, which is the main contribution of this work. In Section 5 we then formulate fully discrete variants for which we show existence and uniqueness. In Section 6 we present several numerical computations which show convergence rates as well as the robustness of the approach. Finally, in the appendix we show the consistency of the weak formulations presented in Section 3.

## 2 Mathematical formulations

Let be the periodic interval . Let be a parameterization of a closed curve . On assuming that on , we introduce the arclength of the curve, i.e. , and set

 →τ=→xsand→ν=−→τ⊥, (1)

where denotes a clockwise rotation by . For the curvature of it holds that

 ϰ→ν=→ϰ=→τs=→xss=1|→xρ|[→xρ|→xρ|]ρ. (2)

Let be an open set with metric tensor

 [(→v,→w)g](→z)=g(→z)→v.→w∀ →v,→w∈R2 for →z∈H, (3)

where is the standard Euclidean inner product, and where is a smooth positive weight function. The length induced by (3) is defined as

 [|→v|g](→z)=([(→v,→v)g](→z))12=g12(→z)|→v|∀ →v∈R2 for →z∈H. (4)

For , we define the generalized elastic energy as

 Wg,λ(→x)=12∫I(ϰ2g+2λ)|→xρ|gdρ, (5)

where

 ϰg=g−12(→x)[ϰ−12→ν.∇lng(→x)] (6)

is the curvature of the curve with respect to the metric , see [6] for details. Generalized elastic flow is defined as the –gradient flow of (5), and it was established in [6] that a strong formulation is given by

 Vg=g12(→x)→xt.→ν=−(ϰg)sgsg−12ϰ3g−S0ϰg+λϰg, (7)

where and

 S0=−Δlng2g (8)

is the sectional curvature of . We refer to [6] for further details.

The two weak formulations of (7), for , introduced in [6] are based on the equivalent equation

 g(→x)→xt.→ν=−1|→xρ|⎛⎜⎝[ϰg]ρg12(→x)|→xρ|⎞⎟⎠ρ−12g12(→x)ϰ3g−g12(→x)S0(→x)ϰg.

The first uses as a variable, while the second uses as a variable.
: Let . For find and such that \cref@addtoresetequationparentequation

 ∫Ig(→x)→xt.→νχ|→xρ|dρ=∫Ig−12(→x)(g−12(→x)[ϰ−12→ν.∇lng(→x)])ρχρ|→xρ|−1dρ −12∫Ig−1(→x)[ϰ−12→ν.∇lng(→x)]3χ|→xρ|dρ −∫IS0(→x)[ϰ−12→ν.∇lng(→x)]χ|→xρ|dρ∀ χ∈H1(I), (9a) ∫Iϰ→ν.→η|→xρ|dρ+∫I(→xρ.→ηρ)|→xρ|−1dρ=0∀ →η∈[H1(I)]2. (9b)

: Let . For find and such that \cref@addtoresetequationparentequation

 ∫Ig(→x)→xt.→νχ|→xρ|dρ=∫Ig−12(→x)[ϰg]ρχρ|→xρ|−1dρ−12∫Ig12(→x)ϰ3gχ|→xρ|dρ −∫IS0(→x)g12(→x)ϰgχ|→xρ|dρ∀ χ∈H1(I), (10a) ∫Ig(→x)ϰg→ν.→η|→xρ|dρ+∫I[∇g12(→x).→η+g12(→x)→xρ.→ηρ|→xρ|2]|→xρ|dρ=0∀ →η∈[H1(I)]2. (10b)

For the numerical approximations based on and it does not appear possible to prove stability results that show that discrete analogues of (5), for , decrease monotonically in time. Based on the techniques in [5], it is possible to introduce alternative weak formulations, for which semidiscrete continuous-in-time finite element approximations admit such a stability result.

We end this section with some example metrics that are of particular interest in differential geometry. Two families of metrics are given by \cref@addtoresetequationparentequation

 g(→z)=(→z.→e2)−2μ, μ∈R,withH=H2={→z∈R2:→z.→e2>0}, (11a) and g(→z)=4(1−α|→z|2)2,withH={Dα={→z∈R2:|→z|<α−12}α>0,R2α≤0. (11b) The case (11a) with μ=1 models the hyperbolic plane, while μ=0 corresponds to the Euclidean case. The case (11b) with α=1 gives a model for the hyperbolic disk, while α=−1 models the geometry of the elliptic plane. Of course, α=0 collapses to the Euclidean case. Further metrics of interest are induced by conformal parameterizations →Φ:H→Rd, d≥3, of the two-dimensional Riemannian manifold M⊂Rd, i.e. M=→Φ(H) and |∂→e1→Φ(→z)|2=|∂→e2→Φ(→z)|2 and ∂→e1→Φ(→z).∂→e2→Φ(→z)=0 for all →z∈H. Here examples include the Mercator projection of the unit sphere without the north and the south pole, →Φ(→z)=cosh−1(→z.→e1)(cos(→z.→e2),sin(→z.→e2),sinh(→z.→e1))T, so that g(→z)=cosh−2(→z.→e1),withH=R2, (11c) as well as the catenoid parameterization →Φ(→z)=(cosh(→z.→e1)cos(→z.→e2),cosh(→z.→e1) sin(→z.→e2),→z.→e1)T, so that g(→z)=cosh2(→z.→e1),withH=R2. (11d) Based on [20, p. 593] we also recall the following conformal parameterization of a torus with large radius R>1 and small radius r=1 from [6]. In particular, we let s=[R2−1]12 and define , so that g(→z)=s2([s2+1]12−cos(→z.→e2))−2,withH=R2. (11e)

## 3 Weak formulations

We define the first variation of a quantity depending in a differentiable way on , in the direction as

 [δδ→xA(→x)](→χ)=limε→0A(→x+ε→χ)−A(→x)ε, (1)

and observe that, for sufficiently smooth,

 [δδ→xA(→x)](→xt)=ddtA(→x). (2)

For later use, on noting (1), (1) and (4), we observe that \cref@addtoresetequationparentequation

 [δδ→xgβ(→x)](→χ) =βgβ−1(→x)→χ.∇g(→x)=βgβ(→x)→χ.∇lng(→x)∀ β∈R, (3a) [δδ→x∇lng(→x)](→χ) =(D2lng(→x))→χ, (3b) [δδ→x|→xρ|](→χ) =→xρ.→χρ|→xρ|=→τ.→χρ=→τ.→χs|→xρ|, (3c) [δδ→x|→xρ|g](→χ) =(→τ.→χs+12→χ.∇lng(→x))|→xρ|g, (3d) [δδ→x→τ](→χ) =[δδ→x→xρ|→xρ|](→χ)=→χρ|→xρ|−→xρ|→xρ|2→xρ.→χρ|→xρ|=→χs−→τ(→χs.→τ) =(→χs.→ν)→ν, (3e) [δδ→x→ν](→χ) =−[δδ→x→τ⊥](→χ)=−(→χs.→ν)→ν⊥=−(→χs.→ν)→τ, (3f) [δδ→x→ν|→xρ|](→χ) =−[δδ→x→x⊥ρ](→χ)=−→χ⊥ρ=−→χ⊥s|→xρ|, (3g)

where we always assume that is sufficiently smooth so that all the quantities are defined almost everywhere. E.g.  for (3a), (3b), and for (3c)–(3g). In addition, on recalling (1), we have for all that \cref@addtoresetequationparentequation

 →a.→b⊥=−→a⊥.→b, (4a) →a⊥=(→a⊥.→τ)→τ+(→a⊥.→ν)→ν=(→a⊥.→ν⊥)→τ−(→a⊥.→τ⊥)→ν=(→a.→ν)→τ−(→a.→τ)→ν. (4b)

Let denote the –inner product on . In the following we will discuss the –gradient flow of the energy

 Wg,λ(→x)=(12ϰ2g+λ,|→xρ|g)=(12g−12(→x)(ϰ−12→ν.∇lng(→x))2+λg12(→x),|→xρ|), (5)

treating either or formally as an independent variable that has to satisfy the side constraint (9b), or (10b), respectively. For the weak formulations of the –gradient flow obtained in this way, it can be shown that they are consistent with the strong formulation (7), see the appendix. Moreover, we will formally establish that solutions to these weak formulations are indeed solutions to the –gradient flow of (5). Mimicking these stability proofs on the discrete level will yield the main results of this paper.

### 3.1 Based on ϰ

We define the Lagrangian

 L(→x,ϰ⋆,→y) =12(g−12(→x)(ϰ⋆−12→ν.∇lng(→x))2+2λg12(→x),|→xρ|) −(ϰ⋆→ν,→y|→xρ|)−(→xs,→ys|→xρ|), (6)

which is obtained on combining (5) and the side constraint

 (ϰ⋆→ν,→η|→xρ|)+(→xs,→ηs|→xρ|)=0∀ →η∈[H1(I)]2, (7)

recall (9b) and (1). Taking variations in , and setting we obtain (7). Combining (7) and (9b) yields, on recalling (1), that , and we are going to use this identity from now. Taking variations in and setting leads to

 (g−12(→x)(ϰ−12→ν.∇lng(→x))−→y.→ν,χ|→xρ|)=0∀ χ∈L2(I), (8)

which implies that

 →y.→ν=g−12(→x)(ϰ−12→ν.∇lng(→x))⟺ϰ=g12(→x)→y.→ν+12→ν.∇lng(→x). (9)

Taking variations in , and then setting , where we have noted (7) and (4), yields, on recalling (1), that

 (g32(→x)→xt.→ν,→χ.→ν|→xρ|)=−12((ϰ−12→ν.∇lng(→x))2,[δδ→xg−12(→x)|→xρ|](→χ)) +12(g−12(→x)(ϰ−12→ν.∇lng(→x)),[δδ→x→ν.∇lng(→x)](→χ)|→xρ|) +(ϰ→y,[δδ→x→ν|→xρ|](→χ))+(→yρ,[δδ→x→τ](→χ))−λ(1,[δδ→xg12(→x)|→xρ|](→χ)), (10)

for all . On choosing in (10) we obtain, on noting (2), that

 (g32(→x)(→xt.→ν)2,|→xρ|)=−12((ϰ−12→ν.∇lng(→x))2,[g−12(→x)|→xρ|]t) +12(g−12(→x)(ϰ−12→ν.∇lng(→x)),[→ν.∇lng(→x)]t|→xρ|) +(ϰ→y,[→ν|→xρ|]t)+(→yρ,→τt)−λ(1,[g12(→x)|→xρ|]t). (11)

Differentiating (7) with respect to time, and then choosing yields, on recalling that , that

 (ϰt,→y.→ν|→xρ|)+(ϰ→y,(→ν|→xρ|)t)+(→τt,→yρ)=0. (12)

Combining (11), (12) and (9) gives, on noting (5), that

 (g32(→x)(→xt.→ν)2,|→xρ|)=−12((ϰ−12→ν.∇lng(→x))2,[g−12(→x)|→xρ|]t) +12(g−12(→x)(ϰ−12→ν.∇lng(→x)),[→ν.∇lng(→x)]t|→xρ|) −(ϰt,g−12(→x)(ϰ−12→ν.∇lng(→x))|→xρ|)−λ(1,[g12(→x)|→xρ|]t) =−ddtWg,λ(→x). (13)

The above yields the gradient flow property of the new weak formulation, on noting from (7) and (4) that the left hand side of (3.1) can be equivalently written as .

In order to derive a suitable weak formulation, we now return to (10). Combining (10), (3) and (4a) yields that

 (g32(→x)→xt.→ν,→χ.→ν|→xρ|)=−12(g−12(→x)(ϰ−12→ν.∇lng(→x))2+2λg12(→x),→χs.→τ|→xρ|) +14(g−12(→x)(ϰ−12→ν.∇lng(→x))2−2λg12(→x),→χ.(∇lng(→x))|→xρ|) +12(g−12(→x)(ϰ−12→ν.∇lng(→x))→ν,(D2lng(→x))→χ|→xρ|) −12(g−12(→x)(ϰ−12→ν.∇lng(→x))[lng(→x)]s,→ν.→χs|→xρ|)+(→ys.→ν,→χs.→ν|→xρ|) +(ϰ→y⊥,→χs|→xρ|)∀ →χ∈[H1(I)]2. (14)

Overall we obtain the following weak formulation.
: Let . For find and such that (14), (8) and

 (ϰ→ν,→η|→xρ|)+(→xs,→ηs|→xρ|)=0∀ →η∈[H1(I)]2 (15)

hold. We remark that in the Euclidean case (8) collapses to , and so on eliminating from (14) and (15), and on noting (4b), we obtain that the formulation collapses to [5, (2.4a,b)] for the Euclidean elastic flow.

### 3.2 Based on ϰg

We recall that was inspired by the formulation , which is based on acting as a variable. In order to derive an alternative formulation, we now start from , where the curvature is a variable.

We begin by equivalently rewriting the side constraint (10b) as

 (g12(→x)ϰg→ν,→η|→xρ|g)+(→xs,→ηs|→xρ|g)+12(∇lng(→x),→η|→xρ|g)=0∀ →η∈[H1(I)]2, (16)

where we have noted (1), (4) and . Combining (5) and (16) leads to the Lagrangian

 Lg(→x,ϰ⋆g,→yg)