An accurate and efficient local one-dimensional method for the 3D acoustic wave equation

Mengling Wu; Yunzhi Jiang; Yongbin Ge

doi:10.1515/dema-2022-0148

Open Access Published by De Gruyter Open Access September 22, 2022

An accurate and efficient local one-dimensional method for the 3D acoustic wave equation

Mengling Wu , Yunzhi Jiang and Yongbin Ge

From the journal Demonstratio Mathematica

https://doi.org/10.1515/dema-2022-0148

Abstract

We establish an accurate and efficient scheme with four-order accuracy for solving three-dimensional (3D) acoustic wave equation. First, the local one-dimensional method is used to transfer the 3D wave equation into three one-dimensional wave equations. Then, a new scheme is obtained by the Padé formulas for computation of spatial second derivatives and the correction of the truncation error remainder for discretization of temporal second derivative. It is compact and can be solved directly by the Thomas algorithm. Subsequently, the Fourier analysis method and the Lax equivalence theorem are employed to prove the stability and convergence of the present scheme, which shows that it is conditionally stable and convergent, and the stability condition is superior to that of most existing numerical methods of equivalent order of accuracy in the literature. It allows us to reduce computational cost with relatively large time step lengths. Finally, numerical examples have demonstrated high accuracy, stability, and efficiency of our method.

Keywords: local one-dimensional; wave equation; variable coefficient; high accuracy compact difference scheme; stability and convergence

MSC 2010: 65M06; 65M12; 35L05

1 Introduction

The numerical solution of wave equation plays an important role in seismic wave propagation, full waveform inversion, and seismic imaging [1,2,3]. Therefore, the research on numerical solutions of initial values (or initial-boundary values) problems of wave equation has extremely important theoretical value and practical significance.

Due to its high accuracy, low memory, and fast computing speed [4,5, 6,7], the finite difference method is a powerful tool for acoustic or seismic wave simulations, especially for models with complex geological structures. The temporally second-order and spatially fourth-order accuracy schemes are commonly applied to reduce the numerical dispersion caused by the discretization of time and space derivatives. For instance, based on the Taylor expansion, both Liang et al. [1] and Liu and Sen [2] used the time-space domain method of plane wave to suppress the numerical dispersion. To balance the spatial and temporal accuracy of the schemes, many researchers have developed fourth-order finite difference schemes in both time and space [8,9,10, 11,12,13, 14,15]. For instance, Yang et al. [9] introduced an approximate analytic central difference method for solving large-scale problems. The polynomial interpolation method was also adopted to construct high-order compact (HOC) difference schemes [10,11]. Abdulkadir [14] introduced a noncompact implicit difference scheme. Compared with the compact difference schemes, the noncompact schemes have a wider stencil and the boundary conditions are difficult to be dealt with. In addition, for three-dimensional (3D) problems, high-order implicit difference schemes usually result in a series of large-scale sparse linear systems, which need iteration on each time step. To improve computational efficiency, two typical splitting methods, i.e., alternating direction implicit (ADI) and local one-dimensional (LOD) methods were proposed. Through these two methods, the multidimensional problems are simplified to a series of one-dimensional (1D) problems.

To efficiently solve parabolic and elliptic equations, Peaceman and Rachford [16] first proposed the ADI method. Until now, quite a lot ADI schemes have been utilized to study the wave equations [17,18,19, 20,21,22]. Liao et al. used the ADI method and Padé approximation to derive a low numerical dispersion scheme with fourth-order accuracy for solving two-dimensional (2D) [17] and 3D [18] wave equations with constant coefficients. Their Courant-Friedrichs-Lewy (CFL) condition numbers are 0.7657 [17] and 0.6079 [18], respectively. Later, Liao et al. extended this method to the 2D and 3D wave equations with variable coefficients [19,22]. Their CFL condition numbers are 0.7321 [19] and 0.5770 [22], respectively. Up to now, the ADI methods have successfully solved various wave propagation, but the stability condition of this method is still quite strict.

The LOD method was originally introduced by Samarskii [23] for solving the hyperbolic equations. The LOD technique has many advantages, such as less computational cost, low computer memory, and good stability [24,25]. For the nonhomogeneous wave equation with variable coefficient, the high-order accuracy ADI methods have been extensively studied, while the high-order LOD methods have not been developed to the same extent. We notice that many researchers mainly study homogeneous wave equations [26,27, 28,29]. Zhang and Jiang presented some fourth-order LOD schemes to solve the 2D [26] and 3D wave equations [27]. They are conditionally stable, and the CFL condition numbers are 0.7321 [26] and 0.6160 [27], respectively. Later, for the nonhomogeneous wave equation, Zhang [29] considered an implicit fourth-order accuracy scheme by setting an unknown parameter θ in front of the source term, which increases computational complexity and reduces computational efficiency. In this paper, on the basis of the LOD technique, we present a compact difference method for the 3D nonhomogeneous problem. The scheme uses only a three-point stencil to approximate the second derivatives of spaces. Moreover, it has the advantages of high spatial-temporal accuracy, better stability, and convenience in dealing with the nonhomogeneous term.

In Section 2, the derivation of a novel HOC method based on the LOD method is given. In Section 3, the stability, consistency, and convergence of the novel scheme are analyzed. In Section 4, numerical examples verify our theoretical results. Finally, the conclusion is given in Section 5.

2 Deduction of the HOC-LOD scheme

We consider the 3D nonhomogeneous variable velocity acoustic wave equation

(1) ∂ 2 u ∂ t 2 = v 2 ( x , y , z ) ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 + ∂ 2 u ∂ z 2 + f ( x , y , z , t ) , ( x , y , z , t ) ∈ Ω × ( 0 , T ] ,

with the initial conditions

(2) u ( x , y , z , 0 ) = φ ( x , y , z ) , ∂ u ( x , y , z , 0 ) ∂ t = ψ ( x , y , z ) , ( x , y , z ) ∈ Ω ,

and the boundary conditions

(3) u ( x , y , z , t ) = g ( x , y , z , t ) , ( x , y , z , t ) ∈ ∂ Ω × ( 0 , T ] ,

where Ω = [ c , d ] 3 ∈ R 3 and c < d . ∂ Ω is the boundary of Ω . u ( x , y , z , t ) ∈ C Ω , t 5 , 5 is unknown function. f ( x , y , z , t ) ∈ C Ω , t 3 , 3 is the source term. v ( x , y , z ) ∈ C 4 ( Ω ) is nonzero function that represents wave velocity. ψ ( x , y , z ) ∈ C 4 ( Ω ) , φ ( x , y , z ) ∈ C 4 ( Ω ) , and g ( x , y , z , t ) ∈ C Ω , t 5 , 5 are known functions.

To build a difference scheme, first, we divide [ c , d ] into N subintervals by equidistant grid, and the space step length is represented by h = d − c N . ( 0 , T ] is the time region, and it is equally divided into M subintervals, and the time step length is τ = T M . The grid points are denoted as ( x i , y j , z k , t n ) , where x i = c + i h , y j = c + j h , z k = c + k h , t n = n τ , n = 0 , … M , i , j , k = 0 , … N .

The LOD method is used to ideally split the 3D wave equation (1) into three 1D wave equations [23] as follows:

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

(5) 1 3 ∂ 2 u ∂ t 2 = v 2 ( x , y , z ) ∂ 2 u ∂ y 2 + 1 3 f ( x , y , z , t ) ,

(6) 1 3 ∂ 2 u ∂ t 2 = v 2 ( x , y , z ) ∂ 2 u ∂ z 2 + 1 3 f ( x , y , z , t ) .

To advance the solution from u ( x i , y j , z k , t n ) to u ( x i , y j , z k , t n + 1 ) , we assume that equation (4) holds from u ( x i , y j , z k , t n ) to u x i , y j , z k , t n + 1 3 , equation (5) holds from u x i , y j , z k , t n + 1 3 to u x i , y j , z k , t n + 2 3 , and equation (6) holds from u x i , y j , z k , t n + 2 3 to u ( x i , y j , z k , t n + 1 ) , respectively.

By analyzing the splitting error of equations (4)–(6) to equation (1) (refer to Appendix A for the corresponding detailed derivation), we can obtain the following conclusions:

Case I: when ∂ 3 u ∂ t ∂ x 2 − ∂ 3 u ∂ t ∂ z 2 ≠ 0 , the splitting error is O ( τ ) .
Case II: when ∂ 3 u ∂ t ∂ x 2 − ∂ 3 u ∂ t ∂ z 2 = 0 5 ∂ 4 u ∂ t 2 ∂ x 2 + 2 ∂ 4 u ∂ t 2 ∂ y 2 − 7 ∂ 4 u ∂ t 2 ∂ z 2 ≠ 0 , the splitting error is O ( τ 2 ) .
Case III: when ∂ 3 u ∂ t ∂ x 2 − ∂ 3 u ∂ t ∂ z 2 = 0 5 ∂ 4 u ∂ t 2 ∂ x 2 + 2 ∂ 4 u ∂ t 2 ∂ y 2 − 7 ∂ 4 u ∂ t 2 ∂ z 2 = 0 , the splitting error is O ( τ 4 ) .

Next, we construct a HOC difference scheme to solve the system (4)–(6). To this end, first, we consider the first equation as follows:

(7) 1 3 ∂ 2 u ∂ t 2 i , j , k n = v i , j , k 2 ∂ 2 u ∂ x 2 i , j , k n + 1 3 f i , j , k n .

For the time second derivative ∂ 2 u ∂ t 2 in equation (7), we adopt the following expression:

(8) ∂ 2 u ∂ t 2 i , j , k n = 9 u i , j , k n + 1 3 − 2 u i , j , k n + u i , j , k n − 1 3 τ 2 − τ 2 108 ∂ 4 u ∂ t 4 i , j , k n + O ( τ 4 ) .

According to equation (4), we have

(9) ∂ 4 u ∂ t 4 = ∂ 2 ∂ t 2 ∂ 2 u ∂ t 2 = 3 v 2 ∂ 2 ∂ x 2 ∂ 2 u ∂ t 2 + ∂ 2 f ∂ t 2 .

We define the central difference operator for the spatial second derivative as follows:

(10) δ x 2 u i , j , k n = u i + 1 , j , k n − 2 u i , j , k n + u i − 1 , j , k n h 2 .

Substituting equations (9) and (10) into equation (8),

∂ 2 u ∂ t 2 i , j , k n = 9 u i , j , k n + 1 3 − 2 u i , j , k n + u i , j , k n − 1 3 τ 2 − v i , j , k 2 τ 2 36 ∂ 2 ∂ x 2 ∂ 2 u ∂ t 2 i , j , k n − τ 2 108 ∂ 2 f ∂ t 2 i , j , k n + O ( τ 4 ) = 9 u i , j , k n + 1 3 − 2 u i , j , k n + u i , j , k n − 1 3 τ 2 − v i , j , k 2 τ 2 36 δ x 2 ∂ 2 u ∂ t 2 i , j , k n − τ 2 108 ∂ 2 f ∂ t 2 i , j , k n + O ( τ 4 + τ 2 h 2 ) = 9 u i , j , k n + 1 3 − 2 u i , j , k n + u i , j , k n − 1 3 τ 2 − v i , j , k 2 τ 2 36 h 2 ∂ 2 u ∂ t 2 i + 1 , j , k n + ∂ 2 u ∂ t 2 i − 1 , j , k n + v i , j , k 2 τ 2 18 h 2 ∂ 2 u ∂ t 2 i , j , k n − τ 2 108 ∂ 2 f ∂ t 2 i , j , k n + O ( τ 4 + τ 2 h 2 ) .

Multiplying by τ 2 on the both sides of aforementioned formula, it can be rewritten as follows:

(11) u i , j , k n + 1 3 = 2 u i , j , k n − u i , j , k n − 1 3 + v i , j , k 2 τ 4 324 h 2 ∂ 2 u ∂ t 2 i + 1 , j , k n − 2 ∂ 2 u ∂ t 2 i , j , k n + ∂ 2 u ∂ t 2 i − 1 , j , k n + τ 2 9 ∂ 2 u ∂ t 2 + τ 2 108 ∂ 2 f ∂ t 2 i , j , k n + O ( τ 6 + τ 4 h 2 ) .

Defining the grid ratio λ = τ / h , substituting equation (7) into equation (11) and omitting the truncation error term, we obtain

(12) u i , j , k n + 1 3 = 2 u i , j , k n − u i , j , k n − 1 3 + v i , j , k 2 τ 2 λ 2 108 v i + 1 , j , k 2 ∂ 2 u ∂ x 2 i + 1 , j , k n + v i − 1 , j , k 2 ∂ 2 u ∂ x 2 i − 1 , j , k n + v i , j , k 2 τ 2 3 − v i , j , k 4 τ 2 λ 2 54 ∂ 2 u ∂ x 2 i , j , k n + v i , j , k 2 τ 2 λ 2 324 ( f i + 1 , j , k n + f i − 1 , j , k n ) + τ 2 9 − v i , j , k 2 τ 2 λ 2 162 f i , j , k n + τ 4 972 ∂ 2 f ∂ t 2 i , j , k n .

The function f and its derivative ∂ 2 f ∂ t 2 are known at each grid point. Next, we use the fourth-order Padé formula in [30] to solve ∂ 2 u ∂ x 2 n in equation (12) as follows:

(13) 1 12 ∂ 2 u ∂ x 2 i + 1 , j , k n + 5 6 ∂ 2 u ∂ x 2 i , j , k n + 1 12 ∂ 2 u ∂ x 2 i − 1 , j , k n = u i + 1 , j , k n − 2 u i , j , k n + u i − 1 , j , k n h 2 + O ( h 4 ) , i = 1 , … , N − 1 , j , k = 0 , … N .

The boundaries of the ∂ 2 u ∂ x 2 i , j , k n can be obtained by equations (3) and (7):

(14) ∂ 2 u ∂ x 2 0 , j , k n = 1 3 v 0 , j , k 2 ∂ 2 u ∂ t 2 − f 0 , j , k n = 1 3 v 0 , j , k 2 ∂ 2 g ∂ t 2 − f 0 , j , k n , j , k = 0 , … N ,

(15) ∂ 2 u ∂ x 2 N , j , k n = 1 3 v N , j , k 2 ∂ 2 u ∂ t 2 − f N , j , k n = 1 3 v N , j , k 2 ∂ 2 g ∂ t 2 − f N , j , k n , j , k = 0 , … N .

Then, we consider the second equation as follows:

(16) 1 3 ∂ 2 u ∂ t 2 i , j , k n + 1 3 = v i , j , k 2 ∂ 2 u ∂ y 2 i , j , k n + 1 3 + 1 3 f i , j , k n + 1 3 .

By the similar treatment process, we can obtain

(17) u i , j , k n + 2 3 = 2 u i , j , k n + 1 3 − u i , j , k n + v i , j , k 2 τ 2 λ 2 108 v i , j + 1 , k 2 ∂ 2 u ∂ y 2 i , j + 1 , k n + 1 3 + v i , j − 1 , k 2 ∂ 2 u ∂ y 2 i , j − 1 , k n + 1 3 + v i , j , k 2 τ 2 3 − v i , j , k 4 τ 2 λ 2 54 ∂ 2 u ∂ y 2 i , j , k n + 1 3 + v i , j , k 2 τ 2 λ 2 324 f i , j + 1 , k n + 1 3 + f i , j − 1 , k n + 1 3 + τ 2 9 − v i , j , k 2 τ 2 λ 2 162 f i , j , k n + 1 3 + τ 4 972 ∂ 2 f ∂ t 2 i , j , k n + 1 3 .

The internal nodes of the ∂ 2 u ∂ y 2 i , j , k n + 1 3 can also be computed by the fourth-order Padé formula [30]:

(18) 1 12 ∂ 2 u ∂ y 2 i , j + 1 , k n + 1 3 + 5 6 ∂ 2 u ∂ y 2 i , j , k n + 1 3 + 1 12 ∂ 2 u ∂ y 2 i , j − 1 , k n + 1 3 = u i , j + 1 , k n + 1 3 − 2 u i , j , k n + 1 3 + u i , j − 1 , k n + 1 3 h 2 + O ( h 4 ) , j = 1 , … , N − 1 , i , k = 0 , … N .

The boundaries of the ∂ 2 u ∂ y 2 i , j , k n + 1 3 can be obtained by equations (3) and (16):

(19) ∂ 2 u ∂ y 2 i , 0 , k n + 1 3 = 1 3 v i , 0 , k 2 ∂ 2 u ∂ t 2 − f i , 0 , k n + 1 3 = 1 3 v i , 0 , k 2 ∂ 2 g ∂ t 2 − f i , 0 , k n + 1 3 , i , k = 0 , … N ,

(20) ∂ 2 u ∂ y 2 i , N , k n + 1 3 = 1 3 v i , N , k 2 ∂ 2 u ∂ t 2 − f i , N , k n + 1 3 = 1 3 v i , N , k 2 ∂ 2 g ∂ t 2 − f i , N , k n + 1 3 , i , k = 0 , … N .

Finally, we consider the last equation as follows:

(21) 1 3 ∂ 2 u ∂ t 2 i , j , k n + 2 3 = v i , j , k 2 ∂ 2 u ∂ z 2 i , j , k n + 2 3 + 1 3 f i , j , k n + 2 3 .

By the similar treatment process, we obtain

(22) u i , j , k n + 1 = 2 u i , j , k n + 2 3 − u i , j , k n + 1 3 + v i , j , k 2 τ 2 λ 2 108 v i , j , k + 1 2 ∂ 2 u ∂ z 2 i , j , k + 1 n + 2 3 + v i , j , k − 1 2 ∂ 2 u ∂ z 2 i , j , k − 1 n + 2 3 + v i , j , k 2 τ 2 3 − v i , j , k 4 τ 2 λ 2 54 ∂ 2 u ∂ z 2 i , j , k n + 2 3 + v i , j , k 2 τ 2 λ 2 324 f i , j , k + 1 n + 2 3 + f i , j , k − 1 n + 2 3 + τ 2 9 − v i , j , k 2 τ 2 λ 2 162 f i , j , k n + 2 3 + τ 4 972 ∂ 2 f ∂ t 2 i , j , k n + 2 3 .

The internal nodes of the ∂ 2 u ∂ z 2 i , j , k n + 2 3 can also be computed by the fourth-order Padé formula [30]:

(23) 1 12 ∂ 2 u ∂ z 2 i , j , k + 1 n + 2 3 + 5 6 ∂ 2 u ∂ z 2 i , j , k n + 2 3 + 1 12 ∂ 2 u ∂ z 2 i , j , k − 1 n + 2 3 = u i , j , k + 1 n + 2 3 − 2 u i , j , k n + 2 3 + u i , j , k − 1 n + 2 3 h 2 + O ( h 4 ) , k = 1 , … , N − 1 , i , j = 0 , … N .

The boundaries of the ∂ 2 u ∂ z 2 i , j , k n + 2 3 can be obtained by equations (3) and (21):

(24) ∂ 2 u ∂ z 2 i , j , 0 n + 2 3 = 1 3 v i , j , 0 2 ∂ 2 u ∂ t 2 − f i , j , 0 n + 2 3 = 1 3 v i , j , 0 2 ∂ 2 g ∂ t 2 − f i , j , 0 n + 2 3 , i , j = 0 , … N ,

(25) ∂ 2 u ∂ z 2 i , j , N n + 2 3 = 1 3 v i , j , N 2 ∂ 2 u ∂ t 2 − f i , j , N n + 2 3 = 1 3 v i , j , N 2 ∂ 2 g ∂ t 2 − f i , j , N n + 2 3 , i , j = 0 , … N .

We can adopt equations (12)–(15), (17)–(20), and (22)–(25) to solve equations (4)–(6), respectively. According to the derivation process, it is an HOC scheme based on the LOD method, which can be denoted as HOC-LOD scheme. Moreover, the HOC-LOD method is a four-order accuracy five-step scheme. Therefore, two start-up time steps need to be computed at t = 2 τ 3 and t = τ . To approximate the start-up time steps u i , j , k 2 3 and u i , j , k 1 , by Taylor expansions, we have

(26) u i , j , k 2 3 = u i , j , k 0 + 2 τ 3 ∂ u ∂ t i , j , k 0 + 2 τ 2 9 ∂ 2 u ∂ t 2 i , j , k 0 + 4 τ 3 81 ∂ 3 u ∂ t 3 i , j , k 0 + 2 τ 4 243 ∂ 4 u ∂ t 4 i , j , k 0 + O ( τ 5 ) ,

(27) u i , j , k 1 = u i , j , k 0 + τ ∂ u ∂ t i , j , k 0 + τ 2 2 ∂ 2 u ∂ t 2 i , j , k 0 + τ 3 6 ∂ 3 u ∂ t 3 i , j , k 0 + τ 4 24 ∂ 4 u ∂ t 4 i , j , k 0 + O ( τ 5 ) .

For simplicity, we denote the time-independent functions as the following notations:

∂ ζ ∂ x = ζ x , ∂ ζ ∂ y = ζ y , ∂ ζ ∂ z = ζ z , ∂ 2 ζ ∂ x 2 = ζ x x , ∂ 2 ζ ∂ y 2 = ζ y y , ∂ 2 ζ ∂ z 2 = ζ z z , ∂ 3 ζ ∂ x 3 = ζ x x x , ∂ 3 ζ ∂ x ∂ y 2 = ζ x y y , ∂ 3 ζ ∂ x ∂ z 2 = ζ x z z , ∂ 3 ζ ∂ y ∂ x 2 = ζ x x y , ∂ 3 ζ ∂ z ∂ x 2 = ζ x x z , ∂ 3 ζ ∂ y 3 = ζ y y y , ∂ 3 ζ ∂ y ∂ z 2 = ζ y z z , ∂ 3 ζ ∂ z ∂ y 2 = ζ y y z , ∂ 3 ζ ∂ z 3 = ζ z z z , ∂ 4 ζ ∂ x 4 = ζ x x x x , ∂ 4 ζ ∂ y 4 = ζ y y y y , ∂ 4 ζ ∂ z 4 = ζ z z z z , ∂ 4 ζ ∂ x 2 ∂ y 2 = ζ x x y y , ∂ 4 ζ ∂ x 2 ∂ z 2 = ζ x x z z , ∂ 4 ζ ∂ y 2 ∂ z 2 = ζ y y z z ,

where ζ represents φ , ψ , and v 2 .

By equations (1) and (2), we can obtain

(28) ∂ 2 u ∂ t 2 i , j , k 0 = v i , j , k 2 ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 + ∂ 2 u ∂ z 2 i , j , k 0 + f i , j , k 0 = v i , j , k 2 ( φ x x + φ y y + φ z z ) i , j , k + f i , j , k 0 ,

(29) ∂ 3 u ∂ t 3 i , j , k 0 = ∂ ∂ t ∂ 2 u ∂ t 2 i , j , k 0 = ∂ ∂ t v 2 ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 + ∂ 2 u ∂ z 2 + f i , j , k 0 = v i , j , k 2 ∂ 2 ∂ x 2 + ∂ 2 ∂ y 2 + ∂ 2 ∂ z 2 ∂ u ∂ t i , j , k 0 + ∂ f ∂ t i , j , k 0 = v i , j , k 2 ( ψ x x + ψ y y + ψ z z ) i , j , k + ∂ f ∂ t i , j , k 0 ,

(30) ∂ 4 u ∂ t 4 i , j , k 0 = ∂ 2 ∂ t 2 ∂ 2 u ∂ t 2 i , j , k 0 = ∂ 2 ∂ t 2 v 2 ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 + ∂ 2 u ∂ z 2 + f i , j , k 0 = v i , j , k 2 ( v x x 2 + v y y 2 + v z z 2 ) i , j , k ( φ x x + φ y y + φ z z ) i , j , k + 2 v i , j , k 2 [ ( v x 2 ) i , j , k ( φ x x x + φ x y y + φ x z z ) i , j , k + ( v y 2 ) i , j , k ( φ y y y + φ x x y + φ y z z ) i , j , k + ( v z 2 ) i , j , k ( φ z z z + φ x x z + φ y y z ) i , j , k ] + v i , j , k 4 ( φ x x x x + φ y y y y + φ z z z z + 2 φ x x y y + 2 φ x x z z + 2 φ y y z z ) i , j , k + v i , j , k 2 ∂ 2 f ∂ x 2 + ∂ 2 f ∂ y 2 + ∂ 2 f ∂ z 2 i , j , k 0 + ∂ 2 f ∂ t 2 i , j , k 0 .

Substituting equations (28)–(30) into equation (26) and omitting the truncation error, we obtain

(31) u i , j , k 2 3 = φ i , j , k + 2 τ 3 ψ i , j , k + 2 τ 2 9 v i , j , k 2 ( φ x x + φ y y + φ z z ) i , j , k + 4 τ 3 81 v i , j , k 2 ( ψ x x + ψ y y + ψ z z ) i , j , k + 2 τ 4 243 v i , j , k 2 [ ( v x x 2 + v y y 2 + v z z 2 ) i , j , k ( φ x x + φ y y + φ z z ) i , j , k + 2 ( v x 2 ) i , j , k ( φ x x x + φ x y y + φ x z z ) i , j , k + 2 ( v y 2 ) i , j , k ( φ y y y + φ x x y + φ y z z ) i , j , k + 2 ( v z 2 ) i , j , k ( φ z z z + φ x x z + φ y y z ) i , j , k ] + 2 τ 4 243 [ v i , j , k 4 ( φ x x x x + φ y y y y + φ z z z z + 2 φ x x y y + 2 φ x x z z + 2 φ y y z z ) i , j , k ] + 2 τ 2 9 f i , j , k 0 + 4 τ 3 81 ∂ f ∂ t i , j , k 0 + 2 τ 4 243 v i , j , k 2 ∂ 2 f ∂ x 2 + ∂ 2 f ∂ y 2 + ∂ 2 f ∂ z 2 i , j , k 0 + ∂ 2 f ∂ t 2 i , j , k 0 .

Similarly, substituting equations (28)–(30) into equation (27), we have

(32) u i , j , k 1 = φ i , j , k + τ ψ i , j , k + τ 2 2 v i , j , k 2 ( φ x x + φ y y + φ z z ) i , j , k + τ 3 6 v i , j , k 2 ( ψ x x + ψ y y + ψ z z ) i , j , k + τ 4 24 v i , j , k 2 [ ( v x x 2 + v y y 2 + v z z 2 ) i , j , k ( φ x x + φ y y + φ z z ) i , j , k + 2 ( v x 2 ) i , j , k ( φ x x x + φ x y y + φ x z z ) i , j , k + 2 ( v y 2 ) i , j , k ( φ y y y + φ x x y + φ y z z ) i , j , k + 2 ( v z 2 ) i , j , k ( φ z z z + φ x x z + φ y y z ) i , j , k ] + τ 4 24 [ v i , j , k 4 ( φ x x x x + φ y y y y + φ z z z z + 2 φ x x y y + 2 φ x x z z + 2 φ y y z z ) i , j , k ] + τ 2 2 f i , j , k 0 + τ 3 6 ∂ f ∂ t i , j , k 0 + τ 4 24 v i , j , k 2 ∂ 2 f ∂ x 2 + ∂ 2 f ∂ y 2 + ∂ 2 f ∂ z 2 i , j , k 0 + ∂ 2 f ∂ t 2 i , j , k 0 .

Based on the aforementioned description, an algorithm for solving equations (4)–(6) with equations (2)–(3) by the HOC-LOD scheme is given as follows:

Step 1: Letting n = 1 and giving the initial time step u i , j , k 0 , the start-up time steps u i , j , k 2 3 and u i , j , k 1 are computed by equations (31) and (32).
Step 2: The values of ∂ 2 u ∂ x 2 i , j , k n are computed by equations (13)–(15). Then, the function values u i , j , k n + 1 3 are computed explicitly by equation (12).
Step 3: The values of ∂ 2 u ∂ y 2 i , j , k n + 1 3 are computed by equations (18)–(20). Then, the function values u i , j , k n + 2 3 are computed explicitly by equation (17).
Step 4: The values of ∂ 2 u ∂ z 2 i , j , k n + 2 3 are computed by equations (23)–(25). Then, the function values u i , j , k n + 1 are computed explicitly by equation (22).
Step 5: Let n ← n + 1 , repeat Steps 2–4 until time reaches the final moment, and the computation is terminated.

3 Theory analysis

3.1 Consistency analysis

Theorem 1

The HOC-LOD scheme in Section 2 is consistent with equation (1) under three cases as follows:

Case I: when ∂ 3 u ∂ t ∂ x 2 − ∂ 3 u ∂ t ∂ z 2 ≠ 0 , the truncation error of the scheme is O τ + τ 2 h 2 + τ 3 h + τ 4 + τ 5 h + τ 6 h 2 .
Case II: when ∂ 3 u ∂ t ∂ x 2 − ∂ 3 u ∂ t ∂ z 2 = 0 5 ∂ 4 u ∂ t 2 ∂ x 2 + 2 ∂ 4 u ∂ t 2 ∂ y 2 − 7 ∂ 4 u ∂ t 2 ∂ z 2 ≠ 0 , the truncation error of the scheme is O τ 2 + τ 2 h 2 + τ 3 h + τ 4 + τ 5 h + τ 6 h 2 .
Case III: when ∂ 3 u ∂ t ∂ x 2 − ∂ 3 u ∂ t ∂ z 2 = 0 5 ∂ 4 u ∂ t 2 ∂ x 2 + 2 ∂ 4 u ∂ t 2 ∂ y 2 − 7 ∂ 4 u ∂ t 2 ∂ z 2 = 0 , the truncation error of the scheme is O τ 2 h 2 + τ 3 h + τ 4 + τ 5 h + τ 6 h 2 .

Proof

First, for equation (12) using Taylor expansions of u i , j , k n ± 1 3 , v 2 ∂ 2 u ∂ x 2 i ± 1 , j , k n , and f i ± 1 , j , k n at point ( x i , y j , z k , t n ) , we can obtain

(33) u i , j , k n ± 1 3 = u i , j , k n ± τ 3 ∂ u ∂ t i , j , k n + τ 2 18 ∂ 2 u ∂ t 2 i , j , k n ± τ 3 162 ∂ 3 u ∂ t 3 i , j , k n + τ 4 1944 ∂ 4 u ∂ t 4 i , j , k n ± τ 5 29160 ∂ 5 u ∂ t 5 i , j , k n + O ( τ 6 ) ,

(34) v 2 ∂ 2 u ∂ x 2 i ± 1 , j , k n = v 2 ∂ 2 u ∂ x 2 i , j , k n ± h ∂ v 2 ∂ x ∂ 2 u ∂ x 2 + v 2 ∂ 3 u ∂ x 3 i , j , k n + h 2 2 ∂ 2 v 2 ∂ x 2 ∂ 2 u ∂ x 2 + 2 ∂ v 2 ∂ x ∂ 3 u ∂ x 3 + v 2 ∂ 4 u ∂ x 4 i , j , k n ± h 3 6 ∂ 3 v 2 ∂ x 3 ∂ 2 u ∂ x 2 + 3 ∂ 2 v 2 ∂ x 2 ∂ 3 u ∂ x 3 + 3 ∂ v 2 ∂ x ∂ 4 u ∂ x 4 + v 2 ∂ 5 u ∂ x 5 i , j , k n + O ( h 4 ) ,

(35) f i ± 1 , j , k n = f i , j , k n ± h ∂ f ∂ x i , j , k n + h 2 ∂ 2 f ∂ x 2 i , j , k n ± h 2 6 ∂ 3 f ∂ x 3 i , j , k n + O ( h 4 ) .

Substituting equations (33)–(35) into equation (12), we obtain

(36) τ 2 9 ∂ 2 u ∂ t 2 i , j , k n + τ 4 972 ∂ 4 u ∂ t 4 i , j , k n = v i , j , k 2 τ 4 108 h 2 2 v 2 ∂ 2 u ∂ x 2 i , j , k n + h 2 ∂ 2 v 2 ∂ x 2 ∂ 2 u ∂ x 2 + 2 ∂ v 2 ∂ x ∂ 3 u ∂ x 3 + v 2 ∂ 4 u ∂ x 4 i , j , k n + v i , j , k 2 τ 2 3 − v i , j , k 4 τ 4 54 h 2 ∂ 2 u ∂ x 2 i , j , k n + v i , j , k 2 τ 4 324 h 2 2 f i , j , k n + h 2 ∂ 2 f ∂ x 2 i , j , k n + τ 2 9 − v i , j , k 2 τ 4 162 h 2 f i , j , k n + τ 4 972 ∂ 2 f ∂ x 2 i , j , k n + O ( τ 6 + τ 4 h 2 ) .

By using the relationship of equation (4), equation (36) is equivalent to

(37) τ 2 9 ∂ 2 u ∂ t 2 i , j , k n + v i , j , k 2 τ 4 324 ∂ 2 ∂ x 2 3 v 2 ∂ 2 u ∂ x 2 + f i , j , k n + τ 4 972 ∂ 2 f ∂ t 2 i , j , k n = v i , j , k 2 τ 2 3 ∂ 2 u ∂ x 2 i , j , k n + v i , j , k 2 τ 4 108 ∂ 2 v 2 ∂ x 2 ∂ 2 u ∂ x 2 + 2 ∂ v 2 ∂ x ∂ 3 u ∂ x 3 + v 2 ∂ 4 u ∂ x 4 i , j , k n + τ 2 9 f i , j , k n + v i , j , k 2 τ 4 324 ∂ 2 f ∂ x 2 i , j , k n + τ 4 972 ∂ 2 f ∂ x 2 i , j , k n + O ( τ 6 + τ 4 h 2 ) .

It can be rewritten as follows:

(38) τ 2 9 ∂ 2 u ∂ t 2 i , j , k n = v i , j , k 2 τ 2 3 ∂ 2 u ∂ x 2 i , j , k n + τ 2 9 f i , j , k n + O ( τ 6 + τ 4 h 2 ) .

Multiplying by 3 τ 2 on the both sides of equation (38), i.e.,

(39) 1 3 ∂ 2 u ∂ t 2 i , j , k n = v i , j , k 2 ∂ 2 u ∂ x 2 i , j , k n + 1 3 f i , j , k n + O ( τ 4 + τ 2 h 2 ) .

Second, for equation (17) using the Taylor expansions of u i , j , k n + 2 3 , u i , j , k n + 1 3 , v 2 ∂ 2 u ∂ y 2 i , j ± 1 , k n + 1 3 , v 2 ∂ 2 u ∂ y 2 i , j , k n + 1 3 , f i , j ± 1 , k n + 1 3 , f i , j , k n + 1 3 , and ∂ 2 f ∂ t 2 i , j , k n + 1 3 at point ( x i , y j , z k , t n ) , we can obtain

(40) u i , j , k n + 2 3 = u i , j , k n + 2 τ 3 ∂ u ∂ t i , j , k n + 2 τ 2 9 ∂ 2 u ∂ t 2 i , j , k n + 4 τ 3 81 ∂ 3 u ∂ t 3 i , j , k n + 2 τ 4 243 ∂ 4 u ∂ t 4 i , j , k n + 4 τ 5 3654 ∂ 5 u ∂ t 5 i , j , k n + O ( τ 6 ) ,

(41) v 2 ∂ 2 u ∂ y 2 i , j , k n + 1 3 = v 2 ∂ 2 u ∂ y 2 i , j , k n + τ 3 v 2 ∂ 3 u ∂ t ∂ y 2 i , j , k n + τ 2 18 v 2 ∂ 4 u ∂ t 2 ∂ y 2 i , j , k n + τ 3 162 v 2 ∂ 5 u ∂ t 3 ∂ y 2 i , j , k n + O ( τ 4 ) ,

(42) v 2 ∂ 2 u ∂ y 2 i , j ± 1 , k n + 1 3 = v 2 ∂ 2 u ∂ y 2 i , j , k n ± h ∂ v 2 ∂ y ∂ 2 u ∂ y 2 + v 2 ∂ 3 u ∂ y 3 i , j , k n + h 2 2 ∂ 2 v 2 ∂ y 2 ∂ 2 u ∂ y 2 + 2 ∂ v 2 ∂ y ∂ 3 u ∂ y 3 + v 2 ∂ 4 u ∂ y 4 i , j , k n ± h 3 6 ∂ 3 v 2 ∂ y 3 ∂ 2 u ∂ y 2 + 3 ∂ 2 v 2 ∂ y 2 ∂ 3 u ∂ y 3 + 3 ∂ v 2 ∂ y ∂ 4 u ∂ y 4 + v 2 ∂ 5 u ∂ y 5 i , j , k n + τ 3 v 2 ∂ 3 u ∂ t ∂ y 2 i , j , k n ± τ h 3 ∂ v 2 ∂ y ∂ 3 u ∂ t ∂ y 2 + v 2 ∂ 4 u ∂ t ∂ y 3 i , j , k n + τ h 2 6 ∂ 2 v 2 ∂ y 2 ∂ 3 u ∂ t ∂ y 2 + 2 ∂ v 2 ∂ y ∂ 4 u ∂ t ∂ y 3 + v 2 ∂ 5 u ∂ t ∂ y 4 i , j , k n + τ 2 18 v 2 ∂ 4 u ∂ t 2 ∂ y 2 i , j , k n ± τ 2 h 18 ∂ v 2 ∂ y ∂ 4 u ∂ t 2 ∂ y 2 + v 2 ∂ 5 u ∂ t 2 ∂ y 3 i , j , k n + τ 3 162 v 2 ∂ 5 u ∂ t 3 ∂ y 2 i , j , k n + O ( τ 4 + h 4 + τ h 3 + τ 2 h 2 + τ 3 h ) ,

(43) f i , j ± 1 , k n + 1 3 = f i , j , k n ± h ∂ f ∂ x i , j , k n + h 2 2 ∂ 2 f ∂ x 2 i , j , k n ± h 3 6 ∂ 3 f ∂ x 3 i , j , k n + τ 3 ∂ f ∂ t i , j , k n ± τ h 3 ∂ 2 f ∂ t ∂ x i , j , k n + τ h 2 6 ∂ 3 f ∂ t ∂ x 2 i , j , k n + τ 2 18 ∂ 2 f ∂ t 2 i , j , k n ± τ 2 h 18 ∂ 3 f ∂ t 2 ∂ x i , j , k n + τ 3 162 ∂ 3 f ∂ t 3 i , j , k n + O ( τ 4 + h 4 + τ h 3 + τ 2 h 2 + τ 3 h ) ,

(44) f i , j , k n + 1 3 = f i , j , k n + τ 3 ∂ f ∂ t i , j , k n + τ 2 18 ∂ 2 f ∂ t 2 i , j , k n + τ 3 162 ∂ 3 f ∂ t 3 i , j , k n + O ( τ 4 ) ,

(45) ∂ 2 f ∂ t 2 i , j , k n + 1 3 = ∂ 2 f ∂ t 2 i , j , k n + τ 3 ∂ 3 f ∂ t 3 i , j , k n + O ( τ 2 ) .

Substituting equations (40)–(45) and equation (33) into equation (17), we obtain

τ 2 9 ∂ 2 u ∂ t 2 i , j , k n + τ 3 27 ∂ 3 u ∂ t 3 i , j , k n + 7 τ 4 972 ∂ 4 u ∂ t 4 i , j , k n + 4 τ 5 972 ∂ 5 u ∂ t 5 i , j , k n = v i , j , k 2 τ 4 108 h 2 2 v 2 ∂ 2 u ∂ y 2 + 2 τ 3 v 2 ∂ 3 u ∂ t ∂ y 2 + τ 2 9 v 2 ∂ 4 u ∂ t 2 ∂ y 2 + τ 3 9 v 2 ∂ 5 u ∂ t 3 ∂ y 2 + h 2 ∂ 2 v 2 ∂ y 2 ∂ 2 u ∂ y 2 + 2 ∂ v 2 ∂ y ∂ 3 u ∂ y 3 + v 2 ∂ 4 u ∂ y 4 + τ h 2 3 ∂ 2 v 2 ∂ y 2 ∂ 3 u ∂ t ∂ y 2 + 2 ∂ v 2 ∂ y ∂ 4 u ∂ t ∂ y 3 + v 2 ∂ 5 u ∂ t ∂ y 4 i , j , k n

(46) + v i , j , k 2 τ 2 3 − v i , j , k 4 τ 4 54 h 2 ∂ 2 u ∂ y 2 + τ 3 ∂ 3 u ∂ t ∂ y 2 + τ 2 18 ∂ 4 u ∂ t 2 ∂ y 2 + τ 3 162 ∂ 5 u ∂ t 3 ∂ y 2 i , j , k n + v i , j , k 2 τ 4 324 h 2 2 f i , j , k n + h 2 ∂ 2 f ∂ x 2 i , j , k n + 2 τ 3 ∂ f ∂ t i , j , k n + τ h 2 3 ∂ 3 f ∂ t ∂ x 2 i , j , k n + τ 2 9 ∂ 2 f ∂ t 2 i , j , k n + τ 3 81 ∂ 3 f ∂ t 3 i , j , k n + τ 2 9 − v i , j , k 2 τ 4 162 h 2 f i , j , k n + τ 3 ∂ f ∂ t i , j , k n + τ 2 18 ∂ 2 f ∂ t 2 i , j , k n + τ 3 162 ∂ 3 f ∂ t 3 i , j , k n + τ 4 972 ∂ 2 f ∂ t 2 i , j , k n + τ 3 ∂ 3 f ∂ t 3 i , j , k n + O τ 6 + τ 4 h 2 + τ 5 h + τ 7 h + τ 8 h 2 .

By using the relationship of equation (5), equation (46) is equivalent to

(47) τ 2 9 ∂ 2 u ∂ t 2 i , j , k n + v i , j , k 2 τ 3 9 ∂ 3 u ∂ t ∂ y 2 i , j , k n + 7 v i , j , k 2 τ 4 324 ∂ 2 ∂ y 2 3 v 2 ∂ 2 u ∂ y 2 + f i , j , k n + v i , j , k 2 τ 5 324 ∂ 3 ∂ t ∂ y 2 3 v 2 ∂ 2 u ∂ y 2 + f i , j , k n + τ 3 27 ∂ f ∂ t i , j , k n + 7 τ 4 972 ∂ 2 f ∂ t 2 i , j , k n + τ 5 972 ∂ 3 f ∂ t 3 i , j , k n = τ 2 3 v 2 ∂ 2 u ∂ y 2 i , j , k n + v i , j , k 2 τ 3 9 ∂ 3 u ∂ t ∂ y 2 i , j , k n + 7 v i , j , k 2 τ 4 108 + v i , j , k 2 τ 5 108 ∂ ∂ t ∂ 2 v 2 ∂ y 2 ∂ 2 u ∂ y 2 + 2 ∂ v 2 ∂ y ∂ 3 u ∂ y 3 + v 2 ∂ 4 u ∂ y 4 i , j , k n + τ 2 9 f i , j , k n + τ 3 27 ∂ f ∂ t i , j , k n + 7 τ 4 972 ∂ 2 f ∂ t 2 i , j , k n + τ 5 972 ∂ 3 f ∂ t 3 i , j , k n + 7 v i , j , k 2 τ 4 324 ∂ 2 f ∂ y 2 i , j , k n + v i , j , k 2 τ 5 324 ∂ 3 f ∂ t ∂ y 2 i , j , k n + O τ 6 + τ 4 h 2 + τ 5 h + τ 7 h + τ 8 h 2 .

It can be rewritten as follows:

(48) τ 2 9 ∂ 2 u ∂ t 2 i , j , k n = τ 2 3 v 2 ∂ 2 u ∂ y 2 i , j , k n + τ 2 9 f i , j , k n + O τ 6 + τ 4 h 2 + τ 5 h + τ 7 h + τ 8 h 2 .

Multiplying by 3 τ 2 on the both sides of equation (48), i.e.,

(49) 1 3 ∂ 2 u ∂ t 2 i , j , k n = v i , j , k 2 ∂ 2 u ∂ y 2 i , j , k n + 1 3 f i , j , k n + O τ 4 + τ 2 h 2 + τ 3 h + τ 5 h + τ 6 h 2 .

Similarly, we can obtain

(50) 1 3 ∂ 2 u ∂ t 2 i , j , k n = v i , j , k 2 ∂ 2 u ∂ z 2 i , j , k n + 1 3 f i , j , k n + O τ 4 + τ 2 h 2 + τ 3 h + τ 5 h + τ 6 h 2 .

Combining equation (39) with equations (49) and (50), noticing the splitting error between equation (1) and equations (4)–(6) (Appendix A). We obtain the following conclusions about the consistency of the HOC-LOD scheme:

Case I: when ∂ 3 u ∂ t ∂ x 2 − ∂ 3 u ∂ t ∂ z 2 ≠ 0 , the truncation error of the scheme is O τ + τ 2 h 2 + τ 3 h + τ 4 + τ 5 h + τ 6 h 2 .
Case II: when ∂ 3 u ∂ t ∂ x 2 − ∂ 3 u ∂ t ∂ z 2 = 0 5 ∂ 4 u ∂ t 2 ∂ x 2 + 2 ∂ 4 u ∂ t 2 ∂ y 2 − 7 ∂ 4 u ∂ t 2 ∂ z 2 ≠ 0 , the truncation error of the scheme is O τ 2 + τ 2 h 2 + τ 3 h + τ 4 + τ 5 h + τ 6 h 2 .
Case III: when ∂ 3 u ∂ t ∂ x 2 − ∂ 3 u ∂ t ∂ z 2 = 0 5 ∂ 4 u ∂ t 2 ∂ x 2 + 2 ∂ 4 u ∂ t 2 ∂ y 2 − 7 ∂ 4 u ∂ t 2 ∂ z 2 = 0 , the truncation error of the scheme is O τ 2 h 2 + τ 3 h + τ 4 + τ 5 h + τ 6 h 2 .□

3.2 Stability analysis

Lemma 1

[31] The sufficient and necessary condition for the roots of the quadratic equation μ 2 − b μ − c = 0 with real coefficients to be not greater than one is ∣ c ∣ ≤ 1 , ∣ b ∣ ≤ 1 − c .

Theorem 2

The HOC-LOD method is stable if

(51) max 1 ≤ i , j , k ≤ N v i , j , k ⋅ τ h = v max λ ≤ 0.7290 ,

where v max = max 1 ≤ i , j , k ≤ N ∣ v i , j , k ∣ .

Proof

Letting u i , j , k n = ξ n e I σ 1 x i e I σ 2 y j e I σ 3 z k , ∂ 2 u ∂ x 2 i , j , k n = η n e I σ 1 x i e I σ 2 y j e I σ 3 z k , ∂ 2 u ∂ y 2 i , j , k n = γ n e I σ 1 x i e I σ 2 y j e I σ 3 z k , ∂ 2 u ∂ z 2 i , j , k n = ς n e I σ 1 x i e I σ 2 y j e I σ 3 z k , where η n , ς n , γ n , and ξ n are amplitudes at the ( n ) t h time step; σ 1 , σ 2 , and σ 3 are wave number in x -, y - and z -directions, respectively; and I = − 1 is the imaginary unit. Equations (13), (18), and (23) can be concluded that

(52) η n e I σ 1 x i e I σ 2 y j e I σ 3 z k ( e I σ 1 h + 10 + e − I σ 1 h ) = 12 h 2 ξ n e I σ 1 x i e I σ 2 y j e I σ 3 z k ( e I σ 1 h − 2 + e − I σ 1 h ) ,

(53) γ n + 1 3 e I σ 1 x i e I σ 2 y j e I σ 3 z k ( e I σ 2 h + 10 + e − I σ 2 h ) = 12 h 2 ξ n + 1 3 e I σ 1 x i e I σ 2 y j e I σ 3 z k ( e I σ 2 h − 2 + e − I σ 2 h ) ,

(54) ς n + 2 3 e I σ 1 x i e I σ 2 y j e I σ 3 z k ( e I σ 3 h + 10 + e − I σ 3 h ) = 12 h 2 ξ n + 2 3 e I σ 1 x i e I σ 2 y j e I σ 3 z k ( e I σ 3 h − 2 + e − I σ 3 h ) .

By Euler formula, e I σ h = cos σ h + I sin σ h , e − I σ h = cos σ h − I sin σ h , and equations (52)–(54) turn to

(55) η n ξ n = 12 ( cos σ 1 h − 1 ) h 2 ( cos σ 1 h + 5 ) ,

(56) γ n + 1 3 ξ n + 1 3 = 12 ( cos σ 2 h − 1 ) h 2 ( cos σ 2 h + 5 ) ,

(57) ς n + 2 3 ξ n + 2 3 = 12 ( cos σ 3 h − 1 ) h 2 ( cos σ 3 h + 5 ) .

Letting v i , j , k n + 1 3 = u i , j , k n , ( v i , j , k 2 ) max = a , we assume that f ( x , y , z , t ) is exact and does not affect the stability. Thus, the homogeneous form of equation (12) can be written in the matrix form as follows:

(58) u i , j , k n + 1 3 v i , j , k n + 1 3 = 2 − 1 1 0 u i , j , k n v i , j , k n + a τ 2 3 − a 2 τ 2 λ 2 54 0 0 0 ( u x x ) i , j , k n ( v x x ) i , j , k n + a 2 τ 2 λ 2 108 0 0 0 ( u x x ) i + 1 , j , k n ( v x x ) i + 1 , j , k n + a 2 τ 2 λ 2 108 0 0 0 ( u x x ) i − 1 , j , k n ( v x x ) i − 1 , j , k n .

Substituting [ u i , j , k n , v i , j , k n ] T = ξ n e I σ 1 x i e I σ 2 y j e I σ 3 z k and [ ( u x x ) i , j , k n , ( v x x ) i , j , k n ] T = η n e I σ 1 x i e I σ 2 y j e I σ 3 z k into the aforementioned formula, we obtain

(59) ξ n + 1 3 = 2 − 1 1 0 ξ n + a τ 2 3 + a 2 τ 2 λ 2 ( cos σ 1 − 1 ) 54 0 0 0 η n .

Substituting equation (55) into equation (59), we have

(60) ξ n + 1 3 = 2 − 1 1 0 ξ n + a τ 2 3 + a 2 τ 2 λ 2 ( cos σ 1 h − 1 ) 54 0 0 0 12 ( cos σ 1 h − 1 ) h 2 ( cos σ 1 h + 5 ) ξ n .

Similarly, equation (17) can be treated as follows:

(61) ξ n + 2 3 = 2 − 1 1 0 ξ n + 1 3 + a τ 2 3 + a 2 τ 2 λ 2 ( cos σ 2 h − 1 ) 54 0 0 0 12 ( cos σ 2 h − 1 ) h 2 ( cos σ 2 h + 5 ) ξ n + 1 3 ,

Similarly, equation (22) can be treated as follows:

(62) ξ n + 1 = 2 − 1 1 0 ξ n + 2 3 + a τ 2 3 + a 2 τ 2 λ 2 ( cos σ 3 h − 1 ) 54 0 0 0 12 ( cos σ 3 h − 1 ) h 2 ( cos σ 3 h + 5 ) ξ n + 2 3 .

Substituting equations (60)–(61) into equation (62), the error propagation matrix can be obtained as follows:

(63) G = C − 1 1 0 B − 1 1 0 A − 1 1 0 = C B A − A − C 1 − C B B A − 1 − B ,

where,

(64) A = 2 + 4 a λ 2 ( cos σ 1 h − 1 ) ( cos σ 1 h + 5 ) + 2 a 2 λ 4 ( cos σ 1 h − 1 ) 2 9 ( cos σ 1 h + 5 ) , B = 2 + 4 a λ 2 ( cos σ 2 h − 1 ) ( cos σ 2 h + 5 ) + 2 a 2 λ 4 ( cos σ 2 h − 1 ) 2 9 ( cos σ 2 h + 5 ) , C = 2 + 4 a λ 2 ( cos σ 3 h − 1 ) ( cos σ 3 h + 5 ) + 2 a 2 λ 4 ( cos σ 3 h − 1 ) 2 9 ( cos σ 3 h + 5 ) .

The characteristic equation can be obtained as follows:

(65) ∣ μ I − G ∣ = μ − ( C B A − A − C ) C B − 1 1 − B A μ + B = μ 2 − ( C B A − A − B − C ) μ + 1 = 0 .

According to the Lemma 1, ∣ b ∣ = ∣ C B A − A − B − C ∣ ≤ 2 . Since A , B , and C have the same range of values, we only need to solve − 2 ≤ A , B , C ≤ − 1 or − 1 ≤ A , B , C ≤ 1 or 1 ≤ A , B , C ≤ 2 . Here, B and C are similar to A , we only analyze the value range of A .

Letting cos σ 1 h = ω , ω ∈ [ − 1 , 1 ] , we assume that A = F ( ω ) = 2 + 4 a λ 2 ( ω − 1 ) ω + 5 + 2 a 2 λ 4 ( ω − 1 ) 2 9 ( ω + 5 ) is a quadratic function of ω . So, F ′ ( ω ) = 18 a 2 λ 4 ω 2 + 180 a λ 2 ω − 198 a 2 λ 4 + 1944 a λ 2 81 ( ω + 5 ) 2 , and the two roots of F ′ ( ω ) are ω 1 = − 5 − 36 − 108 a λ 2 and ω 2 = − 5 + 36 − 108 a λ 2 .

First, when 36 − 108 a λ 2 ≤ 0 , i.e., a λ 2 ≤ 3 , ω 1 and ω 2 do not exist, F ′ ( ω ) > 0 always holds. The value range of the function F ( ω ) is [ F ( − 1 ) , F ( 1 ) ] = 2 + 2 a 2 λ 4 − 18 a λ 2 9 , 2 . When 1 ≤ F ( ω ) ≤ 2 , we obtain 2 + 2 a 2 λ 4 − 18 a λ 2 9 ≥ 1 and 2 + 2 a 2 λ 4 − 18 a λ 2 9 < 2 , i.e., a λ 2 ≤ 9 − 3 7 2 or 9 + 3 7 2 ≤ a λ 2 < 9 . Due to a λ 2 ≤ 3 , so, a λ 2 ≤ 9 − 3 7 2 ≈ 0.5314 .

Second, when ω 2 ≤ − 1 , i.e., 3 ≤ a λ 2 ≤ 5.4 , the value range of the function F ( ω ) is [ F ( − 1 ) , F ( 1 ) ] = 2 + 2 a 2 λ 4 − 18 a λ 2 9 , 2 . When 1 ≤ F ( ω ) ≤ 2 , i.e., 2 + 2 a 2 λ 4 − 18 a λ 2 9 ≥ 1 and 2 + 2 a 2 λ 4 − 18 a λ 2 9 < 2 , we have a λ 2 ≤ 9 − 3 7 2 or 9 + 3 7 2 ≤ a λ 2 < 9 . Because 3 ≤ a λ 2 ≤ 5.4 , there is no solution.

Third, when − 1 ≤ ω 2 ≤ 1 , i.e., 5.4 ≤ a λ 2 ≤ 9 , the value range of the function F ( ω ) is [ F ( ω 2 ) , F ( 1 ) ] = F − 5 + 36 − 108 a λ 2 , 2 . When 1 ≤ F ( ω ) ≤ 2 , we have F − 5 + 36 − 108 a λ 2 ≥ 1 . Due to 5.4 ≤ a λ 2 ≤ 9 , F − 5 + 36 − 108 a λ 2 − 1 ≤ 0 always holds. So there is no solution.

Finally, when a λ 2 > 9 , the value range of the function F ( ω ) is [ F ( ω 2 ) , F ( − 1 ) ] = F − 5 + 36 − 108 a λ 2 , 2 + 2 a 2 λ 4 − 18 a λ 2 9 . Due to F max ( ω ) = F ( − 1 ) ≥ 2 , so it does not satisfy this condition − 2 ≤ A , B , C ≤ − 1 or − 1 ≤ A , B , C ≤ 1 or 1 ≤ A , B , C ≤ 2 .

We can conclude that when 0 ≤ a λ 2 ≤ 9 − 3 7 2 , i.e., v max λ = a λ ≤ 9 − 3 7 2 ≈ 0.7290 , the HOC-LOD scheme is stable. This completes the proof.□

Table 1 presents the stability conditions of various schemes in [18,20,22,27] and the present HOC-LOD scheme for the 3D wave equations. CPD-ADI [18] means compact Padé difference ADI scheme. RE [20] stands for Richardson Extrapolation. The conventional fourth-order Runge-Kutta method is represented by RK4 [20]. Among them, some methods aim at the wave equation with constant wave velocity such as [18,27], in which, ∣ v ∣ represents the magnitude when the wave velocity is constant. The others deal with the wave equation with variable wave coefficients such as [20,22], where v max = ∥ v ∥ = max α ≤ x , y , z ≤ β ∣ v ( x , y , z ) ∣ represents the maximum value of the variable coefficients. We can conclude that the HOC-LOD scheme has better stability than other methods in the literature.

Table 1

The stability conditions of various difference schemes for the 3D wave equations

Schemes	Stability conditions
CPD-ADI [18]	∣ v ∣ λ ∈ ( 0 , 0.6079 )
LOD-1 [27]	∣ v ∣ λ ∈ ( 0 , 0.6157 )
LOD-2 [27]	∣ v ∣ λ ∈ ( 0 , 0.6160 )
RE [20]	v max λ ∈ ( 0 , 0.4714 )
ADI [22]	v max λ ∈ ( 0 , 0.5770 )
RK4 [20]	v max λ ∈ ( 0 , 0.7071 )
HOC-LOD	v max λ ∈ ( 0 , 0.7290 )

3.3 Convergence analysis

Theorem 3

The HOC-LOD scheme in Section 2 is conditionally convergent, and the convergent condition is v max λ ∈ ( 0 , 0.7290 ) .

Proof

On the one hand, from Theorem 1, we obtain the HOC-LOD scheme is consistent with equation (1). On the other hand, from Theorem 2, we obtain the HOC-LOD scheme is conditionally stable, and stable condition is v max λ ∈ ( 0 , 0.7290 ) . By the Lax equivalence theorem [32], we complete the proof.□

4 Illustrative experiments

We utilize six numerical examples to illustrate the computational efficiency and stability of the present HOC-LOD method. First four examples originate from the literature and their solutions satisfy Case III, so we can obtain fourth-order accurate numerical solutions in both time and space for these examples. In addition, we construct Examples 5 and 6 to make their solutions to satisfy Case I and Case II, respectively. For the latter two examples, the HOC-LOD scheme only has the first and second-order temporal accuracy, but the fourth-order spatial accuracy in theory. All programs are written in Fortran 90 language by using double-precision arithmetic and run on an Intel Core i5-4210U CPU@1.70GHz 2.40GHz desk computer with 4GB memory.

L 2 , L ∞ errors and the Rate of convergence are defined as follows:

L ∞ = max 0 ≤ i , j , k ≤ N ∣ u i , j , k M − u ( x i , y j , z k , t M ) ∣ , L 2 = h 3 ∑ i , j , k = 0 N [ u i , j , k M − u ( x i , y j , z k , t M ) ] 2 ,

Rate = log [ L ∞ ( h 1 ) / L ∞ ( h 2 ) ] log ( h 1 / h 2 ) or Rate = log [ L 2 ( h 1 ) / L 2 ( h 2 ) ] log ( h 1 / h 2 ) .

Example 1 [27]

∂ 2 u ∂ t 2 = ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 + ∂ 2 u ∂ z 2 , ( x , y , z ) ∈ [ 0 , 1 ] 3 , t > 0 u ( x , y , z , 0 ) = 0 , ∂ u ( x , y , z , 0 ) ∂ t = 3 π sin ( π x ) sin ( π y ) sin ( π z ) , u ( 0 , y , z , t ) = 0 , u ( 1 , y , z , t ) = 0 , u ( x , 0 , z , t ) = 0 , u ( x , 1 , z , t ) = 0 , u ( x , y , 0 , t ) = 0 , u ( x , y , 1 , t ) = 0 .

The exact solution is u ( x , y , z , t ) = sin ( 3 π t ) sin ( π x ) sin ( π y ) sin ( π z ) .

First, to test the time accuracy of the HOC-LOD scheme, we take various τ at T = 1 and h = 0.008 , which means that the approximation error in space direction is negligible. L 2 and L ∞ errors and rate are presented in Table 2. It is obvious that the HOC-LOD method can achieve fourth-order time convergence. Moreover, since the CFL condition number of the HOC-LOD scheme is 0.7290, when τ ≤ 0.7290 h ≈ 0.0058 , the HOC-LOD scheme is convergent, while τ = 0.0060 > 0.0058 , the HOC-LOD diverges because it exceeds the range of the stability condition.

Table 2

L 2 , L ∞ errors and rate at T = 1 and h = 0.008 by the HOC-LOD scheme for variable τ for Example 1

τ	L ∞	Rate	L 2	Rate
0.0060	6.9207 + 004		8.9000 + 003
0.0058	1.3903( − 8 )		4.9166( − 9 )
0.0055	1.0383( − 8 )	5.4968	3.6718( − 9 )	5.4968
0.0045	2.9193( − 9 )	6.3229	1.0324( − 9 )	6.3228
0.0035	1.1307( − 9 )	3.7742	3.9984( − 10 )	3.7745

Next, the L 2 and L ∞ errors computed by the HOC-LOD scheme are listed in Table 3. Here, we take τ = 0.2 h and T = 0.2 . To make a comparison, we also list the numerical results of the LOD-1 and LOD-2 schemes developed in [27]. The calculation errors of the method in this paper are obviously smaller than that of the LOD-1 and LOD-2 schemes. To test the calculation accuracy of the HOC-LOD scheme, we draw log–log plots for the L 2 and L ∞ errors for the LOD-1 [27], LOD-2 [27], and HOC-LOD schemes. From Figure 1, it can be seen that the slope of the line of the HOC-LOD scheme is close to 4, which represents that the HOC-LOD scheme achieves fourth-order convergence overall. Then, we use this example to confirm the stability of the HOC-LOD scheme. The L ∞ error with h = 1 / 32 for different T and τ is presented in Table 4. As time T grows, we notice that when τ ≤ 1 44 , v max λ ≤ 0.7273 , the HOC-LOD scheme is convergent. When τ ≥ 1 43 , v max λ ≥ 0.7442 , the HOC-LOD scheme is divergent. These computed results are consistent with the theoretical analysis.

Table 3

L 2 , L ∞ errors when τ = 0.2 h and T = 0.2 for Example 1

h	LOD-1 scheme [27]		LOD-2 scheme [27]		HOC-LOD scheme
	L ∞	L 2	L ∞	L 2	L ∞	L 2
1/10	1.9134( − 5 )	6.7648( − 6 )	1.9017( − 5 )	6.7237( − 6 )	4.8506( − 6 )	1.7149( − 6 )
1/20	1.1970( − 6 )	4.2321( − 7 )	1.1897( − 6 )	4.2063( − 7 )	3.0509( − 7 )	1.0787( − 7 )
1/40	7.4830( − 8 )	2.6456( − 8 )	7.4372( − 8 )	2.6294( − 8 )	1.9054( − 8 )	6.7365( − 9 )
1/60	1.4782( − 8 )	5.2262( − 9 )	1.4690( − 8 )	5.1936( − 9 )	3.7605( − 9 )	1.3295( − 9 )
1/80	4.6774( − 9 )	1.6537( − 9 )	4.6460( − 9 )	1.6426( − 9 )	1.1892( − 9 )	4.2045( − 10 )
1/100	1.9159( − 9 )	6.7736( − 10 )	1.8998( − 9 )	6.7168( − 10 )	4.8691( − 10 )	1.7215( − 10 )
1/120	9.2408( − 10 )	3.2670( − 10 )	9.1212( − 10 )	3.2248( − 10 )	2.3478( − 10 )	8.3000( − 11 )
1/140	4.9895( − 10 )	1.7639( − 10 )	4.8731( − 10 )	1.7229( − 10 )	1.2669( − 10 )	4.4793( − 11 )

$Figure 1 Log–log plots for the L ∞ {L}_{\infty } (a) and L 2 {L}_{2} (b) errors.$

Figure 1

Log–log plots for the L ∞ (a) and L 2 (b) errors.

Table 4

The L ∞ error with h = 1 / 32 by the HOC-LOD scheme for different T and τ for Example 1

τ	v max λ	T = 1	T = 5	T = 10
1/50	0.6400	1.7386( − 6 )	6.3249( − 7 )	6.7520( − 6 )
1/49	0.6531	1.9506( − 6 )	8.8928( − 7 )	6.9107( − 6 )
1/48	0.6667	2.1853( − 6 )	1.1727( − 6 )	7.0885( − 6 )
1/47	0.6809	2.4458( − 6 )	1.4863( − 6 )	7.2880( − 6 )
1/46	0.6957	2.7355( − 6 )	1.8341( − 6 )	7.5121( − 6 )
1/45	0.7111	3.0585( − 6 )	2.2206( − 6 )	7.7644( − 6 )
1/44	0.7273	3.4195( − 6 )	2.6514( − 6 )	8.0491( − 6 )
1/43	0.7442	3.8240( − 6 )	1.8841 + 007	5.1103 + 031
1/42	0.7619	4.2787( − 6 )	3.4990 + 017	3.6057 + 052

Example 2 [22]

∂ 2 u ∂ t 2 = 1 + x π 2 + y π 2 + z π 2 ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 + ∂ 2 u ∂ z 2 + 2 + 3 x π 2 + 3 y π 2 + 3 z π 2 cos ( t ) sin ( x ) sin ( y ) sin ( z ) , ( x , y , z ) ∈ [ 0 , π ] 3 , t > 0 u ( x , y , z , 0 ) = sin ( x ) sin ( y ) sin ( z ) , ∂ u ( x , y , z , 0 ) ∂ t = 0 , u ( 0 , y , z , t ) = 0 , u ( π , y , z , t ) = 0 , u ( x , 0 , z , t ) = 0 , u ( x , π , z , t ) = 0 , u ( x , y , 0 , t ) = 0 , u ( x , y , π , t ) = 0 .

The exact solution is

u ( x , y , z , t ) = cos ( t ) sin ( x ) sin ( y ) sin ( z ) .

Example 2 is a 3D nonhomogeneous wave equation with variable coefficients. When T = 1 , τ = 0.0025 with various h , Table 5 presents the L ∞ error and rate by the ADI scheme [22] and the HOC-LOD scheme. The results in Table 5 show that the spatial accuracy of the HOC-LOD scheme reaches four, while the accuracy of the ADI scheme [22] is lower than four. In order to confirm the time accuracy, when T = 1 with various h and τ , we compute the L ∞ error and rate in Table 6. It is obviously shown that the ADI scheme [22] has higher temporal accuracy than the HOC-LOD scheme. However, the numerical results of the HOC-LOD scheme are better than that of the ADI scheme [22].

Table 5

L ∞ error and rate when T = 1 and τ = 0.0025 for various h for Example 2

h	ADI scheme [22]		HOC-LOD scheme
	L ∞	Rate	L ∞	Rate
π / 10	4.3196( − 4 )		6.1785( − 5 )
π / 16	5.4662( − 5 )	3.5842	9.4064( − 6 )	4.0048
π / 20	2.1191( − 5 )	3.4368	3.8507( − 6 )	4.0025
π / 32	2.5748( − 6 )	3.7659	5.8719( − 7 )	4.0014
π / 40	9.9448( − 7 )	3.7659	2.4046( − 7 )	4.0010

Table 6

L ∞ error and rate when T = 1 for different h and τ for Example 2

( τ , h )	ADI scheme [22]		HOC-LOD scheme
	L ∞	Rate	L ∞	Rate
( 1 / 16 , π / 10 )	4.0768( − 4 )		5.7613( − 5 )
( 1 / 32 , π / 20 )	2.0790( − 5 )	4.2934	3.7141( − 6 )	3.9553
( 1 / 64 , π / 40 )	9.5059( − 7 )	4.4509	2.3569( − 7 )	3.9781
( 1 / 128 , π / 80 )	4.2791( − 8 )	4.4734	1.4859( − 8 )	3.9875

Example 3 [22]

∂ 2 u ∂ t 2 = [ 1 + sin 2 ( x ) + sin 2 ( y ) + sin 2 ( z ) ] ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 + ∂ 2 u ∂ z 2 + [ 4 + 3 sin 2 ( x ) + 3 sin 2 ( y ) + 3 sin 2 ( z ) ] e − t cos ( x ) cos ( y ) cos ( z ) , ( x , y , z ) ∈ [ 0 , π ] 3 , t > 0 u ( x , y , z , 0 ) = cos ( x ) cos ( y ) cos ( z ) , ∂ u ∂ t ( x , y , z , 0 ) = − cos ( x ) cos ( y ) cos ( z ) , u ( 0 , y , z , t ) = e − t cos ( y ) cos ( z ) , u ( π , y , z , t ) = − e − t cos ( y ) cos ( z ) , u ( x , 0 , z , t ) = e − t cos ( x ) cos ( z ) , u ( x , π , z , t ) = − e − t cos ( x ) cos ( z ) , u ( x , y , 0 , t ) = e − t cos ( x ) cos ( y ) , u ( x , t , π , t ) = − e − t cos ( x ) cos ( y ) .

The exact solution is

u ( x , y , z , t ) = e − t cos ( x ) cos ( y ) cos ( z ) .

We compute the L ∞ error and rate by the HOC-LOD scheme and compare it with the ADI scheme [22] in Table 7 with various h and τ at T = 1 . The calculation results and the accuracy of the HOC-LOD scheme are much better than that of the ADI scheme in [22].

Table 7

When T = 1 for different h and τ , L ∞ error and rate for Example 3

( τ , h )	ADI scheme [22]		HOC-LOD scheme
	L ∞	Rate	L ∞	Rate
( 1 / 20 , π / 16 )	5.1391( − 5 )		3.1773( − 7 )
( 1 / 40 , π / 32 )	4.2849( − 6 )	3.5842	2.6940( − 8 )	3.4342
( 1 / 80 , π / 64 )	3.9569( − 7 )	3.4368	1.8710( − 9 )	3.7946
( 1 / 160 , π / 128 )	2.9088( − 8 )	3.7659	1.2232( − 10 )	3.9351

For this example, v max = max ( x , y , z ) ∈ [ 0 , π ] × [ 0 , π ] × [ 0 , π ] [ 1 + sin 2 ( x ) + sin 2 ( y ) + sin 2 ( z ) ] = 2 ; therefore, the grid ratio allowed by the HOC-LOD scheme is τ h ≤ 0.7290 2 = 0.3645 . The grid ratio allowed by the ADI scheme [22] is τ h < 0.5770 2 ≈ 0.2885 . First, we choose τ h = 9 10 π ≈ 0.2865 ≤ 0.2885 , v max τ h ≤ 0.5770 , the numerical solutions are tabulated in Table 8. As can be seen, when τ h = 9 10 π , both the ADI scheme [22] and the HOC-LOD scheme are convergent. Then, we choose τ h = 19 20 π ≈ 0.3024 ≥ 0.2885 , v max τ h ≥ 0.5770 , the numerical results are shown in Table 9. When τ h = 19 20 π , the ADI scheme [22] is divergent, while the HOC-LOD scheme is still convergent. Therefore, the stability condition of the scheme in this paper is superior to that of the ADI scheme [22].

Table 8

L ∞ error when T = 1 , v max τ h ≤ 0.5770 with various h and τ for Example 3

( τ , h )	ADI scheme [22]		HOC-LOD scheme
	L ∞	Rate	L ∞	Rate
( 1 / 20 , π / 18 )	3.3689( − 5 )		2.0887( − 7 )
( 1 / 40 , π / 36 )	2.7867( − 6 )	3.5842	1.7211( − 8 )	3.4837
( 1 / 60 , π / 54 )	6.9049( − 7 )	3.4410	3.6951( − 9 )	3.7945
( 1 / 80 , π / 72 )	2.7001( − 8 )	3.2638	1.1912( − 9 )	3.9351

Table 9

L ∞ error when T = 1 , v max τ h ≥ 0.5770 with various h and τ for Example 3

( τ , h )	ADI scheme [22]		HOC-LOD scheme
	L ∞	Rate	L ∞	Rate
( 1 / 20 , π / 19 )	8.5884( − 5 )		1.7655( − 7 )
( 1 / 40 , π / 38 )	1.1864( − 5 )	2.8557	1.4205( − 8 )	3.6356
( 1 / 60 , π / 57 )	7.2660(−1)	—	3.0316( − 9 )	3.8092
( 1 / 80 , π / 76 )	1.3165 + 006	—	9.7333( − 10 )	3.9492

To further verify the stability condition of the HOC-LOD scheme, Table 10 presents that the L ∞ error of the HOC-LOD scheme with h = π / 80 for different T and τ . As time T grows, we notice that when τ ≤ 1 70 , v max λ ≤ 0.7276 , the HOC-LOD scheme is convergent. When τ ≥ 1 69 , v max λ ≥ 0.7381 , the HOC-LOD scheme is divergent. The numerical results are consistent with the theoretical ones.

Table 10

When h = π / 80 for different T and τ , the L ∞ error for Example 3

τ	v max λ	T = 1	T = 5	T = 8
1/90	0.5659	7.9002( − 10 )	3.1846( − 10 )	3.6333( − 10 )
1/80	0.6366	8.0560( − 10 )	3.1785( − 10 )	3.7681( − 10 )
1/75	0.6791	8.1877( − 10 )	3.1656( − 10 )	3.8885( − 10 )
1/72	0.7074	8.2980( − 10 )	3.1711( − 10 )	3.9998( − 10 )
1/71	0.7173	8.3421( − 10 )	3.1739( − 10 )	4.0479( − 10 )
1/70	0.7276	8.3898( − 10 )	3.1765( − 10 )	4.1000( − 10 )
1/69	0.7381	8.4414( − 10 )	2.2032(-2)	5.7085 + 008
1/68	0.7490	8.4998( − 10 )	5.5407 + 017	2.0297 + 040
1/67	0.7601	8.5631( − 10 )	5.4253 + 029	2.4194 + 059

Example 4 [22]

Example 4 is a 3D domain [ 0 , 1280 m ] × [ 0 , 1280 m ] × [ 0 , 800 m ] model with a wave source. The Ricker wavelet source that generates the wave is expressed as follows:

f ( x , y , z , t ) = δ ( x − x 0 , y − y 0 , z − z 0 ) [ 1 − 2 π 2 f p 2 ( t − d r ) 2 ] e − π 2 f p 2 ( t − d r ) 2 ,

where ( x 0 , y 0 , z 0 ) = ( 640 m , 640 m , 400 m ) . d r = 0.5 / f p represents temporal delay, which can guarantee the initial conditions is zero. f p = 10 Hz is the peak frequency.

v ( x , y , z ) = 1,200 + 400 x x max 2 + 100 y y max 2 + 800 z z max 2 .

The uniform grid size h = h x = h y = h z = 10 m and τ = 0.001 s are chosen to solve Example 4.

We draw the snapshots of wave field in different times. Figure 2 describes wave field snapshots computed at z = z max / 2 by the HOC-LOD scheme. In Figure 2(a) and (b), the wave forms a complete circle before reaching the interface, which verify the accuracy of the HOC-LOD scheme. When the time is 0.4 s, the wave reaches the right side of the interface in Figure 2(c), with the increase of time, the reflection phenomenon appears near the interface. Finally, in Figure 2(d), the left and right sides of the wave are slightly impacted, resulting in reflection. The simulation results are consistent with that in the literature and the real physical process.

$Figure 2 Wave fields snapshots: (a) Z-section at t = 0.2 s t=0.2\hspace{0.33em}{\rm{s}} , (b) Z-section at t = 0.3 s t=0.3\hspace{0.33em}{\rm{s}} , (c) Z-section at t = 0.4 s t=0.4\hspace{0.33em}{\rm{s}} , and (d) Z-section at t = 0.5 s t=0.5\hspace{0.33em}{\rm{s}} .$

Figure 2

Wave fields snapshots: (a) Z-section at t = 0.2 s , (b) Z-section at t = 0.3 s , (c) Z-section at t = 0.4 s , and (d) Z-section at t = 0.5 s .

Example 5

∂ 2 u ∂ t 2 = ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 + ∂ 2 u ∂ z 2 , ( x , y , z ) ∈ [ 0 , 1 ] 3 , t > 0 , u ( x , y , z , 0 ) = sin ( π x ) + sin ( π y ) + sin ( π z ) , ∂ u ∂ t ( x , y , z , 0 ) = 0 , u ( 0 , y , z , t ) = cos ( π t ) [ sin ( π y ) + sin ( π z ) ] , u ( 1 , y , z , t ) = cos ( π t ) [ sin ( π y ) + sin ( π z ) ] , u ( x , 0 , z , t ) = cos ( π t ) [ sin ( π x ) + sin ( π z ) ] , u ( x , 1 , z , t ) = cos ( π t ) [ sin ( π x ) + sin ( π z ) ] , u ( x , y , 0 , t ) = cos ( π t ) [ sin ( π x ) + sin ( π y ) ] , u ( x , y , 1 , t ) = cos ( π t ) [ sin ( π x ) + sin ( π y ) ] .

The exact solution is

u ( x , y , z , t ) = cos ( π t ) [ sin ( π x ) + sin ( π y ) + sin ( π z ) ] .

For this example, the error of the HOC-LOD scheme is O τ + τ 2 h 2 + τ 3 h + τ 4 + τ 5 h + τ 6 h 2 . First, to test the time accuracy of the HOC-LOD scheme, we take various τ at T = 0.1 and h = 0.008 , which means that the approximation error in space direction is negligible. The L ∞ , L 2 errors and rate are listed in Table 11. Obviously, the HOC-LOD scheme has only first-order accuracy in time direction. Then, we take τ = h 4 and T = 0.001 to further test the spatial accuracy of the HOC-LOD scheme. It can be seen from the numerical results in Table 12 that the spatial accuracy of the HOC-LOD scheme is fourth-order. This is consistent with the theoretical analysis results.

Table 11

L 2 , L ∞ errors and rate at T = 0.1 and h = 0.008 by the HOC-LOD scheme for variable τ for Example 5

τ	L ∞	Rate	L 2	Rate
0.0020	4.3055( − 4 )		1.7010( − 4 )
0.0030	6.4095( − 4 )	0.9813	2.5269( − 4 )	0.9761
0.0040	8.5737( − 4 )	1.0112	3.3648( − 4 )	0.9955
0.0050	1.0735( − 3 )	1.0075	4.1824( − 4 )	0.9748
0.0060	1.2997( − 3 )	1.0487	5.0350( − 4 )	1.0176

Table 12

L 2 , L ∞ errors when τ = h 4 and T = 0.001 by the HOC-LOD scheme at variable h for Example 5

h	L ∞	Rate	L 2	Rate
1/10	2.0497( − 7 )		9.0026( − 8 )
1/15	5.0569( − 8 )	3.4517	2.2242( − 8 )	3.4482
1/20	1.7264( − 8 )	3.7358	7.5245( − 9 )	3.7674
1/25	7.3736( − 9 )	3.8124	3.2075( − 9 )	3.8212
1/30	4.0191( − 9 )	3.3284	1.5914( − 9 )	3.8442

Example 6

∂ 2 u ∂ t 2 = ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 + ∂ 2 u ∂ z 2 + π 2 cos ( π t ) sin ( π x ) sin ( π z ) , ( x , y , z ) ∈ [ 0 , 1 ] 3 , t > 0 u ( x , y , z , 0 ) = sin ( π x ) sin ( π z ) + sin ( π y ) , ∂ u ∂ t ( x , y , z , 0 ) = 0 , u ( 0 , y , z , t ) = cos ( π t ) sin ( π y ) , u ( 1 , y , t ) = cos ( π t ) sin ( π y ) , u ( x , 0 , z , t ) = cos ( π t ) sin ( π x ) sin ( π z ) , u ( x , 1 , z , t ) = cos ( π t ) sin ( π x ) sin ( π z ) , u ( x , y , 0 , t ) = cos ( π t ) sin ( π y ) , u ( x , y , 1 , t ) = cos ( π t ) sin ( π y ) .

The exact solution is

u ( x , y , z , t ) = cos ( π t ) [ sin ( π x ) sin ( π z ) + sin ( π y ) ] .

When we use the HOC-LOD scheme to solve this example, the error of the HOC-LOD scheme is O τ 2 + τ 2 h 2 + τ 3 h + τ 4 + τ 5 h + τ 6 h 2 . First, for variable τ , we take T = 0.1 and h = 0.008 to compute the L ∞ , L 2 errors and rate in Table 13. It clearly shows that the HOC-LOD scheme has second-order time accuracy. Then, in Table 14, we take τ = h 2 and T = 0.1 to test the spatial accuracy of the HOC-LOD scheme. It represents that the convergence rate of the HOC-LOD scheme is fourth-order. This is consistent with the theoretical analysis results.

Table 13

L 2 , L ∞ errors and rate at T = 0.1 and h = 0.008 by the HOC-LOD scheme for variable τ for Example 6

τ	L ∞	Rate	L 2	Rate
0.0020	4.6057( − 6 )		9.0462( − 7 )
0.0030	1.0071( − 5 )	1.9296	2.0173( − 6 )	1.9780
0.0040	1.7743( − 5 )	1.9686	3.5876( − 6 )	2.0012
0.0050	2.7165( − 5 )	1.9088	5.5812( − 6 )	1.9804
0.0060	3.9058( − 5 )	1.9916	8.0787( − 6 )	2.0284

Table 14

L 2 , L ∞ errors when τ = h 2 and T = 0.1 by the HOC-LOD scheme at variable h for Example 6

h	L ∞	Rate	L 2	Rate
1/20	6.8531( − 6 )		1.1903( − 6 )
1/30	1.3753( − 6 )	3.9610	2.5843( − 7 )	3.7669
1/40	4.3637( − 7 )	3.9903	8.4990( − 8 )	3.8657
1/50	1.7955( − 7 )	3.9797	3.5710( − 8 )	3.8859
1/60	8.8353( − 8 )	3.8894	1.7446( − 8 )	3.9289

5 Conclusion

In this paper, a novel accurate LOD method to solve the 3D acoustic wave equation is available. First, the correction of the truncation error remainder and the Padé difference approximation are employed to obtain overall fourth-order accuracy. The LOD technique and explicit time marching process are used to reduce computational cost. Moreover, the novel method is five step in time, the computations of the two initial time steps are explicitly conducted. Then, the Fourier analysis method shows that the maximum value of the CFL condition number is 0.7290, which is more flexible compared to those schemes with equivalent accuracy in the literature. The consistency and convergence of the new scheme are also proved. Finally, through several numerical examples, we can conclude that the proposed method are accurate and efficient. It needs being pointed out the HOC-LOD scheme has the fourth-order accuracy in both time and space for some special cases, while only the first- or second-order temporal accuracy but fourth-order spatial accuracy for other cases.

Recently, some authors have developed efficient and accurate difference schemes to solve the nonlinear 3D wave equations [33,34,35] and some complex fluid models [36,37,38]. To generalize the present HOC-LOD scheme to these kinds of problems is our present research interest.

Funding information: This work was partially supported by National Natural Science Foundation of China (12161067, 11772165, 11961054, and 11902170), National Natural Science Foundation of Ningxia (2022AAC02023 and 2020AAC03059), the Key Research and Development Program of Ningxia (2018BEE03007), National Youth Top-notch Talent Support Program of Ningxia, and the First Class Discipline Construction Project in Ningxia Universities: Mathematics, and Graduate Innovation Program of Ningxia University (2021043).
Conflict of interest: The authors state no conflict of interest.
Data availability statement: The data that support this study are available from the corresponding author on reasonable request.

Appendix

The LOD method is used to split the 3D wave equation (1) as follows:

(A1) 1 3 ∂ 2 u ∂ t 2 i , j , k n = v i , j , k 2 ∂ 2 u ∂ x 2 i , j , k n + 1 3 f i , j , k n ,

(A2) 1 3 ∂ 2 u ∂ t 2 i , j , k n + 1 3 = v i , j , k 2 ∂ 2 u ∂ y 2 i , j , k n + 1 3 + 1 3 f i , j , k n + 1 3 ,

(A3) 1 3 ∂ 2 u ∂ t 2 i , j , k n + 2 3 = v i , j , k 2 ∂ 2 u ∂ z 2 i , j , k n + 2 3 + 1 3 f i , j , k n + 2 3 .

First, for equation (A2), we use the Taylor expansions of ∂ 2 u ∂ t 2 i , j , k n + 1 3 , ∂ 2 u ∂ y 2 i , j , k n + 1 3 and f i , j , k n + 1 3 at point ( x i , y j , z k , t n ) , we obtain

(A4) ∂ 2 u ∂ t 2 i , j , k n + 1 3 = ∂ 2 u ∂ t 2 i , j , k n + τ 3 ∂ 3 u ∂ t 3 i , j , k n + τ 2 18 ∂ 4 u ∂ t 4 i , j , k n + τ 3 162 ∂ 5 u ∂ t 5 i , j , k n + O ( τ 4 ) ,

(A5) ∂ 2 u ∂ y 2 i , j , k n + 1 3 = ∂ 2 u ∂ y 2 i , j , k n + τ 3 ∂ 3 u ∂ t ∂ y 2 i , j , k n + τ 2 18 ∂ 4 u ∂ t 2 ∂ y 2 i , j , k n + τ 3 162 ∂ 5 u ∂ t 3 ∂ y 2 i , j , k n + O ( τ 4 ) ,

(A6) f i , j , k n + 1 3 = f i , j , k n + τ 3 ∂ f ∂ t i , j , k n + τ 2 18 ∂ 2 f ∂ t 2 i , j , k n + τ 3 162 ∂ 3 f ∂ t 3 i , j , k n + O ( τ 4 ) .

Substituting equations (A4)–(A6) into equation (A2), we obtain

(A7) 1 3 ∂ 2 u ∂ t 2 i , j , k n + τ 9 ∂ 3 u ∂ t 3 i , j , k n + τ 2 54 ∂ 4 u ∂ t 4 i , j , k n + τ 3 486 ∂ 5 u ∂ t 5 i , j , k n = v i , j , k 2 ∂ 2 u ∂ y 2 i , j , k n + τ 3 ∂ 3 u ∂ t ∂ y 2 i , j , k n + τ 2 18 ∂ 4 u ∂ t 2 ∂ y 2 i , j , k n + τ 3 162 ∂ 5 u ∂ t 3 ∂ y 2 i , j , k n + 1 3 f i , j , k n + τ 9 ∂ f ∂ t i , j , k n + τ 2 54 ∂ 2 f ∂ t 2 i , j , k n + τ 3 486 ∂ 3 f ∂ t 3 i , j , k n + O ( τ 4 ) .

Similarly, for equation (A3), we can obtain

(A8) 2 3 ∂ 2 u ∂ t 2 i , j , k n + 2 τ 9 ∂ 3 u ∂ t 3 i , j , k n + 2 τ 2 27 ∂ 4 u ∂ t 4 i , j , k n + 4 τ 3 243 ∂ 5 u ∂ t 5 i , j , k n = v i , j , k 2 ∂ 2 u ∂ z 2 i , j , k n + 2 τ 3 ∂ 3 u ∂ t ∂ z 2 i , j , k n + 2 τ 2 9 ∂ 4 u ∂ t 2 ∂ z 2 i , j , k n + 4 τ 3 81 ∂ 5 u ∂ t 3 ∂ z 2 i , j , k n + 2 3 f i , j , k n + 2 τ 9 ∂ f ∂ t i , j , k n + 2 τ 2 27 ∂ 2 f ∂ t 2 i , j , k n + 4 τ 3 243 ∂ 3 f ∂ t 3 i , j , k n + O ( τ 4 ) .

Equations (A1), (A7), and (A8) are added together, and then, by the relationships in equation (1), we obtain

(A9) ∂ 2 u ∂ t 2 i , j , k n = v i , j , k 2 ∂ 2 u ∂ x 2 i , j , k n + ∂ 2 u ∂ y 2 i , j , k n + ∂ 2 u ∂ z 2 i , j , k n + f i , j , k n + v i , j , k 2 τ 3 ∂ 3 u ∂ t ∂ x 2 − ∂ 3 u ∂ t ∂ z 2 i , j , k n + v i , j , k 2 τ 2 54 5 ∂ 4 u ∂ t 2 ∂ x 2 + 2 ∂ 4 u ∂ t 2 ∂ y 2 − 7 ∂ 4 u ∂ t 2 ∂ z 2 i , j , k n + v i , j , k 2 τ 3 162 3 ∂ 5 u ∂ t 3 ∂ x 2 + 2 ∂ 5 u ∂ t 3 ∂ y 2 − 5 ∂ 5 u ∂ t 3 ∂ z 2 i , j , k n + O ( τ 4 ) .

Through equation (A9), we can obtain the following conclusions:

Case I: when ∂ 3 u ∂ t ∂ x 2 − ∂ 3 u ∂ t ∂ z 2 ≠ 0 , the splitting error is O ( τ ) .
Case II: when ∂ 3 u ∂ t ∂ x 2 − ∂ 3 u ∂ t ∂ z 2 = 0 5 ∂ 4 u ∂ t 2 ∂ x 2 + 2 ∂ 4 u ∂ t 2 ∂ y 2 − 7 ∂ 4 u ∂ t 2 ∂ z 2 ≠ 0 , the splitting error is O ( τ 2 ) .
Case III: when ∂ 3 u ∂ t ∂ x 2 − ∂ 3 u ∂ t ∂ z 2 = 0 5 ∂ 4 u ∂ t 2 ∂ x 2 + 2 ∂ 4 u ∂ t 2 ∂ y 2 − 7 ∂ 4 u ∂ t 2 ∂ z 2 = 0 , we can obtain ∂ 3 u ∂ t ∂ x 2 = ∂ 3 u ∂ t ∂ z 2 ∂ 4 u ∂ t 2 ∂ x 2 = ∂ 4 u ∂ t 2 ∂ z 2 = ∂ 4 u ∂ t 2 ∂ y 2 ∂ 5 u ∂ t 3 ∂ x 2 = ∂ 5 u ∂ t 3 ∂ y 2 = ∂ 5 u ∂ t 3 ∂ z 2 . At this time, 3 ∂ 5 u ∂ t 3 ∂ x 2 + 2 ∂ 5 u ∂ t 3 ∂ y 2 − 5 ∂ 5 u ∂ t 3 ∂ z 2 = 0 always holds. Therefore, the splitting error is O ( τ 4 ) .

References

[1] W. Liang, C. Yang, Y. Wang, and H. Liu, Acoustic wave equation modeling with new time-space domain finite difference operators, Chin. J. Geophys. 56 (2013), no. 10, 3497–3506, https://doi.org/10.6038/cjg20131024. Search in Google Scholar

[2] Y. Liu and M. Sen, Acoustic VTI modeling with a time-space domain dispersion-relation-based finite-difference scheme, Geophys. 75 (2010), no. 3, 11–17, https://doi.org/10.1190/1.3374477. Search in Google Scholar

[3] Y. Jiang and Y. Ge, An explicit fourth-order compact difference scheme for solving the 2D wave equation, Adv. Differ. Equ. 2020 (2020), 415, https://doi.org/10.1186/s13662-020-02870-z. Search in Google Scholar

[4] Y. Liu and K. Mrinal, A new time-space domain high-order finite-difference method for the acoustic wave equation, J. Comput. Phys. 228 (2009), no. 23, 8779–8806, https://doi.org/10.1016/j.jcp.2009.08.027. Search in Google Scholar

[5] B. Finkelstein and R. Kastner, Finite difference time domain dispersion reduction schemes, J. Comput. Phys. 221 (2017), 422–438, https://doi.org/10.1016/j.jcp.2006.06.016. Search in Google Scholar

[6] H. L. Liao and Z. Z. Sun, Maximum norm error estimates of efficient difference schemes for second-order wave equations, J. Comput. Appl. Math. 235 (2011), no. 8, 2217–2233, https://doi.org/10.1016/j.cam.2010.10.019. Search in Google Scholar

[7] M. Zapata, R. Balam, and J. Urquizo, High-order implicit staggered grid finite differences methods for the acoustic wave equation, Numer Methods Partial Differ. Equ. 34 (2018), no. 2, 602–625, https://doi.org/10.1002/num.22217. Search in Google Scholar

[8] S. Britt, E. Turkel, and S. Tsynkov, A high order compact time/space finite difference scheme for the wave equation with variable speed of sound, J. Sci. Comput. 76 (2018), no. 2, 777–811, https://doi.org/10.1007/s10915-017-0639-9. Search in Google Scholar

[9] D. Yang, P. Tong, and X. Deng, A central difference method with low numerical dispersion for solving the scalar wave equation, Geophys. Prospect. 60 (2012), no. 5, 885–905, https://doi.org/10.1111/j.1365-2478.2011.01033.x. Search in Google Scholar

[10] J. Shragge and B. Tapley, Solving the tensorial 3D acoustic wave equation, Geophys. 82 (2017), no. 4, 1–58, https://doi.org/10.1190/GEO2016-0691.1. Search in Google Scholar

[11] R. Shukla and X. Zhong, Derivation of high-order compact finite difference schemes for non-uniform grid using polynomial interpolation, J. Comput. Phys. 204 (2005), no. 2, 404–429, https://doi.org/10.1016/j.jcp.2004.10.014. Search in Google Scholar

[12] S. Britt, S. Tsynkov, and E. Turkel, Numerical solution of the wave equation with variable wave speed on nonconforming domains by high-order difference potentials, J. Comput. Phys. 354 (2017), 26–42, https://doi.org/10.1016/j.jcp.2017.10.049. Search in Google Scholar

[13] K. Lipnikov, G. Manzini, and M. Shashkov, Mimetic finite difference method, J. Comput. Phys. 257 (2014), 1163–1227, https://doi.org/10.1016/j.jcp.2013.07.031. Search in Google Scholar

[14] Y. Abdulkadir, Comparison of finite difference schemes for the wave equation based on dispersion, J. Appl. Math. Phys. 3 (2015), no. 11, 1544–1562, https://doi.org/10.4236/jamp.2015.311179. Search in Google Scholar

[15] C. Chu and P. Stoffa, Implicit finite-difference simulation of seismic wave propagation, Geophys. 77 (2012), no. 2, 57–67, https://doi.org/10.1190/GEO2011-0180.1. Search in Google Scholar

[16] D. Peaceman and H. Rachford, The numerical solution of parabolic and elliptic differential equations, J. Soc. Indust. Appl. Math. 3 (1955), no. 1, 28–41, https://doi.org/10.1137/0103003. Search in Google Scholar

[17] S. Das, W. Liao and A. Gupta, An efficient fourth-order low dispersive finite difference scheme for a 2D acoustic wave equation, J. Comput. Appl. Math. 258 (2014), no. 3, 151–167, https://doi.org/10.1016/j.cam.2013.09.006. Search in Google Scholar

[18] W. Liao, On the dispersion, stability and accuracy of a compact higher-order finite difference scheme for 3D acoustic wave equation, J. Comput. Appl. Math. 270 (2014), 571–583, https://doi.org/10.1016/j.cam.2013.08.024. Search in Google Scholar

[19] W. Liao, P. Yong, H. Dastour, and J. Huang, Efficient and accurate numerical simulation of acoustic wave propagation in a 2D heterogeneous media, Appl. Math. Comput. 321 (2018), 385–400, https://doi.org/10.1016/j.amc.2017.10.052. Search in Google Scholar

[20] K. Li and W. Liao, An efficient and high accuracy finite-difference scheme for the acoustic wave equation in 3D heterogeneous media, J. Comput. Sci. 40 (2020), 101063, https://doi.org/10.1016/j.jocs.2019.101063. Search in Google Scholar

[21] H. L. Liao and Z. Z. Sun, A two-level compact ADI method for solving second-order wave equations, Int. J. Comput. Math. 90 (2013), no. 7, 1471–1488, https://doi.org/10.1080/00207160.2012.754016.Search in Google Scholar

[22] K. Li, W. Liao, and Y. Lin, A compact high-order alternating direction implicit method for three-dimensional acoustic wave equation with variable coefficient, J. Comput. Appl. Math. 361 (2019), 113–129, https://doi.org/10.1016/j.cam.2019.04.013. Search in Google Scholar

[23] A. Samarskii, Local one-dimensional difference schemes for multi-dimensional hyperbolic equations in an arbitrary region, USSR Comput. Math. Math. Phys. 4 (1964), no. 4, 21–35, https://doi.org/10.1016/0041-5553(64)90002-3. Search in Google Scholar

[24] S. Kim and H. Lim, High-order schemes for acoustic waveform simulation, Appl. Numer. Math. 57 (2007), no. 4, 402–414, https://doi.org/10.1016/j.apnum.2006.05.003. Search in Google Scholar

[25] N. Yun, C. Sun, and C. Sim, An optimal nearly analytic splitting method for solving 2D acoustic wave equations, J. Appl. Geophys. 177 (2020), 104029, https://doi.org/10.1016/j.jappgeo.2020.104029. Search in Google Scholar

[26] W. Zhang, L. Tong, and E. Chung, A new high accuracy locally one-dimensional scheme for the wave equation, J. Comput. Appl. Math. 236 (2011), no. 6, 1343–4353, https://doi.org/10.1016/j.cam.2011.08.022. Search in Google Scholar

[27] W. Zhang and J. Jiang, A new family of fourth-order locally one-dimensional schemes for the three-dimensional wave equation, J. Comput. Appl. Math. 311 (2017), 130–147, https://doi.org/10.1016/j.cam.2016.07.020. Search in Google Scholar

[28] C. Sim, C. Sun, and N. Yun, A nearly analytic symplectic partitioned Runge-Kutta method based on a locally one-dimensional technique for solving two-dimensional acoustic wave equations, Geophys. Prospec. 68 (2020), no. 4, 1253–1269, https://doi.org/10.1111/1365-2478.12924. Search in Google Scholar

[29] W. Zhang, A new family of fourth-order locally one-dimensional schemes for the 3D elastic wave equation, J. Comput. Appl. Math. 348 (2018), 246–260, https://doi.org/10.1016/j.cam.2018.08.056. Search in Google Scholar

[30] S. Lele, Compact finite difference schemes with spectral-like resolution, J. Comput. Phys. 103 (1992), no. 1, 16–42, https://doi.org/10.1016/0021-9991(92)90324-R. Search in Google Scholar

[31] D. Yang, Iterative Solution for Large Linear System, Academic Press, New York, 1991. Search in Google Scholar

[32] D. Smith, Numerical solution of partial differential equations, finite difference methods, Math. Comp. 48 (1985), no. 178, 834–835, https://doi.org/10.2307/2007849. Search in Google Scholar

[33] G. Zhang, Two conservative and linearly-implicit compact difference schemes for the nonlinear fourth-order wave equation, Appl. Math. Comput. 401 (2021), 126055, https://doi.org/10.1016/j.amc.2021.126055. Search in Google Scholar

[34] R. Mohanty and V. Gopal, A new off-step high order approximation for the solution of three-space dimensional nonlinear wave equations, Appl. Math. Model. 37 (2013), no. 5, 2802–2815, https://doi.org/10.1016/j.apm.2012.06.021. Search in Google Scholar

[35] D. Li and W. Sun, Linearly implicit and high-order energy-conserving schemes for nonlinear wave equations, J. Sci. Comput. 83 (2020), 65, https://doi.org/10.1007/s10915-020-01245-6. Search in Google Scholar

[36] M. Fardi, I. Pishkar, J. Alidousti, and Y. Khan, Numerical investigation of the MHD suction-injection model of viscous fluid using a kernel-based method, Arch. Appl. Mech. 91 (2021), 4205–4221, https://doi.org/10.1007/s00419-021-02003-2. Search in Google Scholar

[37] Y. Khan and N. Faraz, Simple use of the Maclaurin series method for linear and non-linear differential equations arising in circuit analysis, Int. J. Comput. Math. Electr. Electron. Eng. 40 (2021), no. 3, 593–601, https://doi.org/10.1108/COMPEL-08-2020-0286. Search in Google Scholar

[38] Y. Khan, A Hausdorff fractal Nizhnik-Novikov-Veselov model arising in the incompressible fluid, Int. J. Numer. Methods Heat Fluid Flow 32 (2022), no. 5, 1674–1685, https://doi/org/10.1108/HFF-03-2021-0232. Search in Google Scholar

Received: 2022-04-07

Revised: 2022-06-19

Accepted: 2022-07-15

Published Online: 2022-09-22

This work is licensed under the Creative Commons Attribution 4.0 International License.

An accurate and efficient local one-dimensional method for the 3D acoustic wave equation

Abstract

1 Introduction

2 Deduction of the HOC-LOD scheme

3 Theory analysis

3.1 Consistency analysis

Theorem 1

Proof

3.2 Stability analysis

Lemma 1

Theorem 2

Proof

3.3 Convergence analysis

Theorem 3

Proof

4 Illustrative experiments

Example 1 [27]

Example 2 [22]

Example 3 [22]

Example 4 [22]

Example 5

Example 6

5 Conclusion

Appendix

References

Journal and Issue

Articles in the same Issue