The new alternating direction implicit difference methods for the wave equations

Journal of Computational and Applied Mathematics 230 (2009) 213–223

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

The new alternating direction implicit difference methods for the wave equations Jinggang Qin Institute of Plasma Physics, Chinese Academy of Sciences, Hefei, Anhui Prov., 230031, China

article

abstract

info

Article history: Received 28 May 2008 Received in revised form 18 October 2008 Keywords: Wave equation ADI method Compact difference scheme Error analysis

A new second-order alternating direction implicit (ADI) scheme, based on the idea of the operator splitting, is presented for solving two-dimensional wave equations. The scheme is also extended to a high-order compact difference scheme. Both of them have the advantages of unconditional stability, less impact of the perturbing terms on the accuracy, and being convenient to compute the boundary values of the intermediates. Besides this, the compact scheme has high-order accuracy and costs less in computational time. Numerical examples are presented and the results are very satisfactory. © 2008 Elsevier B.V. All rights reserved.

1. Introduction In this paper, we consider the following two-dimensional nonhomogeneous wave differential equation with initial and boundary conditions on domain Ω = [0, 1]2

 2  ∂ 2u ∂ u ∂ 2u − a + b + cu = g (x, y, t ), ∂t2 ∂ x2 ∂ y2 u(x, y, 0) = u0 (x, y), u(x, y, t ) = ϕ(x, y, t ),

ut (x, y, 0) = ψ(x, y)

(x, y) ∈ Ω , t ∈ (0, T ], (x, y) ∈ Ω ,

(x, y) ∈ ∂ Ω , t ∈ (0, T ],

(1.1) (1.2) (1.3)

where a and b are positive constants, c is a non-negative constant, g , u0 , ψ , and ϕ are given functions with sufficient smoothness, and ∂ Ω is the boundary of the domain Ω . Because of the importance of differential equations, the research on their numerical algorithms is always an active area in numerical computation. Since it is simple to use, the oldest finite difference method (FDM) is still used extensively in practical computations. Today, new difference schemes are constantly presented, and for multidimensional problems, alternating direction implicit schemes get much attention for their unconditional stability and high efficiency. This paper is concerned with numerical solutions to nonhomogeneous wave problems. The wave equation is often solved by explicit time-stepping schemes, which require to choose a time step size which is sufficiently small to satisfy the stability condition and to reduce numerical dispersion as well. As is well known, the alternating direction implicit (ADI) schemes [1–11] are unconditionally stable and only need to solve a sequence of tridiagonal linear systems. The ADI methods first introduced by Peaceman, Rachford [1] and Douglas [2] cost less CPU time to solve the heat equation. In the context of high order finite difference methods, compact schemes feature high-order accuracy and small stencils. Recently, there has been a renewed interest in the development and the application of compact finite difference methods for the numerical solution of the differential equations [4–11]. To obtain satisfactory higher order numerical results with reasonable computational cost, there have been attempts to develop higher order compact ADI methods [5,8,11].

E-mail address: [email protected]. 0377-0427/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2008.11.001

214

J. Qin / Journal of Computational and Applied Mathematics 230 (2009) 213–223

In this paper, we propose a new set of ADI methods for nonhomogeneous wave problems. In Section 2, we obtain a new scheme with truncation error O(1t 2 + h2 ) for two-dimensional wave equations and analyze the error estimate by discrete energy method. In Section 3, we further extend the method to compact ADI scheme with truncation error O(1t 2 + h4 ). In Section 4, we provide two numerical examples to illustrate the effectiveness of the scheme. We also compare the two methods with some other schemes and the results show they have much higher efficiency. 2. The new second-order ADI scheme for two-dimensional problems In this section, we derive a new second-order ADI method for the numerical solution of the wave equation (1.1)–(1.3). 2.1. Construction of the new ADI scheme Introducing v = ∂∂ut , we can rewrite (1.1) as follows

 2  ∂ u ∂v ∂ 2u − a 2 + b 2 + cu = g (x, y, t ), ∂t ∂x ∂y

(x, y) ∈ Ω , t ∈ (0, T ],

(2.1)

∂u , (x, y) ∈ Ω , t ∈ (0, T ], ∂t u(x, y, 0) = u0 (x, y), v(x, y, 0) = ψ(x, y) (x, y) ∈ Ω , v=

u(x, y, t ) = ϕ(x, y, t ),

(2.2) (2.3)

(x, y) ∈ ∂ Ω , t ∈ (0, T ].

(2.4)

The domain Ω is divided into a mesh by points xi = ih, yj = jh(i, j = 0, . . . , N ), denoted by Ωh , where h = mesh size in both x and y directions. Let t = n1t , t n

n+ 12

n+ 12

1 N

is the spatial

= (n + )1t, with 1t being the time increment. Denote the values 1 2

n+ 21

of u(xi , yj , t ) and g (xi , yj , t ) by and gi,j respectively. Also denote the difference solutions to uni,j and vin,j by Uin,j and n Vi,j , respectively. Applying Crank–Nicolson implicit discretization to (2.1), we have n

Vin,j+1 − Vin,j

− (aδx2 + bδy2 )

4t Vin,j+1 + Vin,j

=

2

Uin,j+1 − Uin,j

1t

uni,j

Uin,j+1 + Uin,j

+c

2

Uin,j+1 + Uin,j 2

1t 2 4

= Vin,j +

(2.6)

Uin,j+1 = Let A = 1 +

4

1t 2

c 1t 2 , 4

Vin,j+1 −

(aδx2 + bδy2 )Vin,j+1 +

1t 2

(2.5)

,

where δx2 Uin,j = h12 (Uin+1,j − 2Uin,j + Uin−1,j ) and δy2 Uin,j = From (2.5) and (2.6), we obtain Vin,j+1 −

n+ 1

= gi,j 2 ,

(aδx2 + bδy2 )Vin,j −

c 1t 2 4

1 h2

(Uin,j+1 − 2Uin,j + Uin,j−1 ).

Vin,j+1

c 1t 2 4

n+ 21

Vin,j + (aδx2 + bδy2 − c 1t )Uin,j + 1tgi,j

,

(2.7)

(Vin,j+1 + Vin,j ) + Uin,j .

(2.8)

Eq. (2.7) can be reorganized as

1t 2 4A

(aδx2 + bδy2 )Vin,j+1 = Vin,j +

1t 2 4A

(aδx2 + bδy2 )Vin,j −

c 1t 2 2A

Vin,j +

1t A

(aδx2 + bδy2 − c )Uin,j +

1t A

n+ 12

gi,j

.

(2.9)

1t 4 ab x2 y2 16A2

δ δ (Vin,j+1 − Vin,j ) to the left hand side of (2.9), we factor (2.9) as    1t 2 2 1t 2 2 1− aδ x 1− bδy Vin,j+1

Adding perturbing term

4A

 = 1+

4A

1t

2

4A

aδx2

 1+

1t 2 4A

bδy2



Vin,j −

c 1t 2 2A

Vin,j +

1t A

(aδx2 + bδy2 − c )Uin,j +

1t A

n+ 12

gi,j

.

(2.10)

J. Qin / Journal of Computational and Applied Mathematics 230 (2009) 213–223

215

n+ 12

, we obtain the PR ADI-like scheme    1t 2 2 1t 2 2 c 1t 2 n 1t 4A 1t n+ 12 n+ 1 2 1− aδx Vi,j 2 = 1 + bδy Vin,j − V i ,j + + b δ − c Uin,j + g , y 2 4A 4A 4A 2A 1t 2A i,j       1t 2 2 1t 2 2 c 1t 2 n 1t 4A 1t n+ 12 n+ 21 n+1 2 bδy Vi,j = 1 + aδx Vi,j − Vi,j − − bδy + c Uin,j + g 1− , 4A 4A 4A 2A 1t 2 2A i,j

Introducing the intermediate variable Vi,j





Uin,j+1 =

1t 2



(2.11)

(2.12)

(Vin,j+1 + Vin,j ) + Uin,j .

(2.13)

From (2.11) and (2.12), we get the boundary equation n+ 12

=

V i ,j

Uin,j+1 + Uin,j

1t

−

1t 2

bδy2 (Vin,j+1 − Vin,j ).

8A

(2.14)

But we usually compute the boundary values of the intermediate variable using the following simple equality n+ 12

Uin,j+1 + Uin,j

, 1t where i = 0, N , j = 0, . . . , N − 1. =

V i ,j

(2.15)

We see that schemes (2.11)–(2.15) can solve nonhomogeneous wave problems. In addition, following the idea of Douglas [2,3], a Douglas-like scheme follows:

 1−

 1−

a1t 2 4A b1t 2 4A

δ

2 x

δy2



n+ 21

Vi,j



−

Vin,j

1t 2

 =

n+ 21

The intermediate value of Vi,j




2A n+ 21

Vin,j+1 − Vin,j = Vi,j




aδx2 + bδy2 Vin,j −

c 1t 2 2A

Vin,j +

1t A

(aδx2 + bδy2 − c )Uin,j +

1t

− Vin,j .

A

n+ 12

gi,j

,

(2.16)

(2.17)

on the boundary is easily obtained from (2.17).

2.2. Error estimate We have derived a new kind of ADI scheme in Section 2.1. We further analyze the error of the scheme starting from (2.10). Expanding (2.10), we get Vin,j+1 − Vin,j

1t

− aδx2 + bδy2


Uin,j+1 + Uin,j 2

+c

Uin,j+1 + Uin,j 2

+

1t 3 16A

n+ 12

abδx2 δy2 (Vin,j+1 − Vin,j ) = gi,j

.

(2.18)

Discretizing (2.1) and (2.2) as above and letting ξ = u − U , η = v − V , we obtain the error equations 1 ηin,+ − ηin,j j

1t 1 ηin,+ + ηin,j j

2

− aδx2 + bδy2


ξin,j+1 + ξin,j

ξin,j+1 − ξin,j

=

1t

2

+c

ξin,j+1 + ξin,j 2

+

1t 3

n+ 21

1 − ηin,j ) = Ri,j abδx2 δy2 (ηin,+ j 16A

,

(2.19)

n+ 1

+e Ri,j 2 ,

(2.20)

where the truncation errors n+ 12

Ri,j

= 1t 2



1 ∂ 3v

−

1



∂ 4u ∂ 4u + b 2 2 ∂ t 2 ∂ y2 ∂ t ∂ x

a

 +

c ∂ 2u

+

∂ 5v 16A ∂ t ∂ x2 ∂ y2

1t 4

24 ∂ 8 8∂  ∂ 4u ∂ 4u n+ 12 4 − a 4 + b 4 (xi , yj , t ) + O(1t + 1t 2 h2 + h4 ), 12 ∂x ∂y  2  1∂ v 1 ∂ 3u 1 n+ 12 2 e Ri,j = 1t − (xi , yj , t n+ 2 ) + O(1t 4 ) 8 ∂t2 24 ∂ t 3 1t 2 ∂ 3 u 1 = (xi , yj , t n+ 2 ) + O(1t 4 ). 3 12 ∂ t Hence there exist positive constants C1 , C2 such that t3

h2

n+ 1

|Ri,j 2 | ≤ C1 (1t 2 + h2 ),

1

n+ |e Ri,j 2 | ≤ C2 1t 2 .

t2



1

(xi , yj , t n+ 2 ) (2.21)

(2.22)

216

J. Qin / Journal of Computational and Applied Mathematics 230 (2009) 213–223

Assume that Vh = {w | w = {wi,j } ∈ Ωh , and w|∂ Ω = 0}. For w = {wi,j } ∈ Vh and v = {vi,j } ∈ Vh , denote by 1 vin,+ − vin,j j

δt vin,j =

1t

,

δx vin,j =

vin+1,j − vin,j h

,

δy vin,j =

vin,j+1 − vin,j h

.

For ∀v, w ∈ Vh , define their discrete inner product and norms as follows

(v, w) =

N −1 X

vi,j wi,j h2 ,

kvk =

p (v, v),

i,j=1

v u N −1 N −1 uX X (δx vi,j )2 h2 , kδx vk = t

v u N −1 N −1 uX X kδy vk = t (δy vi,j )2 h2 ,

i=0 j=1

i=1 j=0

v u N −1 N −1 uX X kδx δy vk = t (δx δy vi,j )2 h2 ,

|v|H 1 =

q

kδx vk2 + kδy vk2 .

i=0 j=0

Lemma 1 ([12,13]). For w ∈ Vh , the following equalities hold

−

N −1 X N −1 X

wi,j (δx2 wi,j )h2 =

i=0 j=1

i=1 j=1

−

N −1 X N −1 X

N −1 X N −1 X (δx wi,j )2 h2 = kδx wk2 ,

N −1 X N −1 X wi,j (δy2 wi,j )h2 = (δy wi,j )2 h2 = kδy wk2 .

i=1 j=1

i=1 j=0

From Lemma 1 and Cauchy–Schwarz inequality, multiplying easily follows that

 I1 =

ηn+1 − ηn ηn+1 + ηn , 1t 2

1 ηin,+ +ηin,j 2 j

2

h to both sides of (2.17), computing the inner product, it



1

=

(kηn+1 k2 − kηn k2 ),   ξ n +1 + ξ n η n +1 + η n −a δx2 , 2 2   n +1 + ξ n ξ n+1 − ξ n en+ 1 2ξ −a δx , +R 2 2 1t   n +1  a ξ + ξ n en+ 1 n +1 2 (kδx ξ k − kδx ξ n k2 ) − −a δx , δx R 2 21t

n+1 2n

ξ a +ξ n+ 1

kδxe (kδx ξ n+1 k2 − kδx ξ n k2 ) − a

δx

R 2k 21t 2 a a 1 (kδx ξ n+1 k2 − kδx ξ n k2 ) − (kδx ξ n+1 k2 + kδx ξ n k2 + kδxe Rn+ 2 k2 ), 21t 4

21t

I21 =

= = ≥ ≥

b

b

1

(kδy ξ n+1 k2 − kδy ξ n k2 ) − (kδy ξ n+1 k2 + kδy ξ n k2 + kδye Rn+ 2 k2 ), 4  n +1   n +1  ξ + ξ n η n +1 + η n ξ + ξ n ξ n+1 − ξ n en+ 1 I3 = c , =c , +R 2 , 2 2 2 1t

n+1 n

ξ c + ξ 1 n +1 2 n 2 n +

ke ≥ (kξ k − kξ k ) − c

R 2k 21t 2 c c 1 ≥ (kξ n+1 k2 − kξ n k2 ) − (kξ n+1 k2 + kξ n k2 + ke Rn+ 2 k2 ), 21t 4   3 1t ηn+1 + ηn ab1t 3 2 2 n +1 n I4 = abδx δy (η − η ), = (kδx δy ηn+1 k2 − kδx δy ηn k2 ), I22 ≥

21t

16A

 I5 =

1

Rn+ 2 ,

2

η n +1 + η n 2

 ≤

1 4

1

32A

1

kRn+ 2 k2 + (kηn+1 k2 + kηn k2 ). 4

(2.23)

(2.24) (2.25)

(2.26) (2.27)

(2.28)

J. Qin / Journal of Computational and Applied Mathematics 230 (2009) 213–223

217

Multiplying 21t to both sides of the result derived from (2.23)–(2.28), we have

(kηn+1 k2 − kηn k2 ) + a(kδx ξ n+1 k2 − kδx ξ n k2 ) + b(kδy ξ n+1 k2 − kδy ξ n k2 ) + c (kξ n+1 k2 − kξ n k2 ) + ≤

ab1t 4

1t 2

+

16A

(kδx δy ηn+1 k2 − kδx δy ηn k2 ) 1

1

1

1

(akδxe Rn+ 2 k2 + bkδye Rn+ 2 k2 + c ke Rn+ 2 k2 + kRn+ 2 k2 )

1t 2

(kηn+1 k2 + kηn k2 + a(kδx ξ n+1 k2 + kδx ξ n k2 ) + b(kδy ξ n+1 k2 + kδy ξ n k2 )

+ c (kξ n+1 k2 + kξ n k2 )).

(2.29)

Summing for n to both sides of (2.29), we obtain

kηn k2 + akδx ξ n k2 + bkδy ξ n k2 + c kξ n k2 ≤

n−1   1t X 1 1 1 1 Rl+ 2 k2 + bkδye Rl+ 2 k2 + c ke kRl+ 2 k2 + akδxe R l + 2 k2

2

+

l=1 n−1 1t X

2


kηl k2 + akδx ξ l k2 + bkδy ξ l k2 + c kξ l k2 ,

(2.30)

l=1

noticing that ξ 0 = η0 = 0, by Gronwall Lemma, we obtain the following theorem. Theorem 1. Assume that u, v are the accurate solutions of (2.1)–(2.4) with sufficient smoothness and U , V are the difference solutions of (2.11)–(2.13). Let ξ = u − U , η = v − V , then there exists a positive constant K independent of 1t and h, such that max (kηn k + akδx ξ n k + bkδy ξ n k + c kξ n k) ≤ K (1t 2 + h2 ).

n≤[T /1t ]

3. Generalized compact ADI scheme In this section, we extend the method to compact ADI scheme for the numerical solution of the wave Eqs. (1.1)–(1.3). 3.1. Construction of the compact ADI scheme Introducing fx =

∂ 2u , ∂ x2

fy =

∂ 2u . ∂ y2

(3.1)

We can rewrite (2.1) as follows:

∂v − (afx + bfy ) + cu = g (x, y, t ), ∂t ∂ 2u fx = , ∂ x2 ∂ 2u . fy = ∂ y2

(3.2) (3.3) (3.4)

Using the fourth-order compact finite difference discretization [6,8], we discretize (3.3), (3.4) as follows Ax (fx )ni,j := Ay (fy )ni,j :=

1 12 1 12

5

1

6 5

12 1

6

12

(fx )ni−1,j + (fx )ni,j + (fy )ni,j−1 + (fy )ni,j +

(fx )ni+1,j = δx2 Uin,j ,

(3.5)

(fy )ni,j+1 = δy2 Uin,j ,

(3.6)

where δx2 Uin,j = h12 (Uin+1,j − 2Uin,j + Uin−1,j ), and δy2 Uin,j = h12 (Uin,j+1 − 2Uin,j + Uin,j−1 ). The Crank–Nicolson implicit discretization of (3.2) reads Vin,j+1 − Vin,j

4t

−

1 a((fx )ni,+ + (fx )ni,j ) + b((fy )ni,+j 1 + (fy )ni,j ) j

2

+c

Uin,j+1 + Uin,j 2

n+ 1

= gi,j 2 ,

(3.7)

218

J. Qin / Journal of Computational and Applied Mathematics 230 (2009) 213–223

where i, j = 1, 2, . . . , N − 1. Multiplying (3.7) by 1t, applying to both of its sides with Ax Ay , and using the fact that the two operators Ax and Ay commute with each other, we obtain the implicit approximation of (2.1)–(2.4)

1t

Ax Ay (Vin,j+1 − Vin,j ) −

+

c 1t




n+ 21

,

Ax Ay (Uin,j+1 + Uin,j ) = 1tAx Ay gi,j

2

Vin,j+1 + Vin,j c 1t 2 , 4

 Ax Ay −

Uin,j+1 − Uin,j

=

2 Let A = 1 +

1 a(Ay Ax (fx )ni,+ + Ay Ax (fx )ni,j ) + b(Ax Ay (fy )ni,+j 1 + Ax Ay (fy )ni,j ) j

2

1t

(3.8)

.

(3.9)

(3.8) can be conveniently written as

1t 2 4A

δ −

aAy x2

1t 2 4A

δ

bAx y2



Vin,j+1

 =

Ax Ay +

+ ab1t 4 16A2

Adding perturbing term n+ 12

Vi,j

1t A

1t 2 4A

δ +

aAy x2

1t 2 4A

δ

bAx y2



Vin,j −

(aAy δx2 + bAx δy2 − cAx Ay )Uin,j +

1t A

c 1t 2 2A

Ax Ay Vin,j n+ 12

Ax Ay gi,j

.

(3.10)

δx2 δy2 (Vin,j+1 − Vin,j ) to the left hand side of (3.10), and introducing the intermediate variable

, we obtain the compact PR ADI-like scheme a1 t 2

   b1t 2 2 1t 4A c 1t 2 1t n+ 1 n+ 1 δx2 Vi,j 2 = Ay + δy Vin,j − Ay Vin,j + ( 2 Ay + bδy2 − cAy )Uin,j + Ay gi,j 2 , 4A 4A 4A 2A 1t 2A     2 2 2 1 a1t 2 b1t 2 c 1t 1t 4A 1t n+ n+ 1 Ay − δy Vin,j+1 = Ax + δx Vi,j 2 − Ay Vin,j − ( 2 Ay − bδy2 + cAy )Uin,j + Ay gi,j 2 , 4A 4A 4A 2A 1t 2A 

Ax −

Uin,j+1 =

1t 2

(Vin,j+1 + Vin,j ) + Uin,j .

(3.11)

(3.12) (3.13)

From (3.11) and (3.12), we obtain the boundary equation n+ 12

= Ay

Ax Vi,j

Uin,j+1 + Uin,j

−

1t

b1t 2 4A

δy2 (Vin,j+1 − Vin,j ).

But usually we compute the boundary values of the intermediate using the following simple equality n+ 12

Vi,j

=

Uin,j+1 + Uin,j

1t

,

(3.14)

where i = 0, N , j = 0, . . . , N − 1. In addition, following the idea of Douglas [2,3], we can get a compact Douglas scheme

 Ax −

a1t 2 4A

δ

2 x



n+ 12

Vi,j

−

Vin,j

 =

1t 2

Ay −

b1t 2 4A

δy2



n+ 21

Vin,j+1 − Vin,j = Vi,j




n+ 21

The intermediate value of Vi,j




2A

+ 

aAy δx2 + bAx δy2 Vin,j −

1t A

c 1t 2 2A

Ax Ay Vin,j

(aAy δx2 + bAx δy2 − cAx Ay )Uin,j +

1t A

n+ 12

Ax Ay gi,j

,

(3.15)

− Vin,j .

(3.16)

on the boundary is easily obtained from (3.16).

3.2. Error estimate We have derived a kind of compact PR ADI-like scheme in Section 3.1. In the following, we further analyze the error of the scheme. Expanding (3.11)–(3.12), we get

Ax Ay

Vin,j+1 − Vin,j

1t

− aAy δx2 + bAx δy2


Uin,j+1 + Uin,j 2

+ cAx Ay

Uin,j+1 + Uin,j 2

+

1t 4 16A

abδx2 δy2

Vin,j+1 − Vin,j

1t

n+ 1

= Ax Ay gi,j 2 .

(3.17)

J. Qin / Journal of Computational and Applied Mathematics 230 (2009) 213–223 h2

Noting that Ax Uin,j = (1 + h2

 1+

δ

2 x

δ )Uin,j , Ay Uin,j = (1 +

2 12 x

 1+

h2

δ

2 y

Vin,j+1 − Vin,j



h2

δ )Uin,j , (3.17) can be written as

2 12 y

  − a 1+

h2

δ

2 y





δ +b 1+ 2 x

h2

1t 12 12  n + 1 n + 1 n n Vi,j − Vi,j h2 h2 2 Ui,j + Ui,j 1t 4 + c 1 + δx2 1+ δy + abδx2 δy2 12 12 2 16A 1t    2 2 1 h 2 h n+ 1+ δy gi,j 2 . = 1 + δx2 12



12

219

δ

2 x



δ

2 y



Uin,j+1 + Uin,j 2



12

(3.18)

12

Discretizing (2.1), (2.2) as above and letting ξ = u − U , η = v − V , we obtain the error equations h2

 1+

12

δx2



2

 +c 1 + 1 + ηin,j ηin,+ j

2

1+ h

12

=

δx2

h2 12

δy2



 1+

ξin,j+1 − ξin,j 1t

1 ηin,+ − ηin,j j

h

2

12

      n +1 ξi,j + ξin,j h2 h2 − a 1 + δy2 δx2 + b 1 + δx2 δy2

1t  n +1 ξi,j + ξin,j 2

δy

2

12

+

1t

4

16A

abδx2 δy2

12

η

n +1 i,j

−η

n i ,j

1t

2

n+ 1

= Ri,j 2 ,

n+ 21

+e Ri,j

(3.19)

(3.20)

where the truncation errors n+ 12

Ri,j

1

n+ e R 2 i,j

  ∂ 2v 1 ∂ 4u ∂ 4u c ∂ 2u 1t 4 ∂ 5v = 1t Ax Ay 3 − aAy 2 2 + bAx 2 2 + + ab 24 ∂t 8 ∂t ∂x ∂t ∂y 8 ∂t2 16A2 ∂ t ∂ x2 ∂ y2   h4 ∂ 6u ∂ 6u 1 − aAy 6 + bAx 6 (xi , yj , t n+ 2 ) + O(1t 4 + 1t 2 h4 + h6 ), 240 ∂x ∂y   2 1 ∂ 3u 1∂ v 1 − (xi , yj , t n+ 2 ) + O(1t 4 ) = 1t 2 8 ∂t2 24 ∂ t 3 1t 2 ∂ 3 u 1 = (xi , yj , t n+ 2 ) + O(1t 4 ). 12 ∂ t 3 2



1

(3.21)

(3.22)

Hence, there exist positive constants C 0 , C 00 such that 1

n+ 1

n+ |e Ri,j 2 | ≤ C 00 1t 2 .

|Ri,j 2 | ≤ C 0 (1t 2 + h4 ),

The definitions of discrete inner product and norm are the same as those in Section 2.2. Lemma 2 ([12]). For w ∈ Vh , the following inequalities hold 1 4 1 4

1

h2 kδx wk2 ≤ kwk2 ,

4

h2 kδy wk2 ≤ kwk2 , 1

h2 kδx δy wk2 ≤ kδx wk2 ,

4

h2 kδx δy wk2 ≤ kδy wk2 .

From Lemmas 1 and 2 and Cauchy–Schwarz inequality, multiplying product, we have

1 ηin,+ +ηin,j 2 j

2

h to both sides of (3.19) and computing the inner

  η n +1 + η n δy2 δt ηn , 12 2   12 2       n +1 n n +1 n +1 n +1 η + η h η + ηn h2 + ηn h4 + ηn n 2 n 2 n η 2 2 n η = δt η , + δx δt η , + δy δt η , + δx δy δt η , 2 12 2 12 2 144 2  2 2 1 h h kηn+1 k2 − kηn k2 − (kδx ηn+1 k2 − kδx ηn k2 ) − (kδy ηn+1 k2 − kδy ηn k2 ) = 21t 241t 241t  h4 + (kδx δy ηn+1 k2 − kδx δy ηn k2 ) , (3.23) 2881t 

I1 =

1+

h2

δx2



1+

h2

220

J. Qin / Journal of Computational and Applied Mathematics 230 (2009) 213–223

 ξ n+1 + ξ n ηn+1 + ηn −a 1+ δ δ , 12 2 2     n +1 n n +1 n 2 ξ + e η + η ah ξ n+1 + ξ n ηn+1 + ηn −a δx2 , − δx2 δy2 , 2 12 2   n +1  n+1 2 n  2 ξ n +1 − ξ n ξ + ξ n en+ 1 ξ +ξ , δx + a δx , δx R 2 a δx 2 1t 2     ξ n +1 + ξ n ξ n +1 − ξ n ah2 ξ n+1 + ξ n ah2 1 δx δy , δx δy − δx δy , δx δye Rn+ 2 − 12 2 1t 12 2
ah2 1
a n +1 2 n 2 n+1 2 kδx ξ k − kδx ξ k − kδx δy ξ k − kδx δy ξ n k2 21t 12 21t a ah2 1 1 n +1 2 n 2 − (kδx ξ k + kδx ξ k + kδxe Rn+ 2 k2 ) − Rn+ 2 k2 (kδx δy ξ n+1 k2 + kδx δy ξ n k2 ) + kδx δye h2



I21 =

= =

≥

4 a

≥

kδx ξ

21t a

−

3

2 x

k − kδx ξ k

n +1 2

n 2




48

ah2

−

1

kδx δy ξ n+1 k2 − kδx δy ξ n k2

12 21t 1 Rn+ 2 k2 ), (kδx ξ n+1 k2 + kδx ξ n k2 + kδxe

 I22 = −b

1+

b

(kδy e

≥

2 y



21t b

h2 12


(3.24)

  ξ n+1 + ξ n ηn+1 + ηn δx2 δy2 , 2

2

k − kδy e k ) −

n +1 2

n 2

bh2

1

12 21t

(kδx δy en+1 k2 − kδx δy en k2 )

1 Rn+ 2 k2 ), − (kδy ξ n+1 k2 + kδy ξ n k2 + kδye

(3.25)

3

h2

 I3 = c

1+

12

 1+

= c ≥

h

c

2

12

δx2 δx2

h2

 1+

12

 1+

h

2

12

δy2 δy2

(kξ n+1 k2 − kξ n k2 ) −

21t ch2

1

−



ξin,j+1 + ξin,j ηn+1 + ηn ,



ξ

2

ch2

n +1 i ,j

2

+ξ

n i,j

2 1

12 21t

!

,

ξ

+ ξ n en+ 1 +R 2 1t

n +1

!

(kδx ξ n+1 k2 − kδx ξ n k2 ) ch4

(kδy ξ n+1 k2 − kδy ξ n k2 ) +

1

(kδx δy ξ n+1 k2 − kδx δy ξ n k2 )

12 21t 144 21t c ch4 1 n +1 2 n 2 n+ 21 2 e − (kξ k + kξ k + kR k ) − (kδx δy ξ n+1 k2 + kδx δy ξ n k2 + ke Rn+ 2 k2 ), 3 576 I4 =

1t 3



16A2 ab1t 3

abδx2 δy2 (ηn+1 − ηn ),

ηn+1 + η 2

(kδx δy ηn+1 k2 − kδx δy ηn k2 ),   n+1 + ηn n+ 12 η I5 = R , =

(3.26)

 n

(3.27)

32A2

2

≤

1 4

n+ 12

(kR

k + kηn+1 k2 + kηn k2 ). 2

(3.28)

Multiplying 21t to both sides of the result derived from (3.23)–(3.28), we have h2

h2

h4

kηn+1 k2 − kδx ηn+1 k2 − kδy ηn+1 k2 + kδx δy ηn+1 k2 12 12 144     h2 h2 n +1 2 n+1 2 n+1 2 n +1 2 + a kδx ξ k − kδx δy ξ k + b kδy ξ k − kδx δy ξ k 12 12   2 2 4
h h h ab1t 4 n +1 2 n+1 2 n+1 2 n+1 2 + c kξ k − kδx ξ k − kδy ξ k + kδx δy ξ k + kδx δy ηn+1 k2 − kδx δy ηn k2 2 12

≤ kηn k2 −

h

2

12

kδx ηn k2 −

12

h

2

12

kδy ηn k2 +

144

h

4

144

kδx δy ηn k2

16A

J. Qin / Journal of Computational and Applied Mathematics 230 (2009) 213–223

h2





h2



+ a kδx ξ n k2 − kδx δy ξ n k2 + b kδy ξ n k2 − kδx δy ξ n k2 12 12   2 2 4 h h h + c kξ n k2 − kδx ξ n k2 − kδy ξ n k2 + kδx δy ξ n k2 12

2a1t

+

3 2c 1t

+

3

1t

+

2

12

221



144

n+ 21

R (kδx ξ n+1 k2 + kδx ξ n k2 + kδxe

k2 ) +

2b1t 3

1 Rn+ 2 k2 ) (kδy ξ n+1 k2 + kδy ξ n k2 + kδye

ch4 1t 1 1 (kξ n+1 k2 + kξ n k2 + ke Rn+ 2 k2 ) + (kδx δy ξ n+1 k2 + kδx δy ξ n k2 + kδx δye Rn+ 2 k2 ) 288

(kR

n+ 21

2

k + kη

k + kη k ).

n +1 2

n 2

(3.29)

Summing for n to both sides of (3.29), and by Lemma 2, we obtain 1 3

2

ch4

3

144

kηn k2 + (akδx ξ n k2 + bkδy ξ n k2 + c kξ n k2 ) + n −1 X

≤ 1t

kηn k2 +

l =1 n −1 1t X

+

2

kδx δy ξ n k2

n −1 41t X

3

1

kRn+ 2 k2 +

l =1

(akδx ξ n k2 + bkδy ξ n k2 + c kξ n k2 ) + 1t

n−1 X ch4

l =1

144 l=1

kδx δy ξ n k2

n−1 21t X

3

1 1 1 (akδxe Rn+ 2 k2 + bkδye Rn+ 2 k2 + c ke Rn+ 2 k2 ) + 1t

l=1

n −1 X ch4

144 l =1

1 kδx δye Rn+ 2 k2 .

(3.30)

Noticing that ξ 0 = η0 = 0, by Gronwall Lemma, we obtain the following theorem. Theorem 2. Assume that u, v are the accurate solutions of (2.1)–(2.4) with sufficient smoothness and U , V are the difference solutions of (3.11)–(3.13). Let ξ = u − U , η = v − V , then there exists a positive constant K 0 independent of 1t and h, such that

 max

n≤[T /1t ]

kηn k + akδx ξ n k + bkδy ξ n k + c kξ n k +

ch4 144

 kδx δy ξ n k ≤ K 0 (1t 2 + h4 ).

4. Numerical experiment In this section, we exemplify two numerical examples to illustrate the effectiveness of the present ADI scheme. The computer language used for the programming is Fortran, and the programs are performed on Lenovo M5400(CPU is AMD Athlon64 X2 Dual Core Processor 4400+, and Memory is 2.00 G). The new set of ADI schemes obtained in this article are compared to the schemes (2.16)–(2.17), (3.15)–(3.16) and the following ADI schemes [13]: the second-order ADI scheme:

1t 2

 1−

 1−

2

1t 2 2

aδ

2 x

bδ

2 y





U i,j = (aδx2 + bδy2 )Uik,j + cUik,j + gik,j ,

δt2 Uik,j = U i,j , 1



∂ 2 u0 ∂ 2 u0 + ∂ x2 ∂ y2



1

= u0 (xi , yj ) + 1t ψ(xi , yj ) + 1t (xi , yj ) + 1t 2 g (xi , yj , 0) 0 ≤ i, j ≤ N . 2 2   2 1t U i,j = 1 − bδy2 δt2 Uik,j i = 0, N , j = 1, . . . , N − 1, Ui1,j

2

2

δt2 Uik,j =

1 21t

(ϕ(xi , yj , t k+1 ) − ϕ(xi , yj , t k−1 )) i = 1, . . . , N , j = 0, N

the compact ADI scheme:

 Ax −

 Ay −

1t 2 2

1t 2 2

aδx2 bδy2





U i,j = (aAy δx2 + bAx δy2 )Uik,j + cAx Ay Uik,j + Ax Ay gik,j ,

δt2 Uik,j = U i,j ,

(4.1)

222

J. Qin / Journal of Computational and Applied Mathematics 230 (2009) 213–223

Table 1 The L2 -norm and CPU time for the second-order schemes. N = 40, 1t = h

N = 60, 1t = h

SCHEME

L2 -norm

CPU time (s)

L2 -norm

CPU time (s)

PR DOU SchemeI

4.0266 × 10−4 4.1016 × 10−4 1.6469 × 10−3

0.12500 0.12500 0.12500

1.7203 × 10−4 1.7335 × 10−4 7.3023 × 10−4

0.406250 0.453125 0.406250

Table 2 The L2 -norm and CPU time for the compact schemes. N = 40, 1t = h2

N = 60, 1t = h2

SCHEME

L2 -norm

CPU time (s)

L2 -norm

CPU time (s)

CPR CDOU CSchemeI

3.3075 × 10−7 3.3101 × 10−7 1.0870 × 10−6

5.343750 6.296875 5.812500

6.5916 × 10−8 6.5886 × 10−8 2.1494 × 10−7

26.65625 31.67188 28.95313

Table 3 The L2 -norm and CPU time for all schemes with h = 0.025, 1t = h2 . SCHEME

L2 -norm

PR CPR Dou CDou SchemeI CSchemeI

3.5947 × 10 1.2506 × 10−7 3.5947 × 10−5 1.2513 × 10−7 3.6201 × 10−5 4.4343 × 10−7

CPU time (s) −5

Ui1,j = u0 (xi , yj ) + 1t ψ(xi , yj ) +

 U i ,j =

Ay −

1t 2 2

bδy2



1 2

1t 2



∂ 2 u0 ∂ 2 u0 + 2 ∂x ∂ y2



2.078125 2.234375 2.031250 2.671875 2.013256 2.453125

1

(xi , yj ) + 1t 2 g (xi , yj , 0) 0 ≤ i, j ≤ N . 2

δt2 Uik,j i = 0, N , j = 1, . . . , N − 1,

1 (ϕ(xi , yj , t k+1 ) − ϕ(xi , yj , t k−1 )) i = 1, . . . , N , j = 0, N . (4.2) 21t For convenience, we denote the present second-order ADI scheme (2.11)–(2.13) by PR, the present compact ADI scheme (3.11)–(3.13) by CPR, the second-order ADI scheme (2.16)–(2.17) by Dou, the compact ADI scheme (3.15)–(3.16) by CDou, the second-order ADI scheme (4.1) by SchemeI and the compact ADI scheme (4.2) by CSchemeI.

δt2 Uik,j =

Example 1. Let a = 1.0, b = 1.0, c = 0, and g = 3π 2 e−π t sin π (x + y) in (1.1), then the true solution is u = e−π t sin π (x + y). We chose different grids and compared the accuracy of the computed solutions with other schemes. The solutions were obtained for t = 2.0. The quantities that we compared are the maximum L2 -norm and the total elapsed time (CPU) in seconds. Tables 1 and 2 display the L2 -norm errors and the CPU time. From Table 1, one can easily see that the error from the present second-order scheme is smaller than that from other schemes, and it converges as fast as other schemes. From Table 2, we can find the superiority of the present compact ADI scheme, which has the highest accuracy and fastest computational efficiency. √

Example 2. Let a = 1.0, b = 1.0, c = 0, and g = 0 in (1.1). The exact solution is u = ex+y− 2t , the initial and boundary conditions can be obtained based on the exact solution. We chose the grids 40 × 40, and the time step with 1t = h2 . The solutions were obtained for t = 1.0. We compared the accuracy of the computed solutions from the present ADI scheme and other four ADI schemes. The quantities that we compared are the L2 -norm and the total elapsed time (CPU) in seconds. Table 3 lists the total elapsed time (CPU) in seconds and the L2 -norm error. One can see that the error and the elapsed time from the same order scheme is very close. From the table, we can see that the present fourth-order scheme not only has high accuracy but also converges much faster than other compact schemes. So, the present compact scheme is more efficient than other methods. 5. Conclusions In this paper, we proposed two new ADI schemes for solving nonhomogeneous wave problems. It has been proved that they are unconditionally stable. The compact scheme proposed in this paper is fourth-order accurate in space and

J. Qin / Journal of Computational and Applied Mathematics 230 (2009) 213–223

223

second-order accurate in time, and allows a considerable saving in computing time. Numerical examples are given to test their high accuracy and to show their superiority over some other schemes in terms of accuracy and computational costs. Acknowledgements The author is very grateful to Professor Wang Tongke of Tianjin Normal University for his encouragement and to referees for their valuable remarks, which helped to improve the presentation substantially. References [1] D.W. Peaceman, H. Rachford, The numerical solution of parabolic and elliptic differential equations, J. Soc. Indust. Appl. Math. 3 (1959) 28–41. [2] Jim Douglas Jr, D.W. Peaceman, Numerical solution for two-dimensional heat flow problems, Am. Inst. Chem. Eng. J. 1 (1955) 505–512. [3] Jim Douglas Jr, J.E. Gunn, A general formulation of alternating direction methods: Part I, parabolic and hyperbolic problems, Numer. Math. 6 (1964) 428–453. [4] Jim Douglas Jr, Seongjai Kim, Improved accuracy for locally one-dimensional methods for parabolic equations, Math. Models Methods Appl. Sci. 11 (9) (2001) 1563–1579. [5] Seongjai Kim, Hyeona Lim, High-order schemes for acoustic waveform simulation, Appl. Numer. Math. 57 (2007) 402–414. [6] S.K. Lele, Compact finite difference schemes with spectral-like solution, J. Comput. Phys. 103 (1992) 16–42. [7] J. Zhang, An explicit fourth-order compact finite difference scheme for three dimensional convection-diffusion equation, Commun. Numer. Meth. Eng. 14 (1998) 209–218. [8] Weidong Dai, Raja Nassar, Compact ADI method for solving parabolic differential equations, Numer. Methods Partial Differential Equations 18 (2) (2002) 129–142. [9] M. Li, T Tang, B. Fornberg, A compact fourth-order finite difference scheme for the steady incompressible Navier–Stokes equations, Internat. J. Numer. Methods Fluids 20 (1995) 1137–1151. [10] Zhang Jun, Multigrid method and fourth-order compact scheme for 2D Poisson equation with unequal mesh-size discretization, J. Comput. Phys. 179 (2002) 170–179. [11] Jichun Li, Yitung Chen, Guoqing Liu, High-order compact ADI methods for parabolic equations, Comput. Math. Appl. 52 (2006) 1343–1356. [12] Cheng aijie, Improvement of stability and convergence for Douglas scheme in two space variables, Numer. Math. J. Chinese Univ. 20 (3) (1998) 265–272 (in Chinese). [13] Sun Zhizhong, Numerical Methods for Differential Equations, China Science Press, 2005 (in Chinese).