# On properties of an explicit in time fourth-order vector compact scheme for the multidimensional wave equation

An initial-boundary value problem for the n-dimensional wave equation is considered. A three-level explicit in time and conditionally stable 4th-order compact scheme constructed recently for n=2 and the square mesh is generalized to the case of any n≥ 1 and the rectangular uniform mesh. Another approach to approximate the solution at the first time level (not exploiting high-order derivatives of the initial functions) is suggested. New stability bounds in the mesh energy norms and the discrete energy conservation laws are given, and the 4th order error bound is rigorously proved. Generalizations to the cases of the non-uniform meshes in space and time as well as of the wave equation with variable coefficients are suggested.

## Authors

• 4 publications
11/28/2020

### On compact 4th order finite-difference schemes for the wave equation

We consider compact finite-difference schemes of the 4th approximation o...
01/26/2021

### On Properties of Compact 4th order Finite-Difference Schemes for the Variable Coefficient Wave Equation

We consider an initial-boundary value problem for the n-dimensional wave...
05/06/2021

### A mollifier approach to regularize a Cauchy problem for the inhomogeneous Helmholtz equation

The Cauchy problem for the inhomogeneous Helmholtz equation with non-uni...
10/12/2019

### On the numerical approximations of the periodic Schrödinger equation

We consider semidiscrete finite differences schemes for the periodic Scr...
12/02/2020

### A compact higher-order finite-difference scheme for the wave equation can be strongly non-dissipative on non-uniform meshes

We study necessary conditions for stability of a Numerov-type compact hi...
09/17/2011

### A KdV-like advection-dispersion equation with some remarkable properties

We discuss a new non-linear PDE, u_t + (2 u_xx/u) u_x = epsilon u_xxx, i...
11/04/2019

### Uniform observability of the one-dimensional wave equation for non-cylindrical domains. Application to the control's support optimization

This work is concerned with the distributed controllability of the one-d...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The higher-order compact finite-difference schemes form an important class of numerical methods to solve PDEs due to their much greater computational efficiency compared to more standard 2nd order schemes. In the last decade, the construction and study of higher-order compact schemes for initial-boundary value problems for the 2D and 3D wave or telegraph equations, including the case of variable coefficients, have attracted great interest. Such compact schemes are implicit and most often conditionally stable. Among them, the first type of schemes is constituted by implicit schemes which require application of FFT (for constant coefficients) or iterative methods (for variable coefficients) for their efficient implementation, see, in particular, [2, 8, 18, 24, 22], where additional references are contained.

The second type of compact schemes belongs to the alternating direction implicit (ADI)-like methods (for example, see [15, 17]) which implementation is direct, reduces to solving collections of 1D linear algebraic tridiagonal systems and is much faster. The three-level in time compact ADI schemes of several types with the approximation order were suggested for the 2D wave or telegraph equations with constant coefficients, in particular, in [9, 11, 14]; for the 3D and multidimensional cases, see also [9, 16]. Here and are mesh steps in time and space. More important compact ADI schemes with the 4th approximation order, i.e. , for such 2D and 3D equations were constructed and studied, in particular, in [4, 5, 12, 19, 24]. For the variable coefficient 2D wave equation, the compact ADI scheme having the approximation order has recently been constructed in [3], and the 4th order compact ADI-like schemes were proposed and studied in the 2D and 3D cases in [10, 13, 23].

In the recent paper [7], a new third type of compact schemes has been suggested for the 2D wave equation in a square, with the non-homogeneous Dirichlet boundary conditions. The scheme is three-level explicit in time and conditionally stable; the uniform mesh in time and the square spatial mesh are taken. Its construction begins with applying the usual three-point 1D Numerov scheme in time but then, importantly, not only the solution to the wave equation is approximated as usual but its second spatial derivatives are also approximated independently as auxiliary sought functions. For these auxiliary functions, the three-point 1D Numerov schemes in the corresponding space directions are used. Thus the scheme is not completely explicit (as the classical explicit 2nd order scheme) since the auxiliary functions need to be computed by solving collections of 1D linear algebraic systems with tridiagonal matrices. In this respect, the scheme implementation and its cost are similar to the ADI methods, but the difference is that the auxiliary functions are completely independent and can be computed in parallel. Theoretically, only the spectral stability analysis on harmonics was accomplished in the case of zero free terms in the equations of the scheme, but stability bounds and error bounds were not given. The successful results of numerical experiments were presented. Note that other two-level vector scheme with the first time derivative as the additional sought function are well-known, for example, see [1, 8, 25, 26].

In this paper, the scheme from [7] is generalized to the case of the -dimensional wave equation,

, and the uniform rectangular mesh (recall that the square mesh cannot be constructed in any rectangular domain), and the scheme is classified as a vector one more systematically. The new explicit two-level vector equations for the approximate solution at the first time level are constructed too. They are similar to the main equations and exploit only second order finite-differences of the initial functions. Notice that we do not use derivatives of the free term and initial functions in contrast to

[7] that is essential for applications of the scheme in the non-smooth case, for example, see [8, 10, 24, 22]. The scheme is new even in the simplest 1D case.

Next, new results on the conditional stability of the scheme in the mesh energy norm and the discrete energy conservation law are derived from the corresponding general results in [24, 26] (proved by the energy method) by excluding the auxiliary sought functions. The case of general free terms in the main and auxiliary equations of the scheme is treated that is essential for the completeness of the stability analysis. The stability bound and the energy law contain the operators inverse to the Numerov averages in th space direction, . The imposed stability conditions are shown to be very close to those arising in the spectral analysis, and the sufficient conditions on their validity in the standard terms of the ratios of the time and spatial steps are given. Alternative forms of the stability bound and the energy conservation law not exploiting the operators are presented too. The 4th order error bound in the mesh energy norm is also rigorously proved for the first time based on the derived stability theorem.

Finally, the generalizations of the scheme to the cases of the non-uniform rectangular mesh in space and the non-uniform mesh in time are presented following [24], and the important case of the wave equation with variable coefficients is covered too. For general non-uniform meshes, the scheme has the 3th approximation order as other compact schemes. The non-uniform mesh in time is important to accomplish a dynamic choice of the time step. The scheme becomes non-local in time in this case but it still remains explicit in time and simply implemented. In the case of variable coefficients, the algorithm of the scheme implementation and its cost do not change, and the cost is similar to the ADI methods once again.

The paper is organized as follows. Section 2 is devoted to the construction of the scheme. In Section 3, the main theorems are proved. Section 4 deals with the generalizations of the scheme.

## 2 Construction of an explicit in time 4th-order compact vector scheme for the n-dimensional wave equation

We deal with the following initial-boundary value problem (IBVP) for the -dimensional wave equation under the nonhomogeneous Dirichlet boundary condition

 ∂2tu(x,t)−Lu(x,t)=f(x,t),  L:=a21∂21+…+a2n∂2n,  in  QT=Ω×(0,T); (2.1) u|ΓT=g(x,t);  u|t=0=u0(x),  ∂tu|t=0=u1(x),  x∈Ω=(0,X1)×…×(0,Xn). (2.2)

Here are constants (we take them different for uniformity with [24, 22]), , , . Also is the boundary of and is the lateral surface of .

We first rewrite formally the wave equation (2.1) as a system of equations, with a unique second order partial derivative in or in each of them:

 ∂2tu(x,t)−(a21u11(x,t)+…+a2nunn(x,t))=f(x,t)  in  QT, (2.3) ukk(x,t):=∂2ku(x,t),  1⩽k⩽n,  in  QT, (2.4)

where are the auxiliary sought functions. We consider smooth in solutions to the IBVP (2.1)-(2.2), then (2.3) implies

 ∂4tu=∂2t(Lu+f)=L∂2tu+∂2tf=L(a21u11+…+a2nunn+f)+∂2tf. (2.5)

Define the uniform mesh on a segment , with the step and . Let be the internal part of . Introduce the mesh averages (including the Numerov one) and difference operators in

 ¯sty=12(ˇy+y),  sty=12(y+^y),  stNy=112(ˇy+10y+^y), ¯δty=y−ˇyht, δty=^y−yht, ˚δty=^y−ˇy2ht, Λty=δt¯δty=^y−2y+ˇyh2t

with , and , as well as the operator of summation with the variable upper limit

 Imhty=htm∑l=1yl  for  1⩽m⩽M,  I0hty=0.

Due to Taylor’s formula in , the wave equation (2.1) and formula (2.5), we get

 Λtu=∂2tu+112h2t∂4tu+O(h4t)=a21u11+…+a2nunn+112h2tL(a21u11+…+a2nunn) +f+112h2t(L+Λt)f+O(h4t)  on  ωht, (2.6)

where we prefer to use instead of in [7].

Define the uniform mesh in with the step . Let be its internal part. Introduce the standard difference approximation to and the Numerov average in :

 (Λkw)l:=1h2k(wl+1−2wl+w\l−1),  skNwl:=112(wl−1+10wl+wl+1)=(I+112h2kΛk)wl (2.7)

on , where .

Define the uniform rectangular mesh in with . Let and be the internal part and boundary of . Define also the meshes in and on , where .

Let be the Euclidean space of functions defined on and equal 0 on and endowed with the inner product

 (v,w)h=h1…hn∑xk∈ωhv(xk)w(xk),  k=(k1,…,kn).

Any operator in generates the norm in . Recall that

 0<−Λk⩽λmaxkI  with  λmaxk=4h2ksinπ(Nk−1)2Nk<4h2k,  23I

Let be a difference operator such that on for sufficiently smooth in . Various can be used but below we confine ourselves by the simplest operator

 Lh:=a21Λ1+…+a2nΛn  %suchthat  0<−Lh<4(a21h21+…+a2nh2n)I  in  Hh. (2.9)

Then from expansion (2.6) we immediately get

 Λtu=(I+112h2tLh)(a21u11+…+a2nunn)+fh+O(|h|4) on ωh, (2.10)

where is the identity operator and

 fh:=f+112h2t(Lh+Λt)f=(stN+112h2tLh)f.

Using the Numerov approximation in for equation (2.4

) treated as an ordinary differential equation (ODE) in

, we have

 skNukk−Λku=O(h4k)  on  ωh,  1⩽k⩽n. (2.11)

Omitting the residual terms in the last two expansions (2.10)-(2.11), we pass to the three-level explicit in time vector compact scheme for the wave equation

 (2.12) skNvkk=Λkv  on  ωh,  1⩽k⩽n. (2.13)

Here the sought functions and are defined on and , respectively. We can set the discrete boundary conditions

 v|∂ωh=g,  a2kvkk|∂ωh=gk,  1⩽k⩽n, (2.14)

where in accordance with the wave equation in and the boundary condition one can calculate on and thus set

 gk:=∂2tg−∑1⩽l⩽n,l≠ka2l∂2lg−f for xk=0,Xk,  gk:=a2k∂2kg for xl=0,Xl, 1⩽l⩽n, l≠k,

on . Recall that equation (2.13) is the usual Numerov scheme for the ODE (2.4) in .

The constructed scheme is explicit in time. The reason is that is found explicitly on from (2.12), since means that , provided that the functions are known on . The scheme is not completely explicit since are (easily) computed from 1D difference equations (2.13) with the boundary conditions from (2.14) by solving the collections of corresponding linear algebraic systems with tridiagonal matrices. In this respect, the implementation is similar to the ADI methods. But the difference is that the functions are independent and can be computed not sequentially but in parallel when is known. Moreover, there is no need to store all of them simultaneously since only their weighted sum is used to compute , for .

We need also to define at the first time level with the 4th order of accuracy. This can be done explicitly by using Taylor’s formula and the wave equation, for example, see [2, 5, 7, 18]. But these formulas involve higher-order derivatives (or differences) of the initial functions and that is inconvenient in the non-smooth case. We prefer to avoid this and proceed alternatively following [24, 25] and construct the explicit equation for in the form close to (2.12)-(2.13):

 (δtv)0=12ht(I+112h2tLh)(a21v011+…+a2nv0nn)+u1h+12htf0h  on  ωh, (2.15) skNv0kk=Λkv0  on  ωh,  1⩽k⩽n, (2.16)

with suitable and . The values of on can be taken as in (2.14) (or according to the values of on ). See additional arguments in favor of these equations for after the proof of Theorem 3.2 below.

###### Lemma 2.1.

For the given functions

 u1h=(I+16h2tLh)u1  on  ωh, (2.17) f0h=f(0)dht+112h2tLhf0,  f(0)dht=f(0)d+O(h3t)  on  ωh, (2.18)

where with for any (see explicit forms of in Remark 2.2

below), the following estimate for the approximation error of equation (

2.15) holds

 ψ0:=(δtu)0−12ht(I+112h2tLh)(a21u110+…+a2nunn0)−u1h−12htf0h=O(|h|4) (2.19)

on .

###### Proof.

Using Taylor’s formula at , the initial conditions (2.2), the wave equation and formula (2.5), we get

 (δtu)0=u1+12ht(Lu0+f0)+16h2t[Lu1+(∂tf)0]+124h3t[L(Lu0+f0)+(∂2tf)0]+O(h4t).

Since , we derive

 ψ0=u1+16h2tLu1−u1h+12ht(f(0)d+112h2tLf0−f0h)+O(h4t).

For functions (2.17)-(2.18), this implies formula (2.19). ∎

###### Remark 2.1.

If is smooth and known analytically, then equation (2.15) can be simplified as

 (δtv)0=12ht(I+112h2tLh)Lu0+u1h+12htf0h  on  ωh, (2.20)

with omitting equations (2.16), and formula (2.19) remains valid. But, otherwise, equation (2.20) cannot be used and equations (2.15)-(2.16) become preferable and are more general.

###### Remark 2.2.

Let . If is sufficiently smooth in in (or ), then (see (2.18)) for the following three- and two-level approximations

 f(0)dht=712f0+12f1−112f2,  f(0)dht=13f0+23f1/2  with  f1/2:=f|t=ht/2 (2.21)

(or with ). See [24, Remark 3]. These formulas are checked easily applying Taylor’s formula at .

## 3 Stability and error bounds for the explicit 4th-order compact scheme

To study stability in more detail, we need to take the inhomogeneous version of equations (2.13) and (2.16):

 skNvmkk=Λkvm+bmk  on  ωh,  0⩽m⩽M−1,  1⩽k⩽n, (3.1)

where are given functions. Note that in practice these functions are never identically zero due to the round-off errors in computations, so their influence on the solution have to be studied. This is also necessary to derive an error bound below.

To state the stability theorem, we introduce the following self-adjoint operator

 (3.2)

We impose the following condition between the steps in time and space

 13h2t(a21h21+…+a2nh2n)⩽1−ε  for some  0⩽ε<1. (3.3)

The second inequality (2.8) implies in , and thus inequalities (2.9) and this condition ensure that

 εI

We also define the mesh norm .

###### Theorem 3.1.

Let and . Then, for scheme (2.12), (2.14)-(2.16) and (3.1), under the conditions (3.3) and

 14h2tAh⩽(1−ε20)I  for some  0<ε0<1, (3.5)

the following stability bound

 max1⩽m⩽M(ε20∥¯δtvm∥2h+∥¯stvm∥2Ah)1/2⩽(∥v0∥2Ah+ε−20∥u1h∥2h)1/2+2ε−10∥∥fh+βh∥∥L1ht(Hh) (3.6)

and the corresponding discrete energy conservation law

 ∥¯δtvm∥2h−14h2t∥¯δtvm∥2Ah+∥¯stvm∥2Ah =∥¯δtv1∥2h−14h2t∥¯δtv1∥2h+∥¯stv1∥2Ah+2Im−1ht(fh+βh,˚δtv)h =(Ahv0,stv0)h+(u1,δtv0)h+12ht(f0h+β0h,δtv0)h+2Im−1ht(fh+βh,˚δtv)h, (3.7)

for , are valid. Here

 βh:=(I+112h2tLh)(a21s−11Nb1+…+a2ns−1nNbn)  in  Hh  on  {tm}M−1m=0, (3.8)

and functions : and are any ( and are not only those specific defined above). The first equality (3.7) holds for any (not only defined by equations (2.15)-(2.16)).

###### Proof.

1. First, we recall the general three-level scheme

 BhΛtv+Ahv=φ  in  Hh  % on  ωht, (3.9) Bh(δtv)0+12htAhv0=u1+12htφ0  in  Hh (3.10)

(with the weight ), in particular, see [25, 26, 22, 24]. Here and are any operators in related by the inequality

 14h2tAh⩽(1−ε20)Bh  for some  0<ε0<1. (3.11)

Then the following stability bound

 max1⩽m⩽M(ε20∥¯δtvm∥2Bh+∥¯stvm∥2Ah)1/2 ⩽(∥v0∥2Ah+ε−20∥B−1/2hu1h∥2h)1/2+2ε−10∥∥B−1/2hf∥∥L1ht(Hh) (3.12)

and the discrete energy conservation law (that entails the stability bound)

 ∥¯δtvm∥2Bh−14h2t∥¯δtvm∥2Ah+∥¯stvm∥2Ah =∥¯δtv1∥2Bh−14h2t∥¯δtv1∥2h+∥¯stv1∥2Ah+2Im−1ht(φ,˚δtv)h =(Ahv0,stv0)h+(u1,δtv0)h+12ht(φ0,δtv0)h+2Im−1ht(φ,˚δtv)h,  1⩽m⩽M, (3.13)

are valid according to [24, Theorem 1] and the proof of Theorem 1 in [26] (see also [17]). Notice that the first and second equalities (3.13) are valid respectively for any and defined by equation (3.10). Also inequality (3.11) ensures that

 ε20∥¯δtvm∥2Bh⩽∥¯δtvm∥2Bh−14h2t∥¯δtvm∥2Ah  ∀vm−1,vm∈Hh.

2. Since and , clearly : and : , . From equation (3.1), we can express through :

 vkk=s−1kN(Λkv+bk)  in  Hh  on  {tm}M−1m=0,  1⩽k⩽n.

Inserting these formulas to equations (2.12) and (2.15), we get the closed equations for :

 Λtv=(I+112h2tLh)(a21s−11NΛ1+…+a2ns−1nNΛn)v+fh+βh  in  Hh  on  ωht, (3.14) (δtv)0=12ht(I+112h2tLh)(a21s−11NΛ1+…+a2ns−1nNΛn)v0+u1h+12ht(f0h+β0h)  in  Hh, (3.15)

with defined in (3.8).

These equations present a particular case of equations (3.9)-(3.10) with the operators and given in (3.2). Then the stability condition (3.11) takes the form (3.5), and the stability bound (3.6) and energy conservation law (3.7) follow from respective general relations (3.12) and (3.13). In addition, the following inequality holds

 ε20∥¯δtvm∥2h⩽∥¯δtvm∥2h−14h2t∥¯δtvm∥2Ah

that is essential on the left in (3.7). ∎

Note that the norms of and in (3.6) are bounded uniformly in by the same norms of , and due to the stability condition (3.3).

Let us compare the above stability conditions with the corresponding conditions arising in the frequently used spectral method (it was used in [7]

as well). In this method, one can consider the system of eigenvectors of the operator

in :

 −Lheℓ=λℓ(−Lh)eℓ on ωh,  eℓ=sinπl1x1X1…sinπlnxnXn,  1⩽l1⩽N1−1,…,1⩽ln⩽Nn−1.

with . Hereafter

are the corresponding eigenvalues of the operator

. Inserting the solutions in the form ,