Numerical approximation of bi-harmonic wave maps into spheres

Abstract.

We construct a structure preserving non-conforming finite element approximation scheme for the bi-harmonic wave maps into spheres equation. It satisfies a discrete energy law and preserves the non-convex sphere constraint of the continuous problem. The discrete sphere constraint is enforced at the mesh-points via a discrete Lagrange multiplier. This approach restricts the spatial approximation to the (non-conforming) linear finite elements. We show that the numerical approximation converges to the weak solution of the continuous problem in spatial dimension $d=1$. The convergence analysis in dimensions $d>1$ is complicated by the lack of a discrete product rule as well as the low regularity of the numerical approximation in the non-conforming setting. Hence, we show convergence of the numerical approximation in higher-dimensions by introducing additional stabilization terms in the numerical approximation. We present numerical experiments to demonstrate the performance of the proposed numerical approximation and to illustrate the regularizing effect of the bi-Laplacian which prevents the formation of singularities.

1. Introduction

Let $T>0$, $\Omega\subset\mathbb{R}^{d}$ be an open, bounded, (convex) polyhedral domain, $d\leq 3$, and $\Omega_{T}:=(0,T)\times\Omega$. Further, let $\mathbb{S}^{2}\subset\mathbb{R}^{3}$ denote the unit sphere. We consider bi-harmonic wave maps $\mathbf{u}:(0,T)\times\Omega\to\mathbb{S}^{2}$. Formally, these are critical points of the action functional

$$\Phi(\mathbf{u}):=\frac{1}{2}\int_{\Omega_{T}}|\partial_{t}\mathbf{u}|^{2}-|\Delta\mathbf{u}|^{2}\mathrm{d}(t,x)$$

under the sphere target constraint. More precisely, we consider the Cauchy problem defined by the Euler-Lagrange equation, i.e.

(1.1) $$\partial_{t}^{2}\mathbf{u}+\Delta^{2}\mathbf{u} \perp T_{\textbf{u}}(\mathbb{S}^{2})\qquad \Omega_{T},\qquad\qquad$$ in

(1.2) $$\frac{\partial\mathbf{u}}{\partial{\bf n}} =\frac{\partial\Delta\mathbf{u}}{\partial{\bf n}}=0\qquad \partial\Omega,$$ on

(1.3) $$\mathbf{u}(0) =\mathbf{u}_{0},\quad\partial_{t}\mathbf{u}(0)=\mathbf{v}_{0}\quad \Omega,$$ in

where ${\bf n}$ is the outward unit normal vector field to $\partial\Omega$ and $T_{\textbf{u}(t,x)}(\mathbb{S}^{2})$ denotes the tangent space to $\mathbf{u}(t,x)\in\mathbb{S}^{2}$. (1.1) can be equivalently written as

(1.4) $$\partial_{t}^{2}\mathbf{u}+\Delta^{2}\mathbf{u}=\lambda_{w}u,\;\text{with}\;\lambda_{w}=|\Delta\mathbf{u}|^{2}-|\partial_{t}\mathbf{u}|^{2}-\Delta|\nabla\mathbf{u}|-2\mathrm{div}\langle\Delta\mathbf{u},\nabla\mathbf{u}\rangle.$$

There are a number of analytical results on this problem in the setting $\Omega=\mathbb{R}^{d}$ without boundary condition. In [14], weak solutions in the energy space are constructed by a suitable Ginzburg-Landau-type approximation. Furthermore, in [15] local well-posedness for initial data of sufficiently high Sobolev regularity has been proven (for more general targets), including a blow-up criterion. In dimension $d=1,2$, the energy conservation has been used to prove global existence [18].

As far as we are aware the present paper is the first one addressing the numerical approximation of the bi-harmonic wave maps into spheres equation.

There are a number of results on the numerical analysis for (second order) wave maps, for instance in [6, 7, 16, 3] for the sphere target and with more general target manifolds in [2], [10]. The numerical approximation of wave maps shares many features with the numerical approximation of its parabolic counterparts, the harmonic heat flow equation and with the Landau-Lifshitz-Gilbert equation. Concerning the numerical analysis and further references of these well-studied problems we refer to [1, 9, 5].

A conforming finite element discretisation of problems involving (fourth order) bi-harmonic operators requires at least $C^{1}$ regularity of the employed polynomial approximation spaces, cf. [12]. Apart from the fact that the $H^{2}$-conforming finite element methods are rather complicated to implement, cf. e.g. [17], the requirement to satisfy the sphere constraint on the discrete level restricts the spatial discretization to a piecewise linear finite element setting. Non-conforming methods for fourth order problems based on piecewise linear finite elements have been studied in [11, 13].

Here, in order to show convergence of the numerical approximation we require strong convergence of the gradient of the numerical approximation, which is complicated by the fact the the discrete solution is only $H^{1}$-regular. This prohibits the application of the standard Aubin-Lions-Simon compactness result for the compact embedding $H^{2}\subset H^{1}$. We derive a new discrete compactness result for the non-conforming setting.

The paper is organized as follows. Notation and preliminary results are introduced in Section 2. The numerical scheme is introduced in Section 3 and its convergence towards a weak solution for $d=1$ is shown in Section 4. In Section 5 we discuss a stabilized variant of the numerical scheme that is convergent in higher spatial dimensions. In Section 4 we present a variant of the numerical scheme which preserves the discrete energy. Computational experiments are given in Section 7.

2. Notation and preliminaries

We employ the following notation throughout this paper. By $\mathbf{L}^{2}$, $\mathbf{H}^{1}$ we denote the $\mathbb{R}^{3}$-valued function spaces $L^{2}(\Omega,\mathbb{R}^{3})$, $H^{1}(\Omega,\mathbb{R}^{3})$, respectively. The standard norm in Lebesgue space $L^{p}(\Omega)$, $p>0$ is denoted as $\|\cdot\|_{L^{p}}$. For $p=2$ we denote the inner product on $L^{2}(\Omega)$ as $(\cdot,\cdot):=(\cdot,\cdot)_{L^{2}(\Omega)}$ and the corresponding norm by $\|\cdot\|:=\|\cdot\|_{L^{2}}$; the inner product on $L^{2}(\omega)$ for $\omega\subset\Omega$ is denoted as $(\cdot,\cdot)_{\omega}:=(\cdot,\cdot)_{L^{2}(\omega)}$.

Let $V_{h}\equiv V_{h}(\mathcal{T}_{h})\subset H^{1}(\Omega)$ be the lowest order $H^{1}$-conforming finite element space subordinate to a simplicial partition $\mathcal{T}_{h}$ of $\Omega$, and ${\bf V}_{h}:=[V_{h}]^{3}$. By $\mathcal{N}_{h}$, we denote the set of all nodes, i.e., vertices of all elements $T\in\mathcal{T}_{h}$. Furthermore, we consider the nodal basis $\{\varphi_{\bf z};\,{\bf z}\in\mathcal{N}_{h}\}$ of $V_{h}$ consisting of piecewise linear continuous functions that satisfy $\varphi_{{\bf z}_{k}}({\bf z}_{l})=\delta_{kl}$ for ${\bf z}_{k},{\bf z}_{l}\in\mathcal{N}_{h}$. For ${\bf z}\in\mathcal{N}_{h}$ we denote $\omega_{\bf z}=\mathrm{supp}(\varphi_{\bf z})$. For $T\in\mathcal{T}_{h}$ we denote $h_{T}=\mathrm{diam}(T)$ and the mesh size is denoted as $h=\max_{T\in\mathcal{T}}h_{T}$.

We define the nodal interpolation operator $\mathcal{I}_{h}:C(\overline{\Omega})\rightarrow V_{h}$ by ${\mathcal{I}}_{h}w:=\sum_{{\bf z}\in\mathcal{N}_{h}}w({\bf z})\varphi_{{\bf z}}$, for all $w\in C(\overline{\Omega})$. The following interpolation estimate holds for $p\in(1,\infty]$, [12]:

(2.1) $$\|w-{\mathcal{I}}_{h}w\|_{L^{p}}+h\|\nabla(w-{\mathcal{I}}_{h}w)\|_{L^{p}}\leq Ch^{2}\|D^{2}w\|_{L^{p}}\quad\text{for }w\in W^{2,p}(\Omega),$$

where $h=\max_{T\in\mathcal{T}_{h}}\mathrm{diam}(T)$. We also note that the interpolation operator is $W^{1,\infty}$-stable, cf. e.g. [4, Example 3.7], i.e.: $\|\mathcal{I}_{h}w\|_{W^{1,\infty}}\leq C\|w\|_{W^{1,\infty}}$ for $w\in W^{1,\infty}$.

We denote the inner product on $\mathbb{R}^{3}$ by $\langle\cdot,\cdot\rangle$ and define the discrete version of the $L^{2}$-inner product on ${\bf V}_{h}$ as

$$(\mathbf{v},\mathbf{w})_{h}:=\int_{\Omega}{\mathcal{I}}_{h}\langle\mathbf{v}(x),\mathbf{w}(x)\rangle\,{\rm d}{\bf x}=\sum_{{\bf z}\in\mathcal{N}_{h}}\beta_{{\bf z}}\langle\mathbf{v}({\bf z}),\mathbf{w}({\bf z})\rangle\qquad\forall\,\mathbf{v},\mathbf{w}\in[C(\overline{\Omega})]^{3},$$

where $\beta_{{\bf z}}=\frac{|\omega_{\bf z}|}{d+1}=\int_{\Omega}\varphi_{{\bf z}}\,{\rm d}{\bf x}$, for all ${\bf z}\in\mathcal{N}_{h}$.

It is well known the discrete inner product satisfies

(2.2) $$|(\mathbf{v},\mathbf{w})_{h}-(\mathbf{v},\mathbf{w})| \leq Ch\|\mathbf{v}\|\|\nabla\mathbf{w}\|\,,$$

and the corresponding norm $\|\mathbf{v}\|_{h}^{2}:=(\mathbf{v},\mathbf{v})_{h}$ is equivalent to the $L^{2}$-norm, i.e.,

(2.3) $$\|\mathbf{v}\|_{h}^{2}\leq\|\mathbf{v}\|^{2} \leq(d+2)\|\mathbf{v}\|_{h}^{2}\,,$$

for $\mathbf{v},\mathbf{w}\in{\bf V}_{h}$.

We partition the time interval $[0,T]$ into equidistant subintervals $(t_{n-1},t_{n})$, $n=1,\dots,N$, $t_{n}=n\tau$ with a time-step size $\tau=\frac{T}{N}>0$. For a sequence $\{\psi^{n}\}_{n\geq 0}$ we denote $d_{t}\psi^{n}:=\tau^{-1}\{\psi^{n}-\psi^{n-1}\}$, and $\psi^{n+1/2}:=\frac{1}{2}\{\psi^{n+1}+\psi^{n}\}$ with $n\geq 1$. By $C,\tilde{C}>0$ we denote generic bounded constants which may depend on $\Omega$ and $T$ but are independent of $\tau$, $h$.

For $u\in H^{1}(\Omega)$ we define the discrete Laplacian $\Delta_{h}:H^{1}(\Omega)\rightarrow V_{h}$ as

(2.4) $$-(\Delta_{h}u,\varphi_{h})_{h}=(\nabla u,\nabla\varphi_{h})\qquad\forall\varphi_{h}\in V_{h}\,.$$

The discrete Laplacian satisfies an inverse estimate

(2.5) $$\|\Delta_{h}\phi_{h}\|_{h}^{2}\leq\frac{1}{h^{2}}\|\nabla\phi_{h}\|^{2}\qquad\forall\phi_{h}\in V_{h}.$$

To show convergence of the numerical apprpximation we require the following conditions to be satisfied.

Assumption 1.

We assume that the partition $\mathcal{T}_{h}$ satisfies:

• $\mathcal{T}_{h}$ is quasiuniform and consists of right-angled simplices, i.e. that, one of the following holds: each simplex $T\in\mathcal{T}_{h}$ has one vertex such that all edges intersect at this vertex at right angles (Type-$1$ triangulation, $d=2,3$), each simplex $T\in\mathcal{T}_{h}$ has two vertices at which two edges intersect at right angles (Type-$2$ triangulation, $d=3$).

• $\mathcal{T}_{h}$ is consistent with the discrete Laplace operator (2.4) in the sense that that for all $\boldsymbol{\Phi}\in C^{\infty}_{0}\bigl{(}[0,T);C^{\infty}_{0}(\Omega,\mathbb{R}^{3})\bigr{)}$, the discrete Laplacian $\Delta_{h}:C(\overline{\Omega})\rightarrow V_{h}(\mathcal{T}_{h})$ converges strongly for $h\rightarrow 0$, i.e.,

  $$\|\Delta_{h}\mathcal{I}_{h}\boldsymbol{\Phi}-\Delta\boldsymbol{\Phi}\|_{h}\rightarrow 0\quad\mathrm{for}\quad h\rightarrow 0.$$

Remark 2.1.

The conditions of Assumption 1 can be verified explicitly for certain uniform partitions $\mathcal{T}_{h}$ of $\Omega$ on which the discrete Laplacian agrees with the finite difference approximation of the Laplace operator. For $d=2$ the triangulation consisting of a union of halved squares with side $h$ (Type-$1$ triangulation) agrees with the $5$-point finite difference approximation and for $d=3$ a partition that consists of cubes with side $h$ divided into six tetrahedra (Type-$2$ triangulation) agrees with the $7$-point finite difference approximation. In these two cases the standard theory of finite difference approximation implies that that

$$\Delta_{h}\phi({\bf z})-\Delta\phi({\bf z})=\mathcal{O}(h^{2})\qquad\forall\phi\in C^{4}(\overline{\Omega})$$

for all ${\bf z}\in\mathcal{N}_{h}$. We note that on general meshes only weak convergence of the discrete Laplace operator can be expected, cf. [13].

2.1. Weak solution

Definition 2.1.

Let $(\,\mathbf{u}_{0},{\bf v}_{0}\,)\in H^{2}(\Omega,\mathbb{R}^{3})\times L^{2}(\Omega,\mathbb{R}^{3})$ be such that $|\mathbf{u}_{0}|=1$, and $\langle\mathbf{u}_{0},{\bf v}_{0}\rangle=0$ a.e. in $\Omega\subset\mathbb{R}^{d}$. We call a function $\mathbf{u}:\Omega_{T}\rightarrow\mathbb{R}^{3}$, such that $\partial_{t}\mathbf{u},\Delta\mathbf{u}\in L^{2}(\Omega_{T},\mathbb{R}^{3})$ a weak solution of (1.1), if there holds:

1. $\mathbf{u}(0)=\mathbf{u}_{0}$, and $\partial_{t}\mathbf{u}(0)={\bf v}_{0}$ a.e. in $\Omega$,

2. $|\mathbf{u}|=1$ a.e. in $\Omega_{T}$,

3. for all $\boldsymbol{\Phi}\in C^{\infty}_{0}\bigl{(}[0,T);W^{2,2}(\Omega,\mathbb{R}^{3})\bigr{)}$ there holds

   (2.6) $$-\int_{0}^{T}\bigl{(}\partial_{t}\mathbf{u},\partial_{t}(\mathbf{u}\times\boldsymbol{\Phi})\bigr{)}\,{\rm d}t+\int_{0}^{T}\bigl{(}\Delta\mathbf{u},\Delta(\mathbf{u}\times\boldsymbol{\Phi})\bigr{)}\,{\rm d}t=\bigl{(}{\bf v}_{0}\wedge\mathbf{u}_{0},\boldsymbol{\Phi}(0)\bigr{)}.$$

Recently, (global) existence of weak solutions ${\bf u}:[0,T]\times\Omega\rightarrow\mathbb{S}^{2}$ (for $\Omega\equiv\mathbb{R}^{d}$) has been established in [14]. Note that for $\mathbf{u},\boldsymbol{\Phi}$ as in Definition 2.1 we have $\bigl{(}\partial_{t}\mathbf{u},\partial_{t}(\mathbf{u}\times\boldsymbol{\Phi})\bigr{)}=\bigl{(}\partial_{t}\mathbf{u},\mathbf{u}\times\partial_{t}\boldsymbol{\Phi}\bigr{)}$ and that $|\nabla u|^{2}\leq|\Delta u|$ a.e..

We remark that our definition of a weak solution differs slightly from [14, (1.4)] in the sense that in Definition 2.1 we restrict to a (dense) class of more regular test functions, but this is inessential. Similar to [14, Lemma 2.1] we observe that any sufficiently regular mapping $\mathbf{u}:\Omega_{T}\rightarrow\mathbb{R}^{m}$ satisfying

$$\int_{\Omega_{T}}\Bigl{[}-\langle\partial_{t}\mathbf{u},\partial_{t}\boldsymbol{\Phi}\rangle+\langle\Delta\mathbf{u},\Delta\boldsymbol{\Phi}\rangle-\big{(}|\Delta\mathbf{u}|^{2}-|\partial_{t}\mathbf{u}|^{2}\big{)}\langle\mathbf{u},\boldsymbol{\Phi}\rangle+|\nabla\mathbf{u}|^{2}\Delta\langle\mathbf{u},\boldsymbol{\Phi}\rangle-2\langle\Delta\mathbf{u},\nabla\mathbf{u}\rangle\nabla\langle\mathbf{u},\boldsymbol{\Phi}\rangle\Bigr{]}\,{\rm d}{\bf x}\,{\rm d}{t}=0\,,$$

for all $\boldsymbol{\Phi}\in C^{\infty}_{0}\bigl{(}(0,T);W^{2,2}(\Omega,\mathbb{R}^{3})\bigr{)}$ is a weak solution to (1.1).

3. Numerical approximation

Let $\mathbf{U}^{0}$, $\mathbf{V}^{0}$, s.t. $\langle\mathbf{U}^{0}({\bf z}),\mathbf{V}^{0}({\bf z})\rangle=0$ and $|{\bf U}({\bf z})|=1$ for ${\bf z}\in\mathcal{N}_{h}$, be given and set $\mathbf{U}^{-1}=\mathbf{U}^{0}-\tau\mathbf{V}^{0}$. The numerical approximation of (1.1) is then defined as follows: for $n=0,\dots,N-1$ determine $\mathbf{U}^{n+1},\mathbf{W}^{n+1}\in{\bf V}_{h}$ such that for all $\boldsymbol{\Phi},\boldsymbol{\Psi}\in{\bf V}_{h}$ there holds

(3.1) $$\bigl{(}d_{t}^{2}\mathbf{U}^{n+1},\boldsymbol{\Phi}\bigr{)}_{h}+(\nabla\mathbf{W}^{n+1/2},\nabla\boldsymbol{\Phi}) =\bigl{(}\lambda_{w}^{n+1}\mathbf{U}^{n+1/2},\boldsymbol{\Phi}\bigr{)}_{h}\,,$$

$$(\mathbf{W}^{n+1/2},\boldsymbol{\Psi})_{h} =(\nabla\mathbf{U}^{n+1/2},\nabla\boldsymbol{\Psi})\,,$$

where the discrete Lagrange multiplier $\lambda_{w}^{n+1}\in V_{h}$ is defined for ${\bf z}\in\mathcal{N}_{h}$ as

(3.5) $$\lambda_{w}^{n+1}({\bf z})=\left\{\begin{array}[]{ll}0&\mbox{if }\mathbf{U}^{n+1/2}({\bf z})={\bf 0}\,,\\ \frac{\frac{1}{2}{\langle\mathbf{V}^{0}({\bf z}),\mathbf{V}^{1}({\bf z})\rangle}}{|\mathbf{U}^{1/2}({\bf z})|^{2}}+\frac{\bigl{(}\nabla\mathbf{W}^{1/2},\mathbf{U}^{1/2}({\bf z})\otimes\nabla\varphi_{{\bf z}}\bigr{)}}{\beta_{{\bf z}}|\mathbf{U}^{1/2}({\bf z})|^{2}}&\text{if }n=0,\,\mathbf{U}^{1/2}({\bf z})\neq{\bf 0},\\ -\frac{\bigl{\langle}d_{t}\mathbf{U}^{n}({\bf z}),{d_{t}}\mathbf{U}^{n+1/2}({\bf z})\bigr{\rangle}}{|\mathbf{U}^{n+1/2}({\bf z})|^{2}}+\frac{\bigl{(}\nabla\mathbf{W}^{n+1/2},\mathbf{U}^{n+1/2}({\bf z})\otimes\nabla\varphi_{{\bf z}}\bigr{)}}{\beta_{{\bf z}}|\mathbf{U}^{n+1/2}({\bf z})|^{2}}&\mbox{else}\,.\end{array}\right.$$

The above explicit formula for the discrete Lagrange multiplier $\lambda_{w}^{n+1}$ guarantees that $|{\bf U}^{n+1}({\bf z})|^{2}=|{\bf U}^{n}({\bf z})|^{2}$, see Lemma 3.1 below. The formula for $n>1$ can be deduced by setting $\boldsymbol{\Phi}=\mathbf{U}^{n+1/2}({\bf z})\varphi_{{\bf z}}$ in (3.1) from the identity

$$\frac{1}{\tau}\Bigl{\langle}d_{t}{\bf U}^{n+1}({\bf z})-d_{t}{\bf U}^{n}({\bf z}),{\bf U}^{n+1/2}({\bf z})\Bigr{\rangle}$$

$$\qquad=\frac{1}{2\tau^{2}}\Bigl{[}|{\bf U}^{n+1}({\bf z})|^{2}-|{\bf U}^{n}({\bf z})|^{2}\Bigr{]}-\frac{1}{\tau}\Bigl{\langle}d_{t}{\bf U}^{n}({\bf z}),{\bf U}^{n-1/2}({\bf z})+\tau d_{t}{\bf U}^{n+1/2}({\bf z})\Bigr{\rangle}$$

$$\qquad=-\frac{1}{2\tau^{2}}\Bigl{[}|{\bf U}^{n}({\bf z})|^{2}-|{\bf U}^{n-1}({\bf z})|^{2}\Bigr{]}-\Bigl{\langle}d_{t}{\bf U}^{n}({\bf z}),d_{t}{\bf U}^{n+1/2}({\bf z})\Bigr{\rangle}$$

$$\qquad=-\Bigl{\langle}d_{t}{\bf U}^{n}({\bf z}),d_{t}{\bf U}^{n+1/2}({\bf z})\Bigr{\rangle}\,,$$

where we used that $|{\bf U}^{n+1}({\bf z})|^{2}=|{\bf U}^{n}({\bf z})|^{2}=|{\bf U}^{n-1}({\bf z})|^{2}=1$. The case $n=0$ is considered separately since $|{\bf U}^{-1}({\bf z})|\neq 1$ if $|{\bf V}_{0}({\bf z})|>0$ and the formula can be derived analogously using the orthogonality $\langle\mathbf{U}^{0},\mathbf{V}^{0}\rangle=0$.

For $n\geq 1$ we denote $\mathbf{V}^{n+1}:=d_{t}\mathbf{U}^{n+1}$ and for $n=0,\dots,N$ we define the discrete energy as

$$\mathcal{E}_{h}\bigl{(}{\bf V}^{n+1},{\bf U}^{n+1}\bigr{)}:=\frac{1}{2}\|\mathbf{V}^{n+1}\|_{h}^{2}+\frac{1}{2}\|\Delta_{h}\mathbf{U}^{n+1}\|_{h}^{2}\,.$$

Lemma 3.1.

Let $\mathcal{T}_{h}$ be a quasiuniform triangulation of $\Omega\subset\mathbb{R}^{d}$, and $(\,{\bf U}^{0},{\bf V}^{0}\,)\in{\bf V}_{h}\times{\bf V}_{h}$ such that $|{\bf U}^{0}({\bf z})|=1$, and $\bigl{\langle}{\bf U}^{0}({\bf z}),{\bf V}^{0}({\bf z})\bigr{\rangle}=0$ for all ${\bf z}\in\mathcal{N}_{h}$. For $n\geq 0$, and sufficiently small $\tilde{C}\equiv\tilde{C}(\Omega,\mathcal{T}_{h})>0$ independent of $k,h>0$ such that $\tau\leq\tilde{C}h^{2}$, there exist ${\bf U}^{n+1},\,{\bf W}^{n+1}\in{\bf V}_{h}$, which solve (3.1). Furthermore $|{\bf U}^{n+1}({\bf z})|=1$ for all ${\bf z}\in\mathcal{N}_{h}$, and

(3.6) $$\max_{n=1,\dots,N}\mathcal{E}_{h}\bigl{(}{\bf V}^{n},{\bf U}^{n}\bigr{)}+\frac{\tau^{2}}{2}\sum_{n=1}^{N}\|d_{t}{\bf V}^{n}\|^{2}_{h}=\mathcal{E}_{h}\bigl{(}{\bf V}^{0},{\bf U}^{0}\bigr{)}.$$

Proof.

i) We note that the second equation in (3.1) and (2.4) imply that ${\bf W}^{n+1}=-\Delta_{h}{\bf U}^{n+1}\in{\bf V}_{h}$, $n\geq 0$; the existence of ${\bf U}^{n+1}\in{\bf V}_{h}$ is shown in Lemma 3.2 below.

ii) The property $|{\bf U}^{n+1}({\bf z})|^{2}=|{\bf U}^{n}({\bf z})|^{2}=1$ for all ${\bf z}\in\mathcal{N}_{h}$ follows by choosing $\boldsymbol{\Phi}=\mathbf{U}^{n+1/2}({\bf z})\varphi_{z}$ in (3.1)₁ and noting the definition of the discrete Lagrange multiplier (3.5).

iii) Energy law: we choose $\boldsymbol{\Phi}=d_{t}\mathbf{U}^{n+1}$, $\boldsymbol{\Psi}=d_{t}\Delta_{h}\mathbf{U}^{n+1}$ in (3.1). Noting that

$$\bigl{(}\lambda_{w}^{n+1}\mathbf{U}^{n+1/2},d_{t}\mathbf{U}^{n+1}\bigr{)}_{h}=\bigl{(}\lambda_{w}^{n+1},|\mathbf{U}^{n+1}|^{2}-|\mathbf{U}^{n}|^{2}\bigr{)}_{h}=0,$$

(since $\langle\mathbf{U}^{n+1/2},d_{t}\mathbf{U}^{n+1}\rangle(\mathbf{z})=(|\mathbf{U}^{n+1}|^{2}-|\mathbf{U}^{n}|^{2})(\mathbf{z})=0$ by part i)) we deduce that

(3.7) $$\bigl{(}d_{t}\mathbf{V}^{n+1},\mathbf{V}^{n+1}\bigr{)}_{h}+(\nabla\mathbf{W}^{n+1/2},\nabla d_{t}\mathbf{U}^{n+1}) =0\,,$$

$$(\mathbf{W}^{n+1/2},d_{t}\Delta_{h}\mathbf{U}^{n+1})_{h} =(\nabla\mathbf{U}^{n+1/2},\nabla d_{t}\Delta_{h}\mathbf{U}^{n+1})\,.$$

From the definition of $\Delta_{h}$ it follows that $(\nabla\mathbf{U}^{n+1/2},\nabla d_{t}\Delta_{h}\mathbf{U}^{n+1})=-(\Delta_{h}\mathbf{U}^{n+1/2},d_{t}\Delta_{h}\mathbf{U}^{n+1})_{h}$ and $(\mathbf{W}^{n+1/2},d_{t}\Delta_{h}\mathbf{U}^{n+1})_{h}=-(\nabla\mathbf{W}^{n+1/2},\nabla d_{t}\mathbf{U}^{n+1})$. Hence, we add the respective equations in (3.7) to cancel the mixed terms and obtain after using the identitis $2(a-b)a=a^{2}+(a-b)^{2}-b^{2}$, $(a+b)(a-b)=a^{2}-b^{2}$ that

$$\frac{1}{2}\|\mathbf{V}^{n+1}\|_{h}^{2}+\frac{1}{2}\|\mathbf{V}^{n+1}-\mathbf{V}^{n}\|_{h}^{2}-\frac{1}{2}\|\mathbf{V}^{n}\|_{h}^{2}+\frac{1}{2}\|\Delta_{h}\mathbf{U}^{n+1}\|_{h}^{2}-\frac{1}{2}\|\Delta_{h}\mathbf{U}^{n}\|_{h}^{2}=0\,.$$

The discrete energy law (3.6) then follows directly after summation of the above identity for $n=0,\dots,N-1$.

∎

The solvability of the nonlinear system given by (3.1) is guaranteed by the following lemma.

Lemma 3.2.

Given $\mathbf{U}^{n}$, $\mathbf{U}^{n-1}\in\mathbf{V}_{h}$, for sufficiently small $\tau\leq Ch^{2}$, there exist $\mathbf{U}^{n+1}$, $\mathbf{W}^{n+1}\in\mathbf{V}_{h}$ that satisfy (3.1).

Proof.

$i)$ From (3.1)₂ and (2.4) we observe that $\mathbf{W}^{n}=-\Delta_{h}\mathbf{U}^{n}$, $\mathbf{W}^{n+1}=-\Delta_{h}\mathbf{U}^{n+1}$. Hence, it is enough to show the existence of $\mathbf{U}^{n+1}$.

For $\frac{1}{8}<\varepsilon\leq\frac{1}{4}$ we define the continuous map $\mathbf{F}_{\varepsilon}:\mathbf{V}_{h}\rightarrow\mathbf{V}_{h}$ via

$$(\mathbf{F}_{\varepsilon}(\mathbf{u}),\boldsymbol{\Phi}):=\frac{2}{\tau^{2}}\bigl{(}\mathbf{u}-2\mathbf{U}^{n}+\mathbf{U}^{n-1/2},\boldsymbol{\Phi}\bigr{)}_{h}+\big{(}\nabla\mathbf{w},\nabla\boldsymbol{\Phi}\big{)}-\bigl{(}\lambda_{w,\varepsilon}^{n+1}\mathbf{u},\boldsymbol{\Phi}\bigr{)}_{h},$$

where $\mathbf{w}:=-\Delta_{h}\mathbf{u}$ and

$$\lambda_{w,\varepsilon}^{n+1}({\bf z})=-\frac{\bigl{\langle}d_{t}\mathbf{U}^{n}({\bf z}),{k^{-1}\big{(}\mathbf{u}({\bf z})-\mathbf{U}^{n-1/2}\big{)}}\bigr{\rangle}}{|\mathbf{u}({\bf z})|_{\varepsilon}^{2}}+\frac{\bigl{(}\nabla\mathbf{w},\mathbf{u}({\bf z})\otimes\nabla\varphi_{{\bf z}}\bigr{)}}{\beta_{{\bf z}}|\mathbf{u}({\bf z})|_{\varepsilon}^{2}}:=\lambda_{1}^{n+1}(\mathbf{z})+\lambda_{2}^{n+1}(\mathbf{z})\qquad{\bf z}\in\mathcal{N}_{h}.$$

with $|\cdot|_{\varepsilon}^{2}:=\max\{|\cdot|^{2},\varepsilon\}$.

We estimate the first term above as

$$\bigl{(}\lambda_{1}^{n+1}\mathbf{u},\mathbf{u}\bigr{)}_{h} ={\sum_{{\bf z}\in\mathcal{N}_{h}}}\bigl{(}\lambda_{1}^{n+1}\mathbf{u},\mathbf{u}({\bf z}){\varphi_{{\bf z}}}\bigr{)}_{h}\leq{\sum_{{\bf z}\in\mathcal{N}_{h}}}\frac{|\mathbf{u}({\bf z})|^{2}}{|\mathbf{u}({\bf z})|_{\varepsilon}^{2}}\frac{1}{\tau}\|\mathbf{u}-\mathbf{U}^{n-1/2}\|_{h}\|d_{t}\mathbf{U}^{n}\|_{h}$$

(3.8) $$\leq\frac{1}{{2}\tau^{2}}\|\mathbf{u}-\mathbf{U}^{n-1/2}\|_{h}^{2}+C\|d_{t}\mathbf{U}^{n}\|_{h}^{2}\leq\frac{1}{\tau^{2}}\|\mathbf{u}\|^{2}+\frac{1}{\tau^{2}}\|\mathbf{U}^{n-1/2}\|_{h}^{2}+C\|d_{t}\mathbf{U}^{n}\|_{h}^{2}.$$

Noting that $\|\varphi_{{\bf z}}\|_{L^{\infty}}\leq Ch^{-1}$ we estimate

$$\begin{split}\bigl{(}\lambda_{2}^{n+1}\mathbf{u},\mathbf{u}({\bf z})\varphi_{{\bf z}}\bigr{)}_{h}&\leq\frac{|\mathbf{u}({\bf z})|^{2}}{|\mathbf{u}({\bf z})|_{\varepsilon}^{2}}\bigl{(}\nabla\mathbf{w},\mathbf{u}({\bf z})\otimes\nabla\varphi_{{\bf z}}\bigr{)}_{\omega_{{\bf z}}}\leq C\|\nabla\mathbf{w}\|_{\omega_{{\bf z}}}\|\mathbf{u}({\bf z})\nabla\varphi_{{\bf z}}\|_{\omega_{{\bf z}}}\\ &\leq\frac{C}{h}\|\nabla\mathbf{w}\|_{\omega_{\bf z}}\|\mathbf{u}({\bf z})\|_{\omega_{\bf z}}.\end{split}$$

Consequently, we estimate the second term using the above estimate, the inverse inequality $\|\nabla\mathbf{v}\|\leq Ch^{-1}\|\mathbf{v}\|$ for $\mathbf{v}\in\mathbf{V}_{h}$ and (2.3) as

(3.9) $$\begin{split}\bigl{(}\lambda_{2}^{n+1}\mathbf{u},\mathbf{u}\bigr{)}_{h}&={\sum_{{\bf z}\in\mathcal{N}_{h}}}\bigl{(}\lambda_{2}^{n+1}\mathbf{u},\mathbf{u}({\bf z})\varphi_{{\bf z}}\bigr{)}_{h}\leq\frac{C}{h}\sum_{{\bf z}\in\mathcal{N}_{h}}\|\nabla\mathbf{w}\|_{\omega_{\bf z}}\|\mathbf{u}({\bf z})\|_{\omega_{\bf z}}\\ &\leq\frac{C}{h}\|\nabla\mathbf{w}\|\|\mathbf{u}\|_{h}\leq\frac{C}{h^{2}}\|\mathbf{w}\|\|\mathbf{u}\|_{h}\leq\frac{1}{2}\|\Delta_{h}\mathbf{u}\|_{h}^{2}+\frac{C}{h^{4}}\|\mathbf{u}\|_{h}^{2},\end{split}$$

where we also used the fact that $\sum_{{\bf z}\in\mathcal{N}_{h}}\|\mathbf{v}\|_{\omega_{\bf z}}^{2}\leq C(\mathcal{T}_{h})\|\mathbf{v}\|^{2}$ (where the constant $C(\mathcal{T}_{h})>0$ is $h$-independent due to the regularity of the partition $\mathcal{T}_{h}$).

We use (3), (3.9) along with (2.4) to conclude

$$(\mathbf{F}_{\varepsilon}(\mathbf{u}),\mathbf{u}) \geq\frac{2}{\tau^{2}}\left(\|\mathbf{u}\|_{h}^{2}-2|\bigl{(}\mathbf{U}^{n},\mathbf{u}\bigr{)}_{h}|-|\bigl{(}\mathbf{U}^{n-1/2},\mathbf{u}\bigr{)}_{h}|\right)+\|\Delta_{h}\mathbf{u}\|_{h}^{2}$$

$$\quad-\frac{1}{\tau^{2}}\|\mathbf{u}\|_{h}^{2}-\frac{1}{\tau^{2}}\|\mathbf{U}^{n-1/2}\|_{h}^{2}-C\|d_{t}\mathbf{U}^{n}\|_{h}^{2}-\frac{1}{2}\|\Delta_{h}\mathbf{u}\|_{h}^{2}-\frac{C}{h^{4}}\|\mathbf{u}\|_{h}^{2}$$

$$\geq\frac{1}{\tau^{2}}\left(\left(1-\frac{C\tau^{2}}{h^{4}}\right)\|\mathbf{u}\|_{h}^{2}-{4}\|\mathbf{U}^{n}\|_{h}\|\mathbf{u}\|_{h}-2\|\mathbf{U}^{n-1/2}\|_{h}\|\mathbf{u}\|_{h}\right)$$

$$\quad+\frac{1}{2}\|\Delta_{h}\mathbf{u}\|_{h}^{2}-\frac{1}{2\tau^{2}}\|\mathbf{U}^{n-1/2}\|_{h}^{2}-C\|d_{t}\mathbf{U}^{n}\|_{h}^{2}$$

$$\geq\frac{1}{\tau^{2}}\|\mathbf{u}\|_{h}\left(\left(1-\frac{C\tau^{2}}{h^{4}}\right)\|\mathbf{u}\|_{h}-4\|\mathbf{U}^{n}\|_{h}-2\|\mathbf{U}^{n-1/2}\|_{h}\right)$$

$$\quad+\frac{1}{2}\|\Delta_{h}\mathbf{u}\|_{h}^{2}-\frac{1}{2\tau^{2}}\|\mathbf{U}^{n-1/2}\|_{h}^{2}-C\|d_{t}\mathbf{U}^{n}\|_{h}^{2}.$$

For $\tau\leq\tilde{C}h^{2}$ for sufficiently small $\tilde{C}\equiv\tilde{C}(\Omega)>0$ it is possible to find $R>0$ such that for all $\mathbf{u}\in\mathbf{V}_{h}$ with $\|\mathbf{u}\|_{h}\geq R$ it holds

$$(\mathbf{F}_{\varepsilon}(\mathbf{u}),\mathbf{u})\geq\frac{1}{\tau^{2}}\|\mathbf{u}\|_{h}-\frac{1}{2\tau^{2}}\|\mathbf{U}^{n-1/2}\|_{h}^{2}-C\|d_{t}\mathbf{U}^{n}\|_{h}^{2}>0.$$

Then the Brouwer’s fixed point theorem implies the existence of $\mathbf{u}\equiv\mathbf{U}^{n+1/2}\in\mathbf{V}_{h}$ such that $\mathbf{F}_{\varepsilon}(\mathbf{U}^{n+1/2})=0$.

$ii)$ We show that for sufficiently small $\tau\leq Ch^{2}$ the solution $\mathbf{U}^{n+1/2}$ from step $i)$ also satisfies $\mathbf{F}_{\varepsilon}(\mathbf{U}^{n+1/2})=0$ for $\varepsilon=0$ which implies that $\mathbf{U}^{n+1}=2\mathbf{U}^{n+1/2}-\mathbf{U}^{n}$ is a solution of (3.1). It suffices to show that for $\tau\leq Ch^{2}$ it holds

(3.10) $$|\mathbf{U}^{n+1/2}({\bf z})|\geq 1-\frac{\tau}{2}|\mathrm{d}_{t}\mathbf{U}^{n+1}({\bf z})|>\frac{1}{2}\qquad\forall{\bf z}\in\mathcal{N}_{h}.$$

We proceed by induction. Suppose that for $0\leq j\leq n$ its holds that $\mathbf{U}^{j}\in\mathbf{V}_{h}$ satisfies $|\mathbf{U}^{j}({\bf z})|=1$ for all ${\bf z}\in\mathcal{N}_{h}$ and that

$$\mathcal{E}_{h}\bigl{(}{\bf V}^{n},{\bf U}^{n}\bigr{)}+\frac{\tau^{2}}{2}\sum_{j=1}^{n}\|d_{t}{\bf V}^{j}\|^{2}_{h}=\mathcal{E}_{h}\bigl{(}{\bf V}^{0},{\bf U}^{0}\bigr{)}.$$

The solution $\mathbf{U}^{n+1}$ satisfies

$$\bigl{(}d_{t}^{2}\mathbf{U}^{n+1},\boldsymbol{\Phi}\bigr{)}_{h}+\big{(}\Delta_{h}\mathbf{U}^{n+1/2},\Delta_{h}\boldsymbol{\Phi}\big{)}_{h}=\bigl{(}\lambda_{w,\varepsilon}^{n+1}\mathbf{U}^{n+1/2},\boldsymbol{\Phi}\bigr{)}_{h}=\bigl{(}(\lambda_{1}^{n+1}+\lambda_{2}^{n+1})\mathbf{U}^{n+1/2},\boldsymbol{\Phi}\bigr{)}_{h}.$$

We set $\boldsymbol{\Phi}=d_{t}\mathbf{U}^{n+1}$ in the above and get

$$\frac{1}{2}d_{t}\left[\|d_{t}\mathbf{U}^{n+1}\|_{h}^{2}+\|\Delta_{h}\mathbf{U}^{n+1}\|_{h}^{2}\right]+\frac{\tau}{2}\|d_{t}^{2}\mathbf{U}^{n+1}\|_{h}^{2}=\bigl{(}(\lambda_{1}^{n+1}+\lambda_{2}^{n+1})\mathbf{U}^{n+1/2},d_{t}\mathbf{U}^{n+1}\bigr{)}_{h}.$$

Noting that $\tau d_{t}\mathbf{U}^{n+1}=2\mathbf{U}^{n+1/2}-2\mathbf{U}^{n}$ we estimate

$$(\lambda_{1}^{n+1}\mathbf{U}^{n+1/2},d_{t}\mathbf{U}^{n+1})_{h}\leq\frac{1}{2}\left(\frac{|d_{t}\mathbf{U}^{n}|(|d_{t}\mathbf{U}^{n+1}|+|d_{t}\mathbf{U}^{n}|)}{\max\{|\mathbf{U}^{n+1/2}|,\varepsilon\}},|d_{t}\mathbf{U}^{n+1}|\right)_{h}$$

$$\leq\frac{1}{\tau}\left(\frac{|d_{t}\mathbf{U}^{n}|(|d_{t}\mathbf{U}^{n+1}|+|d_{t}\mathbf{U}^{n}|)}{\max\{|\mathbf{U}^{n+1/2}|,\varepsilon\}},|\mathbf{U}^{n+1/2}|+|\mathbf{U}^{n}|\right)_{h}$$

$$\leq\frac{C}{\tau}\left(1+\frac{\|\mathbf{U}^{n}\|_{\mathbf{L}^{\infty}}}{\varepsilon}\right)\|d_{t}\mathbf{U}^{n}\|_{h}(\|d_{t}\mathbf{U}^{n+1}\|_{h}+\|d_{t}\mathbf{U}^{n}\|_{h}).$$

Noting(3.9) we estimate

$$(\lambda_{2}^{n+1}\mathbf{U}^{n+1/2},d_{t}\mathbf{U}^{n+1})_{h}=\sum_{{\bf z}\in\mathcal{N}_{h}}\beta_{\bf z}\frac{\bigl{(}\nabla\mathbf{W}^{n+1/2},\mathbf{U}^{n+1/2}({\bf z})\otimes\nabla\varphi_{{\bf z}}\bigr{)}}{\beta_{{\bf z}}\max\{|\mathbf{U}^{n+1/2}({\bf z})|^{2},\varepsilon\}}\mathbf{U}^{n+1/2}({\bf z})d_{t}\mathbf{U}^{n+1}({\bf z})$$

$$\leq C\sum_{{\bf z}\in\mathcal{N}_{h}}\frac{\big{(}|\nabla\mathbf{W}^{n+1/2}|,|\nabla\varphi_{{\bf z}}|\bigr{)}_{\omega_{{\bf z}}}}{\max\{|\mathbf{U}^{n+1/2}({\bf z})|^{2},\varepsilon\}}|\mathbf{U}^{n+1/2}({\bf z})|^{2}|d_{t}\mathbf{U}^{n+1}({\bf z})|\leq\frac{C}{h^{2}}\|\Delta_{h}\mathbf{U}^{n+1/2}\|_{h}\|d_{t}\mathbf{U}^{n+1}\|_{h},$$

where we used

$$\big{(}|\nabla\mathbf{W}^{n+1/2}|,|d_{t}\mathbf{U}^{n+1}({\bf z})||\nabla\varphi_{{\bf z}}|\bigr{)}_{\omega_{{\bf z}}} \leq\|\nabla\mathbf{W}^{n+1/2}\|_{L^{2}(\omega_{\bf z})}\||d_{t}\mathbf{U}^{n+1}({\bf z})||\nabla\varphi_{{\bf z}}|\|_{L^{2}(\omega_{\bf z})}$$

$$\leq\frac{1}{h^{2}}\|\mathbf{W}^{n+1/2}\|_{L^{2}(\omega_{\bf z})}\||d_{t}\mathbf{U}^{n+1}({\bf z})||\varphi_{{\bf z}}|\|_{L^{2}(\omega_{\bf z})}.$$

Hence, we conclude that

$$\frac{1}{2\tau}\left(1-\frac{1}{4}-\frac{C\tau}{h^{2}}\right)\|d_{t}\mathbf{U}^{n+1}\|_{h}^{2}+\frac{1}{2\tau}\left(1-\frac{C\tau}{h^{2}}\right)\|\Delta_{h}\mathbf{U}^{n+1}\|_{h}^{2}+\frac{\tau}{2}\|d_{t}^{2}\mathbf{U}^{n+1}\|_{h}^{2}$$

$$\leq\frac{1}{2\tau}\left(1+C\left(1+\frac{\|\mathbf{U}^{n}\|_{\mathbf{L}^{\infty}}}{\varepsilon}\right)^{2}\right)\|d_{t}\mathbf{U}^{n}\|_{h}^{2}+\frac{1}{2\tau}\left(1+\frac{C\tau}{h^{2}}\right)\|\Delta_{h}\mathbf{U}^{n}\|_{h}^{2}.$$

By the inverse estimate $\|v_{h}\|_{L^{\infty}}\leq h^{-d/2}\|v_{h}\|$, recalling the induction assumption, we deduce for sufficiently small $\tau\leq Ch^{2}$ that

$$\frac{\tau^{2}}{2}\|d_{t}\mathbf{U}^{n+1}\|_{L^{\infty}}^{2} \leq C\tau^{2}h^{-d}\|d_{t}\mathbf{U}^{n+1}\|_{h}^{2}$$

$$\leq\frac{C\tau^{2}h^{-d}}{4}\left[\left(1+C\left(1+\frac{\|\mathbf{U}^{n}\|_{\mathbf{L}^{\infty}}}{\varepsilon}\right)^{2}\right)\|d_{t}\mathbf{U}^{n}\|_{h}^{2}+2\|\Delta_{h}\mathbf{U}^{n}\|_{h}^{2}\right]\leq\tilde{C}\tau^{2}h^{-d}<\frac{1}{2}.$$

It follows that (3.10) is satisfied for sufficiently small $\tau\leq Ch^{d/2}\leq Ch^{2}$ for $d\leq 3$. ∎

Remark 3.1.

Compared to the (second order) wave map equation, which requires a time step restriction $\tau\leq Ch^{\max\{1,d/2\}}$ for the existence of the discrete solution [7, 10], the bi-harmonic wave map equation requires a stronger restriction $\tau\leq Ch^{2}$. Nevertheless, the fixed-point algorithm introduced in Section 7.1 requires a weaker condition $\tau\leq Ch$ for convergence in practical computations; a similar condition is also required for the convergence of the corresponding fixed-point algorithm for the wave map equation [10].

4. Convergence of the numerical approximation

Given the solution $\{\mathbf{U}^{n}\}_{n=0}^{N}\subset\mathbf{V}_{h}$ we consider the interpolants ${\bf u}_{\tau,h},{\bf u}^{-}_{\tau,h},\overline{{\bf u}}_{\tau,h}:\Omega_{T}\rightarrow\mathbb{R}^{3}$ such that for $t\in[t_{n},t_{n+1})\times\Omega$

$${\bf u}_{\tau,h}(t)=\frac{t-t_{n}}{\tau}{\bf U}^{n+1}+\frac{t_{n+1}-t}{\tau}{\bf U}^{n}\,,$$

(4.1) $${\bf u}^{-}_{\tau,h}(t)={\bf U}^{n}\,,\qquad{\bf u}^{+}_{\tau,h}(t)={\bf U}^{n+1}\,,$$

$$\overline{{\bf u}}_{\tau,h}(t):={\bf U}^{n+1/2},$$

and analogously we define the corresponding interpolants for $\{\mathbf{V}^{n}\}_{n=0}^{N}$ ($\mathbf{V}^{n}=\mathrm{d}_{t}\mathbf{U}^{n}$ for $n\geq 1$), and $\{\mathbf{W}^{n}\}_{n=0}^{N}$ ($\mathbf{W}^{n}=-\Delta_{h}\mathbf{U}^{n}$ for $n\geq 1$). In the sequel we replace the subscript $\tau,h$ by $h$ to simplify the notation. The scheme (3.1) can be equivalently written in terms of the above interpolants for $\boldsymbol{\Phi}\in{C^{\infty}_{0}}([0,T];\mathbf{V}_{h})$ as

(4.2) $$\begin{split}\int_{0}^{T}\left[\big{(}\partial_{t}{\bf v}_{h},\boldsymbol{\Phi}\big{)}_{h}+\big{(}\nabla\overline{{\bf w}}_{h},\nabla\boldsymbol{\Phi}\big{)}-\big{(}\lambda_{h}^{+}\overline{{\bf u}}_{h},\boldsymbol{\Phi}\big{)}_{h}\right]\mathrm{d}t=0,\\ \int_{0}^{T}\left[\big{(}\overline{{\bf w}}_{h},\boldsymbol{\Phi}\big{)}-\big{(}\nabla\overline{{\bf u}}_{h},\nabla\boldsymbol{\Phi}\big{)}_{h}\right]\mathrm{d}t=0.\end{split}$$

In the next lemma we restate (4.2) in a form that is suitable for showing convergence to the weak formulation (2.6).

Lemma 4.1.

For all $\boldsymbol{\Psi}\in C^{\infty}_{0}([0,T);C^{\infty}(\Omega;\mathbb{R}^{3}))$ it holds that

(4.3) $$\left|\int_{0}^{T}\left[-\big{(}\partial_{t}{\bf u}_{h},\partial_{t}[{\bf u}_{h}\times\boldsymbol{\Psi}]\big{)}_{h}+\big{(}\Delta_{h}\overline{{\bf u}}_{h},\Delta_{h}\mathcal{I}_{h}[\overline{{\bf u}}_{h}\times\boldsymbol{\Psi}]\big{)}_{h}\right]\mathrm{d}t-\big{(}\mathbf{V}^{0},\mathbf{U}^{0}\times\boldsymbol{\Psi}(0)\big{)}_{h}\right|$$

$$\qquad\leq\left|\int_{0}^{T}\big{(}[{\bf v}_{h}-{\bf v}_{h}^{+}],\partial_{t}[{\bf u}_{h}\times\boldsymbol{\Psi}]\big{)}_{h}\mathrm{d}t\right|+C\tau^{1/2}\mathcal{E}_{h}(\mathbf{V}^{0},\mathbf{U}^{0})^{1/2}\|\boldsymbol{\Psi}\|_{L^{\infty}(\Omega_{T})}.$$

Proof.

We take $\boldsymbol{\Phi}=\mathcal{I}_{h}[\overline{{\bf u}}_{h}\times\boldsymbol{\Psi}]$ with $\boldsymbol{\Psi}\in C^{\infty}_{0}([0,T);C^{\infty}(\Omega;\mathbb{R}^{3})$ in (4.2) and note that the term containing $\lambda_{h}^{+}$ disappears by orthogonality. Next, we rewrite the first term as

(4.4) $$\big{(}\partial_{t}{\bf v}_{h},\overline{{\bf u}}_{h}\times\boldsymbol{\Psi}\big{)}_{h}=\big{(}\partial_{t}{\bf v}_{h},{\bf u}_{h}\times\boldsymbol{\Psi}\big{)}_{h}+\big{(}\partial_{t}{\bf v}_{h},(\overline{{\bf u}}_{h}-{\bf u}_{h})\times\boldsymbol{\Psi}\big{)}_{h}.$$

Noting that

$$\overline{{\bf u}}_{h}-{\bf u}_{h}=\Big{(}\frac{1}{2}(t_{n+1}+t_{n})-t\Big{)}{\bf v}_{h}^{+}\qquad\text{for }t\in[t_{n},t_{n+1}),$$

we estimate $\big{(}\partial_{t}{\bf v}_{h},(\overline{{\bf u}}_{h}-{\bf u}_{h})\times\boldsymbol{\Psi}\big{)}_{h}\leq\frac{\tau}{2}\|{\partial_{t}{\bf v}_{h}}\|_{h}\|{\bf v}_{h}^{+}\|_{h}\|\boldsymbol{\Psi}\|_{L^{\infty}}$. Hence, we use the discrete energy estimate (3.6) to bound the second term in (4.4) as

$$\int_{0}^{T}\big{(}\partial_{t}{\bf v}_{h},(\overline{{\bf u}}_{h}-{\bf u}_{h})\times\boldsymbol{\Psi}\big{)}_{h}\mathrm{d}t \leq C\tau^{1/2}\left(\int_{0}^{T}\tau\|{\partial_{t}{\bf v}_{h}}\|_{h}^{2}\mathrm{d}t\right)^{1/2}\|\boldsymbol{\Psi}\|_{L^{\infty}(\Omega_{T})}$$

$$\leq C\tau^{1/2}\mathcal{E}_{h}(\mathbf{V}^{0},\mathbf{U}^{0})^{1/2}\|\boldsymbol{\Psi}\|_{L^{\infty}(\Omega_{T})}.$$

We use integration by parts and that $\partial_{t}{\bf u}_{h}={\bf v}_{h}^{+}$ to rewrite the the first term in (4.4) as

$$\int_{0}^{T}\big{(}\partial_{t}{\bf v}_{h},{\bf u}_{h}\times\boldsymbol{\Psi}\big{)}_{h}\mathrm{d}t=-\int_{0}^{T}\big{(}{\bf v}_{h},\partial_{t}[{\bf u}_{h}\times\boldsymbol{\Psi}]\big{)}_{h}\mathrm{d}t-\big{(}\mathbf{V}^{0},\mathbf{U}^{0}\times\boldsymbol{\Psi}(0)\big{)}_{h}$$

$$\qquad=\int_{0}^{T}\big{(}-\partial_{t}{\bf u}_{h}-[{\bf v}_{h}-{\bf v}_{h}^{+}],\partial_{t}[{\bf u}_{h}\times\boldsymbol{\Psi}]\big{)}_{h}\mathrm{d}t-\big{(}\mathbf{V}^{0},\mathbf{U}^{0}\times\boldsymbol{\Psi}(0)\big{)}_{h}.$$

By the definition of the discrete Laplacian (2.4) we note that the second equation in (4.2) implies that $\overline{{\bf w}}_{h}=-\Delta_{h}\overline{{\bf u}}_{h}$. Hence, using again (2.4) we rewrite

$$\big{(}\nabla\overline{{\bf w}}_{h},\nabla\mathcal{I}_{h}[\overline{{\bf u}}_{h}\times\boldsymbol{\Psi}]\big{)}=\big{(}\overline{{\bf w}}_{h},-\Delta_{h}\mathcal{I}_{h}[\overline{{\bf u}}_{h}\times\boldsymbol{\Psi}]\big{)}_{h}=\big{(}\Delta_{h}\overline{{\bf u}}_{h},\Delta_{h}\mathcal{I}_{h}[\overline{{\bf u}}_{h}\times\boldsymbol{\Psi}]\big{)}_{h}.$$

The inequality (4.3) then follows from the last three identities above. ∎

Next, we state a discrete analogue of the Aubin-Lions-Simin compactness lemma for the embedding $H^{2}\subset H^{1}$, which provides strong $H^{1}$ convergence of the numerical approximation. This strong convergence is required to identify the limits of the nonlinear terms in the numerical scheme.

Lemma 4.2.

Let $\{u_{h}\}_{h>0}$ be a sequence in $\{L^{2}(0,T,V_{h})\}_{h>0}$ such that

$$\int_{0}^{T}\|u_{h}\|_{h}^{2}\mathrm{d}t\leq C,\quad\int_{0}^{T}\|\Delta_{h}u_{h}\|_{h}^{2}\mathrm{d}t\leq C\text{ and }\int_{0}^{T}\|\partial_{t}u_{h}\|_{h}^{2}\mathrm{d}t\leq C.$$

Then also $\displaystyle\int_{0}^{T}\|\Delta_{h}u_{h}\|_{h}^{2}\mathrm{d}t\leq C$ for $h>0$ and there exists a limit $u\in L^{2}(0,T;H^{2}(\Omega))$ and a subsequence such that

$$u_{h} \rightarrow u \text{ in }L^{2}(0,T;H^{1}(\Omega)),$$

$$\Delta_{h}u_{h} \rightharpoonup\Delta u \text{in }L^{2}(0,T;L^{2}(\Omega)).$$

Proof.

The equivalence of norms (2.3) implies that $\{u_{h}\}$, $\{\Delta_{h}u_{h}\}$ are bounded in $L^{2}(0,T;L^{2}(\Omega))$. Moreover, noting (2.2) we bound by the Cauchy-Schwarz and Young inequalities

(4.5) $$\|\nabla u_{h}\|^{2} =(\nabla u_{h},\nabla u_{h})=-(\Delta_{h}u_{h},u_{h})_{h}\leq|(\Delta_{h}u_{h},u_{h})|+Ch\|\Delta_{h}u_{h}\|\|\nabla u_{h}\|$$

$$\leq(1+C)\|\Delta_{h}u_{h}\|\|u_{h}\|\leq C\left(\|\Delta_{h}u_{h}\|^{2}+\|u_{h}\|^{2}\right),$$

where we also used the inverse inequality $\|\nabla u_{h}\|\leq Ch^{-1}\|u_{h}\|$ and (2.2).

The above inequality along with the bounds for $\{u_{h}\}$, $\{\Delta_{h}u_{h}\}$ implies the sequence $\{u_{h}\}$ is also bounded in $L^{2}(0,T;H^{1}(\Omega))$. Hence, we deduce that there exist a subsequence (not relabeled) such that $u_{h}\rightharpoonup u$ in $L^{2}(0,T;H^{1}(\Omega))$. Moreover, since $\partial_{t}u_{h}$ is bounded in $L^{2}(0,T;L^{2}(\Omega))$ a standard Aubin-Lions-Simon result, see, e.g., [19, Theorem 2.1, Ch. III], implies strong convergence $u_{h}\rightarrow u$ in $L^{2}(0,T;L^{2}(\Omega))$.

The bound $\|\Delta_{h}u_{h}\|_{L^{2}(0,T;L^{2})}\leq C$ implies the existence of a limit $w\in L^{2}(0,T;L^{2}(\Omega))$ such that $\Delta_{h}u_{h}\rightharpoonup w$ in $L^{2}(0,T;L^{2}(\Omega))$. By the weak convergence of $\{u_{h}\}$ in $L^{2}(0,T;H^{1}(\Omega))$ we deduce that $w=\Delta u$ since (taking a strongly converging sequence $\varphi_{h}\rightarrow\varphi$, e.g., $\varphi_{h}=I_{h}\varphi$)

$$\int_{0}^{T}(\nabla u,\nabla\phi)\mathrm{d}t=\lim_{h\rightarrow 0}\int_{0}^{T}(\nabla u_{h},\nabla\phi_{h})\mathrm{d}t=\lim_{h\rightarrow 0}\int_{0}^{T}(-\Delta_{h}u_{h},\phi_{h})_{h}\mathrm{d}t$$

$$=\lim_{h\rightarrow 0}\int_{0}^{T}\big{[}(-\Delta_{h}u_{h},\phi_{h})+\mathcal{O}(h)\big{]}\mathrm{d}t=\lim_{h\rightarrow 0}\int_{0}^{T}(-\Delta_{h}u_{h},\phi)\mathrm{d}t=\int_{0}^{T}(-w,\phi)\mathrm{d}t,$$

where we used that (2.2) implies

$$|(-\Delta_{h}u_{h},\phi_{h})_{h}-(-\Delta_{h}u_{h},\phi_{h})|\leq Ch\|\Delta_{h}u_{h}\|\|\nabla\phi_{h}\|.$$

Since $\nabla u_{h}\rightharpoonup\nabla u$ the strong convergence of the gradient follows if we show the convergence of the norm $\lim_{h\rightarrow 0}\int_{0}^{T}\|\nabla u_{h}\|^{2}\mathrm{d}t=\int_{0}^{T}\|\nabla u\|^{2}\mathrm{d}t$. Consequently, using that $\Delta_{h}u_{h}\rightharpoonup\Delta u$, $u_{h}\rightarrow u$ in $L^{2}(0,T;L^{2}(\Omega))$ and (2.2) we conclude that

$$\lim_{h\rightarrow 0}\int_{0}^{T}\|\nabla u_{h}\|^{2}\mathrm{d}t=\lim_{h\rightarrow 0}\int_{0}^{T}(\nabla u_{h},\nabla u_{h})\mathrm{d}t=\lim_{h\rightarrow 0}\int_{0}^{T}(-\Delta u_{h},u_{h})_{h}\mathrm{d}t$$

$$=\int_{0}^{T}(-\Delta u,u)\mathrm{d}t=\int_{0}^{T}\|\nabla u\|^{2}\mathrm{d}t.$$

∎

Next we formulate a (perturbed) discrete product rule for the discrete Laplace operator. The lemma enables us to use discrete orthogonality property in the sequel, which is required to identify the limit of the second term in on the left-had side of (4.3).

Lemma 4.3.

The following discrete product rule holds for the discrete Laplacian (2.4) on right-angled triangulations:

$$\left|(-\Delta_{h}\mathcal{I}_{h}[v_{h}w],u_{h})_{h}-\big{[}(-\Delta_{h}v_{h})\mathcal{I}_{h}w,u_{h})_{h}+(v_{h}(-\Delta_{h}\mathcal{I}_{h}w),u_{h})_{h}-2(\nabla v_{h}\nabla\mathcal{I}_{h}w,u_{h})\big{]}\right|$$

$$\qquad\leq Ch\|\nabla u_{h}\|\|\nabla v_{h}\|\|\nabla w\|_{L^{\infty}(\Omega)}.$$

for any $u_{h},v_{h}\in V_{h}$ and $w\in C(\overline{\Omega})\cap W^{1,\infty}(\Omega)$.

Proof.

$i)$ We first consider Type-$1$ triangulations. We employ the transformation on the reference simplex $\hat{T}\subset\mathbb{R}^{d}$ with vertices $\hat{p}_{0}={\bf 0}$, $\hat{p}_{i}=\alpha_{i}\mathbf{e}_{i}$, $i=1,\dots,d$. Since $\mathcal{T}_{h}$ is right-angled, to each element $T\in\mathcal{T}_{h}$, there is a rotation matrix $R_{T}$ such that the map $x=\mathcal{R}_{T}(\hat{x})=p_{0}+R_{T}\hat{x}$ maps $\hat{T}$ onto $T$ (for suitable $\alpha_{i}$). We denote $\hat{u}(\hat{x})=u(\mathcal{R}_{T}(\hat{x}))=u(x)$ and note the identity $\nabla_{x}u=R_{T}\nabla_{\hat{x}}\hat{u}$, since $R_{T}^{T}=R_{T}^{-1}$. We note that for $\hat{u}\in\mathcal{P}^{1}(\hat{T})$ it holds that

$$\partial_{\hat{x}_{i}}\hat{u}|_{\hat{T}}=\frac{\hat{u}(\hat{p}_{i})-\hat{u}(\hat{p}_{0})}{\alpha_{i}}.$$

Using the definition of the discrete Laplacian we also note that

(4.6) $$(-\Delta_{h}\mathcal{I}_{h}[v_{h}w],\varphi_{\bf z})_{h}=(\nabla\mathcal{I}_{h}[v_{h}w],\nabla\varphi_{\bf z})=\sum_{T\subset\omega_{z}}(\nabla\mathcal{I}_{h}[v_{h}w],\nabla\varphi_{\bf z})_{T}.$$

Using that $\nabla_{x}u=R_{T}\nabla_{\hat{x}}\hat{u}=:R_{T}\hat{\nabla}\hat{u}$ and that $R_{T}$ is a rotation we observe that

$$(\nabla\mathcal{I}_{h}[v_{h}w],\nabla\varphi_{\bf z})_{T}=(R_{T}\hat{\nabla}\mathcal{I}_{h}[\hat{v}\hat{w}],R_{T}\hat{\nabla}\hat{\varphi}_{\hat{{\bf z}}})_{\hat{T}}=(\hat{\nabla}\mathcal{I}_{h}[\hat{v}\hat{w}],\hat{\nabla}\hat{\varphi}_{\hat{{\bf z}}})_{\hat{T}}.$$

Hence, we may proceed by considering the the reference element $\hat{T}$, below we drop the $\hat{}$ to simplify the notation.

On each interior triangle there are $d+1$ non-zero basis functions associated with its nodes $p_{0},\dots,p_{d}$ which we denote as $\varphi_{\bf z}$, ${\bf z}=p_{0},\dots,p_{d}$. We note that if ${\bf z}=p_{0}$

$$\nabla\varphi_{\bf z}=-(\alpha_{1}^{-1},\dots,\alpha_{d}^{-1})^{T},$$

and for ${\bf z}=p_{i}$ with $i>0$ it holds

$$\nabla\varphi_{\bf z}=(0,\dots,0,\alpha_{i}^{-1},0\dots,0)^{T}.$$

Below we consider the two above cases separately.

For ${\bf z}=p_{i}$, $i>0$ we rewrite

$$\{\nabla\mathcal{I}_{h}[v_{h}w]\}_{i}=\partial_{x_{i}}[v_{h}w]=\frac{[v_{h}w](p_{i})-[v_{h}w](p_{0})}{\alpha_{i}}=\frac{v_{h}(p_{i})(w(p_{i})-w(p_{0}))}{\alpha_{i}}+\frac{(v_{h}(p_{i})-v_{h}(p_{0}))w(p_{0})}{\alpha_{i}}$$

$$=\frac{v_{h}(p_{i})(w(p_{i})-w(p_{0}))}{\alpha_{i}}+\frac{(v_{h}(p_{i})-v_{h}(p_{0}))w(p_{i})}{\alpha_{i}}-\frac{(v_{h}(p_{i})-v_{h}(p_{0}))(w(p_{i})-w(p_{0}))}{\alpha_{i}}$$

Hence, we deduce that

(4.7) $$\langle\nabla\mathcal{I}_{h}[v_{h}w],\nabla\varphi_{\bf z}\rangle=v_{h}({\bf z})\langle\nabla\mathcal{I}_{h}w,\nabla\varphi_{\bf z}\rangle+w({\bf z})\langle\nabla v_{h},\nabla\varphi_{\bf z}\rangle-\partial_{x_{i}}v_{h}\partial_{x_{i}}\mathcal{I}_{h}w.$$

In the case that ${\bf z}=p_{0}$ we proceed similarly as in the previous case and arrive at

(4.8) $$\langle\nabla\mathcal{I}_{h}[v_{h}w],\nabla\varphi_{\bf z}\rangle=v_{h}({\bf z})\langle\nabla\mathcal{I}_{h}w,\nabla\varphi_{\bf z}\rangle+w({\bf z})\langle\nabla v_{h},\nabla\varphi_{\bf z}\rangle-\langle\nabla v_{h},\nabla\mathcal{I}_{h}w\rangle.$$

For each $T\in\mathcal{T}_{h}$, $\overline{T}=\mathrm{conv}\{p_{0},\dots,p_{d}\}$, noting that $u_{h}|_{T}=\sum_{i=0}^{d}u_{p_{i}}\varphi_{p_{i}}:=\sum_{i=0}^{d}u_{i}\varphi_{i}$, we deduce from (4.7), (4.8) that

$$(\nabla\mathcal{I}_{h}[v_{h}w],\nabla u_{h})_{T}=\sum_{i=0}^{d}u_{i}(\nabla\mathcal{I}_{h}[v_{h}w],\nabla\varphi_{i})_{T}$$

$$=\sum_{i=0}^{d}\Big{(}v_{i}u_{i}\langle\nabla\mathcal{I}_{h}w,\nabla\varphi_{i}\rangle|T|+w_{i}u_{i}\langle\nabla v_{h},\nabla\varphi_{i}\rangle|T|\Big{)}-u_{0}\langle\nabla v_{h},\nabla\mathcal{I}_{h}w\rangle|T|-\sum_{i=1}^{d}u_{i}\partial_{x_{i}}v_{h}\partial_{x_{i}}\mathcal{I}_{h}w|T|$$

$$:=I_{T}-II_{T}-III_{T}.$$

By the definition of the discrete Laplacian we deduce that

$$\sum_{T\in\mathcal{T}_{h}}I_{T} =\sum_{T\in\mathcal{T}_{h}}\left(\big{(}\nabla\mathcal{I}_{h}w,\nabla[u_{h}v_{h}]\big{)}_{T}+\big{(}\nabla v_{h},\nabla\mathcal{I}_{h}[u_{h}w]\big{)}_{T}\right)=\big{(}\nabla\mathcal{I}_{h}w,\nabla[u_{h}v_{h}]\big{)}+\big{(}\nabla v_{h},\nabla\mathcal{I}_{h}[u_{h}w]\big{)}$$

(4.9) $$=(v_{h}(-\Delta_{h}\mathcal{I}_{h}w),u_{h})_{h}+((-\Delta_{h}v_{h})\mathcal{I}_{h}w,u_{h})_{h}.$$

Next, we consider the term $III_{T}$

$$III_{T}=\sum_{i=1}^{d}u_{i}\partial_{x_{i}}v_{h}\partial_{x_{i}}\mathcal{I}_{h}w|T|=u_{0}\langle\nabla v_{h},\nabla\mathcal{I}_{h}w\rangle|T|+\sum_{i=1}^{d}(u_{i}-u_{0})\partial_{x_{i}}v_{h}\partial_{x_{i}}\mathcal{I}_{h}w|T|:=III_{T,1}+III_{T,2}.$$

We estimate the second term as

$$III_{T,2}=\sum_{i=1}^{d}\alpha_{i}\partial_{x_{i}}u_{h}\partial_{x_{i}}v_{h}\partial_{x_{i}}\mathcal{I}_{h}w|T|\leq Ch\|\nabla u_{h}\|_{T}\|\nabla v_{h}\|_{T}\|\nabla w\|_{L^{\infty}(T)}.$$

We add $III_{T,1}$ to $II_{T}$ and get

$$II_{T}+III_{T,1}=2u_{0}\langle\nabla v_{h},\nabla\mathcal{I}_{h}w\rangle|T|=2\sum_{i=0}^{d}\frac{u_{0}}{d+1}\langle\nabla v_{h},\nabla\mathcal{I}_{h}w\rangle|T|$$

$$=2\sum_{i=0}^{d}\frac{u_{i}}{d+1}\langle\nabla v_{h},\nabla\mathcal{I}_{h}w\rangle|T|+2\sum_{i=1}^{d}\frac{(u_{0}-u_{i})}{d+1}\langle\nabla v_{h},\nabla\mathcal{I}_{h}w\rangle|T|$$

$$=2(\nabla v_{h}\nabla\mathcal{I}_{h}w,u_{h})_{T}+2\sum_{i=1}^{d}\frac{(u_{0}-u_{i})}{d+1}\langle\nabla v_{h},\nabla\mathcal{I}_{h}w\rangle|T|,$$

where we used $\varphi_{i}\equiv\varphi_{p_{i}}(p_{i})=1$ to deduce the last equality.

We estimate the last term on the right-hand side in the above equality as

$$\sum_{i=1}^{d}(u_{0}-u_{i})\langle\nabla v_{h},\nabla\mathcal{I}_{h}w\rangle|T|=\sum_{i=1}^{d}\alpha_{i}\partial_{x_{i}}u_{h}\langle\nabla v_{h},\nabla\mathcal{I}_{h}w\rangle|T|\leq Ch\|\nabla u_{h}\|_{T}\|\nabla v_{h}\|_{T}\|\nabla w\|_{L^{\infty}(T)}$$

Collecting the above estimates for $II_{T}$, $III_{T}$, summing over $T\in\mathcal{T}_{h}$, noting (4) and (4.6) along with the identity $(-\Delta_{h}\mathcal{I}_{h}[v_{h}w],u_{h})_{h}=\sum_{{\bf z}\in\mathcal{N}_{h}}(\nabla\mathcal{I}_{h}[v_{h}w],u_{h}({\bf z})\nabla\varphi_{\bf z})$ concludes the proof.

$ii)$ For Type-$2$ triangulations we may proceed analogously as in the first part (recall that in this case $d=3$). For simplicity we assume a uniform partition with mesh size $h$, a generalization to non-uniform partitions is straightforward. As in the first part, by rotation and translation, all tetrahedra can be mapped onto a reference simplex $\hat{T}=\mathrm{conv}\{\hat{p}_{0},\dots,\hat{p}_{3}\}$ where $\hat{p}_{0}={\bf 0}$, $\hat{p}_{1}=h\mathbf{e}_{1}$, $\hat{p}_{2}=h(\mathbf{e}_{1}+\mathbf{e}_{2})$ and $\hat{p}_{3}=h\mathbf{e}_{3}$. Consequently we observe that for $\hat{u}\in\mathcal{P}^{1}(\hat{T})$ it holds that

$$\partial_{\hat{x}_{1}}\hat{u}|_{\hat{T}}=\frac{\hat{u}(\hat{p}_{1})-\hat{u}(\hat{p}_{0})}{h},\quad\partial_{\hat{x}_{2}}\hat{u}|_{\hat{T}}=\frac{\hat{u}(\hat{p}_{2})-\hat{u}(\hat{p}_{1})}{h},\quad\partial_{\hat{x}_{3}}\hat{u}|_{\hat{T}}=\frac{\hat{u}(\hat{p}_{1})-\hat{u}(\hat{p}_{0})}{h}.$$

The rest of the proof follows as in part $i)$ with minor modifications, therefore we omit it. ∎

Remark 4.1.

A closer inspection of the proof of Lemma 4.3 reveals that in the one dimensional case $\Omega\subset\mathbb{R}^{1}$ the discrete product rule simplifies to an equality which is analogous to the continuous product rule:

(4.10) $$(-\Delta_{h}\mathcal{I}_{h}[v_{h}w],u_{h})_{h}=(-\Delta_{h}v_{h})\mathcal{I}_{h}w,u_{h})_{h}+(v_{h}(-\Delta_{h}\mathcal{I}_{h}w),u_{h})_{h}-2(\nabla v_{h}\nabla\mathcal{I}_{h}w,u_{h}).$$

The discrete product rule (4.10) is used in the next theorem to conclude convergence of the numerical approximation for $d=1$.

Theorem 4.1 (Convergence).

Let Assumption 1 hold and let $\Omega\subset\mathbb{R}^{1}$. Furthermore, let $\mathbf{U}^{0}\rightarrow\mathbf{u}_{0}$, $\mathbf{V}^{0}\rightarrow\mathbf{v}_{0}$ for $h\rightarrow 0$, $\mathcal{E}_{h}(\mathbf{U}^{0},\mathbf{V}^{0})\leq C$ and $|{\bf U}^{0}({\bf z})|=1$, $\langle\mathbf{U}^{0}({\bf z}),\mathbf{V}^{0}({\bf z})\rangle=0$ for ${\bf z}\in\mathcal{N}_{h}$. Then there exists a weak solution ${\bf u}:[0,T]\times\Omega\rightarrow\mathbb{S}^{2}$ of (1.1) according to Definition 2.1 and a subsequence $\{{\bf u}_{h}\}_{h,\tau}$ such that for $(\tau,h)\rightarrow 0$

$${\bf u}_{h} \overset{\ast}{\rightharpoonup}{\bf u} \text{ in }L^{\infty}(0,T;{\bf H}^{1}),$$

$$\partial_{t}{\bf u}_{h} \overset{\ast}{\rightharpoonup}\partial_{t}{\bf u} \text{ in }L^{\infty}(0,T;{\bf L}^{2}).$$

Proof.

Estimate (3.6) and Lemma 4.2 implies the boundedness of the sequence $\{{\bf u}_{h}\}_{h,\tau>0}$ in $L^{\infty}(0,T;{\bf H}^{1})$. This together with the bound on $\{\partial_{t}{\bf u}_{h}\}_{h,\tau>0}$ implies the (sub)convergence

(4.11) $${\bf u}_{h},{\bf u}_{h}^{+},\overline{{\bf u}}_{h} \overset{\ast}{\rightharpoonup}{\bf u} \text{ in }L^{\infty}(0,T;{\bf H}^{1}),$$

(4.12) $${\bf u}_{h},{\bf u}_{h}^{+},\overline{{\bf u}}_{h} \rightarrow{\bf u} \text{ in }L^{2}(0,T;{\bf L}^{2}),$$

(4.13) $$\partial_{t}{\bf u}_{h},{\bf v}_{h},{\bf v}_{h}^{+} \overset{\ast}{\rightharpoonup}\partial_{t}{\bf u} \text{ in }L^{\infty}(0,T;{\bf L}^{2}),$$

(4.14) $$\Delta_{h}{\bf u}_{h},\Delta_{h}\overline{{\bf u}}_{h},\Delta_{h}{\bf u}_{h}^{+} \overset{\ast}{\rightharpoonup}\Delta{\bf u} \text{ in }L^{\infty}(0,T;{\bf L}^{2}),$$

where to show the convergence of ${\bf v}_{h}$, ${\bf v}_{h}^{+}$ to the same limit for $\tau,h\rightarrow 0$ can be concluded thanks to the the numerical dissipation term in (3.6), since,

(4.15) $$\|{\bf v}_{h}-{\bf v}_{h}^{+}\|_{L^{2}(0,T;{\bf L}^{2})}^{2}\leq C\tau^{3}\sum_{n=1}^{N}\|d_{t}\mathbf{V}^{n}\|_{h}^{2}\leq C\tau\rightarrow 0.$$

Similarly we also conclude that the limits of $\Delta_{h}{\bf u}_{h},\Delta_{h}\overline{{\bf u}}_{h},\Delta_{h}{\bf u}_{h}^{+}$ coincide.

Furthermore, Lemma 4.2 implies the strong convergence

(4.16) $$\nabla{\bf u}_{h},\nabla\overline{{\bf u}}_{h},\nabla{\bf u}_{h}^{+}\rightarrow\nabla{\bf u}\qquad\text{ in }L^{2}(0,T;{\bf L}^{2}).$$

Since $|{\bf u}_{h}^{+}({\bf z})|=1$ for any $T\in\mathcal{T}_{h}$ and ${\bf z}\in\mathcal{N}_{h}$ we have for ${\bf x},{\bf z}\in T$ that $|{\bf u}_{h}^{+}({\bf x})|-1|=||{\bf u}_{h}^{+}({\bf x})|-|{\bf u}_{h}^{+}({\bf z})||\leq\mathrm{diam}(T)|\nabla{\bf u}_{h}^{+}|_{T}|$. Hence we obtain

$$\big{\|}|{\bf u}_{h}^{+}|-1\big{\|}_{L^{2}(T)}\leq Ch\|\nabla{\bf u}_{h}^{+}\|_{L^{2}(T)}.$$

This implies that $|{\bf u}_{h}^{+}|\rightarrow 1$ almost everywhere in $\Omega_{T}$ for $(\tau,h)\rightarrow 0$ and consequently $|{\bf u}|=1$ almost everywhere in $\Omega_{T}$.

To show that the limiting function ${\bf u}$ is a weak solution we pass to the limit in (4.3) for $(\tau,h)\rightarrow 0$. Noting the orthogonality of $\Delta_{h}\overline{{\bf u}}_{h}$ and $\Delta_{h}\overline{{\bf u}}_{h}\times\mathcal{I}_{h}\boldsymbol{\Psi}$ we deduce by the discrete product rule (4.10) that

(4.17) $$\big{(}\Delta_{h}\overline{{\bf u}}_{h},\Delta_{h}\mathcal{I}_{h}[\overline{{\bf u}}_{h}\times\boldsymbol{\Psi}]\big{)}_{h}=\big{(}\Delta_{h}\overline{{\bf u}}_{h},\overline{{\bf u}}_{h}\times\Delta_{h}\mathcal{I}_{h}\boldsymbol{\Psi}\big{)}_{h}+2\big{(}\Delta_{h}\overline{{\bf u}}_{h},\nabla\overline{{\bf u}}_{h}\times\nabla\mathcal{I}_{h}\boldsymbol{\Psi}\big{)}.$$

To estimate the first term in (4.17) we bound

$$|\big{(}\Delta_{h}\overline{{\bf u}}_{h},\overline{{\bf u}}_{h}\times\Delta_{h}\mathcal{I}_{h}\boldsymbol{\Psi}\big{)}_{h}-\big{(}\Delta{\bf u},{\bf u}\times\Delta\boldsymbol{\Psi}\big{)}|\leq|\big{(}\Delta_{h}\overline{{\bf u}}_{h},\overline{{\bf u}}_{h}\times\Delta_{h}\mathcal{I}_{h}\boldsymbol{\Psi}\big{)}_{h}-\big{(}\Delta_{h}\overline{{\bf u}}_{h},\overline{{\bf u}}_{h}\times\Delta\boldsymbol{\Psi}\big{)}_{h}|$$

$$\qquad+|\big{(}\Delta_{h}\overline{{\bf u}}_{h},\overline{{\bf u}}_{h}\times\Delta\boldsymbol{\Psi}\big{)}_{h}-\big{(}\Delta_{h}\overline{{\bf u}}_{h},\overline{{\bf u}}_{h}\times\mathcal{I}_{h}\Delta\boldsymbol{\Psi}\big{)}|+|\big{(}\Delta_{h}\overline{{\bf u}}_{h},\overline{{\bf u}}_{h}\times\mathcal{I}_{h}\Delta\boldsymbol{\Psi}\big{)}-\big{(}\Delta_{h}\overline{{\bf u}}_{h},\overline{{\bf u}}_{h}\times\Delta\boldsymbol{\Psi}\big{)}|$$

$$\qquad+|\big{(}\Delta_{h}\overline{{\bf u}}_{h},\overline{{\bf u}}_{h}\times\Delta\boldsymbol{\Psi}\big{)}-\big{(}\Delta_{h}\overline{{\bf u}}_{h},{\bf u}\times\Delta\boldsymbol{\Psi}\big{)}|+|\big{(}\Delta_{h}\overline{{\bf u}}_{h},{\bf u}\times\Delta\boldsymbol{\Psi}\big{)}-\big{(}\Delta{\bf u},{\bf u}\times\Delta\boldsymbol{\Psi}\big{)}|$$

$$:=A_{1}+A_{2}+A_{3}+A_{4}+A_{5}.$$

We estimate the first term by the Cauchy-Schwarz inequality and $|\mathbf{u}_{h}|\leq 1$ as

$$A_{1}\leq C\|\Delta_{h}\overline{{\bf u}}_{h}\|_{h}\|\Delta_{h}\mathcal{I}_{h}\boldsymbol{\Psi}-\Delta\boldsymbol{\Psi}\|_{h}.$$

Hence, noting Assumption 1 and Lemma 3.1 we deduce that $\lim_{h\rightarrow 0}\int_{0}^{T}A_{1}\,\mathrm{d}t=0.$

To estimate the second term we use (2.2), $|{\bf u}_{h}|\leq 1$ and $W^{1,\infty}$-stability of the interpolation operator

$$A_{2}\leq Ch\|\Delta_{h}\overline{{\bf u}}_{h}\|\|\nabla(\overline{{\bf u}}_{h}\times\mathcal{I}_{h}\Delta\boldsymbol{\Psi})\|\leq Ch\|\Delta_{h}\overline{{\bf u}}_{h}\|\left(\|\nabla\overline{{\bf u}}_{h}\|+1\right)\|\Delta\boldsymbol{\Psi}\|_{W^{1,\infty}}.$$

Consequently, noting Lemma 3.1 and $\|\nabla\overline{{\bf u}}_{h}\|_{L^{\infty}(0,T;\mathbf{L}^{2})}\leq C$ we get that $\lim_{h\rightarrow 0}\int_{0}^{T}A_{2}\,\mathrm{d}t=0$ .

Next, Lemma 3.1, $|\mathbf{u}_{h}|\leq 1$ and (2.1) imply that

$$A_{3}\leq Ch^{2}\|\Delta_{h}\overline{{\bf u}}_{h}\||\Delta\boldsymbol{\Psi}|_{W^{2,2}}\leq Ch^{2},$$

and consequently $\lim_{h\rightarrow 0}\int_{0}^{T}A_{3}\,\mathrm{d}t=0$.

Using the Cauchy-Schwarz inequality we estimate

$$A_{4}\leq\|\Delta_{h}\overline{{\bf u}}_{h}\|\|\overline{{\bf u}}_{h}-{\bf u}\|\|\Delta\boldsymbol{\Psi}\|_{\mathbf{L}^{\infty}},$$

and consequently Lemma 3.1 and (4.12) imply that $\lim_{h\rightarrow 0}\int_{0}^{T}A_{4}\,\mathrm{d}t=0.$ For the last term we conclude by the weak convergence (4.14) that $\lim_{\tau,h\rightarrow 0}\int_{0}^{T}A_{5}\,\mathrm{d}t=0$.

Collecting the above limits for $A_{1}$, $\dots$, $A_{5}$ we conclude that

(4.18) $$\lim_{\tau,h\rightarrow 0}\int_{0}^{T}\big{(}\Delta_{h}\overline{{\bf u}}_{h},\overline{{\bf u}}_{h}\times\Delta_{h}\mathcal{I}_{h}\boldsymbol{\Psi}\big{)}_{h}\mathrm{d}t=\int_{0}^{T}\big{(}\Delta{\bf u},{\bf u}\times\Delta\boldsymbol{\Psi}\big{)}\mathrm{d}t.$$

The second term in (4.17) is estimated as

$$|\big{(}\Delta_{h}\overline{{\bf u}}_{h},\nabla\overline{{\bf u}}_{h}\times\nabla\mathcal{I}_{h}\boldsymbol{\Psi}\big{)}-\big{(}\Delta{{\bf u}},\nabla{{\bf u}}\times\nabla\boldsymbol{\Psi}\big{)}|\leq|\big{(}\Delta_{h}\overline{{\bf u}}_{h},\nabla\overline{{\bf u}}_{h}\times\nabla\mathcal{I}_{h}\boldsymbol{\Psi}\big{)}-\big{(}\Delta_{h}\overline{{\bf u}}_{h},\nabla\overline{{\bf u}}_{h}\times\nabla\boldsymbol{\Psi}\big{)}|$$

$$\qquad+|\big{(}\Delta_{h}\overline{{\bf u}}_{h},\nabla\overline{{\bf u}}_{h}\times\nabla\boldsymbol{\Psi}\big{)}-\big{(}\Delta_{h}\overline{{\bf u}}_{h},\nabla{\bf u}\times\nabla\boldsymbol{\Psi}\big{)}|+|\big{(}\Delta_{h}\overline{{\bf u}}_{h},\nabla{\bf u}\times\nabla\boldsymbol{\Psi}\big{)}-\big{(}\Delta{\bf u},\nabla{\bf u}\times\nabla\boldsymbol{\Psi}\big{)}|$$

$$:=B_{1}+B_{2}+B_{3}$$

Using Lemma 3.1 along with the bound $\|\nabla{\bf u}_{h}\|_{L^{\infty}(0,T;{\bf L}^{2})}\leq C$ and (2.1) we estimate

$$B_{1}\leq\|\Delta_{h}\overline{{\bf u}}_{h}\|\|\nabla\overline{{\bf u}}_{h}\|\|\mathcal{I}_{h}\boldsymbol{\Psi}-\boldsymbol{\Psi}\|_{\mathbf{W}^{1,\infty}}\leq Ch|\boldsymbol{\Psi}|_{\mathbf{W}^{2,\infty}}.$$

The second term can be estimate by Lemma 3.1 as

$$B_{2}\leq\|\Delta_{h}\overline{{\bf u}}_{h}\|\|\nabla\overline{{\bf u}}_{h}-\nabla{\bf u}\|\|\boldsymbol{\Psi}\|_{\mathbf{W}^{1,\infty}}.$$

Consequently, noting the convergences (4.14), (4.16) we conclude that

(4.19) $$\lim_{\tau,h\rightarrow 0}\int_{0}^{T}\big{(}\Delta_{h}\overline{{\bf u}}_{h},\nabla\overline{{\bf u}}_{h}\times\nabla\mathcal{I}_{h}\boldsymbol{\Psi}\big{)}\mathrm{d}t=\int_{0}^{T}\big{(}\Delta{{\bf u}},\nabla{{\bf u}}\times\nabla\boldsymbol{\Psi}\big{)}\mathrm{d}t.$$

The remaining terms can be treated as in [7, Theorem 4.1]. I.e., we have that

$$\lim_{\tau,h\rightarrow 0}\int_{0}^{T}\big{(}\partial_{t}{\bf u}_{h},\partial_{t}[{\bf u}_{h}\times\boldsymbol{\Psi}]\big{)}_{h}\mathrm{d}t=\int_{0}^{T}\big{(}\partial_{t}{\bf u},\partial_{t}[{\bf u}\times\boldsymbol{\Psi}]\big{)}\mathrm{d}t,$$

and

$$\lim_{\tau,h\rightarrow 0}\big{(}\mathbf{V}^{0},\mathbf{U}^{0}\times\boldsymbol{\Psi}(0)\big{)}_{h}=\big{(}\mathbf{v}^{0},\mathbf{u}^{0}\times\boldsymbol{\Psi}(0)\big{)},$$

for all $\boldsymbol{\Psi}\in C^{\infty}_{0}([0,T);C^{\infty}(\Omega;\mathbb{R}^{3}))$. Furthermore, owing to (4.15), we get for the first term on the right-hand side of (4.3) that

$$\lim_{\tau\rightarrow 0}\int_{0}^{T}\big{(}[{\bf v}_{h}-{\bf v}_{h}^{+}],\partial_{t}[{\bf u}_{h}\times\boldsymbol{\Psi}]\big{)}_{h}\mathrm{d}t=0.$$

Collecting last three eqalities along with (4.18), (4.19) we conclude from (4.3) that the limit ${\bf u}$ satisfies (2.6). ∎

5. Higher-dimensional case

In the higher-dimensional case $d>1$ we employ Lemma 4.3 (instead of (4.17) for $d=1$) and get that

(5.1) $$\begin{split}\lim_{\tau,h\rightarrow 0}\big{(}\Delta_{h}\overline{{\bf u}}_{h},\Delta_{h}\mathcal{I}_{h}[\overline{{\bf u}}_{h}\times\boldsymbol{\Psi}]\big{)}_{h}\leq&\lim_{\tau,h\rightarrow 0}\left(\big{(}\Delta_{h}\overline{{\bf u}}_{h},\overline{{\bf u}}_{h}\times\Delta_{h}\mathcal{I}_{h}\boldsymbol{\Psi}\big{)}_{h}+2\big{(}\Delta_{h}\overline{{\bf u}}_{h},\nabla\overline{{\bf u}}_{h}\times\nabla\mathcal{I}_{h}\boldsymbol{\Psi}\big{)}\right)\\ &+\lim_{\tau,h\rightarrow 0}Ch\|\nabla\Delta_{h}\overline{{\bf u}}_{h}\|\|\nabla\overline{{\bf u}}_{h}\|\|\nabla\boldsymbol{\Psi}\|_{{\bf L}^{\infty}}.\end{split}$$

When passing to the limit in the discrete formulation we need to show that the perturbation term $h\|\nabla\Delta_{h}\overline{{\bf u}}_{h}\|\|\nabla\overline{{\bf u}}_{h}\|\|\nabla\boldsymbol{\Psi}\|_{{\bf L}^{\infty}}$ disappears for $h\rightarrow 0$. However we lack an $h$-independent estimate for $\nabla\Delta_{h}\overline{{\bf u}}_{h}$. Estimating the term $\|\nabla\Delta_{h}\overline{{\bf u}}_{h}\|$ by an inverse estimate is not viable since it would introduce a factor $h^{-1}$ which would cancel the factor $h$ in the last term in (5.1). Under the assumption that $h^{\alpha/2}\|\nabla\Delta_{h}\overline{{\bf u}}_{h}\|<C$ for some $0<\alpha<2$ (cf. (5.10) below) we can control the perturbation term as follows

(5.2) $$h\int_{0}^{T}\|\nabla\Delta_{h}\overline{{\bf u}}_{h}\|\|\nabla\overline{{\bf u}}_{h}\|\|\nabla\boldsymbol{\Psi}\|_{{\bf L}^{\infty}}\mathrm{d}t\leq Ch^{1-\alpha/2}\|\nabla\boldsymbol{\Psi}\|_{L^{\infty}(0,T;{\bf L}^{\infty})}\rightarrow 0\text{ for }h\rightarrow 0.$$

To guarantee a $h$-independent bound of $\nabla\Delta_{h}\overline{{\bf u}}_{h}$ we introduce an artificial stabilization term in the numerical formulation (3.1). The corresponding numerical approximation is defined as follows: given $\mathbf{U}^{0}$, $\mathbf{V}^{0}$ set $\mathbf{U}^{-1}=\mathbf{U}^{0}-\tau\mathbf{V}^{0}$ and for $n=0,\dots,N-1$ determine $\mathbf{U}^{n+1},\mathbf{W}^{n+1}\in{\bf V}_{h}$ such that for all $\boldsymbol{\Phi},\boldsymbol{\Psi}\in{\bf V}_{h}$ there holds

(5.3) $$\bigl{(}d_{t}^{2}\mathbf{U}^{n+1},\boldsymbol{\Phi}\bigr{)}_{h}+(\nabla\mathbf{W}^{n+1/2},\nabla\boldsymbol{\Phi}){+}h^{\alpha}\left(\Delta_{h}\mathbf{W}^{n+1/2},\Delta_{h}\boldsymbol{\Phi}\right) =\bigl{(}\lambda_{w}^{n+1}\mathbf{U}^{n+1/2},\boldsymbol{\Phi}\bigr{)}_{h}\,,$$

$$(\mathbf{W}^{n+1/2},\boldsymbol{\Psi})_{h} =(\nabla\mathbf{U}^{n+1/2},\nabla\boldsymbol{\Psi})\,,$$

where $0<\alpha<2$ and the discrete Lagrange multiplier $\lambda_{w}^{n+1}\in V_{h}$ is defined for ${\bf z}\in\mathcal{N}_{h}$ as

(5.9) $$\lambda_{w}^{n+1}({\bf z})=\left\{\begin{array}[]{ll}0&\mbox{ if }\mathbf{U}^{n+1/2}({\bf z})={\bf 0},\\ \frac{\frac{1}{2}\langle\mathbf{V}^{0}({\bf z}),\mathbf{V}^{1}({\bf z})\rangle}{|\mathbf{U}^{1/2}({\bf z})|^{2}}+\frac{\bigl{(}\nabla\mathbf{W}^{1/2},\mathbf{U}^{1/2}({\bf z})\otimes\nabla\varphi_{{\bf z}}\bigr{)}}{\beta_{{\bf z}}|\mathbf{U}^{1/2}({\bf z})|^{2}}\\ +\frac{h^{\alpha}\bigl{(}\Delta_{h}\mathbf{W}^{1/2},\mathbf{U}^{1/2}({\bf z})\Delta_{h}\varphi_{{\bf z}}\bigr{)}}{\beta_{{\bf z}}|\mathbf{U}^{1/2}({\bf z})|^{2}}&\text{if }n=0,\,\mathbf{U}^{1/2}({\bf z})\neq{\bf 0},\\ -\frac{\bigl{\langle}d_{t}\mathbf{U}^{n}({\bf z}),{d_{t}}\mathbf{U}^{n+1/2}({\bf z})\bigr{\rangle}}{|\mathbf{U}^{n+1/2}({\bf z})|^{2}}+\frac{\bigl{(}\nabla\mathbf{W}^{n+1/2},\mathbf{U}^{n+1/2}({\bf z})\otimes\nabla\varphi_{{\bf z}}\bigr{)}}{\beta_{{\bf z}}|\mathbf{U}^{n+1/2}({\bf z})|^{2}}\\ +\frac{h^{\alpha}\bigl{(}\Delta_{h}\mathbf{W}^{n+1/2},\mathbf{U}^{n+1/2}({\bf z})\Delta_{h}\varphi_{{\bf z}}\bigr{)}}{\beta_{{\bf z}}|\mathbf{U}^{n+1/2}({\bf z})|^{2}}&\mbox{ else}\,.\end{array}\right.$$

Setting $\boldsymbol{\Phi}=\mathbf{U}^{n+1}-\mathbf{U}^{n}$, $\boldsymbol{\Psi}=\Delta_{h}(\mathbf{U}^{n+1}-\mathbf{U}^{n})$ in (5.3) and noting that $\mathbf{W}^{n}=-\Delta_{h}\mathbf{U}^{n}$ we deduce that the following stability estimate holds

(5.10) $$\max_{n=1,\dots,N}\left[\mathcal{E}_{h}\bigl{(}{\bf V}^{n},{\bf U}^{n}\bigr{)}+\frac{h^{\alpha}}{2}\|\nabla\mathbf{W}^{n}\|^{2}\right]\leq\mathcal{E}_{h}\bigl{(}{\bf V}^{0},{\bf U}^{0}\bigr{)}+\frac{h^{\alpha}}{2}\|\nabla\mathbf{W}^{0}\|^{2}\,,$$

5.1. Existence

The existence of $\mathbf{U}^{n+1}$, $\mathbf{W}^{n+1}=-\Delta_{h}\mathbf{U}^{n+1}$ which solve (5.3) follow as in Lemma 3.2 under the assumption that $\tau\leq Ch^{3-\alpha/2}$ is sufficiently small. The additional term (5.9) (which imposes the stronger assumption $\tau\leq Ch^{3-\alpha/2}$) can be estimated analogously to (3.9) as

$$\frac{h^{\alpha}\bigl{(}\Delta_{h}\mathbf{W}^{n+1/2},\mathbf{U}^{n+1/2}({\bf z})\Delta_{h}\varphi_{{\bf z}}\bigr{)}}{\beta_{{\bf z}}|\mathbf{U}^{n+1/2}({\bf z})|^{2}} \leq C\frac{h^{\alpha}}{h^{2}}\|\Delta_{h}\mathbf{W}^{n+1/2}\|_{h}\|\mathbf{U}^{n+1/2}\|_{h}\leq C\frac{h^{\alpha}}{2h^{3}}\|\nabla\mathbf{W}^{n+1/2}\|\|\mathbf{U}^{n+1/2}\|_{h}$$

$$\leq\frac{h^{\alpha}}{2}\|\nabla\mathbf{W}^{n+1/2}\|^{2}+Ch^{\alpha-6}\|\mathbf{U}^{n+1/2}\|_{h}^{2}.$$

where we used the inverse estimate (2.5) and that $\|\Delta_{h}\varphi_{{\bf z}}\|_{L^{\infty}}\leq Ch^{-2}$ .

Moreover to show (3.10) we estimate the additional term

$$\sum_{{\bf z}\in\mathcal{N}_{h}} \beta_{{\bf z}}\frac{h^{\alpha}\bigl{(}\Delta_{h}\mathbf{W}^{n+1/2},\mathbf{U}^{n+1/2}({\bf z})\Delta_{h}\varphi_{{\bf z}}\bigr{)}}{\beta_{{\bf z}}|\mathbf{U}^{n+1/2}({\bf z})|^{2}}d_{t}\mathbf{U}^{n+1}({\bf z})\leq Ch^{\alpha-3}\|\nabla\mathbf{W}^{n+1/2}\|\|d_{t}\mathbf{U}^{n+1}\|_{h}$$

$$\leq C{h^{\alpha/2-3}}h^{\alpha}\|\nabla\mathbf{W}^{n+1/2}\|^{2}+C{h^{\alpha/2-3}}\|d_{t}\mathbf{U}^{n+1}\|_{h}^{2},$$

which can be absorbed in the corresponding terms for sufficiently small $\tau\leq C{h^{3-\alpha/2}}$

5.2. Convergence

Under the additional assumption that $h^{\alpha}\|\nabla\mathbf{W}^{0}\|^{2}\leq C$ the convergence of the numerical approximation (5.3) follows as in Theorem 4.1 with minor modifications: instead of (4.17) we consider (5.1) and the perturbation term disappears for $h\rightarrow 0$ thanks to (5.2); moreover thanks to (5.10) the stabilization term vanishes in the limit, i.e., $h^{\alpha}\left(\Delta_{h}\mathbf{W}^{n+1/2},\Delta_{h}\boldsymbol{\Phi}\right)\leq h^{\alpha}\|\Delta_{h}\mathbf{W}^{n+1/2}\|\|\Delta_{h}\boldsymbol{\Phi}\|\rightarrow 0$ for $h\rightarrow 0$. The remainder of the proof is the same as in the case $d=1$.

6. Conservative scheme

We present an alternative approximation scheme that preserves the discrete energy, i.e, the scheme satisfies a (discrete) energy equality (6.2) below.

Let $\mathbf{U}^{0},\mathbf{V}^{0}\in\mathbf{V}_{h}$, $\langle\mathbf{V}^{0}({\bf z}),\mathbf{U}^{0}({\bf z})\rangle=0$, $\mathbf{U}^{-1}=\mathbf{U}^{0}-\tau\mathbf{V}^{0}$, for $n\geq 1$ define $\mathbf{V}^{n}=d_{t}\mathbf{U}^{n}$ and $\mathbf{U}^{n^{*}+\frac{1}{2}}=\frac{1}{2}\left(\mathbf{U}^{n+1}+\mathbf{U}^{n-1}\right)$, $\mathbf{W}^{n^{*}+\frac{1}{2}}=\frac{1}{2}\left(\mathbf{W}^{n+1}+\mathbf{W}^{n-1}\right)$. For $n=0,\dots,N-1$ determine $\mathbf{U}^{n+1}$, $\mathbf{W}^{n+1}$ as the solution of

(6.1) $$\bigl{(}d_{t}\mathbf{V}^{n+1},\boldsymbol{\Phi}\bigr{)}_{h}+(\nabla\mathbf{W}^{n^{*}+\frac{1}{2}},\nabla\boldsymbol{\Phi}) =\bigl{(}\widetilde{\lambda}_{w}^{n+1}\mathbf{U}^{n^{*}+\frac{1}{2}},\boldsymbol{\Phi}\bigr{)}_{h}\,,$$

$$(\mathbf{W}^{n^{*}+\frac{1}{2}},\boldsymbol{\Psi})_{h} =(\nabla\mathbf{U}^{n^{*}+\frac{1}{2}},\nabla\boldsymbol{\Psi})\,,$$

with

$$\widetilde{\lambda}_{w}^{n+1}({\bf z})=\left\{\begin{array}[]{ll}0&\text{if }\mathbf{U}^{n^{*}+1/2}({\bf z})={\bf 0},\\ -\frac{\bigl{\langle}\mathbf{V}^{0}({\bf z}),\mathbf{V}^{1}({\bf z})\bigr{\rangle}-\frac{1}{2}|\mathbf{V}^{0}({\bf z})|^{2}}{|\mathbf{U}^{0^{*}+\frac{1}{2}}({\bf z})|^{2}}+\frac{\bigl{(}\nabla\mathbf{W}^{0^{*}+\frac{1}{2}},\mathbf{U}^{0^{*}+\frac{1}{2}}({\bf z})\otimes\nabla\varphi_{{\bf z}}\bigr{)}}{\beta_{{\bf z}}|\mathbf{U}^{0^{*}+\frac{1}{2}}({\bf z})|^{2}}&\text{if }n=0,\,\mathbf{U}^{0^{*}+\frac{1}{2}}({\bf z})\neq{\bf 0}\\ -\frac{\bigl{\langle}\mathbf{V}^{n}({\bf z}),\mathbf{V}^{n+1}({\bf z})\bigr{\rangle}}{|\mathbf{U}^{n^{*}+\frac{1}{2}}({\bf z})|^{2}}+\frac{\bigl{(}\nabla\mathbf{W}^{n^{*}+\frac{1}{2}},\mathbf{U}^{n^{*}+\frac{1}{2}}({\bf z})\otimes\nabla\varphi_{{\bf z}}\bigr{)}}{\beta_{{\bf z}}|\mathbf{U}^{n^{*}+\frac{1}{2}}({\bf z})|^{2}}&\text{else},\end{array}\right.$$

where the modification for $n=0$, accounts for the fact that possibly $|{\bf U}^{-1}|\neq 1$ if $\mathbf{V}_{0}\neq{\bf 0}$, cf. [7], [8].

Setting $\boldsymbol{\Phi}=\mathbf{U}^{n^{*}+\frac{1}{2}}({\bf z})\varphi_{{\bf z}}$ in (6.1) then implies that $|{\bf U}^{n+1}({\bf z})|^{2}=|{\bf U}^{n}({\bf z})|^{2}$ $\forall{\bf z}\in\mathcal{N}_{h}$, $n=0,\dots,N-1$.

Setting $\boldsymbol{\Phi}=\mathbf{U}^{n+1}-\mathbf{U}^{n-1}=\mathbf{V}^{n+1}+\mathbf{V}^{n}$, $\boldsymbol{\Psi}=\Delta_{h}(\mathbf{U}^{n+1}-\mathbf{U}^{n-1})$ in (6.1) implies that above scheme conserves the discrete energy $\widetilde{\mathcal{E}}_{h}\bigl{(}{\bf V}^{n},{\bf U}^{n},{\bf U}^{n-1}\bigr{)}:=\frac{1}{2}\left(\|\mathbf{V}^{n}\|_{h}^{2}+\frac{1}{2}\left[\|\Delta_{h}\mathbf{U}^{n}\|^{2}+\|\Delta_{h}\mathbf{U}^{n-1}\|_{h}^{2}\right]\right)$, i.e., for $n=1,\dots,N$ the following equality holds:

(6.2) $$\widetilde{\mathcal{E}}_{h}\bigl{(}{\bf V}^{n},{\bf U}^{n}\bigr{)}=\frac{1}{2}\left(\|\mathbf{V}^{0}\|_{h}^{2}+\frac{\tau}{2}\left(\widetilde{\lambda}_{w}^{1}\mathbf{V}^{0},\mathbf{V}^{0}\right)_{h}+\frac{1}{2}\left[\|\Delta_{h}\mathbf{U}^{0}\|^{2}+\|\Delta_{h}\mathbf{U}^{-1}\|_{h}^{2}\right]\right).$$

Compared to (3.6) the numerical damping term $\frac{\tau^{2}}{2}\sum_{n=1}^{N}\|d_{t}{\bf V}^{n}\|^{2}_{h}$ is not present in (6.2) and the energy of the initial data is exactly preserved at all times. Due to the lack of numerical damping term $\frac{\tau^{2}}{2}\sum_{n=1}^{N}\|d_{t}{\bf V}^{n}\|^{2}_{h}$ the convergence of the conservative scheme is not clear. In particular we are not able to show that the right hand side in the (counterpart of) Lemma 4.1 vanishes for $\tau\rightarrow 0$.

7. Numerical experiments

In this section we present numerical simulation results using the secheme (3.1). The results computed with the conservative scheme (6.1) were simular on finite time intervals. In the long-time limit the results will inevitably diverge, since the numerical solution of (3.1) will eventually become stationary, due to the effects of the numerical dissipation.

7.1. Solution of the nonlinear systems

For every $n\geq 1$, we solve the nonlinear problem (3.1)-(3.5) by a fixed point algorithm with projection. Given a tolerance $\epsilon>0$ we proceed for $n=1,\dots,N$ as follows: set ${\bf u}^{0}={\bf u}^{0}_{*}={\bf U}^{n}$, ${\bf w}^{0}={\bf w}^{0}_{*}=-\Delta_{h}{\bf U}^{n}$ and iterate for $l=0,\dots$

1. • Find ${\bf u}^{l+1},\mathbf{w}^{l+1}=-\Delta_{h}{\bf u}^{l+1}\in{\bf V}_{h}$ such that for all $\boldsymbol{\Phi}\in{\bf V}_{h}$

     $$\begin{array}[]{rcl}\frac{1}{\tau}\Bigl{(}{\bf u}^{l+1},\boldsymbol{\Phi}\Bigr{)}_{h}+\frac{\tau}{2}\Bigl{(}\nabla{\bf w}^{l+1},\nabla\boldsymbol{\Phi}\Bigr{)}-\frac{\tau}{2}\Bigl{(}\lambda^{l}_{w*}{\bf u}^{l+1},\boldsymbol{\Phi}\Bigr{)}_{h}&=&\frac{1}{\tau}\Bigl{(}2{\bf U}^{n}-{\bf U}^{n-1},\boldsymbol{\Phi}\Bigr{)}_{h}-\frac{\tau}{2}\Bigl{(}\nabla{\bf W}^{n},\nabla\boldsymbol{\Phi}\Bigr{)}\\ &&+\frac{\tau}{2}\Bigl{(}{\lambda}^{l}_{w*}{\bf U}^{n},\boldsymbol{\Phi}\Bigr{)}_{h}\,.\end{array}$$

   • For all ${\bf z}\in\mathcal{N}_{h}$ set

     $${\bf u}^{l+1}_{*}({\bf z}) =\frac{{\bf u}^{l+1}({\bf z})}{|{\bf u}^{l+1}({\bf z})|},$$

     $${\bf w}^{l+1}_{*}({\bf z}) =-\Delta_{h}{\bf u}^{l+1}_{*}({\bf z})\,.$$

2. Stop the iterations once

   $$\max_{{\bf z}\in\mathcal{N}_{h}}|{\bf u}^{l+1}({\bf z})-{\bf u}_{*}^{l}({\bf z})|\leq\epsilon\,.$$

3. Set ${\bf U}^{n+1}={\bf u}_{*}^{l+1}$ and proceed to the next time level.

The Lagrange multiplier $\lambda_{w*}^{l}$ is a counterpart of (3.5) which replaces the unknown value ${\bf U}^{n+1}$, ${\bf W}^{n+1}$ by the respective solutions ${\bf u}^{l}_{*}$, ${\bf w}^{l}_{*}$ which are known from the previous iteration of the nonlinear solver, i.e., we denote $\mathbf{U}^{n+1/2}_{*}=\frac{1}{2}\big{(}{\bf u}_{*}^{l+1}+{\bf U}^{n}\big{)}$, $d_{t}\mathbf{U}^{n+1}_{*}=\frac{1}{\tau}\big{(}{\bf u}_{*}^{l+1}-{\bf U}^{n}\big{)}$, and set

(7.4) $$\lambda_{w*}^{l}({\bf z})=\left\{\begin{array}[]{ll}0&\mbox{if }\mathbf{U}_{*}^{n+1/2}({\bf z})={\bf 0}\,,\\ \frac{\frac{1}{2}\langle\mathbf{V}^{0}({\bf z}),d_{t}\mathbf{U}^{1}_{*}({\bf z})\rangle}{|\mathbf{U}_{*}^{1/2}({\bf z})|^{2}}+\frac{\bigl{(}\nabla\mathbf{W}_{*}^{1/2},\mathbf{U}_{*}^{1/2}({\bf z})\otimes\nabla\varphi_{{\bf z}}\bigr{)}}{\beta_{{\bf z}}|\mathbf{U}_{*}^{1/2}({\bf z})|^{2}}&\text{if }n=0,\,\mathbf{U}_{*}^{1/2}({\bf z})\neq{\bf 0},\\ -\frac{\bigl{\langle}d_{t}\mathbf{U}^{n}({\bf z}),{d_{t}}\mathbf{U}_{*}^{n+1/2}({\bf z})\bigr{\rangle}}{|\mathbf{U}_{*}^{n+1/2}({\bf z})|^{2}}+\frac{\bigl{(}\nabla\mathbf{W}_{*}^{n+1/2},\mathbf{U}_{*}^{n+1/2}({\bf z})\otimes\nabla\varphi_{{\bf z}}\bigr{)}}{\beta_{{\bf z}}|\mathbf{U}_{*}^{n+1/2}({\bf z})|^{2}}&\mbox{else}\,.\end{array}\right.$$

We choose the stopping tolerance $\epsilon=10^{-8}$ in all experiments below.

Remark 7.1 (Projection step).

The use of an additional projection step in the fixed-point iteration was introduced in [10] for the related second-order problems. Without the projection step (step $(1)-(ii)$ in the fixed-point algorithm above) the fixed-point algorithm requires to choose $\tau=\mathcal{O}(h^{2})$ to guarantee convergence. Similarly as in [10], the algorithm with the projection step requires a weaker time-step restriction $\tau=\mathcal{O}(h)$ for convergence. Compared to the second order harmonic wave map equation [10], the bi-harmonic problem requires roughly $100$-times smaller time-step to achieve similar convergence of the fixed-point algorithm, however, we also observe a faster evolution than in the case of second order problems.

7.2. Numerical experiments

We consider the initial condition on $\Omega=\left(-\frac{1}{2},\frac{1}{2}\right)^{2}$

(7.7) $${\bf u}_{0}({\bf x})=\left\{\begin{array}[]{ll}(4{\bf x}A,A^{2}-4|{\bf x}|^{2})/(A^{2}+4|{\bf x}|^{2})&\mbox{for }|{\bf x}|\leq 1/2,\\ (0,0,-1)&\mbox{for }|{\bf x}|\geq 1/2,\\ \end{array}\right.$$

where $A=(1-4|{\bf x}|)^{4}$. In the case of second order problems, i.e., the harmonic heat flow and wave map equations, respectively, the evolution starting from the above initial condition leads to a (discrete) finite-time blow-up, cf., [10] and the references therein.

We set $\mathbf{U}^{0}=\mathcal{I}_{h}{\bf u}_{0}$, $\mathbf{W}^{0}=-\Delta_{h}\mathbf{U}^{0}$ and $\mathbf{V}^{0}=\mathbf{0}$. We compute the problem for different mesh sizes $h=2^{k}$ with corresponding time-step steps $\tau=2^{6-k}\times 10^{-4}$ for $k=5,6,7$. In Figure 1 we display the evolution of the computed energies as well as of the gradient $\|\nabla\mathbf{u}_{h}\|_{L^{\infty}(\Omega)}$. Compared to the case of wave map into spheres, cf. for instance [10, Example 4.3] we do not observe formation of singularities (discrete blow-up), i.e., the gradient of the numerical solution remains bounded independently of the mesh size. Snapshots of the solution computed with $h=2^{-5}$ at different times are displayed in Figure 2 (the figure is colored by the magnitude of the $z$-component of the solution).

Figure 1. Evolution of the discrete energy (left) and of the gradient $\|\nabla\mathbf{u}_{h}\|_{L^{\infty}(\Omega)}$ (right) for $h=2^{k}$, $k=5,6,7$.

Figure 2. Numerical solution computed with $h=2^{-5}$ at time $t=0,0.02,0.045,0.095,0.13,0.24$.