A CP2 TIME-STEPPING VIRTUAL ELEMENT METHOD FOR LINEAR WAVE EQUATIONS ON POLYGONAL MESHES

This paper is concerned with a CP2 time-stepping virtual element method (VEM) for solving linear wave equations on polygonal meshes. The spatial discretization is carried out by the VEM while the temporal discretization is obtained based on the CP2 time-stepping approach, leading to a fully discrete method. The error estimates in the H1 semi-norm and L2 norm are derived after detailed derivation. Finally, the numerical performance and efficiency of the proposed method is illustrated by several numerical experiments.


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 0continuous 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 seminorm and L 2 norm by using the energy method, which are equal to O(h k + τ 3 ) and O(h k+1 +τ 3 ), respectively, where h and τ 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 0 P 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. 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 ∈ L 2 (Ω × I) and the initial data Ψ 0 , Ψ 1 ∈ H 1 0 (Ω). Then the variational formulation of the problem (2.1) reads where (· , ·) denotes the L 2 (Ω) inner product, a(u, v) := (∇u, ∇v) for u and v in H 1 0 (Ω), and the partial derivative u tt should be understood in the sense of distribution with respect to t.
Define the H 1 (Ω) semi-norm via |v| 2 1 = a(v, v), which is actually a norm on H 1 0 (Ω) by the Poincaré-Friedrichs inequality and is equivalent to the usual H 1 (Ω) norm. The bilinear form a(· , ·) is continuous and coercive, i.e. there exist two uniform positive constants M and α such that . As shown in [24,27], Problem (2.2) has a unique solution u(t) satisfying a(u(t), u(t)) + u t (t) 2 for all t in (0, T ). Next, introduce some notation for later requirements. As in [1,10], given a bounded domain D in 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 (·, ·) ω the L 2 scalar product on an open subset ω ⊆ D.
In our forthcoming analysis, the symbol C stands for a generic positive constant independent of the mesh size h and the time step size τ , which may take different values at different occurrences. And for any two quantities a and b, a b indicates a ≤ Cb.

2.2.
The virtual semi-discrete method. Let {T h } h be a family of conforming partitions of Ω into general polygonal polygons E, and For each polygonal element E with n E edges, denote by {P i } n E i=1 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 semidiscrete problem based on the VEM given in [28,27], with the family of polygonal meshes {T h } h satisfying the following condition (cf. [12,16]): A1. For each E ∈ T h , there exists a "virtual triangulation" T E of E such that T E is uniformly shape regular and quasi-uniform. The corresponding mesh size of T E is proportional to h E . Each edge of E is a side of a certain triangle in 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 ≥ 1 and E ∈ T h , introduce the following function spaces: • P k (E) is the set of all polynomials on E with the total degree no more than k; For a broken Sobolev space For all E ∈ T h , introduce an augmented local virtual element spaceṼ k (E) bỹ Next, we consider the dual space X k (E) of W k (E) defined as

Then define a local modified virtual element space
where the functional vectors χ v , χ k−2 e , χ k−2 E are obtained via introducing the following set of linear functionals of W k (E). For all v ∈ W k (E) we take (see Figure  1): • χ v : Values of v at the n E internal vertices.
• χ k−2 e : 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]). • χ k−2 E : 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)) = X k (E). So the triple (E, W k (E), X k (E)) really forms a finite element. Relabel the d.o.f.s by a single index i = 1, 2, · · · , n dof := dim W k (E). Therefore, with these d.o.f.s we can associate a canonical basis {φ 1 , φ 2 , · · · , φ n dof } of W k (E) such that χ i (φ j ) = δ ij , i, j = 1, 2, · · · , n dof . The basis {φ j } n dof j=1 do not need to be written explicitly and this is the reason of the word "virtual" in VEM. Then every function v ∈ W k (E) can be expanded as and in numerical computation it can be identified as a vector v ∈ R n dof in the form Thus we can establish an isomorphism between functions in W k (E) and vectors in R n dof : Next, define the local H 1 projection operator: where P 0,E : H 1 (E) → R is taken as In order to handle the time term (u tt (t), v) in (2.2), we also introduce the local As in [2], we conclude that Π ∇,E k = Π 0,E k for the case k ≤ 2. Using an integration by parts, the right hand side of the first equation of (2.4) can be written as are not computable, we should modify the two local bilinear forms accordingly. So we define two bilinear forms a E h (·, ·) and does not lead to a stable method, a stabilization term S E (·, ·) (resp. R E (·, ·)) 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 R E (·, ·) and S E (·, ·) (cf. [3]).
• k-consistency: For all q ∈ P k (E) and v ∈ W k (E), • Stability: There exist positive constants α * , α * independent of h and E, such that, for allṽ 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 S E (·, ·) and R E (·, ·). 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 T h , and N dof the number of degrees of freedom of W h . The global degrees of freedom for W h are chosen as follows.
• G1: Values of v at the n E vertexes of the polygon E.
• G2: Values of v at the k − 1 internal points of the (k + 1)-Gauss-Lobatto quadrature rule on e.
• G3: All moments up to order k − 2 of v on each element E.
Also, define the global approximated bilinear forms a h (·, ·) : It is evident that Thus we introduce the approximated H 1 semi-norm and the approximated L 2 norm by As usual, in order to impose the computable initial conditions, introduce the modified H 1 projection P ∇ : Then the suitable discrete initial data Ψ 0,h and Ψ 1,h can be chosen as (2.14) 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.
, since the latter one is not computable.
For developing error analysis later on, we also introduce the L 2 -projection operator P 0 : L 2 (Ω) → W h defined by (2. 16) Referring to [28,27], the following approximation results hold.
3. The C 0 P 2 time-stepping VEM. Consider a non-uniform subdivision for the time interval I with the nodes: 0 = t 0 < t 1 < · · · < t N = T . For 0 ≤ n ≤ N − 1, write I n := (t n , t n+1 ) with time length τ n = t n+1 − t n , τ := max 0≤n≤N −1 τ n , and introduce the following notation: and denote by S hτ n the restriction of S hτ 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.
We call the scheme (3.2) as the C 0 P 2 time-stepping VEM. Next, we present the numerical implementation method (3.2). We know U ∈ S hτ is continuous at t = t n and is a second order polynomial in the variable t on the subinterval I n , so U | In is uniquely determined by U n , U n+1 andU n + . Hence, by some direct computation, for t ∈ (t n , t n+1 ), the function U andU can be expressed as 3) 2) and takingv to be φ and (t − t n )η respectively on I n where φ, η ∈ W h , we are able to derive the explicit formulation for the full discrete scheme (3.2), i.e. the following linear system: (3.5) 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 andU 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.  By introducing a time transformationt = T − t, we can rewrite (4.4) as the same form for (2.1). We then use the C 0 P 2 time stepping VEM (3.2) to solve this problem, and the numerical scheme is equivalent to finding g ∈ S hτ such that where Φ 0,h and Φ 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 ∈ S hτ given by (4.5) satisfies ġ ∞,0,h (|Φ 0,h | 1,h + Φ 1,h 0,h ). (4.7) 4.3. The time-direction interpolation operators Q τ . Along the time direction, introduce the local quadratic interpolation operators Q n : C 1 (Ī n ) → P 2 (I n ) by for n = 0, 1, · · · , N − 1. It is easy to check that With these operators, we then denote the global interpolation operator Q τ : The following lemma presents an error bound for the interpolation operator Q τ (cf. [15]).

5.
Error estimates for the C 0 P 2 time-stepping VEM. In this section, we will develop error analysis for the C 0 P 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 u tt ∈ L 1 (0, T ; H k+1 (Ω)), u ttt ∈ L 1 (0, T ; H 2 (Ω)).
where e 1 := u − Q τ (P ∇ u) and e 2 := U − Q τ (P ∇ 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 P ∇ , 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 P ∇ (resp. 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 [u] n = 0 for 0 ≤ n ≤ N − 1, we obtain Furthermore, we write the above equation associated with the initial conditions as follows. Inserting (5.3) into (5.8) and summing these equations over n from 0 to N − 1, we obtain the following identity: In We first simplify the LHS of (5.9). We have by integration by parts and a direct manipulation that The definition of Q τ and (2.14) implies that Then (5.10) can be reformulated as Now, in order to simplify the expression, we choose v as g, where g ∈ S hτ is the solution of (4.5). Then by using (4.6) we know the LHS of (5.9) satisfies In Next, we try to simplify the RHS of (5.9): Since P 0 u ∈ C 1 (0, T ; W h ) and P ∇ u ∈ C 1 (0, T ; W h ), we have by (4.8) and (4.9) that Then using integration by parts and observing the factg 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ġ N + = 0, from the stability of a h (·, ·) and (5.16), we have n=0 In 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 h k+1 utt L 1 (0,T ; H k+1 (Ω)) |(e2) N |1. (5.18) For the second term, since the time derivative is commutative with P ∇ (resp. P 0 ), we have by integration by parts that For the last term, we find by the Cauchy-Schwarz inequality and (2.12) that In In where · stands for the norm | · | 1,Ω or · 0,Ω , 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, · · · , N , there hold and + τ 3 u ttt L 1 (0,T ; H 2 (Ω)) . 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 0 P 2 time-stepping VEM (3.2). The codes are designed based on the references [5,25] and polygonal Voroni meshes 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 τ . In particular, let {U n } N n=0 be the function sequence generated by the method (3.2), where U n = U (·, t n ) for n = 0, 1, . . . , N . Define where K and 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 τ and τ , 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 τ , respectively. Example 6.1. Consider the problem (2.1) defined on the unit square Ω = (0, 1) 2 and choose the time interval I = (0, 1). We take the load term f , the initial data Ψ 0 and Ψ 1 to be in accordance with the exact solution u(x 1 , x 2 , t) = sin(πx 1 )sin(πx 2 )sin(t 2 ).
We will choose the orders of the VEM approximation as k = 1 and k = 2. For simplicity, we adopt the C 0 P 2 element with the uniform division in the time discretization.
Firstly, in order to observe the convergence order of error E H 1 , Et L 2 in space direction, we fix time-division size τ = 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 1st-order and Et L 2 is 2nd-order in space direction. Then we consider the order of error E H 1 , Et L 2 with respect to τ . In this case, we take h = τ 3 and h 2 = τ 3 so that we can decrease the heavy computational cost but still illustrate the need results. For E H 1 , we choose τ = 1/2, 1/3, 1/4, 1/5, 1/6 with the space mesh size h = τ 3 in scheme (3.2). Similarly, for Et L 2 , we set τ = 1/4, 1/6, 1/8, 1/10, 1/12, 1/14, 1/16 with the space mesh size h 2 = τ 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 0 P 2 time stepping VEM proposed in (3.2). We take h = 1/2 j+1 , j = 1, ..., 6 with h = τ to generate different meshes, and show the numerical results in Table 4 and Figure 6. From the 4th and 7th column of Table 4, we observe that the ratios between E H 1 (resp. Et L 2 ) and h + τ 3 (resp. h 2 + τ 3 ) can be controlled a bounded quantity, which indicates that there exists some absolute positive constant C such that  In addition, from the 5th and 8th column of Table 4 and Figure 6, we can see that the relative errors are also decreasing when mesh sizes are refined.
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 timedivision size τ = 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 2nd-order and Et L 2 is 3rd-order in space direction.
In addition, from the 5th and 8th 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 , Ψ 0 and Ψ 1 in the problem (2.1) are given such that it has the following exact solution: u(x 1 , x 2 , t) = e −t (x 2 1 − x 1 )(x 2 2 − x 2 ). We use the Voronoi meshes for the solution domain Ω = (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 + τ 3 ) for the C 0 P 2 time-stepping VEM while O(h k+1 + τ 2 ) for the Newmark trapezoidal VEM. Therefore, our method performs better than the latter one in the time direction.