High Efficient Numerical Methods for Viscous and Nonviscous Wave Problems

Abstract

This paper is concerned with accurate and efficient numerical methods for solving viscous and nonviscous wave problems. The paper first introduces a new second-order PR-ADI like scheme. For an efficient simulation, the scheme is also extended to a high-order compact PRADI like method. 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 results are presented to show the accuracy and efficiency of the new algorithms.

1. Introduction

In this paper, we consider the following differential problem with initial and boundary conditions on domain Ω=[0,1]²: 𝛾₁𝜕𝑢𝜕𝑡+𝛾₂𝜕²𝑢𝜕𝑡²𝜕−𝑎▵𝑢−𝑞,(],𝑢𝜕𝑡▵𝑢=𝑆𝑥,𝑦,𝑡𝑥,𝑦∈Ω,𝑡∈0,𝑇𝑥,𝑦,0=𝑢₀𝑥,𝑦,𝜕𝑢,𝑢,(],𝜕𝑡𝑥,𝑦,0=𝜓𝑥,𝑦𝑥,𝑦∈Ω,𝑥,𝑦,𝑡=𝜑𝑥,𝑦,𝑡𝑥,𝑦∈𝜕Ω,𝑡∈0,𝑇 (1)(2)(3) where 𝛾₁,𝛾₂,𝑎, and 𝑞 are nonnegative constants, and 𝑎𝛾₂>𝛾₁𝑞, 𝑆 denotes the source, and 𝑢₀,𝜓, and 𝜑 are given functions, and 𝜕Ω is the boundary of the domain Ω.

Equations (1)–(3) govern various physical phenomena by choosing the coefficients. When 𝛾₂=0 and 𝑞=0, the equation represents standard heat transfer, and when 𝛾₁=0 and 𝑞=0, the equation represents standard wave equations. Since there are many efficient algorithms known for them, we exclude their consideration.

Some conventional numerical approaches for solving wave equations introduce auxiliary variables to rewrite the equations as first-order hyperbolic systems [1–3]; however, these approaches introduce new unknowns which results in an increase in the number of variables in the discrete problems. Thus, there are advantages in keeping the formulation (1)–(3) involving the second time-derivative and a scalar unknown. It has been the case that it is hard to construct methods combining good stability with high accuracy with this formulation [4].

This paper is concerned with numerical solutions to viscous and nonviscous nonhomogeneous wave problems. The wave equation is often solved by explicit time-stepping schemes, which require to choose a time step size sufficiently small to satisfy the stability condition and to reduce numerical dispersion as well. As is well known, the alternating direction implicit (ADI) schemes [3–14] are unconditionally stable and only need to solve a sequence of tridiagonal linear systems. In 1955, Peaceman and Rachford Jr. [15] presented a method to solve two-dimensional parabolic equations. When constructing the difference scheme, they first evaluated the derivative with respect to 𝑥 implicitly, the derivative with respect to 𝑦 explicitly and used a time step of Δ𝑡/2 to get a scheme, then took another half step with explicit in 𝑥 and implicit in 𝑦 to get another scheme. The final two-step scheme, called Peaceman and Rachford (PR) scheme, is unconditionally stable and second-order accurate. Alternative conventional approaches for solving wave equations introduce an auxiliary variable to rewrite the equation as a first-order hyperbolic system. With these approaches one introduces new unknowns, which result in an increase in the number of variables in the discrete problems. Thus, there are good reasons to try to keep the formulation involving the second time-derivative and a scalar unknown. Lim et al. [11] introduce a stable three-level locally one-dimensional (LOD) method for solving (1)–(3), which is second-order in space.

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 [8–10, 12, 14]. To obtain satisfactory higher order numerical results with reasonable computational cost, there have been attempts to develop higher order compact ADI methods.

In this paper, we propose a new set of ADI methods for viscous and nonviscous wave problems. In Section 3, we obtain a new PR ADI-like scheme with truncation error 𝑂(Δ𝑡²+ℎ²) and analyze the error estimate by discrete energy method. In Section 4, we further extend the method to compact PR ADI-like scheme with truncation error 𝑂(Δ𝑡²+ℎ⁴). In Section 5, 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 that they have much higher efficiency.

2. Notation

The domain Ω is divided into a mesh by points 𝑥_𝑖=𝑖ℎ,𝑦_𝑗=𝑗ℎ(𝑖,𝑗=0,…,𝑁), denoted by Ω_ℎ, where ℎ=1/𝑁 is the spatial mesh size in both 𝑥 and 𝑦 directions. Let 𝑡^𝑛=𝑛Δ𝑡,𝑡^𝑛+1/2=(𝑛+1/2)Δ𝑡, with Δ𝑡 being the time increment. Denote the values of 𝑢(𝑥_𝑖,𝑦_𝑗,𝑡^𝑛) and 𝑆(𝑥_𝑖,𝑦_𝑗,𝑡^𝑛+1/2) by 𝑢^𝑛_𝑖,𝑗 and 𝑆^𝑛+1/2_𝑖,𝑗, respectively. Denote also the difference solutions to 𝑢^𝑛_𝑖,𝑗 and 𝑣^𝑛_𝑖,𝑗 by 𝑈^𝑛_𝑖,𝑗 and 𝑉^𝑛_𝑖,𝑗, respectively.

Assume that 𝑉_ℎ={𝑤∣𝑤={𝑤_𝑖,𝑗}∈Ω_ℎ, and 𝑤|_𝜕Ω=0}. For 𝑤={𝑤_𝑖,𝑗}∈𝑉_ℎ and 𝑣={𝑣_𝑖,𝑗}∈𝑉_ℎ, denote by 𝛿_𝑡𝑣^𝑛_𝑖,𝑗=𝑣^𝑛+1_𝑖,𝑗−𝑣^𝑛_𝑖,𝑗Δ𝑡,𝛿_𝑥𝑣^𝑛_𝑖,𝑗=𝑣^𝑛_𝑖+1,𝑗−𝑣^𝑛_𝑖,𝑗ℎ,𝛿_𝑦𝑣^𝑛_𝑖,𝑗=𝑣^𝑛_𝑖,𝑗+1−𝑣^𝑛_𝑖,𝑗ℎ,𝛿²_𝑥𝑣^𝑛_𝑖,𝑗=𝑣^𝑛_𝑖+1,𝑗−2𝑣^𝑛_𝑖,𝑗+𝑣^𝑛_{𝑖−1,𝑗}ℎ²,𝛿²_𝑦𝑣^𝑛_𝑖,𝑗=𝑣^𝑛_𝑖,𝑗+1−2𝑣^𝑛_𝑖,𝑗+𝑣^𝑛_{𝑖,𝑗−1}ℎ². (4) For ∀𝑣,𝑤∈𝑉_ℎ, define their discrete inner product and norms as follows: (𝑣,𝑤)=^𝑁−1_𝑖,𝑗=1𝑣_𝑖,𝑗𝑤_𝑖,𝑗ℎ²,‖𝑣‖=(),𝛿𝑣,𝑣_𝑥𝑣=^𝑁−1_𝑖=0^𝑁−1_𝑗=1𝛿_𝑥𝑣_𝑖,𝑗²ℎ²,𝛿_𝑦𝑣=^𝑁−1_𝑖=1^𝑁−1_𝑗=0𝛿_𝑦𝑣_𝑖,𝑗²ℎ²,𝛿_𝑥𝛿_𝑦𝑣=^𝑁−1_𝑖=0^𝑁−1_𝑗=0𝛿_𝑥𝛿_𝑦𝑣_𝑖,𝑗²ℎ²,|𝑣|_𝐻¹=𝛿_𝑥𝑣²+𝛿_𝑦𝑣². (5)

3. 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 differential problems (1)–(3).

3.1. Construction of the New ADI Scheme

Introducing 𝑣=𝛾₁𝑢+𝛾₂(𝜕𝑢/𝜕𝑡), we can rewrite (1) as follows: 𝜕𝑣(],𝜕𝑡−𝛼▵𝑣−𝛽▵𝑢=𝑆,𝑥,𝑦∈Ω,𝑡∈0,𝑇𝑣=𝛾₁𝑢+𝛾₂𝜕𝑢,𝑢𝜕𝑡𝑥,𝑦,0=𝑢₀,𝑣𝑥,𝑦𝑥,𝑦,0=𝛾₁𝑢₀𝑥,𝑦+𝛾₂𝜓,𝑢,𝑣𝑥,𝑦𝑥,𝑦∈Ω,𝑥,𝑦,𝑡=𝜑𝑥,𝑦,𝑡𝑥,𝑦,𝑡=𝛾₁𝜑𝑥,𝑦,𝑡+𝛾₂𝜑_𝑡,(],𝑥,𝑦,𝑡𝑥,𝑦∈𝜕Ω,𝑡∈0,𝑇 (6)(7)(8)(9) where 𝛼=𝑞/𝛾₂,𝛽=𝑎−(𝛾₁𝑞/𝛾₂).

Applying Crank-Nicolson implicit discretization to (6), we have 𝑉^𝑛+1_𝑖,𝑗−𝑉^𝑛_𝑖,𝑗𝛿Δ𝑡−𝛼²_𝑥+𝛿²_𝑦𝑉^𝑛+1_𝑖,𝑗+𝑉^𝑛_𝑖,𝑗2−𝛽𝛿²_𝑥+𝛿²_𝑦𝑈^𝑛+1_𝑖,𝑗+𝑈^𝑛_𝑖,𝑗2=𝑆^𝑛+1/2_𝑖,𝑗,𝑉^𝑛+1_𝑖,𝑗+𝑉^𝑛_𝑖,𝑗2=𝛾₁𝑈^𝑛+1_𝑖,𝑗+𝑈^𝑛_𝑖,𝑗2+𝛾₂𝑈^𝑛+1_𝑖,𝑗−𝑈^𝑛_𝑖,𝑗.Δ𝑡 (10)

From (10), we obtain 𝑉^𝑛+1_𝑖,𝑗−Δ𝑡2𝛼𝛿²_𝑥+𝛿²_𝑦𝑉^𝑛+1_𝑖,𝑗−Δ𝑡2𝛽𝛿²_𝑥+𝛿²_𝑦𝑈^𝑛+1_𝑖,𝑗=𝑉^𝑛_𝑖,𝑗+Δ𝑡2𝛼𝛿²_𝑥+𝛿²_𝑦𝑉^𝑛_𝑖,𝑗+Δ𝑡2𝛽𝛿²_𝑥+𝛿²_𝑦𝑈^𝑛_𝑖,𝑗+Δ𝑡𝑆^𝑛+1/2_𝑖,𝑗,𝑈^𝑛+1_𝑖,𝑗=Δ𝑡𝛾₁Δ𝑡+2𝛾₂𝑉^𝑛+1_𝑖,𝑗+𝑉^𝑛_𝑖,𝑗−𝛾₁Δ𝑡−2𝛾₂𝛾₁Δ𝑡+2𝛾₂𝑈^𝑛_𝑖,𝑗. (11)(12) Let 𝜆₁=Δ𝑡𝛾₁Δ𝑡+2𝛾₂,𝜆₂=𝛾₁Δ𝑡−2𝛾₂𝛾₁Δ𝑡+2𝛾₂,𝑐=𝛼+𝛽𝜆₁1,𝛾=2𝛽1−𝜆₂, (13) then (11) can be reorganized as 𝑉^𝑛+1_𝑖,𝑗−𝑐Δ𝑡2𝛿²_𝑥+𝛿²_𝑦𝑉^𝑛+1_𝑖,𝑗=𝑉^𝑛_𝑖,𝑗+𝑐Δ𝑡2𝛿²_𝑥+𝛿²_𝑦𝑉^𝑛_𝑖,𝑗+𝛾Δ𝑡𝛿²_𝑥+𝛿²_𝑦𝑈^𝑛_𝑖,𝑗+Δ𝑡𝑆^𝑛+1/2_𝑖,𝑗. (14) Adding perturbing term (𝑐²Δ𝑡²/4)𝛿²_𝑥𝛿²_𝑦(𝑉^𝑛+1_𝑖,𝑗−𝑉^𝑛_𝑖,𝑗) to the left hand side of (14), we factor (14) as 1−𝑐Δ𝑡2𝛿²_𝑥1−𝑐Δ𝑡2𝛿²_𝑦𝑉^𝑛+1_𝑖,𝑗=1+𝑐Δ𝑡2𝛿²_𝑥1+𝑐Δ𝑡2𝛿²_𝑦𝑉^𝑛_𝑖,𝑗+𝛾Δ𝑡𝛿²_𝑥+𝛿²_𝑦𝑈^𝑛_𝑖,𝑗+Δ𝑡𝑆^𝑛+1/2_𝑖,𝑗. (15) Introducing the intermediate variable 𝑉^𝑛+1/2_𝑖,𝑗, we obtain the PR ADI-like scheme as follows: 1−𝑐Δ𝑡2𝛿²_𝑥𝑉^𝑛+1/2_𝑖,𝑗=1+𝑐Δ𝑡2𝛿²_𝑦𝑉^𝑛_𝑖,𝑗+𝛾Δ𝑡2𝛿²_𝑦+1𝑐𝑈^𝑛_𝑖,𝑗+Δ𝑡2𝑆^𝑛+1/2_𝑖,𝑗,1−𝑐Δ𝑡2𝛿²_𝑦𝑉^𝑛+1_𝑖,𝑗=1+𝑐Δ𝑡2𝛿²_𝑥𝑉^𝑛+1/2_𝑖,𝑗+𝛾Δ𝑡2𝛿²_𝑦−1𝑐𝑈^𝑛_𝑖,𝑗+Δ𝑡2𝑆^𝑛+1/2_𝑖,𝑗,𝑈^𝑛+1_𝑖,𝑗=𝜆₁𝑉^𝑛+1_𝑖,𝑗+𝑉^𝑛_𝑖,𝑗−𝜆₂𝑈^𝑛_𝑖,𝑗. (16)(17)(18) From (16), (17), we get the boundary equation as follows: 𝑉^𝑛+1/2_𝑖,𝑗=𝑉^𝑛+1_𝑖,𝑗+𝑉^𝑛_𝑖,𝑗2+𝛾𝑐𝑈^𝑛_𝑖,𝑗−𝑐Δ𝑡4𝛿²_𝑦𝑉^𝑛+1_𝑖,𝑗−𝑉^𝑛_𝑖,𝑗. (19) But usually we compute the boundary values of the intermediate variable using the following simple equality: 𝑉^𝑛+1/2_𝑖,𝑗=𝑉^𝑛+1_𝑖,𝑗+𝑉^𝑛_𝑖,𝑗2+𝛾𝑐𝑈^𝑛_𝑖,𝑗, (20) where 𝑖=0,𝑁,𝑗=0,…,𝑁−1.

We see that schemes (16)–(20) can solve differential problems (1)–(3). In addition, following the idea of Douglas ([2, 3]), a Douglas-like scheme is as follows: 1−𝑐Δ𝑡2𝛿²_𝑥𝑉^𝑛+1/2_𝑖,𝑗−𝑉^𝑛_𝑖,𝑗=𝑐Δ𝑡𝛿²_𝑥+𝛿²_𝑦𝑉^𝑛_𝑖,𝑗+𝛾Δ𝑡𝛿²_𝑥+𝛿²_𝑦𝑈^𝑛_𝑖,𝑗+Δ𝑡𝑆^𝑛+1/2_𝑖,𝑗,1−𝑐Δ𝑡2𝛿²_𝑦𝑉^𝑛+1_𝑖,𝑗−𝑉^𝑛_𝑖,𝑗=𝑉^𝑛+1/2_𝑖,𝑗−𝑉^𝑛_𝑖,𝑗. (21)(22) The intermediate value of 𝑉^𝑛+1/2_𝑖,𝑗 on the boundary is easy to be obtained from (22). PR ADI-like scheme (16)–(18) and Douglas-like scheme (21)-(22) are all derived from equation (15), but PR ADI-like scheme has more efficiency, which will be found in numerical experiment.

3.2. Error Estimate

We have derived a new kind of ADI scheme in Section 3.1. We further analyze the error of the scheme starting from (15). Expanding (15), we get 𝑉^𝑛+1_𝑖,𝑗−𝑉^𝑛_𝑖,𝑗𝛿Δ𝑡−𝛼²_𝑥+𝛿²_𝑦𝑉^𝑛+1_𝑖,𝑗+𝑉^𝑛_𝑖,𝑗2−𝛽𝛿²_𝑥+𝛿²_𝑦𝑈^𝑛+1_𝑖,𝑗+𝑈^𝑛_𝑖,𝑗2+𝑐²Δ𝑡²4𝛿²_𝑥𝛿²_𝑦𝑉^𝑛+1_𝑖,𝑗−𝑉^𝑛_𝑖,𝑗Δ𝑡=𝑆^𝑛+1/2_𝑖,𝑗. (23) Discretizing (6) as above and letting 𝜉=𝑢−𝑈,𝜂=𝑣−𝑉, we obtain the error equations as follows: 𝜂^𝑛+1_𝑖,𝑗−𝜂^𝑛_𝑖,𝑗𝛿Δ𝑡−𝛼²_𝑥+𝛿²_𝑦𝜂^𝑛+1_𝑖,𝑗+𝜂^𝑛_𝑖,𝑗2−𝛽𝛿²_𝑥+𝛿²_𝑦𝜉^𝑛+1_𝑖,𝑗+𝜉^𝑛_𝑖,𝑗2+𝑐²Δ𝑡²4𝛿²_𝑥𝛿²_𝑦𝜂^𝑛+1_𝑖,𝑗−𝜂^𝑛_𝑖,𝑗Δ𝑡=𝑅^𝑛+1/2_𝑖,𝑗,𝜂^𝑛+1_𝑖,𝑗+𝜂^𝑛_𝑖,𝑗2=𝛾₁𝜉^𝑛+1_𝑖,𝑗+𝜉^𝑛_𝑖,𝑗2+𝛾₂𝜉^𝑛+1_𝑖,𝑗−𝜉^𝑛_𝑖,𝑗+Δ𝑡^∼𝑅^𝑛+1/2_𝑖,𝑗, (24)(25) where the truncation errors 𝑅^𝑛+1/2_𝑖,𝑗=Δ𝑡²1𝜕24³𝑣𝜕𝑡³−𝛼8𝜕⁴𝑣𝜕𝑡²𝜕𝑥²+𝜕⁴𝑣𝜕𝑡²𝜕𝑦²−𝛽8𝜕⁴𝑢𝜕𝑡²𝜕𝑥²+𝜕⁴𝑢𝜕𝑡²𝜕𝑦²+𝑐²Δ𝑡²4𝜕⁵𝑣𝜕𝑡𝜕𝑥²𝜕𝑦²𝑥_𝑖,𝑦_𝑗,𝑡^𝑛+1/2−𝛼ℎ²𝜕12⁴𝑣𝜕𝑥⁴+𝜕⁴𝑣𝜕𝑦⁴𝑥_𝑖,𝑦_𝑗,𝑡^𝑛+1/2−𝛽ℎ²𝜕12⁴𝑢𝜕𝑥⁴+𝜕⁴𝑢𝜕𝑦⁴𝑥_𝑖,𝑦_𝑗,𝑡^𝑛+1/2+𝑂Δ𝑡⁴+Δ𝑡²ℎ²+ℎ⁴,_∼𝑅^𝑛+1/2_𝑖,𝑗=Δ𝑡²18𝜕²𝑣𝜕𝑡²−𝛾₁8𝜕²𝑢𝜕𝑡²−𝛾₂𝜕24³𝑢𝜕𝑡³𝑥_𝑖,𝑦_𝑗,𝑡^𝑛+1/2+𝑂Δ𝑡⁴. (26) Hence there exist positive constants 𝐶₁,𝐶₂ such that 𝑅^𝑛+1/2_𝑖,𝑗≤𝐶₁Δ𝑡²+ℎ²,𝑅^∼𝑛+1/2_𝑖,𝑗≤𝐶₂Δ𝑡². (27)

Lemma 1 (see [1]). For 𝑤∈𝑉_ℎ, the following equalities hold: −^𝑁−1_𝑖=1^𝑁−1_𝑗=1𝑤_𝑖,𝑗𝛿²_𝑥𝑤_𝑖,𝑗ℎ²=^𝑁−1_𝑖=0^𝑁−1_𝑗=1𝛿_𝑥𝑤_𝑖,𝑗²ℎ²=𝛿_𝑥𝑤²,−^𝑁−1_𝑖=1^𝑁−1_𝑗=1𝑤_𝑖,𝑗𝛿²_𝑦𝑤_𝑖,𝑗ℎ²=^𝑁−1_𝑖=1^𝑁−1_𝑗=0𝛿_𝑦𝑤_𝑖,𝑗²ℎ²=𝛿_𝑦𝑤². (28)

From Lemma 1 and Cauchy-Schwarz inequality, multiplying ((𝜂^𝑛+1_𝑖,𝑗+𝜂^𝑛_𝑖,𝑗)/2)ℎ² to both sides of (24), computing the inner product, it is easy to follow that 𝐼₁=𝜂^𝑛+1−𝜂^𝑛,𝜂Δ𝑡^𝑛+1+𝜂^𝑛2=1𝜂2Δ𝑡^𝑛+12−𝜂^𝑛2,𝐼₂₁=−𝛼𝛿²_𝑥𝜂^𝑛+1+𝜂^𝑛2,𝜂^𝑛+1+𝜂^𝑛2=𝛼𝛿_𝑥𝜂^𝑛+1+𝜂^𝑛2²,𝐼₂₂=−𝛼𝛿²_𝑦𝜂^𝑛+1+𝜂^𝑛2,𝜂^𝑛+1+𝜂^𝑛2=𝛼𝛿_𝑦𝜂^𝑛+1+𝜂^𝑛2²,𝐼₃₁=−𝛽𝛿²_𝑥𝜉^𝑛+1+𝜉^𝑛2,𝜂^𝑛+1+𝜂^𝑛2=−𝛽𝛿²_𝑥𝜉^𝑛+1+𝜉^𝑛2,𝜂^𝑛+1+𝜂^𝑛2=−𝛽𝛿²_𝑥𝜉^𝑛+1+𝜉^𝑛2,𝛾₁𝜉^𝑛+1+𝜉^𝑛2+𝛾₂𝜉^𝑛+1−𝜉^𝑛+Δ𝑡^∼𝑅^𝑛+1/2=𝛽𝛾₁𝛿_𝑥𝜉^𝑛+1+𝜉^𝑛2²+𝛽𝛾₂𝛿2Δ𝑡^𝑥𝜉^𝑛+12−𝛿_𝑥𝜉^𝑛2−−𝛽𝛿_𝑥𝜉^𝑛+1+𝜉^𝑛2,𝛿𝑅_𝑥^∼𝑛+1/2≥𝛽𝛾₂𝛿2Δ𝑡^𝑥𝜉^𝑛+12−𝛿_𝑥𝜉^𝑛2−𝛽𝛿_𝑥𝜉^𝑛+1+𝜉^𝑛2𝛿𝑅_𝑥^∼𝑛+1/2≥𝛽𝛾₂𝛿2Δ𝑡^𝑥𝜉^𝑛+12−𝛿_𝑥𝜉^𝑛2−𝛽42𝛿_𝑥𝜉^𝑛+12+2𝛿_𝑥𝜉^𝑛2+𝛿𝑅_𝑥^∼𝑛+1/22,𝐼₃₂≥𝛽𝛾₂𝛿2Δ𝑡^𝑦𝜉^𝑛+12−𝛿_𝑦𝜉^𝑛2−𝛽42𝛿_𝑦𝜉^𝑛+12+2𝛿_𝑦𝜉^𝑛2+𝛿𝑅_𝑦^∼𝑛+1/22,𝐼₄=𝑐²Δ𝑡²4𝛿²_𝑥𝛿²_𝑦𝜂^𝑛+1−𝜂^𝑛,𝜂Δ𝑡^𝑛+1+𝜂^𝑛2=12Δ𝑡𝑐Δ𝑡²4𝛿_𝑥𝛿_𝑦𝜂^𝑛+12−𝛿_𝑥𝛿_𝑦𝜂^𝑛2,𝐼₅=𝑅^𝑛+1/2,𝜂^𝑛+1+𝜂^𝑛2≤14𝑅^𝑛+1/22+12𝜂^𝑛+12+𝜂^𝑛2. (29) Multiplying 2Δ𝑡 to both sides of the result derived from (29), we have 𝜂^𝑛+12−𝜂^𝑛2+2𝛼Δ𝑡𝛿_𝑥𝜂^𝑛+1+𝜂^𝑛2+𝛿_𝑦𝜂^𝑛+1+𝜂^𝑛2+𝛽𝛾₂𝛿_𝑥𝜉^𝑛+12+𝛿_𝑦𝜉^𝑛+12−𝛿_𝑥𝜉^𝑛2−𝛿_𝑦𝜉^𝑛2+𝑐Δ𝑡²4𝛿_𝑥𝛿_𝑦𝜂^𝑛+12−𝛿_𝑥𝛿_𝑦𝜂^𝑛2≤Δ𝑡2𝛽𝛿𝑅_𝑥^∼𝑛+1/22+𝛽𝛿𝑅_𝑦^∼𝑛+1/22+𝑅^𝑛+1/22+Δ𝑡𝜂^𝑛+12+𝜂^𝑛2+𝛽𝛿_𝑥𝜉^𝑛+12+𝛿_𝑥𝜉^𝑛2+𝛿_𝑦𝜉^𝑛+12+𝛿_𝑦𝜉^𝑛2. (30) Summing for 𝑛 to both sides of (30), we obtain 𝜂^𝑛2+𝛽𝛾₂𝛿_𝑥𝜉^𝑛2+𝛿_𝑦𝜉^𝑛2≤Δ𝑡2^𝑛−1_𝑙=1𝑅^𝑙+1/22+𝛽𝛿𝑅_𝑥^∼𝑙+1/22+𝛽𝛿𝑅_𝑦^∼𝑙+1/22+Δ𝑡^𝑛−1_𝑙=12𝜂^𝑙2+𝛽𝛿_𝑥𝜉^𝑙2+𝛽𝛿_𝑦𝜉^𝑙2, (31) noticing that 𝜉⁰=𝜂⁰=0, by Gronwall Lemma, we obtain the following theorem.

Theorem 2. Assume that 𝑢,𝑣 are the accurate solutions of (6)–(9) with sufficient smoothness and 𝑈,𝑉 are the difference solutions of (16)–(18). Let 𝜉=𝑢−𝑈,𝜂=𝑣−𝑉, then there exists a positive constant 𝐾 independent of Δ𝑡 and ℎ, such that max_{[]𝑛≤𝑇/Δ𝑡}𝜂^𝑛+𝛿_𝑥𝜉^𝑛+𝛿_𝑦𝜉^𝑛≤𝐾Δ𝑡²+ℎ². (32)

4. Generalized Compact ADI Scheme

In this section, we extend the method to compact ADI scheme for the numerical solution of the wave problems (1)–(3).

4.1. Construction of the Compact ADI Scheme

Introducing 𝑓_𝑥𝑣=𝜕²𝑣𝜕𝑥²,𝑓_𝑦𝑣=𝜕²𝑣𝜕𝑦²,𝑓_𝑥𝑢=𝜕²𝑢𝜕𝑥²,𝑓_𝑦𝑢=𝜕²𝑢𝜕𝑦². (33) We can rewrite (6) as follows: 𝜕𝑣𝑓𝜕𝑡−𝛼_𝑥𝑣+𝑓_𝑦𝑣−𝛽𝑓_𝑥𝑢+𝑓_𝑦𝑢=𝑆𝑥,𝑦,𝑡. (34)

Using the fourth-order compact finite difference discretization [6, 8], we discretize (33) as follows: 𝐴_𝑥𝑓_𝑥^𝑛_𝑖,𝑗1∶=𝑓12^{𝑥𝑛𝑖−1,𝑗}+56𝑓_𝑥^𝑛_𝑖,𝑗+1𝑓12^{𝑥𝑛𝑖+1,𝑗}=𝛿²_𝑥𝑈^𝑛_𝑖,𝑗,𝐴_𝑦𝑓_𝑦^𝑛_𝑖,𝑗1∶=𝑓12^{𝑦𝑛𝑖,𝑗−1}+56𝑓_𝑦^𝑛_𝑖,𝑗+1𝑓12^{𝑦𝑛𝑖,𝑗+1}=𝛿²_𝑦𝑈^𝑛_𝑖,𝑗. (35)

The Crank-Nicolson implicit discretization of (34) reads 𝑉^𝑛+1_𝑖,𝑗−𝑉^𝑛_𝑖,𝑗−𝛼𝑓▵𝑡^{𝑥𝑣𝑛+1𝑖,𝑗}+𝑓_𝑥𝑣^𝑛_𝑖,𝑗+𝑓_𝑦𝑣^𝑛+1_𝑖,𝑗+𝑓_𝑦𝑣^𝑛_𝑖,𝑗2−𝛽𝑓_𝑥𝑢^𝑛+1_𝑖,𝑗+𝑓_𝑥𝑢^𝑛_𝑖,𝑗+𝑓_𝑦𝑢^𝑛+1_𝑖,𝑗+𝑓_𝑦𝑢^𝑛_𝑖,𝑗2=𝑔^𝑛+1/2_𝑖,𝑗, (36) where 𝑖,𝑗=1,2,…,𝑁−1. Multiplying (36) by Δ𝑡, applying to its both sides with 𝐴_𝑥𝐴_𝑦, and using the fact that the two operators 𝐴_𝑥 and 𝐴_𝑦 commute with each other, we obtain the implicit approximation of (6)–(9) 𝐴_𝑥𝐴_𝑦𝑉^𝑛+1_𝑖,𝑗−𝑉^𝑛_𝑖,𝑗−𝛼Δ𝑡2𝐴_𝑦𝐴_𝑥𝑓_𝑥𝑣^𝑛+1_𝑖,𝑗+𝐴_𝑦𝐴_𝑥𝑓_𝑥𝑣^𝑛_𝑖,𝑗+𝐴_𝑥𝐴_𝑦𝑓_𝑦𝑣^𝑛+1_𝑖,𝑗+𝐴_𝑥𝐴_𝑦𝑓_𝑦𝑣^𝑛_𝑖,𝑗−𝛽Δ𝑡2𝐴_𝑦𝐴_𝑥𝑓_𝑥𝑢^𝑛+1_𝑖,𝑗+𝐴_𝑦𝐴_𝑥𝑓_𝑥𝑢^𝑛_𝑖,𝑗+𝐴_𝑥𝐴_𝑦𝑓_𝑦𝑢^𝑛+1_𝑖,𝑗+𝐴_𝑥𝐴_𝑦𝑓_𝑦𝑢^𝑛_𝑖,𝑗=Δ𝑡𝐴_𝑥𝐴_𝑦𝑆^𝑛+1/2_𝑖,𝑗,𝑉^𝑛+1_𝑖,𝑗+𝑉^𝑛_𝑖,𝑗2=𝛾₁𝑈^𝑛+1_𝑖,𝑗+𝑈^𝑛_𝑖,𝑗2+𝛾₂𝑈^𝑛+1_𝑖,𝑗−𝑈^𝑛_𝑖,𝑗.Δ𝑡 (37)(38) Let 𝜆₁=Δ𝑡𝛾₁Δ𝑡+2𝛾₂,𝜆₂=𝛾₁Δ𝑡−2𝛾₂𝛾₁Δ𝑡+2𝛾₂,𝑐=𝛼+𝛽𝜆₁1,𝛾=2𝛽1+𝜆₂, (39) then (37) can be conveniently written as 𝐴_𝑥𝐴_𝑦−𝑐Δ𝑡2𝐴_𝑦𝛿²_𝑥−𝑐Δ𝑡2𝐴_𝑥𝛿²_𝑦𝑉^𝑛+1_𝑖,𝑗=𝐴_𝑥𝐴_𝑦+𝑐Δ𝑡2𝐴_𝑦𝛿²_𝑥+𝑐Δ𝑡2𝐴_𝑥𝛿²_𝑦𝑉^𝑛_𝑖,𝑗+𝛾Δ𝑡𝐴_𝑦𝛿²_𝑥+𝐴_𝑥𝛿²_𝑦𝑈^𝑛_𝑖,𝑗+Δ𝑡𝐴_𝑥𝐴_𝑦𝑆^𝑛+1/2_𝑖,𝑗. (40) Adding perturbing term (𝑐²Δ𝑡²/4)𝛿²_𝑥𝛿²_𝑦(𝑉^𝑛+1_𝑖,𝑗−𝑉^𝑛_𝑖,𝑗) to the left hand side of (40) and introducing the intermediate variable 𝑉^𝑛+1/2_𝑖,𝑗, we obtain the compact PR ADI-like scheme as follows: 𝐴_𝑥−𝑐Δ𝑡2𝛿²_𝑥𝑉^𝑛+1/2_𝑖,𝑗=𝐴_𝑦+𝑐Δ𝑡2𝛿²_𝑦𝑉^𝑛_𝑖,𝑗+𝛾Δ𝑡2𝛿²_𝑦+1𝑐𝐴_𝑦𝑈^𝑛_𝑖,𝑗+Δ𝑡2𝐴_𝑦𝑆^𝑛+1/2_𝑖,𝑗,𝐴_𝑦−𝑐Δ𝑡2𝛿²_𝑦𝑉^𝑛+1_𝑖,𝑗=𝐴_𝑥+𝑐Δ𝑡2𝛿²_𝑥𝑉^𝑛+1/2_𝑖,𝑗+𝛾Δ𝑡2𝛿²_𝑦−1𝑐𝐴_𝑦𝑈^𝑛_𝑖,𝑗+Δ𝑡2𝐴_𝑦𝑆^𝑛+1/2_𝑖,𝑗,𝑈^𝑛+1_𝑖,𝑗=𝜆₁𝑉^𝑛+1_𝑖,𝑗+𝑉^𝑛_𝑖,𝑗−𝜆₂𝑈^𝑛_𝑖,𝑗. (41)(42)(43) From (41) and (42), we obtain the boundary equation as follows: 𝐴_𝑥𝑉^𝑛+1/2_𝑖,𝑗=𝐴_𝑦𝑉^𝑛+1_𝑖,𝑗+𝑉^𝑛_𝑖,𝑗2+𝛾𝑐𝐴_𝑦𝑈^𝑛_𝑖,𝑗−𝑐Δ𝑡4𝛿²_𝑦𝑉^𝑛+1_𝑖,𝑗−𝑉^𝑛_𝑖,𝑗. (44) But usually we compute the boundary values of the intermediate using the following simple equality: 𝑉^𝑛+1/2_𝑖,𝑗=𝑉^𝑛+1_𝑖,𝑗+𝑉^𝑛_𝑖,𝑗2+𝛾𝑐𝑈^𝑛_𝑖,𝑗, (45) where 𝑖=0,𝑁,𝑗=1,𝑁−1. In addition, following the idea of Douglas [2, 3], we can get a compact Douglas scheme as follows: 𝐴_𝑥−𝑐Δ𝑡2𝛿²_𝑥𝑉^𝑛+1/2_𝑖,𝑗−𝑉^𝑛_𝑖,𝑗=𝑐Δ𝑡𝐴_𝑦𝛿²_𝑥+𝐴_𝑥𝛿²_𝑦𝑉^𝑛_𝑖,𝑗+𝛾Δ𝑡𝐴_𝑦𝛿²_𝑥+𝐴_𝑥𝛿²_𝑦𝑈^𝑛_𝑖,𝑗+Δ𝑡𝐴_𝑥𝐴_𝑦𝑔^𝑛+1/2_𝑖,𝑗,𝐴_𝑦−𝑐Δ𝑡2𝛿²_𝑦𝑉^𝑛+1_𝑖,𝑗−𝑉^𝑛_𝑖,𝑗=𝑉^𝑛+1/2_𝑖,𝑗−𝑉^𝑛_𝑖,𝑗. (46)(47) The intermediate value of 𝑉^𝑛+1/2_𝑖,𝑗 on the boundary is easy to be obtained from (47).

4.2. Error Estimate

We have derived a kind of compact PR ADI-like scheme in Section 4.1. In the following, we further analyze the error of the scheme. Expanding (41)-(42), we get 𝐴_𝑥𝐴_𝑦𝑉^𝑛+1_𝑖,𝑗−𝑉^𝑛_𝑖,𝑗𝐴Δ𝑡−𝛼_𝑦𝛿²_𝑥+𝐴_𝑥𝛿²_𝑦𝑉^𝑛+1_𝑖,𝑗+𝑉^𝑛_𝑖,𝑗2−𝛽𝐴_𝑦𝛿²_𝑥+𝐴_𝑥𝛿²_𝑦𝑈^𝑛+1_𝑖,𝑗+𝑈^𝑛_𝑖,𝑗2+𝑐²Δ𝑡²4𝛿²_𝑥𝛿²_𝑦𝑉^𝑛+1_𝑖,𝑗−𝑉^𝑛_𝑖,𝑗Δ𝑡=𝐴_𝑥𝐴_𝑦𝑆^𝑛+1/2_𝑖,𝑗. (48) Noting that 𝐴_𝑥𝑈^𝑛_𝑖,𝑗=(1+(ℎ²/12)𝛿²_𝑥)𝑈^𝑛_𝑖,𝑗,𝐴_𝑦𝑈^𝑛_𝑖,𝑗=(1+(ℎ²/12)𝛿²_𝑦)𝑈^𝑛_𝑖,𝑗, (48) can be written as ℎ1+²𝛿12^2𝑥ℎ1+²𝛿12^2𝑦𝑉^𝑛+1_𝑖,𝑗−𝑉^𝑛_𝑖,𝑗ℎΔ𝑡−𝛼1+²𝛿12^2𝑦𝛿²_𝑥+ℎ1+²𝛿12^2𝑥𝛿²_𝑦𝑉^𝑛+1_𝑖,𝑗+𝑉^𝑛_𝑖,𝑗2ℎ−𝛽1+²𝛿12^2𝑦𝛿²_𝑥+ℎ1+²𝛿12^2𝑥𝛿²_𝑦𝑈^𝑛+1_𝑖,𝑗+𝑈^𝑛_𝑖,𝑗2+𝑐²Δ𝑡²4𝛿²_𝑥𝛿²_𝑦𝑉^𝑛+1_𝑖,𝑗−𝑉^𝑛_𝑖,𝑗=ℎΔ𝑡1+²𝛿12^2𝑥ℎ1+²𝛿12^2𝑦𝑆^𝑛+1/2_𝑖,𝑗. (49) Discretizing (6) as above and letting 𝜉=𝑢−𝑈,𝜂=𝑣−𝑉, we obtain the error equations as follows: ℎ1+²𝛿12^2𝑥ℎ1+²𝛿12^2𝑦𝜂^𝑛+1_𝑖,𝑗−𝜂^𝑛_𝑖,𝑗ℎΔ𝑡−𝛼1+²𝛿12^2𝑦𝛿²_𝑥+ℎ1+²𝛿12^2𝑥𝛿²_𝑦𝜂^𝑛+1_𝑖,𝑗+𝜂^𝑛_𝑖,𝑗2ℎ−𝛽1+²𝛿12^2𝑦𝛿²_𝑥+ℎ1+²𝛿12^2𝑥𝛿²_𝑦𝜉^𝑛+1_𝑖,𝑗+𝜉^𝑛_𝑖,𝑗2+𝑐²Δ𝑡²4𝛿²_𝑥𝛿²_𝑦𝜂^𝑛+1_𝑖,𝑗−𝜂^𝑛_𝑖,𝑗Δ𝑡=𝑅^𝑛+1/2_𝑖,𝑗,𝜂^𝑛+1_𝑖,𝑗+𝜂^𝑛_𝑖,𝑗2=𝛾₁𝜉^𝑛+1_𝑖,𝑗+𝜉^𝑛_𝑖,𝑗2+𝛾₂𝜉^𝑛+1_𝑖,𝑗−𝜉^𝑛_𝑖,𝑗+Δ𝑡^∼𝑅^𝑛+1/2_𝑖,𝑗, (50)(51) where the truncation errors 𝑅^𝑛+1/2_𝑖,𝑗=Δ𝑡²1𝐴24^𝑥𝐴_𝑦𝜕²𝑣𝜕𝑡³−𝛼8𝐴_𝑦𝜕⁴𝑣𝜕𝑡²𝜕𝑥²+𝐴_𝑥𝜕⁴𝑣𝜕𝑡²𝜕𝑦²−𝛽8𝐴_𝑦𝜕⁴𝑢𝜕𝑡²𝜕𝑥²+𝐴_𝑥𝜕⁴𝑢𝜕𝑡²𝜕𝑦²+𝑐²Δ𝑡⁴4𝜕⁵𝑣𝜕𝑡𝜕𝑥²𝜕𝑦²−𝛼ℎ⁴𝐴240^𝑦𝜕⁶𝑣𝜕𝑥⁶+𝐴_𝑥𝜕⁶𝑣𝜕𝑦⁶−𝛽ℎ⁴𝐴240^𝑦𝜕⁶𝑢𝜕𝑥⁶+𝐴_𝑥𝜕⁶𝑢𝜕𝑦⁶𝑥_𝑖,𝑦_𝑗,𝑡^𝑛+1/2+𝑂Δ𝑡⁴+Δ𝑡²ℎ⁴+ℎ⁶,_∼𝑅^𝑛+1/2_𝑖,𝑗=Δ𝑡²18𝜕²𝑣𝜕𝑡²−𝛾₁8𝜕²𝑢𝜕𝑡²−𝛾₂𝜕24³𝑢𝜕𝑡³𝑥_𝑖,𝑦_𝑗,𝑡^𝑛+1/2+𝑂Δ𝑡⁴. (52) Hence there exist positive constants 𝐶^′,𝐶^′′ such that 𝑅^𝑛+1/2_𝑖,𝑗≤𝐶^′Δ𝑡²+ℎ⁴,𝑅^∼𝑛+1/2_𝑖,𝑗≤𝐶^′′Δ𝑡². (53)

Lemma 3 (see [1]). For 𝑤∈𝑉_ℎ, the following inequalities hold: 14ℎ²𝛿_𝑥𝑤²≤‖𝑤‖²,14ℎ²𝛿_𝑦𝑤²≤‖𝑤‖²,14ℎ²𝛿_𝑥𝛿_𝑦𝑤²≤𝛿_𝑥𝑤²,14ℎ²𝛿_𝑥𝛿_𝑦𝑤²≤𝛿_𝑦𝑤². (54)

From Lemmas 1 and 3 and Cauchy-Schwarz inequality, multiplying ((𝜂^𝑛+1_𝑖,𝑗+𝜂^𝑛_𝑖,𝑗)/2)ℎ² to both sides of (50) and computing the inner product, we have 𝐼₁=ℎ1+²𝛿12^2𝑥ℎ1+²𝛿12^2𝑦𝛿_𝑡𝜂^𝑛,𝜂^𝑛+1+𝜂^𝑛2=𝛿_𝑡𝜂^𝑛,𝜂^𝑛+1+𝜂^𝑛2+ℎ²𝛿12^2𝑥𝛿_𝑡𝜂^𝑛,𝜂^𝑛+1+𝜂^𝑛2+ℎ²𝛿12^2𝑦𝛿_𝑡𝜂^𝑛,𝜂^𝑛+1+𝜂^𝑛2+ℎ⁴𝛿144^2𝑥𝛿²_𝑦𝛿_𝑡𝜂^𝑛,𝜂^𝑛+1+𝜂^𝑛2=1𝜂2Δ𝑡^𝑛+12−𝜂^𝑛2−ℎ²𝛿24Δ𝑡^𝑥𝜂^𝑛+12−𝛿_𝑥𝜂^𝑛2−ℎ²𝛿24Δ𝑡^𝑦𝜂^𝑛+12−𝛿_𝑦𝜂^𝑛2+ℎ⁴𝛿288Δ𝑡^𝑥𝛿_𝑦𝜂^𝑛+12−𝛿_𝑥𝛿_𝑦𝜂^𝑛2,𝐼₂₁𝛼ℎ=−1+²𝛿12^2𝑦𝛿²_𝑥𝜂^𝑛+1+𝜂^𝑛2,𝜂^𝑛+1+𝜂^𝑛2=𝛼𝛿_𝑥𝜂^𝑛+1+𝜂^𝑛2²−𝛼ℎ²𝛿12^𝑥𝛿_𝑦𝜂^𝑛+1+𝜂^𝑛2²≥𝛼ℎ²6𝛿_𝑥𝛿_𝑦𝜂^𝑛+1+𝜂^𝑛2²𝐼≥0,₂₂𝛼ℎ=−1+²𝛿12^2𝑥𝛿²_𝑦𝜂^𝑛+1+𝜂^𝑛2,𝜂^𝑛+1+𝜂^𝑛2=𝛼𝛿_𝑦𝜂^𝑛+1+𝜂^𝑛2²−𝛼ℎ²𝛿12^𝑥𝛿_𝑦𝜂^𝑛+1+𝜂^𝑛2²≥𝛼ℎ²6𝛿_𝑥𝛿_𝑦𝜂^𝑛+1+𝜂^𝑛2²𝐼≥0,₃₁𝛽ℎ=−1+²𝛿12^2𝑦𝛿²_𝑥𝜉^𝑛+1+𝜉^𝑛2,𝜂^𝑛+1+𝜂^𝑛2=−𝛽𝛿²_𝑥𝜉^𝑛+1+𝜉^𝑛2,𝜂^𝑛+1+𝜂^𝑛2−𝛽ℎ²𝛿12^2𝑥𝛿²_𝑦𝜉^𝑛+1+𝜉^𝑛2,𝜂^𝑛+1+𝜂^𝑛2=𝛽𝛾₁𝛿_𝑥𝜉^𝑛+1+𝜉^𝑛2,𝛿_𝑥𝜉^𝑛+1+𝜉^𝑛2+𝛽𝛾₂𝛿_𝑥𝜉^𝑛+1+𝜉^𝑛2,𝛿_𝑥𝜉^𝑛+1−𝜉^𝑛𝛿Δ𝑡+𝛽_𝑥𝜉^𝑛+1+𝜉^𝑛2,𝛿𝑅_𝑥^∼𝑛+1/2−𝛽𝛾₁ℎ²𝛿12^𝑥𝛿_𝑦𝜉^𝑛+1+𝜉^𝑛2,𝛿_𝑥𝛿_𝑦𝜉^𝑛+1+𝜉^𝑛2−𝛽𝛾₂ℎ²𝛿12^𝑥𝛿_𝑦𝜉^𝑛+1+𝜉^𝑛2,𝛿_𝑥𝛿_𝑦𝜉^𝑛+1−𝜉^𝑛−Δ𝑡𝛽ℎ²𝛿12^𝑥𝛿_𝑦𝜉^𝑛+1+𝜉^𝑛2,𝛿_𝑥𝛿𝑅_𝑦^∼𝑛+1/2≥𝛽𝛾₂𝛿2Δ𝑡^𝑥𝜉^𝑛+12−𝛿_𝑥𝜉^𝑛2−ℎ²12𝛽𝛾₂𝛿2Δ𝑡^𝑥𝛿_𝑦𝜉^𝑛+12−𝛿_𝑥𝛿_𝑦𝜉^𝑛2−𝛽4𝛿_𝑥𝜉^𝑛+12+𝛿_𝑥𝜉^𝑛2+𝛿𝑅_𝑥^∼𝑛+1/22−𝛽ℎ²𝛿48^𝑥𝛿_𝑦𝜉^𝑛+12+𝛿_𝑥𝛿_𝑦𝜉^𝑛2+𝛿_𝑥𝛿𝑅_𝑦^∼𝑛+1/22≥𝛽𝛾₂𝛿2Δ𝑡^𝑥𝜉^𝑛+12−𝛿_𝑥𝜉^𝑛2−ℎ²12𝛽𝛾₂𝛿2Δ𝑡^𝑥𝛿_𝑦𝜉^𝑛+12−𝛿_𝑥𝛿_𝑦𝜉^𝑛2−𝛽32𝛿_𝑥𝜉^𝑛+12+2𝛿_𝑥𝜉^𝑛2+𝛿𝑅_𝑥^∼𝑛+1/22,𝐼₃₂𝛽ℎ=−1+²𝛿12^2𝑥𝛿²_𝑦𝜉^𝑛+1+𝜉^𝑛2,𝜂^𝑛+1+𝜂^𝑛2≥𝛽𝛾₂𝛿2Δ𝑡^𝑦𝑒^𝑛+12−𝛿_𝑦𝑒^𝑛2−ℎ²12𝛽𝛾₂𝛿2Δ𝑡^𝑥𝛿_𝑦𝑒^𝑛+12−𝛿_𝑥𝛿_𝑦𝑒^𝑛2−𝛽32𝛿_𝑦𝜉^𝑛+12+2𝛿_𝑦𝜉^𝑛2+𝛿𝑅_𝑦^∼𝑛+1/22,𝐼₄=𝑐²Δ𝑡4𝛿²_𝑥𝛿²_𝑦𝜂^𝑛+1−𝜂^𝑛,𝜂^𝑛+1+𝜂^𝑛2=𝑐²Δ𝑡8𝛿_𝑥𝛿_𝑦𝜂^𝑛+12−𝛿_𝑥𝛿_𝑦𝜂^𝑛2,𝐼₅=𝑅^𝑛+1/2,𝜂^𝑛+1+𝜂^𝑛2≤14𝑅^𝑛+1/22+2𝜂^𝑛+12+2𝜂^𝑛2. (55)

Multiplying 2Δ𝑡 to both sides of the result derived from (55), we have 𝜂^𝑛+12−ℎ²𝛿12^𝑥𝜂^𝑛+12−ℎ²𝛿12^𝑦𝜂^𝑛+12+ℎ⁴𝛿144^𝑥𝛿_𝑦𝜂^𝑛+12+𝛽𝛾₂𝛿_𝑥𝜉^𝑛+12−ℎ²𝛿12^𝑥𝛿_𝑦𝜉^𝑛+12+𝛿_𝑦𝜉^𝑛+12−ℎ²𝛿12^𝑥𝛿_𝑦𝜉^𝑛+12+𝑐²Δ𝑡²4𝛿_𝑥𝛿_𝑦𝜂^𝑛+12≤𝜂^𝑛2−ℎ²𝛿12^𝑥𝜂^𝑛2−ℎ²𝛿12^𝑦𝜂^𝑛2+ℎ⁴𝛿144^𝑥𝛿_𝑦𝜂^𝑛2+𝛽𝛾₂𝛿_𝑥𝜉^𝑛2−ℎ²𝛿12^𝑥𝛿_𝑦𝜉^𝑛2+𝛿_𝑦𝜉^𝑛2−ℎ²𝛿12^𝑥𝛿_𝑦𝜉^𝑛2+2𝛽Δ𝑡32𝛿_𝑥𝜉^𝑛+12+2𝛿_𝑥𝜉^𝑛2+𝛿𝑅_𝑥^∼𝑛+1/22+2𝛿_𝑦𝜉^𝑛+12+2𝛿_𝑦𝜉^𝑛2+𝛿𝑅_𝑦^∼𝑛+1/22+𝑐²Δ𝑡²4𝛿_𝑥𝛿_𝑦𝜂^𝑛2+Δ𝑡2𝑅^𝑛+1/22+2𝜂^𝑛+12+2𝜂^𝑛2. (56) Summing for 𝑛 to both sides of (56), and by Lemma 3, we obtain 13𝜂^𝑛2+2𝛽𝛾₂3𝛿_𝑥𝜉^𝑛2+𝛿_𝑦𝜉^𝑛2+ℎ⁴𝛿144^𝑥𝛿_𝑦𝜉^𝑛2≤2Δ𝑡^𝑛−1_𝑙=1𝜂^𝑛2+8𝛽Δ𝑡3^𝑛−1_𝑙=1𝛿_𝑥𝜉^𝑛2+𝛿_𝑦𝜉^𝑛2+Δ𝑡2^𝑛−1_𝑙=1𝑅^𝑛+1/22+2𝛽Δ𝑡3^𝑛−1_𝑙=1𝛿𝑅_𝑥^∼𝑛+1/22+𝛿𝑅_𝑦^∼𝑛+1/22. (57) Noticing that 𝜉⁰=𝜂⁰=0, by Gronwall Lemma, we obtain the following theorem.

Theorem 4. Assume that 𝑢,𝑣 are the accurate solutions of (6)–(9) with sufficient smoothness and 𝑈,𝑉 are the difference solutions of (41)–(43). Let 𝜉=𝑢−𝑈,𝜂=𝑣−𝑉, then there exists a positive constant 𝐾^′ independent of Δ𝑡 and ℎ, such that max_{[]𝑛≤𝑇/Δ𝑡}𝜂^𝑛+𝛿_𝑥𝜉^𝑛+𝛿_𝑦𝜉^𝑛+ℎ⁴𝛿144^𝑥𝛿_𝑦𝜉^𝑛≤𝐾^′Δ𝑡²+ℎ⁴. (58)

5. Numerical Experiment

In this section, we exemplify two numerical examples to illustrate the effectiveness of the present ADI scheme. Lele [9] discussed some finite differencing errors for first and second derivatives by establishing the relationship of wave number (for differential equation) and numerical wave number (for difference scheme) and indicated that standard Pade scheme (which is known as compact finite difference scheme in this paper) is a good choice to discrete differential equations. In this paper, we focus on the efficient alternating direction technique for solving second-order wave equations. Because the compact method in spatial direction is standard, the scheme is effective.

The computer language used for the programming is Fortran, and the programs are performed on professional computer (CPU is 3.10 GHz, and Memory is 3.47 G). For convenience, we denote the present second-order ADI scheme (16)–(18) by PR, the present compact ADI scheme (41)–(43) by CPR, the second-order ADI scheme (21)–(22) by Dou, the compact ADI scheme (46)–(47) by CDou, and the three-level LOD scheme by LOD.

Example 5. Let 𝛾₁=1.0,𝛾₂=1.0,𝑎=1.0,𝑞=1.0, and 𝑆=(𝜋²−3𝜋+2)e^−𝜋𝑡sin(𝑥+𝑦) in (1), then the true solution is 𝑢=e^−𝜋𝑡sin(𝑥+𝑦). We chose different grids and compared the accuracy of the computed solutions with other schemes. The solutions were obtained for 𝑡_max=2.0. The quantities that we compared are the maximum 𝐿²-norm and the total elapsed time (CPU) in seconds.

Tables 1 and 2 display the 𝐿²-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. In Figures 1 and 2, we plot the errors and cost time at various grids. The figures show the superiority of the compact ADI over other schemes (second-order scheme). The CPR has much higher and more efficiency than CDou, which may be the best scheme to solve above questions.

Table 1 

The 𝐿²-norm and CPU time for the second-order schemes.

Table 2 

The 𝐿²-norm and CPU time for the fourth-order schemes.

Figure 1 

The error of all schemes for various grid sizes when Δ𝑡=ℎ2.

Figure 2 

The CPU time of all schemes for various grid sizes when Δ𝑡=ℎ2.

Example 6. Let 𝛾₁=1.0,𝛾₂=1.0,𝑎=1.0,𝑞=1.0, and 𝑆=𝜋²sin5𝜋𝑥sin2𝜋𝑦[34+(0.5+17𝜋²)cos𝜋²𝑡+(17−0.5𝜋²)sin𝜋²𝑡] in (1). The exact solution is 𝑢=(1+0.5sin𝜋²𝑡)sin5𝜋𝑥cos3𝜋𝑦, and the initial and boundary conditions can be obtained based on the exact solution. We chose the grids 25×25 and 50×50, and the time step with Δ𝑡=ℎ². The solutions were obtained for 𝑡_max=5.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 𝐿²-norm and the total elapsed time (CPU) in seconds.

Table 3 lists 𝑡 the total elapsed time (CPU) in seconds and the 𝐿²-norm error are listed in Table 3 and shown in Figures 3 and 4. One can see that the error from the same order scheme is very close. From the table and figures, we can see that the present fourth-order scheme not only has high accuracy but also converges as fast as other schemes, especially for CPR scheme.

Table 3 

The 𝐿²-norm and CPU time for all schemes.

Figure 3 

The error of all schemes for various grid sizes when Δ𝑡=ℎ2.

Figure 4 

The CPU time of all schemes for various grid sizes when Δ𝑡=ℎ2.

From Figures 1 and 3, it is easily to find that the tested methods degrade to first-order method as the grid gets coarse. In all schemes, we add perturbing term (𝑐²Δ𝑡/4)𝛿²_𝑥𝛿²_𝑦(𝑉^𝑛+1_𝑖,𝑗−𝑉^𝑛_𝑖,𝑗) that is, an approximation to (𝑐²Δ𝑡²/4)(𝜕5𝑣/𝜕𝑡𝜕𝑥²𝜕𝑦²). When step ℎ becomes large, the time step Δ𝑡 varies correspondingly because of Δ𝑡=ℎ for the second-order schemes or Δ𝑡=ℎ² for the compact schemes. Hence, the effect of perturbing term is enhanced, which results in the schemes seem to degrade to first-order methods as the grid gets coarse. We can choose suitable ℎ and get required accuracy.

6. 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 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.

The present scheme can be applied to the following variable coefficient differential models on domain Ω=[0,1]²: 𝛾₁𝜕𝑢𝜕𝑡²+𝛾₂𝜕²𝑢𝜕𝑡²()−𝜕−∇⋅𝐴∇𝑢()(],𝑢𝜕𝑡∇⋅𝑄∇𝑢=𝑆,𝑥,𝑦∈Ω,𝑡∈0,𝑇𝑥,𝑦,0=𝑢₀𝑥,𝑦,𝑢_𝑡,𝑢,(],𝑥,𝑦,0=𝜓𝑥,𝑦𝑥,𝑦∈Ω,𝑥,𝑦,𝑡=𝜑𝑥,𝑦,𝑡𝑥,𝑦∈𝜕Ω,𝑡∈0,𝑇 (59) where 𝐴=𝐴(𝑥,𝑦) and 𝑄=𝑄(𝑥,𝑦) are diagonal, nonnegative diffusion tensors.

Acknowledgment

This work is partly supported by National Magnetic Confinement Fusion Science Program, Grant no. 2011GB112004.