Numerical resolution of the wave equation using the spectral method

We present a new procedure for the numerical study of the wave equation. We use the spectral discretization method associated with the Euler scheme for spatial and temporal discretization. A detailed numerical analysis leads to an a priori error estimate. We confirm the high precision of the method presented by a numerical study.

The time and space approximation for the wave equation has been studied in [9] using the finite element method. We present a new study of this problem using the spectral method associated with backward Euler discretization scheme. The spectral method is known for its high precision [10,11]. For the space discretization, the discrete spaces are constructed from spaces of polynomials of high degree. Then the discrete problem is obtained using the Galerkin method combined with numerical integration.
The outline of the paper is as follows. In Sect. 2, we study the continuous problem and present some energy estimate properties. In Sect. 3, we are interested to the timesemidiscrete problem. We discretize the second time derivative by using a second difference quotient of the solution on a nonuniform temporal grid. Then we obtain an optimal a priori error estimate. In Sect. 4, we study the fully discrete problem and establish an optimal a priori error estimate. Finally, in Sect. 5, we present a numerical study. sure in the space H 1 ( ) of infinitely differentiable functions with compact support in , and H -1 ( ) is its dual space. We denote by (·, ·) and · , respectively, the scalar product and associated norm in L 2 ( ). By H 1 2 (∂ ) we denote the trace space of functions in H 1 ( ). Let γ ⊂ ∂ , H 1 2 00 (γ ) be the space of functions in H 1 2 (γ ) such that their extension by zero to ∂ /γ belongs to H 1 2 (∂ ). In the following, we define • u(x, t) on × ]0, T[ as where X is a separable Banach space. • C j (0, T; X) represents the set of functions of time of class C j with values in X. It is a Banach space with norm where ∂ l t u is the time partial derivative of u of order l.
• H s (0, T; X) = {v ∈ L 2 (0, T; X); ∂ k v ∈ L 2 (0, T; X); k ≤ s} is the Hilbert space with scalar product • W m,1 (0, T, X) is the space of functions in L 1 (0, T, X) such that all their derivatives up to the order m belong to L 1 (0, T, X). Consider the wave equation problem where the wave u is the unknown defined on × ]0, T[, and (u 0 , v 0 ) are the data functions defined on .
This problem can be written in a more general form, Proof 1 Taking the inner product of the first equation of problem (2) andu v and integrating by parts the second term, we obtain By integrating this inequality between 0 and t we get estimate (3).
To formulate the time semidiscrete problem, we apply the Euler implicit method to system (1). Then it consists in finding the sequence of functions (u k ) 0≤k≤K in the space We suppose that the data (u 0 , v 0 ) ∈ H 1 0 ( ) × L 2 ( ). Then, if u 0 and v 0 are known, then we easily show that u k+1 , k ≥ 1, is a solution of the following weak problem: Proof 2 Using the Lax-Milgram theorem, we show that problem (6) has a unique solution. Then, by iteration on k, we deduce that problem (5) has a unique solution.
Taking the inner product of u k+1 -u k δt k and the first equation in (5), we obtain Applying the Cauchy-Schwarz inequality leads to Then by iteration on k we have Finally, we conclude by using the third and fourth equations of system (5).
Remark 2 1) We notice that for k ≥ 1, the solution u k+1 of problem (6) belongs to H s+1 ( ), s ≥ 1 2 . When the domain is convex or of dimension 1, s ≥ 1 is explicitly known. For 1 2 ≤ s ≤ 1, from condition (7) we derive the inequality where the constant C is independent of the step δt. This inequality is not optimal since u k+1 2 is not bounded independently of the step δt.
2) Using the implicit Euler scheme for the time discretization, problem (2) is written as follows: where F k+1 = f k+1 g k+1 . For n = 0, systems (11) and (5) coincide if F k+1 = 0, k ≥ 1. When n = 0, we propose the following two cases where the two systems completely coincide: i) We replace the fourth equation of system (5) by the following implicit equation: ii) We replace the third equation of system (11): Multiplying the first equation in (13) byu k+1 u k+1 , we obtain the following stability condition: However, if we take the inner product of the first equation of system (13) and u k+1 ( ) -1 u k+1 , we obtain the following stability condition in terms of weaker norms: To obtain the error estimate between the solutions u of (1) and (u) k 0≤k≤K of (5), we define the error We can easily show that (ϒ) k 0≤k≤K is a solution of (13), where the two components of F k+1 are the following consistency errors: where C is a positive constant independent of δt. (13), where the second member is F k+1 . Then applying the stability condition (14) leads to thanks to the Taylor theorem with integral remainder to bound the terms ε(v) j , ∇ε(u) j , e(v) 0 , and ∇e(u) 0 . Then we conclude the desired estimate (17).
We can find the following error estimate in weaker norms by using the same technique as the proof of Theorem 1 and replacing condition (14) by (15).

Corollary 1
Suppose that the solution u of system (1) belongs to W 3,1 (0, T; L 2 ( )) ∩ W 2,1 (0, T; H 1 0 ( )). Then the following a priori error estimate between the solution u and the solution (u) k 0≤k≤K of system (5) holds for 0 ≤ k ≤ K : where C is a positive constant independent of δt.
We remark that the obtained estimates (18) and (19) are optimal of order 1 in time, since the discretization is based on the implicit Euler scheme, which is of order 1.

The spectral discretization
We further suppose that is a rectangle for d = 2 or a parallelepiped rectangle for d = 3.
Let P N ( ) the space of polynomials of degree ≤ N (N ≥ 2) for each variable, and let P 0 We recall the following property (see [10,13]): The reference domain ]-1, 1[ d (d = 2, 3) is transformed to the domain using the affine mapping T, and the scalar product is defined on continuous functions u and v by Remark 3 For simplicity of analysis, we suppose that the spectral discretization is fixed over time.
We suppose that u 0 and v 0 are continuous on¯ . The discrete problem is deduced from (5) by applying the Galerkin method combined with numerical integration: where I N is the interpolating operator from L 2 ( ) into P N ( ), find u k As in (6), u k+1 N , 1 ≤ k ≤ K , is the solution of the discrete weak problem Proposition 3 Let the data (u 0 , v 0 ) ∈ H 1 0 ( ) × L 2 ( ). If u 0 N and v 0 N are known, then problem (25) has a unique solution u k+1 N , k ≥ 1, in H 1 0 ( ). Moreover, the solution (u k N ) 0≤k≤K of problem (23)-(24) satisfies for 0 ≤ k ≤ K the following stability condition: Proof 4 We show that problem (25) has a unique solution using the Lax-Milgram theorem and property (21).
To prove the stability condition (26), we define · d the discrete norm deduced from the discrete scalar product (·, ·) N . Now letting v N = Using the Cauchy-Schwarz inequality and (21), we have Then iterating over k, we obtain Finally, estimate (26) is deduced from (23).

Proposition 4
Let u 0 , v 0 be continuous on , and let u 0 N , v 0 N be known. The error estimate between solutions u k+1 , k ≥ 1, and u k+1 N , k ≥ 1, of problems (6) and (25), respectively, is where and C is a positive constant independent of N .
Proof 5 Consider χ k+1 N ∈ P 0 N ( ). By the triangle inequality we have To estimate u k+1 Nχ k+1 N , we begin by writing problems (5) and (25) for v N ∈ P 0 N ( ). Then we consider τ k = δt k δt k-1 and doing the difference term by term, we obtain Since K k is linear and continuous on P 0 N ( ), by the Riesz theorem there exists ϑ k N in P 0 N ( ) such that Applying the result proved in [14,Prop. 4.1] and [15], we get where C is a positive constant independent of N . So we conclude (27), since where C is a positive constant independent of N .
To find the order of convergence as a function of N , it is necessary to estimate each of the terms of the second member of inequality (27).
• Estimation of T 1,j We consider j+1 = u j+1u j and χ where 1,0 N is the orthogonal projection operator from H 1 0 ( ) into P 0 N ( ) related to the inner product defined by the semi norm | · | 1, . (See ( [13], Lemma VI.2.5) and [10] for all the properties of this operator.) • Estimation of T 2,j Since the Gauss-Lobatto quadrature formula is exact for a polynomial of degree ≤ 2N -1, we have Thanks to the triangle and Cauchy-Schwarz inequalities, we have Thus we conclude using the properties of 1,0 N-1 .
• Estimation of T 3,j Let θ j = u ju j-1 . We use for this estimation N-1 the orthogonal projection from L 2 ( ) into P N-1 ( ). By the same argument above for the Gauss-Lobatto formula, we have Using inequality (21) in each direction leads to Thanks to the approximation properties of operator N-1 (see [10,Thm. 7.1]) and I N (see [10,Thm. 14.2]), for θ j ∈ H s ( ); s > 1, we obtain Finally, to estimate where C is a positive constant independent of N .

Numerical results
Consider the interpolating Lagrange polynomial ϕ j defined by where δ ij is the Kronecker symbol. The solution u k+1 N of problem (25) is written as Let U k+1 be the admissible solution vector of components u k+1 N (ζ i , ζ j ). The matrix system deduced form of the discrete problem (25) is written where D k+1 is a diagonal matrix of coefficients r s , 1 ≤ r, s ≤ N -1, A k+1 is the matrix with coefficients (∇(ϕ i ϕ j ); ∇(ϕ r ϕ s )), 1 ≤ i, j, r, s ≤ N -1, and F k is the vector with components We remark that the matrix D k+1 + δt 2 k A k+1 is symmetric and positive definite. Then system (34) is solved using the gradient conjugate method at each iteration.

Iterative algorithm
Step 1: Let Step 2: Suppose u k-1 N and u k N are known The linear system (34) is solved by the gradient conjugate method We do the iterations until the following condition is satisfied: where ξ = 10 -10 for all the following numeric tests.
(37) Figure 3 (respectively, Fig. 4) deals with log 10 uu Nδt H 1 ( ) (in blue) and log 10 uu Nδt L 2 ( ) (in red) as functions of N (respectively, log 10 (N)). We remark that the error norms log 10 uu Nδt H 1 ( ) decrease exponentially until N = 10 and stagnate for N > 10. The errors log 10 uu Nδt L 2 ( ) decrease until N = 10 and stagnate for N > 10. We remark that convergence stagnates. This is due to the fact that the time order of convergence is less than the order of the spectral method. Finally, the isovalues of the exact and discrete solutions (37) are presented in Fig. 5.

Conclusion and future work
This work concerns the numerical analysis of the implicit Euler scheme in time and the spectral discretization in space of the second-order wave equation. We prove an optimal error estimate in time and space. These estimates depend only on the regularity of the solution. Although the spectral methods are known as high-order methods in space, we remark that they have the disadvantage of losing part of this accuracy due to lower order of temporal discretization (often of order 1 or 2). The numerical analysis and implementation of the more difficult case where the spectral discretization depends on time using the second-order BDF method for the time discretization will be the subject of our forthcoming work.