Emigration, Migration und Kultur

Chen and Ge Advances in Difference Equations ( 2018) 2018:361 https://doi.org/10.1186/s13662-018-1825-2 RESEARCH Open Access High order locally one-dimensional methods for solving two-dimensional parabolic equations Jianhua Chen1 and Yongbin Ge1* * Correspondence: gyb@nxu.edu.cn Institute of Applied Mathematics and Mechanics, Ningxia University, Yinchuan, China 1 Abstract Based on the locally one-dimensional strategy, we propose two high order finite difference schemes for solving two-dimensional linear parabolic equations. In the first method, fourth order approximation in space and (2, 2) Padé formula in time are considered. These lead to a fourth order finite difference scheme in both space and time. For the second method, we employ sixth order approximation in space and (3, 3) Padé formula in time. This yields a novel sixth order scheme in both space and time. The methods are proved to be unconditionally stable, and the Sheng–Suzuki barrier is successfully avoided. Numerical experiments are given to illustrate our conclusions as well as computational effectiveness. Keywords: Parabolic equations; Locally one-dimensional strategy; Padé approximations; High order methods; Unconditional stability 1 Introduction Consider the following two-dimensional linear parabolic equation:   2 ∂u ∂ u ∂ 2u , =a + ∂t ∂x2 ∂y2 (x, y, t) ∈  × [0, T], (1) together with the initial condition u(x, y, 0) = φ0 (x, y) and boundary conditions u(0, y, t) = g0 , u(1, y, t) = g1 , u(x, 0, t) = d0 , u(x, 1, t) = d1 , where  = [0, l] × [0, l] is a spatial domain and l is a positive real number. φ0 is a sufficiently smooth function and g0 , g1 , d0 , d1 are constants. Many efforts have been made to the development of accurate and stable methods for the numerical solution of (1) [1–19]. Various high order methods [1–4, 6–9, 11–19] have been © The Author(s) 2018. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Chen and Ge Advances in Difference Equations ( 2018) 2018:361 Page 2 of 17 proposed. Among them, splitting strategies including alternating direction implicit (ADI) and locally one-dimensional (LOD) methods have been extensively explored for high order difference schemes [4, 6–9, 11–17]. These methods are extremely efficient for solving multi-dimensional equations by converting multi-dimensional equations to successions of one-dimensional equations. Subsequently, only sequences of linear tri-diagonal systems need to be solved. Recently, Dai and Nassar [15], Karra [16], developed a high order ADI difference scheme and a high order LOD difference scheme, respectively, for solving the two-dimensional parabolic equations with Dirichlet boundary conditions. Zhao et al. [17] proposed two sets of high order LOD difference schemes to solve the two- and three-dimensional heat problems with Neumann boundary conditions. Qin [7] proposed a family of ADI methods for the three-dimensional nonhomogeneous parabolic equation. All these schemes obtain fourth order accuracy in space, but only second order accuracy in time, since the Crank–Nicolson implicit method is employed for time discretization. Very recently, a semi-discrete method and Padé approximations or the Runge–Kutta methods were exploited to increase the temporal accuracy [19–24]. Vu and Alexander [19] developed a series of explicit exponential Runge–Kutta methods of high order for parabolic problems. For the one-dimensional hyperbolic equation, Gao and Chi [20] used semi-discrete method to transform it into a system consisting of ordinary differential equations with respect to time, whose exact solution containing an infinite matrix series was approximated by (1, 1) and (2, 2) Padé approximations. They obtained two schemes with third and fifth order accuracy in time, respectively. Liu et al. [21] used a similar strategy, (2, 2) and (3, 3) Padé approximations for the time discretization and C 3 quartic spline approximation for space discretization, to get two higher order difference schemes for the one-dimensional linear hyperbolic equation. Zhang [23] provided a (3, 3) Padé approximation method for solving space fractional Fokker–Planck equations. Liu et al. [25] developed a series of compact implicit schemes of fourth and sixth orders for solving differential equations involved in geodynamics simulations. And Liu et al. [26] proposed a sixth order accuracy solution to a system of nonlinear differential equations with coupled compact method. This paper targets at the development of two high order LOD finite difference schemes for solving two-dimensional parabolic equations. Sheng–Suzuki accuracy barrier [27] is avoided. To achieve high order accuracies in both time and space, we successfully combine high order approximations in the spatial discretization with techniques of semi-discrete and high order Padé approximations in the temporal discretization [20, 21, 23]. The outline of this paper is as follows: Sect. 2 presents two splitting methods for solving (1). Their stabilities are analyzed in Sect. 3. Numerical verification is carried out in Sect. 4. Further comments and conclusions are given in Sect. 5. 2 High order difference methods We superimpose the space-time domain  × [0, T] by an N × N × M mesh, where N , M are positive integers. Let h = l/N denote the step size of space and τ = T/M for step size of time. We further define xi = ih, yj = jh, i, j = 0, 1, 2, . . . , N, tk = kτ , k = 0, 1, 2, . . . , M. Chen and Ge Advances in Difference Equations ( 2018) 2018:361 Page 3 of 17 Applying an LOD strategy, we split Eq. (1) to the following one-dimensional equations ideally: ∂ 2u 1 ∂u =a 2, 2 ∂t ∂x (2) ∂ 2u 1 ∂u =a 2. 2 ∂t ∂y (3) 1 To advance the solution from t k to t k+1 , we assume that Eq. (2) holds from t k to t k+ 2 , and Eq. (3) holds from t k+ 21 to t k+1 , respectively. We will investigate several high order semi- discretization methodologies for dealing with (2) and (3), respectively. 2.1 O(τ 4 + h4 ) finite difference method Firstly, we build the fourth order scheme in space for Eq. (2). To this end, we employ a fourth order compact formula in [28] to discrete the second spatial derivative in x:  ∂ 2u ∂x2 –1    h2 δx2 uij + O h4 , = 1 + δx2 12 ij  (4) where δx2 is the second order central difference operator in the x-direction, which is defined as δx2 uij = (ui+1j –2uij +ui–1j )/h2 . Substituting Eq. (4) into Eq. (2) and dropping the truncated error O(h4 ), we acquire the following relation:       5 ∂u 1 ∂u 1 ∂u + + 24 ∂t i+1,j 12 ∂t i,j 24 ∂t i–1,j = 2a a a ui+1,j – 2 ui,j + 2 ui–1,j . h2 h h (5) This semi-discrete scheme is of fourth order accuracy in space. Consider the initial and boundary conditions. The above can be rewritten into the following: ⎧ ⎨A dU(t) = BU(t) + G, dt ⎩U(0) = φ0 , (6) which is a system of ordinary differential equations with an unknown function at each spatial grid point, where U(t) = [u1j , u2j , . . . , uN–1j ]T , G = 1, and the matrices A and B of order N – 1 are given by ⎡ 5 ⎢ 12 1 ⎢ 24 ⎢ ⎢ A=⎢ ⎢ ⎢ ⎣ 1 24 5 12 1 24 .. .. . . 1 24 ⎤ .. . 5 12 1 24 ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ 1 ⎥ 24 ⎦ 5 12 (N–1)×(N–1) a [g , 0, . . . , 0, g1 ]T , j h2 0 = 1, 2, . . . , N – Chen and Ge Advances in Difference Equations ( 2018) 2018:361 Page 4 of 17 and ⎡ –2 ⎢ ⎢1 a2 ⎢ ⎢ B= 2 ⎢ h ⎢ ⎢ ⎣ 1 –2 ⎤ 1 .. . .. . 1 .. . –2 1 ⎥ ⎥ ⎥ ⎥ , ⎥ ⎥ ⎥ 1⎦ –2 (N–1)×(N–1) Eq. (6) can be comprised to ⎧ ⎨ dU(t) = A–1 BU(t) + A–1 G, dt ⎩U(0) = φ0 . (7) The exact solution to (7) is –1 B U(t) = –B–1 G + etA   U(0) + B–1 G . (8) Discretizing Eq. (8) in temporal variable t, we get  1  1 –1  U t k+ 2 = –B–1 G + e(k+ 2 )τ A B U(0) + B–1 G . (9) After rearranging the above, we arrive at  1  τ –1  –1 –1 U t k+ 2 = –B–1 G + e 2 A B ekτ A B U(0) + ekτ A B B–1 G . (10) Namely   1   τ –1 τ –1 U t k+ 2 = e 2 A B – I B–1 G + e 2 A B U(tk ), (11) τ –1 where I is an identify matrix of order N – 1. Now, approximate e 2 A B to get the numerical solution of fourth order in time is the key issue. The (2, 2) Padé approximation [29] is an efficient approximation to eZ of fourth order, i.e., –1      12 + 6Z + Z2 + O Z4 . eZ = 12 – 6Z + Z2 τ –1 B Replacing the e 2 A e τ A–1 B 2 (12) by the (2, 2) Padé approximation, we get         τ –1 2 τ –1 2 –1 –1 –1 12 + 3τ A B + . A B A B = 12 – 3τ A B + 2 2 (13) Substituting Eq. (13) into Eq. (11), we acquire the following difference scheme:  1 U t k+ 2 =          τ –1 2 –1 τ –1 2 –1 –1 12 – 3τ A B + A B 12 + 3τ A B + A B – I B–1 G 2 2         τ –1 2 τ –1 2 –1 –1 –1 A B 12 + 3τ A B + A B + 12 – 3τ A B + 2 2   (14) × U tk . Chen and Ge Advances in Difference Equations ( 2018) 2018:361 Page 5 of 17 1 This recurrence relation is used to calculate U from t k to t k+ 2 . We can easily find that Eq. (14) gets fourth order accuracy in both time and space. A similar approach is used to 1 tackle Eq. (3) to obtain a recurrence relation which can calculate from t k+ 2 to t k+1 . Fourth order approximation for the second spatial derivative of y is given by  ∂ 2u ∂y2 –1    h2 δy2 uij + O h4 , = 1 + δy2 12 ij  (15) where δy2 is the second order central difference operator in the y-direction which is defined as δy2 uij = (uij+1 – 2uij + uij–1 )/h2 . Substituting Eq. (15) into Eq. (3) and dropping the truncated errors O(h4 ), we obtain the following semi-discrete scheme of fourth order accuracy in space:       1 ∂u 5 ∂u 1 ∂u + + 24 ∂t i,j+1 12 ∂t i,j 24 ∂t i,j–1 = a 2a a ui,j+1 – 2 ui,j + 2 ui,j–1 . h2 h h (16) Consider the initial and boundary conditions, they can be rewritten as a system of ordinary differential equations: ⎧ ⎨ dU(t) = A–1 BU(t) + A–1 F, dt ⎩U(k + 1 ) = φ 1 , (17) k+ 2 2 where U(t) = [ui1 , ui2 , . . . , uiN–1 ]T , F = ha2 [d0 , 0, . . . , 0, d1 ]T , i = 1, 2, . . . , N – 1, and the matrices A and B are defined as before. The exact solution of Eq. (17) is tA–1 B –1 U(t) = –B F + e     1 –1 +B F . U k+ 2 (18) We discretize Eq. (18) in time and receive      k+1  1 –1 (k+1)τ A–1 B –1 U k+ = –B F + e +B F . U t 2 (19) Rearranging Eq. (19) leads to   τ –1   1  τ –1 U t k+1 = e 2 A B – I B–1 F + e 2 A B U t k+ 2 . τ –1 B Replacing the e 2 A   U t k+1 = (20) by the (2, 2) Padé approximation, we obtain        τ –1 2 τ –1 2 –1 12 + 3τ A–1 B + – I B–1 F A B A B 2 2         τ –1 2 –1 τ –1 2 12 + 3τ A–1 B + + 12 – 3τ A–1 B + A B A B 2 2  k+ 1  (21) ×U t 2 .  12 – 3τ A–1 B +  Chen and Ge Advances in Difference Equations ( 2018) 2018:361 Page 6 of 17 Combining Eq. (14) with Eq. (21), we obtain a fourth-order difference scheme in both time and space for solving Eq. (1) as follows: ⎧ 1 ⎪ U(t k+ 2 ) = {[12 – 3τ A–1 B + ( τ2 A–1 B)2 ]–1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ × [12 + 3τ A–1 B + ( τ2 A–1 B)2 ] – I}B–1 G ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ + {[12 – 3τ A–1 B + ( τ2 A–1 B)2 ]–1 ⎪ ⎪ ⎪ ⎪ ⎨ × [12 + 3τ A–1 B + ( τ A–1 B)2 ]}U(t k ), 2 ⎪ ⎪ U(t k+1 ) = {[12 – 3τ A–1 B + ( τ2 A–1 B)2 ]–1 ⎪ ⎪ ⎪ ⎪ ⎪ × [12 + 3τ A–1 B + ( τ2 A–1 B)2 ] – I}B–1 F ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ + {[12 – 3τ A–1 B + ( τ2 A–1 B)2 ]–1 ⎪ ⎪ ⎪ ⎪ 1 ⎩ × [12 + 3τ A–1 B + ( τ A–1 B)2 ]}U(t k+ 2 ). (22) 2 Applying Eq. (22), we may complete the entire calculation from t k to t k+1 . Due to the use of fourth order scheme for discretizing the space variables and (2, 2) Padé approximation for the temporal variable, it is not difficult to find that equations in (22) are of fourth order accuracy in both time and space. 2.2 O(τ 6 + h6 ) finite difference method In [2], the authors presented a sixth order approximation for the second order derivative together with the constant boundary conditions. To this end, we have  u0 – 2u1 + u2 a h2     121 127 1 1 1 d 3 u0 + u1 + u2 – u3 + u4 + O h6 , = 2 dt 32 150 1200 150 2400   ui–1 – 2ui + ui+1 a h2     1 d 1 1 97 1 1 = – ui–2 + ui–1 + ui + ui+1 – ui+2 + O h6 2 dt 240 10 120 10 240  for i = 2, 3, 4, . . . , N – 2 (23) (24) and a  uN–2 – 2uN–1 + uN h2     121 127 1 1 1 d 3 uN – uN–1 + uN–2 – uN–3 + uN–4 + O h6 . = 2 dt 32 150 1200 150 2400  (25) Based on the above, we discrete the second order spatial derivative of x in Eq. (2). This gives the following system of ordinary differential equations: ⎧ ⎨C dU(t) = BU(t) + P, dt ⎩U(0) = φ0 . (26) Chen and Ge Advances in Difference Equations ( 2018) 2018:361 Page 7 of 17 Matrix B is defined as before. And matrix C is defined as follows: ⎡ 121 ⎢ 300 ⎢ 1 ⎢ 20 ⎢ 1 ⎢– 480 ⎢ ⎢ C=⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 127 2400 97 240 1 20 1 – 300 1 20 97 240 1 4800 1 – 480 1 20 .. .. .. . . 1 – 480 . 1 20 1 – 480 1 4800 ⎤ 1 – 480 .. . .. 97 240 1 20 1 – 300 1 20 97 240 127 2400 . ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ 1 ⎥ – 480 ⎥ ⎥ 1 ⎥ 20 ⎦ 121 300 , (N–1)×(N–1) where U(t) = [u1j , u2j , . . . , uN–1j ]T , P = ha2 [g0 , 0, . . . , 0, g1 ]T , and j = 1, 2, . . . , N – 1. The exact solution to the system of ordinary differential equations can be formed as follows: U(t) = –B–1 g + etC –1 B   U(0) + B–1 P . (27) We discretize Eq. (27) in time and rearrange it, we obtain   1   τ –1 τ –1 U t k+ 2 = e 2 C B – I B–1 P + e 2 C B U(tk ). (28) The (3, 3) Padé approximation [29] is employed to get sixth order accuracy for temporal variable –1      120 + 60Z + 12Z2 + Z3 + O Z6 . eZ = 120 – 60Z + 12Z2 – Z3 τ So e 2 C –1 B τ e2C (29) can be written as –1 B      2 τ –1 3 –1 = 120 – 30τ C –1 B + 3 τ C –1 B – C B 2      2 τ –1 3 C B . × 120 + 30τ C –1 B + 3 τ C –1 B + 2 (30) Substituting Eq. (30) into Eq. (28) and rearranging it, we gain  1 U t k+ 2 =      –1 2 τ –1 3 –1 –1 120 – 30τ C B + 4 τ C B – C B 2       –1 2 τ –1 3 –1 C B – I B–1 P × 120 + 30τ C B + 3 τ C B + 2      –1 2 τ –1 3 –1 –1 + 120 – 30τ C B + 4 τ C B – C B 2        –1 2 τ –1 3 –1 × 120 + 30τ C B + 3 τ C B + C B U tk . 2 (31) Similarly, combining the O(h6 ) approximation method for discretizing the space variable with y the (3, 3) Padé approximation for discretizing the temporal variable t to solve Eq. (3), Chen and Ge Advances in Difference Equations ( 2018) 2018:361 Page 8 of 17 1 we can get a recurrence relation which can accomplish the calculation from t k+ 2 to t k+1 , i.e.,   U t k+1 =     –1 2 τ –1 3 –1 C B 120 – 30τ C B + 3 τ C B – 2       –1 2 τ –1 3 –1 – I B–1 Q × 120 + 30τ C B + 3 τ C B + C B 2      –1 2 τ –1 3 –1 –1 C B + 120 – 30τ C B + 3 τ C B – 2      –1 2  1 τ –1 3 –1 C B U t k+ 2 . × 120 + 30τ C B + 3 τ C B + 2  –1 (32) The matrices B and C are defined as before, and the vector Q can be easily got. To achieve recurrence calculation from t k to t k+1 , combining Eq. (31) with Eq. (32), we obtain ⎧ 1 ⎪ ⎪U(t k+ 2 ) = {[120 – 30τ C –1 B + 3(τ C –1 B)2 – ( τ2 C –1 B)3 ]–1 ⎪ ⎪ ⎪ ⎪ ⎪ × [120 + 30τ C –1 B + 3(τ C –1 B)2 + ( τ2 C –1 B)3 ] – I}B–1 P ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ + {[120 – 30τ C –1 B + 3(τ C –1 B)2 – ( τ2 C –1 B)3 ]–1 ⎪ ⎪ ⎪ ⎪ ⎨ × [120 + 30τ C –1 B + 3(τ C –1 B)2 + ( τ C –1 B)3 ]}U(t k ), 2 ⎪ ⎪ U(t k+1 ) = {[120 – 30τ C –1 B + 3(τ C –1 B)2 – ( τ2 C –1 B)3 ]–1 ⎪ ⎪ ⎪ ⎪ ⎪ × [120 + 30τ C –1 B + 3(τ C –1 B)2 + ( τ2 C –1 B)3 ] – I}B–1 Q ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ + {[120 – 30τ C –1 B + 3(τ C –1 B)2 – ( τ2 C –1 B)3 ]–1 ⎪ ⎪ ⎪ ⎪ 1 ⎩ × [120 + 30τ C –1 B + 3(τ C –1 B)2 + ( τ C –1 B)3 ]}U(t k+ 2 ). (33) 2 Because of the use of the O(h6 ) approximation method in space and the (3, 3) Padé approximation in time, it is easy to see local truncation errors of the two equations in (33) to be O(τ 6 + h6 ). 3 Stability and convergence analysis Proposition 1 Assume that λ is an eigenvalue of matrix A–1 B, and x, a vector of dimension N – 1, is a corresponding eigenvector. Then λ is real, and furthermore, λ ≤ 0. Proof Let λ and x be eigenvalues and corresponding eigenvector of matrix A–1 B, respectively. They satisfy the following condition: A–1 Bx = λx or xT Bx = λxT Ax. Since  a2  2 x + x22 + · · · + x2N–1 – x1 x2 – · · · – xN–2 xN–1 h2 1  a2  = – 2 x21 + x2N–1 + (x1 – x2 )2 + (x2 – x3 )2 + · · · (xN–2 – xN–1 )2 , h xT Bx = –2 Chen and Ge Advances in Difference Equations ( 2018) 2018:361 Page 9 of 17 hence xT Bx < 0 and  1 5 2 x + x22 + · · · + x2N–1 + (x1 x2 + x2 x3 + · · · + xN–2 xN–1 ) 12 1 12     1 1 2 = x21 + x2N–1 + x2 + x23 + · · · + x2N–2 8 12  1 (x1 + x2 )2 + (x2 + x3 )2 · · · + (xN–2 + xN–1 )2 + 24  1 + x21 + x22 + x23 + · · · + x2N–2 + x2N–1 . 4 xT Ax = Hence xT Ax > 0. The above two results indicate that λ is real and λ ≤ 0.  Proposition 2 Assume that μ is an eigenvalue of matrix C –1 B, and x, a vector of dimension N – 1, is a corresponding eigenvector. Then μ is real and satisfies μ ≤ 0. Proof For matrix C, we get xT Cx = 1 1 1 97 2 121 2 127 x1 + x1 x2 – x1 x3 + x1 x4 + x1 x2 + x 300 2400 300 4800 20 240 1 1 1 1 97 2 1 + x2 x3 – x2 x4 – x1 x3 + x2 x3 + x + x3 x4 20 480 480 20 240 3 20 1 1 97 2 1 1 1 x3 x5 – x2 x4 + x3 x4 + x + x4 x5 – x4 x6 – 480 480 20 240 4 20 480 1 97 2 1 1 xN–5 xN–3 + xN–4 xN–3 + x + xN–3 xN–2 + ··· – 480 20 240 N–3 20 1 1 1 97 2 – xN–3 xN–1 – xN–4 xN–2 + xN–3 xN–2 + x 480 480 20 240 N–2 1 1 127 1 xN–4 xN–1 – xN–3 xN–1 + xN–2 xN–1 + xN–2 xN–1 + 20 4800 300 2400 121 2 x . + 300 N–1 Applying the inequalities 2xy ≥ –x2 – y2 and 2xy ≤ x2 + y2 , Chen and Ge Advances in Difference Equations ( 2018) 2018:361 Page 10 of 17 we obtain      1 1 1 1 1 121 1 127 x Cx ≥ + – x21 – + + – 300 2 2400 4800 20 2 300 480     97 1 127 1 3 2 + – × x2 – + 240 2 2400 20 2 480 2    97 1 3 4 1 1 + – x23 – × + – 240 2 20 2 480 300      N–5   97 1 4 1 1 4 4 4 97 – × x24 + x2i – + – – + 240 2 20 4800 2 480 240 20 480 i=5     1 1 4 97 1 4 – × x2 – + + 240 2 20 4800 2 480 N–4    4 1 1 3 97 1 + – × + – – x2N–3 240 2 20 2 480 300     1 3 2 97 1 127 + – × x2 – + 240 2 2400 20 2 480 N–2    1 1 121 1 127 – + + + 300 2 2400 4800 20   1 1 1 + – x2N–1 > 0, – 2 300 480 T and we also have proved that xT Bx < 0. According to the two results, we obtain that the eigenvalue μ of matrix C –1 B is real and satisfies μ ≤ 0.  Theorem 1 Finite difference schemes (22) and (33) are unconditionally stable, respectively. Proof Let λi (i = 1, 2, . . . , N – 1) be eigenvalues of matrix A–1 B, then 12 + 3τ λi + ( τ2 λi )2 ≤ 1, 12 – 3τ λi + ( τ2 λi )2 and thus   12 + 3τ λi + ( τ2 λi )2 ρ(H1 ) = max ≤ 1. i 12 – 3τ λi + ( τ2 λi )2 Let μi (i = 1, 2, . . . , N – 1) be eigenvalues of matrix C –1 B, then 120 + 30τ μi + 3(τ μi )2 + ( τ2 μi )3 ≤ 1, 120 – 30τ μi + 3(τ μi )2 – ( τ2 μi )3 we get   120 + 30τ μi + 3(τ μi )2 + ( τ2 μi )3 ρ(H2 ) = max ≤ 1, i 120 – 30τ μi + 3(τ μi )2 – ( τ2 μi )3 Chen and Ge Advances in Difference Equations ( 2018) 2018:361 Page 11 of 17 where         τ –1 2 –1 τ –1 2 A B A B 12 + 3τ A–1 B + , H1 = 12 – 3τ A–1 B + 2 2      2 τ –1 3 –1 C B H2 = 120 – 30τ C –1 B + 3 τ C –1 B – 2      2 τ –1 3 , C B × 120 + 30τ C –1 B + 3 τ C –1 B + 2 ρ(H1 ) and ρ(H2 ) are the spectral radii of H1 and H2 , respectively. This shows that the new developed schemes (22) and (33) are unconditionally stable bypassing the accuracy barrier  theorem [27]. Theorem 2 Difference schemes (22) and (33) are unconditionally convergent, respectively. Proof From the derivation process of the schemes, it is readily seen that local truncation errors of (22) and (33) are O(τ 4 + h4 ) and O(τ 6 + h6 ), respectively. So, they are consistent with the initial-boundary problem considered. And in Theorem 1, we have proved that the two schemes are unconditionally stable, therefore we can naturally declare that the difference schemes (22) and (33) are unconditionally convergent by the Lax equivalence  theorem [30] regardless of the Courant number. 4 Numerical experiments In this section, we consider two problems whose exact solutions are given to numerically illustrate the validity and effectiveness of our schemes as compared with those of the Peaceman–Rachford (P–R) method [11]. All computations are completed on a P4/2.4G private computer using double precision arithmetic. We estimate the rate of convergence of each method through the asymptotic formula Rate = log(Error(h1 )/ Error(h2 )) , log(h1 /h2 ) in which Error(h1 ) and Error(h2 ) are L2 -norm errors based on different mesh sizes h = h1 and h = h2 , respectively. Problem 1 ∂u ∂ 2 u ∂ 2 u = + , ∂t ∂x2 ∂y2 0 ≤ x, y ≤ 1, t > 0, together with proper initial-boundary conditions. The exact solution of this problem is 2 u(x, y, t) = e–2π t sin(πx) sin(πy). Table 1 gives the L2 -norm errors (Error) and the convergence rate (Rate) by using scheme (22), scheme (33), and the P–R scheme with h = τ = 1/10, 1/20, 1/30, 1/40 at T = 1. The data in Table 1 show that scheme (22) and scheme (33) have overall fourth order and sixth order accuracy, respectively. Chen and Ge Advances in Difference Equations ( 2018) 2018:361 Page 12 of 17 Table 1 L2 -norm error and convergence rate for Problem 1 with h = τ at T = 1 P–R scheme h Error 1/10 1/20 1/30 1/40 1.09(–9) 4.16(–10) 2.01(–10) 1.16(–10) Fourth order scheme (22) Sixth order scheme (33) Rate Error Rate Error Rate 1.39 1.79 1.91 3.84(–11) 2.28(–12) 4.46(–13) 1.41(–13) 4.07 4.02 4.00 2.39(–13) 3.62(–15) 3.15(–16) 5.55(–17) 6.04 6.02 6.04 Table 2 L2 -error and convergence rate for Problem 1 with h = 0.02 at T = 1 τ P–R scheme 1/10 1/20 1/30 1/40 Fourth order scheme (22) Sixth order scheme (33) Error Rate Error Rate Error Rate 7.89(–9) 1.16(–9) 4.00(–10) 1.93(–10) 2.76 2.62 2.53 3.37(–11) 2.21(–12) 4.34(–13) 1.38(–13) 3.93 4.01 3.98 2.51(–13) 3.83(–15) 3.48(–16) 5.77(–17) 6.03 5.92 6.25 Table 3 L2 -error and convergence rate for Problem 1 with τ = 0.001 at T = 1 P–R scheme h Error 1/10 1/20 1/30 1/40 2.39(–10) 5.62(–11) 2.46(–11) 1.37(–11) Fourth order scheme (22) Sixth order scheme (33) Rate Error Rate Error Rate 2.09 2.04 2.03 1.08(–12) 6.70(–14) 1.32(–14) 4.19(–15) 4.01 4.00 3.99 1.26(–14) 2.02(–16) 1.78(–17) 3.16(–18) 5.96 5.99 6.00 To verify the accuracy in time, we fix spatial grid size h = 0.02 and decrease the temporal sizes from 1/10 to 1/40 in Table 2. It shows that scheme (22) and scheme (33) achieve the expected fourth order and sixth order accuracy in time, respectively. Table 3 shows the L2 norm error at T = 1 when we fix temporal grid size τ = 0.001 and decrease the spatial grid sizes from 1/10 to 1/40. The results in Table 3 confirm that scheme (22) and scheme (33) are fourth order and sixth order accuracy in space, respectively. These computed results are in full agreement with the theoretical accuracy order. Figure 1 depicts the absolute error obtained from scheme (22) and scheme (33) with h = τ = 1/20 at T = 1. It indicates that the present schemes indeed achieve a very high accuracy on comparably coarse mesh, as compared to conventional P–R splitting methods. Problem 2  2  ∂u ∂ u ∂ 2u 1 , = + ∂t 17π 2 ∂x2 ∂y2 0 ≤ x, y ≤ 4, t > 0. The exact solution of this problem is u(x, y, t) = e–t sin(4πx) sin(πy). We design this problem to let the solution u change much faster in x direction than in y direction. To evaluate the overall convergence rate of the present two schemes for solving Problem 2, we choose h = τ = 1/10, 1/20, 1/30, 1/40 at T = 1 in Table 4. The data show that the present scheme (22) has an overall fourth order convergence rate and scheme (33) has Chen and Ge Advances in Difference Equations ( 2018) 2018:361 Page 13 of 17 Figure 1 The exact solution (a), the absolute error obtained from P–R scheme (b), present fourth order scheme (c), and present sixth order scheme (d), with h = τ = 1/20 at T = 1, Problem 1 Table 4 L2 -norm error and convergence rate for Problem 2 with h = τ at T = 1 h P–R scheme Error 1/10 1/20 1/30 1/40 1.08(–2) 1.08(–2) 5.03(–3) 2.84(–3) Fourth order scheme (22) Sixth order scheme (33) Rate Error Rate Error Rate 0.00 1.88 1.99 7.36(–3) 4.57(–4) 8.95(–5) 2.82(–5) 4.01 4.02 4.01 1.31(–3) 2.15(–5) 1.91(–6) 3.40(–7) 5.92 5.97 5.99 Table 5 L2 -norm error for Problem 2 with τ = 0.2 at T = 1 h r = τ /h 1/40 1/80 1/160 8 16 32 P–R scheme Fourth order scheme (22) Sixth order scheme (33) Error Error Error 1.85(–3) 3.07(–4) 8.47(–4) 2.94(–5) 2.97(–6) 1.32(–7) 3.39(–7) 5.02(–9) 1.46(–10) an overall sixth order convergence rate. In order to verify the unconditional stability and unconditional convergence of the present methods, we choose the different mesh ratios r, which is defined as r = τ /h at T = 1 in Table 5. It is easy to observe that the present methods are unconditionally stable. The exact solution and the absolute error obtained from the P– R scheme, the fourth order scheme, and the sixth order scheme with h = τ = 1/20 at T = 1 Chen and Ge Advances in Difference Equations ( 2018) 2018:361 Page 14 of 17 Figure 2 The exact solution (a), the absolute error obtained from P–R scheme (b), present fourth order scheme (c), and present sixth order scheme (d), with h = τ = 1/20 at T = 1, Problem 2 are plotted in Fig. 2 for Problem 2. It shows that our splitting methods can still produce very accurate solution via significantly coarse grid in the situation bypassing the accuracy barrier successfully. 5 Conclusion In this paper, we propose two types of high order splitting finite difference schemes for solving the two-dimensional parabolic equation. Based on the LOD strategy, we separate the two-dimensional equation into two one-dimensional equations to construct the new schemes. Traditional splitting accuracy barrier is bypassed naturally. Semi-discretized formulas are utilized to get high order finite difference schemes without any influence of Courant numbers. We apply the fourth order and the sixth order approximation to discrete the spatial variables, (2, 2) Padé and (3, 3) Padé approximation to discrete the temporal variable, respectively. We obtain two splitting schemes being fourth order and sixth order accuracy in both time and space, respectively. By rigorous matrix analysis, we show that all finite difference schemes involved are unconditionally stable. Two testing experiments are carried out to demonstrate the high accuracies and efficiencies of our new LOD Chen and Ge Advances in Difference Equations ( 2018) 2018:361 Page 15 of 17 schemes. It is worthy of being pointed out that the present methods can be straightforwardly extended to the three, or more, dimensional linear parabolic equations. Higher order and stable splitting methods can also be constructed in similar ways. We plan to report results from forthcoming research in this aspect in the near future. Appendix A sixth order finite difference for the one-dimensional steady diffusion equation ∂ 2u = f (x), ∂x2 0<x<1 with the boundary condition u(0) = α, u(1) = β. Let h = N1 denote the step size of space, and define xi = ih, i = 0, 1, 2, . . . , N . Based on Taylor’s series expansion of continuous functions, we get f ′′ (xi )h2 =   –fi–2 + 16fi–1 – 30fi + 16fi+1 – fi+2 + O h6 . 12 For i = 2, 3, . . . , N – 1, we obtained   ui–1 – 2ui + ui+1 1 1 = (fi–1 + 28fi + fi+1 ) + f ′′ (xi )h2 + O h6 h2 30 20   1 –fi–2 + 16fi–1 – 30fi + 16fi+1 – fi+2 = (fi–1 + 28fi + fi+1 ) + + O h6 30 240   1 97 1 1 1 fi–2 + fi–1 + fi + fi – fi+2 + O h6 . =– 240 10 120 10 240 At the grid point 1, the sixth order formula is expressed as follows:    h2  ′′ 1 1 6 u0 – 2u1 + u2 ′′ f f – 2f + f – (x ) + f (x ) = + 28f + f ) + (f 0 1 2 0 2 2 1 0 h2 30 20 5 10  6 +O h       1 h2 ′′ 3 28 3 1 3 = f0 + f1 + f2 – + – + f (x0 ) 30 50 30 25 30 50 200   h2 ′′ f (x2 ) + O h6 – 200     3 1 16 28 3 1 = f0 + f1 + + – – 30 50 2400 30 25 2400     16 3 30 1 h2 ′′ 1 + f2 – + + f3 + f4 – f (x0 ) + O h6 30 50 2400 2400 2400 200 =   121 127 1 1 h2 ′′ 3 f0 + f1 + f2 – f3 + f4 – f (x0 ) + O h6 . 32 150 1200 150 2400 200 Chen and Ge Advances in Difference Equations ( 2018) 2018:361 Page 16 of 17 Similarly, we can construct the sixth order formula at grid point N – 1: 3 121 127 1 uN–2 – 2uN–1 + uN = fN – fN–1 + fN–2 – fN–3 h2 32 150 1200 150   h2 ′′ 1 fN–4 – f (xN ) + O h6 . + 2400 200 We can see more details in Ref. [2]. For the sake of completeness, we repeated the above deduction here. Funding This work has been partially supported by the National Natural Science Foundation of China under Grant 11772165, the National Natural Science Foundation of Ningxia under Grant 2018AAC02003, and the Key Research and Development Program of Ningxia under Grant 2018BEE03007. Availability of data and materials Not applicable. Competing interests The authors declare that they have no competing interests. Authors’ contributions JC deduced the two high order difference schemes in this paper, analyzed their stability and convergence, and conducted the numerical experiments. YG presented the idea of this work and wrote this manuscript. All authors read and approved the final manuscript. Authors’ information Not applicable. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Received: 19 July 2018 Accepted: 30 September 2018 References 1. Sun, H.W., Zhang, J.: A high-order compact boundary value method for solving one-dimensional heat equations. Numer. Methods Partial Differ. Equ. 9, 846–857 (2003) 2. Lin, Y.X., Gao, J., Xiao, M.Q.: A high-order finite difference method for 1D nonhomogeneous heat equation. Numer. Methods Partial Differ. Equ. 25, 327–346 (2009) 3. Zhao, J., Dai, W., Niu, T.C.: Fourth order compact schemes of a heat conduction problem with Neumann boundary conditions. Numer. Methods Partial Differ. Equ. 23, 949–959 (2007) 4. Li, J., Chen, Y., Liu, G.: High order compact ADI methods for parabolic equations. Comput. Math. Appl. 52, 1343–1356 (2006) 5. Bialecki, B., De Frutos, J.: ADI spectral collocation methods for parabolic problems. J. Comput. Phys. 229, 5182–5193 (2010) 6. Zhou, H., Wu, Y.J., Tian, W.: Extrapolation algorithm of compact ADI approximation for two-dimensional parabolic equations. Appl. Math. Comput. 219, 2875–2884 (2012) 7. Qin, J.: The new alternating direction implicit difference methods for solving three-dimensional parabolic equations. Appl. Math. Model. 34, 890–897 (2010) 8. Wang, Y.M.: Error and extrapolation of a compact LOD method for parabolic differential equations. J. Comput. Appl. Math. 235, 1367–1382 (2011) 9. Wang, T., Wang, Y.M.: A higher-order compact LOD method and its extrapolations for nonhomogeneous parabolic differential equations. Appl. Math. Comput. 237, 512–530 (2014) 10. Wang, T.: Alternating direction finite volume element methods for 2D parabolic partial difference differential equations. Numer. Methods Partial Differ. Equ. 24, 24–40 (2008) 11. Peaceman, D.W., Rachford, H.H.: The numerical solution of parabolic and elliptic differential equation. J. Soc. Ind. Appl. Math. 3, 28–41 (1955) 12. Dehghan, M.: A new ADI technique for two-dimensional parabolic equation with an integral condition. Comput. Math. Appl. 43, 1477–1488 (2002) 13. Daoud, D.S.: On the numerical solution of multi-dimensional parabolic problem by the additive splitting up method. Appl. Math. Comput. 162, 197–210 (2005) 14. Douglas, J., Kimy, S.: Improved accuracy for locally one-dimensional methods for parabolic equations. Math. Models Methods Appl. Sci. 11, 1563–1579 (2001) 15. Dai, W., Nasar, R.: Compact finite difference scheme for solving parabolic differential equations. Numer. Methods Partial Differ. Equ. 18, 129–141 (2002) 16. Karaa, S.: An accurate LOD scheme for two-dimensional parabolic problems. Appl. Math. Comput. 27, 209–224 (2005) Chen and Ge Advances in Difference Equations ( 2018) 2018:361 Page 17 of 17 17. Zhao, J., Dai, W.Z., Zhang, S.Y.: Fourth order compact schemes for solving multi-dimensional heat conduction problems with Neumann boundary conditions. Numer. Methods Partial Differ. Equ. 24, 165–178 (2008) 18. Hiroaki, N.: First-, second-, and third-order finite-volume schemes for diffusion. J. Comput. Phys. 256, 791–805 (2014) 19. Vu, T.L., Alexander, O.: Explicit exponential Runge–Kutta methods of high order for parabolic problems. J. Comput. Appl. Math. 259, 262–361 (2014) 20. Gao, F., Chi, C.M.: Unconditionally stable difference schemes for a one-space-dimensional linear hyperbolic equation. Appl. Math. Comput. 187, 1272–1276 (2007) 21. Liu, H.W., Liu, L.B., Chen, Y.P.: A semi-discretization method based on quartic splines for solving one-space-dimensional hyperbolic equations. Appl. Math. Comput. 210, 508–514 (2009) 22. Piao, X.F., Choi, H., Kim, S.D.: A fast singly diagonally implicit Runge–Kutta method for solving 1D unsteady convection-diffusion equations. Numer. Methods Partial Differ. Equ. 30, 788–812 (2013) 23. Zhang, Y.X.: [3, 3] Padé approximation method for solving space fractional Fokker–Planck equations. Appl. Math. Lett. 35, 109–114 (2014) 24. Hao, W., Zhu, S.: Domain decomposition schemes with high-order accuracy and unconditional stability. Appl. Math. Comput. 219, 6170–6181 (2013) 25. Liu, D., Kuang, W., Tangborn, A.: High-order compact implicit difference methods for parabolic equations in geodynamo simulation. Adv. Math. Phys. 2009, Article ID 568296 (2009). https://doi.org/10.1155/2009/568296 26. Liu, D., Chen, Q., Wang, Y.: A sixth order accuracy solution to a system of nonlinear differential equations with coupled compact method. J. Comput. Eng. 2013, Article ID 432192 (2013). https://doi.org/10.1155/2013/432192 27. Sheng, Q.: Solving linear partial differential equations by exponential splitting. IMA J. Numer. Anal. 9, 199–212 (1989) 28. Hirsh, R.S.: High order accurate difference solutions of fluid mechanics problem by a compact differencing technique. J. Comput. Phys. 19, 90–109 (1975) 29. Cuyt, A.: Padé Approximants for Operators: Theory and Applications. Lecture Notes in Mathematics. Springer, Berlin (1984) 30. Smith, G.D.: Numerical Solution of Partial Differential Equations: Finite Difference Methods, 3rd edn. Oxford University Press, London (1985)