A new fourth-order numerical algorithm for a class of nonlinear wave equations

Applied Numerical Mathematics 62 (2012) 1864–1879

Contents lists available at SciVerse ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

A new fourth-order numerical algorithm for a class of nonlinear wave equations ✩ Dingwen Deng a,b , Chengjian Zhang a,∗ a b

School of Mathematics and Statistics, Huazhong University of Science and Technology, Wuhan 430074, People’s Republic of China College of Mathematics and Information Science, Nanchang Hangkong University, Nanchang 330063, People’s Republic of China

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 3 March 2011 Received in revised form 22 March 2012 Accepted 30 July 2012 Available online 28 August 2012 Keywords: Wave equation ADI method Compact finite difference Convergence

In this paper, a new three-level compact alternating direction implicit (ADI) difference scheme is derived for solving a kind of nonlinear wave equations. Basing on a fourthorder approximation to the exact solution at the first time level, it is shown by the energy method that the numerical solution is conditionally convergent with an order of O (t 2 + h4x + h4y ) in H 1 - and L ∞ -norms. A new Richardson extrapolation formula based on three time-grid parameters is given to get numerical solution of fourth-order accuracy in both time and space. The performance of the new algorithm is illustrated by numerical experiments. © 2012 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction In this article, we are concerned with the numerical approximation to the following initial boundary value problem (IBVP) of a class of nonlinear wave equations

∂ 2u ∂u ∂ 2u ∂ 2u +ρ + β 2 u = a 2 + b 2 + f (u , x, y , t ), 2 ∂t ∂t ∂x ∂y u (x, y , 0) = φ(x, y ),

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

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

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

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

¯ (x, y ) ∈ Ω,

(1.1) (1.2) (1.3)

where Ω is a rectangular domain, namely Ω = (a1 , b1 ) × (a2 , b2 ), ∂Ω is the boundary of Ω and Ω¯ = Ω ∪ ∂Ω . φ(x, y ), ϕ (x, y ), ψ(x, y , t ) and f (u , x, y , t ) are sufficiently smooth functions to keep convergence rate and consistency of the difference scheme under consideration. ρ  0, β, a > 0 and b > 0 are all constants. Special cases of Eq. (1.1) are frequently encountered in a vast array of fields. For instance, taking a = b = 1, β > 0 and f (u , x, y , t ) = 0, Eq. (1.1) reduces to telegraph equation (TE) which describes a wide range of phenomena (cf. [23]), such as propagation of electromagnetic waves in superconducting media, propagation of pressure waves occurring in pulsatile blood flow in arteries and a planar random motion of the particle in fluid flowing. When a = b = 1 and ρ = β = 0, Eq. (1.1) becomes a typical nonlinear Klein–Gordon equation (KGE). Taking a = b = 1, β = 0 and f (u , x, y , t ) = − (x, y ) sin u, Eq. (1.1) belongs to the type of the damped sineGordon equation (SGE). KGE and SGE also model many physical phenomena (cf. [2,4,12]) including propagation of fluxion in Josephson junctions between two superconductors, interaction of solitons in a collisionless plasma and so on. Analytical ✩

*

This work is supported by NSFC (Grant No. 11171125, 91130003) and HSFH (Grant No. 2011CDB289), SRF of NHU (No. Ec200907255). Corresponding author. E-mail addresses: [email protected] (D. Deng), [email protected] (C. Zhang).

0168-9274/$36.00 © 2012 IMACS. Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.apnum.2012.07.004

D. Deng, C. Zhang / Applied Numerical Mathematics 62 (2012) 1864–1879

1865

solutions to KGE and SGE were investigated in [12,13,17], respectively. Recently, numerical studies on TE (cf. [5,10,23,24]), KGE (cf. [1,6,16,26,25]) and SGE (cf. [2–4,8,28,30]) have also attracted much attention. Over the years, researchers have devoted much attention to the development and application of higher-order compact (HOC) schemes. Compared with standard difference method, they have advantages of their own, such as high accuracy, compactness and better resolution for high-frequency waves [18] and have been more extensively applied in many fields, such as finance [22], fluid dynamics [15,27], quantum mechanics [3,4,11], biology [7,19]. Operator-splitting methods such as ADI and local one-dimensional (LOD) methods [9,14,21,29] have proven to be very useful in the approximation of the solutions of multi-dimensional hyperbolic problems. More recently, HOC ADI methods, which preserve the high accuracy of HOC schemes and the high efficiency of ADI methods, have been successfully applied to solve hyperbolic problems. For example, in [4], Cui has developed a HOC ADI method for two-dimensional (2D) generalized SGE. This scheme is of order 2 in time and order 4 in space. A class of unconditionally stable HOC ADI methods have been derived for multi-dimensional TEs in [23]. These methods have fourth-order accuracy in space, but only second-order accuracy in time. To further raise computational efficiency, the applications of Richardson extrapolations to HOC schemes are really good alternatives. For example, a combination of HOC ADI method with a global Richardson extrapolation was effectively utilized to solve linear parabolic equation in [20]. However, the study on combinations of HOC ADI schemes with Richardson extrapolations for hyperbolic problems is quite limited. In this paper, enlightened by the work of [4], a new three-level, second-order in time and fourth-order in space, HOC ADI scheme for solving IBVP (1.1)–(1.3) is designed. Then a Richardson extrapolation based on three time-grid parameters is established to make final solution fourth-order accurate in both time and space. By the discrete energy method, various error estimates are given. We should mention that a two-grid Richardson extrapolation in time on a second-order scheme cannot yield fourth-order accuracy, even though the truncation error of the new ADI scheme has a temporal truncation error in the form of O (t 2 ) + O (t 4 ). In fact, because asymptotic expansion of approximate scheme at the first time level includes odd powers of time-grid size, a new Richardson extrapolation formula based on three time-grid parameters is introduced. The paper is organized as follows. In Section 2, we define notations associated to spatial and temporal grids. Section 3 is devoted to the construction of HOC ADI method. Convergence analysis is discussed in Section 4. A new Richardson extrapolation formula based on three time-grid parameters is established in Section 5. Three numerical examples are presented to test the performance of our new algorithms in Section 6. Conclusions are given in Section 7. 2. Partition and notations For temporal discretization, t is temporal increment and there exist two positive integers N, n such that t = T / N, 1

tn = nt, 0  n  N. ∀ v n ∈ { v n | 0  n  N } and letting δt v n−1/2 = ( v n − v n−1 )/t and v n+ 2 = ( v n + v n+1 )/2, we write: δt2 v n = (δt v n+1/2 − δt v n−1/2 )/t, μt v n = ( v n+1 + v n−1 )/2, D t v n = (δt v n+1/2 + δt v n−1/2 )/2. For spatial approximation, let h x = (b1 − a1 )/ L 1 , h y = (b2 − a2 )/ L 2 be grid spacing in x- and y-dimensions, respectively, where L 1 and L 2 are positive integers. Write xi = a1 + ih x , 0  i  L 1 ; y j = a2 + jh y , 0  j  L 2 . The discrete grid Ωh = {(xi , y j ) | 1  i  L 1 − 1, 1  j  L 2 − 1} whose discrete boundary ∂Ωh = {(xi , y j ) | i = 0, L 1 ; j = 0, L 2 }, and Ω¯ h = Ωh ∪ ∂Ωh . ¯ h }, setting δx v 1 = ( v i , j − v i −1, j )/h x and δx2 v i , j = (δx v 1 − δx v 1 )/h x , For grid function v = { v i , j | (xi , y j ) ∈ Ω i− , j i+ , j i− , j 2

2

2

we define: δ y δx v i − 1 , j − 1 = (δx v i − 1 , j − δx v i − 1 , j −1 )/h y , δ y δx2 v i , j − 1 = (δx2 v i , j − δx2 v i , j −1 )/h y . Notations δ y v i , j − 1 , δ 2y v i , j , 2

2

2

2

2

2

δx δ y v i − 1 , j − 1 , δx δ 2y v i − 1 , j and δx2 δ 2y v i , j can be defined similarly as well. 2

2

2

We define a vector space Sh0 as follows: Sh0 = {v | v = ( v 0,0 , v 1,0 , . . . , v L 1 −1,0 , v L 1 ,0 , v 0,1 , v 1,1 , . . . , v L 1 −2, L 2 , v L 1 −1, L 2 ,

v L 1 , L 2 ) T , and v i , j = 0 if (xi , y j ) ∈ ∂Ωh }. That is, a vector in Sh0 can be viewed as a grid function with zero values on ∂Ωh . Obviously, Sh0 ⊂ R( L 1 +1)×( L 2 +1) . ∀v, w ∈ Sh0 , we define inner products as follows: L 1 −1 L 2 −1

(v, w) =

  i =1

(δx v, δx w) =

v i, j w i, j hxh y ,

j =1

i =1 j =1

1 −1 L 2 −1  L  2   2 δx v i , j w i , j h x h y , δ x v, w =

i =1

Similarly,

(δ 2y v, w)

L 1 L 2 −1  

2

2

1 −1 L 2 −1  L  2 2   2 2 δx δ y v i , j w i , j h x h y . δ x δ y v, w =

j =1

i =1

j =1

and (δ y v, δ y w) can be defined as well. Besides, we introduce

v2 = (v, v),

δx v2 = (δx v, δx v),

1 −1 L 2 −1  2 2 L  2 2 δ v = δx v i , j h x h y , x

i =1

δx δ y v2 =

(δx v i − 1 , j )(δx w i − 1 , j )h x h y ,

L1 L2   (δx δ y v i − 1 , j − 1 )2 h x h y , i =1 j =1

2

2

j =1

L2 1 −1   2 2 L  2 2 δ y δ v = δ y δx v i , j − 1 h x h y . x i =1 j =1

2

1866

D. Deng, C. Zhang / Applied Numerical Mathematics 62 (2012) 1864–1879

Similar symbols δ y v2 , δ 2y v2 , δ y δx v2 and δx δ 2y v2 are also defined. We further denote

|v|2H 1 = δx v2 + δ y v2 ,

v21h = |v|2H 1 + v2 ,

 2  2 |v|2H 2 = δx2 v + 2δx δ y v2 + δ 2y v ,

v∞ =

max

(xi , y j )∈Ωh

| v i , j |,

L 1 −1 L 2 −1

h v2 =

  i =1

(h v i , j )2 h x h y

j =1

where h v i , j = δx2 v i , j + δ 2y v i , j , (xi , y j ) ∈ Ωh . 3. Three-level compact ADI difference scheme Using methods similar to those used in [4], a compact ADI difference scheme is derived for IBVP (1.1)–(1.3) in this ¯ h . Then, grid functions un = {un | (xi , y j ) ∈ Ωh } and u (t ) = section. Write u (xi , y j , tn ) = uni, j , u (xi , y j , t ) = u i , j (t ), (xi , y j ) ∈ Ω i, j n n {u i , j (t ) | (xi , y j ) ∈ Ωh } defined on Ωh form vectors u = [u 1,1 , un2,1 , . . . , unL −1,1 , un1,2 , . . . , unL −2, L −1 , unL −1, L −1 ] T and u(t ) 1

1

2

1

2

= [u 1,1 (t ), u 2,1 (t ), . . . , u L 1 −1,1 (t ), u 1,2 (t ), . . . , u L 1 −2, L 2 −1 (t ), u L 1 −1, L 2 −1 (t )] T , respectively. The approximation to un is denoted by Un = [U 1n,1 , U 2n,1 , . . . , U nL −1,1 , U 1n,2 , . . . , U nL −2, L −1 , U nL −1, L −1 ] T . Besides, letting L and v be a linear operator and a grid 1 1 2 1 2 function defined on Ωh , respectively, we denote Lv = [L v 1,1 , L v 2,1 , . . . , L v L 1 −1,1 , L v 1,2 , . . . , L v L 1 −2, L 2 −1 , L v L 1 −1, L 2 −1 ] T . 2 Defining the difference operators Aχ = 1 + (h2χ /12)δχ (χ = x or y) and applying Taylor series expansion with MacLaurin remainder, we have

Aχ u χ χ (xi , y j , t ) = δχ2 u i , j (t ) + R χ (xi , y j , t ),

(xi , y j ) ∈ Ωh ,

(3.1)

that is 1 2 −1 u χ χ ( x i , y j , t ) = A− χ δχ u i , j (t ) + Aχ R χ (xi , y j , t ),

(xi , y j ) ∈ Ωh ,

(3.2)

where 1 A− χ =1









2 , 1 + h2χ /12 δχ

R χ ( xi , y j , t ) = −

h4χ ∂ 6 u 240 ∂ χ 6

  (xi , y j , t ) + O h6χ .

1 2 −1 2 Let us define operator A = aA− x δx + b A y δ y . Applying (3.2) to Eq. (1.1) on Ωh , we obtain a second-order initial value problem as follows:



u (t ) + ρ u (t ) + β 2 u(t ) = Au(t ) + f(t ) + η(t ), u(t 0 ) = φ, u (t 0 ) = ϕ ,

(3.3)

1 −1 where f i , j (t ) = f (u (xi , y j , t ), xi , y j , t ), and ηi , j (t ) = aA− x R x (xi , y j , t ) + b A y R y (xi , y j , t ), (xi , y j ) ∈ Ωh . Applying Taylor series expansion with MacLaurin remainder yields

  δt2 un + ρ D t un + β 2 μt un = A μt un + f n + gn ,

(3.4)

where f in, j = f (uni, j , xi , y j , tn ) and

g ni, j

=η

n i, j

+ − 

+ O t

6

5t 2 ∂ 4 u 12

∂t4

+ t 2 h4x

−

 n t 4 d4 f + − − + 2 dt 2 180 ∂ t 6 30 ∂ t 5 24 dt 4 i , j ∂t3

ρ t 2 ∂ 3 u 3

+ t 2 h4y



,

t 2 d2 f

7t 4 ∂ 6 u

ρ t 4 ∂ 5 u

(xi , y j ) ∈ Ωh , 1  n  N − 1.

Since

D t un =

t 2

δt2 un + δt un−1/2 ,

μt un =

t 2 2

δt2 un + un ,

(3.5)

after substitution of (3.5) into (3.4) and rearrangement, we have



σ2 − where

t 2 2

     A δt2 un = Aun − ρ δt un−1/2 + β 2 un + f n + gn ,

σ = 1 + (ρ /2)t + (β 2 /2)t 2 . Adding the following perturbation term abt 4  −1 2  −1 2  2 n Ax δx A y δ y δt u = Rnxy 4σ 2

to (3.6) and according to the definition of A, we find that

(3.6)

D. Deng, C. Zhang / Applied Numerical Mathematics 62 (2012) 1864–1879





at 2  −1 2  Ax δx 2σ

σ−

1867



σ−

bt 2  −1 2   2 n  A y δ y δt u 2σ

    n  1 2 −1 2 n n−1/2 + β 2 un + f n + R˜ (u) , = a A− x δx + b A y δ y u − ρ δt u

(3.7)

˜ (u)]n = gn + Rnxy . As δx2 and δ 2y depend only on i and j, respectively, δx2 and δ 2y commute with each other. Similarly, where [R two operators Ax and A y also commute with each other, i.e., Ax A y = A y Ax . Consequently, applying Ax A y to (3.7) leads to



σ Ax −

δx2

2σ



aA y δx2

=



at 2

σ Ay −

+ bAx δ 2y



bt 2 2σ

δ 2y

 

δt2 un

 

n

u + Ax A y −

 





n

ρ δt un−1/2 + β 2 un + f n + R(u) ,

(3.8)

where truncation error



R(u)

n

   n = Ax A y R˜ (u) = O t 2 + h4x + h4y .

(3.9)

Omitting small term [R(u)]n in (3.8) and replacing un by its approximation Un yields



σ Ax −



at 2

δx2

2σ

σ Ay −

bt 2 2σ

δ 2y

 

δt2 Un



      = aA y δx2 + bAx δ 2y Un + Ax A y − ρ δt Un−1/2 + β 2 Un + fˆn ,

(3.10)

where ˆf in, j = f (U in, j , xi , y j , tn ), (xi , y j ) ∈ Ωh . Thus, multiplying t 2 to (3.10) and introducing two intermediate variables Gunn type [9,20] is derived as follows:



σ Ax −



σ Ay −



at 2

δx2

2σ

bt 2 2σ



ωi∗, j and ωi∗∗ , j , an ADI scheme of the Douglas–



w ∗i , j = t 2 aA y δx2 + bAx δ 2y U in, j

     n− 1 + t 2 Ax A y − ρ δt U i , j 2 + β 2 U in, j + f U in, j , xi , y j , tn ,



(3.11)

∗ δ 2y w ∗∗ i, j = w i, j ,

n −1 n U in,+j 1 = w ∗∗ i , j + 2U i , j − U i , j ,

(3.12)

(xi , y j ) ∈ Ωh , 1  n  N − 1.

(3.13)

When Eqs. (3.11)–(3.12) are being solved, we need the following boundary conditions







w ∗∗ = t 2 δt2 ψin,0 , i ,0 w ∗0, j



= t

2

σ Ay −

bt 2 2σ

 δ 2y



w ∗∗ = t 2 δt2 ψin, L , i,L



2



δt2 ψ0n, j



∗

w L 1 , j = t

,

2



2

σ Ay −

bt 2 2σ

δ 2y



 δt2 ψ Ln1 , j ,

(3.14)

which can be obtained from (1.2)–(1.3) and (3.12)–(3.13). ¯ h } and {U 1 | Noting that Eqs. (3.11)–(3.13) is a three-level ADI difference scheme, we need both {U i0, j | (xi , y j ) ∈ Ω i, j

(xi , y j ) ∈ Ω¯ h } to start computation. {U i0, j | (xi , y j ) ∈ Ω¯ h } and {U i1, j | (xi , y j ) ∈ ∂Ωh } can be exactly solved by (1.2)–(1.3). An

application of Taylor expansion with integral remainder gives

u 1i , j = u 0i , j +

1 4 3  t 4 t k ∂ 3 u ∂ u ( x , y , 0 ) + (xi , y j , st )(1 − s)3 ds. i j k! ∂ t 3 6 ∂t4 k =1

(3.15)

0

By (1.1)–(1.3), we can exactly compute utt (xi , y j , 0) and uttt (xi , y j , 0). Inserting them, φ(xi , y j , 0) and and then neglecting small term in (3.15) leads to

 U i1, j

=

β 2 t 2

1−

 + +

2

bt 2

t 3 6

2

−

+

ρ β 2 t 3

bρ t 6

6

3

     ρ t 2 ρ 2 t 3 β 2 t 3 at 2 aρ t 3 φ + t − φxx + − ϕ+ − 2

φyy +

f u (φ, x, y , 0)ϕ +

ϕ (xi , y j , 0) into (3.15)



at

3

6

t 2 2

−

ϕxx + ρ t 3 6

bt



6

6

3

ϕyy +

6

t 6

2

6

3

f t (φ, x, y , 0)



f (φ, x, y , 0) , i, j

(xi , y j ) ∈ Ωh .

(3.16)

1868

D. Deng, C. Zhang / Applied Numerical Mathematics 62 (2012) 1864–1879

¯ h } and {U 1 | (xi , y j ) ∈ Ω¯ h } are first obtained. Then, we run the x-sweep Thus, by (1.2)–(1.3) and (3.16), {U i0, j | (xi , y j ) ∈ Ω i, j

¯ h } and { w ∗∗ | (xi , y j ) ∈ Ω¯ h } according to (3.11) and (3.12), respectively. and y-sweep procedures to get { w ∗i , j | (xi , y j ) ∈ Ω i, j

¯ h } by (3.13). It should be mentioned that the {U n+1 | (xi , y j ) ∈ Ω¯ h } is solvable Finally, we determine {U in,+j 1 | (xi , y j ) ∈ Ω i, j solely because of symmetric positive-definite properties of coefficient matrices. 4. Convergence analysis

In this section, various error estimates are given by the discrete energy method. We firstly introduce several lemmas used later. Lemma 1. (Cf. [11].) For any grid function v ∈ Sh0 , we have that

v2 

(b1 − a1 )2 4

δx v2 ,

v2 

(b2 − a2 )2 4

v21h  c˜ |v|2H 1 ,

δ y v2 ,

where c˜ = 1 + |Ω|/8 and |Ω| = (b1 − a1 )(b2 − a2 ). Lemma 2. (Cf. [4,11].) For any grid function v ∈ Sh0 , we have

2

(1)

3 4

(2)

9

2

v2  (Ax v, v)  v2 ,

3

v2  (Ax A y v, v)  v2 ,

v2  (A y v, v)  v2 , 1 3

h2x δx v2  4v2 ;

|v|2H 1  (Ax A y v, −h v),

h2y δ y v2  4v2 .

Lemma 3. (Cf. [20,19].) For any grid function v ∈ Sh0 , we have that |v| H 2 = h v and there exists a positive constant cˆ such that v∞  cˆ h v. Lemma 4. (Cf. [20].) For any time sequences { v 0 , v 1 , . . . , v N } and { g 0 , g 1 , . . . , g N }, we have

2t

l 

g

n



Dt v

n =1

n







1 

ε

v

1 2

2

+ t

l −1  

v

 n+ 12 2



+ v

 l+ 12 2



 +ε



g

 1 2

+ t

n =1

l −1  

δt g

 n+ 12 2

 l 2

+ g

 ,

n =1

for any ε > 0. Write Aˆ h = Ax A y , Ah = aA y δx2 + bAx δ 2y and Bh = δx2 δ 2y . Let eni, j = uni, j − U in, j , G ni, j = f (uni, j , xi , y j , tn ) − f (U in, j , xi , y j , tn ), (xi , y j ) ∈ Ω¯ h , 0  n  N. Then corresponding vectors en , Gn ∈ Sh0 . Subtracting (3.10) from (3.8) and (3.16) from (3.15), respectively, one gets

⎧   abt 4  2 n     n ⎪ ⎪ Bh δt e = Ah μt en + Aˆh Gn + R (u) , ⎪ Aˆh δt2 en + ρ D t en + β 2 μt en + ⎪ 2 ⎪ 4 σ ⎪ ⎨ 1  n  N − 1 ; e0 = 0 ; 1 4 ⎪ ⎪ t 4 ∂ u ⎪ 1 ⎪ e = (xi , y j , st )(1 − s)3 ds, (xi , y j ) ∈ Ωh . ⎪ i, j ⎪ ⎩ 6 ∂t4

(4.1)

0

4.1. H 1 -norm error

¯ × [0, T ] and |˜ε |  ε˜ 0 , it holds that We firstly assume there are positive constants c 1 , ε˜ 0 , so that, for any (x, y , t ) ∈ Ω

f (u + ε˜ , x, y , t ) − f (u , x, y , t )  c 1 |˜ε |.

(4.2)

Meanwhile, setting h = max{h x , h y }, we also assume that there exists constant c 2 > 0 such that c 2 h  h x , h y  h. Noting that 1/ 2

δt e i , j = (e 1i , j − e 0i , j )/t, we assume that there are four positive constants c 3 , c 4 , c 5 and c 6 , such that

  1  1/2       R(u) n 2  c 2 t 4 + h8 + h8 , e   c 4 t 4 , δt e   c 4 t 3 , x y 3  1 1/2  e 1  c 5 t 4 , δt e 1  c 5 t 3 , δx δ y δt e1/2   c 6 t 3 . H H In what follows, we apply mathematical induction to prove Theorem 1.

(4.3)

D. Deng, C. Zhang / Applied Numerical Mathematics 62 (2012) 1864–1879

1869

¯ h , 0  K  N } be the solution of the difference scheme (3.11)–(3.13) and (3.16) at Theorem 1. Let grid function {U iK, j | (xi , y j ) ∈ Ω

time level K , {u iK, j | (xi , y j ) ∈ Ω¯ h , 0  K  N } be exact solution of IBVP (1.1)–(1.3) at time t = K t. Then under assumptions of (4.2) 1

and t  ch 2 + ( > 0), we arrive at



      δt e K + 12 2 + e K 2  C t 4 + h8 + h8 x y 1h

4

max

9

0 K  N

(4.4)

for N t = T , where

 C = max

1+

2

+



c 7 + c 32 T exp(c 8 T ),

min{a, b}

β 2 c 42 t 4

c 7 = c 42 t 2 +

 

3c˜

abc 62 t 6 4 2

σ

+

4c 42 t 2 9

max{a, b}c 52 t 4 2

   + c 42 + c 52 t 4 ,

,

c8 =

3c 12 (a1 − b1 )2 4a

9

+ . 2

Proof. Obviously, (4.4) is valid for K = 0, 1. Now suppose that (4.4) is valid for K = 0, 1, . . . , l (2  l  N − 1). Next, it is shown that (4.4) will be valid for K = l + 1 as well. Noting that t  ch1/2+ , it follows from induction hypothesis that 1 1 n  1 − 12  n  − 12   2 2 e  h− e   h − C t 4 + h8x + h8y 2  ε˜ 0 , x hy x hy i, j

1  n  l,

(4.5)

only if t, h x and h y are enough small. So, a combination of assumption (4.2) with (4.5) yields that |G ni, j |  c 1 |eni, j |, 1  n  l, which implies that



2 1  2 c 2  2  1 Aˆh Gn , D t en  δt en+1/2  + δt en−1/2  + 1 en  ,

1  n  l.

(4.6)

2     δt en+1/2 2 + 1 δt en−1/2 2 + c 3 t 4 + h8 + h8 . x y

(4.7)

4

4

2

It follows from Cauchy–Schwarz inequality and (4.3) that





n

R(u) , D t en 

1 4

4

2

Writing





F n = Aˆh δt en+1/2 , δt en+1/2 +

+

β 2  2

   a     Aˆh en+1 , en+1 + Aˆh en , en + A y δx en+1 , δx en+1 + A y δx en , δx en 2

abt 4 

     δx δ y δt en+1/2 2 + b Ax δ y en+1 , δ y en+1 + Ax δ y en , δ y en ,

4σ 2

(4.8)

2

and applying Lemmas 1 and 2, we find that

 n 2 (b1 − a1 )2  n 2 e   δx e   3 (b1 − a1 )2 F n , 4

4 9

 δt en+1/2 2  F n ,

 n +1  2 e   1h

4a



3c˜ min{a, b}

F n,

(4.9)



F 0  c 7 t 4 + h8x + h8y .

(4.10)

Taking inner product with the first equation of (4.1) by D t en h x h y , and then using the discrete Green formula and (4.6)–(4.8) gives

4ρ 

2    D t en 2 + 1 F n  1 F n−1 + c 1 en 2 + 9 2t 2t 2

1 2

2     δt en+1/2 2 + 1 δt en−1/2 2 + c 3 t 4 + h8 + h8 . x y

2

2

(4.11)

Multiplying both sides of (4.11) by 2t, summing n from 1 to l and using the discrete Gronwall’s lemma yields









F l  c 7 + c 32 T exp(c 8 T ) t 4 + h8x + h8y ,

(4.12)

where (4.9) and (4.10) have been used. Furthermore, it follows from (4.9)–(4.10) that

4 9

     δt el+1/2 2 + el+1 2  C t 4 + h8 + h8 . x y 1h

It is proved that (4.4) is valid for K = l + 1. The proof is completed.

(4.13)

2

1870

D. Deng, C. Zhang / Applied Numerical Mathematics 62 (2012) 1864–1879

4.2. Maximum norm error To derive the maximum norm error, we should give another two hypotheses.

¯ × [0, T ]), we assume that there exists constant (1) Setting ∂ u (u (x, y , t ), x, y , t ) ∈ L ∞ (Ω ∂f

μ1 > 0, such that

 ∂f u (x, y , t ), x, y , t  μ1 . ¯ 0, T ] ∂ u (x, y ,t )∈Ω×[ ess sup

(4.14)

μ2 and ε0 , so that, for any (x, y , t ) ∈ Ω¯ × [0, T ] and |εk |  ε0 , k = 1, 2,   f u (u + ε1 , x, y , t + ε2 ) − f u (u , x, y , t )  μ2 |ε1 | + |ε2 | . (4.15)

(2) Also, suppose that there are positive constants

Moreover, let us assume that there are two positive constants

  h e1   μ3 t 4 ,

  h e 12   0.5μ3 t 4 ,

μ3 and μ4 , such that

     1 δt R(u) n+ 2 2  μ2 t 4 + h8 + h8 . x

4

y

(4.16)

Now, we give the mathematical proof of the following theorem.

¯ h , 0  K  N } be the solution of the difference scheme (3.11)–(3.13) and (3.16) at Theorem 2. Let grid function {U iK, j | (xi , y j ) ∈ Ω ¯ h , 0  K  N } be exact solution of IBVP (1.1)–(1.3) at time t = K t. Then under assumptions of (4.2) time level K , {u iK, j | (xi , y j ) ∈ Ω 1

(4.14), (4.15) and t  ch 2 + ( > 0), the following error estimate holds

max

0 K  N

     h e K 2  M t 4 + h8 + h8 x

(4.17)

y

for N t = T , where

M=



exp( T )



72

c 12 C + c 32 +

min{a, b} min{a, b}

+

10β 2 c 52 t 4

+

3

7 2

μ24 T 2

max{a, b}μ23 t 4 +

+

9C μ25 T

6abc

 +

4

20c 52 t 2

2 2 4 3 t h 2 +4  c2

4

μ

3

 ,

and μ5 > 0 is only dependent on μ1 , μ2 , ε0 and |u (x, y , t )| (see (4.18)). 1

Proof. By Theorem 1 and t  ch 2 + , and applying the same method as (4.5), we get

n e  ε0 , i, j

n +1 e  ε0 , i, j

 n +1    u − uni, j − eni ,+j 1 − eni, j s  ε0 , i, j

n+ 12  ε0 δt e i, j

provided t, h x and h y are small enough. Combining (4.14) with (4.15) gives n+ 1 δt G i , j 2

=

n+ 1  δt e i , j 2



1





f u uni ,+j 1 − eni ,+j 1 s, xi , y j , tn+1 ds

0

+

eni, j

1

t









f u uni ,+j 1 − eni ,+j 1 s, xi , y j , tn+1 − f u uni, j − eni, j s, xi , y j , tn



ds

0

 n+ 1   μ5 δt e i , j 2 + eni, j , where

μ5 is only dependent on μ1 , μ2 , ε0 and |u (x, y , t )|. Thus, Theorem 1 implies that   n + 1 2       δt G 2   2μ2 δt en+ 12 2 + en 2  9 C μ2 t 4 + h8 + h8 . x y 5 5 2

(4.18)

A combination of Theorem 1 with assumption (4.2) leads to

   n 2   G   c 2 en 2  c 2 C t 4 + h8 + h8 , 1

Next, we introduce

1

x

y

1  n  N.

(4.19)

D. Deng, C. Zhang / Applied Numerical Mathematics 62 (2012) 1864–1879



2

2

T n+1 = en+1 H 1 + en H 1 −

−

    δ 2 en+1 2 + δ 2 en 2 − x x

h2x  12

 2 n+1 2  2 n 2  δ e  + δ e  +

h2y

y

12

y

h2x h2y 144



h2x 12



+

1871

    δx δ y en+1 2 + δx δ y en 2

h2y  12

 2 n+1 2  2 n 2  2 n+1 2  2 n 2  δ y δ e  + δ y δ e  + δx δ e  + δx δ e  , x

x

y

y

(4.20)

   h2       δ 2 δt en+ 12 2 + δx δ y δt en+ 12 2 − y δx δ y δt en+ 12 2 + δ 2 δt en+ 12 2

h2x 

S n +1 = −

x

12

y

12

    h2x h2y  1 2 δ y δ 2 δt en+ 12 2 + δx δ 2 δt en+ 12 2 , + δt en+ 2 H 1 + x y

(4.21)

144





 









 







 

2 2 2 2 2 2 F n+1 = a δx2 en+1  + δx2 en  + b δ 2y en+1  + δ 2y en  + (a + b) δx δ y en+1  + δx δ y en 

 − S˜ n+1 =



ah2y

+

12

        δ y δ 2 en+1 2 + δ y δ 2 en 2 + δx δ 2 en+1 2 + δx δ 2 en 2 ,

bh2x 

x

12

x

y

(4.22)

y

    δ y δ 2 δt en+ 12 2 + δx δ 2 δt en+ 12 2 . x y

abt 4  4σ 2

(4.23)

Clearly, it follows from Lemma 2 that

   1  n+1 2 e 1 + en 2 1  T n+1  10 en+1 2 1 + en 2 1 , H H H H 3 9

(4.24)

1 3 2 3

δt en+ 12 2 1  S n+1  10 δt en+ 12 2 1 , H H

(4.25)

9







 







 

2 2 2 2 min{a, b} h en+1  + h en   F n+1  max{a, b} h en+1  + h en  .

(4.26)

Taking inner products between the first equation of (4.1) and −4(t )h x h y h D t en , 1  n  K < N, summing up over n, and using discrete Green formula, Lemma 2 and Lemma 4 with ε = 3/ min{a, b}, we deduce that

2S K +1 +

K 4ρ t 

D t en 2 1 + β 2 T K +1 + F K +1 + 2 S˜ K +1

3

1

H

n =1

2

1



˜1

1

 2S + β T + F + 2 S +

+

3



    h e K +1 2 + h e K 2 +

1  2

K −1        h e 12 2 + t h en+1 2 + h en 2

2 min{a, b} 

2



12 min{a, b}

n =1

 1 2     G  +  R(u) 1 2

  n+ 1 2   n+ 12 2   K 2   K 2  + G  +  R(u)  . δt G 2  + δt R(u)

K −1 

+ t

(4.27)

n =1

Therefore, by (4.24)–(4.27), we have that



1 min{a, b}



 K  K + 1 2 2 2  2 K 2    K +1 2 n 2 K + 1 D t e 1 + β e 1 + e 1 + 6 S˜ 2 δt e 2 H 1 + 4ρ t + h e K +1  + h e K  H H H

1 min{a, b}

n =1

20 3

2        δt e 12 2 1 + 10β e1 2 1 + e0 2 1 + 3 max{a, b} h e1 2 + h e0 2 + 6 S˜1 H H H



3

  K −1      1 2 t   n +1  2 n 2     2 + 2 h e + + h e h e 2

+

36

(min{a, b})2

n =1

  K −1   1 2    n+ 1 2      1      G  +  R(u) 1 2 + t δt G 2  + δt R(u) n+ 2 2 + G K 2 +  R(u) K 2 n =1

1872

D. Deng, C. Zhang / Applied Numerical Mathematics 62 (2012) 1864–1879



K −1 

     h en+1 2 + h en 2 t +

n =1

20c 52 t 2

+

3

10β 2 c 52 t 4

+

3

+

7 2

1



72

min{a, b} min{a, b}

 c 12 C

+ c 32

+ 

max{a, b}μ23 t 4 +

μ24 T 2

+

9C μ25 T



4

 t 4 + h8x + h8y .

6abc 4 μ23 t 2 h4  c 22+4

(4.28)

Hence, an application of the discrete Gronwall’s lemma to (4.28) yields

      h e K +1 2 + h e K 2  M t 4 + h8 + h8 . x

(4.29)

y

We obtain Theorem 2 because K  N − 1 is an arbitrary positive integer.

2

By Lemma 3, we immediately obtain Theorem 3. Theorem 3. Let u (x, y , t ) ∈ C (6,4) (Ω¯ × [0, T ]) be the exact solution of IBVP (1.1)–(1.3). Then under conditions of Theorem 2, the numerical solution of the new ADI method (3.11)–(3.13) and (3.16) is convergent with an order of O (t 2 + h4x + h4y ) in L ∞ -norm. 5. Improvement of accuracy in temporal dimension In fact, a more precise bound in (4.29) is









˜ t 2 + t 3 + h4x + h4y , max h e K   M

0 K  N

which implies









max e K ∞ = O t 2 + t 3 + h4x + h4y ,

0 K  N

˜ > 0. Therefore, to get the numerical solution of order 4 in time, a three-grid Richardson extrapolation where constant M must be established. ¯ × [0, T ]) be the exact solution of IBVP (1.1)–(1.3), and {U K (t , h x , h y ) | Theorem 4. Let smooth function u (x, y , t ) ∈ C (6,6) (Ω i, j

(xi , y j ) ∈ Ω¯ h , 0  K  N } be the numerical solution of ADI method (3.11)–(3.13) and (3.16) at time t = K t. Then, under conditions of Theorem 2, we find that









max u K − (U E ) K ∞ = O t 4 + h4x + h4y ,

(5.1)

1 K  N

¯ h , and the extrapolation solution at time level K is defined by where u iK, j = u (xi , y j , K t ), (xi , y j ) ∈ Ω



(U E )iK, j

=

1 U K (t , h x , h y ) − 12 U 2K ( t , h x , h y ) + 32 U 4K ( t , h x , h y ), 21 i , j 21 i , j 2 21 i , j 4

ψ(xi , y j , t K ),

(xi , y j ) ∈ ∂Ωh , 1  K  N .

(xi , y j ) ∈ Ωh , 1  K  N ,

(5.2)

Proof. Setting

r (x, y , t ) = −

ρ ∂ 3u 3 ∂t3

−

1 ∂ 4u 2 ∂t4

+

1 d2 f 2 dt 2

,

hˆ (x, y , t ) =

1 ∂ 4u 12 ∂ t 4

+ r (x, y , t ),

it is clear from (3.9) that



n





R (u ) i , j = t 2 Aˆh hˆ (xi , y j , tn ) + O t 4 + h4x + h4y ,

(xi , y j ) ∈ Ωh .

(5.3)

We assume that functions p (x, y , t ) and q(x, y , t ) satisfy the following two IBVPs

⎧ 2 ˆ ⎪ ⎨ ptt + ρ pt + β p = ap xx + bp y y + H (x, y , t ), (x, y ) ∈ Ω, 0 < t  T , ¯ p (x, y , 0) = 0, pt (x, y , 0) = 0, (x, y ) ∈ Ω, ⎪ ⎩ p (x, y , t ) = 0, (x, y ) ∈ ∂Ω, 0 < t  T , ⎧ 2 ˜ ⎪ ⎨ qtt + ρ qt + β q = aq xx + bq y y + H (x, y , t ), (x, y ) ∈ Ω, 0 < t  T , ¯ q(x, y , 0) = 0, qt (x, y , 0) = −0.5r (x, y , 0), (x, y ) ∈ Ω, ⎪ ⎩ q(x, y , t ) = 0, (x, y ) ∈ ∂Ω, 0 < t  T ,

(5.4)

(5.5)

D. Deng, C. Zhang / Applied Numerical Mathematics 62 (2012) 1864–1879

1873

ˆ (x, y , t ) = f u (u , x, y , t ) p + hˆ (x, y , t ) and H˜ (x, y , t ) = f u (u , x, y , t )q. Define grid function pn = respectively, where H i, j n ¯ h , 0  n  N. Then corresponding vectors pn , qn ∈ S0 . Similar to (3.8), we can p ( xi , y j , t n ) , q i , j = q ( xi , y j , t n ) , ( xi , y j ) ∈ Ω h easily develop their equivalent difference equations as follows:

⎧   abt 4  2 n      ⎪ ˆ n + R(p) n , ⎪ Aˆh δt2 pn + ρ D t pn + β 2 μt pn + Bh δt p = Ah μt pn + Aˆh H ⎪ 2 ⎪ 4σ ⎪ ⎪ ⎨ 1  n  N − 1; p0 = 0; 1 ⎪ ⎪ t 2 ˆ t 3 ⎪ 1 ⎪ ⎪ p = h ( x i , y j , 0) + p (3) (xi , y j , st )(1 − s)2 ds, (xi , y j ) ∈ Ωh , ⎪ ⎩ i, j 2 2

(5.6)

0

⎧   abt 4  2 n      ⎪ ˜ n + R(q) n , ⎪ Aˆh δt2 qn + ρ D t qn + β 2 μt qn + Bh δt q = Ah μt qn + Aˆh H ⎪ 2 ⎪ 4σ ⎪ ⎪ ⎨ 1  n  N − 1; q0 = 0; 1 ⎪ ⎪ ⎪ 1 2 ⎪ ⎪ q = − 0 . 5r ( x , y , t ) t +  t q(2) (xi , y j , t 0 + st )(1 − s) ds, (xi , y j ) ∈ Ωh , i j 0 ⎪ ⎩ i, j

(5.7)

0

ˆ n = f u (un , xi , y j , tn ) pn + hˆ n , H˜ n = f u (un , xi , y j , tn )qn , (xi , y j ) ∈ Ωh . Obviously, [ R ( p )]n = respectively, in which H i, j i, j i, j i, j i, j i, j i, j i, j 2 4 4 O(t + h x + h y ), [ R (q)]ni, j = O(t 2 + h4x + h4y ), (xi , y j ) ∈ Ωh . Define vector sn = en − t 2 pn − t 3 qn , 0  n  N. Clearly, vector sn ∈ Sh0 . Multiplying (5.6) and (5.7) by t 2 and t 3 , respectively, and subtracting the resulting systems from (4.1), it is easy to find that

⎧   abt 4  2 n     n ⎪ ⎪ ⎪ Aˆh δt2 sn + ρ D t sn + β 2 μt sn + Bh δt s = Ah μt sn + Aˆh Hn + R(s) , ⎨ 4σ 2 0 ⎪ 1  n  N − 1; s = 0; ⎪   ⎪ ⎩ s1 = O t 5 , (x , y ) ∈ Ω , i j h i, j

(5.8)

where H ni, j = f (uni, j − t 2 pni, j − t 3 qni, j , xi , y j , tn ) − f (U in, j , xi , y j , tn ), (xi , y j ) ∈ Ωh . By (5.3), we find that



n





n



n

n









R (s) i , j = R (u ) i , j − t 2 Aˆh hˆ ni, j − t 2 R ( p ) i , j − t 3 R (q) i , j + O t 4 = O t 4 + h4x + h4y ,

(xi , y j ) ∈ Ωh .

Hence, using the same technique in the proof of Theorem 2, we obtain easily









max s K ∞  C ∗ t 4 + h4x + h4y ,

1 K  N

where C ∗ is a positive constant and independent of grid parameters t, h x and h y . Namely, there exist two functions p (x, y , t ) and q(x, y , t ) such that





u iK, j = U iK, j (t , h x , h y ) + t 2 p iK, j + t 3 qiK, j + O t 4 + h4x + h4y ,

(xi , y j ) ∈ Ωh .

(5.9)

When the values of time-step size are equal to t /2 and t /4, respectively, we similarly obtain

 u iK, j

=

U i2K ,j

t 2

 u iK, j = U i4K ,j

t 4





t

2



t

3





qiK, j + O t 4 + h4x + h4y ,

(5.10)

  2  3   t t , hx , h y + p iK, j + qiK, j + O t 4 + h4x + h4y .

(5.11)

, hx , h y +

p iK, j

2

4

+

2

4

Multiplying Eqs. (5.9)–(5.11) by 1/21, −12/21, and 32/21, respectively, and adding the results gives





u iKj − (U E )iK, j = O t 4 + h4x + h4y ,

(xi , y j ) ∈ Ωh .

2

(5.12)

Remark. By Theorem 4, Richardson extrapolation algorithm based on two time-grid parameters



(U˜ E )iK, j

=

4 2K t U ( , h x , h y ) − 13 U iK, j (t , h x , h y ), 3 i, j 2

(xi , y j ) ∈ Ωh , 1  K  N ,

ψ(xi , y j , t K ),

(xi , y j ) ∈ ∂Ωh , 1  K  N ,

(5.13)

cannot achieve the convergence order of O (t 4 + h4x + h4y ). This is illustrated numerically in the first example in the next section as well.

1874

D. Deng, C. Zhang / Applied Numerical Mathematics 62 (2012) 1864–1879

Table 1 Computational results for Example 1 at T = 1 with t = h. h

1/4

1/8

1/16

1/32

1/64

Eu1 RE(I) Ru1 t1

1.1 × 10−2 1.8 × 10−2 3.8 × 10−5

1.8 × 10−3 2.61 3.0 × 10−3 2.7 × 10−4

3.1 × 10−4 2.55 5.1 × 10−4 2.2 × 10−3

6.4 × 10−5 2.27 1.1 × 10−4 1.7 × 10−2

1.5 × 10−5 2.12 2.4 × 10−5 1.4 × 10−1

Eu2 RE(II) Ru2 t2

6.6 × 10−4 ∗ 1.1 × 10−3 2.7 × 10−4

3.0 × 10−5 4.45 5.0 × 10−5 1.9 × 10−3

1.7 × 10−6 4.12 2.9 × 10−6 1.5 × 10−2

9.8 × 10−8 4.15 1.6 × 10−7 1.2 × 10−1

2.7 × 10−9 5.20 4.4 × 10−9 9.8 × 10−1

Eu3 RE(III) Ru3 t3

2.8 × 10−3

2.6 × 10−4 3.47 4.2 × 10−4 8.0 × 10−4

2.1 × 10−5 3.63 3.4 × 10−5 6.5 × 10−3

1.9 × 10−6 3.48 3.1 × 10−6 5.2 × 10−2

1.8 × 10−7 3.38 3.0 × 10−7 4.2 × 10−1

∗

∗ 4.7 × 10−3 1.2 × 10−4

6. Numerical experiments In this section, three IBVPs whose analytical solutions are known to us are used to illustrate the practicability of our algorithms. For brevity, the present ADI method is called Algorithm I. Richardson extrapolation algorithm (5.2), which is used in conjunction with the present ADI method to solve IBVP (1.1)–(1.3), is called Algorithm II. Moreover, Richardson extrapolation algorithm (5.13), which is applied along with the present ADI method to solve IBVP (1.1)–(1.3), is called Algorithm III. In addition, we only adopt uniform spatial grids (h x = h y = h) with different meshsizes h. The L ∞ -, L 2 - and maximum relative errors at T = N t, which are defined by u N − U N ∞ , u N − U N  and (u N − UN )/UN ∞ , respectively, are applied to measure the performance of numerical methods. For simplicity, we denote L ∞ -, L 2 - and maximum relative errors of numerical solution and computational cost in seconds from Algorithm I by Eu1 , Lu1 , Ru1 and t 1 , respectively. Similarly, corresponding errors of numerical solution and time cost in seconds from Algorithm II are represented by Eu2 , Lu2 , Ru2 and t 2 , respectively. And, Eu3 , Lu3 , Ru3 and t 3 denote corresponding errors of numerical solution and computational cost in seconds from Algorithm III, respectively. To test the accuracy of the new algorithms, we also report the orders of L ∞ - and L 2 -errors of numerical solution in the following tables. We define them as follows: RE(σ ) = log2 [Euδ (h)/Euδ (h/2)], RL(σ ) = log2 [Luδ (h)/Luδ (h/2)], (σ , δ) = (I, 1) or (II, 2) or (III, 3) when t = h and t = h2 ; and RE(I) = log2 [Eu1 (t , h)/Eu1 (t /4, h/2)], RL(I) = log2 [Lu1 (t , h)/Lu1 (t /4, h/2)] when t = h. All computer programs were coded in Fortran 95. Example 1. The first IBVP to be solved is

∂ 2u ∂ 2u ∂ 2u 1 ∂ u + + 4u = 2 + 2 + f (x, y , t ), 2 2 ∂t ∂t ∂x ∂y

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

where Ω = (0, 1) × (0, 1) and f (x, y , t ) = 2(π 2 + 2) exp(−0.5t ) sin(π x) sin(π y ). The theoretical solution is u (x, y , t ) = exp(−0.5t ) sin(π x) sin(π y ). To accurately assess the performance of three methods, for a fixed grid, we run 100 times, then take the average computational time as time cost (CPU time) in Tables 1–2. The computational results listed in Tables 1–2 show the following conclusions. (1) From Table 1, we can see that Algorithm II has the convergence rate of O (h4 ) in L ∞ -norm, while Algorithm I has only O(h2 ). Meanwhile, we also find that Algorithm III cannot attain convergence order of O(h4 ). These coincide well with our analysis. (2) Table 1 also shows the superiority of Algorithm II over Algorithm I and Algorithm III. For example, with the same meshsize h, Algorithm II provides the most accurate solution. For achieving nearly the same accuracy, computational cost of Algorithm II is also the lowest. For instance, we can see from Table 1 that Eu1 = 1.5 × 10−5 obtained by Algorithm I with h = 1/64 needs 0.14 s; Eu3 = 1.9 × 10−6 obtained by Algorithm III with h = 1/32 requires 5.2 × 10−2 s; however, Eu2 = 1.7 × 10−6 gained by Algorithm II with h = 1/16 only costs 1.5 × 10−2 s. (3) Table 2 shows that Algorithm I is temporally second-order accurate and spatially fourth-order accurate in L ∞ - and L 2 norms. There seems to be a little degradation of the convergence rate for h = 1/64, t = 1/256, where round-off errors probably start to play a role. Also, the same kind of phenomenon exists in other tables. Example 2. Consider the following sine-Gordon equation

∂ 2u ∂ 2u ∂ 2u = 2 + 2 − sin u , ∂t2 ∂x ∂y

−7  x, y  7, 0 < t  T ,

which was studied in [4,8]. The exact solution is u (x, y , t ) = 4 arctan(exp(x + y − t )). Numerical results are displayed in Tables 3–6 and Fig. 1.

D. Deng, C. Zhang / Applied Numerical Mathematics 62 (2012) 1864–1879

1875

Table 2 Numerical results for Example 1 at T = 2. h

1/4

1/8

1/16

1/32

1/64

t

1/2 3.1 × 10−2

1/8 1.8 × 10−3 4.13 8.9 × 10−4 4.13 4.8 × 10−3 5.4 × 10−4

1/32 1.1 × 10−4 3.98 5.6 × 10−5 3.98 3.1 × 10−4 8.6 × 10−3

1/128 6.9 × 10−6 4.03 3.4 × 10−6 4.03 1.9 × 10−5 1.4 × 10−1

1/256 4.5 × 10−7 3.95 2.2 × 10−7 3.95 1.2 × 10−6 2.22

Eu1 RE(I) Lu1 RL(I) Ru1 t1

∗ 1.5 × 10−2 ∗ 8.4 × 10−2 4.0 × 10−5

Table 3 Numerical results for Example 2 at T = 5 with t = h. h

1/4

1/8

1/16

1/32

1/64

Eu1 RE(I) Lu1 RL(I) t1

2.3 × 10−1

4.9 × 10−2 2.23 1.7 × 10−1 2.18 0.18

1.1 × 10−2 2.12 4.0 × 10−2 2.06 1.55

2.8 × 10−3 2.01 9.9 × 10−3 2.01 12.89

7.0 × 10−4 2.00 2.5 × 10−3 2.00 103.83

5.3 × 10−3

2.5 × 10−4 4.39 5.2 × 10−4 4.45 1.37

1.1 × 10−5 4.56 2.5 × 10−5 4.39 11.12

6.0 × 10−7 4.16 1.5 × 10−6 4.08 95.33

4.2 × 10−8 3.82 1.3 × 10−7 3.44 776.79

Eu2 RE(II) Lu2 RL(II) t2

∗ 7.6 × 10−1 ∗ 0.03 ∗ 1.1 × 10−2 ∗ 0.19

Table 4 Numerical results for Example 2 at T = 5 for with t = h2 . h

1/4

1/8

1/16

1/32

1/64

Eu1 RE(I) Lu1 RL(I) t1

1.2 × 10−2

7.1 × 10−4 4.01 2.5 × 10−3 4.02 1.55

4.5 × 10−5 4.00 1.6 × 10−4 4.01 24.77

2.8 × 10−6 4.01 9.8 × 10−6 4.02 410.86

1.5 × 10−7 4.24 5.0 × 10−7 4.29 6634.24

∗ 4.1 × 10−2 ∗ 0.10

Table 5 Numerical convergence of Algorithm I in L ∞ -norm and L 2 -norm for Example 2 at T = 1.

t

h

Eu1

RE(I)

Lu1

RL(I)

t1

1/4 1/16 1/64 1/256 1/1024

1/4 1/8 1/16 1/32 1/64

3.2 × 10−2 2.6 × 10−3 1.7 × 10−4 1.1 × 10−5 6.5 × 10−7

∗ 3.64 3.95 3.99 4.01

1.4 × 10−1 1.0 × 10−2 6.4 × 10−4 4.0 × 10−5 2.5 × 10−6

∗ 3.80 3.96 3.99 4.01

0.00 0.07 1.25 20.72 333.56

The data in Table 3 once again demonstrate that Algorithm I has only convergence rate of O (h2 ), while the convergence rate of Algorithm II is O (h4 ) when t = h. When t = h2 , we can see from Table 4 that convergence orders RE(I) and LE(I) are almost equal to four. However, compared with the data in Table 3, with the same meshsize h, Algorithm II not only provides more accurate solution, but also costs much less CPU time. Table 5 further shows that Algorithm I is of order 2 in time and order 4 in space with respect to L ∞ - and L 2 -norms. To exhibit propagation of pointwise errors (i.e. e iN, j = u iN, j − U iN, j , (xi , y j ) ∈ Ω¯ h ), and behavior of numerical solution, we use Algorithm II with h = t = 0.125 to solve Example 2 and find that Eu2 = 2.7 × 10−4 and Lu2 = 6.8 × 10−4 at T = 7. Meanwhile, Fig. 1 presents numerical solution at T = 7 and corresponding pointwise errors. In Table 6, Eu4 and Lu4 denote L ∞ - and L 2 -errors of numerical solution from the algorithm proposed in [2], respectively; corresponding errors of numerical solution from the algorithm presented in [4] are represented by Eu5 and Lu5 , respectively. The data of Eu5 and Lu5 in Table 6 come from Table 4 of paper [4]. Table 6 shows that Algorithm II generates the best solution. It is worth mentioning that although the algorithm proposed in [4] is also second-order accurate in time and fourth-order accurate in space, numerical results show Algorithm I can give more accurate solution. This meets our anticipation because the splitting error introduced in this paper is O (t 4 ), while that introduced in [4] is O (t 3 ).

1876

D. Deng, C. Zhang / Applied Numerical Mathematics 62 (2012) 1864–1879

Table 6 L ∞ - and L 2 -errors at different time levels provided by four different methods for Example 2, with t = 0.01 and h = 0.1. T

2

3

4

5

6

7

Eu4 Lu4 Eu5 Lu5 Eu1 Lu1 Eu2 Lu2

2.7 × 10−3 8.5 × 10−3 4.0 × 10−4 1.1 × 10−3 1.5 × 10−4 4.8 × 10−4 6.5 × 10−6 1.6 × 10−5

4.0 × 10−3 1.2 × 10−2 5.2 × 10−4 1.5 × 10−3 2.2 × 10−4 6.9 × 10−4 7.9 × 10−6 2.0 × 10−5

4.8 × 10−3 1.6 × 10−2 5.4 × 10−4 1.8 × 10−3 2.7 × 10−4 9.0 × 10−4 8.0 × 10−6 2.2 × 10−5

5.3 × 10−3 1.9 × 10−2 5.7 × 10−4 2.0 × 10−3 2.9 × 10−4 1.0 × 10−3 8.2 × 10−6 2.4 × 10−5

5.7 × 10−3 2.0 × 10−2 6.2 × 10−4 2.0 × 10−3 3.2 × 10−4 1.1 × 10−3 8.8 × 10−6 2.5 × 10−5

6.4 × 10−3 1.9 × 10−2 7.0 × 10−4 1.9 × 10−3 3.6 × 10−4 1.1 × 10−3 9.6 × 10−6 2.5 × 10−5

Fig. 1. Example 2 (solved by algorithm II with h = t = 0.125): numerical solution at T = 7 (left column) and corresponding pointwise errors (right column).

Example 3. The final problem to be solved is

 2  ∂ 2u ∂ 2u 2 ∂ u + γ u3 + α u = α + ∂t2 ∂ x2 ∂ y2

√

on

the region −7  x, y  7 and 0 < t  T . Its exact solution is u (x, y , t ) = α /γ tanh(κ (x + y − ct )), where κ = α /[2(c 2 − 2α 2 )] with α , γ > 0 and any c such that c 2 > 2α 2 . This equation is called kink wave. A vital property of Example 3 is the conservation of energy,

E = E (t ) =

1



2

   γ (ut )2 + α 2 (u x )2 + (u y )2 + α u 2 − u 4 dx dy , 2

R2

namely E (t ) = E (0), ∀t > 0. Therefore, the conserved quantity can be applied to verify the efficiency of numerical method. Let Ω = (−7, 7) × (−7, 7), then

E = E (t ) =

1 2



   γ (ut )2 + α 2 (u x )2 + (u y )2 + α u 2 − u 4 dx dy ,

(6.1)

2

Ω¯

may be regarded roughly as approximate conservation. To confirm the conservation of energy, we introduce average operator μ¯ uni, j = (uni, j + uni+1, j + uni, j+1 + uni+1, j+1 )/4 and apply the following composite trapezoidal rule

 L −1 L −1

 L 1 −1 L 2 −1 1 2       n 2 γ  n 4 2   n− 12 2 E = μ¯ δt u i, j + α u i, j − u i, j hxh y + α2 μ¯ δx uni− 1 , j h x h y n

1 2

+

+

i =0

L 2 −1  α2 

2

j =0

L 1 −1  α2 

2

2

j =0

i =0

δx unL

2 1 1− 2 , j

2 δ y uni, L − 1 2 2

 + δx unL +



i =1

L 1 −1 L 2 −1

2  1 1 − 2 , j +1

hxh y + α2

2  δ y uni+1, L − 1 h x h y 2 2



  i =0

j =1

2

j =0



μ¯ δ y uni, j− 1 2

2 

hxh y

D. Deng, C. Zhang / Applied Numerical Mathematics 62 (2012) 1864–1879

1877

Table 7 Numerical results for Example 3 with different parameters at T = 2 with t = h. h

Eu1

α = 0.1, γ = 1, c = 0.3

3.3 × 10−4 9.4 × 10−5 2.5 × 10−5 6.4 × 10−6

1/2 1/4 1/8 1/16

α = 0.1, 1/2 1/4 1/8 1/16

γ = 10, c = 0.3

α = 0.2, 1/2 1/4 1/8 1/16

γ = 10, c = 0.4

1.0 × 10−4 3.0 × 10−5 7.8 × 10−6 2.3 × 10−6 1.5 × 10−3 4.3 × 10−4 1.2 × 10−4 3.0 × 10−5

Table 8 Numerical results for Example 3 with

RE(I)

Lu1

Eu2

RE(II)

Lu2

∗ 1.80 1.92 1.94

1.5 × 10−3 4.7 × 10−4 1.6 × 10−4 3.3 × 10−5

4.7 × 10−5 3.0 × 10−6 1.9 × 10−7 1.2 × 10−8

∗ 4.00 3.98 3.91

1.8 × 10−4 1.1 × 10−5 6.7 × 10−7 4.6 × 10−8

∗ 1.80 1.92 1.74

4.9 × 10−5 1.5 × 10−4 4.0 × 10−5 1.0 × 10−5

1.5 × 10−5 9.4 × 10−7 5.9 × 10−8 3.9 × 10−9

∗ 4.00 3.98 3.91

5.7 × 10−5 3.4 × 10−6 2.1 × 10−7 1.4 × 10−8

∗ 1.81 1.86 1.96

6.2 × 10−3 1.8 × 10−3 4.8 × 10−4 1.2 × 10−4

7.7 × 10−4 5.4 × 10−5 3.2 × 10−6 2.0 × 10−7

∗ 3.85 4.05 4.00

2.7 × 10−3 1.6 × 10−4 9.6 × 10−6 6.0 × 10−7

α = 1, γ = 1 and c = 2 at different time levels, with t = 0.01 and h = 0.05.

T

1

2

3

4

5

6

Eu1 Lu1 Eu2 Lu2

3.3 × 10−5

8.3 × 10−5

2.6 × 10−4

9.1 × 10−4

3.3 × 10−3

1.2 × 10−2 8.8 × 10−2 3.2 × 10−6 2.0 × 10−5

1.9 × 10−4 3.2 × 10−8 1.5 × 10−7

5.6 × 10−4 6.1 × 10−8 2.7 × 10−7

1.7 × 10−3 8.4 × 10−8 4.2 × 10−7

6.1 × 10−3 1.3 × 10−7 9.2 × 10−7

2.3 × 10−2 5.1 × 10−7 3.5 × 10−6

Fig. 2. Example 3 with α = 0.1, γ = 1 and c = 0.3 (solved by algorithm II with h = t = 0.125): numerical solution at T = 7 (left column) and numerical solution at y = 0 (right column).

to approximate E (nt ). The notations E 1 and E 2 represent numerical energies obtained by Algorithm I and Algorithm II, respectively. Computational results are listed in Tables 7–8 and Figs. 2–3. Table 7 further testifies the convergence orders of Algorithm I and Algorithm II. Table 7 also indicates that Algorithm II gives very good results for a coarser grid. To describe behavior of approximation solution, we apply Algorithm II with t = h = 0.125 to solve Example 3 with α = 0.1, γ = 1 and c = 0.3 and obtain max1n56 un − Un ∞ = 7.8 × 10−7 . Graphs of approximate solution at T = 7 and at y = 0 are plotted in left and right columns of Fig. 2, respectively. With the help of Matlab 7.0, E (0) = 0.52037069969504 when α = 0.1, γ = 1 and c = 0.3; E (0) = 0.05203706996951 when α = 0.1, γ = 10 and c = 0.3; and E (0) = 0.224194312196306 when α = 0.2, γ = 10 and c = 0.4. For these cases, Algorithm II with t = h = 0.05 is used and yields numerical results as follows:





max un − Un ∞ = 2.3 × 10−8 ,

when α = 0.1,

γ = 1 and c = 0.3;

max un − Un ∞ = 7.3 × 10−9 ,

when α = 0.1,

γ = 10 and c = 0.3;

1n140



1n140



1878

D. Deng, C. Zhang / Applied Numerical Mathematics 62 (2012) 1864–1879

Fig. 3. Numerical energies from 0 to 7 for Example 3 with different parameters.





max un − Un ∞ = 1.0 × 10−7 ,

1n140

when α = 0.2,

γ = 10 and c = 0.4.

Meanwhile, their numerical energies at each time level are plotted in Fig. 3, from which we find they almost keep a constant as time increases. Finally, Algorithms I and II are used to solve Example 3 with α = γ = 1 and c = 2. Numerical results listed in Table 8, from which we find that Algorithm II can yields better solution. 7. Conclusions Motivated by the work of [4], we have designed a new ADI method for IBVP (1.1)–(1.3) in this article. By the discrete energy method, it has been shown that this method has convergence order of O (t 2 + h4x + h4y ) in H 1 - and L ∞ -norms. Moreover, a novel algorithm, which is a combination of the present ADI method with Richardson extrapolation (5.2), has been proposed. By introducing two auxiliary IBVPs, we have strictly proved that the algorithm can make final solution convergent with an order of O (t 4 + h4x + h4y ) in L ∞ -norm. Numerical results support theoretical analysis and exhibit computational efficiency of the algorithms. Extension of the new algorithm to (n + 1)-dimensional equation (1.1) is straightforward. Acknowledgements We are deeply indebted to Prof. Mingrong Cui whose work presented in [4] is very useful for us in this paper. We are also grateful to Principal Editor Martin Berzins and two referees for their valuable comments and suggestions, which have greatly improved the manuscript. References [1] K.C. Basak, P.C. Ray, R.K. Bera, Solution of non-linear Klein–Gordon equation with a quadratic non-linear term by Adomian decomposition method, Commun. Nonlinear Sci. Numer. Simul. 13 (2009) 718–723. [2] A.G. Bratsos, The solution of the two-dimensional sine-Gordon equation using the method of lines, J. Comput. Appl. Math. 206 (2007) 251–277. [3] M.R. Cui, Fourth-order compact scheme for the one-dimensional sine-Gordon equation, Numer. Methods Partial Differential Equations 25 (2009) 685– 711. [4] M.R. Cui, High order compact alternating direction implicit method for the generalized sine-Gordon equation, J. Comput. Appl. Math. 235 (2010) 837–849. [5] M. Dehghan, A. Mohebbi, High order implicit collocation method for the solution of two-dimensional linear hyperbolic equation, Numer. Methods Partial Differential Equations 25 (2009) 232–243. [6] M. Dehghan, A. Mohebbi, Z. Asgari, Fourth-order compact solution of the nonlinear Klein–Gordon equation, Numer. Algor. 52 (2009) 523–540. [7] D.W. Deng, Z.Y. Zhang, A new high-order algorithm for a class of nonlinear evolution equation, J. Phys. A: Math. Theor. 41 (2008) 015202. [8] K. Djidjeli, W.G. Price, E.H. Twizell, Numerical solutions of a damped sine-Gordon equation in two space variables, J. Engrg. Math. 29 (1995) 347–369. [9] J. Douglas Jr., J.E. Gunn, A general formulation of alternating direction methods. Part I. Parabolic and hyperbolic problems, Numer. Math. 6 (1964) 428–453. [10] F. Gao, C. Chi, Unconditionally stable difference schemes for a one-space-dimensional linear hyperbolic equation, Appl. Math. Comput. 187 (2007) 1272–1276. [11] Z. Gao, S. Xie, Fourth-order alternating direction implicit compact finite difference schemes for two-dimensional Schrödinger equations, Appl. Numer. Math. 61 (2011) 593–614. [12] D.R. Gulevich, F.V. Kusmartsev, S. Savel’ev, V.A. Yampol’skii, F. Nori, Shape and wobbling wave excitations in Josephson junctions: Exact solutions of the (2 + 1)-dimensional sine-Gordon model, Phys. Rev. B 80 (2009) 094509. [13] H. Hosseinzadeh, H. Jafari, M. Roohani, Application of Laplace decomposition method for solving Klein–Gordon equation, World Appl. Sci. J. 8 (2010) 809–813.

D. Deng, C. Zhang / Applied Numerical Mathematics 62 (2012) 1864–1879

1879

[14] K.J. in ’t Hout, B.D. Welfert, Unconditional stability of second-order ADI schemes applied to multi-dimensional diffusion equations with mixed derivative terms, Appl. Numer. Math. 59 (2009) 677–692. [15] J.C. Kalita, D.C. Dalal, A.K. Dass, A class of higher order compact schemes for the unsteady two-dimensional convection–diffusion equation with variable convection coefficients, Int. J. Numer. Meth. Fluids 38 (2002) 1111–1131. [16] S.A. Khuri, A. Sayfy, A spline collocation approach for the numerical solution of generalized nonlinear Klein–Gordon equation, Appl. Math. Comput. 216 (2010) 1047–1056. [17] G. Leibbrandt, New exact solutions of the classical sine-Gordon equation in 2 + 1 and 3 + 1 dimensions, Phys. Rev. Lett. 41 (1978) 435–438. [18] R.K. Lele, Compact finite difference schemes with spectral-like resolution, J. Comput. Phys. 103 (1992) 16–42. [19] H.L. Liao, Z.Z. Sun, Maximum norm error bounds of ADI and compact ADI methods for solving parabolic equations, Numer. Methods Partial Differential Equations 26 (2010) 37–60. [20] H.L. Liao, Z.Z. Sun, Maximum norm error estimates of efficient difference schemes for second-order wave equations, J. Comput. Appl. Math. 235 (2011) 2217–2233. [21] H. Lim, S. Kim, J. Douglas Jr., Numerical methods for viscous and nonviscous wave equations, Appl. Numer. Math. 57 (2007) 194–212. [22] A. Mayo, Higher-order accurate implicit finite difference method for evaluating American options, The European Journal of Finance 10 (2004) 212–237. [23] R.K. Mohanty, New unconditionally stable difference schemes for the solution of multi-dimensional telegraphic equations, Int. J. Comput. Math. 86 (2009) 2061–2071. [24] R.K. Mohanty, M.K. Jain, U. Arora, An unconditionally stable ADI method for the linear hyperbolic equation in three space dimensions, Int. J. Comput. Math. 79 (2002) 133–142. [25] J. Rashidinia, R. Mohammadi, Tension spline approach for the numerical solution of nonlinear Klein–Gordon equation, Comput. Phys. Comm. 181 (2010) 78–91. [26] J. Rashidinia, M. Ghasemi, R. Jalilian, Numerical solution of the nonlinear Klein–Gordon equation, J. Comput. Appl. Math. 233 (2010) 1866–1878. [27] A. Shah, L. Yuan, A. Khan, Upwind compact finite difference scheme for time-accurate solution of the incompressible Navier–Stokes equation, Appl. Math. Comput. 215 (2010) 3201–3213. [28] Q. Sheng, A.Q.M. Khaliq, D.A. Voss, Numerical simulation of two-dimensional sine-Gordon solitons via a split cosine scheme, Math. Comput. Simulation 68 (2005) 355–373. [29] Z.Y. Zhang, D.W. Deng, A new alternating-direction finite element method for hyperbolic equation, Numer. Methods Partial Differential Equations 23 (2007) 1530–1559. [30] C.X. Zheng, Numerical solution to the sine-Gordon equation defined on the whole real axis, SIAM J. Sci. Comput. 29 (2007) 2494–2506.