Abstract
In this paper, we develop a highly accurate and efficient finite difference scheme for solving the two-dimensional (2D) wave equation. Based on the local one-dimensional (LOD) method and Padé difference approximation, a fourth-order accuracy explicit compact difference scheme is proposed. Then, the Fourier analysis method is used to analyze the stability of the scheme, which shows that the new scheme is conditionally stable and the Courant-Friedrichs-Lewy (CFL) condition is superior to most existing methods of equivalent order of accuracy in the literature. Finally, numerical experiments demonstrate the high accuracy, stability, and efficiency of the proposed method.
1. Introduction
Wave equation is mainly used to describe various wave phenomena in nature, such as sound wave, light wave, and water wave. Due to complex physical background, it is very difficult to obtain the exact solutions of practical problems. Therefore, researches on numerical solutions of initial value (or initial-boundary value) problems for wave equation have extremely important theoretical value and practical significance.
The finite difference method is one of the most important and popular methods for solving the wave equation. It has many advantages, such as it is easy to be implemented, has low memory requirement, and has a high computing speed [1–8]. The second-order central difference scheme is the most commonly used, but when the number of spatial sampling points per wavelength is too small, the method has serious dispersion. In order to effectively suppress numerical dispersion, Yang et al. [9] proposed a local interpolation strategy based on the Runge-Kutta method. Later, for solving large-scale problems, Yang et al. [10] proposed an approximate analytic central difference method. However, both methods are conditionally stable; their maximum values of CFL conditions are 0.5080 [9] and 0.8440 [10], respectively. Most of the existing schemes with fourth-order accuracy in both space and time have very strict stability conditions. For instance, Wang et al. [11] developed a new finite difference stencil where the range of the CFL number is 0.7071. Feo et al. [12] combined novel tensor mimetic discretizations in space and a leapfrog approximation in time to produce an explicit scheme; the range of the CFL number is 0.8660. In addition, many implicit finite difference schemes are also used to solve the 2D wave equations [13–16]. For instance, Wang et al. [13] introduced two efficient fourth-order implicit time-space-domain finite difference schemes by using linear and nonlinear optimization methods. The high-order compact (HOC), compact Padé difference, and noncompact Padé difference implicit schemes were proposed in Ref. [14], all of which have a truncation error of . Since these implicit difference methods require iterations at each time step, this greatly increases the computational cost. In order to improve the computational efficiency, a combination of the finite difference method and splitting method is used by many researchers [17, 18]. Two types of classical splitting methods are the alternating direction implicit (ADI) method and the LOD method. Both of these methods simplify the multidimensional problems into a series of one-dimensional (1D) problems that form the tridiagonal systems and are easy to be solved. The ADI method was originally introduced by Peaceman and Rachford [19] to solve parabolic and elliptic equations. It was also used to solve hyperbolic equations [20, 21]. Up to now, many ADI methods have been studied to solve the wave equations [22–26]. For instance, Das et al. [22] and Liao et al. [23] used a combination of ADI and the Padé approximation to derive a difference scheme with fourth-order accuracy and low numerical dispersion. Their maximum values of CFL conditions are 0.7657 [22] and 0.7321 [23], respectively. Qin [24] constructed a compact Douglas scheme based on the ADI method, which has fourth-order accuracy in space, but only second-order accuracy in time.
On the other hand, the LOD method was originally introduced by Samarskii [27] for solving the 2D hyperbolic equations. The LOD technique has been proven to be very useful in reducing computer memory and computational cost and improving stability [28, 29]. We notice that many researchers mainly study the homogeneous wave equation. For instance, Yun et al. [30] applied a new fourth-order accuracy optimal nearly analytic splitting method to solve the 2D acoustic wave equation based on the LOD method. Sim et al. [31] proposed a nearly analytic symplectic partitioned Runge-Kutta method based on the LOD technique, which has second-order temporal accuracy and fourth-order spatial accuracy. Very recently, Zhang [29] considered the nonhomogeneous wave equation and set an unknown parameter in front of the source term. Then, the value of was calculated by solving a series of large matrices, so as to obtain an implicit scheme with fourth-order accuracy in both time and space. This method increases computational complexity and reduces computational efficiency. Therefore, in this paper, we are aiming at developing an explicit compact difference method for solving the 2D nonhomogeneous wave equation based on the LOD method. The merits of the present method include high-order accuracy, better stability, and easy treatment for the nonhomogeneous term.
The structure of this paper is arranged as follows. In Section 2, the derivation of the new explicit HOC difference scheme based on the LOD method is given and the corresponding numerical algorithm is introduced. In Section 3, the stability of the new scheme is analyzed. In Section 4, the numerical examples are employed to verify our theoretical analysis results. Finally, the conclusion is given in Section 5.
2. Derivation of the New Scheme
The 2D wave equation with variable wave velocity is as follows: with the initial conditions and the Dirichlet boundary conditions
where represents pressure and is the wave velocity. is the source term. , , and are known smooth functions. , , and represent the boundary of . are constants, and . is the time region.
is divided into subintervals by equidistant grid, and the space step length is represented by . is also equally divided into subintervals, and the time step length is . We use which represents the mesh point, .
The LOD method is used to ideally split the 2D wave equation (1) into two 1D wave equations (27) as follows:
To advance the solution from to , we assume that Equation (4) holds from to and Equation (5) holds from to .
Firstly, we consider Equation (4) at the grid point , i.e.,
For the time second derivative on the left side of Equation (6), we adopt the following expression:
According to the Equation (4), we have
We define the central difference operator for the spatial second derivative as follows:
Substituting Equations (8) and (9) into Equation (7), we obtain
It can be rewritten as
Defining the grid ratio , substituting Equation (4) into Equation (11), and omitting the truncation error term, we obtain
The function and its derivative are known at each grid point. Next, we use the fourth-order Padé scheme given in Ref. [32] to calculate the unknown term in Equation (12) as follows:
in which the nodes of the on the boundaries can be obtained by using Equations (3) and (6):
Then, we consider Equation (5) at the grid point , i.e.,
Using a similar treatment process as for Equation (6), omitting the truncation error terms, we can obtain
The internal nodes of the can also be computed by the fourth-order Padé scheme [32].
in which, the nodes of the on the boundaries can be obtained by using Equations (3) and (16):
Applying Equations (12)–(20), we can complete the entire calculation from to for Equation (1). The scheme is an explicit HOC difference scheme based on the LOD method, which can be denoted as the EHOC-LOD scheme. According to the derivation process, we know that the EHOC-LOD scheme is a four-level scheme with fourth-order accuracy in both time and space. Therefore, two start-up time steps need to be calculated at and . To approximate the start-up time steps and , using Taylor’s expansions, we have
For simplicity, we denote the time-independent functions as the following notations:
in which can represent , etc.
By using Equations (1) and (2), we can get
Substituting Equations (24)–(26) into Equations (21) and (22), respectively, and omitting the truncation error, we obtain
Based on the description above, an algorithm for solving Equation (1) by using the EHOC-LOD scheme is given as follows:
Step 1. Let , give the initial time step , and the start-up time steps and are computed by using Equations (27) and (28) .
Step 2. The values of are computed by Equations (13)–(15) . Then, the function values are computed by using Equation (12) . The function values and are computed by Equation (3) .
Step 3. The values of are computed by Equations (18)–(20) . Then, the function values are computed by using Equation (17) . The function values and are computed by Equation (3) .
Step 4. Let , repeat Steps 2 and 3 until time reaches the final moment, and the computation is terminated.
3. Stability Analysis
In this section, we analyze the stability of the EHOC-LOD scheme by using the Fourier analysis method.
Lemma 1 (see [33]). The sufficient and necessary condition for the roots of the quadratic equation with real coefficients to be not greater than one is , .
Theorem 1. The EHOC-LOD scheme in Section 2 is stable if
in which, .
Proof. Let , , and , where are amplitudes, are phase angles, and . We assume that is exact and does not produce error. Equations (13) and (18) can conclude that
By using Euler’s formula, , Equations (30) and (31) turn to
Let and , and Equation (12) is written in matrix form:
Substituting into the formula above, we get
Substituting Equation (32) into Equation (35), we have
Similarly, Equation (17) can be treated as
Substituting Equation (36) into Equation (37), the error propagation matrix can be obtained as where
The characteristic equation can be obtained as
According to Lemma 1, . Due to and having the same range of values, we only need to solve , or .
Case 1. When , we only need to solve , i.e.,
Firstly, for the right side of the inequality (41) we have
Then, for the left side of the inequality (41), Letting , the inequality (44) can be expressed as
Letting , we assume that is a quadratic function of . Then, the function image is a parabola with an upward opening, the axis of symmetry is , the ordinate of the vertex is . In this case, we only need to determine the position of the symmetry axis, so that on .
Firstly, since , the axis of symmetry always holds.
Secondly, when the symmetry axis , we have . is an increasing function; we only need , so or . But , so .
Finally, when the axis of symmetry , we have . We only need , so . But , so there is no solution.
In summary, when , i.e., , the EHOC-LOD scheme is stable.
Case 2. When , we only need to solve
Firstly, for the right side of the inequality (46), According the result of Case 1, we can get that
Then, for the left side of the inequality (46), Letting , the inequality (49) be expressed as
Letting , we assume that is a quadratic function of . Then, the function image is a parabola with an upward opening, the axis of symmetry is , and the ordinate of the vertex is . In this case, we only need to determine the position of the symmetry axis, so that on .
Firstly, since , the axis of symmetry always holds.
Secondly, when the axis of symmetry , we have . is an increasing function; we only need , so or . But , so .
Finally, when the axis of symmetry , we have . We only need , so . But , so .
In summary, due to , there is no solution for this case.
Combining Cases 1 and 2, we can conclude that when , the EHOC-LOD scheme is stable. This completes the proof.
Table 1 gives the stability conditions of various schemes in Refs. [9–12, 15, 16, 22, 23, 28, 30] and the present EHOC-LOD scheme for the 2D wave equation. Among them, some methods are aimed at the constant wave number wave equation such as Refs. [9, 10, 12, 16, 22, 28, 30]; represents the absolute value of the constant wave coefficient. The others deal with the wave equation with variable wave coefficients such as Refs. [11, 15, 23]; represents the maximum value of the variable wave coefficients. The theoretical analysis result of the stability condition of the EHOC-LOD scheme in this paper is . So, it is more flexible than all the other schemes listed in the literature.
4. Numerical Experiments
This section uses some numerical examples to confirm the accuracy and stability of the proposed EHOC-LOD method. All programs are written in Fortran 90 language by using double-precision arithmetic and run on an Intel Core i5-4210U [email protected] GHz 2.40 GHz desk computer with 4 GB memory.
and norm errors and the rate of convergence are obtained by the following definitions:
The CPU time for all the computation is given in seconds.
Example 1 (see [28]). The exact solution is
Example 1 is a 2D homogeneous wave equation with constant coefficients. Table 2 shows the and norm errors for the EHOC-LOD scheme with and . In order to compare, we give the calculated results of the LOD scheme proposed by Zhang et al. in Ref. [28] in Table 2. The results show that the and norm errors of the EHOC-LOD scheme are smaller than those of the LOD scheme [28]. In order to verify the computed accuracy of the EHOC-LOD scheme in space direction, in Figure 1, we draw log-log plots for the and norm errors in the space direction for the LOD scheme [28] and the EHOC-LOD scheme. That displays that the slope of the line of the EHOC-LOD scheme is close to 4, which means that the EHOC-LOD scheme can converge with the fourth order in the space direction. When the CFL number is fixed to 0.2 and , the and norm errors are listed in Table 3. Obviously, the numerical solution of the EHOC-LOD scheme is a bit better than that of the LOD scheme [28]. Then, two log-log plots for the and norm errors in the time direction for the LOD scheme [28] and the EHOC-LOD scheme are shown in Figure 2. Figure 2 shows that the slope of the line of the EHOC-LOD scheme exceeds 4, which means that the EHOC-LOD scheme can converge with the fourth order in the time direction.
Example 2 (see [24]). The exact solution is
(a)
(b)
(a)
(b)
Example 2 is a 2D nonhomogeneous wave equation with constant coefficients. In order to make a better comparison, we program the algorithms of the CPR, CDOU, and CSchemeI schemes with second-order accuracy in time and fourth-order accuracy in space in Ref. [24] also in Fortran 90 language. We run the programs on the same computer to make the comparison fair. When , for different spatial grid step length, in Table 4, the norm errors and CPU time are shown by using the present EHOC-LOD scheme with , the CPR, CDOU, and CSchemeI schemes with in Ref. [24]. We can notice that the EHOC-LOD scheme obtains fourth-order accuracy in both time and space, but those schemes in Ref. [24] must take to achieve fourth-order accuracy in the space direction. Moreover, the present method gets a more accurate solution. When we define the time as and the number of grids as , the total number of time advance steps of the method in this paper is , while the total number of time advance steps of those methods in Ref. [24] is ; therefore, the EHOC-LOD scheme expends much less CPU time than those methods in Ref. [24]. With the increase of time and grid number , the method in this paper has more significant advantages.
Then, we test the stability condition of the EHOC-LOD method. The norm error of the EHOC-LOD scheme with for different and is shown in Table 5. As the time increases, when , the EHOC-LOD scheme is convergent. When , the EHOC-LOD scheme is divergent. These computed results are in good agreement with the theoretical analysis.
Example 3 (see [23]). The exact solution is
When , for various and , we compute the norm error and by using the NCV-CPD-ADI scheme in Ref. [23] and the EHOC-LOD scheme in Table 6. The EHOC-LOD scheme produces the numerical solution with a smaller error than the NCV-CPD-ADI scheme does [23]. In addition, The EHOC-LOD scheme can achieve the fourth-order accuracy in both time and space, while the accuracy of the NCV-CPD-ADI scheme [23] is lower than that of the fourth order.
The stability condition is given by , where for this example. Therefore, the grid ratio allowed by the EHOC-LOD scheme is . The grid ratio allowed by the NCV-CPD-ADI scheme [23] is . When we fix , then and . We take several numbers near ; Table 7 shows the results. When , the norm error of the NCV-CPD-ADI scheme [23] and the present scheme are both quite accurate. When , the norm error of the NCV-CPD-ADI scheme [23] increases faster, and eventually, it is divergent due to exceeding the stability range. However, the norm error of the EHOC-LOD scheme is still quite accurate and convergent. Therefore, the CFL condition of the scheme in this paper is superior to that of the NCV-CPD-ADI scheme [23].
To further verify the stability range of the EHOC-LOD scheme, Table 8 shows the norm error of the EHOC-LOD scheme with for different and . We notice that as time grows, when , the EHOC-LOD scheme is convergent. When , the EHOC-LOD scheme is divergent. The numerical experiment results are consistent with the theoretical ones.
Example 4 (see [23]). Example 4 is a