# The Splitting Crank–Nicolson Scheme with Intrinsic Parallelism for Solving Parabolic Equations

Authors: Guanyu Xue; Yunjie Gong; Hui Feng
Journal: Journal of Function Spaces
 (2020)

---

## Abstract

In this paper, a splitting Crank–Nicolson (SC-N) scheme with intrinsic parallelism is proposed for parabolic equations. The new algorithm splits the Crank–Nicolson scheme into two domain decomposition methods, each one is applied to compute the values at (n + 1)th time level by use of known numerical solutions at n-th time level, respectively. Then, the average of the above two values is chosen to be the numerical solutions at (n + 1)th time level. The new algorithm obtains accuracy of the Crank–Nicolson scheme while maintaining parallelism and unconditional stability. This algorithm can be extended to solve two-dimensional parabolic equations by alternating direction implicit (ADI) technique. Numerical experiments illustrate the accuracy and efficiency of the new algorithm.


---

## Body

## 1. Introduction

The parallel numerical schemes for parabolic equations have been studied rapidly. In 1985, Evans [1] proposed the Alternating Group Explicit (AGE) scheme for the diffusion equation by Saul’yev asymmetric schemes [2]. The Alternating Segment Crank–Nicolson (ASC-N) scheme was designed in [3]. Afterwards, the alternating segment algorithms (AGE scheme and ASC-N scheme) above became very effective methods for parabolic equations, such as convection-diffusion equation [4], dispersive equation [5–8], and fourth-order parabolic equation [9, 10]. Meanwhile, domain decomposition methods (DDMs) for the partial differential equations have been developed. Since the advantages of the low computation and communication costs at each time step, the nonoverlapping DDMs have been studied extensively [11–20]. The concept called “intrinsic parallelism” was presented in [21–23]. In 1999, the alternating difference schemes were presented, and the unconditional stability analysis was given in [24]. The unconditionally stable domain decomposition method was obtained by the alternating technique in [25–30]. Recently, the numerical methods for parabolic equations have attracted great attention of scholars [31–42]. Meanwhile, the finite volume element methods for elliptic problems are studied by Bi et al. [43, 44], and the finite volume element method for second-order hyperbolic equations is proposed by Chen et al. [45]. Obviously, the Crank–Nicolson scheme is one of the most classical methods for PDEs [46, 47]. But, there is litter work on parallel algorithms that satisfy the accuracy and stability of the Crank–Nicolson scheme.On this basis, a new parallel algorithm called splitting Crank–Nicolson scheme for parabolic equations will be presented in this paper. The idea of the new algorithm is to divide the classical Crank–Nicolson scheme into two parts, which are DDMs, each of which is used in computations at (n + 1)th time level utilizing the numerical solution at time level n. Then, the average of the above two values is chosen to be the numerical solutions at (n + 1)th time level. It can be described as follows:(1)
Split the Crank–Nicolson scheme into DDM I and DDM II by Saul’yev asymmetric schemes.(2)
DDM I is applied to compute the values at (n + 1)th time level noted as Vn+1 by use of known value Un at n-th time level. DDM II is also used to compute the values at (n + 1)th time level noted as Wn+1 by use of Un at n-th time level.(3)
The valueUn+1 = (Vn+1 + Wn+1)/2 are set as numerical solutions at (n + 1)th time level.The advantage of the SC-N scheme is spitting the C-N scheme into two parallel algorithms, which can be computed by parallel computers. Then, the average value obtained is restored to C-N scheme approximately. This paper is organized as follows: in Section2, we introduce a SC-N scheme for parabolic equations. For simplicity of presentation, we focus on a model problem, namely, one-dimensional parabolic equations. The new algorithm and detailed presentations are given. Then, we extend the new algorithm to solve two-dimensional parabolic equations by ADI technique. Finally, numerical experiments illustrated the accuracy of SC-N scheme is approximate to the C-N scheme and the new algorithm is efficient.


## 2. Algorithm Presentation

Considering the model problem of one-dimensional parabolic equations,(1)∂u∂t−a∂2u∂x2=fx,t,x∈0,l,t∈0,T,with the initial and boundary conditions(2)ux,0=u0x,x∈0,l,u0,t=g0t,ul,t=glt,t∈0,T,where a > 0 is a constant.Leth and τ be the spatial and temporal step sizes, respectively. Denote xj = jh, j = 0, 1, …, m and tn = nτ, n = 0, 1, …, N. Let ujn be the approximate solution at (xj, tn). u (x, t) represents the exact solution of (1).The well-known Crank–Nicolson scheme can be written as(3)−ar2uj−1n+1+1+arujn+1−ar2uj+1n+1=ar2uj−1n+1−arujn+ar2uj+1n+τfjn,where r = τ/h2.Letμ = ar, and (3) can be written as the matrix form(4)AUn+1=BUn+Fn,where Un=u1n,u2n,…,um−1nT and Fn=μ/2u0n+u0n+1+τf1n,τf2n,…,τfm−2n,μ/2umn+umn+1+τfm−1nT.The matricesA and B are as follows:(5)A=1+μ−μ2−μ21+μ−μ2⋱⋱⋱−μ21+μ−μ2−μ21+μ,B=1−μμ2μ21−μμ2⋱⋱⋱μ21−μμ2μ21−μ.Splitting Crank–Nicolson scheme (4), we obtain(6)A1+A2Un+1=B1+B2Un+2Fn,where A1 and A2 are block diagonal matrices, respectively.(7)A=A1+A22,B=B1+B22.V
n+1 and Wn+1 can be defined as follows:(8)A1Vn+1≜B1Un+Fn,(9)A2Wn+1≜B2Un+Fn,where Vn+1 and Wn+1 approximate to Un+1, respectively. Then,(10)Un+1≈Vn+1+Wn+12.Scheme (8) is DDM I, and scheme (9) is called DDM II. The two DDMs aforementioned are suitable for parallel computing.In order to construct two domain decomposition methods, DDM I and DDM II, at (n + 1)th time level, we need to consider four forms of Saul’yev asymmetric difference schemes corresponding to (3) as follows:(11)−ar2uj−1n+1+1+3ar2ujn+1−aruj+1n+1=ar2uj−1n+1−ar2ujn+τfjn,(12)−aruj−1n+1+1+3ar2ujn+1−ar2uj+1n+1=1−ar2ujn+ar2uj+1n+τfjn,(13)1+ar2ujn+1−ar2uj+1n+1=aruj−1n+1−3ar2ujn+ar2uj+1n+τfjn,(14)−ar2uj−1n+1+1+ar2ujn+1=ar2uj−1n+1−3ar2ujn+aruj+1n+τfjn.The flow chart of schemes (11)–(14) are displayed in Figure 1.Figure 1
The diagram of Saul’yev asymmetric difference schemes. (a) Scheme (11). (b) Scheme (12). (c) Scheme (13). (d) Scheme (14).
(a)
(b)
(c)
(d)Assumem − 1 = 6K, where K is a positive integer.


### 2.1. DDM I

For the valuesu1n+1,u2n+1,andu3n+1, we use the formulas as follows:(15)−ar2u0n+1+1+aru1n+1−ar2u2n+1=ar2u0n+1−aru1n+ar2u2n+τf1n,−ar2u1n+1+1+aru2n+1−ar2u3n+1=ar2u1n+1−aru2n+ar2u3n+τf2n,−ar2u2n+1+1+ar2u3n+1=ar2u2n+1−3ar2u3n+aru4n+τf3n.Find the valuesu6k−2n+1,u6k−1n+1,u6kn+1,u6k+1n+1,u6k+2n+1,u6k+3n+1 by using the following formulas (k = 1, 2, …, K − 1):(16)1+ar2u6k−2n+1−ar2u6k−1n+1=aru6k−3n+1−3ar2u6k−2n+ar2u6k−1n+τf6k−2n,−ar2u6k−2n+1+1+aru6k−1n+1−ar2u6kn+1=ar2u6k−2n+1−aru6k−1n+ar2u6kn+τf6k−1n,−ar2u6k−1n+1+1+3ar2u6kn+1−aru6k+1n+1=ar2u6k−1n+1−ar2u6kn+τf6kn,−aru6kn+1+1+3ar2u6k+1n+1−ar2u6k+2n+1=1−ar2u6k+1n+ar2u6k+2n+τf6k+1n,−ar2u6k+1n+1+1+aru6k+2n+1−ar2u6k+3n+1=ar2u6k+1n+1−aru6k+2n+ar2u6k+3n+τf6k+2n,−ar2u6k+2n+1+1+ar2u6k+3n+1=ar2u6k+2n+1−3ar2u6k+3n+aru6k+4n+τf6k+3n.Obviously, each subdomain contains 6 nodes which can be computed by (16) independently.For the valuesum−3n+1,um−2n+1,um−1n+1 by using the formulas as follows:(17)1+ar2um−3n+1−ar2um−2n+1=arum−4n+1−3ar2um−3n+ar2um−2n+τfm−3n,−ar2um−3n+1+1+arum−2n+1−ar2um−1n+1=ar2um−3n+1−arum−2n+ar2um−1n+τfm−2n,−ar2um−2n+1+1+arum−1n+1−ar2umn+1=ar2um−2n+1−arum−1n+ar2umn+τfm−1n.DDM I can be written as the matrix form(18)A1Un+1=B1Un+Fn,where Un=u1n,u2n,⋯,um−1nT, A1 = (I + μG1), and B1 = (I − μG2).The matricesG1 and G2 are block diagonal matrices which are as follows:(19)G1=P1Q⋱QP2,G2=Q1Q⋱QQ2,where(20)P1=1−120−121−120−1212,P2=12−120−121−120−121,Q=12−12−121−12−1232−1−132−12−121−12−1212,Q1=1−12−121−12−1232−1−132−12−121−12−1212,Q2=12−12−121−12−1232−1−132−12−121−12−121.Each block matrix system (i.e., each subdomain) can be solved independently. It is evident that DDM I (18) has intrinsic parallelism.


### 2.2. DDM II

For the valuesu1n+1,u2n+1,…,u6n+1 by using the following formulas:(21)−ar2u0n+1+1+aru1n+1−ar2u2n+1=ar2u0n+1−aru1n+ar2u2n+τf1n,−ar2u1n+1+1+aru2n+1−ar2u3n+1=ar2u1n+1−aru2n+ar2u3n+τf2n,−ar2u2n+1+1+3ar2u3n+1−aru4n+1=ar2u2n+1−ar2u3n+τf3n,−aru3n+1+1+3ar2u4n+1−ar2u5n+1=1−ar2u4n+ar2u5n+τf4n,−ar2u4n+1+1+aru5n+1−ar2u6n+1=ar2u4n+1−aru5n+ar2u6n+τf5n,−ar2u5n+1+1+ar2u6n+1=ar2u5n+1−3ar2u6n+aru7n+τf6n.Find the valuesu6k+1n+1,u6k+2n+1,u6k+3n+1,u6k+4n+1,u6k+5n+1,u6k+6n+1 by using the formulas as follows (k = 1, 2, …, K − 2):(22)1+ar2u6k+1n+1−ar2u6k+2n+1=aru6kn+1−3ar2u6k+1n+ar2u6k+2n+τf6k+1n,−ar2u6k+1n+1+1+aru6k+2n+1−ar2u6k+3n+1=ar2u6k+1n+1−aru6k+2n+ar2u6k+3n+τf6k+2n,−ar2u6k+2n+1+1+3ar2u6k+3n+1−aru6k+4n+1=ar2u6k+2n+1−ar2u6k+3n+τf6k+3n,−aru6k+3n+1+1+3ar2u6k+4n+1−ar2u6k+5n+1=1−ar2u6k+4n+ar2u6k+5n+τf6k+4n,−ar2u6k+4n+1+1+aru6k+5n+1−ar2u6k+6n+1=ar2u6k+4n+1−aru6k+5n+ar2u6k+6n+τf6k+5n,−ar2u6k+5n+1+1+ar2u6k+6n+1=ar2u6k+5n+1−3ar2u6k+6n+aru6k+7n+τf6k+6n.Obviously, each subdomain contains 6 nodes which can be computed by (22) independently.Find the valuesum−6n+1,um−5n+1,…,um−1n+1 by using the formulas as follows:(23)1+ar2um−6n+1−ar2um−5n+1=arum−7n+1−3ar2um−6n+ar2um−5n+τfm−6n,−ar2um−6n+1+1+arum−5n+1−ar2um−4n+1=ar2um−6n+1−arum−5n+ar2um−4n+τfm−5n,−ar2um−5n+1+1+3ar2um−4n+1−arum−3n+1=ar2um−5n+1−ar2um−4n+τfm−4n,−arum−4n+1+1+3ar2um−3n+1−ar2um−2n+1=1−ar2um−3n+ar2um−2n+τfm−3n,−ar2um−3n+1+1+arum−2n+1−ar2um−1n+1=ar2um−3n+1−arum−2n+ar2um−1n+τfm−2n,−ar2um−2n+1+1+arum−1n+1−ar2umn+1=ar2um−2n+1−arum−1n+ar2umn+τfm−1n.DDM II can be written as the matrix form(24)A2Un+1=B2Un+Fn,where A2 = (I + μG2) and B2 = (I − μG1).The matricesG1 and G2 are block diagonal matrices which are as follows:(25)G1=P1Q⋱QP2,G2=Q1Q⋱QQ2.Each block matrix system (i.e., each subdomain) can be solved independently. It is evident that DDM II (24) has intrinsic parallelism.Schemes (15)–(17) and schemes (21)–(23) construct two domain decomposition methods (18) and (24), respectively. The corresponding algorithm can be described in Algorithm 1. The matrix form of Algorithm 1 can be written as follows:(26)A1Vn+1=B1Un+Fn,A2Wn+1=B2Un+Fn,Un+1=12Vn+1+Wn+1,n=0,1,2,…,where Un=u1n,u2n,…,um−1nT, Vn=v1n,v2n,…,vm−1nT, and Wn=w1n,w2n,…,wm−1nT.Algorithm 1:The SC-N scheme for one-dimensional parabolic equations.
Require: Initialization U0 (xj) ⟵ u0 (xj).
forn = 0, 1, ⋯, Ndo
forj = 0, 1, ⋯, mdo
Solve the valuesVjn+1 by using DDM I (18).
Solve the valuesWjn+1 by using DDM II (24).
The average of two values will be calculated, i.e.Ujn+1=1/2Vjn+1+Wjn+1.
end for
end for
Ensure: Output UN (xj).Remark 1.
Obviously, C-N scheme (4) is equal to A1+A2/2Un+1=B1+B2/2Un+Fn.Remark 2.
Because the SC-N scheme is derived from the C-N scheme, it also has the properties of the C-N scheme, i.e., the SC-N scheme is unconditionally stable and maintains the second-order numerical accuracyO (τ2 + h2).


## 2.1. DDM I

For the valuesu1n+1,u2n+1,andu3n+1, we use the formulas as follows:(15)−ar2u0n+1+1+aru1n+1−ar2u2n+1=ar2u0n+1−aru1n+ar2u2n+τf1n,−ar2u1n+1+1+aru2n+1−ar2u3n+1=ar2u1n+1−aru2n+ar2u3n+τf2n,−ar2u2n+1+1+ar2u3n+1=ar2u2n+1−3ar2u3n+aru4n+τf3n.Find the valuesu6k−2n+1,u6k−1n+1,u6kn+1,u6k+1n+1,u6k+2n+1,u6k+3n+1 by using the following formulas (k = 1, 2, …, K − 1):(16)1+ar2u6k−2n+1−ar2u6k−1n+1=aru6k−3n+1−3ar2u6k−2n+ar2u6k−1n+τf6k−2n,−ar2u6k−2n+1+1+aru6k−1n+1−ar2u6kn+1=ar2u6k−2n+1−aru6k−1n+ar2u6kn+τf6k−1n,−ar2u6k−1n+1+1+3ar2u6kn+1−aru6k+1n+1=ar2u6k−1n+1−ar2u6kn+τf6kn,−aru6kn+1+1+3ar2u6k+1n+1−ar2u6k+2n+1=1−ar2u6k+1n+ar2u6k+2n+τf6k+1n,−ar2u6k+1n+1+1+aru6k+2n+1−ar2u6k+3n+1=ar2u6k+1n+1−aru6k+2n+ar2u6k+3n+τf6k+2n,−ar2u6k+2n+1+1+ar2u6k+3n+1=ar2u6k+2n+1−3ar2u6k+3n+aru6k+4n+τf6k+3n.Obviously, each subdomain contains 6 nodes which can be computed by (16) independently.For the valuesum−3n+1,um−2n+1,um−1n+1 by using the formulas as follows:(17)1+ar2um−3n+1−ar2um−2n+1=arum−4n+1−3ar2um−3n+ar2um−2n+τfm−3n,−ar2um−3n+1+1+arum−2n+1−ar2um−1n+1=ar2um−3n+1−arum−2n+ar2um−1n+τfm−2n,−ar2um−2n+1+1+arum−1n+1−ar2umn+1=ar2um−2n+1−arum−1n+ar2umn+τfm−1n.DDM I can be written as the matrix form(18)A1Un+1=B1Un+Fn,where Un=u1n,u2n,⋯,um−1nT, A1 = (I + μG1), and B1 = (I − μG2).The matricesG1 and G2 are block diagonal matrices which are as follows:(19)G1=P1Q⋱QP2,G2=Q1Q⋱QQ2,where(20)P1=1−120−121−120−1212,P2=12−120−121−120−121,Q=12−12−121−12−1232−1−132−12−121−12−1212,Q1=1−12−121−12−1232−1−132−12−121−12−1212,Q2=12−12−121−12−1232−1−132−12−121−12−121.Each block matrix system (i.e., each subdomain) can be solved independently. It is evident that DDM I (18) has intrinsic parallelism.


## 2.2. DDM II

For the valuesu1n+1,u2n+1,…,u6n+1 by using the following formulas:(21)−ar2u0n+1+1+aru1n+1−ar2u2n+1=ar2u0n+1−aru1n+ar2u2n+τf1n,−ar2u1n+1+1+aru2n+1−ar2u3n+1=ar2u1n+1−aru2n+ar2u3n+τf2n,−ar2u2n+1+1+3ar2u3n+1−aru4n+1=ar2u2n+1−ar2u3n+τf3n,−aru3n+1+1+3ar2u4n+1−ar2u5n+1=1−ar2u4n+ar2u5n+τf4n,−ar2u4n+1+1+aru5n+1−ar2u6n+1=ar2u4n+1−aru5n+ar2u6n+τf5n,−ar2u5n+1+1+ar2u6n+1=ar2u5n+1−3ar2u6n+aru7n+τf6n.Find the valuesu6k+1n+1,u6k+2n+1,u6k+3n+1,u6k+4n+1,u6k+5n+1,u6k+6n+1 by using the formulas as follows (k = 1, 2, …, K − 2):(22)1+ar2u6k+1n+1−ar2u6k+2n+1=aru6kn+1−3ar2u6k+1n+ar2u6k+2n+τf6k+1n,−ar2u6k+1n+1+1+aru6k+2n+1−ar2u6k+3n+1=ar2u6k+1n+1−aru6k+2n+ar2u6k+3n+τf6k+2n,−ar2u6k+2n+1+1+3ar2u6k+3n+1−aru6k+4n+1=ar2u6k+2n+1−ar2u6k+3n+τf6k+3n,−aru6k+3n+1+1+3ar2u6k+4n+1−ar2u6k+5n+1=1−ar2u6k+4n+ar2u6k+5n+τf6k+4n,−ar2u6k+4n+1+1+aru6k+5n+1−ar2u6k+6n+1=ar2u6k+4n+1−aru6k+5n+ar2u6k+6n+τf6k+5n,−ar2u6k+5n+1+1+ar2u6k+6n+1=ar2u6k+5n+1−3ar2u6k+6n+aru6k+7n+τf6k+6n.Obviously, each subdomain contains 6 nodes which can be computed by (22) independently.Find the valuesum−6n+1,um−5n+1,…,um−1n+1 by using the formulas as follows:(23)1+ar2um−6n+1−ar2um−5n+1=arum−7n+1−3ar2um−6n+ar2um−5n+τfm−6n,−ar2um−6n+1+1+arum−5n+1−ar2um−4n+1=ar2um−6n+1−arum−5n+ar2um−4n+τfm−5n,−ar2um−5n+1+1+3ar2um−4n+1−arum−3n+1=ar2um−5n+1−ar2um−4n+τfm−4n,−arum−4n+1+1+3ar2um−3n+1−ar2um−2n+1=1−ar2um−3n+ar2um−2n+τfm−3n,−ar2um−3n+1+1+arum−2n+1−ar2um−1n+1=ar2um−3n+1−arum−2n+ar2um−1n+τfm−2n,−ar2um−2n+1+1+arum−1n+1−ar2umn+1=ar2um−2n+1−arum−1n+ar2umn+τfm−1n.DDM II can be written as the matrix form(24)A2Un+1=B2Un+Fn,where A2 = (I + μG2) and B2 = (I − μG1).The matricesG1 and G2 are block diagonal matrices which are as follows:(25)G1=P1Q⋱QP2,G2=Q1Q⋱QQ2.Each block matrix system (i.e., each subdomain) can be solved independently. It is evident that DDM II (24) has intrinsic parallelism.Schemes (15)–(17) and schemes (21)–(23) construct two domain decomposition methods (18) and (24), respectively. The corresponding algorithm can be described in Algorithm 1. The matrix form of Algorithm 1 can be written as follows:(26)A1Vn+1=B1Un+Fn,A2Wn+1=B2Un+Fn,Un+1=12Vn+1+Wn+1,n=0,1,2,…,where Un=u1n,u2n,…,um−1nT, Vn=v1n,v2n,…,vm−1nT, and Wn=w1n,w2n,…,wm−1nT.Algorithm 1:The SC-N scheme for one-dimensional parabolic equations.
Require: Initialization U0 (xj) ⟵ u0 (xj).
forn = 0, 1, ⋯, Ndo
forj = 0, 1, ⋯, mdo
Solve the valuesVjn+1 by using DDM I (18).
Solve the valuesWjn+1 by using DDM II (24).
The average of two values will be calculated, i.e.Ujn+1=1/2Vjn+1+Wjn+1.
end for
end for
Ensure: Output UN (xj).Remark 1.
Obviously, C-N scheme (4) is equal to A1+A2/2Un+1=B1+B2/2Un+Fn.Remark 2.
Because the SC-N scheme is derived from the C-N scheme, it also has the properties of the C-N scheme, i.e., the SC-N scheme is unconditionally stable and maintains the second-order numerical accuracyO (τ2 + h2).


## 3. Extension to Two-Dimensional Parabolic Equations

In this section, we will extend Algorithm1 to solve two-dimensional parabolic equations:(27)∂u∂t−∂∂xa∂u∂x−∂∂yb∂u∂y=fx,y,t,x,y∈Ω,t∈0,T,(28)ux,y,0=u0x,y,x,y∈Ω,(29)ux,y,t=0,x,y∈∂Ω,t∈0,T,where domain Ω ∈ (0, Lx) × (0, Ly) and a > 0 and b > 0 are diffusion coefficients.Letui,jn be the approximate solution at (xi, yj, tn), and u (x, y, t) represents the exact solution of (27). With the same time and space discretization of Algorithm 1, we obtain its extended algorithm by alternating direction implicit (ADI) technique [48] for equations (27)–(29).


### 3.1.x-Direction

Letr1 = aτ/(2h2) and r2 = bτ/(2h2), and the matrix form of the SC-N scheme in x-direction can be written as follows:(30)I+r1G1V1n+1/2=I−r1G2Un+b1n,I+r1G2W1n+1/2=I−r1G1Un+b1n,Un+1/2=12V1n+1/2+W1n+1/2,where(31)b1n=r2u0,jn+r1u0,jn+1/2/2+r2u1,j−1n−2u1,jn+u1,j+1n+τf1,jn/2r2u2,j−1n−2u2,jn+u2,j+1n+τf2,jn/2⋮r2um−2,j−1n−2um−2,jn+um−2,j+1n+τfm−2,jn/2r2um,jn+r1um,jn+1/2/2+r2um−1,j−1n−2um−1,jn+um−1,j+1n+τfm−1,jn/2,Un+1/2=u1,jn+1/2,u2,jn+1/2,…,um−1,jn+1/2T,V1n+1/2=v1,jn+1/2,v2,jn+1/2,…,vm−1,jn+1/2T,W1n+1/2=w1,jn+1/2,w2,jn+1/2,…,wm−1,jn+1/2T,j=1,2,…,m−1.


### 3.2.y-Direction

The matrix form of the SC-N scheme iny-direction can be written as follows:(32)I+r2G1V2n+1=I−r2G2Un+1/2+b2n+1/2,I+r2G2W2n+1=I−r2G1Un+1/2+b2n+1/2,Un+1=12V2n+1+W2n+1,where(33)b2n+1/2=r2ui,0n+1/2+r1ui,0n+1/2+r1ui−1,1n+1−2ui,1n+1+ui+1,1n+1+τfi,1n+1/2/2r1ui−1,2n+1−2ui,2n+1+ui+1,2n+1+τfi,2n+1/2/2⋮r1ui−1,m−2n+1−2ui,m−2n+1+ui+1,m−2n+1+τfi,m−2n+1/2/2r2ui,mn+1/2+r1ui,mn+1/2+r1ui−1,m−1n+1−2ui,m−1n+1+ui+1,m−1n+1+τfi,m−1n+1/2/2,Un+1=ui,1n+1,ui,2n+1,…,ui,m−1n+1T,V2n+1=vi,1n+1,vi,2n+1,…,vi,m−1n+1T,W2n+1=wi,1n+1,wi,2n+1,…,wi,m−1n+1T,i=1,2,…,m−1.The corresponding algorithm can be described in Algorithm2.Algorithm 2:The SC-N scheme for two-dimensional parabolic equations.
Require: Initialization U0 (xi, yj) ⟵ u0 (xi, yj).
forn = 0, 1, ⋯, Ndo
fori = 0, 1, ⋯, mdo
forj = 0, 1, ⋯, mdo
Solve the valuesUi,jn+1/2 by using the scheme (30).
Solve the valuesUi,jn+1 again by using the scheme (32).
end for
end for
end for
Ensure: Output UN (xi, yj).Similar to Algorithm1, it is obvious that Algorithm 2 has unconditional stability and parallelism.Remark 3.
In fact,I+r1G1+I+r1G2/2Un+1/2=I−r1G1+I−r1G2/2Un+b1n is the C-N scheme in x-direction, and I+r2G1+I+r2G2/2Un+1=I−r2G1+I−r2G2/2Un+1/2+b2n+1/2 is the C-N scheme in y-direction.Remark 4.
In Algorithm2, the domain is divided into many subdomains by using two DDMs. In each time interval, we first solve the values along x-direction by (30) at the half-time step and then solve the values along y-direction by (32) at the next half-time step. Schemes (30) and (32) lead to block diagonal algebraic systems that can be solved independently. So, Algorithm 2 not only is suitable for parallel computation but also maintains the accuracy. Based on the advantage of ADI technique, Algorithm 2 reduces computational complexities. Though it is developed for two-dimensional problems, Algorithm 2 can be easily extended to solve high-dimensional parabolic equations.


## 3.1.x-Direction

Letr1 = aτ/(2h2) and r2 = bτ/(2h2), and the matrix form of the SC-N scheme in x-direction can be written as follows:(30)I+r1G1V1n+1/2=I−r1G2Un+b1n,I+r1G2W1n+1/2=I−r1G1Un+b1n,Un+1/2=12V1n+1/2+W1n+1/2,where(31)b1n=r2u0,jn+r1u0,jn+1/2/2+r2u1,j−1n−2u1,jn+u1,j+1n+τf1,jn/2r2u2,j−1n−2u2,jn+u2,j+1n+τf2,jn/2⋮r2um−2,j−1n−2um−2,jn+um−2,j+1n+τfm−2,jn/2r2um,jn+r1um,jn+1/2/2+r2um−1,j−1n−2um−1,jn+um−1,j+1n+τfm−1,jn/2,Un+1/2=u1,jn+1/2,u2,jn+1/2,…,um−1,jn+1/2T,V1n+1/2=v1,jn+1/2,v2,jn+1/2,…,vm−1,jn+1/2T,W1n+1/2=w1,jn+1/2,w2,jn+1/2,…,wm−1,jn+1/2T,j=1,2,…,m−1.


## 3.2.y-Direction

The matrix form of the SC-N scheme iny-direction can be written as follows:(32)I+r2G1V2n+1=I−r2G2Un+1/2+b2n+1/2,I+r2G2W2n+1=I−r2G1Un+1/2+b2n+1/2,Un+1=12V2n+1+W2n+1,where(33)b2n+1/2=r2ui,0n+1/2+r1ui,0n+1/2+r1ui−1,1n+1−2ui,1n+1+ui+1,1n+1+τfi,1n+1/2/2r1ui−1,2n+1−2ui,2n+1+ui+1,2n+1+τfi,2n+1/2/2⋮r1ui−1,m−2n+1−2ui,m−2n+1+ui+1,m−2n+1+τfi,m−2n+1/2/2r2ui,mn+1/2+r1ui,mn+1/2+r1ui−1,m−1n+1−2ui,m−1n+1+ui+1,m−1n+1+τfi,m−1n+1/2/2,Un+1=ui,1n+1,ui,2n+1,…,ui,m−1n+1T,V2n+1=vi,1n+1,vi,2n+1,…,vi,m−1n+1T,W2n+1=wi,1n+1,wi,2n+1,…,wi,m−1n+1T,i=1,2,…,m−1.The corresponding algorithm can be described in Algorithm2.Algorithm 2:The SC-N scheme for two-dimensional parabolic equations.
Require: Initialization U0 (xi, yj) ⟵ u0 (xi, yj).
forn = 0, 1, ⋯, Ndo
fori = 0, 1, ⋯, mdo
forj = 0, 1, ⋯, mdo
Solve the valuesUi,jn+1/2 by using the scheme (30).
Solve the valuesUi,jn+1 again by using the scheme (32).
end for
end for
end for
Ensure: Output UN (xi, yj).Similar to Algorithm1, it is obvious that Algorithm 2 has unconditional stability and parallelism.Remark 3.
In fact,I+r1G1+I+r1G2/2Un+1/2=I−r1G1+I−r1G2/2Un+b1n is the C-N scheme in x-direction, and I+r2G1+I+r2G2/2Un+1=I−r2G1+I−r2G2/2Un+1/2+b2n+1/2 is the C-N scheme in y-direction.Remark 4.
In Algorithm2, the domain is divided into many subdomains by using two DDMs. In each time interval, we first solve the values along x-direction by (30) at the half-time step and then solve the values along y-direction by (32) at the next half-time step. Schemes (30) and (32) lead to block diagonal algebraic systems that can be solved independently. So, Algorithm 2 not only is suitable for parallel computation but also maintains the accuracy. Based on the advantage of ADI technique, Algorithm 2 reduces computational complexities. Though it is developed for two-dimensional problems, Algorithm 2 can be easily extended to solve high-dimensional parabolic equations.


## 4. Numerical Experiments

To illustrate the efficiency of the SC-N scheme for parabolic equations, we will compare the accuracy of the new algorithm with the existing method.Example 1.
(34)
∂
u
∂
t
−
∂
2
u
∂
x
2
=
1
+
π
2
e
t
sin
π
x
,
x
∈
0,1
,
t
∈
0,1
,
u
x
,
0
=
sin
π
x
,
x
∈
0,1
,
u
0
,
t
=
u
1
,
t
=
0
,
t
∈
0,1
.
The exact solution of Example1 is(35)ux,t=etsinπx.
Firstly, we examine the convergence rate of Algorithm1. We divide the mesh point into many segments, such as K = 3, K = 4, K = 5, and K = 6. We calculate errors L∞=uxj,tn−ujn∞ taking τ = 0.001, and the following rate of convergence in the space is(36)rate≈logLh1∞/Lh2∞logh1/h2.
Clearly, the errors appear to be of orderO (h2) in Table 1.
Next, we present the error results of the SC-N scheme in terms of the absolute errors and the relative errors, where the absolute error (A. E.) is defined by(37)ejn=ujn−uxj,tn,and the relative error (R. E.) is defined by(38)Ejn=ejnuxj,tn×100%.
Tables2 and 3 display the absolute errors and the relative errors obtained by presented Algorithm 1 for h = 1/19 (i.e., K = 3) and h = 1/25 (i.e., K = 4) at t = 0.2, t = 0.4, and t = 0.8 when taking r = 1.5 (r = τ/h2). From Tables 2 and 3, it is obvious that our algorithm has high accuracy.
We compare Algorithm1 with the ASC-N scheme in [3] and C-N scheme (3) by the maximum errors for h = 1/19 (i.e., K = 3) and h = 1/25 (i.e., K = 4) at different times t = 0.2, t = 0.4, t = 0.6, and t = 0.8. With the increasing of computation time, the errors of the ASC-N scheme in [3] increase more than those of Algorithm 1 for different r (r = τ/h2) in Tables 4 and 5. Moreover, the results show that Algorithm 1 can achieve the same accuracy as the classic C-N scheme while maintaining parallelism. We consider an example for h = 1/121 (i.e., K = 20) with large grid ratio r = 15 (r = τ/h2). Table 6 shows that Algorithm 1 has a better accuracy than two others, and it is indicated that Algorithm 1 is stable.Table 1
Convergence rate of Algorithm1 for h at t = 0.4.
h
L
∞ error
Rate
1/19 (K = 3)
4.998746e − 4
—
1/25 (K = 4)
3.085232e − 4
2.0167
1/31 (K = 5)
2.002844e − 4
1.9955
1/37 (K = 6)
1.412289e − 4
1.9903Table 2
The absolute errors and relative errors of numerical solutions to Example1 for h = 1/19 (i.e., K = 3).
x
j
Algorithm 1 (t = 0.2)
Algorithm 1 (t = 0.4)
Algorithm 1 (t = 0.8)
A.E.
R. E.
A. E.
R. E.
A. E.
R. E.
0.11
0.3731e − 3
1.7559e − 2
0.5021e − 3
1.9312e − 2
0.7577e − 3
1.9534e − 2
0.21
0.6184e − 3
1.5420e − 2
0.8428e − 3
1.7176e − 2
1.2739e − 3
1.7399e − 2
0.31
0.8750e − 3
1.5999e − 2
1.1881e − 3
1.7754e − 2
1.7950e − 3
1.7976e − 2
0.42
0.9753e − 3
1.5409e − 2
1.3294e − 3
1.7165e − 2
2.0093e − 3
1.7388e − 2
0.52
1.0066e − 3
1.5469e − 2
1.3715e − 3
1.7225e − 2
2.0729e − 3
1.7447e − 2
0.63
0.9135e − 3
1.5280e − 2
1.2463e − 3
1.7037e − 2
1.8839e − 3
1.7259e − 2
0.73
0.7582e − 3
1.5778e − 2
1.0310e − 3
1.7534e − 2
1.5579e − 3
1.7757e − 2
0.84
0.5469e − 3
1.7563e − 2
0.7361e − 3
1.9316e − 2
1.1109e − 3
1.9538e − 2
0.95
0.1891e − 3
1.7556e − 2
0.2545e − 3
1.9309e − 2
0.3841e − 3
1.9532e − 2Table 3
The absolute errors and relative errors of numerical solutions to Example1 for h = 1/25 (i.e., K = 4).
x
j
Algorithm 1 (t = 0.2)
Algorithm 1 (t = 0.4)
Algorithm 1 (t = 0.8)
A.E.
R. E.
A. E.
R. E.
A. E.
R. E.
0.04
0.5579e − 4
9.0292e − 3
0.7499e − 4
9.9277e − 3
1.1317e − 4
1.0041e − 2
0.12
1.6391e − 4
9.0313e − 3
2.2031e − 4
9.9298e − 3
3.3247e − 4
1.0043e − 2
0.24
2.7890e − 4
8.2703e − 3
3.7803e − 4
9.1694e − 3
5.7102e − 4
0.9283e − 2
0.36
3.5674e − 4
8.0054e − 3
4.8512e − 4
8.9048e − 3
7.3305e − 4
0.9018e − 2
0.48
3.8737e − 4
7.8818e − 3
5.2760e − 4
8.7813e − 3
7.9739e − 4
0.8895e − 2
0.60
3.6642e − 4
7.8242e − 3
4.9945e − 4
8.7237e − 3
7.5491e − 4
0.8837e − 2
0.72
2.9685e − 4
7.8239e − 3
4.0462e − 4
8.7235e − 3
6.1158e − 4
0.8837e − 2
0.84
1.8851e − 4
7.9457e − 3
2.5655e − 4
8.8452e − 3
3.8769e − 4
0.8959e − 2
0.96
0.5579e − 4
9.0292e − 3
0.7499e − 4
9.9277e − 3
1.1317e − 4
1.0041e − 2Table 4
Comparisons by the maximum errors to Example1 for h = 1/19 (i.e., K = 3).
Methods
r
t = 0.2
t = 0.4
t = 0.6
t = 0.8
Algorithm1
0.36
3.5539e − 4
4.8387e − 4
5.9795e − 4
7.3131e − 4
ASC-N scheme [3]
4.7468e − 4
7.1743e − 4
1.0705e − 3
1.5969e − 3
C-N scheme
3.5468e − 4
4.7723e − 4
5.7998e − 4
7.0232e − 4
Algorithm1
0.72
5.7866e − 4
7.8813e − 4
9.7397e − 4
1.1912e − 3
ASC-N scheme [3]
7.3439e − 4
1.1100e − 3
1.6563e − 3
2.4709e − 3
C-N scheme
5.7499e − 4
7.8239e − 4
9.6534e − 4
1.0928e − 3
Algorithm1
1.45
1.0066e − 3
1.3715e − 3
1.6949e − 3
2.0729e − 3
ASC-N scheme [3]
1.0999e − 3
1.6591e − 3
2.4754e − 3
3.6929e − 3
C-N scheme
9.9823e − 4
1.3421e − 3
1.6034e − 3
1.8807e − 3
Algorithm1
1.80
1.2153e − 3
1.6560e − 3
2.0465e − 3
2.5029e − 3
ASC-N scheme [3]
1.2843e − 3
1.9356e − 3
2.8880e − 3
4.3084e − 3
C-N scheme
1.2031e − 3
1.6423e − 3
1.9034e − 3
2.2222e − 3Table 5
Comparisons by the maximum errors to Example1 for h = 1/25 (i.e., K = 4).
Methods
r
t = 0.2
t = 0.4
t = 0.6
t = 0.8
Algorithm1
0.63
2.2660e − 4
3.0852e − 4
3.8124e − 4
4.6627e − 4
ASC-N scheme [3]
2.8997e − 4
4.3826e − 4
6.5392e − 4
9.7553e − 4
C-N scheme
2.2632e − 4
3.0324e − 4
3.7834e − 4
4.5342e − 4
Algorithm1
0.94
3.0760e − 4
4.1963e − 4
5.1826e − 4
6.3353e − 4
ASC-N scheme [3]
4.5949e − 4
6.9190e − 4
8.3264e − 4
1.5416e − 3
C-N scheme
3.0642e − 4
4.1823e − 4
5.1654e − 4
6.1278e − 4
Algorithm1
1.35
3.8737e − 4
5.2761e − 4
6.5199e − 4
7.9740e − 4
ASC-N scheme [3]
4.2810e − 4
6.4710e − 4
9.6552e − 4
1.4404e − 3
C-N scheme
3.8535e − 4
5.2252e − 4
6.4721e − 4
7.7627e − 4
Algorithm1
1.76
4.6541e − 4
6.3397e − 4
7.8343e − 4
9.5815e − 4
ASC-N scheme [3]
4.6674e − 4
7.0553e − 4
1.0527e − 3
1.5704e − 3
C-N scheme
4.5121e − 4
6.2937e − 4
7.6320e − 4
9.3729e − 4Table 6
Maximum error comparison forr = 15 and h = 1/121 (i.e., K = 20).
Algorithm1
ASC-N method [3]
AGE method [1]
t = 0.2
1.3994e − 4
1.4759e − 4
1.2717e − 3
t = 0.4
3.2238e − 4
2.8034e − 4
3.2242e − 3
t = 0.6
5.0639e − 4
6.5004e − 4
6.2125e − 3
t = 0.8
6.9237e − 4
8.3824e − 4
9.2284e − 3Example 2.
(39)
∂
u
∂
t
−
∂
2
u
∂
x
2
−
∂
2
u
∂
y
2
=
f
x
,
y
,
t
,
x
,
y
∈
Ω
,
t
∈
0,1
,
u
x
,
y
,
0
=
sin
π
x
sin
π
y
,
x
,
y
∈
Ω
,
u
0
,
y
,
t
=
u
1
,
y
,
t
=
0
,
y
∈
0,1
,
t
∈
0,1
,
u
x
,
0
,
t
=
u
x
,
1
,
t
=
0
,
x
∈
0,1
,
t
∈
0,1
,where the domain is Ω = (0, 1) × (0, 1), and the right side function is(40)fx,y,t=1+2π2etsinπxsinπy.
The corresponding exact solution is(41)ux,y,t=etsinπxsinπy.
Leth = 1/19 (i.e., K = 3), h = 1/25 (i.e., K = 4), h = 1/31 (i.e., K = 5), and h = 1/37 (i.e., K = 6), and the maximum errors of Algorithm 2 for τ = 0.0004 at different times are displayed in Table 7. It is obvious that the accuracy of Algorithm 2 is satisfied. Table 8shows the comparison of the CPU calculation time among Algorithm 2, the classical C-N scheme, and the ASC-N scheme. Algorithm 2 not only has high accuracy but also has good parallel efficiency.Table 7
The maximum errors of Algorithm2 to Example 2 for τ = 0.0004.
h = 1/19
h = 1/25
h = 1/31
h = 1/37
t = 0.2
6.0275e − 4
1.5451e − 4
8.0677e − 5
4.9096e − 5
t = 0.4
4.4976e − 4
1.9272e − 4
1.0052e − 4
6.1127e − 5
t = 0.5
5.1593e − 4
2.1305e − 4
1.1113e − 4
6.7574e − 5
t = 0.6
5.7284e − 4
2.3546e − 4
1.2282e − 4
7.4684e − 5
t = 0.8
7.0014e − 4
2.8760e − 4
1.5001e − 4
9.1219e − 5Table 8
Comparison of three schemes’ calculation time forr = 1.2 at t = 0.5.
h
Algorithm2
C-N scheme
ASC-N scheme
1/37 (K = 6)
6.4121s
20.2625s
16.7089s
1/49 (K = 8)
21.4226s
49.9863s
41.0513s
1/61 (K = 10)
53.2623s
175.2632s
146.0665s
1/91 (K = 15)
341.5628s
1167.2047s
966.9273s


## 5. Conclusion

We have proposed and analyzed the SC-N scheme for parabolic equations. This algorithm consists of two DDMs, which split the Crank–Nicolson scheme; each one is used to solve the values at the same time level. Then, the average of two values is calculated. The SC-N scheme maintains the properties of the C-N scheme, i.e., unconditional stability and second-order numerical accuracy. Then, we extend the new algorithm to two-dimensional parabolic equations by ADI technique, which means that high-dimensional parabolic equations can be solved by the proposed algorithm in this paper. Numerical experiments illustrate the good performance of the new algorithm.



---

*Source: 8571625-2020-03-30.xml*