High-order compact LOD methods for solving high-dimensional advection equations

In this paper, by using the local one-dimensional (LOD) method, Taylor series expansion and correction for the third derivatives in the truncation error remainder, two high-order compact LOD schemes are established for solving the two- and three- dimensional advection equations, respectively. They have the fourth-order accuracy in both time and space. By the von Neumann analysis method, it shows that the two schemes are unconditionally stable. Besides, the consistency and convergence of them are also proved. Finally, numerical experiments are given to confirm the accuracy and efficiency of the present schemes.


Introduction
Consider the high-dimensional advection equations as follows with the initial conditions Communicated by Cassio Oishi.
B Yongbin Ge gyb@nxu.edu.cn Bo Hou houbo950211@163.com and the following periodic boundary condition For the two-dimensional (2D) advection equation, X = (x, y) are spatial variables, v = (a, b) are the advection coefficients, u = u (x, y, t) is the unknown function, D is a 2D spatial domain. And for the three-dimensional (3D) advection equation, X = (x, y, z), v = (a, b, c), u = u(x, y, z, t), D is a 3D spatial domain. Advection equation is one of the very important hyperbolic partial differential equations (PDEs), which is widely used in many fields, such as biomathematics, energy development, and aerodynamics, etc. However, because practical problems are usually very complicated, it is often difficult to obtain the exact solutions of them, therefore, it is very necessary and valuable to solve them with numerical methods.
In the past few decades, many scholars have devoted to the development of accurate and stable numerical methods for solving advection equations. Vabishchevich (2019), Vabishchevich (2018) established an unconditionally stable implicit three-level scheme with the secondorder accuracy and proposed several two-level schemes for solving the advection equations. Wei and Zhang (2010) proposed a two-parameter display difference scheme for the advection equation and discussed the stability and the calculation accuracy by using Miller theorem, the scheme is stable when the two parameters are non-negative. Erdogan (2014) used exponential fitting method to present a difference scheme. Bourchtein and Bourchtein (2012) constructed a three-level, five-point central leapfrog scheme for the advection equation, but this scheme is conditionally stable. Kim (2003) proposed an upwind leapfrog approach for solving the advection equations, which is non-dissipative and very accurate, and then was extended to higher-order and multiple dimensions. Abarbanel and Chertock (2000) developed a class of finite difference schemes for the hyperbolic initial boundary value problems in one and two space dimensions. They possess the property of strict stability. Sheu and Lee (2003) presented a class of Taylor-Galerkin (TG) finite-element models for solving the first-order hyperbolic equation which admits discontinuities, in which, five parameters were introduced for purposes of controlling stability, monotonicity and accuracy. Sarra (2008) explored the accuracy and eigenvalue stability of symmetric, asymmetric radial basis functions (RBF) collocation methods for solving some model hyperbolic initial boundary value problems in one and two dimensions. Hou and Ge (2019) constructed an unconditionally stable high-order compact difference scheme with the truncation error is O(τ 4 +τ 2 h 2 +h 4 ) to solve the one-dimensional advection equation. Fan (2015) proposed a box scheme of the variable coefficient advection equation and the energy analysis method was used to verify the unconditional stability of the scheme. Montmollin (2000) employed the space-time integrated least squares (STILS) method to solve the pure advection equation and proposed a time-marching method.
For high-dimensional differential equations, using fully implicit difference methods to get a solution, usually, a large linear system at each time step has to be resolved, which would be too expensive to be used in practice. To overcome this difficulty, operator splitting methods were proposed, such as the alternating direction implicit (ADI) (Peaceman and Rechford 1955;Douglas and Kimy 2011) and the local one-dimensional (LOD) (Samarskii 1964) methods. These methods can transfer the high-dimensional problems to some lowdimensional problems, which makes the algorithm much simpler and reduces the amount of calculation. Among splitting methods, various LOD schemes have been developed. For instance, Wang (2011) used the LOD method to obtain a scheme for the 2D non-homogeneous parabolic equations with the second-order accuracy in time direction and the fourth-order accuracy in space direction and temporal accuracy was improved to the fourth-order by using the Richardson extrapolation technique. Wang and Wang (2014) used the Richardson extrapolation algorithms and LOD method to propose two schemes for solving the 2D nonhomogeneous diffusion reaction equations. Based on the LOD strategy, Chen and Ge (2018) proposed two kinds high-order finite difference schemes for solving the 2D linear parabolic equation, whose truncation errors are O(τ 4 + h 4 ) and O(τ 6 + h 6 ), respectively. Xu (2016) constructed two improved LOD schemes for the nonlinear parabolic equations, and gave out the error and stability analysis. Karaa (2005) constructed a high-order LOD difference scheme for solving the 2D parabolic equations with Dirichlet boundary conditions. Zhao et al. (2017) developed two sets of high-order LOD difference schemes to solve the 2D and 3D heat problems with Neumann boundary conditions. On the other hand, ADI method is another commonly used operator splitting method, which has been developed extensively by many authors to solve the parabolic and hyperbolic equations. For more details, readers are referred to Tian (2011, Sun and Sun (2020), Karaa and Zhang (2004), Karaa (2006), You (2006), Liao and Sun (2013), Dai and Nassar (1998) and the references therein.
In this paper, based on the LOD method, a class of high-order compact difference schemes are introduced for solving the high-dimensional advection equations. In Sect. 2, we propose a high-order compact LOD difference scheme for solving the 2D advection equation, the consistency, stability and convergence are proved. In Sect. 3, we extend the present method to the 3D case. In Sect. 4, numerical experiments are given to demonstrate the accuracy and reliability of the present method. In Sect. 5, some conclusions are given.

2D advection equation
Now we consider the 2D advection equation as following with initial condition and periodic boundary condition where u(x, y, t) is the unknown function and is assumed to have sufficient smoothness, a and b are the advection coefficients in x and y directions, respectively, ϕ 2 (x, y) is a smooth known function, D 2 = { x, y| c ≤ x, y ≤ d}, ∂ D 2 is the boundary of D 2 .

High-order compact LOD scheme
We divide the domain D 2 × (0, T ] into an N 2 × M mesh, with h = d−c N and τ = T M , respectively. x i = c + ih, y j = c + jh, for i, j = 0, 1, . . . N are the ith and jth mesh node. In temporal dimension, the uniform step size τ is used, and t n = nτ is the nth time level. The quantity u(x i , y j , t n ) represents the exact solution, while u n i, j denotes its approximate solution.
Applying the LOD strategy, we split Eq. (4) into the following one-dimensional equations ideally: To advance the solution from t n to t n+1 , we assume that Eq. (7) holds from t n to t n+ 1 2 , and Eq. (8) holds from t n+ 1 2 to t n+1 .
Firstly, the points u n+ 1 2 i, j and u n i, j are expanded with Taylor series at the grid point (x i , y j , t n+ 1 4 ) as follows: Combining Eqs. (9) and (10), we can get ∂u ∂t Similarly, we can obtain Substituting Eqs. (11) and (12) into Eq. (7), we can obtain In which, δ t u . In order to achieve the fourth-order accuracy in both time and space directions, the third derivatives of t and x in Eq. (13) must be discretized with the second-order. And at the same time, in order to ensure the compactness of the scheme, we make the following deformations of Eq. (7): Substituting the expressions of ∂ 3 u ∂t 3 and ∂ 3 u ∂ x 3 above into Eq. (13), we have , it will degenerate accuracy of the scheme to be the second-order in time direction. In order to improve the temporal accuracy to be the fourth-order, we use Eqs. (9) and (10) to get Then, we can get Substituting Eq. (17) into Eq. (15), we have Using the definitions of δ t and δ x , and the Crank-Nicolson method for time derivative in the mixed derivative term ∂ 3 u ∂ x 2 ∂t , we get Discarding O(τ 4 + τ 2 h 2 + h 4 ) and multiplying Eq. (21) by τ , we can get the high-order compact scheme in x direction, In which, λ = τ h . Similarly, we can obtain the high-order compact scheme for Eq. (8) in y direction So the high-order compact LOD scheme for the 2D advection equation (4) is written as Remark 1 The scheme (24) is the fourth-order accuracy in both time and space and is a twostep scheme. Only the value of the initial time step requires to be known, so it does not need to construct the scheme for the computation of start-up time step, like the three-level schemes in Vabishchevich (2018Vabishchevich ( , 2019 and Bourchtein and Bourchtein (2012).

Remark 2
The scheme (24) has only three unknown points and the coefficient matrices of the scheme are tridiagonal, so we can solve the scheme (24) by using the efficient tridiagonal matrix algorithm (TDMA). However, this method requires that the main diagonal elements of the coefficient matrices can not be zero, so we let aλ, bλ = 2.

Theory analysis
Theorem 2.1 High-order compact LOD scheme (24) is the fourth-order consistent.
Proof For the first equation in the scheme (24), using the Taylor series expansions of Substituting Eqs. (25)-(28) into Eq. (24), we can obtain Making use of the relationships in Eqs. (14), (30) is equivalent to Similarly, we can get Combining Eqs. (31) and (32), we can get ∂u ∂t Thus, the scheme (24) is the fourth-order consistent to the original differential equation (4). In other words, the scheme is the fourth-order accurate in both time and space. Therefore, Theorem 1 is proved.

Theorem 2.2 High-order compact LOD scheme (24) is unconditionally stable.
Proof In the following, we use the von Neumann method to analyze the stability of the scheme (24).
Let u n i, j = η n e I σ 1 x i e I σ 2 y j ,u n+ 1 2 i, j = η n+ 1 2 e I σ 1 x i e I σ 2 y j , where I = √ −1 is imaginary unit, η is the amplitude. σ 1 and σ 2 are the wave numbers in x and y directions, respectively. For x direction, substituting u n i, j = η n e I σ 1 x i e I σ 2 y j ,u n+ 1 2 i, j = η n+ 1 2 e I σ 1 x i e I σ 2 y j into Eq. (24), and eliminating e I σ 1 x i e I σ 2 y j on both sides of the equation, we have Using the Euler formulas, e I σ 1 h = cos σ 1 h + I sin σ 1 h, e −I σ 1 h = cos σ 1 h − I sin σ 1 h, we get Similarly, we can draw Then, we can get the error amplification factor According to the von Neumann stability theorem, when |G| ≤ 1, the scheme is stable. We can draw the error amplification factor of the scheme is |G| = 1, so the scheme is unconditionally stable. (24) is unconditionally convergent.

Theorem 2.3 High-order compact LOD scheme
Proof From Theorem 2.1 and Theorem 2.2, we can get the scheme (24) is unconditionally convergent by the Lax equivalence theorem (Smith 1985).

Extension to the 3D case
Next, we consider the 3D advection equation as following: with initial condition and periodic boundary condition u (d, y, z, t) = u(e, y, z, t), u(x, d, z, t) = u(x, e, z, t), u(x, y, d, t Where u(x, y, z, t) is the unknown function, which is assumed to have sufficient smoothness. a, b and c are the advection coefficients in x, y and z directions, respectively. ϕ 3 (x, y, z) is a smooth known function, Applying the LOD strategy, we split Eq. (39) into the following one-dimensional equations ideally: To advance the solution from t n to t n+1 , we assume that Eq. (42) holds from t n to t n+ 1 3 , Eq. (43) holds from t n+ 1 3 to t n+ 2 3 and Eq. (44) holds from t n+ 2 3 to t n+1 . Using the method similar to the 2D advection equation, a high-order compact LOD scheme for the 3D advection equation can be obtained as follows: u n i+1, j,k + 2 3 − a 2 λ 2 6 u n i, j,k + 1 6 + aλ 4 + a 2 λ 2 12 u n i−1, j,k  (45) is the fourth-order accuracy in both time and space and is a twostep scheme. It involves three unknowns for each time step. We can solve three tridiagonal linear equations by the TDMA algorithm for each time step according to Eq. (45). Because it requires that the main diagonal elements of the coefficient matrices can not be zero, we set aλ, bλ, cλ = 2.
Remark 4 In Sect. 2, we have proved the consistency, stability and convergence of the 2D high-order compact LOD scheme (24). Similarly, we can proved the scheme (45) is also the fourth-order consistent, unconditionally stable and convergent by using the same method.

Numerical experiments
In order to verify the accuracy and reliability of the present method in this article, we select the following five examples to conduct numerical experiments and compare them with the calculated results in the existing literature. All calculations in the article are programmed in Fortran 77 language. We conduct our computations with double precision arithmetic and perform them on a Intel IV/Dual-Core/2.6 GHz private computer with 8 GB memory. The convergence rate can be got by the following formula:

Problem 1
We consider the 2D advection equation with the initial and periodic boundary conditions The exact solution is Table 1 gives the maximum absolute errors and the convergence rate of the Beam-Warming scheme (Beam and Warming 1976), the RHOC-ADI scheme (Tian 2011) and the LOD scheme (24) with τ = 0.5h, t = 0.2, a = b = 1 for Problem 1. Because the RHOC-ADI scheme is the second-order in time and the fourth-order in space, we also compute it with τ = h 2 . The data in Table 1 show that the present LOD scheme (24) is the fourth-order in both space and time, which is consistent with our theoretical analysis. While the RHOC-ADI scheme achieves spatial fourth-order accuracy only if τ = O(h 2 ). Otherwise, spatial accuracy is just Table 1 Maximum absolute error and convergence rate when t = 0.2, a = 1, b = 1 for Problem 1 N Beam-Warming scheme (Beam and Warming 1976) RHOC-ADI scheme (Tian 2011) RHOC-ADI scheme (Tian 2011) LOD scheme (24)  (24) is better (smaller errors) than that of the Beam-Warming scheme and the RHOC-ADI scheme. In Fig. 1(a) gives the exact solution, b-d give the absolute error of the LOD scheme (24), the RHOC-ADI scheme and the Beam-Warming scheme when N = 20, t = 1, τ = 0.05 for Problem 1. We can see that the absolute error of the LOD scheme (24) is smaller than that of the RHOC-ADI scheme and the Beam-Warming scheme. Figure 1(e) gives comparison of the numerical solutions of the three schemes with the exact solution when x = 0.5, N = 20, t = 1, τ = 0.05 for Problem 1. It shows that the LOD scheme (24) is the closest to the exact solution.

Problem 2
We consider the 2D advection equation with the initial and periodic boundary conditions The exact solution is Table 2 gives the maximum absolute errors, the convergence rate and the CPU time of the Beam-Warming scheme, the RHOC-ADI scheme and the LOD scheme (24) when τ = 0.5, a = 1, b = 1 2 for Problem 2. It is easy to find the computed results of the LOD scheme (24) is superior to the Beam-Warming scheme when τ = 0.5h and the RHOC-ADI scheme with τ = h 2 . The CPU time of the LOD scheme (24) is less than that of the RHOC-ADI scheme. Because the Beam-Warming scheme is second-order in time and space, its CPU time is less than that of the LOD scheme (24). Figure 2 gives the exact solution (a), the absolute error (b,c,d) (24) is the minimal under the same condition. Besides, we can also see more directly that the numerical solution of the LOD scheme (24) is the closest to the exact solution compared to other two schemes from Fig. 2(e).

Problem 3
We consider the 2D advection equation The exact solution is The problem 3 is a Gauss wave case. In fact, it is an initial value problem, so we set the boundary value to be zero in a enough big spatial domain. Table 3 gives the maximum absolute error of the Beam-Warming scheme, the RHOC-ADI scheme and the LOD scheme(24) when t = 5, a = 0.1, b = 0.1, p = 0, q = 20, σ = 1, τ = 0.025, (x, y) = (4, 4) for Problem 3. It is easy to find the computed results of the RHOC-ADI scheme and the LOD scheme (24) is better than that of the Beam-Warming scheme. Figure 3 gives the numerical solution of the LOD scheme(24) on t = 10 (a, b), t = 50 (c, d), t = 80 (e, f), t = 120 (g, h) when N = 100, τ = 0.0025, a = 0.1, b = 0.1, p = 0, q = 20, σ = 1, (x, y) = (4, 4) and compared with the exact solution for Problem 3. Figure 4 gives the numerical solution of the LOD scheme (24) on t = 100 (a, b), t = 200 (c, d), t = 400 (e, f) when N = 200, τ = 0.1, a = 0.1, b = 0.1, c = 0, d = 100, σ = 10, (x, y) = (30, 30) and compared with the exact solution for Problem 3. It is not difficult to see that the movement track of the wave. It also can be seen from these figures that the simulation result is well in terms of phase and amplitude. In other words, the numerical solution of the LOD scheme (24) is consistent with the exact solution.

Problem 4
We consider the 3D advection equation with the initial and periodic boundary conditions u(x, y, z, 0) = sin π x a + y b + z c u (0, y, z, t) = u(2, y, z, t), The exact solution is Table 4 gives the maximum absolute errors and the convergence rate of the Beam-Warming scheme, the RHOC scheme and the LOD scheme(45) with τ = 0.5h and the RHOC scheme with τ = h 2 when t = 0.2, a = 1 2 , b = 1 3 , c = 1 4 for Problem 4. It is known that the LOD scheme(45) is the fourth-order accurate in time and space, so it has higher-order accuracy in time direction compared to the RHOC scheme and the Beam-Warming scheme. In addition, the computed results of the LOD scheme (45) are also better than that of the RHOC scheme and the Beam-Warming scheme. Figure 5 gives the exact solution (a), the absolute error and the contours of Problem 4 from the Beam-Warming scheme (d, g), the RHOC scheme (c, f) and the LOD scheme (45) (b, e) with N = 32, t = 1, τ = 0.1. It is easy to obverse that the error of the scheme (45) is less than that of the Beam-Warming scheme and the RHOC scheme from Fig. 5(b-d). From Fig. 5(e-g), we can find that the numerical solution of the LOD scheme (45) is consistent with the exact solution, while the RHOC scheme and the Beam-Warming scheme produce relatively big errors.

Problem 5
We consider the 3D advection equation with the initial and periodic boundary conditions u (0, y, z, t) = u(2, y, z, t), t > 0 u(x, 0, z, t) = u(x, 2, z, t), t > 0 u(x, y, 0, t) = u(x, y, 2, t), t > 0 The exact solution is  Table 5 gives the maximum absolute errors, convergence rate and CPU time by using the Beam-Warming scheme, the RHOC scheme and the LOD scheme (45) with t = 0.5, a = 1, b = 1, c = 1. It is easy to see that the accuracy and the CPU time of the scheme (45) are better than that of the RHOC scheme. Figure 6 gives the exact solution (a), numerical solution, absolute error and contours of the LOD scheme (b, e, h), the RHOC scheme (c, f, i) and the Beam-Warming scheme (d, g, j) with N = 20, t = 1, τ = 0.05. It can find that the numerical solution of the LOD scheme (45) is in good agreement with the exact solution  through obverse Fig. 6(b-d, h-j). For Fig. 6(e-g), the error of the LOD scheme (45) can achieve 10 −2 , but for the Beam-Warming scheme, and the RHOC scheme, the errors just achieve 10 0 and 10 −1 . Therefore, the LOD scheme (45) is better in the accuracy compared with the RHOC scheme and the Beam-Warming scheme. Besides, the CPU time of the LOD scheme (45) is less than that of the RHOC scheme.

Concluding remarks
In this study, we have proposed two high-order compact LOD difference schemes for solving the 2D and 3D advection equations, respectively. The scheme are two-step, unconditionally stable and have the fourth-order accuracy in both time and space. Firstly, based on the LOD strategy, we separate the 2D equations into two 1D equations. Then, for the 1D equations, the Taylor series expansion and correction for the third derivative in the truncation error remainder of the central difference scheme are used for discretization of time and space. Besides, the consistency, stability and convergence are proved. Afterwards, we extend the present method to the 3D case. At last, the accuracy, convergence and efficiency of the present schemes are validated by some numerical experiments. It shows that the present schemes are more accurate and efficient than those schemes from Tian (2011)  During the establishment of the present high-order compact schemes, for the constant coefficient problems, time derivative is approximated with the Crank-Nicolson method combined with its correction of the truncation error remainder. Objective differential equations, e.g., Eqs. (7), (8), (42), (43), are used to treat the higher derivatives in the truncation error terms. However, if we extend the present high-order scheme to the variable coefficient problems, e.g., 1 2 ∂u ∂t + a(x, t) ∂u ∂ x = 0, we need to approximate the ( ∂ 2 a ∂t 2 ) n+ 1 2 i, j to the fourth-order. This cannot be achieved with the two-level and three-point stencil. In another words, the scheme would be the second-order instead of the fourth-order in time. Therefore, the method cannot be extended to the variable coefficient problems. On the other hand, the present method requires the unknown function to be sufficiently smooth, at least the third derivatives are continuous. But for the shock wave problems, we know, it does not meet this condition, So the present method also can not be used to catch shocks waves. For this problem, the conservation law should be considered and some special numerical methods should be selected and used. Concerning these two problems, we will continue to study in our follow-up work through using other effective methods.