A $C^0P_2$ time-stepping virtual element method for linear wave equations on polygonal meshes

1.   Introduction

In the past few years, a new numerical method called the virtual element method (VEM) has received much attention in the field of scientific computing. The pioneering works can be found in [3,2,5,11]. VEMs can be regarded as a generalization of usual finite element methods that can easily handle polytopal meshes, which have several significant advantages over standard finite element methods, such as (1) they are convenient for solving problems on complex geometric domains; (2) they are suitable for attacking problems associated with high regularity admissible spaces (see [8]). Until now, such methods have been used to solve a variety of mathematical physical problems. We refer to [4,7,9,13,14,17,18,19,20,23,27,28,29,30,31] and references therein for details about their applications for steady variational problems in primal formulations.

However, there are few works on unsteady problems using VEMs. Some relevant results can be found in [27,28]. In [28], after temporal discretization by the usual backward Euler method and spatial discretization by the VEM in [2], a fully discrete method was proposed and analyzed for linear parabolic equations. Similar ideas were further extended to handle linear wave equations in [27] with the focus on discussing the semidiscrete method. A fully discrete method is developed by the Newmark method and the Bathe method for temporal discretization, and the error analysis is studied briefly. It is remarked that all these methods are equal-sized in time direction. In this paper, we are devoted to proposing a fully discrete method for linear wave equations with the temporal discretization carried out by the $C^0$-continuous time-stepping methods as used in [21,22] while the spatial discretization by the VEM in [2], respectively. We develop optimal error estimates in the $H^1$ semi-norm and $L^2$ norm by using the energy method, which are equal to $O(h^{k}+\tau^3)$ and $O(h^{k+1}+\tau^3)$, respectively, where $h$ and $\tau$ denote respectively the mesh sizes in space and time, and $k$ denotes the order of the VEM approximation. Our method can use a non-uniform subdivision in time, more flexible in numerically solving linear wave equations. We also offer several numerical experiments to illustrate the performance of the proposed method. Moreover, the numerical comparison indicates that our method proposed performs better than the one in [27] in accuracy.

The rest of this paper is organized as follows. In Section 2, we introduce the model continuous problem and present the VEM semi-discrete method. In Section 3, we propose the $C^0P_2$ time-stepping VEM. Some fundamental results are given in Section 4, followed by error estimates for the method in Section 5. Finally we present some numerical tests in Section 6 to show the performance of the method proposed.

2.   The continuous and discrete problems

2.1.   The continuous problem

Let $\Omega \subset \mathbb{R}^{2}$ be a polygonal domain of interest. The linear wave equation can be described as

where $I: = (0,T)$, $u$ represents the unknown function to be sought, $u_t$ and $u_{tt}$ denote respectively its first and second order time derivative. Assume the source term $f \in {L}^2(\Omega \times {I})$ and the initial data $\Psi_0$, $\Psi_1 \in H^1_0(\Omega)$.

Then the variational formulation of the problem (2.1) reads

where $( \cdot \, , \, \cdot)$ denotes the $L^2(\Omega)$ inner product, $a(u, \, v): = (\nabla u, \, \nabla v)$ for $u$ and $v$ in $H^1_0(\Omega)$, and the partial derivative $u_{tt}$ should be understood in the sense of distribution with respect to $t$.

Define the $H^1(\Omega)$ semi-norm via $| v |^2_1 = a(v,v)$, which is actually a norm on $H^1_0(\Omega)$ by the Poincaré-Friedrichs inequality and is equivalent to the usual $H^1(\Omega)$ norm. The bilinear form $a( \cdot \, , \, \cdot)$ is continuous and coercive, i.e. there exist two uniform positive constants $M$ and $\alpha$ such that

As shown in [24,27], Problem (2.2) has a unique solution $u(t)$ satisfying

for all t in $(0,T)$.

Next, introduce some notation for later requirements. As in [1,10], given a bounded domain $D$ in $\mathbb{R}^d\; (d = 1,2)$ and a non-negative integer $s$, denote by $H^{s}(D)$ the standard Sobolev spaces with the following norm and seminorm:

Moreover, denote by $(\cdot,\cdot)_{\omega}$ the $L^2$ scalar product on an open subset $\omega \subseteq {D}$.

Let $B$ be a Banach space equipped with the norm $\Vert\cdot\Vert_{B}$. For $1\leq p\leq\infty$, define

where

and

We also define

equipped with the norm

To simplify the presentation, write

In our forthcoming analysis, the symbol $C$ stands for a generic positive constant independent of the mesh size $h$ and the time step size $\tau$, which may take different values at different occurrences. And for any two quantities $a$ and $b$, $a \lesssim{b}$ indicates $a \leq {C}b$.

2.2.   The virtual semi-discrete method

Let $\{\mathcal{T}_h\}_h$ be a family of conforming partitions of $\Omega$ into general polygonal polygons $E$, and

For each polygonal element $E$ with $n_E$ edges, denote by $\{P_i\}_{i = 1}^{n_E}$ its vertices which are numbered in a counter-clockwise order and by $e_i$ the edge connecting $P_i$ and $P_{i+1}$, where we identify $P_{n_E+1}$ as $P_1$. The dependence on $E$ will be always omitted when no confusion can arise. In this subsection, we will obtain the VEM semi-discrete problem based on the VEM given in [28,27], with the family of polygonal meshes ${\left\{ {\mathcal{T}_h} \right\}_h}$ satisfying the following condition (cf. [12,16]):

A1. For each $E\in {\mathcal T}_h$, there exists a "virtual triangulation" ${\mathcal T}_E$ of $E$ such that ${\mathcal T}_E$ is uniformly shape regular and quasi-uniform. The corresponding mesh size of ${\mathcal T}_E$ is proportional to $h_E$. Each edge of $E$ is a side of a certain triangle in ${\mathcal T}_E$.

As shown in [16], this condition covers the usual geometric assumptions frequently used in the context of VEMs.

We will adopt two main steps to get the virtual element discretization of Problem (2.2). Firstly, we consider how to construct the local virtual element space and the local discrete bilinear form on each element $E$. Then we obtain the global virtual element space and the global discrete bilinear form based on their local counterparts.

For all natural number $k \geq 1$ and $E \in \mathcal{T}_h$, introduce the following function spaces:

● $\mathbb{P}_k(E)$ is the set of all polynomials on $E$ with the total degree no more than $k$;

● $\mathbb{B}_k(\partial E) : = \{v \in C^0(\partial{E}) ; \, v_{|_e}\in \mathbb{P}_k(e) \; \, \forall \, {e} \subset \partial {E} \, \}$,

where conventionally $\mathbb{P}_{-1}(E) = \{ 0 \}$.

For a broken Sobolev space $H^1(\mathcal{T}_h) : = \Pi_{E \in \mathcal{T}_h} H^1(E)$, define a broken $H^1$ semi-norm by

For all $E \in\mathcal{T}_h$, introduce an augmented local virtual element space $\tilde{V}_k(E)$ by

Then define a local modified virtual element space $W_k(E)$ (which is isomorphic to $\tilde{V}_k(E)$) as

where $\Pi_k^{\nabla,E}$ stands for the $H^1$-projection to $\mathbb{P}_k(E)$, defined later in this subsection, and $\mathbb{M}_{\ell}(D)$ denotes the scaled monomial on a $d$-dimensional domain $D$, i.e.

with $\boldsymbol{x}_D$ the barycenter of $D$, $h_D$ the diameter of $D$ or $h_D = |D|^{1/d}$, and $\ell$ a non-negative integer. Next, we consider the dual space $\mathcal{X}_k(E)$ of $W_k(E)$ defined as

where the functional vectors $\chi_v, \, \chi_e^{k-2}, \, \chi_E^{k-2}$ are obtained via introducing the following set of linear functionals of $W_k(E)$. For all $v \in W_k(E)$ we take (see Figure 1):

● $\chi_v$: Values of $v$ at the $n_E$ internal vertices.

● $\chi_e^{k-2}$: All moments on edges of $E$ up to order $k-2$, i.e.

(alternatively, we can also choose the function values at the $k-1$ internal points of the $(k+1)$-Gauss-Lobatto quadrature rule in $e$, as suggested in [5]).

● $\chi_E^{k-2}$: All moments on elements $E$ up to order $k-2$, i.e.

As shown in [3,6], the set of degrees of freedom (d.o.f.s) given above are unisolvent for $W_k(E)$, i.e. $(W_k(E))' = \mathcal{X}_k(E)$. So the triple $(E, W_k(E), \mathcal{X}_k(E))$ really forms a finite element. Relabel the d.o.f.s by a single index $i = 1,2,\cdots, n_{\rm{dof}} : = \dim {W}_k(E)$. Therefore, with these d.o.f.s we can associate a canonical basis $\{\phi_1, \, \phi_2, \, \cdots, \, \phi_{n_{\rm{dof}}}\}$ of ${W}_k(E)$ such that

The basis $\{ \phi_j \}_{j = 1}^{n_{\rm{dof}}}$ do not need to be written explicitly and this is the reason of the word "virtual" in VEM. Then every function $v \in {W}_k(E)$ can be expanded as

and in numerical computation it can be identified as a vector $\boldsymbol{v} \in \mathbb{R}^{n_{\rm{dof}}}$ in the form

Thus we can establish an isomorphism between functions in $W_k(E)$ and vectors in $\mathbb{R}^{n_{\rm{dof}}}$:

Next, define the local $H^1$ projection operator: $\Pi_k^{\nabla,E} \colon H^1(E) \to \mathbb{P}_k(E)$ as the solution of

where $P^{0,E} \colon H^1(E) \to \mathbb{R}$ is taken as

In order to handle the time term $(u_{tt}(t), \, v)$ in (2.2), we also introduce the local $L^2(E)$ projection operator $\Pi^{0,E}_{k}$ from $L^2(E)$ to $\mathbb{P}_{k}(E)$, defined by

As in [2], we conclude that $\Pi_k^{\nabla, E} = \Pi_k^{0, E}$ for the case $k\leq{2}$.

Using an integration by parts, the right hand side of the first equation of (2.4) can be written as

therefore $\Pi_k^{\nabla,E}v$ is computable with the d.o.f.s of $v\in W_k(E)$. So is the quantity $\Pi^{0, E}_{k}v$.

Since for all $u_h, v_h \in W_k(E)$, the quantities $a^E(u_h,v_h): = (\nabla{u_h},\, \nabla{v_h})_E$ and $(u_h,v_h)_E$ are not computable, we should modify the two local bilinear forms accordingly. So we define two bilinear forms $a_h^E(\cdot, \, \cdot)$ and $m_h^E(\cdot, \, \cdot)$ such that

to approximate the previous bilinear forms.

Note that $a^E(\Pi_k^{\nabla,E}u, \, \Pi_k^{\nabla,E}v)$ (resp. $(\Pi_k^{0,E}u, \Pi_k^{0,E}v)_E$) alone does not lead to a stable method, a stabilization term $\mathcal{S}^E(\cdot, \, \cdot)$ (resp. $\mathcal{R}^E(\cdot, \, \cdot)$) should be added to overcome the difficulty. To ensure the stability while maintaining the accuracy, we should impose the following assumptions on the element-wise stabilization terms $\mathcal{R}^E(\cdot, \, \cdot)$ and $\mathcal{S}^E(\cdot, \, \cdot)$ (cf.[3]).

● $\boldsymbol{k}$-consistency: For all $q\in \mathbb{P}_k(E)$ and $v \in {W}_k(E)$,

● Stability: There exist positive constants $\alpha_*$, $\alpha^*$ independent of $h$ and $E$, such that, for all $\tilde{v} \in (I-\Pi_k^{\nabla,E})W_k(E)$, there holds

Then following the VEM framework, for all $u$ and $v$ in $W_k(E)$, define

Remark 1. VEMs are in fact a family of schemes different in the choice of the stabilization terms. For our problem discussed here, one of the key points is the construction of $\mathcal{S}^E(\cdot, \, \cdot)$ and $\mathcal{R}^E(\cdot, \, \cdot)$. The typical choice for the two terms was given in [2,3,5], which is also suitable under the geometric assumption A1 by using the norm equivalence given in [16]. Concretely speaking, for all $u$ and $v$ in ${W}_k(E)$,

Now we define the global virtual element space by assembling the set of local spaces $W_k(E)$, described as

and its dimension is

where $N_P$ (resp. $N_V$ and $N_e$) is the number of elements (resp. internal vertexes and edges) in $\mathcal{T}_h$, and $N_{\rm{dof}}$ the number of degrees of freedom of $W_h$. The global degrees of freedom for $W_h$ are chosen as follows.

● $\mathbf{G1}$: Values of $v$ at the $n_E$ vertexes of the polygon $E$.

● $\mathbf{G2}$: Values of $v$ at the $k-1$ internal points of the $(k+1)$-Gauss-Lobatto quadrature rule on $e$.

● $\mathbf{G3}$: All moments up to order $k-2$ of $v$ on each element $E$.

Also, define the global approximated bilinear forms $a_h(\cdot, \, \cdot):{W_h}\times{W_h}\to\mathbb{R}$ and $m_h(\cdot, \, \cdot): W_h \times W_h \to \mathbb{R}$ via summing the local counterparts:

It is evident that

Thus we introduce the approximated $H^1$ semi-norm and the approximated $L^2$ norm by

By means of $\Pi_k^{0,E}$, we define the approximated source term $f_h(t)$ for all $t\in (0,T)$ as

Thus we can construct a computable approximation $(f_h(t), \, v)$ for the right-hand side $(f, \, v)$ in (2.2).

As usual, in order to impose the computable initial conditions, introduce the modified $H^1$ projection $\mathcal{P}^{\nabla} \colon H^1_0(\Omega) \to W_h$ defined by

Then the suitable discrete initial data $\Psi_{0,h}$ and $\Psi_{1,h}$ can be chosen as

According to (2.8), (2.9), (2.10), (2.12) and (2.14), we are ready to construct the virtual semi-discrete approximation to Problem (2.2) as follows.

Remark 2. Different from the finite element method, we should use $m_h( u_{h, tt}(t), \, v_h )$ instead of $( u_{h, tt}(t), \, v_h )$, since the latter one is not computable.

For developing error analysis later on, we also introduce the $L^2$-projection operator $\mathcal{P}^{0} \colon L^2(\Omega) \to W_h$ defined by

Referring to [28,27], the following approximation results hold.

Lemma 2.1. For all $u \in H^1_0(\Omega) \cap H^{k+1}(\Omega)$, there holds

Furthermore, if the domain $\Omega$ is convex, we have

Lemma 2.2. If $u \in H^{k+1}(\Omega)$, there holds

3.   The $C^0P_2$ Lg time-stepping VEM

Consider a non-uniform subdivision for the time interval $I$ with the nodes: $0 = t_{0} < t_{1} < \cdots < t_{N} = T$. For $0\leq n\leq N-1$, write $I_{n}: = (t_{n},t_{n+1})$ with time length $\tau_{n} = t_{n+1}-t_{n}$, $\tau: = \max_{0\leq n\leq N-1}\tau_{n}$, and introduce the following notation:

Now, we define the space-time discrete space by

and denote by $S_{n}^{h\tau}$ the restriction of $S^{h\tau}$ to $I_n$.

Based on (2.15) and following some ideas in [21], we propose a space-time virtual element scheme for Problem (2.1) as follows.

where $m_h(\cdot , \, \cdot)$ and $a_h(\cdot, \, \cdot)$ are defined in (2.10) and (2.9), respectively. Here, we also use the notation $\dot{v} : = v_{t} = \partial{v}/\partial{t}$ and $\ddot{v} : = v_{tt} = \partial^{2}v/\partial{t^2}$.

We call the scheme (3.2) as the $C^0P_2$ time-stepping VEM. Next, we present the numerical implementation method (3.2). We know $U\in S^{h\tau}$ is continuous at $t = t_{n}$ and is a second order polynomial in the variable $t$ on the subinterval $I_{n}$, so $U|_{I_{n}}$ is uniquely determined by $U^n$, $U^{n+1}$ and $\dot{U}_{+}^{n}$. Hence, by some direct computation, for $t \in (t_{n},\, t_{n+1})$, the function $U$ and $\dot{U}$ can be expressed as

By (3.4), we have

Then inserting (3.3) and (3.4) into the first equation of (3.2) and taking $\dot{v}$ to be $\phi$ and $(t-t_{n}) \eta$ respectively on $I_n$ where $\phi,\eta \in W_h$, we are able to derive the explicit formulation for the full discrete scheme (3.2), i.e. the following linear system:

Find $\{U\}_{n = 0}^{N}$ and $\{ \dot{U}_{+}^{n} \}_{n = 0}^{N-1}$ in $W_h$ such that

It is easy to see from the method (3.5) that along the time direction, the numerical solution can be computed one interval after another.

Lemma 3.1. In view of similar arguments as given in [21], if $U$ is the solution of (3.2) with $f_h = 0$, and the initial function $U^0$ and $\dot{U}_{-}^0$ arbitrary given, then there hold

and

The first inequality (3.6) ensures the unique solvability of (3.2).

4.   Some fundamental results

We now present some preliminary results useful in the sequel.

4.1.   Error estimates for static problems

Consider the variational formulation of the static problem of (2.1) and its virtual element discretized problem, described as follows.

Following the ideas for deriving error estimates of the VEM in [2] and exploiting some estimates in [16], we can easily obtain the following lemma.

Lemma 4.1. Let $u$ be the solution of (4.1), and let $u_h \in {W}_h$ be the solution of(4.2), with $W_h$ defined in (2.8), $a_h$ defined in (2.9), and $f_h$ defined in (2.12). Assume further that $\Omega$ is convex, the right-hand side $f\in H^k(\Omega)$, and the exact solution $u$ of (4.1) belongs to $H^{k+1}(\Omega)$. Then there holds

4.2.   The backward homogeneous problem

Consider the dual problem of (2.1) with $f = 0$, described as follows.

By introducing a time transformation $\tilde{t} = T - t$, we can rewrite (4.4) as the same form for (2.1). We then use the $C^0P_2$ time stepping VEM (3.2) to solve this problem, and the numerical scheme is equivalent to finding $g\in S^{h\tau}$ such that

where $\Phi_{0,h}$ and $\Phi_{1,h}$ are arbitrary functions in $W_h$. Summing over $n$ from $0$ to ${N-1}$ yields

According to the Lemma 3.1, the function $g\in{S}^{h\tau}$ given by (4.5) satisfies

4.3.   The time-direction interpolation operators $Q_{\tau}$

Along the time direction, introduce the local quadratic interpolation operators $Q_{n}: C^1(\bar{I}_{n}) \to \mathbb{P}_2(I_{n})$ by

for $n = 0, 1, \cdots, N-1$. It is easy to check that

With these operators, we then denote the global interpolation operator $Q_{\tau} : C^1(\bar{I}) \to \prod_{n = 0}^{N-1}(\mathbb{P}_2(I_{n}))$ by

The following lemma presents an error bound for the interpolation operator $Q_{\tau}$ (cf. [15]).

Lemma 4.2. For any function $v\in C^1(\bar{I})$ with $v_{ttt}\in L^p(I)$, there holds

where $1\leq {p} \leq \infty$ and $v_{ttt}: = \partial^{3}v/\partial{t^3}$.

5.   Error estimates for the $C^0P_2$ Lg time-stepping VEM

In this section, we will develop error analysis for the $C^0P_2$ time stepping VEM (3.2). Throughout this section, we assume that the exact solution $u$ of (2.1) satisfies the following regularity conditions

and

Write

where $e_1 : = u - Q_{\tau}(\mathcal{P}^{\nabla}u)$ and $e_2 : = U - Q_{\tau}(\mathcal{P}^{\nabla}u)$. For clarity, we divide our error analysis into three steps.

Step 1. The estimate for $e_1$. By (4.8) and (2.17),

Observing that the time derivative is commutative with the modified $H^1$ projection $\mathcal{P}^{\nabla}$, we have by (4.8) and (2.18) that

Step 2. The estimate for $e_2$. Using (2.1) combined with an integration by parts, we find

Then by the definition of $\mathcal{P}^{\nabla}$ (resp. $\mathcal{P}^{0}$) (2.13) (resp. (2.16)), (5.6) is equivalent to

Subtracting the first equation of (3.2) from (5.7), and using the fact that $[\dot{u}]_{n} = 0$ for $0\leq{n}\leq{N-1}$, we obtain

Furthermore, we write the above equation associated with the initial conditions as follows.

where

Inserting (5.3) into (5.8) and summing these equations over $n$ from $0$ to $N-1$, we obtain the following identity:

We first simplify the LHS of (5.9). We have by integration by parts and a direct manipulation that

The definition of $Q_{\tau}$ and (2.14) implies that

and

Then (5.10) can be reformulated as

Now, in order to simplify the expression, we choose $v$ as $g$, where $g\in S^{h\tau}$ is the solution of (4.5). Then by using (4.6) we know the LHS of (5.9) satisfies

Next, we try to simplify the RHS of (5.9):

Since $\mathcal{P}^{0}u \in C^1(0, T; \,W_h)$ and $\mathcal{P}^{\nabla}u \in C^1(0, T; \,W_h)$, we have by (4.8) and (4.9) that

Then using integration by parts and observing the fact $\ddot{g}$ is piecewise constant in the time direction, we see that

Therefore, in terms of (5.13) and (5.14), we can rewrite the RHS of (5.9) as

Hence, the combination of (5.9), (5.12) and (5.15) together implies

Now, choosing $g$ such that $g^N = (e_2)^N$ and $\dot{g}_{+}^N = 0$, from the stability of $a_h(\cdot, \, \cdot)$ and (5.16), we have

where

Let us study the first term in the sum (5.17). By using the Cauchy-Schwarz inequality, the triangle inequality, (4.7), (2.18), (2.19) and the consistency property of the approximated $H^1$ semi-norm and $L^2$ norm, we have

For the second term, since the time derivative is commutative with $\mathcal{P}^{\nabla}$ (resp. $\mathcal{P}^{0}$), we have by integration by parts that

Hence, according to the Cauchy-Schwarz inequality, the coercivity of $m_h(\cdot, \, \cdot)$, (4.7) and (4.10),

For the last term, we find by the Cauchy-Schwarz inequality and (2.12) that

Finally, inserting the estimates (5.18), (5.19) and (5.20) into (5.17), we arrive at

If we take

in (5.16) and then use the similar arguments for deriving (5.21), we can obtain

Step 3. The estimate for $e$. Combining (5.4), (5.5), (5.21) and (5.22) obtained in Steps 1 and 2 and using the triangle inequality

where $\| \cdot \|$ stands for the norm $| \cdot |_{1,\Omega}$ or $\| \cdot \|_{0,\Omega}$, we can readily derive the desired error estimates for the numerical method (3.2), described as the following theorem.

Theorem 5.1. Let $u$ and $U$ be the solutions of (2.1) and (3.2), respectively. Assume $u$ satisfies the regularity conditions (5.1) and (5.2). Then for $n = 1,2, \cdots, N$, there hold

and

Corollary 1. Under the conditions in Theorem 5.1, for $n = 1,2, \cdots, N$, there holds

Remark 3. For our method (3.2), if choosing $k = 2$ and adopting uniformly temporal discretization, we have from the above estimate that $\| u(t_n) - U^n \|_{0,\Omega} = O(h^3+\tau^3)$.

6.   Numerical results

In this section, we present several numerical experiments on the operating platform (Matlab R2016a on MacOS10.12.6) to show the performance of the $C^0P_2$ time-stepping VEM (3.2). The codes are designed based on the references [5,25] and polygonal Voroni meshes $\mathcal{V}_h$ are produced using the software Polymesher introduced in [26].

Theoretically, the error generated by a fully discrete scheme has two components: the error due to the spatial discretization depending on $h$, and the error created by the time integrator depending on $\tau$. In particular, let $\{U^n\}_{n = 0}^N$ be the function sequence generated by the method (3.2), where $U^n = U(\cdot,t_n)$ for $n = 0,1,\dots, N$.

Define

where $\boldsymbol{K}$ and $\boldsymbol{M}$ are the related stiffness and mass matrices, respectively.

And introduce the experimental order of convergence:

where $Err$ and $Err'$ are the errors on two consecutive grids with space mesh-size $h$ and $h'$ or time division-size $\tau$ and ${\tau}'$, respectively.

We use the approximated $H^1$ semi-norm and $L^2$ norm to study the orders of the error $E_{H^1}$, $Et_{L^2}$ and $E_{L^2}$ with respect to $h$ and $\tau$, respectively.

Example 6.1. Consider the problem (2.1) defined on the unit square $\Omega = (0,1)^2$ and choose the time interval $I = (0,1)$. We take the load term $f$, the initial data $\Psi_0$ and $\Psi_1$ to be in accordance with the exact solution

We will choose the orders of the VEM approximation as $k = 1$ and $k = 2$. For simplicity, we adopt the $C^0P_2$ element with the uniform division in the time discretization.

Case 1. $k = 1$

Firstly, in order to observe the convergence order of error $E_{H^1}$, $Et_{L^2}$ in space direction, we fix time-division size $\tau = 1/40$ and let $h$ vary from $1/5$ to $1/80$. The numerical results shown in Table 1 and Figure 2 in the log scale indicate that $E_{H^1}$ is $1$st-order and $Et_{L^2}$ is $2$nd-order in space direction.

Then we consider the order of error $E_{H^1}$, $Et_{L^2}$ with respect to $\tau$. In this case, we take $h = {\tau}^3$ and $h^2 = {\tau}^3$ so that we can decrease the heavy computational cost but still illustrate the need results. For $E_{H^1}$, we choose $\tau = 1/2, 1/3, 1/4, 1/5, 1/6$ with the space mesh size $h = {\tau}^3$ in scheme (3.2). Similarly, for $Et_{L^2}$, we set $\tau = 1/4, 1/6, 1/8, 1/10, 1/12, 1/14, 1/16$ with the space mesh size $h^2 = {\tau}^3$. From Table 2 and 3, and Figure 3, we can see that both $E_{H^1}$ and $Et_{L^2}$ in our method proposed in (3.2) have $3$-rd order accuracy in time direction.

Finally, we consider the error bound of the $C^0P_2$ time stepping VEM proposed in (3.2). We take $h = 1/2^{j+1}, \; j = 1,...,6$ with $h = \tau$ to generate different meshes, and show the numerical results in Table 4 and Figure 6. From the $4$th and $7$th column of Table 4, we observe that the ratios between $E_{H^1}$ (resp. $Et_{L^2}$) and $h+{\tau}^3$ (resp. $h^2+{\tau}^3$) can be controlled a bounded quantity, which indicates that there exists some absolute positive constant $C$ such that

In addition, from the $5$th and $8$th column of Table 4 and Figure 6, we can see that the relative errors are also decreasing when mesh sizes are refined.

Case 2. $k = 2$

In this case, we perform the similar experiments as for $k = 1$. Firstly, in order to observe the convergence order of error $E_{H^1}$, $Et_{L^2}$ in space direction, we fix time-division size $\tau = 1/80$ and let $h$ vary from $1/5$ to $1/80$. The numerical results shown in Table 5 and Figure 4 in the log scale indicate that $E_{H^1}$ is $2$nd-order and $Et_{L^2}$ is $3$rd-order in space direction.

Then we consider the order of error $E_{H^1}$, $Et_{L^2}$ with respect to $\tau$. In this case, we take $h^2 = {\tau}^3$ and $h = {\tau}$ so that we can decrease the heavy computational cost but still illustrate the needed convergence orders. For $E_{H^1}$, we choose $\tau = 1/4, 1/6, 1/8, 1/10, 1/12,1/14,1/16$ with the space mesh size $h^2 = {\tau}^3$ in scheme (3.2). Similarly, for $Et_{L^2}$, we set $\tau = 1/2, 1/4, 1/8, 1/16, 1/32, 1/64, 1/128$ with the space mesh size $h = {\tau}$. From Table 6 and 7, and Figure 5, we can see that both $E_{H^1}$ and $Et_{L^2}$ in our method proposed in (3.2) have $3$-rd order accuracy in time direction.

Finally, we consider the error bound of the $C^0P_2$ time stepping VEM proposed in (3.2). We take $h = 1/2^{j+1}, \; j = 1,...,6$ with $h = \tau$ to generate different meshes, and show the numerical results in Table 8 and Figure 6. From the $4$th and $7$th column of Table 8, we observe that the radio between $E_{H^1}$ (resp. $Et_{L^2}$) and $h^2+{\tau}^3$ (resp. $h^3+{\tau}^3$) can be controlled a bounded quantity, which indicates that there exists some absolute positive constant $C$ such that

In addition, from the $5$th and $8$th column of Table 8, we can see that the relative errors are also decreasing when mesh sizes are refined.

All the numerical results for $k = 1$ and $k = 2$ are in agreement with the theoretical results in Theorem 5.1.

Example 6.2. In this example, we make a numerical comparison between our method (3.2) and the Newmark trapezoidal VEM (see [27]). Assume the functions $f$, $\Psi_0$ and $\Psi_1$ in the problem (2.1) are given such that it has the following exact solution:

We use the Voronoi meshes for the solution domain $\Omega = (0,1)^2$ and the uniform meshes for the time interval $I: = (0,1)$. The chosen orders of the VEM approximation are $k = 1$ and $k = 2$. The results of the approximated $L^2$-norm error $E_{L^2}$ are shown in Table 9.

From Table 9 we can observe that for both cases, the two methods have the expected order of convergence in $h$, that is, $O(h^{k+1})$. However, if we consider the error effect from temporal discretization together, we can find from the $4$-th, $5$-th, $7$-th and $8$-th columns of Table 9 that the error $E_{L^2}$ is actually $O(h^{k+1} +{\tau}^3)$ for the $C^0P_2$ time-stepping VEM while $O(h^{k+1} +{\tau}^2)$ for the Newmark trapezoidal VEM. Therefore, our method performs better than the latter one in the time direction.

Acknowledgments

The authors would like to thank the referees for their valuable comments leading to an improvement of the early version of the paper.