Energy-preserving mixed finite element methods for the Hodge wave equation

Yongke Wu and Yanhong Bai School of Mathematical Sciences, University of Electronic Science and Technology of China, Chengdu 611731, China. wuyongke1982@sina.com School of Science, Xihua University, Chengdu 610039, China. baiyanhong1982@126.com

Abstract.

Energy-preserving numerical methods for solving the Hodge wave equation is developed in this paper. Based on the de Rham complex, the Hodge wave equation can be formulated as a first-order system and mixed finite element methods using finite element exterior calculus is used to discretize the space. A continuous time Galerkin method, which can be viewed as a modification of the Crank-Nicolson method, is used to discretize the time which results in a full discrete method preserving the energy exactly when the source term is vanished. A projection based operator is used to establish the optimal order convergence of the proposed methods. Numerical experiments are present to support the theoretical results.

Key words and phrases:

the Hodge wave equation, energy conservation, de Rham complex, optimal error estimates

2010 Mathematics Subject Classification:

65M60; 65M12; 65J08

Y. Wu was supported by the National Natural Science Foundation of China (11971094 and 11501088). Y. Bai was supported by the National Natural Science Foundation of China (11701481).

1. Introduction

We consider energy-preserving numerical methods for solving the Hodge wave equation, the hyperbolic equation in $\mathbb{R}^{n}$ associated to the Hodge Laplacian of differential $k$ -forms for $0\leq k\leq n$ . The initial-boundary value problem we study is: Find $u:\ (0,T]\mapsto H_{0}\Lambda^{k}(\Omega)$ satisfying

(1)

\displaystyle u_{tt}+(\,{\rm d}\delta+\delta\,{\rm d})u

\displaystyle=f\qquad\text{in }\Omega\times(0,T],

with homogeneous boundary conditions

(2)

\displaystyle\operatorname*{tr}(u)=0,\quad\operatorname*{tr}(\star\,{\rm d}u)

\displaystyle=0\qquad\text{on }\partial\Omega\times(0,T],

and initial conditions

(3)

\displaystyle u(\cdot,0)=u_{0}(\cdot),\quad u_{t}(\cdot,0)

\displaystyle=u_{1}(\cdot)\qquad\text{in }\Omega.

Here $\Omega\subset\mathbb{R}^{n}$ is a domain homomorphism to a ball with piecewise smooth and Lipschitz boundary. The unknown $u$ is a time dependent differential $k$ -form on $\Omega$ , $u_{t}$ and $u_{tt}$ denote its partial derivatives with respect to time variable, and $\,{\rm d}$ , $\delta$ , $\star$ , and $\operatorname*{tr}$ denote exterior derivative, co-derivative, Hodge star, and the trace operator, respectively; see Section 2 for precise definitions. We assume that $T$ is a finite positive real number denoting the ending time.

Many physical problems can be described by (1), such as the mathematical models of sound waves ( $n=3$ and $k=0$ ), electromagnetic waves ( $n=3$ and $k=1$ ), structural vibration ( $n=3$ and $k=2$ ) and so on. There are many theoretical analyses of finite element method for (1) in the special case $n=2$ or $3$ and $k=0$ or $k=n-1$ ; see [15, 4, 17, 11, 16, 28, 19, 21, 9, 22, 20] and the references therein. The pioneer work on mixed finite element methods [5] for the general form of the Hodge wave equation (1) can be found in Quenneville-Bélair’s Ph. D. thesis [26]; see also [3]. In this work, he has presented (1) the abstract Hodge wave equation in the mixed form, (2) the semi-discretization in space for solving the Hodge wave equation (3) the existence and uniqueness of the solution for the semi-discretization in space, (4) the error estimates in the $\|\cdot\|_{L^{\infty}(L^{2})}$ norm for the semi-discretization in space based on the elliptic projection operator.

In the present work, we shall give more thorough analysis of the mixed finite element method developed in [26, 3]. Introduce a $(k-1)$ -form $\sigma=\delta u$ and a $(k+1)$ -form $\omega=\,{\rm d}u$ with standard modification for $k=0$ or $k=n$ , and a $k$ -form $\mu=u_{t}$ . The first order formulation of (1) reads as: find $\sigma\in H_{0}\Lambda^{-}$ , $\mu\in H_{0}\Lambda$ , and $\omega\in H_{0}\Lambda^{+}$ such that

(4)	$\displaystyle\langle\sigma_{t},\tau\rangle-\langle\,{\rm d}^{-}\tau,\mu\rangle$	$\displaystyle=0\qquad\qquad\forall~{}~{}\tau\in H_{0}\Lambda^{-},$
(5)	$\displaystyle\langle\mu_{t},v\rangle+\langle\,{\rm d}^{-}\sigma,v\rangle+\langle\omega,\,{\rm d}v\rangle$	$\displaystyle=\langle f,v\rangle\,\qquad\forall~{}~{}v\in H_{0}\Lambda,$
(6)	$\displaystyle\langle\omega_{t},\phi\rangle-\langle\,{\rm d}\mu,\phi\rangle$	$\displaystyle=0\qquad\qquad\forall~{}~{}\phi\in H_{0}\Lambda^{+},$

with initial conditions

\sigma_{0}=\delta u_{0},\quad\mu_{0}=u_{1},\quad\omega_{0}=\,{\rm d}u_{0}.

Comparing with [26], the main contributions of this paper are as follows. Firstly, we use the skew-symmetric property of the formulation (4)-(6) to get the following energy estimates

(7)		$\displaystyle\sup\limits_{0\leq t\leq T}E(t)$	$\displaystyle\leq E(0)+2\int_{0}^{T}\\|f(\cdot,s)\\|\,{\rm d}s,$
(8)		$\displaystyle\sup\limits_{0\leq t\leq T}H(t)$	$\displaystyle\leq H(0)+4\\|f\\|_{L^{\infty}(L^{2})}+2\int_{0}^{T}\\|f_{t}(\cdot,s)\\|\,{\rm d}s,$

with

	$\displaystyle E(t)$	$\displaystyle=(\\|\sigma(\cdot,t)\\|^{2}+\\|\mu(\cdot,t)\\|^{2}+\\|\omega(\cdot,t)\\|^{2})^{1/2},$
	$\displaystyle H(t)$	$\displaystyle=(\\|\,{\rm d}^{-}\sigma(\cdot,t)\\|^{2}+\\|\,{\rm d}\mu(\cdot,t)\\|^{2}+\\|\delta\mu(\cdot,t)\\|^{2}+\\|\delta^{+}\omega(\cdot,t)\\|^{2})^{1/2}.$

These energy estimates imply the existence and uniqueness of solution for (4)-(6); see Remark 2.4. When (4)-(6) is self-conserve, i.e., $f=0$ , the inequality (7)-(8) become equalities which implies the energies $E$ and $H$ are preserved exactly; see Remark 2.3. Due to the structure preserving properties of the finite element exterior calculus (FEEC) [1, 2, 3], the semi-discretization in space also inherit the skew-symmetric property of the spatial differential terms, and thus the energy conservation is preserved naturally. We then use the continuous time Galerkin method [16] to give unconditioned energy conservation schemes. Here we follow the approach in [17, 16, 22], where the energy estimates has been derived for scalar wave equations but not for Hodge wave equations. As we know, energy conservation numerical schemes can have a crucial influence on the quality of the numerical simulations. Especially, in long-time simulations, energy-preserving can have a dramatic effect on stability and global error growth.

Secondly, we obtain the optimal convergence order of the error estimates in both $L^{2}$ -norm and $\|\mathcal{A}(\cdot)\|_{L^{\infty}(L^{2})}$ -norm for the semi- and full-discrete mixed finite element methods, where $\mathcal{A}$ is a skew-symmetric operator defined in Section 2. Such result has been derived for scalar wave equation [17, 16, 22] but generalization to general Hodge wave equation is non-trivial. Technically, the canonical interpolation operators $\pi_{h}$ used in [17, 16, 22] cannot be commutated with the discrete co-derivative operator $\delta_{h}$ , and the $L^{2}$ projection operator $Q_{h}$ cannot be commutated with the exterior derivative operators $\,{\rm d}$ . Using these standard operators in the convergence analysis will lead to the lost of the convergence order. To overcome this difficulty, we choose a projection based interpolation operator $I_{h}$ briefly mentioned in [8, Proposition 5.44] and redefine it based on the Hodge decomposition. Such projection based operators has been introduced for $H^{1},H({\rm curl\,})$ and $H(\operatorname{div})$ spaces in [12, 13, 14, 25], where the authors have proved that these projection based operators made the de Rahm diagram commute and had the quasi-optimal interpolation error bound for $hp$ finite element spaces. Although this projection-based quasi-interpolation operator is not new, the properties we are going to prove are not fully explored in the literature. Specifically, we shall prove that (1) $I_{h}$ is commuted with $\delta_{h}$ , (2) $I_{h}$ is stable in both $\|\,{\rm d}(\cdot)\|$ and $\|\delta_{h}(\cdot)\|$ norms, (3) $I_{h}$ is the $L^{2}$ orthogonal projection to the space $\mathfrak{Z}_{0,h}$ , (4) $I_{h}$ is an orthogonal projection operator with respect to the inner-product $\langle\,{\rm d}(\cdot),\,{\rm d}(\cdot)\rangle$ , (5) $I_{h}$ has the same approximation properties as the classical interpolation operators; see Lemma 3.3 - 3.5. By using these properties of the projection-based operator $I_{h}$ , we get the optimal error estimates for both the semi- and full- discretization with respect to both $\|\cdot\|_{L^{\infty}(L^{2})}$ and $\|\mathcal{A}(\cdot)\|_{L^{\infty}(L^{2})}$ norms (the detail definition of these norms can be found in Section 3), while recall that [26] only give the error estimates of $\|\cdot\|_{L^{\infty}(L^{2})}$ for the semi-discretization in space and as the line of Quenneville’s proof, it seems difficulty to get the error estimate of the energy norm $\|\mathcal{A}(\cdot)\|_{L^{\infty}(L^{2})}$ . But the control of the energy norm is very important, since the $L^{2}$ -norm is possible small but the energy norm is larger due to the small oscillation in the error. Furthermore our error estimate, comparing with [26] is robust to $T$ in the sense that the factor $T$ is absent on the error estimates; see Theorem 3.9, 3.13, 4.6 and 4.8. Such error estimates imply that our algorithms are robust for long time problems and the numerical experiment supports this result; see Table 3.

What remains of this paper is organized as follow. In Section 2 we introduce the required background on finite element exterior calculus (FEEC) and the Hodge wave equation. We obtain the mixed formulation of the Hodge wave equation and get the energy conservation estimates. Section 3, we briefly introduce the finite element spaces on $k$ -forms, give the semi-discrete form of the Hodge wave equation, introduce a projection-based quasi-interpolation operator and explore properties of this operator, obtain the energy estimates of the semi-discrete form, and get the optimal error estimates of the semi-discrete form. In section 4, the full-discrete form of the Hodge wave equation is obtained, the energy estimates and the optimal error estimates are obtained. Section 5 give some numerical experiments to confirm our theoretical results.

Throughout this paper, $i$ , $h$ and $\,\Delta t$ denote the time level, the mesh size and the time step size, respectively. The capital $C$ may be different in different places, denotes a positive constant which is independent on $i$ , $h$ and $\,\Delta t$ . We denote by $\|\cdot\|_{m,p}$ the norm of the classical Sobolev spaces $W^{m,p}\Lambda^{k}(\Omega)$ , $1\leq p\leq\infty$ and $0\leq k\leq n$ . If $p=2$ , we write $\|\cdot\|_{m,p}$ simply as $\|\cdot\|_{m}$ and denote by $|\cdot|_{m}$ the semi-norm in $W^{m,2}\Lambda^{k}(\Omega)$ . In addition, for any Sobolev space $Y$ , we define the space $L^{p}([a,b],Y)$ with norm $\|f\|_{L^{p}(Y)}=\left(\int_{a}^{b}\|f(\cdot,t)\|_{Y}^{p}\,{\rm d}t\right)^{1/p}$ , and if $p=\infty$ , the integral is replaced by the essential supremum.

2. Preliminaries

In this section, we follow the convention of [1, 2, 3] to introduce necessary background of finite element exterior calculus. Then, we introduce the Hodge wave equation and its mixed formulation. Finally, we get the energy conservation estimates for this mixed form.

2.1. de Rham complex

Let $\Omega\subset\mathbb{R}^{n}$ ( $n\geq 2$ ) be a bounded Lipschitz domain. For a given integer $0\leq k\leq n$ , $\Lambda^{k}(\Omega)$ represents the linear space of all smooth $k$ -forms on $\Omega$ . For any $\omega\in\Lambda^{k}(\Omega)$ , $\omega$ can be written as

\omega=\sum\limits_{1\leq\sigma_{1}<\cdots<\sigma_{k}\leq n}a_{\sigma}\,{\rm d}x^{\sigma_{1}}\wedge\cdots\wedge\,{\rm d}x^{\sigma_{k}},

with $a_{\sigma}\in C^{\infty}(\Omega)$ and $\wedge$ the wedge product. As $\Omega$ is a flat domain in $\mathbb{R}^{n}$ , we can identify each tangent space of $\Omega$ with $\mathbb{R}^{n}$ . Given an $\omega\in\Lambda^{k}(\Omega)$ and vectors $v_{1},~{}v_{2},\cdots,~{}v_{k}\in\mathbb{R}^{n}$ , we have that the map $\boldsymbol{x}\in\Omega\mapsto\omega_{\boldsymbol{x}}(v_{1},v_{2},\cdots,v_{k})\in\mathbb{R}$ is a smooth map (infinitely differentiable).

We define the exterior derivative $\,{\rm d}^{k}:~{}\Lambda^{k}(\Omega)\rightarrow\Lambda^{k+1}(\Omega)$ as

\,{\rm d}^{k}\omega_{x}(v_{1},v_{2},\cdots,v_{k+1})=\sum\limits_{j=1}^{k+1}(-1)^{j+1}\partial_{v_{j}}\omega_{x}(v_{1},\cdots,\hat{v}_{j},\cdots,v_{k+1}),

where the hat is used to indicate a suppressed argument. By the definition of $\,{\rm d}^{k}$ , it is easy to see that $\,{\rm d}^{k}$ is a sequence of differential operators satisfying that the range of $\,{\rm d}^{k}$ lies in the domain of $\,{\rm d}^{k+1}$ , i.e., $\,{\rm d}^{k+1}\circ\,{\rm d}^{k}=0$ for $k=0,1,\cdots,n-1$ . For convenience of notation, we shall skip the superscript $k$ if there is no confusion.

Let vol be the unique volume form in $\Lambda^{k}(\Omega)$ , define the $L^{2}$ -inner product of any two differential $k$ -forms on $\Omega$ as the integral of their pointwise inner product:

\langle\omega,\mu\rangle=\int_{\Omega}\langle\omega_{x},\mu_{x}\rangle\text{vol}.

The completion of $\Lambda^{k}(\Omega)$ under the corresponding norm defines the Hilbert space $L^{2}\Lambda^{k}(\Omega)$ . The domain of the exterior derivative $\,{\rm d}^{k}$ can be enlarged to

H\Lambda^{k}(\Omega)=\{\omega\in L^{2}\Lambda^{k}(\Omega):~{}\,{\rm d}\omega\in L^{2}\Lambda^{k+1}(\Omega)\}.

$H\Lambda^{k}(\Omega)$ is a Hilbert space with inner product $\langle\omega,\mu\rangle+\langle\,{\rm d}\omega,\,{\rm d}\mu\rangle$ and associated graph norm $\|\cdot\|_{H\Lambda}$ . The de Rham complex

(9)

\begin{CD}H\Lambda^{0}(\Omega)@>{\,{\rm d}}>{}>H\Lambda^{1}(\Omega)@>{\,{\rm d}}>{}>\cdots @>{\,{\rm d}}>{}>H\Lambda^{n-1}(\Omega)@>{\,{\rm d}}>{}>H\Lambda^{n}(\Omega)\end{CD}

is then bounded in the sense that $\,{\rm d}:~{}H\Lambda^{k}(\Omega)\rightarrow H\Lambda^{k+1}(\Omega)$ is a bounded operator.

For any smooth manifold $M$ and any $x\in M$ , we use $T_{x}M$ to denote the tangential space of $M$ at $x$ . For any smooth $k$ -form $\omega\in\Lambda^{k}(\Omega)$ , we define $\operatorname*{tr}\omega\in\Lambda^{k}(\partial M)$ as

\operatorname*{tr}\omega(v_{1},v_{2},\cdots,v_{k})=\omega(v_{1},v_{2},\cdots,v_{k})

for tangential vectors $v_{i}\in T_{x}\partial M\subset T_{x}M$ ( $i=1,2,\cdots,k$ ). This operator can be extended continuous to Lipschitz domain $\Omega$ , also denote by $\operatorname*{tr}:~{}H^{1}\Lambda^{k}(\Omega)\rightarrow H^{1/2}\Lambda^{k}(\partial\Omega)$ and $\operatorname*{tr}:~{}H\Lambda^{k}(\Omega)\rightarrow H^{-1/2}\Lambda^{k}(\partial\Omega)$ . Define

	$\displaystyle H_{0}\Lambda^{k}(\Omega)$	$\displaystyle=\{\omega\in H\Lambda^{k}(\Omega):~{}\operatorname*{tr}\omega=0\text{ on }\partial\Omega\},$
	$\displaystyle H_{0}^{1}\Lambda^{k}(\Omega)$	$\displaystyle=\{\omega\in H^{1}\Lambda^{k}(\Omega):~{}\operatorname*{tr}\omega=0\text{ on }\partial\Omega\}.$

In the following sections, we will focus on the de Rham complex with homogeneous trace

(10)

\begin{CD}H_{0}\Lambda^{0}(\Omega)@>{\,{\rm d}}>{}>H_{0}\Lambda^{1}(\Omega)@>{\,{\rm d}}>{}>\cdots @>{\,{\rm d}}>{}>H_{0}\Lambda^{n-1}(\Omega)@>{\,{\rm d}}>{}>H_{0}\Lambda^{n}(\Omega)\end{CD}

In order to define the dual complex, we start with the Hodge star operator $\star:~{}\Lambda^{k}(\Omega)\rightarrow\Lambda^{n-k}(\Omega)$ ,

\int_{\Omega}\omega\wedge\mu=\langle\star\omega,\mu\rangle,\qquad\forall~{}~{}\omega\in\Lambda^{k}(\Omega),~{}~{}\mu\in\Lambda^{n-k}(\Omega).

The coderivative operator $\delta^{k}:~{}\Lambda^{k}(\Omega)\rightarrow\Lambda^{k-1}(\Omega)$ is defined as

\delta^{k}\omega=(-1)^{k(n-k+1)}\star\,{\rm d}^{n-k}\star\omega.

$\,{\rm d}^{k-1}$ and $\delta^{k}$ are related by the Stokes theorem

\langle\,{\rm d}\omega,\mu\rangle=\langle\omega,\delta\mu\rangle+\int_{\partial\Omega}\operatorname*{tr}\omega\wedge\operatorname*{tr}(\star\mu),\qquad\omega\in\Lambda^{k-1}(\Omega),~{}~{}\mu\in\Lambda^{k}(\Omega).

We define the spaces

	$\displaystyle H^{*}\Lambda^{k}(\Omega)$	$\displaystyle=\{\omega\in L^{2}\Lambda^{k}(\Omega):~{}\delta\omega\in L^{2}\Lambda^{k-1}(\Omega)\},$
	$\displaystyle H_{0}^{*}\Lambda^{k}(\Omega)$	$\displaystyle=\{\omega\in H^{}\Lambda^{k}(\Omega):~{}\operatorname{tr}\star\omega=0\text{ on }\partial\Omega\}.$

Treat $\,{\rm d}:~{}H_{0}\Lambda^{k}(\Omega)\subset L^{2}\Lambda^{k}(\Omega)\rightarrow L^{2}\Lambda^{k+1}(\Omega)$ as an unbounded and densely defined operator. Then Stokes theorem implies that $\delta:~{}H^{*}\Lambda^{k+1}(\Omega)\subset L^{2}\Lambda^{k+1}(\Omega)\rightarrow L^{2}\Lambda^{k}(\Omega)$ is the adjoint of $\,{\rm d}$ as

(11)

\langle\,{\rm d}\omega,\mu\rangle=\langle\omega,\delta\mu\rangle,\qquad\forall~{}~{}\omega\in H_{0}\Lambda^{k}(\Omega),~{}~{}\mu\in H^{*}\Lambda^{k+1}(\Omega).

We have a dual sequence of (10)

(12)

\begin{CD}H^{*}\Lambda^{0}(\Omega)@<{\delta}<{}<H^{*}\Lambda^{1}(\Omega)@<{\delta}<{}<\cdots @<{{\delta}}<{}<H^{*}\Lambda^{n-1}(\Omega)@<{\delta}<{}<H^{*}\Lambda^{n}(\Omega).\end{CD}

Let $\mathfrak{Z}_{0}^{k}$ be the kernel of $\,{\rm d}$ in the space $H_{0}\Lambda^{k}(\Omega)$ , then $\mathfrak{Z}_{0}^{k}$ can be decomposed as $\mathfrak{Z}_{0}^{k}=\mathfrak{B}_{0}^{k}\oplus^{\bot_{L^{2}}}\mathfrak{H}_{0}^{k}$ , where $\mathfrak{B}_{0}^{k}$ is the range of $\,{\rm d}^{k-1}$ , i.e., $\mathfrak{B}_{0}^{k}=\,{\rm d}(H_{0}\Lambda^{k-1}(\Omega))$ and $\mathfrak{H}_{0}^{k}$ is the space of harmonic forms, i.e., $\mathfrak{H}_{0}^{k}=\{\omega\in H_{0}\Lambda^{k}(\Omega)\cap H^{*}\Lambda^{k}(\Omega):~{}~{}\,{\rm d}\omega=0\text{ and }\delta\omega=0\}$ , $\oplus^{\bot_{L^{2}}}$ means that the decomposition is orthogonal in the sense of the $L^{2}$ -inner product. The following Hodge decomposition has been established in [1, page 22]:

L^{2}\Lambda^{k}(\Omega)=\mathfrak{B}_{0}^{k}\oplus^{\bot_{L^{2}}}\mathfrak{H}_{0}^{k}\oplus^{\bot_{L^{2}}}\delta H^{*}\Lambda^{k+1}(\Omega).

Denote $\mathfrak{K}^{k}$ as the $L^{2}$ orthogonal complement of $\mathfrak{Z}_{0}^{k}$ in $H_{0}\Lambda^{k}(\Omega)$ , i.e., $\mathfrak{K}^{k}=H_{0}\Lambda^{k}(\Omega)\cap\delta H^{*}\Lambda^{k+1}(\Omega)$ . Then we have the Hodge decomposition of $H_{0}\Lambda^{k}(\Omega)$ :

(13)

H_{0}\Lambda^{k}(\Omega)=\mathfrak{Z}_{0}^{k}\oplus^{\bot_{L^{2}}}\mathfrak{K}^{k}=\mathfrak{B}_{0}^{k}\oplus^{\bot_{L^{2}}}\mathfrak{H}_{0}^{k}\oplus^{\bot_{L^{2}}}\mathfrak{K}^{k}.

It should be point out that when $k=0$ , we have $\mathfrak{Z}_{0}^{0}=\{0\}$ and $\mathfrak{K}^{-1}=\{0\}$ . When $k=n$ , we have $\mathfrak{K}^{n}=\{0\}$ .

In the following sections, when spaces of the consecutive differential forms are involved, we use the short sequences

(14)

\begin{CD}H_{0}\Lambda^{-}(\Omega)@>{{\,{\rm d}^{-}}}>{}>H_{0}\Lambda(\Omega)@>{\,{\rm d}}>{}>H_{0}\Lambda^{+}(\Omega)\end{CD}

or the one with the Hodge decomposition

(15)

\begin{CD}\mathfrak{B}_{0}^{-}\oplus^{\bot_{L^{2}}}\mathfrak{H}_{0}^{-}\oplus^{\bot_{L^{2}}}\mathfrak{K}^{-}@>{{\,{\rm d}^{-}}}>{}>\mathfrak{B}_{0}\oplus^{\bot_{L^{2}}}\mathfrak{H}_{0}\oplus^{\bot_{L^{2}}}\mathfrak{K}@>{{\,{\rm d}}}>{}>\mathfrak{B}_{0}^{+}\oplus^{\bot_{L^{2}}}\mathfrak{H}_{0}^{+}\oplus^{\bot_{L^{2}}}\mathfrak{K}^{+}.\end{CD}

In this paper, we consider the domain $\Omega$ with zero Betti numbers, namely, we impose the following assumption on the domain $\Omega$ :

(A):: We assume that $\Omega$ is simple in the sense that $\dim\mathfrak{H}_{0}^{k}=0$ for all $1\leq k\leq n-1$ .

2.2. The Hodge wave equation

The Hodge wave equation reads as given $f:~{}(0,T)\mapsto L^{2}\Lambda^{k}$ , find $u\in H^{2}((0,T),H_{0}\Lambda^{k}\cap H^{*}\Lambda^{k})$ such that

(16)

u_{tt}+\mathcal{L}u=f\qquad\text{in }\Omega,

where $\mathcal{L}=\,{\rm d}^{-}\delta+\delta^{+}\,{\rm d}$ is called the Hodge Laplacian operator [2], with the initial conditions

(17)

u(\cdot,0)=u_{0}(\cdot),\qquad u_{t}(\cdot,0)=u_{1}(\cdot).

For easy to preserve the energy exactly, we will use mixed method to discrete (16). Introduce a $(k-1)$ -form $\sigma=\delta u$ and a $(k+1)$ -form $\omega=\,{\rm d}u$ with standard modification for $k=0$ or $k=n$ , and a $k$ -form $\mu=u_{t}$ . The mixed formulation [26] of the Hodge wave equation (16) is: given $f\in L^{2}((0,T),L^{2}\Lambda)$ , find $(\sigma,\mu,\omega):~{}(0,T]\mapsto H_{0}\Lambda^{-}\times H_{0}\Lambda\times H_{0}\Lambda^{+}:=\boldsymbol{W}$ such that

(18)	$\displaystyle\langle\sigma_{t},\tau\rangle-\langle\,{\rm d}^{-}\tau,\mu\rangle$	$\displaystyle=0\qquad\qquad\forall~{}~{}\tau\in H_{0}\Lambda^{-},$
(19)	$\displaystyle\langle\mu_{t},v\rangle+\langle\,{\rm d}^{-}\sigma,v\rangle+\langle\omega,\,{\rm d}v\rangle$	$\displaystyle=\langle f,v\rangle\,\qquad\forall~{}~{}v\in H_{0}\Lambda^{k},$
(20)	$\displaystyle\langle\omega_{t},\phi\rangle-\langle\,{\rm d}\mu,\phi\rangle$	$\displaystyle=0\qquad\qquad\forall~{}~{}\phi\in H_{0}\Lambda^{+},$

with initial conditions

\sigma(\cdot,0)=\delta u_{0},\quad\mu(\cdot,0)=u_{1}(\cdot),\quad\omega(\cdot,0)=\,{\rm d}u_{0}(\cdot).

Denoted by

\mathcal{A}=\begin{pmatrix}0&\delta&0\\ -\,{\rm d}^{-}&0&-\delta^{+}\\ 0&\,{\rm d}&0\end{pmatrix}.

The existence of solutions for the mixed formulation (18)-(20) can be found in [26] and it can also be obtained by Picard Theorem since the operator

\mathcal{A}:~{}H_{0}\Lambda^{-}\times(H_{0}\Lambda\cap H^{*}\Lambda)\times(H_{0}\Lambda^{+}\cap H^{*}\Lambda^{+})\rightarrow L^{2}\Lambda^{-}\times L^{2}\Lambda\times L^{2}\Lambda^{+}

is bounded. To prove the uniqueness of the solution, we need the energy estimates. We introduce a basic inequality.

Lemma 2.1.

([22, Lemma 1]) Suppose that a real number $x$ satisfies the quadratic inequality

x^{2}\leq\gamma^{2}+\beta x

for $\beta,~{}\gamma\geq 0$ and $\beta^{2}+\gamma^{2}>0$ . Then

x\leq\beta+\gamma.

We define two energies of the mixed formulation (18)-(20) as

E(t)=\left(\|\sigma(\cdot,t)\|^{2}+\|\mu(\cdot,t)\|^{2}+\|\omega(\cdot,t)\|^{2}\right)^{1/2}

and

H(t)=\left(\|\,{\rm d}^{-}\sigma(\cdot,t)\|^{2}+\|\,{\rm d}\mu(\cdot,t)\|^{2}+\|\delta\mu(\cdot,t)\|^{2}+\|\delta^{+}\omega(\cdot,t)\|^{2}\right)^{1/2}.

We have the following energy estimates.

Theorem 2.2.

Let $\boldsymbol{u}=(\sigma,\mu,\omega)^{\intercal}\in\boldsymbol{W}$ be the solution of the mixed formulation (18)-(20). Provided $f\in L^{1}((0,T),L^{2}\Lambda)$ , we have the energy bound

(21)

\sup\limits_{0\leq s\leq T}E(t)\leq E(0)+2\int_{0}^{T}\|f(\cdot,s)\|\,{\rm d}s.

Furthermore, if $f\in W^{1,1}((0,T),L^{2}\Lambda)$ , we have the bound

(22)

\sup\limits_{0\leq t\leq T}H(t)\leq H(0)+4\|f\|_{L^{\infty}(L^{2})}+2\int_{0}^{T}\|f_{t}(\cdot,s)\|\,{\rm d}s.

When $f=0$ , the inequalities become equalities and thus we have the energy conservation

E(t)=E(0),H(t)=H(0),\quad\forall t>0.

Proof.

Taking $\tau=\sigma$ , $v=\mu$ and $\phi=\omega$ in (18) - (20) and adding them together, we obtain

\frac{1}{2}\frac{\,{\rm d}}{\,{\rm d}t}E^{2}(t)=\langle f,\mu\rangle.

Integrate the above equation on the interval $(0,s)$ , for any $s\in(0,T]$ , we have

	$\displaystyle E^{2}(s)$	$\displaystyle=E^{2}(0)+2\int_{0}^{s}\langle f,\mu\rangle\,{\rm d}t$
		$\displaystyle\leq E^{2}(0)+2\sup\limits_{0\leq t\leq T}E(t)\int_{0}^{T}\\|f\\|\,{\rm d}t.$

Then (21) follows by Lemma 2.1.

Taking $\tau=\delta\mu_{t}$ , $v=-\,{\rm d}^{-}\sigma_{t}-\delta^{+}\omega_{t}$ and $\phi=\,{\rm d}\mu_{t}$ in (18)-(20), we have

	$\displaystyle\langle\,{\rm d}^{-}\sigma_{t},\mu_{t}\rangle-\frac{1}{2}\frac{\,{\rm d}}{\,{\rm d}t}\\|\delta\mu\\|^{2}$	$\displaystyle=0,$
	$\displaystyle-\langle\mu_{t},\,{\rm d}^{-}\sigma_{t}\rangle-\langle\,{\rm d}\mu_{t},\omega_{t}\rangle-\frac{1}{2}\frac{\,{\rm d}}{\,{\rm d}t}\\|\,{\rm d}^{-}\sigma\\|^{2}-\frac{1}{2}\frac{\,{\rm d}}{\,{\rm d}t}\\|\delta^{+}\omega\\|^{2}$	$\displaystyle=-\langle f,\,{\rm d}^{-}\sigma_{t}+\delta^{+}\omega_{t}\rangle,$
	$\displaystyle\langle\omega_{t},\,{\rm d}\mu_{t}\rangle-\frac{1}{2}\frac{\,{\rm d}}{\,{\rm d}t}\\|\,{\rm d}\mu\\|^{2}$	$\displaystyle=0.$

Add the above equations together, we obtain

\frac{1}{2}\frac{\,{\rm d}}{\,{\rm d}t}H^{2}(t)=\langle f,\,{\rm d}^{-}\sigma_{t}+\delta^{+}\omega_{t}\rangle.

Pick any $0\leq s\leq T$ and integrate from $0$ to $s$ to obtain

	$\displaystyle H^{2}(s)$	$\displaystyle=H^{2}(0)+2\int_{0}^{s}\langle f,\,{\rm d}^{-}\sigma_{t}+\delta^{+}\omega_{t}\rangle\,{\rm d}t$
		$\displaystyle=H^{2}(0)+2\langle f(\cdot,s),\,{\rm d}^{-}\sigma(\cdot,s)+\delta^{+}\omega(\cdot,s)\rangle-2\langle f(\cdot,0),\,{\rm d}^{-}\sigma(\cdot,0)+\delta^{+}\omega(\cdot,0)\rangle$
		$\displaystyle\quad-2\int_{0}^{s}\langle f_{t},\,{\rm d}^{-}\sigma+\delta^{+}\omega\rangle\,{\rm d}t$
		$\displaystyle\leq H^{2}(0)+\sup\limits_{0\leq t\leq T}H(t)\left(4\\|f\\|_{L^{\infty}(L^{2})}+2\int_{0}^{T}\\|f_{t}\\|\,{\rm d}t\right).$

Then the desired inequality (22) follows from Lemma 2.1. ∎

Remark 2.3.

When the source term $f$ of (16) equal $0$ , i.e., (16) is a self-conserve system, Theorem 2.2 implies that the mixed form (18)-(20) preserves the energies $E$ and $H$ exactly.

Remark 2.4.

Theorem 2.2 implies the uniqueness of solution of (18)-(20) in the space $\boldsymbol{W}$ . Together with the existence of solutions in the space $H_{0}\Lambda^{-}\times(H_{0}\Lambda\cap H^{*}\Lambda)\times(H_{0}\Lambda^{+}\cap H^{*}\Lambda^{+})$ , we obtain that (18)-(20) have a unique solution $\boldsymbol{u}=(\sigma,\mu,\omega)^{\intercal}$ in the space $H_{0}\Lambda^{-}\times(H_{0}\Lambda\cap H^{*}\Lambda)\times(H_{0}\Lambda^{+}\cap H^{*}\Lambda^{+})$ and for any $t\in(0,T]$ satisfying

(23)

\langle\boldsymbol{u}_{t},\boldsymbol{v}\rangle+\langle\mathcal{A}\boldsymbol{u},\boldsymbol{v}\rangle=\langle\boldsymbol{F},\boldsymbol{v}\rangle\qquad\forall~{}~{}\boldsymbol{v}\in\boldsymbol{W},

with $\boldsymbol{F}=(0,~{}f,~{}0)^{\intercal}$ .

3. Semi-discretization of the Hodge wave equation

In this section, we will introduce mixed finite element methods developed in [26, 3] for the spatial discretization of the Hodge wave equation (16), and give the energy estimates and optimal error estimates.

3.1. Finite element spaces

Let $\mathcal{T}_{h}$ be a shape regular triangulation of $\Omega$ . For each $n$ -simplex $K\in\mathcal{T}_{h}$ , we define $h_{K}=|K|^{1/n}$ and $h=\max\limits_{K\in\mathcal{T}_{h}}h_{K}$ . For completeness, we briefly introduce the construction of finite element spaces following [1, 3].

Denote $\mathcal{P}_{r}(\mathbb{R}^{n})$ as the space of polynomials in $n$ variables of degree at most $r$ and $\mathcal{H}_{r}(\mathbb{R}^{n})$ as the space of homogeneous polynomial functions of degree $r$ . Spaces of polynomial differential forms $\mathcal{P}_{r}\Lambda^{k}(\mathbb{R}^{n})$ and $\mathcal{H}_{r}(\mathbb{R}^{n})$ can be defined by using the corresponding polynomial as the coefficients. We will suppress $\mathbb{R}^{n}$ from the notation for simplicity. For each integer $r\geq n$ , we have the polynomial subcomplex of the de Rham complex

\setcounter{MaxMatrixCols}{11}\begin{CD}0@>{}>{}>\mathcal{P}_{r}\Lambda^{0}@>{\,{\rm d}}>{}>\mathcal{P}_{r-1}\Lambda^{1}@>{\,{\rm d}}>{}>\cdots @>{\,{\rm d}}>{}>\mathcal{P}_{r-n}\Lambda^{n}@>{}>{}>0.\end{CD}

Given a point $x\in\mathbb{R}^{n}$ , treat $x$ as a vector in the tangential space $T_{x}\mathbb{R}^{n}$ and define the Koszul operator $\kappa:~{}\Lambda^{k}(\mathbb{R}^{n})\rightarrow\Lambda^{k-1}(\mathbb{R}^{n})$ as

(\kappa\omega)_{x}(v_{1},v_{2},\cdots,v_{k-1})=\omega_{x}(x,v_{1},v_{2},\cdots,v_{k-1}).

This $\kappa$ satisfying the identity $\kappa\,{\rm d}+\,{\rm d}\kappa=(k+r){\rm id}$ [1, Theorem 3.1] on the space $\mathcal{H}_{r}\Lambda^{k}$ and there is a direct sum

\mathcal{H}_{r}\Lambda^{k}=\kappa\mathcal{H}_{r-1}\Lambda^{k+1}\oplus\,{\rm d}\mathcal{H}_{r+1}\Lambda^{k-1}.

Based on the decomposition, the incomplete polynomial differential form can be introduced as

\mathcal{P}_{r}^{-}\Lambda^{k}=\mathcal{P}_{r-1}\Lambda^{k}+\kappa\mathcal{H}_{r-1}\Lambda^{k+1}

and, for $r\geq 1$ , have the following subcomplex of the de Rham complex

\setcounter{MaxMatrixCols}{11}\begin{CD}0@>{}>{}>\mathcal{P}_{r}^{-}\Lambda^{0}@>{\,{\rm d}}>{}>\mathcal{P}_{r}^{-}\Lambda^{1}@>{\,{\rm d}}>{}>\cdots @>{\,{\rm d}}>{}>\mathcal{P}_{r}^{-}\Lambda^{n}@>{}>{}>0.\end{CD}

For each simplex $K\in\mathcal{T}_{h}$ , denote $\mathcal{P}_{r}\Lambda^{k}(K)$ or $\mathcal{P}_{r}^{-}\Lambda^{k}(K)$ as the spaces of $k$ forms obtained by restricting the forms $\mathcal{P}_{r}\Lambda^{k}(\mathbb{R}^{n})$ or $\mathcal{P}_{r}^{-}\Lambda^{k}(\mathbb{R}^{n})$ , respectively, to $K$ . We then obtain the finite element spaces

	$\displaystyle\mathcal{P}_{r}\Lambda^{k}(\mathcal{T}_{h})$	$\displaystyle=\{\omega\in H\Lambda^{k}(\Omega):~{}~{}\omega\|_{K}\in\mathcal{P}_{r}\Lambda^{k}(K),~{}~{}\forall~{}~{}K\in\mathcal{T}_{h}\},$
	$\displaystyle\mathcal{P}_{r}^{-}\Lambda^{k}(\mathcal{T}_{h})$	$\displaystyle=\{\omega\in H\Lambda^{k}(\Omega):~{}~{}\omega\|_{K}\in\mathcal{P}_{r}^{-}\Lambda^{k}(K),~{}~{}\forall~{}~{}K\in\mathcal{T}_{h}\}.$

We choose $V_{h}^{k}=\mathcal{P}_{r}\Lambda^{k}(\mathcal{T}_{h})\cap H_{0}\Lambda^{k}$ or $V_{h}^{k}=\mathcal{P}_{r}^{-}\Lambda^{k}(\mathcal{T}_{h})\cap H_{0}\Lambda^{k}$ so that $(V_{h}^{k},\,{\rm d})$ forms a subcomplex of $(H_{0}\Lambda^{k}(\Omega),\,{\rm d})$ . For the consecutive spaces, we shall use short sequence

\begin{CD}V_{h}^{-}@>{\,{\rm d}^{-}}>{}>V_{h}@>{\,{\rm d}}>{}>V_{h}^{+}.\end{CD}

The discrete coderivative $\delta_{h}:~{}V_{h}\rightarrow V_{h}^{-}$ is defined as the $L^{2}$ -adjoint of $\,{\rm d}^{-}:~{}V_{h}^{-}\rightarrow V_{h}$ , i.e., for any given $\omega_{h}\in V_{h}$ , $\delta_{h}\omega_{h}$ is the unique element in $V_{h}^{-}$ such that

(24)

\langle\delta_{h}\omega_{h},v_{h}\rangle=\langle\omega_{h},\,{\rm d}^{-}v_{h}\rangle\qquad\forall~{}~{}v_{h}\in V_{h}^{-}.

The discrete Hodge decomposition of $V_{h}^{k}$ is

(25)

V_{h}^{k}=\mathfrak{Z}_{0,h}\oplus^{\bot_{L^{2}}}\mathfrak{K}_{h},

where $\mathfrak{Z}_{0,h}=\ker(\,{\rm d})\cap V_{h}=\,{\rm d}^{-}\boldsymbol{V}_{h}^{-}\subset\mathfrak{Z}_{0}$ and $\mathfrak{K}_{h}=\delta_{h}^{+}V_{h}^{+}$ is the $L^{2}$ orthogonal complement of $\mathfrak{Z}_{h,0}$ in $V_{h}$ . Generally $\mathfrak{K}_{h}\not\subset\mathfrak{K}$ , since $\delta_{h}$ is not a conforming discretization of $\delta$ . It should be point out that when $k=0$ , we have $\mathfrak{Z}_{0,h}=\{0\}$ and $\mathfrak{K}_{h}^{-}=\{0\}$ . When $k=n$ , we have $\mathfrak{K}_{h}=\{0\}$ .

We have the following discrete Poincaré inequality; cf. [1, Theorem 5.11]

Lemma 3.1 (discrete Poincaré inequality for $\,{\rm d}$ ).

There is a positive constant $C_{p}$ , independent of $h$ , such that

(26)

\|\omega_{h}\|\leq C_{p}\|\,{\rm d}\omega_{h}\|\qquad\forall~{}~{}\omega_{h}\in\mathfrak{K}_{h}.

Since $\delta_{h}$ is the adjoint operator of $\,{\rm d}^{-}:~{}V_{h}^{-}\rightarrow V_{h}$ , we have the following discrete Poincaré inequality for $\delta_{h}$ as well; cf. [7] and [6, Lemma 4].

Lemma 3.2 (discrete Poincaré inequality for $\delta_{h}$ ).

Let $C_{p}$ be the constant in (26). Then we have

\|\omega_{h}\|\leq C_{p}\|\delta_{h}\omega_{h}\|\qquad\forall~{}~{}\omega_{h}\in\mathfrak{Z}_{0,h}.

3.2. A projection-based quasi-interpolation operator

In this section, we introduce a projection-based quasi-interpolation operator briefly mentioned in [8, Proposition 5.44] which is a generalization of projection based operators introduced for $H^{1},H({\rm curl\,})$ and $H(\operatorname{div})$ spaces in [12, 13, 14, 25]. We redefine this operator based on the Hodge decomposition and prove more properties of this operator: it is commuted with $\delta_{h}$ , stable in both $\|\,{\rm d}(\cdot)\|$ and $\|\delta_{h}(\cdot)\|$ norms, a $L^{2}$ orthogonal projection to the space $\mathfrak{Z}_{0,h}$ , an orthogonal projection operator in the inner product $\langle\,{\rm d}(\cdot),\,{\rm d}(\cdot)\rangle$ on the subspace $\mathfrak{K}_{h}$ , and has the same approximation properties as the classical interpolation operators.

For any given $v\in H_{0}\Lambda(\Omega)$ , define $P_{h}v\in\mathfrak{K}_{h}$ such that

(27)

\langle\,{\rm d}P_{h}v,\,{\rm d}\phi_{h}\rangle=\langle\,{\rm d}v,\,{\rm d}\phi_{h}\rangle\qquad\forall~{}~{}\phi_{h}\in\mathfrak{K}_{h}.

Equation (27) determines $P_{h}v\in\mathfrak{K}_{h}$ uniquely since the Poincaré inequality (26) implies $\langle\,{\rm d}(\cdot),\,{\rm d}(\cdot)\rangle$ is an inner product on the subspace $\mathfrak{K}_{h}$ .

For any $v\in H_{0}\Lambda(\Omega)$ , the Hodge decomposition (13) implies that there exist $v_{1}\in\mathfrak{K}^{-}$ and $v_{2}\in\mathfrak{K}$ such that

v=\,{\rm d}^{-}v_{1}\oplus^{\bot_{L^{2}}}v_{2}

The projection-based quasi-interpolation operator $I_{h}:~{}H_{0}\Lambda(\Omega)\rightarrow V_{h}$ is defined as:

(28)

I_{h}v=\,{\rm d}^{-}P_{h}^{-}v_{1}\oplus^{\bot_{L^{2}}}P_{h}v_{2}.

We have the following properties.

Lemma 3.3.

For any $v\in H_{0}\Lambda(\Omega)$ , there hold

\langle I_{h}v,\,{\rm d}^{-}\phi_{h}\rangle=\langle v,\,{\rm d}^{-}\phi_{h}\rangle\qquad\forall~{}~{}\phi_{h}\in V_{h}^{-},

and

\langle\,{\rm d}I_{h}v,\,{\rm d}\psi_{h}\rangle=\langle\,{\rm d}v,\,{\rm d}\psi_{h}\rangle\qquad\forall~{}~{}\psi_{h}\in V_{h}.

Here we denote $V_{h}^{-1}=\{0\}$ .

Proof.

For any $\phi_{h}\in V_{h}^{-}$ , the discrete Hodge decomposition (25) implies that there exists $\phi_{h,1}\in\mathfrak{K}_{h}^{-}$ such that $\,{\rm d}^{-}\phi_{h}=\,{\rm d}^{-}\phi_{h,1}$ , therefore

\langle I_{h}v,\,{\rm d}^{-}\phi_{h}\rangle=\langle\,{\rm d}^{-}P_{h}^{-}v_{1},\,{\rm d}^{-}\phi_{h,1}\rangle=\langle\,{\rm d}^{-}v_{1},\,{\rm d}^{-}\phi_{h,1}\rangle=\langle v,\,{\rm d}^{-}\phi_{1,h}\rangle=\langle v,\,{\rm d}^{-}\phi_{h}\rangle.

For any $\psi_{h}\in V_{h}$ , the discrete Hodge decomposition (25) implies that there exists $\psi_{h,1}\in\mathfrak{K}_{h}$ such that $\,{\rm d}\psi_{h}=\,{\rm d}\psi_{h,1}$ ,

\langle\,{\rm d}I_{h}v,\,{\rm d}\psi_{h}\rangle=\langle\,{\rm d}P_{h}v_{2},\,{\rm d}\psi_{h,1}\rangle=\langle\,{\rm d}v_{2},\,{\rm d}\psi_{h,1}\rangle=\langle\,{\rm d}v,\,{\rm d}\psi_{h}\rangle.

Then, the desired results are obtained. ∎

We have the following stability results of $I_{h}$ .

Lemma 3.4.

We have the following stability results of $I_{h}$ :

(1)

For any $v\in H_{0}\Lambda(\Omega)$ , there holds

$\|\,{\rm d}I_{h}v\|\leq\|\,{\rm d}v\|.$
(2)

For any $v\in H_{0}\Lambda(\Omega)\cap H^{*}\Lambda(\Omega)$ , it holds

$\delta_{h}I_{h}v=Q_{h}^{-}\delta v,$

where $Q_{h}:~{}L^{2}\Lambda(\Omega)\rightarrow V_{h}$ is the $L^{2}$ projection operator. Therefore

$\|\delta_{h}I_{h}v\|\leq\|\delta v\|.$

Proof.

(1) For any $v\in H_{0}\Lambda(\Omega)$ , there exist $v_{1}\in\mathfrak{K}^{-}$ and $v_{2}\in\mathfrak{K}$ such that

v=\,{\rm d}^{-}v_{1}\oplus^{\bot_{L^{2}}}v_{2}.

We have

\displaystyle\|\,{\rm d}I_{h}v\|=\|\,{\rm d}P_{h}v_{2}\|\leq\|\,{\rm d}v_{2}\|=\|\,{\rm d}v\|.

(2) For any $v\in H_{0}\Lambda(\Omega)\cap H^{*}\Lambda(\Omega)$ . Using the facts that $P_{h}v_{2}\in\delta_{h}^{+}V_{h}^{+}$ , $\delta_{h}\delta_{h}^{+}=0$ and $\delta_{h}I_{h}v\in\mathfrak{K}_{h}^{-}$ , for any $\phi_{h}\in\mathfrak{K}_{h}^{-}$ , we have

	$\displaystyle\langle\delta_{h}I_{h}v,\phi_{h}\rangle$	$\displaystyle=\langle\delta_{h}\,{\rm d}^{-}P_{h}^{-}v_{1},\phi_{h}\rangle=\langle\,{\rm d}^{-}P_{h}^{-}v_{1},\,{\rm d}^{-}\phi_{h}\rangle$
		$\displaystyle=\langle\,{\rm d}^{-}v_{1},\,{\rm d}^{-}\phi_{h}\rangle=\langle v,\,{\rm d}^{-}\phi_{h}\rangle$
		$\displaystyle=\langle\delta v,\phi_{h}\rangle=\langle Q_{h}^{-}\delta v,\phi_{h}\rangle.$

Using the orthogonality result $\mathfrak{Z}_{0,h}^{-}\bot\mathfrak{K}_{h}^{-}$ , we get the desired result. ∎

To get approximation properties of the projection-based quasi-interpolation operator $I_{h}$ , we need the de Rham complexes for smooth differential forms established in [10] and the following Sobolev embedding result

(29)

H_{0}\Lambda\cap H^{*}\Lambda\hookrightarrow H^{1}\Lambda,

which holds when $\Omega$ is convex Lipschitz domain.

Lemma 3.5.

Assume that $\Omega$ is smooth enough such that (29) holds, then for any $v\in H_{0}\Lambda(\Omega)\cap H^{r+1}\Lambda(\Omega)$ with $r\geq 1$ , we have

(30)		$\displaystyle\\|v-I_{h}v\\|$	$\displaystyle\lesssim h^{l}\\|v\\|_{l}\qquad\text{for}\quad 1\leq l\leq r,$
(31)		$\displaystyle\\|\,{\rm d}(v-I_{h}v)\\|$	$\displaystyle\lesssim h^{l}\\|\,{\rm d}v\\|_{l}\qquad\text{for}\quad 1\leq l\leq r.$

Proof.

We shall use the de Rham sequence in [10]. Then for any $v\in H_{0}\Lambda(\Omega)\cap H^{r}\Lambda(\Omega)$ , there exist $v_{1}\in\mathfrak{K}^{-}\cap H^{r}\Lambda^{-}(\Omega)$ and $v_{2}\in\mathfrak{K}\cap H^{r}\Lambda(\Omega)$ such that

v=\,{\rm d}^{-}v_{1}\oplus^{\bot_{L^{2}}}v_{2},

and

(32)

\|v_{2}\|_{l}\lesssim\|v\|_{l}.

Therefore,

I_{h}v=\,{\rm d}^{-}P_{h}^{-}v_{1}\oplus^{\bot_{L^{2}}}P_{h}v_{2}.

Note that $v_{2}$ is the solution of the problem

(33)

v_{2}=\delta^{+}s,\quad\,{\rm d}v_{2}=q\quad\text{in}\quad\Omega,\qquad\operatorname*{tr}s=0\quad\text{on}\quad\partial\Omega.

with $q=\,{\rm d}v$ . Then $P_{h}v_{2}$ is the mixed finite element approximation of $v_{2}$ in $V_{h}$ , the standard error estimates of the mixed finite element method [5, 18] implies that

\|v_{2}-P_{h}v_{2}\|\lesssim h^{l}\|v_{2}\|_{l}\lesssim h^{l}\|v\|_{l},

where in the second inequality, we have used (32).

Also note that $v_{1}$ is the solution of the problem

(34)

\delta\,{\rm d}^{-}v_{1}=g,\quad\delta^{-}v_{1}=0,\quad\operatorname*{tr}v_{1}=0\quad\text{on}\quad\partial\Omega,

with $g=\delta v$ . The definition of $P_{h}^{-}$ implies that $P_{h}^{-}v_{1}$ is the mixed finite element approximation of $v_{1}$ in $V_{h}^{-}$ , then the standard error estimates for the mixed finite element methods [5, 18] imply

\|\,{\rm d}^{-}(v_{1}-P_{h}^{-}v_{1})\|\lesssim h^{l}\|\,{\rm d}^{-}v_{1}\|_{l}\lesssim h^{l}\|g\|_{l-1}=h^{l}\|\delta v\|_{l-1}\lesssim\|v\|_{l}.

Therefore,

\|v-I_{h}v\|\leq\|\,{\rm d}^{-}(v_{1}-P_{h}^{-}v_{1})\|+\|v_{2}-P_{h}v_{2}\|\lesssim h^{l}\|v\|_{l}.

We turn to the estimates of (31). Since for any $\phi_{h}\in\mathfrak{K}_{h}$ , it holds

\displaystyle\langle\,{\rm d}P_{h}v_{2},\,{\rm d}\phi_{h}\rangle

\displaystyle=\langle\,{\rm d}v_{2},\,{\rm d}\phi_{h}\rangle=\langle\,{\rm d}\pi_{h}v_{2}+\,{\rm d}(I-\pi_{h})v_{2},\,{\rm d}\phi_{h}\rangle,

where $\pi_{h}:~{}H_{0}\Lambda(\Omega)\cap H^{r}\Lambda(\Omega)\rightarrow V_{h}$ is the classical interpolation operator [1]. Therefore,

\,{\rm d}P_{h}v_{2}=\,{\rm d}\pi_{h}v_{2}+Q_{\mathfrak{Z}_{0,h}^{+}}\,{\rm d}(I-\pi_{h})v_{2},

where $Q_{\mathfrak{Z}_{0,h}^{+}}:~{}L^{2}\Lambda^{+}\rightarrow\mathfrak{Z}_{0,h}^{+}$ is the $L^{2}$ orthogonal projection operator. Then, we have

\displaystyle\|\,{\rm d}(v-I_{h}v)\|=\|\,{\rm d}(v_{2}-P_{h}v_{2})\|\leq\|(I-\pi_{h}^{+})\,{\rm d}v_{2}\|\lesssim h^{l}\|\,{\rm d}v\|_{l}.

∎

3.3. Semi-discretization and error analysis

The semi-discrete formulation [26, 3] of (18)-(20) is: Given $f\in L^{2}((0,T),L^{2}\Lambda)$ , find $\boldsymbol{u}_{h}=(\sigma_{h},\mu_{h},\omega_{h})^{\intercal}:~{}(0,T]\mapsto V_{h}^{-}\times V_{h}\times V_{h}^{+}:=\boldsymbol{W}_{h}$ such that

(35)	$\displaystyle\langle\sigma_{h,t},\tau_{h}\rangle-\langle\,{\rm d}^{-}\tau_{h},\mu_{h}\rangle$	$\displaystyle=0\,\qquad\qquad\forall~{}~{}\tau_{h}\in V_{h}^{-},$
(36)	$\displaystyle\langle\mu_{h,t},v_{h}\rangle+\langle\,{\rm d}^{-}\sigma_{h},v_{h}\rangle+\langle\omega_{h},\,{\rm d}v_{h}\rangle$	$\displaystyle=\langle f,v_{h}\rangle\qquad\forall~{}~{}v_{h}\in V_{h},$
(37)	$\displaystyle\langle\omega_{h,t},\phi_{h}\rangle-\langle\,{\rm d}\mu_{h},\phi_{h}\rangle$	$\displaystyle=0\,\qquad\qquad\forall~{}~{}\phi_{h}\in V_{h}^{+},$

with initial values

\sigma_{h}(\cdot,0)=I_{h}^{-}\delta u_{0},\quad\mu_{h}(\cdot,0)=I_{h}u_{1},\quad\omega_{h}(\cdot,0)=I_{h}^{+}\,{\rm d}u_{0}.

Introduce

\mathcal{A}_{h}=\begin{pmatrix}0&\delta_{h}&0\\ -\,{\rm d}^{-}&0&-\delta_{h}^{+}\\ 0&\,{\rm d}&0\end{pmatrix}

(35) - (37) can be rewritten as

(38)

\langle\boldsymbol{u}_{h,t},\boldsymbol{v}_{h}\rangle+\langle\mathcal{A}_{h}\boldsymbol{u}_{h},\boldsymbol{v}_{h}\rangle=\langle\boldsymbol{F},\boldsymbol{v}_{h}\rangle\qquad\forall~{}~{}\boldsymbol{v}_{h}\in\boldsymbol{W}_{h}.

Following the same line as the proof of Theorem 2.2, we have the energy estimates.

Theorem 3.6.

Let $\boldsymbol{u}_{h}=(\sigma_{h},\mu_{h},\omega_{h})^{\intercal}\in\boldsymbol{W}_{h}$ be the solution of the mixed formulation (35)-(37) or (38). Provided $f\in L^{1}((0,T),L^{2}\Lambda)$ , we have the energy bound

(39)

\sup\limits_{0\leq s\leq T}\|\boldsymbol{u}_{h}(\cdot,t)\|\leq\|\boldsymbol{u}_{h}(\cdot,0)\|+2\int_{0}^{T}\|f(\cdot,s)\|\,{\rm d}s.

Furthermore, if $f\in W^{1,1}((0,T),L^{2}\Lambda)$ , we have the bound

(40)

\sup\limits_{0\leq t\leq T}\|\mathcal{A}_{h}\boldsymbol{u}_{h}(\cdot,t)\|\leq\|\mathcal{A}_{h}\boldsymbol{u}_{h}(\cdot,0)\|+4\|f\|_{L^{\infty}(L^{2})}+2\int_{0}^{T}\|f_{t}\|\,{\rm d}t.

When $f=0$ , these inequalities become equalities and we have the energy conservation

\|\boldsymbol{u}_{h}(\cdot,t)\|=\|\boldsymbol{u}_{h}(\cdot,0)\|,and\ \|\mathcal{A}_{h}\boldsymbol{u}_{h}(\cdot,t)\|=\|\mathcal{A}_{h}\boldsymbol{u}_{h}(\cdot,0)\|\quad\forall t>0.

Remark 3.7.

Since (35)-(37) (or (38)) is a linear system, Theorem 3.6 implies the existence and uniqueness of the solution at any time level $t\in(0,T]$ . Also, when the source term $f=0$ , Theorem 3.6 implies that the energies $\|\boldsymbol{u}_{h}(\cdot,t)\|$ and $\|\mathcal{A}_{h}\boldsymbol{u}_{h}(\cdot,t)\|$ are preserved exactly.

The rest of this section will focus on the error estimates of the semi-discretization (35)-(37) (or its simplified form (38)). We denote

\mathcal{I}_{h}=\begin{pmatrix}I_{h}^{-}&0&0\\ 0&I_{h}&0\\ 0&0&I_{h}^{+}\end{pmatrix}.

Then for any $\boldsymbol{v}_{h}=(\tau_{h},v_{h},\phi_{h})^{\intercal}\in\boldsymbol{W}_{h}$ , (23) is equivalent to

(41)

\langle\mathcal{I}_{h}\boldsymbol{u}_{t},\boldsymbol{v}_{h}\rangle+\langle\mathcal{A}_{h}\mathcal{I}_{h}\boldsymbol{u},\boldsymbol{v}_{h}\rangle=\langle\boldsymbol{F},\boldsymbol{v}_{h}\rangle+\langle\Theta_{h,t},\boldsymbol{v}_{h}\rangle+\langle\mathcal{A}_{h}\mathcal{I}_{h}\boldsymbol{u}-\mathcal{A}\boldsymbol{u},\boldsymbol{v}_{h}\rangle

with

\Theta_{h}=\mathcal{I}_{h}\boldsymbol{u}-\boldsymbol{u}.

Using the properties of the projection-based quasi-interpolation operator $I_{h}$ , we obtain

	$\displaystyle\langle\mathcal{A}_{h}\mathcal{I}_{h}\boldsymbol{u}-\mathcal{A}\boldsymbol{u},\boldsymbol{v}_{h}\rangle$	$\displaystyle=\langle I_{h}\mu-\mu,\,{\rm d}^{-}\tau_{h}\rangle-\langle\,{\rm d}^{-}(I_{h}^{-}-I)\sigma,v_{h}\rangle$
		$\displaystyle\quad-\langle I_{h}^{+}\omega-\omega,\,{\rm d}v_{h}\rangle+\langle\,{\rm d}(I_{h}-I)\mu,\phi_{h}\rangle$
		$\displaystyle=-\langle\,{\rm d}^{-}(I_{h}^{-}-I)\sigma,v_{h}\rangle+\langle\,{\rm d}(I_{h}-I)\mu,\phi_{h}\rangle$
		$\displaystyle=\langle\boldsymbol{G},\boldsymbol{V}_{h}\rangle,$

with $\boldsymbol{G}=(0,-\,{\rm d}^{-}(\pi_{h}^{-}-I)\sigma,\,{\rm d}(I_{h}-I)\mu)^{\intercal}$ . Denote

\mathcal{E}_{h}=\mathcal{I}_{h}\boldsymbol{u}-\boldsymbol{u}_{h}

and subtracting the semi-discrete form (38) from (41), we get

(42)

\displaystyle\langle\mathcal{E}_{h,t},\boldsymbol{v}_{h}\rangle+\langle\mathcal{A}_{h}\mathcal{E}_{h},\boldsymbol{v}_{h}\rangle=\langle\Theta_{h,t}+\boldsymbol{G},\boldsymbol{v}_{h}\rangle\qquad\forall~{}~{}\boldsymbol{v}_{h}\in\boldsymbol{W}_{h}.

We have the following estimate of $\mathcal{E}_{h}$ .

Lemma 3.8.

Suppose the exact solution $\boldsymbol{u}=(\sigma,\mu,\omega)$ of (23) has time derivatives $\sigma_{t}\in L^{1}((0,T),H^{r}\Lambda^{-})$ , $\mu_{t}\in L^{1}((0,T),H^{r}\Lambda)$ and $\omega_{t}\in L^{1}((0,T),H^{r}\Lambda^{+})$ with $r\geq 1$ . Then, for any $1\leq m\leq r$ and $t\in[0,T]$ , we have the bound

\|\mathcal{E}_{h}(\cdot,t)\|\lesssim h^{m}\int_{0}^{T}\left(\|\boldsymbol{u}_{t}\|_{m}+\|\,{\rm d}^{-}\sigma\|_{m}+\|\,{\rm d}\mu\|_{m}\right)\,{\rm d}t.

Proof.

The fact that $\mathcal{E}_{h}(\cdot,0)=0$ and Theorem 3.6 implies

\sup\limits_{0\leq t\leq T}\|\mathcal{E}_{h}(\cdot,t)\|\leq 2\int_{0}^{T}\|\Theta_{h,t}+\boldsymbol{G}\|\,{\rm d}t.

Using the triangle inequality and the approximation properties of $I_{h}$ , the desired result follows. ∎

We then obtain the following estimates by the triangle inequality, Lemma 3.5, and 3.6.

Theorem 3.9.

Suppose the exact solution $\boldsymbol{u}=(\sigma,\mu,\omega)^{\intercal}$ of (23) has time derivatives $\sigma_{t}\in L^{1}((0,T),H^{r}\Lambda^{-})$ , $\mu_{t}\in L^{1}((0,T),H^{r}\Lambda)$ and $\omega_{t}\in L^{1}((0,T),H^{r}\Lambda^{+})$ with $r\geq 1$ . Let $\boldsymbol{u}_{h}$ be the exact solution of (38). Then, for any $1\leq m\leq r$ and $t\in[0,T]$ , we have the bound

\displaystyle\|\boldsymbol{u}(\cdot,t)-\boldsymbol{u}_{h}(\cdot,t)\|

\displaystyle\lesssim h^{m}\left(\|\boldsymbol{u}\|_{L^{\infty}(H^{m})}+\int_{0}^{T}(\|\boldsymbol{u}_{t}\|_{m}+\|\,{\rm d}^{-}\sigma\|_{m}+\|\,{\rm d}\mu\|_{m})\,{\rm d}t\right).

Remark 3.10.

In this theorem, the convergence order $m$ is determined by the polynomial order of the finite element spaces preserved.

Remark 3.11.

It should be point out that in [26], Quenneville has obtained an error estimates for the semi-discretization in the form

	$\displaystyle\\|\boldsymbol{u}-\boldsymbol{u}_{h}\\|_{L^{\infty}(L^{2})}$	$\displaystyle\leq\\|\pi_{h}\boldsymbol{u}_{0}-\boldsymbol{u}_{0,h}\\|+\\|\pi_{h}\boldsymbol{u}-\boldsymbol{u}\\|_{L^{\infty}(L^{2})}$
		$\displaystyle\quad+(1+T)(\\|(\pi_{h}\boldsymbol{u}-\boldsymbol{u})(\cdot,0)\\|+\\|\pi_{h}\boldsymbol{u}_{t}-\boldsymbol{u}_{t}\\|_{L^{1}(L^{2})}),$

where $\pi_{h}$ is an elliptic projection operator. Comparing with this result, ours do not have the factor $1+T$ and thus is more robust to the time variable. ∎

We now give to the error estimates in the energy norm $\|\mathcal{A}\boldsymbol{u}-\mathcal{A}_{h}\boldsymbol{U}_{h}\|$ , which is equivalent to $\|d^{-}\sigma-d^{-}\sigma_{h}\|+\|\,{\rm d}\mu-\,{\rm d}\mu_{h}\|+\|\delta\mu-\delta_{h}\mu\|+\|\delta^{+}\omega-\delta_{h}^{+}\omega_{h}\|$ . Note that it is possible that the $L^{2}$ -norm is small but the energy norm is larger due to the small oscillation in the error. We shall show the energy norm is still of the same order of convergence. We give the estimate of $\|\mathcal{A}_{h}\mathcal{E}_{h}\|$ first. By Lemma 3.5 and 3.6, we have the following estimate.

Lemma 3.12.

Suppose the exact solution $\boldsymbol{u}=(\sigma,\mu,\omega)^{\intercal}$ of (23) has time derivatives $\sigma_{tt}\in L^{1}((0,T),H^{r}\Lambda^{-})$ , $\mu_{tt}\in L^{1}((0,T),H^{r}\Lambda)$ and $\omega_{tt}\in L^{1}((0,T),H^{r}\Lambda^{+})$ with $r\geq 1$ . Then, for any $1\leq m\leq r$ and $t\in[0,T]$ , we have the bound

	$\displaystyle\\|\mathcal{A}_{h}\mathcal{E}_{h}(\cdot,t)\\|$	$\displaystyle\lesssim h^{m}(\\|\boldsymbol{u}_{t}\\|_{L^{\infty}(H^{m})}+\\|\,{\rm d}^{-}\sigma\\|_{L^{\infty}(H^{m})}+\\|\,{\rm d}\mu\\|_{L^{\infty}(H^{m})})$
		$\displaystyle\quad+h^{m}\int_{0}^{T}\left(\\|\boldsymbol{u}_{tt}\\|_{m}+\\|\,{\rm d}^{-}\sigma_{t}\\|_{m}+\\|\,{\rm d}\mu_{t}\\|_{m}\right)\,{\rm d}t$

Proof.

Using the fact that $\mathcal{E}_{h}(\cdot,0)=0$ and Theorem 3.6, we have

\|\mathcal{A}_{h}\mathcal{E}_{h}(\cdot,t)\|\leq 4\|\Theta_{h,t}+\boldsymbol{G}\|_{L^{\infty}(L^{2})}+2\int_{0}^{T}\|\Theta_{h,tt}+\boldsymbol{G}_{t}\|\,{\rm d}t.

Triangle inequality and Lemma 3.5 imply the desired result. ∎

Theorem 3.13.

	$\displaystyle\\|\mathcal{A}_{h}\boldsymbol{u}_{h}(\cdot,t)-\mathcal{A}\boldsymbol{u}(\cdot,t)\\|\lesssim$	$\displaystyle h^{m}{\large[}\,\\|\boldsymbol{u}_{t}\\|_{L^{\infty}(H^{m})}+\\|\,{\rm d}^{-}\sigma\\|_{L^{\infty}(H^{m})}+\\|\,{\rm d}\mu\\|_{L^{\infty}(H^{m})}$
		$\displaystyle\\|\mathcal{A}\boldsymbol{u}\\|_{L^{\infty}(H^{m})}+\int_{0}^{T}(\\|\boldsymbol{u}_{tt}\\|_{m}+\\|\,{\rm d}^{-}\sigma_{t}\\|_{m}+\\|\,{\rm d}\mu_{t}\\|)\,{\rm d}t\,{\large]}.$

Proof.

The triangle inequality and Lemma 3.4 imply that

	$\displaystyle\\|\mathcal{A}_{h}\boldsymbol{u}_{h}(\cdot,t)-\mathcal{A}\boldsymbol{u}(\cdot,t)\\|$	$\displaystyle\lesssim\\|\,{\rm d}^{-}(\sigma-I_{h}^{-}\sigma)\\|+\\|\,{\rm d}(\mu-I_{h}\mu)\\|+\\|\delta\mu-\delta_{h}I_{h}\mu\\|$
		$\displaystyle\quad+\\|\delta^{+}\omega-\delta^{+}_{h}I_{h}^{+}\omega\\|+\\|\mathcal{A}_{h}\mathcal{E}_{h}(\cdot,t)\\|$
		$\displaystyle=\\|\,{\rm d}^{-}(I-I_{h}^{-})\sigma\\|+\\|\,{\rm d}(\mu-I_{h}\mu)\\|+\\|(I-Q_{h}^{-})\delta\mu\\|$
		$\displaystyle\quad+\\|(I-Q_{h})\delta^{+}\omega\\|+\\|\mathcal{A}_{h}\mathcal{E}_{h}(\cdot,t)\\|.$

Using the properties of the $L^{2}$ projection operators, Lemma 3.5 and 3.12, we get the desired results. ∎

4. Full-discretization

In this section, we will consider the full discretization. We will use a second order continuous time Galerkin method [16] to discretize time variable and will obtain the energy estimates and optimal error estimates.

Energy conservation numerical schemes can have a crucial influence on the quality of the numerical simulations. In long-time simulations, energy-preserving can have a dramatic effect on stability and global error growth. The numerical schemes are not automatically inherit from the semi-discretization and a lot of time discretization methods cannot preserve the energies exactly. These led us to pay more attentions on the time discretization.

4.1. Time discretization

Let $\mathcal{T}_{\,\Delta t}$ denote the equispaced partition of the interval $(0,T)$ with $\,\Delta t=T/N$ and $N$ the number of elements in $\mathcal{T}_{\,\Delta t}$ . For $1\leq i\leq N$ , we denote $t_{i}=i\,\Delta t$ and $\tau_{i}=(t_{i-1},t_{i})$ with $t_{0}=0$ . For any quantity $v(t)$ , we denote $v^{i}=v(t_{i})$ . Define $\mathcal{P}_{1}(\mathcal{T}_{\,\Delta t})$ (abbr. $\mathcal{P}_{1}$ ) as the set of continuous piecewise linear polynomials with respect to the time variable $t$ on $\mathcal{T}_{\,\Delta t}$ and $\mathcal{P}_{0}(\mathcal{T}_{\,\Delta t})$ (abbr. $\mathcal{P}_{0}$ ) as the set of piecewise constant with respect to the time variable $t$ on $\mathcal{T}_{\,\Delta t}$ . For any Sobolev space $S$ associates with the spatial variables, we use $\mathcal{P}_{1}(S)$ to denote the set of functions that are continuous piecewise linear polynomials with respect to the time variable $t$ and in the Sobolev space $S$ with respect to the spatial variables. $\mathcal{P}_{0}(S)$ is defined similarly.

The full discrete formulation of the Hodge wave equation (23) can be written as: Find $\boldsymbol{U}_{h}=(\tilde{\sigma}_{h},\tilde{\mu}_{h},\tilde{\omega}_{h})^{\intercal}\in\mathcal{P}_{1}(\boldsymbol{W}_{h})$ such that

(43)

\int_{0}^{T}\left(\langle\boldsymbol{U}_{h,t},\boldsymbol{V}_{h}\rangle+\langle\mathcal{A}_{h}\boldsymbol{U}_{h},\boldsymbol{V}_{h}\rangle\right)\,{\rm d}t=\int_{0}^{T}\langle\boldsymbol{F},\boldsymbol{V}_{h}\rangle\,{\rm d}t\qquad\forall~{}~{}~{}\boldsymbol{V}_{h}\in\mathcal{P}_{0}(\boldsymbol{W}_{h}).

Remark 4.1.

The full discrete formulation (43) is equivalent to

\int_{t_{i-1}}^{t_{i}}(\langle\boldsymbol{U}_{h,t},\boldsymbol{V}_{h}\rangle+\langle\mathcal{A}_{h}\boldsymbol{U}_{h},\boldsymbol{V}_{h}\rangle)\,{\rm d}t=\int_{t_{i-1}}^{t_{i}}\langle\boldsymbol{F},\boldsymbol{V}_{h}\rangle\,{\rm d}t\quad\forall~{}\boldsymbol{V}_{h}\in\mathcal{P}_{0}(\boldsymbol{W}_{h}),\ 1\leq i\leq N.

The fact that

\displaystyle\int_{t_{i-1}}^{t_{i}}\langle\boldsymbol{U}_{h,t},\boldsymbol{V}_{h}\rangle\,{\rm d}t

\displaystyle=\langle\boldsymbol{U}_{h}^{i}-\boldsymbol{U}_{h}^{i-1},\boldsymbol{V}_{h}\rangle,

and

\int_{t_{i-1}}^{t_{i}}\langle\mathcal{A}_{h}\boldsymbol{U}_{h},\boldsymbol{V}_{h}\rangle\,{\rm d}t=\frac{\,\Delta t}{2}\langle\mathcal{A}_{h}(\boldsymbol{U}_{h}^{i}+\boldsymbol{U}_{h}^{i-1}),\boldsymbol{V}_{h}\rangle,

implies the full discrete formulation (43) is essentially a Crank-Nicolson scheme with exact time integration of the right hand side.

We have the following energy estimates for (43).

Theorem 4.2.

Let $\boldsymbol{U}_{h}=(\tilde{\sigma}_{h},\tilde{\mu}_{h},\tilde{\omega}_{h})^{\intercal}\in\mathcal{P}_{1}(\boldsymbol{W}_{h})$ be the solutions of (43). Assume that $f\in L^{\infty}((0,T),L^{2}\Lambda)$ , then there hold the following energy bound

(44)

\max\limits_{0\leq i\leq N}\|\boldsymbol{U}_{h}^{i}\|\leq\|\boldsymbol{U}_{h}^{0}\|+2\int_{0}^{T}\|\boldsymbol{F}\|\,{\rm d}t.

When $F=0$ , the inequality becomes equality and we have the energy conservation

\|\boldsymbol{U}_{h}^{i}\|=\|\boldsymbol{U}_{h}^{0}\|,\quad\forall\ 1\leq i\leq N.

Proof.

Taking $\boldsymbol{V}_{h}$ in (43) as

\boldsymbol{V}_{h}|_{\tau_{i}}=\boldsymbol{U}_{h}^{i}+\boldsymbol{U}_{h}^{i-1}\qquad\text{and}\qquad\boldsymbol{V}_{h}|_{\mathcal{T}_{\,\Delta t}\setminus\tau_{i}}=0,

we obtain

\int_{t_{i-1}}^{t_{i}}\langle\boldsymbol{U}_{h,t},\boldsymbol{U}_{h}^{i}+\boldsymbol{U}_{h}^{i-1}\rangle\,{\rm d}t+\int_{t_{i-1}}^{t_{i}}\langle\mathcal{A}_{h}\boldsymbol{U}_{h},\boldsymbol{U}_{h}^{i}+\boldsymbol{U}_{h}^{i-1}\rangle\,{\rm d}t=\int_{t_{i-1}}^{t_{i}}\langle\boldsymbol{F},\boldsymbol{U}_{h}^{i}+\boldsymbol{U}_{h}^{i-1}\rangle\,{\rm d}t.

The fact that

\displaystyle\int_{t_{i-1}}^{t_{i}}\langle\mathcal{A}_{h}\boldsymbol{U}_{h},\boldsymbol{U}_{h}^{i}+\boldsymbol{U}_{h}^{i-1}\rangle\,{\rm d}t

\displaystyle=\frac{\,\Delta t}{2}\langle\mathcal{A}_{h}(\boldsymbol{U}_{h}^{i}+\boldsymbol{U}_{h}^{i-1}),\boldsymbol{U}_{h}^{i}+\boldsymbol{U}_{h}^{i-1}\rangle=0

and

\displaystyle\int_{t_{i-1}}^{t_{i}}\langle\boldsymbol{U}_{h,t},\boldsymbol{U}_{h}^{i}+\boldsymbol{U}_{h}^{i-1}\rangle\,{\rm d}t

\displaystyle=\|\boldsymbol{U}_{h}^{i}\|^{2}-\|\boldsymbol{U}_{h}^{i-1}\|^{2}

imply

\displaystyle\|\boldsymbol{U}_{h}^{i}\|^{2}-\|\boldsymbol{U}_{h}^{i-1}\|^{2}

\displaystyle=\int_{t_{i-1}}^{t_{i}}\langle\boldsymbol{F},\boldsymbol{U}_{h}^{i}+\boldsymbol{U}_{h}^{i-1}\rangle\,{\rm d}t\leq 2\max\limits_{0\leq i\leq N}\|\boldsymbol{U}_{h}^{i}\|\int_{t_{i-1}}^{t_{i}}\|\boldsymbol{F}\|\,{\rm d}t.

Summing over $i$ from $1$ to $m\leq N$ , we get

\displaystyle\|\boldsymbol{U}_{h}^{m}\|^{2}-\|\boldsymbol{U}_{h}^{0}\|^{2}\leq 2\max\limits_{0\leq i\leq N}\|\boldsymbol{U}_{h}^{i}\|\int_{0}^{T}\|\boldsymbol{F}\|\,{\rm d}t.

Therefore, the desired result follows by Lemma 2.1. ∎

Remark 4.3.

Since (43) is a linear system, therefore Theorem 4.2 implies the existence and uniqueness of the solution for the full discrete form (43).

Theorem 4.4.

Let $\boldsymbol{U}_{h}=(\tilde{\sigma}_{h},\tilde{\mu}_{h},\tilde{\omega}_{h})^{\intercal}\in\mathcal{P}_{1}(\boldsymbol{W}_{h})$ be the solution of (43). Assume that $f\in W^{1,1}((0,T),L^{2}\Lambda)$ , then there holds the following energy bound

(45)

\max\limits_{0\leq i\leq N}\|\mathcal{A}_{h}\boldsymbol{U}_{h}^{i}\|\leq\|\mathcal{A}_{h}\boldsymbol{U}_{h}^{0}\|+4\|\boldsymbol{F}\|_{L^{\infty}(L^{2})}+2\int_{0}^{T}\|\boldsymbol{F}_{t}\|\,{\rm d}t.

When $F=0$ , the inequality becomes equality and we have the energy conservation

\|\mathcal{A}_{h}\boldsymbol{U}_{h}^{i}\|=\|\mathcal{A}_{h}\boldsymbol{U}_{h}^{0}\|\quad\forall\ 1\leq i\leq N.

Proof.

Taking $\boldsymbol{V}_{h}$ in (43) as

\boldsymbol{V}_{h}|_{\tau_{i}}=\mathcal{A}_{h}(\boldsymbol{U}_{h}^{i}-\boldsymbol{U}_{h}^{i-1})\quad\text{and}\quad\boldsymbol{V}_{h}|_{\mathcal{T}_{\,\Delta t}\setminus\tau_{i}}=0,

we have

	$\displaystyle\int_{t_{i-1}}^{t_{i}}\langle\boldsymbol{U}_{h,t},\mathcal{A}_{h}(\boldsymbol{U}_{h}^{i}-\boldsymbol{U}_{h}^{i-1})\rangle\,{\rm d}t$	$\displaystyle+\int_{t_{i-1}}^{t_{i}}\langle\mathcal{A}_{h}\boldsymbol{U}_{h},\mathcal{A}_{h}(\boldsymbol{U}_{h}^{i}-\boldsymbol{U}_{h}^{i-1})\rangle\,{\rm d}t$
		$\displaystyle=\int_{t_{i-1}}^{t_{i}}\langle\boldsymbol{F},\mathcal{A}_{h}(\boldsymbol{U}_{h}^{i}-\boldsymbol{U}_{h}^{i-1})\rangle\,{\rm d}t.$

Using the fact that

\int_{t_{i-1}}^{t_{i}}\langle\boldsymbol{U}_{h,t},\mathcal{A}_{h}(\boldsymbol{U}_{h}^{i}-\boldsymbol{U}_{h}^{i-1})\rangle\,{\rm d}t=\langle\boldsymbol{U}_{h}^{i}-\boldsymbol{U}_{h}^{i-1},\mathcal{A}_{h}(\boldsymbol{U}_{h}^{i}-\boldsymbol{U}_{h}^{i-1})\rangle=0

and

\int_{t_{i-1}}^{t_{i}}\langle\mathcal{A}_{h}\boldsymbol{U}_{h},\mathcal{A}_{h}(\boldsymbol{U}_{h}^{i}-\boldsymbol{U}_{h}^{i-1})\rangle=\frac{\,\Delta t}{2}(\|\mathcal{A}_{h}\boldsymbol{U}_{h}^{i}\|^{2}-\|\mathcal{A}_{h}\boldsymbol{U}_{h}^{i-1}\|^{2}),

we get

\displaystyle\|\mathcal{A}_{h}\boldsymbol{U}_{h}^{i}\|^{2}-\|\mathcal{A}_{h}\boldsymbol{U}_{h}^{i-1}\|^{2}

\displaystyle=2\int_{t_{i-1}}^{t_{i}}\left\langle\boldsymbol{F},\mathcal{A}_{h}\frac{\boldsymbol{U}_{h}^{i}-\boldsymbol{U}_{h}^{i-1}}{\,\Delta t}\right\rangle\,{\rm d}t=2\int_{t_{i-1}}^{t_{i}}\langle\boldsymbol{F},\mathcal{A}_{h}\boldsymbol{U}_{h,t}\rangle\,{\rm d}t.

Summing over $i$ from $1$ to $m\leq N$ , we obtain

	$\displaystyle\\|\mathcal{A}_{h}\boldsymbol{U}_{h}^{m}\\|^{2}$	$\displaystyle=\\|\mathcal{A}_{h}\boldsymbol{U}_{h}^{0}\\|^{2}+2\int_{0}^{t_{m}}\langle\boldsymbol{F},\mathcal{A}_{h}\boldsymbol{U}_{h,t}\rangle\,{\rm d}t$
		$\displaystyle=\\|\mathcal{A}_{h}\boldsymbol{U}_{h}^{0}\\|^{2}+2\langle\boldsymbol{F}^{m},\mathcal{A}_{h}\boldsymbol{U}_{h}^{m}\rangle-2\langle\boldsymbol{F}^{0},\mathcal{A}_{h}\boldsymbol{U}_{h}^{0}\rangle-2\int_{0}^{t_{m}}\langle\boldsymbol{F}_{t},\mathcal{A}_{h}\boldsymbol{U}_{h}\rangle\,{\rm d}t$
		$\displaystyle\leq\\|\mathcal{A}_{h}\boldsymbol{U}_{h}^{0}\\|^{2}+2\max\limits_{0\leq i\leq N}\\|\mathcal{A}_{h}\boldsymbol{U}_{h}^{i}\\|\left(2\\|\boldsymbol{F}\\|_{L^{\infty}(L^{2})}+\int_{0}^{T}\\|\boldsymbol{F}_{t}\\|\,{\rm d}t\right).$

Then the desired result follows by a direct using of Lemma 2.1. ∎

4.2. Error analysis of the full discretization

In this subsection, we turn to the error estimates of the full discrete formulation (43). We bound the error of the full discrete formulation in various norms. Let

\displaystyle\boldsymbol{e}_{h}=\mathcal{I}_{h}\boldsymbol{u}-\boldsymbol{U}_{h}.

Simple caculation shows that for any $\boldsymbol{V}_{h}\in\mathcal{P}_{0}(\boldsymbol{W}_{h})$ , $\boldsymbol{e}_{h}$ satisfies the following equation

(46)

\int_{0}^{T}(\langle\boldsymbol{e}_{h,t},\boldsymbol{V}_{h}\rangle+\langle\mathcal{A}_{h}\boldsymbol{e}_{h},\boldsymbol{V}_{h}\rangle)\,{\rm d}t=\int_{0}^{T}\langle\Theta_{h}+\boldsymbol{G},\boldsymbol{V}_{h}\rangle\,{\rm d}t,\quad\forall~{}~{}\boldsymbol{V}_{h}\in\mathcal{P}_{0}(\boldsymbol{W}_{h}).

Then, we have the following estimates of $\boldsymbol{e}_{h}$ .

Lemma 4.5.

Let $\boldsymbol{u}=(\sigma,\mu,\omega)^{\intercal}$ be the solution of (23) and $\boldsymbol{U}_{h}=(\tilde{\sigma}_{h},\tilde{\mu}_{h},\tilde{\omega}_{h})^{\intercal}\in\mathcal{P}_{1}(\boldsymbol{W}_{h})$ be the solution of (43). Assume that $\sigma_{tt}\in L^{1}((0,T),H^{r}\Lambda^{-})$ , $\mu_{tt}\in L^{1}((0,T),H^{r}\Lambda)$ and $\omega_{tt}\in L^{1}((0,T),H^{r}\Lambda^{+})$ with $r\geq 1$ . Then, for any $1\leq m\leq r$ , we have the bound

\displaystyle\max\limits_{0\leq i\leq N}\|\boldsymbol{e}_{h}^{i}\|\lesssim

\displaystyle h^{m}\int_{0}^{T}(\|\boldsymbol{u}_{h}\|_{m}+\|\,{\rm d}^{-}\sigma\|_{m}+\|\,{\rm d}\mu\|_{m})\,{\rm d}t+\,\Delta t^{2}\int_{0}^{T}\|\mathcal{A}\boldsymbol{u}_{tt}\|\,{\rm d}t.

Proof.

Taking $\boldsymbol{V}_{h}$ in (46) as

\boldsymbol{V}_{h}|_{\tau_{i}}=\boldsymbol{e}_{h}^{i-1}+\boldsymbol{e}_{h}^{i}\quad\text{and}\quad\boldsymbol{V}_{h}|_{\mathcal{T}_{\,\Delta t}\setminus\tau_{i}}=0,

we obtain

\displaystyle\int_{t_{i-1}}^{t_{i}}\langle\boldsymbol{e}_{h,t},\boldsymbol{e}_{h}^{i-1}+\boldsymbol{e}_{h}^{i}\rangle\,{\rm d}t

\displaystyle+\int_{t_{i-1}}^{t_{i}}\langle\mathcal{A}_{h}\boldsymbol{e}_{h},\boldsymbol{e}_{h}^{i}+\boldsymbol{e}_{h}^{i+1}\rangle\,{\rm d}t=\int_{t_{i-1}}^{t_{i}}\langle\Theta_{h}+\boldsymbol{G},\boldsymbol{e}_{h}^{i-1}+\boldsymbol{e}_{h}^{i}\rangle\,{\rm d}t.

Note that

\displaystyle\int_{t_{i-1}}^{t_{i}}\langle\boldsymbol{e}_{h,t},\boldsymbol{e}_{h}^{i-1}+\boldsymbol{e}_{h}^{i}\rangle\,{\rm d}t

\displaystyle=\|\boldsymbol{e}_{h}^{i}\|^{2}-\|\boldsymbol{e}_{h}^{i-1}\|^{2},

and

		$\displaystyle\int_{t_{i-1}}^{t_{i}}\langle\mathcal{A}_{h}\boldsymbol{e}_{h},\boldsymbol{e}_{h}^{i-1}+\boldsymbol{e}_{h}^{i}\rangle\,{\rm d}t=\int_{t_{i-1}}^{t_{i}}\langle\mathcal{A}_{h}(I-J_{\,\Delta t})\boldsymbol{e}_{h},\boldsymbol{e}_{h}^{i-1}+\boldsymbol{e}_{h}^{i}\rangle\,{\rm d}t$
	$\displaystyle\leq$	$\displaystyle 2\max\limits_{0\leq i\leq N}\\|\boldsymbol{e}_{h}^{i}\\|\int_{t_{i-1}}^{t_{i}}\\|\mathcal{A}_{h}(I-J_{\,\Delta t})\boldsymbol{e}_{h}\\|\,{\rm d}t\leq C\,\Delta t^{2}\max\limits_{0\leq ileqN}\\|\boldsymbol{e}_{h}^{i}\\|\int_{t_{i-1}}^{t_{i}}\\|\mathcal{A}_{h}\mathcal{I}_{h}\boldsymbol{u}_{tt}\\|\,{\rm d}t$
	$\displaystyle\leq$	$\displaystyle C\,\Delta t^{2}\max\limits_{0\leq ileqN}\\|\boldsymbol{e}_{h}^{i}\\|\int_{t_{i-1}}^{t_{i}}\\|\mathcal{A}\boldsymbol{u}_{tt}\\|\,{\rm d}t,$

we then have

		$\displaystyle\\|\boldsymbol{e}_{h}^{i}\\|^{2}-\\|\boldsymbol{e}_{h}^{i-1}\\|^{2}\leq C\max\limits_{0\leq i\leq N}\\|\boldsymbol{e}_{h}^{i}\\|\left(\int_{t_{i-1}}^{t_{i}}(\\|\Theta_{h}\\|+\\|\boldsymbol{G}\\|)\,{\rm d}t+\,\Delta t^{2}\int_{t_{i-1}}^{t_{i}}\\|\mathcal{A}\boldsymbol{u}_{tt}\\|\,{\rm d}t\right)$
	$\displaystyle\lesssim$	$\displaystyle\max\limits_{0\leq i\leq N}\\|\boldsymbol{e}_{h}^{i}\\|\left(h^{m}\int_{t_{i-1}}^{t_{i}}(\\|\boldsymbol{u}_{h}\\|_{m}+\\|\,{\rm d}^{-}\sigma\\|_{m}+\\|\,{\rm d}\mu\\|_{m})\,{\rm d}t+\,\Delta t^{2}\int_{t_{i-1}}^{t_{i}}\\|\mathcal{A}\boldsymbol{u}_{tt}\\|\,{\rm d}t\right).$

Summing over $i$ from $1$ to $m\leq N$ , we get

\displaystyle\|\boldsymbol{e}_{h}^{m}\|^{2}

\displaystyle\lesssim\max\limits_{0\leq i\leq N}\|\boldsymbol{e}_{h}^{i}\|\left(h^{m}\int_{0}^{T}(\|\boldsymbol{u}\|_{m}+\|\,{\rm d}^{-}\sigma\|_{m}+\|\,{\rm d}\mu\|_{m})\,{\rm d}t+\,\Delta t^{2}\int_{0}^{T}\|\mathcal{A}\boldsymbol{u}_{tt}\|\,{\rm d}t\right).

Then the desired result follows. ∎

Using triangle inequality, Lemma 4.5 and the properties of $I_{h}$ , we have the following $L^{2}$ error estimate.

Theorem 4.6.

	$\displaystyle\max\limits_{0\leq i\leq N}\\|\boldsymbol{u}^{i}-\boldsymbol{U}_{h}^{i}\\|\lesssim$	$\displaystyle h^{m}\left(\\|\boldsymbol{u}\\|_{L^{\infty}(H^{m})}+\\|\boldsymbol{u}\\|_{L^{1}(H^{m})}+\\|\,{\rm d}^{-}\sigma\\|_{L^{1}(H^{m})}+\\|\,{\rm d}\mu\\|_{L^{1}(H^{m})}\right)$
		$\displaystyle+\,\Delta t^{2}\\|\mathcal{A}\boldsymbol{u}_{tt}\\|_{L^{1}(L^{2})}$

Now, we turn to the estimates of the energy norm $\|\mathcal{A}(\cdot)\|$ error. We first estimate $\|\mathcal{A}_{h}\boldsymbol{e}_{h}\|$ .

Lemma 4.7.

Let $\boldsymbol{u}=(\sigma,\mu,\omega)^{\intercal}$ be the solution of (23) and $\boldsymbol{U}_{h}=(\tilde{\sigma}_{h},\tilde{\mu}_{h},\tilde{\omega}_{h})^{\intercal}\in\mathcal{P}_{1}(\boldsymbol{W}_{h})$ be the solution of (43). Assume that $\sigma_{tt}\in L^{\infty}((0,T),H^{r}\Lambda^{-})$ , $\mu_{tt}\in L^{\infty}((0,T),H^{r}\Lambda)$ and $\omega_{tt}\in L^{\infty}((0,T),H^{r}\Lambda^{+})$ with $r\geq 1$ . Then, for any $1\leq m\leq r$ we have the bound

\displaystyle\|\mathcal{A}_{h}\boldsymbol{e}_{h}\|

\displaystyle\lesssim h^{m}\left(\|\boldsymbol{u}_{t}\|_{L^{1}(H^{m})}+\|\,{\rm d}^{-}\sigma_{t}\|_{L^{1}(H^{m})}+\|\,{\rm d}\mu_{t}\|_{L^{1}(H^{m})}\right)+\,\Delta t^{2}\|\mathcal{A}\boldsymbol{u}_{tt}\|_{L^{\infty}(L^{2})}.

Proof.

Taking $\boldsymbol{V}_{h}$ in (46) as

\boldsymbol{V}_{h}|_{\tau_{i}}=\mathcal{A}_{h}(\boldsymbol{e}_{h}^{i}-\boldsymbol{e}_{h}^{i-1})\quad\text{and}\quad\boldsymbol{V}_{h}|_{\mathcal{T}_{\,\Delta t}\setminus\tau_{i}}=0,

we obtain

	$\displaystyle\int_{t_{i-1}}^{t_{i}}\langle\boldsymbol{e}_{h,t},\mathcal{A}_{h}(\boldsymbol{e}_{h}^{i}-\boldsymbol{e}_{h}^{i-1})\rangle\,{\rm d}t$	$\displaystyle+\int_{t_{i-1}}^{t_{i}}\langle\mathcal{A}_{h}\boldsymbol{e}_{h},\mathcal{A}_{h}(\boldsymbol{e}_{h}^{i}-\boldsymbol{e}_{h}^{i-1})\rangle\,{\rm d}t$
		$\displaystyle=\int_{t_{i-1}}^{t_{i}}\langle\Theta_{h}+\boldsymbol{G},\mathcal{A}_{h}(\boldsymbol{e}_{h}^{i}-\boldsymbol{e}_{h}^{i-1})\rangle\,{\rm d}t.$

Note that

\displaystyle\int_{t_{i-1}}^{t_{i}}\langle\boldsymbol{e}_{h,t},\mathcal{A}_{h}(\boldsymbol{e}_{h}^{i}-\boldsymbol{e}_{h}^{i-1})\rangle\,{\rm d}t

\displaystyle=\langle(\boldsymbol{e}_{h}^{i}-\boldsymbol{e}_{h}^{i-1}),\mathcal{A}_{h}(\boldsymbol{e}_{h}^{i}-\boldsymbol{e}_{h}^{i-1})\rangle=0,

and

	$\displaystyle\int_{t_{i-1}}^{t_{i}}\langle\mathcal{A}_{h}\boldsymbol{e}_{h},\mathcal{A}_{h}(\boldsymbol{e}_{h}^{i}-\boldsymbol{e}_{h}^{i-1})\rangle\,{\rm d}t$	$\displaystyle=\frac{\,\Delta t}{2}\left(\\|\mathcal{A}_{h}\boldsymbol{e}_{h}^{i}\\|^{2}-\\|\mathcal{A}_{h}\boldsymbol{e}_{h}^{i-1}\\|^{2}\right)$
		$\displaystyle\quad+\int_{t_{i-1}}^{t_{i}}\langle\mathcal{A}_{h}(I-J_{\,\Delta t})\boldsymbol{e}_{h},\mathcal{A}_{h}(\boldsymbol{e}_{h}^{i}-\boldsymbol{e}_{h}^{i-1})\rangle\,{\rm d}t.$

Therefore,

		$\displaystyle\\|\mathcal{A}_{h}\boldsymbol{e}_{h}^{i}\\|^{2}-\\|\mathcal{A}_{h}\boldsymbol{e}_{h}^{i-1}\\|^{2}$
	$\displaystyle=$	$\displaystyle\frac{2}{\,\Delta t}\int_{t_{i-1}}^{t_{i}}\langle\Theta_{h}+\boldsymbol{G},\mathcal{A}_{h}(\boldsymbol{e}_{h}^{i}-\boldsymbol{e}_{h}^{i-1})\rangle\,{\rm d}t-\frac{2}{\,\Delta t}\int_{t_{i-1}}^{t_{i}}\langle\mathcal{A}_{h}(I-J_{\,\Delta t})\boldsymbol{e}_{h},\mathcal{A}_{h}(\boldsymbol{e}_{h}^{i}-\boldsymbol{e}_{h}^{i-1})\rangle\,{\rm d}t$
	$\displaystyle=$	$\displaystyle 2\int_{t_{i-1}}^{t_{i}}\left\langle\Theta_{h}+\boldsymbol{G},\mathcal{A}_{h}\frac{\partial}{\partial t}(J_{\,\Delta t}\boldsymbol{e}_{h})\right\rangle\,{\rm d}t-2\int_{t_{i-1}}^{t_{i}}\left\langle\mathcal{A}_{h}(I-J_{\,\Delta t})\boldsymbol{e}_{h},\mathcal{A}_{h}\frac{\partial}{\partial t}(J_{\,\Delta t}\boldsymbol{e}_{h})\right\rangle\,{\rm d}t.$

Summing over $i$ from $1$ to $m\leq N$ , we get

		$\displaystyle\\|\mathcal{A}_{h}\boldsymbol{e}_{h}^{m}\\|^{2}$
	$\displaystyle=$	$\displaystyle 2\int_{0}^{t_{m}}\left\langle\Theta_{h}+\boldsymbol{G},\mathcal{A}_{h}\frac{\partial}{\partial t}(J_{\,\Delta t}\boldsymbol{e}_{h})\right\rangle\,{\rm d}t-2\int_{0}^{t_{m}}\left\langle\mathcal{A}_{h}(I-J_{\,\Delta t})\boldsymbol{e}_{h},\mathcal{A}_{h}\frac{\partial}{\partial t}(J_{\,\Delta t}\boldsymbol{e}_{h})\right\rangle\,{\rm d}t$
	$\displaystyle=$	$\displaystyle 2\langle\Theta_{h}^{m}+\boldsymbol{G}^{m},\mathcal{A}_{h}\boldsymbol{e}_{h}^{m}\rangle-2\int_{0}^{t_{m}}\langle\Theta_{h,t}+\boldsymbol{G}_{t},\mathcal{A}_{h}J_{\,\Delta t}\boldsymbol{e}_{h}\rangle\,{\rm d}t$
		$\displaystyle+2\int_{0}^{t_{m}}\left\langle\mathcal{A}_{h}\frac{\partial}{\partial t}(I-J_{\,\Delta t})\boldsymbol{e}_{h},\mathcal{A}_{h}(J_{\,\Delta t}\boldsymbol{e}_{h})\right\rangle\,{\rm d}t$
	$\displaystyle\leq$	$\displaystyle 2\max\limits_{0\leq i\leq N}\\|\mathcal{A}_{h}\boldsymbol{e}_{h}^{i}\\|\left(\max\limits_{0\leq i\leq N}\\|\Theta_{h}^{i}+\boldsymbol{G}^{i}\\|+\int_{0}^{T}\\|\Theta_{h,t}+\boldsymbol{G}_{t}\\|\,{\rm d}t+\int_{0}^{T}\left\\|\mathcal{A}_{h}\frac{\partial}{\partial t}(I-J_{\,\Delta t})\boldsymbol{e}_{h}\right\\|\,{\rm d}t\right)$
	$\displaystyle\leq$	$\displaystyle C\max\limits_{0\leq i\leq N}\\|\mathcal{A}_{h}\boldsymbol{e}_{h}^{i}\\|\int_{0}^{T}\left(\\|\Theta_{h,t}\\|+\\|\boldsymbol{G}_{t}\\|+\left\\|\mathcal{A}_{h}\frac{\partial}{\partial t}(I-J_{\,\Delta t})\boldsymbol{e}_{h}\right\\|\right)\,{\rm d}t$
	$\displaystyle\leq$	$\displaystyle C\max\limits_{0\leq i\leq N}\\|\mathcal{A}_{h}\boldsymbol{e}_{h}^{i}\\|\left[h^{m}\left(\\|\boldsymbol{u}_{t}\\|_{L^{1}(H^{m})}+\\|\,{\rm d}^{-}\sigma_{t}\\|_{L^{1}(H^{m})}+\\|\,{\rm d}\mu_{t}\\|_{L^{1}(H^{m})}+\,\Delta t^{2}\\|\\|\mathcal{A}\boldsymbol{u}_{tt}\\|_{L^{\infty}(L^{2})}\right)\right],$

where in the last inequality, we have used the properties of the one dimensional interpolation operator. Then the desired result follows. ∎

We summarize the error estimate for the full discretization below.

Theorem 4.8.

Let $\boldsymbol{u}=(\sigma,\mu,\omega)^{\intercal}$ be the solution of (23) and $\boldsymbol{U}_{h}=(\tilde{\sigma}_{h},\tilde{\mu}_{h},\tilde{\omega}_{h})^{\intercal}\in\mathcal{P}_{1}(\boldsymbol{W}_{h})$ be the solution of (43). Assume that $\sigma_{ttt}\in L^{\infty}((0,T),H^{r}\Lambda^{-})$ , $\mu_{ttt}\in L^{\infty}((0,T),H^{r}\Lambda)$ and $\omega_{ttt}\in L^{\infty}((0,T),H^{r}\Lambda^{+})$ with $r\geq 1$ . Then, for any $1\leq m\leq r$ we have the bound

	$\displaystyle\max\limits_{0\leq i\leq N}\\|\mathcal{A}\boldsymbol{u}^{i}-\mathcal{A}_{h}\boldsymbol{U}_{h}^{i}\\|\lesssim$	$\displaystyle\,h^{m}\left[\\|\boldsymbol{u}_{t}\\|_{L^{1}(H^{m})}+\\|\mathcal{A}\boldsymbol{u}\\|_{L^{\infty}(H^{m})}\right.$
		$\displaystyle+\left.\\|\,{\rm d}^{-}\sigma_{t}\\|_{L^{1}(H^{m})}+\\|\,{\rm d}\mu_{t}\\|_{L^{1}(H^{m})}\right]$
	$\displaystyle+$	$\displaystyle\,\Delta t^{2}\\|\mathcal{A}\boldsymbol{u}_{tt}\\|_{L^{\infty}(L^{2})}.$

5. Numerical experiments

In this section, we will give some simple numerical examples to illustrate the theoretical results. We consider the Hodge wave equation on the unit square $\Omega=(0,1)^{2}$ , i.e., $n=2$ and compute the cases $k=0$ , $k=1$ and $k=2$ .

5.1. The case $k=0$

The Hodge wave equation presents in the standard $H^{1}(\Omega)$ or $H({\rm curl\,})$ language reads as (note that $\sigma=\delta\mu=0$ ): Find $\mu\in H_{0}^{1}$ and $\boldsymbol{\omega}\in H_{0}(\operatorname*{rot})$ or $\mu\in H_{0}({\rm curl\,})$ and $\boldsymbol{\omega}\in H_{0}(\operatorname{div})$ such that

(47)		$\displaystyle(\mu_{t},v)+(\boldsymbol{\omega},{\rm grad\,}v)$	$\displaystyle=(f,v)\qquad\forall~{}~{}\boldsymbol{v}\in H_{0}^{1}(\Omega),$
(48)		$\displaystyle(\boldsymbol{\omega}_{t},\boldsymbol{\phi})-({\rm grad\,}\mu,\boldsymbol{\phi})$	$\displaystyle=0\qquad\forall~{}~{}\boldsymbol{\phi}\in H_{0}(\operatorname*{rot}).$

(49)		$\displaystyle(\mu_{t},v)+(\boldsymbol{\omega},{\rm curl\,}v)$	$\displaystyle=(f,v)\qquad\forall~{}~{}\boldsymbol{v}\in H_{0}({\rm curl\,}),$
(50)		$\displaystyle(\boldsymbol{\omega}_{t},\boldsymbol{\phi})-({\rm curl\,}\mu,\boldsymbol{\phi})$	$\displaystyle=0\qquad\forall~{}~{}\boldsymbol{\phi}\in H_{0}(\operatorname{div}).$

We only give the numerical results for (49)-(50). We choose the exact solutions as

	$\displaystyle\mu(x,y,t)$	$\displaystyle=e^{-t}\sin(\pi x)\sin(\pi y),$
	$\displaystyle\boldsymbol{\omega}(x,y,t)$	$\displaystyle=-\pi e^{-t}\begin{pmatrix}\sin(\pi x)\cos(\pi y)\\ -\cos(\pi x)\sin(\pi y)\end{pmatrix}.$

The initial conditions are $\mu_{0}=\mu(x,y,0)$ and $\boldsymbol{\omega}_{0}=\boldsymbol{\omega}(x,y,0)$ . We use piecewise continuous second order polynomial to discrete $\mu$ and use $RT_{1}$ element [27] to discrete $\omega$ , the numerical results are listed in Table 1.

Table 1. Errors and convergence orders in various norms with

\,\Delta t=0.0001

t=0.0004

$h$	$\\|\mu-\mu_{h}\\|$	$\\|{\rm curl\,}(\mu-\mu_{h})\\|$	$\\|\omega-\omega_{h}\\|$
1/4	3.8253e-03	1.3417e-01	1.3164e-01
1/8	5.0314e-04	3.5025e-02	3.3567e-02
1/16	6.7850e-05	9.5850e-03	8.4467e-03
order	2.914	1.904	1.981

From Table 1, we can see that the mixed finite element method is of second order convergence rate for the variables ${\rm curl\,}\mu$ and $\omega$ , and is of third order convergence rate for $\mu$ . All these variables have optimal convergence order.

5.2. The case $k=1$

The Hodge wave equation is: Find $\sigma\in H_{0}({\rm curl\,})$ , $\boldsymbol{\mu}\in H_{0}(\operatorname{div})$ and $\omega\in L_{0}^{2}(\Omega)$ such that

	$\displaystyle(\sigma_{t},\tau)-(\boldsymbol{\mu},{\rm curl\,}\tau)$	$\displaystyle=0\qquad\forall~{}~{}\tau\in H_{0}({\rm curl\,}),$
	$\displaystyle(\boldsymbol{\mu}_{t},\boldsymbol{v})+({\rm curl\,}\sigma,\boldsymbol{v})+(\omega,\operatorname{div}\boldsymbol{v})$	$\displaystyle=(\boldsymbol{f},\boldsymbol{v})\qquad\forall~{}~{}\boldsymbol{v}\in H_{0}(\operatorname{div}),$
	$\displaystyle(\omega_{t},\phi)-(\operatorname{div}\boldsymbol{\mu},\phi)$	$\displaystyle=0\qquad\forall~{}~{}\phi\in L_{0}^{2}(\Omega).$

This formulation can be viewed as the mixed method for the time-harmonic Maxwell’s equations with divergence free constrain on both $\boldsymbol{\mu}$ and $\boldsymbol{f}$ . The formulation can also be viewed as the mixed method for the elastic wave equation. We use continuous piecewise quadratic polynomial to discrete $\sigma$ , use $RT_{1}$ element [27] to discrete $\mu$ and use discontinuous piecewise linear polynomial to discrete $\omega$ . Firstly, we choose the exact solutions as

(51)	$\displaystyle\sigma(x,y,t)$	$\displaystyle=2e^{-t}\left(\pi\sin^{2}(\pi x)\sin(\pi y)\cos(\pi y)-x(x-1)(2x-1)y^{2}(y-1)^{2}\right)$
(52)	$\displaystyle\boldsymbol{\mu}(x,y,t)$	$\displaystyle=e^{-t}\begin{pmatrix}\sin^{2}(\pi x)\sin^{2}(\pi y)\\ x^{2}(x-1)^{2}y^{2}(y-1)^{2}\end{pmatrix},$
(53)	$\displaystyle\begin{split}\omega(x,y,t)&=-2e^{-t}\left(\pi\sin(\pi x)\cos(\pi x)\sin^{2}(\pi y)\right.\\ &\quad\left.+x^{2}(x-1)^{2}y(y-1)(2y-1)\right),\end{split}$

with initial conditions $\sigma_{0}=\sigma(x,y,0)$ , $\boldsymbol{\mu}_{0}=\boldsymbol{\mu}(x,y,0)$ and $\omega_{0}=\omega(x,y,0)$ . The numerical results are listed in Table 2. We also test the long time robustness of our algorithm, the numerical results are listed in Table 3.

Table 2. Errors and convergence orders in various norms with

\,\Delta t=0.0001

t=0.0004

$h$	$\\|\sigma-\sigma_{h}\\|$	$\\|{\rm curl\,}(\sigma-\sigma_{h})\\|$	$\\|\mu-\mu_{h}\\|$	$\\|\operatorname{div}(\mu-\mu_{h})\\|$	$\\|\omega-\omega_{h}\\|$
1/4	4.8563e-02	1.6691e+00	2.9085e-02	0.1391	0.1388
1/8	6.4968e-03	4.4475e-01	7.6216e-03	0.0364	0.0367
1/16	8.2584e-04	1.1298e-01	1.9354e-03	0.0093	0.0093
order	2.949	1.942	1.950	1.950	1.9763

Table 3. Long time problem with

\,\Delta t=0.1

and

h=1/16

$T$	$\\|\sigma-\sigma_{h}\\|$	$\\|{\rm curl\,}(\sigma-\sigma_{h})\\|$	$\\|\mu-\mu_{h}\\|$	$\\|\operatorname{div}(\mu-\mu_{h})\\|$	$\\|\omega-\omega_{h}\\|$
10	5.4138e-04	1.5087e-02	3.7500e-01	1.3603	0.2374
30	4.8684e-04	1.8186e-02	3.7502e-01	1.3604	0.2486
50	6.1267e-04	1.4339e-02	3.7502e-01	1.3604	0.4158

Then, we chose $\boldsymbol{f}=0$ and the initial conditions $\sigma_{0}=\sigma(x,y,0)$ , $\boldsymbol{\mu}_{0}=\boldsymbol{\mu}(x,y,0)$ and $\omega_{0}=\omega(x,y,0)$ with $\sigma$ , $\boldsymbol{\mu}$ and $\omega$ defined as in (51) - (53). We compute the energies $\|\boldsymbol{U}_{h}^{i}\|$ and $\|\mathcal{A}_{h}\boldsymbol{U}_{h}^{i}\|$ on different time levels, the numerical results are showing in Fig. 1.

Refer to caption — Figure 1. Energies $\|\boldsymbol{U}_{h}\|$ and $\|\mathcal{A}_{h}\boldsymbol{U}_{h}\|$ in different times with $h=1/16$ and $\,\Delta t=0.25$ .

From this example, we have the following observations.

(1)

The mix finite element method is of second order convergence rate for the variables ${\rm curl\,}\sigma$ , $\mu$ , $\operatorname{div}\mu$ and $\omega$ , and is of third order convergence rate for the variable $\sigma$ . All of these variables have optimal convergence order.
(2)

From Table 3, we can see that the mixed finite element method is robust for long-time problem.
(3)

From Fig. 1, we can see that the mixed finite element method conserves the energies $\|\boldsymbol{U}_{h}\|$ and $\|\mathcal{A}_{h}\boldsymbol{U}_{h}\|$ exactly.

5.3. The case $k=2$

The Hodge wave equation presents in the $H(\operatorname{div})$ and $L^{2}$ language reads as (note that $\omega=\,{\rm d}\mu=0$ ): Find $\boldsymbol{\sigma}\in H_{0}(\operatorname{div})$ and $\mu\in L_{0}^{2}(\Omega)$ such that

(54)		$\displaystyle(\boldsymbol{\sigma}_{t},\boldsymbol{\tau})+(\mu,\operatorname{div}\boldsymbol{\tau})$	$\displaystyle=0\qquad\forall~{}~{}\boldsymbol{\tau}\in H_{0}(\operatorname{div}),$
(55)		$\displaystyle(\mu_{t},v)-(\operatorname{div}\boldsymbol{\sigma},v)$	$\displaystyle=(f,v)\qquad\forall~{}~{}v\in L_{0}^{2}(\Omega).$

This formulation is the mixed method for acoustic wave equations [22]. We choose the exact solutions as

(56)		$\displaystyle\boldsymbol{\sigma}(x,y,t)$	$\displaystyle=-\pi e^{-t}\begin{pmatrix}\cos(\pi x)\sin(\pi y)\\ \sin(\pi x)\cos(\pi y)\end{pmatrix},$
(57)		$\displaystyle\mu(x,y,t)$	$\displaystyle=e^{-t}\sin(\pi x)\sin(\pi y),$

and pick initial conditions $\boldsymbol{\sigma}_{0}=\boldsymbol{\sigma}(x,y,0)$ and $\mu_{0}=\mu(x,y,0)$ . We use $RT_{1}$ element to discrete $\boldsymbol{\sigma}$ and use discontinuous piecewise linear polynomials to discrete $\mu$ , the numerical results are listed in Table 4.

Table 4. Errors and convergence orders in various norms with

\,\Delta t=0.0001

t=0.0004

$h$	$\\|\sigma-\sigma_{h}\\|$	$\\|\operatorname{div}(\sigma-\sigma_{h})\\|$	$\\|\mu-\mu_{h}\\|$
1/4	5.6058e-02	3.8201e-01	1.9350e-02
1/8	1.4002e-02	9.6901e-02	4.9051e-03
1/16	3.4958e-03	2.4360e-02	1.2312e-03
order	2.002	1.985	1.992

From Table 4, we can see that the mixed finite element method is of second order convergence rate for all the variables.

Acknowledgments

We would like to thank Professor Long Chen from University of California at Irvine for valuable discussion and suggestions.

References

[1] D. N. Arnold, R. S. Falk, and R. Winther. Finite element exterior calculus, homological techniques, and applications. Acta Numerica, 15:1–155, 2006.
[2] D. N. Arnold, R. S. Falk, and R. Winther. Finite element exterior calculus: from hodge theory to numerical stability. Bulletin of the American Mathematical Society, 47(2):281–354, 2010.
[3] D. N. Arnold. Finite element exterior calculus. SIAM, 2018.
[4] G. A. Baker. Error estimates for finite element methods for second order hyperbolic equations. SIAM Journal on Numerical Analysis, 13(4):564–576, 1976.
[5] F. Brezzi and M. Fortin. Mixed and hybrid finite element methods. Springer-Verlag, New York, 1991.
[6] L. Chen and Y. Wu. Convergence of adaptive mixed finite element methods for the hodge laplacian equation: without harmonic forms. SIAM Journal on Numerical Analysis, 55(6):2905–2929, 2017.
[7] L. Chen, Y. Wu, L. Zhong, and J. Zhou. Multigrid preconditioners for mixed finite element methods of the vector laplacian. Journal of Scientific Computing, 77:101–128, 2018.
[8] S. H. Christiansen, H. Z. Munthe-Kaas and B. Owren. Topics in structure-preserving discretization. Acta Numerica, pp 1-119, 2011.
[9] E. T. Chung and B. Engquist. Optimal discontinuous galerkin methods for the acoustic wave equation in higher dimensions. SIAM Journal on Numerical Analysis, 47(5):3820–3848, 2009.
[10] M. Costabel and A. McIntosh. On Bogovskil̆ and regularized Poincaré integral operators for de Rahm complexes on Lipschitz domains. Math. Z., 265: 297–320, 2010.
[11] L. C. Cowsat, T. F. Dupont, and M. F. Wheeler. A priori estimates for mixed finite element methods for the wave equation. Computer Methods in Applied Mechanics and Engineering, 82(1-3):205–222, 1990.
[12] L. Demkowicz, P. Monk, L. Vardapetyan, and W. Rachowicz. De Rham diagram for $hp$ finite element spaces. Math. Comput. Applications, 39(7-8):29-38, 2000.
[13] L. Demkowicz, I. Babuška. Optimal $p$ interpolation error estimates for edge finite elements of variable order in 2D. SIAM J. Numer. Anal., 41(4):1195–1208, 2003.
[14] L. Demkowicz and A. Buffa. $H^{1}$ , $H({\rm curl\,})$ and $H(\operatorname{div})$ -conforming projection-based interpolation in three dimensions Quasi-optimal $p$ -interpolation estimates. Comput. Methods Appl. Mech. Engrg, 194: 267 – 296, 2005.
[15] T. Dupont. L2-estimates for Galerkin methods for second order hyperbolic equations. SIAM Journal on Numerical Analysis, 10(5):880–889, 1973.
[16] D. French and T. Peterson. A continuous space-time finite element method for the wave equation. Mathematics of Computation, 65(214):491–506, 1996.
[17] T. Geveci. On the application of mixed finite element methods to the wave equations. ESAIM: Mathematical Modelling and Numerical Analysis, 22(2):243–250, 1988.
[18] V. Girault and P.-A. Raviart. Finite element methods for Navier-Stokes equations: theory and algorithms, volume 5. Springer Science & Business Media, 2012.
[19] R. Glowinski and T. Rossi. A mixed formulation and exact controllability approach for the computation of the periodic solutions of the scalar wave equation. (i): Controllability problem formulation and related iterative solution. Academie Des Sciences Comptes Rendus Mathematique, 343(7):493–498, 2006.
[20] R. Griesmaier and P. Monk. Discretization of the wave equation using continuous elements in time and a hybridizable discontinuous Galerkin method in space. J Sci Comput, 58: 472-498.
[21] E. W. Jenkins. Numerical solution of the acoustic wave equation using raviart–thomas elements. Journal of Computational and Applied Mathematics, 206(1):420–431, 2007.
[22] R. C. Kirby and T. T. Kieu. Symplectic-mixed finite element approximation of linear acoustic wave equations. Numerische Mathematik, 130(2):257–291, Oct 2014.
[23] D. Mitrea, M. Mitrea, and M. Taylor. Layer potentials, the Hodge Laplacian, and global boundary problems in nonsmooth Riemannian manifolds, volume 713. American Mathematical Soc., 2001.
[24] M. Mitrea. Dirichlet integrals and gaffney-friedrichs inequalities in convex domains. Forum Math, 13:531–567, 2001.
[25] J. T. Oden, L. Demkowicz, W. Rachowicz, and T. A. Westermann. Towards a universal $h$ - $p$ adaptive finite element strategy, Part2. A posteriori error estimation. Comput. Methods Appl. Meth. Eng., 77: 113–180, 1989.
[26] V. Quenneville-Bélair. A new approach to finite element simulations of general relativity. Ph. D. thesis, University of Minnesota, 2015.
[27] P. A. Raviart and J. M. Thomas. A mixed finite element method for 2-nd order elliptic problems. Mathematical aspects of finite element methods, pages 292–315, 1977.
[28] B. Riviére and M. F. Wheeler. A priori error estimates for mixed finite element approximations of the acoustic wave equation. SIAM Journal on Numerical Analysis, 40(5):1698–1715, 2002.

(30)		$\displaystyle\\|v-I_{h}v\\|$	$\displaystyle\lesssim h^{l}\\|v\\|_{l}\qquad\text{for}\quad 1\leq l\leq r,$
(31)		$\displaystyle\\|\,{\rm d}(v-I_{h}v)\\|$	$\displaystyle\lesssim h^{l}\\|\,{\rm d}v\\|_{l}\qquad\text{for}\quad 1\leq l\leq r.$

	$\displaystyle\\|\mathcal{A}_{h}\boldsymbol{u}_{h}(\cdot,t)-\mathcal{A}\boldsymbol{u}(\cdot,t)\\|$	$\displaystyle\lesssim\\|\,{\rm d}^{-}(\sigma-I_{h}^{-}\sigma)\\|+\\|\,{\rm d}(\mu-I_{h}\mu)\\|+\\|\delta\mu-\delta_{h}I_{h}\mu\\|$
		$\displaystyle\quad+\\|\delta^{+}\omega-\delta^{+}_{h}I_{h}^{+}\omega\\|+\\|\mathcal{A}_{h}\mathcal{E}_{h}(\cdot,t)\\|$
		$\displaystyle=\\|\,{\rm d}^{-}(I-I_{h}^{-})\sigma\\|+\\|\,{\rm d}(\mu-I_{h}\mu)\\|+\\|(I-Q_{h}^{-})\delta\mu\\|$
		$\displaystyle\quad+\\|(I-Q_{h})\delta^{+}\omega\\|+\\|\mathcal{A}_{h}\mathcal{E}_{h}(\cdot,t)\\|.$

	$\displaystyle\\|\mathcal{A}_{h}\boldsymbol{U}_{h}^{m}\\|^{2}$	$\displaystyle=\\|\mathcal{A}_{h}\boldsymbol{U}_{h}^{0}\\|^{2}+2\int_{0}^{t_{m}}\langle\boldsymbol{F},\mathcal{A}_{h}\boldsymbol{U}_{h,t}\rangle\,{\rm d}t$
		$\displaystyle=\\|\mathcal{A}_{h}\boldsymbol{U}_{h}^{0}\\|^{2}+2\langle\boldsymbol{F}^{m},\mathcal{A}_{h}\boldsymbol{U}_{h}^{m}\rangle-2\langle\boldsymbol{F}^{0},\mathcal{A}_{h}\boldsymbol{U}_{h}^{0}\rangle-2\int_{0}^{t_{m}}\langle\boldsymbol{F}_{t},\mathcal{A}_{h}\boldsymbol{U}_{h}\rangle\,{\rm d}t$
		$\displaystyle\leq\\|\mathcal{A}_{h}\boldsymbol{U}_{h}^{0}\\|^{2}+2\max\limits_{0\leq i\leq N}\\|\mathcal{A}_{h}\boldsymbol{U}_{h}^{i}\\|\left(2\\|\boldsymbol{F}\\|_{L^{\infty}(L^{2})}+\int_{0}^{T}\\|\boldsymbol{F}_{t}\\|\,{\rm d}t\right).$

	$\displaystyle\max\limits_{0\leq i\leq N}\\|\mathcal{A}\boldsymbol{u}^{i}-\mathcal{A}_{h}\boldsymbol{U}_{h}^{i}\\|\lesssim$	$\displaystyle\,h^{m}\left[\\|\boldsymbol{u}_{t}\\|_{L^{1}(H^{m})}+\\|\mathcal{A}\boldsymbol{u}\\|_{L^{\infty}(H^{m})}\right.$
		$\displaystyle+\left.\\|\,{\rm d}^{-}\sigma_{t}\\|_{L^{1}(H^{m})}+\\|\,{\rm d}\mu_{t}\\|_{L^{1}(H^{m})}\right]$
	$\displaystyle+$	$\displaystyle\,\Delta t^{2}\\|\mathcal{A}\boldsymbol{u}_{tt}\\|_{L^{\infty}(L^{2})}.$

Energy-preserving mixed finite element methods for the Hodge wave equation

Abstract.

Key words and phrases:

2010 Mathematics Subject Classification:

1. Introduction

2. Preliminaries

2.1. de Rham complex

2.2. The Hodge wave equation

Lemma 2.1.

Theorem 2.2.

Proof.

Remark 2.3.

Remark 2.4.

3. Semi-discretization of the Hodge wave equation

3.1. Finite element spaces

Lemma 3.1 (discrete Poincaré inequality for d\,{\rm d}).

Lemma 3.2 (discrete Poincaré inequality for δh\delta_{h}).

3.2. A projection-based quasi-interpolation operator

Lemma 3.3.

Proof.

Lemma 3.4.

Proof.

Lemma 3.5.

Proof.

3.3. Semi-discretization and error analysis

Theorem 3.6.

Remark 3.7.

Lemma 3.8.

Proof.

Theorem 3.9.

Remark 3.10.

Remark 3.11.

Lemma 3.12.

Proof.

Theorem 3.13.

Proof.

4. Full-discretization

4.1. Time discretization

Remark 4.1.

Theorem 4.2.

Proof.

Remark 4.3.

Theorem 4.4.

Proof.

4.2. Error analysis of the full discretization

Lemma 4.5.

Proof.

Theorem 4.6.

Lemma 4.7.

Proof.

Theorem 4.8.

5. Numerical experiments

5.1. The case k=0k=0

5.2. The case k=1k=1

5.3. The case k=2k=2

Acknowledgments

References

Lemma 3.1 (discrete Poincaré inequality for $\,{\rm d}$ ).

Lemma 3.2 (discrete Poincaré inequality for $\delta_{h}$ ).

5.1. The case $k=0$

5.2. The case $k=1$

5.3. The case $k=2$