- Email: [email protected]
Domain decomposition schemes with high-order accuracy and unconditional stability
Applied Mathematics and Computation 219 (2013) 6170–6181
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Domain decomposition schemes with high-order accuracy and unconditional stability Wenrui Hao a, Shaohong Zhu b,⇑ a b
Department of Applied and Computational Mathematics and Statistics, University of Notre Dame, Notre Dame, IN 46556, USA School of Mathematics Science and LPMC, Nankai University, Tianjin 300071, China
a r t i c l e
i n f o
a b s t r a c t Parallel finite difference schemes with high-order accuracy and unconditional stability for solving parabolic equations are presented. The schemes are based on domain decomposition method, i.e., interface values between subdomains are computed by the explicit scheme; interior values are computed by the implicit scheme. The numerical stability and error are derived in the H1 norm in one dimensional case. Numerical results of both one and two dimensions examining the stability, accuracy, and parallelism of the procedure are also presented. Crown Copyright Ó 2012 Published by Elsevier Inc. All rights reserved.
Keywords: Domain decomposition Finite difference Parabolic equation High-order accuracy Unconditional stability
1. Introduction Domain decomposition is a powerful tool for devising parallel methods to solve time-dependent partial differential equations. There is rich literature on domain decomposition finite difference methods for solving parabolic equations on parallel computers. For the non-overlapping domain decomposition methods, the explicit nature of the calculation at the interface of sub-domain leads some domain decomposition schemes to be conditionally stable, which implies that they have to suffer from temporal step-size restrictions (see [1–5]). Schemes with unconditional stability as well as high-order accuracy being desired in the applications, many investigators have turned to improve the stability of the domain decomposition method. For example, the corrected explicit–implicit domain decomposition algorithms were presented in [6,7]. By adding the correction step to explicit–implicit domain decomposition methods, updating the interface solutions at each time level, the corrected methods were proved to have the unconditional stability. Meanwhile, the domain decomposition schemes without any correction steps were presented in [8,9], and proved to be unconditionally stable. All of these methods with unconditional stability reach the second order accuracy at most. The purpose of this paper is to present the domain decomposition finite difference procedure with third-order accuracy and unconditional stability. We first consider the following Dirichlet boundary problem
@U @ 2 U 2 ¼ 0; x 2 ð0; 1Þ; t 2 ð0; T; @t @x Uð0; tÞ ¼ Uð1; tÞ ¼ 0; t 2 ð0; T; Uðx; 0Þ ¼ U 0 ðxÞ;
0 6 x 6 1;
ð1:1Þ ð1:2Þ ð1:3Þ
where the initial function U 0 ðxÞ satisfies the boundary condition, i.e., U 0 ð0Þ ¼ U 0 ð1Þ ¼ 0. Then we extend the method to the problem of two dimensional space. ⇑ Corresponding author. E-mail addresses: [email protected] (W. Hao), [email protected] (S. Zhu). 0096-3003/$ - see front matter Crown Copyright Ó 2012 Published by Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2012.12.001
6171
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181
We will introduce two new finite difference schemes for solving (1.1)–(1.3) in Section 2, and sketch the domain decomposition procedure, the numerical stability and convergence in Section 3. In Section 5, the proof of the unconditional stability and the error estimate will be given. In Section 4, we extend the method to the problem of two dimensional space. Numerical examples and examination of the algorithm will be provided in Section 6. 2. Two new finite difference schemes Taking the usual h; s mesh in x and t and denoting the approximate value of Uðxj ; tn Þ U nj by unj , where xj ¼ jh and t n ¼ ns, we define the following operators
Dþ unj ¼ unjþ1 unj ;
D unj ¼ unj unj1 ;
1 Ds unj ¼ ðunj un1 Þ: j
s
It is well known that the following Taylor expansion resulting in the fully implicit finite difference scheme is valid.
Ds U nþ1 j
Dþ D U nþ1 j
2
h
! 2 h @2U 4 ¼ ðxj ; t nþ1 Þ þ Oðs2 þ h Þ: þ 2 12 @t 2
s
ð2:4Þ
Noticing that nþ1
Uj @2U ðxj ; tnþ1 Þ ¼ @t 2 and substituting
Ds U nþ1 j
@2 U @t 2
s2
þ OðsÞ
ð2:5Þ
ðxj ; tnþ1 Þ into (2.4), we obtain
Dþ D U nþ1 j
2U nj þ U n1 j
2
h
! nþ1 n n1 2 h U j 2U j þ U j 4 2 þ ¼ Oðs2 þ h þ sh Þ: þ 2 12 s2
s
ð2:6Þ
Replacing j by j þ 1 and j 1 in (2.5), we have n n1 2 U nþ1 jþ1 ¼ 2U jþ1 U jþ1 þ s
@2U ðxjþ1 ; t nþ1 Þ þ Oðs3 Þ @t 2
n n1 2 U nþ1 j1 ¼ 2U j1 U j1 þ s
@2U ðxj1 ; t nþ1 Þ þ Oðs3 Þ: @t 2
and
Substituting U nþ1 j1 into (2.6), we can get
Ds U nþ1 j ¼
nþ1 2U njþ1 U n1 þ 2U nj1 U n1 jþ1 2U j j1
h
s2 @ 2 U h
2
@t 2
ðxjþ1 ; tnþ1 Þ þ
2
! nþ1 n n1 2 h U j 2U j þ U j þ þ 2 12 s2
s
! @2U s3 4 2 nþ1 ðx ; t Þ þ Oðs2 þ h þ sh þ 2 Þ: j1 2 @t h
ð2:7Þ
By omitting the high order term, (2.6) and (2.7) yield two new finite difference schemes for (1.1):
! nþ1 2 n n1 h uj 2uj þ uj ¼ 0; 2 2 2 12 s h ! Dþ D unþ1 2unj þ un1 s h2 unþ1 j j j n nþ1 Ds unþ1 þ rðDs unjþ1 Ds unþ1 þ j jþ1 Þ rðDs uj1 Ds uj1 Þ ¼ 0; 2 2 12 s2 h
Ds unþ1 j
Dþ D unþ1 j
þ
s
ð2:8Þ
þ
ð2:9Þ 4
2
wherer ¼ hs2 . From the derivation of (2.8) and (2.9), we know that the truncation errors of (2.8) and (2.9) are Oðs2 þ h þ sh Þ 2 3 4 2 4 2 and O sh2 þ s2 þ h þ sh þ hs2 respectively. If r is any positive real number, the truncation errors become Oðh Þ and Oðh Þ. 3. Domain decomposition procedure and main results Suppose Jh ¼ 1; N s ¼ T. For simplicity, we will consider a domain decomposition which involves in decomposing ð0; 1Þ into only two subdomains, ð0; xÞ and ð x; 1Þ, where x ¼ xk for some integer k ð1 < k < J 1Þ. We use the explicit scheme (2.9) to compute the solution value unþ1 and the implicit scheme (2.8) to compute other solution values unþ1 ðj – kÞ respectively. j k The system can be written as
Lðunþ1 Þ ¼ 0; j unþ1 j
¼ 0;
1 6 n 6 N 1;
j ¼ 0; J;
0 < j < J;
ð3:10Þ
6172
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181
where the linear operator L is
Lðunþ1 Þ¼ j
8 unþ1 2un þun1 Dþ D unþ1 j j j j > nþ1 s h2 > ; > Ds uj h2 þ 2 þ 12 s2 > > < Dþ D unþ1
2
j – k;
unþ1 2un þun1
j j j h Ds unþ1 h2 j þ 2s þ 12 > j s2 > > > > : rðD un D unþ1 Þ rðD un D unþ1 Þ; s jþ1 s jþ1 s j1 s j1
ð3:11Þ
: j ¼ k:
The resulting system of equations decouples into two disjoint sets of equations corresponding to the subdomains. These systems can be solved in parallel. We are now in position to state two main theorems of this paper, which will be proved in the forthcoming sections. For the discrete function un ¼ funj jj ¼ 0; 1; . . . ; J; un0 ¼ unJ ¼ 0g, define
kun k2 ¼
J1 X junj j2 h;
n 2 kDþ un k2 ¼ sumJ1 j¼0 jDþ uj j h:
j¼1
We have the following theorems. Theorem 3.1. (Stability). For any given r > 0, the finite difference solutions of the parallel scheme (3.10)–(3.11) satisfy
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 þ maxkDþ u k 6 þ 1ðkDþ u1 k þ ku1 u0 kÞ: n 2r 12r2 n
Theorem 3.2. (Convergence). Let enj ¼ Uðxj ; tn Þ unj . For any given r > 0, the finite difference solutions of the parallel scheme (3.10)–(3.11) satisfy 4
maxkDþ en k 6 CðkDþ e1 k þ ke1 e0 k þ h Þ; n
where C is a positive constant independent of h and s. Since the finite difference scheme (3.10) has three time levels, besides taking u0j ¼ U 0 ðjhÞ, thus e0j ¼ 0, we need to find other methods to solve u1j ; 8j. In order to match fourth order accuracy, we can use either fourth order explicit schemes, such 4 as the impact scheme, or the high-order parallel iterative method [10] to compute u1j , i.e., let je1j j ¼ Oðh Þ. Then according to Theorem 3.2, scheme (3.10) will reach third-order accuracy. 4. Extension to two dimensional case In this section, Uðx; y; tÞ will be a solution of the following Dirichlet boundary problem on X ¼ ð0; 1Þ ð0; 1Þ,
@U DU ¼ 0; @t
ðx; yÞ 2 X;
t 2 ð0; T;
ð4:12Þ
Uðx; y; tÞ ¼ 0;
ðx; yÞ 2 @ X;
t 2 ð0; T;
ð4:13Þ
Uðx; y; 0Þ ¼ U 0 ðx; yÞ;
ðx; yÞ 2 X;
ð4:14Þ 2
2
where the initial function U 0 ðx; yÞ satisfies the boundary condition, and DU ¼ @@xU2 þ @@yU2 . We extend the domain decomposition method stated in Section 3 to the above problem. Take
X1 ¼ fðx; yÞ 2 X : x < xg; X2 ¼ fðx; yÞ 2 X : x > xg: Let xi ¼ ih same as in Section 3, and let yj ¼ jh. Suppose that there exists an integer k such that x ¼ xk (see Fig. 1). In analogy with Section 3 we call points ðxi ; yj ; tn Þ as boundary points if ðxi ; yj Þ 2 @ X and interface points if i ¼ k. Otherwise, we call them interior points. The values uni;j will approximate Uðxi ; yj ; t n Þ U ni;j . We denote
Dþ D ui;j ¼ uiþ1;j 2ui;j þ ui1;j ; and
Ds uni;j ¼
uni;j un1 i;j
s
:
dþ d ui;j ¼ ui;jþ1 2ui;j þ ui;j1
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181
6173
Fig. 1. Domain decomposition for two dimensional case.
Then
Dþ D U i;j 2
h dþ d U i;j
2
¼
@2U h @4U 4 ðxi ; yj ; tÞ þ ðxi ; yj ; tÞ þ Oðh Þ; 2 @x 12 @x4 2
¼
2
@2U h @4U 4 ðxi ; yj ; tÞ þ ðxi ; yj ; tÞ þ Oðh Þ; 2 @y 12 @y4
h dþ d Dþ D U i;j 4
h
¼
@4U 2 ðxi ; yj ; tÞ þ Oðh Þ: @x2 @y2
Thus we have
Dþ D U i;j þ dþ d U i;j
1 d d D D U 2 þ þ i;j 6h ! 2 h @4U @4U @4U 4 ¼ DUðxi ; yj ; tÞ þ ðxi ; yj ; tÞ þ 2 2 2 ðxi ; yj ; tÞ þ 4 ðxi ; yj ; tÞ þ Oðh Þ @x @y @y 12 @x4 2
h
þ
2
¼ DUðxi ; yj ; tÞ þ
h @2U 4 ðxi ; yj ; tÞ þ Oðh Þ: 12 @t2
ð4:15Þ
Moreover,
Ds U nþ1 ¼ i;j
@U s @2U ðxi ; yj ; t nþ1 Þ ðxi ; yj ; t nþ1 Þ þ Oðs2 Þ: @t 2 @t 2
ð4:16Þ
Subtract (4.16) with (4.15), we obtain
Ds U nþ1 i;j
Dþ D U nþ1 þ dþ d U nþ1 i;j i;j h
2
1
d d D D U nþ1 2 þ þ i;j
6h
! 2 h @2U 4 ¼ ðxi ; yj ; tnþ1 Þ þ Oðs2 þ h Þ: þ 2 12 @t2
s
Therefore, the fourth order finite difference scheme for (4.12) is
Ds unþ1 i;j
Dþ D unþ1 þ dþ d unþ1 i;j i;j h
2
1 6h
d d D D unþ1 2 þ þ i;j
! nþ1 2 n n1 h ui;j 2ui;j þ ui;j þ ¼ 0; þ 2 12 s2
s
ð4:17Þ
which is a nine-point finite difference scheme with three levels (see Fig. 2). Denote
ui;j ¼ ui1;j1 þ ui1;jþ1 þ uiþ1;j1 þ uiþ1;jþ1 4ui;j D 2 1 ¼ Dþ D ui;j þ dþ d ui;j þ dþ d Dþ D ui;j : 2 Then (4.17) is equivalent to
Ds unþ1 i;j
unþ1 þ 2Dþ D unþ1 þ 2dþ d unþ1 D i;j i;j i;j 3h
2
! nþ1 2 n n1 h ui;j 2ui;j þ ui;j þ ¼ 0; þ 2 12 s2
s
n n1 which can be used for computing unþ1 when i – k. When i ¼ k; unþ1 i;j i1;j in (4.17) are replaced by 2ui1;j ui1;j . In this case, we nþ1 need to solve a tridiagonal matrix using Thomas algorithm for computing uk;j .
6174
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181
Fig. 2. Nine point scheme.
It is straightforward to extend the two-subdomain results to many subdomains on x direction by cutting the whole domain into vertical strips. Moreover, we can also extend the domain decomposition scheme to three dimensional case by dividing x axis into many subdomains. 5. Proof of stability and convergence We first state three auxiliary lemmas. The stability and convergence results are then derived. Lemma 5.1 (Discrete Poincare Inequality). For the discrete function un ¼ funj jj ¼ 0; 1; . . . ; J; un0 ¼ unJ ¼ 0g, there exists
kun k 6 k
Dþ un k: h
Lemma 5.2 (Discrete Green Theorem). If uj and
v j are discrete functions on the set fxj jj ¼ 0; 1; . . . ; Jg, then we have
J1 J X X uj Dþ v j ¼ v j D uj u0 v 1 þ uJ v J : j¼1
j¼1
Lemma 5.3. For any given r > 0; f 2 L1 ð0; T; L2 ð0; 1ÞÞ, the finite difference solutions of system Lðunþ1 Þ ¼ fjnþ1 with the Dirichlet j boundary condition satisfy
maxkDþ un k2 6 kDþ u1 k2 þ n
1 1 2 þ þ 1 ku1 u0 k2 þ Th maxkf n k2 : 2 n 2r 12r
Lemma 5.1 and 5.2 are proved in [11]. Before proving the stability and convergence, we will give the proof of Lemma 5.3. 5.1. Proof of Lemma 5.3 Denoting wnþ1 ¼ unþ1 unj , we can rewrite Lðunþ1 Þ ¼ fjnþ1 with the Dirichlet boundary condition as j j j
8 nþ1 nþ1 1 w r Dþ D unþ1 þ 12 þ 12r ðwj wnj Þ ¼ fjnþ1 s; > j > > j > nþ1 > nþ1 nþ1 1 1 rðwnjþ1 wnþ1 s; > jþ1 Þ rðwj1 wj1 Þ ¼ fj > > > : unþ1 ¼ 0; j ¼ 0; J:
j – k;
j ¼ k;
j
Multiplying the above equations by wnþ1 h; j ¼ 1; 2; . . . ; J 1 and summing them up with respect to j, we have j
J1 J1 J1 X X 1 1 X nþ1 n nþ1 nþ1 n nþ1 þ ðwnþ1 Þ2 h r wnþ1 Dþ D unþ1 hþ ðw wnj Þwnþ1 h rwnþ1 j j j j k ðwkþ1 wkþ1 Þh rwk ðwk1 wk1 Þh 2 12r j¼1 j j¼1 j¼1 ¼
J1 X fjnþ1 wnþ1 sh: j j¼1
ð5:18Þ
6175
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181
From Lemma 5.2, we have J1 J J1 J1 J1 X X 1X 1X 1X wnþ1 Dþ D unþ1 ¼ D unþ1 D wnþ1 ¼ ðDþ unþ1 Þ2 þ ðDþ unj Þ2 ðDþ unþ1 Dþ unj Þ2 : j j j j j j 2 2 2 j¼1 j¼1 j¼0 j¼0 j¼0
Then (5.18) becomes
i r XJ1 XJ1 r hXJ1 1 1 nþ1 2 n 2 nþ1 n 2 þ kwnþ1 k2 þ ð D u Þ h ð D u Þ h þ ð D u D u Þ h þ þ þ þ þ j j j j j¼0 j¼0 j¼0 2 2 2 12r X
X X J1 J1 J1 1 1 1 nþ1 nþ1 ðwnþ1 Þ2 h ðwnj Þ2 h þ ðwnþ1 wnj Þ2 h þ rhwk ðwnþ1 j j kþ1 þ wk1 Þ j¼1 j¼1 j¼1 2 2 2 J1 X nþ1 ¼ rhwk ðwnkþ1 þ wnk1 Þ þ wnþ1 fjnþ1 sh j j¼1 2 2 n ðwnþ1 ðwnþ1 Þ2 þ ðwnk1 Þ2 k Þ þ ðwkþ1 Þ 6 rh þ k 2 2 2 6 rhðwnþ1 k Þ þ rh
! þ
J1 X ðwnþ1 Þ2 þ ðfjnþ1 sÞ2 j
2
j¼1
h
ðwnkþ1 Þ2 þ ðwnk1 Þ2 1 nþ1 2 s2 nþ1 2 þ kw k þ kf k : 2 2 2
ð5:19Þ
Noting that
" # 2 nþ1 2 ðwnþ1 r XJ1 r XJ1 nþ1 2 kþ1 Þ þ ðwk1 Þ nþ1 n 2 nþ1 nþ1 nþ1 2 þ rhðwnþ1 ð D D ð D Þ h þ rh þ uj þ uj Þ h þ rhwk ðwkþ1 þ wk1 Þ ¼ þ wj k Þ ; j¼0 j¼0;j–k;k1 2 2 2 we can simplify (5.19) as follows
i r XJ1 1 nþ1 2 r h 1 1 nþ1 2 kw k þ kDþ unþ1 k2 kDþ un k2 þ þ ðkwnþ1 k2 kwn k2 þ kwnþ1 wn k2 Þ ð D w Þ h þ þ j j¼0;j–k;k1 2 2 2 4 24r " # 2 nþ1 2 ðwnþ1 ðwnkþ1 Þ2 þ ðwnk1 Þ2 s2 kþ1 Þ þ ðwk1 Þ þ rh 6 kf nþ1 k2 : 2 2 2 Since kwnþ1 k2 P 0;
PJ1
nþ1 2 Þ h j¼0;j–k;k1 ðDþ wj
P 0 and kwnþ1 wn k2 P 0, then we have
" # 2 nþ1 2 i 1 ðwnþ1 ðwnkþ1 Þ2 þ ðwnk1 Þ2 rh 1 s2 kþ1 Þ þ ðwk1 Þ nþ1 2 n 2 nþ1 2 n 2 6 kf nþ1 k2 : kDþ u k kDþ u k þ þ ðkw k kw k Þ þ rh 2 4 24r 2 2 2 Summing up with respect to n, we get
2 ðwnþ1 Þ2 þ ðwnþ1 ðw1 Þ2 þ ðw1k1 Þ2 r 1 1 r 1 1 k1 Þ kDþ unþ1 k2 þ þ kwnþ1 k2 þ rh kþ1 þ kw1 k2 þ rh kþ1 6 kDþ u1 k2 þ 2 4 24r 2 4 24r 2 2 Ts maxkf n k2 : þ 2 n Thus
kDþ unþ1 k2 6 kDþ u1 k2 þ Lemma 5.3 is proved.
1 1 2 þ 1 ku1 u0 k2 þ Th maxkf n k2 : þ n 2r 12r 2
h
5.2. Proof of stability Let f ¼ 0 in Lemma 5.3, thus Theorem 3.1 is proved.
h
5.3. Proof of convergence According to the parallel scheme (3.10)–(3.11), the errors enþ1 ð1 6 n < NÞ follows the below equations: j
! nþ1 2 n n1 h ej 2ej þ ej þ ¼ Gnþ1 ; þ j 2 2 2 12 s h ! Dþ D enþ1 2enj þ en1 s h2 enþ1 j j j Ds enþ1 þ þ j 2 2 12 s2 h
Ds enþ1 j
Dþ D enþ1 j
s
2
nþ1 n nþ1 rðDs enjþ1 Ds enþ1 h þ Gnþ1 ; j jþ1 Þ rðDs ej1 Ds ej1 Þ ¼ U
j – k;
j ¼ k;
enþ1 ¼ 0; j
j ¼ 0; J;
6176
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181 4
2
where Unþ1 ¼ 2r 2 @@tU2 ðxk ; tnþ1 Þ; jGnþ1 j ¼ Oðh Þ. In order to get the error estimates for enþ1 , we assume that enþ1 ¼ pnþ1 þ qnþ1 , j j j j j nþ1 where pnþ1 and q are the solutions of the following problems respectively. j j problem I
! nþ1 2 n n1 h pj 2pj þ pj ¼ Gnþ1 ; j – k; j 2 2 2 12 s h ! Dþ D pnþ1 2pnj þ pn1 s h2 pnþ1 j j j nþ1 n nþ1 Ds pnþ1 þ rðDs pnjþ1 Ds pnþ1 ; þ j jþ1 Þ rðDs pj1 Ds pj1 Þ ¼ Gj 2 2 12 s2 h
Ds pnþ1 j
Dþ D pnþ1 j
pnþ1 ¼ 0; j
s
þ
þ
j ¼ k;
ð5:20Þ
j ¼ 0; J;
p0j ¼ e0j ;
p1j ¼ e1j ; 8j:
problem II
Ds qnþ1 j Ds qnþ1 j
Dþ D qnþ1 j 2
h
Dþ D qnþ1 j 2
h
s
þ
2
s
þ
2
þ þ
! nþ1 2 n n1 h qj 2qj þ qj ¼ 0; 2 12 s ! nþ1 2 n n1 h qj 2qj þ qj 12
s2 2
nþ1 n nþ1 rðDs qnjþ1 Ds qnþ1 h ; jþ1 Þ rðDs qj1 Ds qj1 Þ ¼ U
qnþ1 j
¼ 0;
j – k;
j ¼ k;
j ¼ 0; J;
q0j ¼ q1j ¼ 0; 8j: From Lemma 5.3, we can obtain the estimate for pnþ1 , i.e., j
1 1 2 þ kDþ pn k2 6 kDþ p1 k2 þ þ 1 kp1 p0 k2 þ Th maxkGn k2 2 n 2r 12r 1 1 10 6 þ 1 ðkDþ e1 k þ ke1 e0 kÞ2 þ Oðh Þ: þ 2r 12r 2
ð5:21Þ
nþ1 For estimating qnþ1 , we first consider q ð1 6 n 6 N 1Þ which satisfies the following equations j j
nþ1 Dþ D q j 2
h nþ1 Dþ D q k
nþ1 q 0
¼ 0;
j – k; 2
¼ Unþ1 h ; 2 h nþ1 ¼q ¼ 0: J
nþ1 Then the formula of q is j
nþ1 q j
8 0; j ¼ 0; J; > > < j PJ nþ1 4 ¼ J kþ1 U h ; 1 6 j 6 k; > P > : Jj k1 Unþ1 h4 ; k < j < J: i¼1 J
Hence
nþ1 k ¼ Oðh4 Þ and kq nþ1 q n k ¼ Oðsh3 Þ: kDþ q ~nþ1 nþ1 Define q ¼ qnþ1 q , then we have j j j
! nþ1 2 ~ j 2q ~nj þ q ~n1 h q j þ ¼ Rnþ1 ; j – k; þ j 2 2 2 12 s h ! ~nþ1 ~nj þ q ~n1 Dþ D q 2q s h2 q~nþ1 j j j nþ1 ~nþ1 ~njþ1 Ds q ~nþ1 ~n ~nþ1 Ds q þ rðDs q ; þ j jþ1 Þ rðDs qj1 Ds qj1 Þ ¼ Rj 2 2 2 12 s h ~nþ1 Ds q j
where
~nþ1 Dþ D q j
s
~nþ1 ¼ 0; q j
j ¼ 0; J;
~0j ¼ q 0j ; q
~1j ¼ q 1j ; 8j; q
j ¼ k;
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181
Rnþ1 j
6177
8 qnþ1 2qn þqn1 j j j h2 > nþ1 > Ds q þ 2s þ 12 ; j – k; > j s2 > < nþ1 2q n þq n1 q 2 ¼ D q j j j nþ1 h njþ1 Ds q nþ1 þ 2s þ 12 rðDs q s j > jþ1 Þ s2 > > > : nþ1 n j1 Ds q j1 Þ; j ¼ k; rðDs q 3
~nþ1 k. and kRnþ1 k ¼ Oðh Þ. From Lemma 5.3, we can obtain the estimate of kDþ q
~nþ1 k2 6 kDþ q ~1 k2 þ kDþ q
1 1 ~1 q ~0 k2 þ Th2 maxkRn k2 ¼ Oðh8 Þ: þ 1 kq þ n 2r 12r 2
Thus, 4
~nþ1 k þ kDþ q nþ1 k ¼ Oðh Þ: kDþ qnþ1 k 6 kDþ q Combining with (5.21), we get 4
kDþ enþ1 k 6 kDþ pnþ1 k þ kDþ qnþ1 k 6 CðkDþ e1 k þ ke1 e0 k þ h Þ; where C is a positive constant independent of h and s. Theorem 3.2 is proved.
h
6. Numerical experiments In this section, some numerical results are presented to show the stability, accuracy, and parallelism of the scheme described above, and the computational costs are also presented. All the experiments are run on a cluster consisting of a manager that uses one core of a Xeon 5410 processor and up to 32 computing nodes, each containing two Xeon 5410 processors running 64-bit Linux, i.e., each node consists of 8 processing cores. 6.1. One dimensional numerical experiments We consider the problem defined in equations (1.1)–(1.3) with U 0 ðxÞ ¼ sinðpxÞ. Obviously the exact solution of the equa2 tions is Uðx; tÞ ¼ ep t sinðpxÞ. First, we verify the stability of the scheme by taking the step size h ¼ 103 ; r ¼ 1; 10; 100; 1000 with two subdomains. Fig. 3 clearly shows that the norm of un doesn’t occur blowing up even if r is large enough. This explains the unconditional stability of the scheme. Second, we examine the numerical errors in the solutions. Table 1 shows that the errors for each case are roughly of the 3 same order of magnitude, and the errors appear to be Oðh Þ in each case. We also compare them with those resulting from
Fig. 3. The infinity norm of un v.s. t.
6178
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181
Table 1 Numerical errors for different grid length. 2 Processors
4 Processors
J
ku U k=h
r ¼ 10; T ¼ 1 1000 2000 4000 8000
J
2.06821 2.06715 2.06316 1.94919
r ¼ 100; T ¼ 1 1000 2000 4000 8000
8 Processors ku U k=h
J
kun U n k=h
1000 2000 4000 8000
4.05448 4.05336 4.04938 3.93590
1000 2000 4000 8000
8.07385 8.07276 8.06881 7.95632
2.06703e2 2.06729e2 2.06690e2 2.06820e2
1000 2000 4000 8000
4.04766e2 4.05299e2 4.05307e2 4.05434e2
1000 2000 4000 8000
8.04137e2 8.06955e2 8.07216e2 8.07356e2
r ¼ 1000; T ¼ 1 1000 2000 4000 8000
2.01883e4 2.03740e4 2.06360e4 2.06634e4
1000 2000 4000 8000
4.04424e4 3.93818e4 4.03957e4 4.05137e4
1000 2000 4000 8000
8.03459e4 8.03843e4 8.01579e4 8.06558e4
r ¼ 10000; T ¼ 1 1000 2000 4000 8000
2.02383e6 2.04768e6 2.5982e6 2.06920e6
1000 2000 4000 8000
4.05348e6 4.10234e6 4.00321e6 4.01203e6
1000 2000 4000 8000
8.00374e6 8.04875e6 8.00829e6 8.06558e6
n
n
3
n
n
3
3
Table 2 Numerical errors for the scheme in [8]. 2 Processors
4 Processors
J
kun U n k=h
r ¼ 10; T ¼ 1 1000 2000 4000 8000
8 Processors
J
kun U n k=h
J
kun U n k=h
1.02839 1.00293 1.12303 1.00290
1000 2000 4000 8000
2.03842 2.20394 2.04095 2.03938
1000 2000 4000 8000
4.09382 4.02910 4.10293 3.98201
r ¼ 100; T ¼ 1 1000 2000 4000 8000
1.09282e1 1.02938e1 1.10293e1 1.00927e1
1000 2000 4000 8000
2.09384e1 2.07261e1 2.01628e1 2.10928e1
1000 2000 4000 8000
4.08271e1 4.08726e1 4.09281e1 4.09817e1
r ¼ 1000; T ¼ 1 1000 2000 4000 8000
1.03928e2 1.09282e2 1.07621e2 1.01822e2
1000 2000 4000 8000
2.00971e2 2.00761e2 2.00127e2 2.02872e2
1000 2000 4000 8000
4.02893e2 4.02934e2 4.00029e2 4.07625e2
r ¼ 10000; T ¼ 1 1000 2000 4000 8000
1.00121e3 1.00273e3 1.02922e3 1.00461e3
1000 2000 4000 8000
2.04759e3 2.07162e3 2.00937e3 2.00183e3
1000 2000 4000 8000
4.01123e3 4.08211e3 4.00182e3 4.09837e3
2
2
Table 3 Comparison of scheme. Processors n
1 2 4 8 16 32 64
The third order scheme
The second order scheme in [8]
Time tn (s)
Speed-up t1 =tn
Time tn (s)
Speed-up t 1 =t n
4862.66 2450.82 1250.29 621.74 325.05 165.85 100.45
– 1.98 3.89 7.82 14.96 29.32 48.40
4564.28 2340.66 1164.36 598.98 300.84 157.82 98.97
– 1.95 3.92 7.62 15.21 28.92 46.12
2
6179
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181 Table 4 Comparison of time for two dimensional case. Processors n
Time tn (s)
Speed-up t1 =tn
1 2 4 8 16 32 64
10 h 23 m 19 s 6 h 01 m 29 s 2 h 52 m 14 s 1 h 23 m 54 s 42 m 47 s 21 m 23 s 12 m 58 s
– 1.72 3.61 7.41 14.53 29.06 47.91
Table 5 Numerical errors for two dimensional cases. 8 Processors
16 Processors
Jx Jy
ku U k=h
r ¼ 10; T ¼ 1 100 100 200 200 400 400
10.07867 10.06592 10.02345
r ¼ 100; T ¼ 1 100 100 200 200 400 400
n
n
1.09837e4 1.05629e4 1.07820e4
3
32 Processors
Jx Jy
ku U k=h
100 100 200 200 400 400
40.09821 40.06710 40.01023
100 100 200 200 400 400
n
n
3
4.09272e4 4.09983e4 4.05730e4
Jx Jy
kun U n k=h
100 100 200 200 400 400
160.09802 160.02893 160.07802
100 100 200 200 400 400
3
1.60029e5 1.60789e5 1.60892e5
Fig. 4. Initial condition for the Lotka–Volterra system.
the second order scheme in [8], which are shown in Table 2. Table 1 and Table 2 show that our domain decomposition algorithm has higher order accuracy. Third, we test the speed-up for the scheme. Here we take h ¼ 105 ; s ¼ 106 , i.e., r ¼ 104 ; T ¼ 0:5 and list the time for computing and the speed-up in Table 3 for our domain decomposition algorithm and the scheme in [8]. Table 3 shows that our domain decomposition algorithm is competitive with the existing method in computing time.
6180
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181 Table 6 Comparison of time for Lotka–Volterra system. Processors n
Time tn (s)
Speed-up t 1 =tn
1 2 4 8 16 32 64 128 256
21 h 51 m 32 s 11 h 45 m 8 s 6 h 12 m 36 s 2 h 58 m 57 s 1 h 33 m 37 s 45 m 28 s 28 m 12 s 14 m 33 s 7 m 28 s
– 1.86 3.52 7.37 14.16 28.85 46.51 90.12 175.61
Fig. 5. Numerical solutions for T ¼ 1.
6.2. Two dimensional numerical experiments We consider the problem defined in equations (4.12)–(4.14) with U 0 ðx; yÞ ¼ sinðpxÞ sinðpyÞ. Obviously the exact solution 2 of the equations is Uðx; y; tÞ ¼ e2p t sinðpxÞ sinðpyÞ. The computing time and the speed-up of the scheme are shown in Table 4 3 2 with h ¼ 10 ; s ¼ 10 and T ¼ 1. Numerical errors are shown in Table 5 and also achieve the third order accuracy.
6.3. Application to Lotka–Volterra system The competitive Lotka–Volterra equations are a simple model of the population dynamics of species competing for some common resource. Given two populations, u and v, with logistic dynamics, the Lotka–Volterra formulation adds an additional term to account for the species’ interactions. Thus the competitive Lotka–Volterra equations are:
ut ¼ Du þ uð1 v Þ;
ðx; yÞ 2 X
v t ¼ Dv v ð1 uÞ; ðx; yÞ 2 X v ¼ u ¼ 0; ðx; yÞ 2 @ X: We choose random values for uðx; yÞ and v ðx; yÞ as initial conditions shown in Fig. 4. We take h ¼ 0:001, s ¼ 0:01; T ¼ 1 and show the computing time in Table 6. Fig. 5 shows the numerical solution for T ¼ 1.
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181
6181
7. Conclusion We presented domain decomposition finite difference schemes with unconditional stability and third-order accuracy for the parabolic system. Error estimate and stability of the numerical solutions have been derived for one dimensional case. The scheme is easy to implement the parallelism and is extended in two and three dimensional case. The numerical results demonstrate the good performance of the parallel scheme, namely, unconditional stability, the third order accuracy and high degree of parallelism. References [1] C.N. Dawson, Q. Du, T.F. Dupont, A finite difference domain decomposition algorithm for numerical solution of the Heat equation, Math. Comp. 57 (1991) 63–71. [2] G. Yuan, S. Zhu, L. Shen, Domain decomposition algorithm based on the group explicit formula for the heat equation, Int. J. Comput. Math 82 (2005) 1295–1306. [3] S. Zhu, Z. Yu, J. Zhao, A high-order parallel finite difference algorithm, Appl. Math. Comp. 183 (2006) 365–372. [4] C.N. Dawson, T.F. Dupont, Explicit/implicit conservative domain decomposition procedures for parabolic problems based on block-centered finite differences, SIAM J. Numer. Anal. 31 (1994) 1045–1061. [5] G. Yuan, L. Shen, Stability and convergence of the explicit–implicit conservative domain decomposition procedure for parabolic problems, Comput. Math. Appl. 47 (2004) 793–801. [6] H. Shi, H. Liao, Unconditional stability of corrected explicit–implicit domain decomposition algorithms for parallel approximation of heat equations, SIAM J. Numer. Anal. 44 (2006) 1584–1611. [7] H. Liao, H. Shi, Z. Sun, Corrected explicit–implicit domain decomposition algorithms for two-dimensional semilinear parabolic equations, Sci. China Series A 52 (2009) 2362–2388. [8] Z. Sheng, G. Yuan, X. Hang, Unconditional stability of parallel difference schemes with second order accuracy for parabolic equation, Appl. Math. Compt. 184 (2007) 1015–1031. [9] S. Zhu, Conservative domain decomposition procedure with unconditional stability and second-order accuracy, Appl. Math. Compt. 216 (2010) 3275– 3282. [10] W. Hao, S. Zhu, Parallel iterative methods for parabolic equations, Int. J. Comput. Math 86 (2009) 431–440. [11] Y. Zhou, Applications of discrete functional analysis to the finite difference method, International Academic publishers (1991).
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Domain decomposition schemes with high-order accuracy and unconditional stability Wenrui Hao a, Shaohong Zhu b,⇑ a b
Department of Applied and Computational Mathematics and Statistics, University of Notre Dame, Notre Dame, IN 46556, USA School of Mathematics Science and LPMC, Nankai University, Tianjin 300071, China
a r t i c l e
i n f o
a b s t r a c t Parallel finite difference schemes with high-order accuracy and unconditional stability for solving parabolic equations are presented. The schemes are based on domain decomposition method, i.e., interface values between subdomains are computed by the explicit scheme; interior values are computed by the implicit scheme. The numerical stability and error are derived in the H1 norm in one dimensional case. Numerical results of both one and two dimensions examining the stability, accuracy, and parallelism of the procedure are also presented. Crown Copyright Ó 2012 Published by Elsevier Inc. All rights reserved.
Keywords: Domain decomposition Finite difference Parabolic equation High-order accuracy Unconditional stability
1. Introduction Domain decomposition is a powerful tool for devising parallel methods to solve time-dependent partial differential equations. There is rich literature on domain decomposition finite difference methods for solving parabolic equations on parallel computers. For the non-overlapping domain decomposition methods, the explicit nature of the calculation at the interface of sub-domain leads some domain decomposition schemes to be conditionally stable, which implies that they have to suffer from temporal step-size restrictions (see [1–5]). Schemes with unconditional stability as well as high-order accuracy being desired in the applications, many investigators have turned to improve the stability of the domain decomposition method. For example, the corrected explicit–implicit domain decomposition algorithms were presented in [6,7]. By adding the correction step to explicit–implicit domain decomposition methods, updating the interface solutions at each time level, the corrected methods were proved to have the unconditional stability. Meanwhile, the domain decomposition schemes without any correction steps were presented in [8,9], and proved to be unconditionally stable. All of these methods with unconditional stability reach the second order accuracy at most. The purpose of this paper is to present the domain decomposition finite difference procedure with third-order accuracy and unconditional stability. We first consider the following Dirichlet boundary problem
@U @ 2 U 2 ¼ 0; x 2 ð0; 1Þ; t 2 ð0; T; @t @x Uð0; tÞ ¼ Uð1; tÞ ¼ 0; t 2 ð0; T; Uðx; 0Þ ¼ U 0 ðxÞ;
0 6 x 6 1;
ð1:1Þ ð1:2Þ ð1:3Þ
where the initial function U 0 ðxÞ satisfies the boundary condition, i.e., U 0 ð0Þ ¼ U 0 ð1Þ ¼ 0. Then we extend the method to the problem of two dimensional space. ⇑ Corresponding author. E-mail addresses: [email protected] (W. Hao), [email protected] (S. Zhu). 0096-3003/$ - see front matter Crown Copyright Ó 2012 Published by Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2012.12.001
6171
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181
We will introduce two new finite difference schemes for solving (1.1)–(1.3) in Section 2, and sketch the domain decomposition procedure, the numerical stability and convergence in Section 3. In Section 5, the proof of the unconditional stability and the error estimate will be given. In Section 4, we extend the method to the problem of two dimensional space. Numerical examples and examination of the algorithm will be provided in Section 6. 2. Two new finite difference schemes Taking the usual h; s mesh in x and t and denoting the approximate value of Uðxj ; tn Þ U nj by unj , where xj ¼ jh and t n ¼ ns, we define the following operators
Dþ unj ¼ unjþ1 unj ;
D unj ¼ unj unj1 ;
1 Ds unj ¼ ðunj un1 Þ: j
s
It is well known that the following Taylor expansion resulting in the fully implicit finite difference scheme is valid.
Ds U nþ1 j
Dþ D U nþ1 j
2
h
! 2 h @2U 4 ¼ ðxj ; t nþ1 Þ þ Oðs2 þ h Þ: þ 2 12 @t 2
s
ð2:4Þ
Noticing that nþ1
Uj @2U ðxj ; tnþ1 Þ ¼ @t 2 and substituting
Ds U nþ1 j
@2 U @t 2
s2
þ OðsÞ
ð2:5Þ
ðxj ; tnþ1 Þ into (2.4), we obtain
Dþ D U nþ1 j
2U nj þ U n1 j
2
h
! nþ1 n n1 2 h U j 2U j þ U j 4 2 þ ¼ Oðs2 þ h þ sh Þ: þ 2 12 s2
s
ð2:6Þ
Replacing j by j þ 1 and j 1 in (2.5), we have n n1 2 U nþ1 jþ1 ¼ 2U jþ1 U jþ1 þ s
@2U ðxjþ1 ; t nþ1 Þ þ Oðs3 Þ @t 2
n n1 2 U nþ1 j1 ¼ 2U j1 U j1 þ s
@2U ðxj1 ; t nþ1 Þ þ Oðs3 Þ: @t 2
and
Substituting U nþ1 j1 into (2.6), we can get
Ds U nþ1 j ¼
nþ1 2U njþ1 U n1 þ 2U nj1 U n1 jþ1 2U j j1
h
s2 @ 2 U h
2
@t 2
ðxjþ1 ; tnþ1 Þ þ
2
! nþ1 n n1 2 h U j 2U j þ U j þ þ 2 12 s2
s
! @2U s3 4 2 nþ1 ðx ; t Þ þ Oðs2 þ h þ sh þ 2 Þ: j1 2 @t h
ð2:7Þ
By omitting the high order term, (2.6) and (2.7) yield two new finite difference schemes for (1.1):
! nþ1 2 n n1 h uj 2uj þ uj ¼ 0; 2 2 2 12 s h ! Dþ D unþ1 2unj þ un1 s h2 unþ1 j j j n nþ1 Ds unþ1 þ rðDs unjþ1 Ds unþ1 þ j jþ1 Þ rðDs uj1 Ds uj1 Þ ¼ 0; 2 2 12 s2 h
Ds unþ1 j
Dþ D unþ1 j
þ
s
ð2:8Þ
þ
ð2:9Þ 4
2
wherer ¼ hs2 . From the derivation of (2.8) and (2.9), we know that the truncation errors of (2.8) and (2.9) are Oðs2 þ h þ sh Þ 2 3 4 2 4 2 and O sh2 þ s2 þ h þ sh þ hs2 respectively. If r is any positive real number, the truncation errors become Oðh Þ and Oðh Þ. 3. Domain decomposition procedure and main results Suppose Jh ¼ 1; N s ¼ T. For simplicity, we will consider a domain decomposition which involves in decomposing ð0; 1Þ into only two subdomains, ð0; xÞ and ð x; 1Þ, where x ¼ xk for some integer k ð1 < k < J 1Þ. We use the explicit scheme (2.9) to compute the solution value unþ1 and the implicit scheme (2.8) to compute other solution values unþ1 ðj – kÞ respectively. j k The system can be written as
Lðunþ1 Þ ¼ 0; j unþ1 j
¼ 0;
1 6 n 6 N 1;
j ¼ 0; J;
0 < j < J;
ð3:10Þ
6172
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181
where the linear operator L is
Lðunþ1 Þ¼ j
8 unþ1 2un þun1 Dþ D unþ1 j j j j > nþ1 s h2 > ; > Ds uj h2 þ 2 þ 12 s2 > > < Dþ D unþ1
2
j – k;
unþ1 2un þun1
j j j h Ds unþ1 h2 j þ 2s þ 12 > j s2 > > > > : rðD un D unþ1 Þ rðD un D unþ1 Þ; s jþ1 s jþ1 s j1 s j1
ð3:11Þ
: j ¼ k:
The resulting system of equations decouples into two disjoint sets of equations corresponding to the subdomains. These systems can be solved in parallel. We are now in position to state two main theorems of this paper, which will be proved in the forthcoming sections. For the discrete function un ¼ funj jj ¼ 0; 1; . . . ; J; un0 ¼ unJ ¼ 0g, define
kun k2 ¼
J1 X junj j2 h;
n 2 kDþ un k2 ¼ sumJ1 j¼0 jDþ uj j h:
j¼1
We have the following theorems. Theorem 3.1. (Stability). For any given r > 0, the finite difference solutions of the parallel scheme (3.10)–(3.11) satisfy
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 þ maxkDþ u k 6 þ 1ðkDþ u1 k þ ku1 u0 kÞ: n 2r 12r2 n
Theorem 3.2. (Convergence). Let enj ¼ Uðxj ; tn Þ unj . For any given r > 0, the finite difference solutions of the parallel scheme (3.10)–(3.11) satisfy 4
maxkDþ en k 6 CðkDþ e1 k þ ke1 e0 k þ h Þ; n
where C is a positive constant independent of h and s. Since the finite difference scheme (3.10) has three time levels, besides taking u0j ¼ U 0 ðjhÞ, thus e0j ¼ 0, we need to find other methods to solve u1j ; 8j. In order to match fourth order accuracy, we can use either fourth order explicit schemes, such 4 as the impact scheme, or the high-order parallel iterative method [10] to compute u1j , i.e., let je1j j ¼ Oðh Þ. Then according to Theorem 3.2, scheme (3.10) will reach third-order accuracy. 4. Extension to two dimensional case In this section, Uðx; y; tÞ will be a solution of the following Dirichlet boundary problem on X ¼ ð0; 1Þ ð0; 1Þ,
@U DU ¼ 0; @t
ðx; yÞ 2 X;
t 2 ð0; T;
ð4:12Þ
Uðx; y; tÞ ¼ 0;
ðx; yÞ 2 @ X;
t 2 ð0; T;
ð4:13Þ
Uðx; y; 0Þ ¼ U 0 ðx; yÞ;
ðx; yÞ 2 X;
ð4:14Þ 2
2
where the initial function U 0 ðx; yÞ satisfies the boundary condition, and DU ¼ @@xU2 þ @@yU2 . We extend the domain decomposition method stated in Section 3 to the above problem. Take
X1 ¼ fðx; yÞ 2 X : x < xg; X2 ¼ fðx; yÞ 2 X : x > xg: Let xi ¼ ih same as in Section 3, and let yj ¼ jh. Suppose that there exists an integer k such that x ¼ xk (see Fig. 1). In analogy with Section 3 we call points ðxi ; yj ; tn Þ as boundary points if ðxi ; yj Þ 2 @ X and interface points if i ¼ k. Otherwise, we call them interior points. The values uni;j will approximate Uðxi ; yj ; t n Þ U ni;j . We denote
Dþ D ui;j ¼ uiþ1;j 2ui;j þ ui1;j ; and
Ds uni;j ¼
uni;j un1 i;j
s
:
dþ d ui;j ¼ ui;jþ1 2ui;j þ ui;j1
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181
6173
Fig. 1. Domain decomposition for two dimensional case.
Then
Dþ D U i;j 2
h dþ d U i;j
2
¼
@2U h @4U 4 ðxi ; yj ; tÞ þ ðxi ; yj ; tÞ þ Oðh Þ; 2 @x 12 @x4 2
¼
2
@2U h @4U 4 ðxi ; yj ; tÞ þ ðxi ; yj ; tÞ þ Oðh Þ; 2 @y 12 @y4
h dþ d Dþ D U i;j 4
h
¼
@4U 2 ðxi ; yj ; tÞ þ Oðh Þ: @x2 @y2
Thus we have
Dþ D U i;j þ dþ d U i;j
1 d d D D U 2 þ þ i;j 6h ! 2 h @4U @4U @4U 4 ¼ DUðxi ; yj ; tÞ þ ðxi ; yj ; tÞ þ 2 2 2 ðxi ; yj ; tÞ þ 4 ðxi ; yj ; tÞ þ Oðh Þ @x @y @y 12 @x4 2
h
þ
2
¼ DUðxi ; yj ; tÞ þ
h @2U 4 ðxi ; yj ; tÞ þ Oðh Þ: 12 @t2
ð4:15Þ
Moreover,
Ds U nþ1 ¼ i;j
@U s @2U ðxi ; yj ; t nþ1 Þ ðxi ; yj ; t nþ1 Þ þ Oðs2 Þ: @t 2 @t 2
ð4:16Þ
Subtract (4.16) with (4.15), we obtain
Ds U nþ1 i;j
Dþ D U nþ1 þ dþ d U nþ1 i;j i;j h
2
1
d d D D U nþ1 2 þ þ i;j
6h
! 2 h @2U 4 ¼ ðxi ; yj ; tnþ1 Þ þ Oðs2 þ h Þ: þ 2 12 @t2
s
Therefore, the fourth order finite difference scheme for (4.12) is
Ds unþ1 i;j
Dþ D unþ1 þ dþ d unþ1 i;j i;j h
2
1 6h
d d D D unþ1 2 þ þ i;j
! nþ1 2 n n1 h ui;j 2ui;j þ ui;j þ ¼ 0; þ 2 12 s2
s
ð4:17Þ
which is a nine-point finite difference scheme with three levels (see Fig. 2). Denote
ui;j ¼ ui1;j1 þ ui1;jþ1 þ uiþ1;j1 þ uiþ1;jþ1 4ui;j D 2 1 ¼ Dþ D ui;j þ dþ d ui;j þ dþ d Dþ D ui;j : 2 Then (4.17) is equivalent to
Ds unþ1 i;j
unþ1 þ 2Dþ D unþ1 þ 2dþ d unþ1 D i;j i;j i;j 3h
2
! nþ1 2 n n1 h ui;j 2ui;j þ ui;j þ ¼ 0; þ 2 12 s2
s
n n1 which can be used for computing unþ1 when i – k. When i ¼ k; unþ1 i;j i1;j in (4.17) are replaced by 2ui1;j ui1;j . In this case, we nþ1 need to solve a tridiagonal matrix using Thomas algorithm for computing uk;j .
6174
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181
Fig. 2. Nine point scheme.
It is straightforward to extend the two-subdomain results to many subdomains on x direction by cutting the whole domain into vertical strips. Moreover, we can also extend the domain decomposition scheme to three dimensional case by dividing x axis into many subdomains. 5. Proof of stability and convergence We first state three auxiliary lemmas. The stability and convergence results are then derived. Lemma 5.1 (Discrete Poincare Inequality). For the discrete function un ¼ funj jj ¼ 0; 1; . . . ; J; un0 ¼ unJ ¼ 0g, there exists
kun k 6 k
Dþ un k: h
Lemma 5.2 (Discrete Green Theorem). If uj and
v j are discrete functions on the set fxj jj ¼ 0; 1; . . . ; Jg, then we have
J1 J X X uj Dþ v j ¼ v j D uj u0 v 1 þ uJ v J : j¼1
j¼1
Lemma 5.3. For any given r > 0; f 2 L1 ð0; T; L2 ð0; 1ÞÞ, the finite difference solutions of system Lðunþ1 Þ ¼ fjnþ1 with the Dirichlet j boundary condition satisfy
maxkDþ un k2 6 kDþ u1 k2 þ n
1 1 2 þ þ 1 ku1 u0 k2 þ Th maxkf n k2 : 2 n 2r 12r
Lemma 5.1 and 5.2 are proved in [11]. Before proving the stability and convergence, we will give the proof of Lemma 5.3. 5.1. Proof of Lemma 5.3 Denoting wnþ1 ¼ unþ1 unj , we can rewrite Lðunþ1 Þ ¼ fjnþ1 with the Dirichlet boundary condition as j j j
8 nþ1 nþ1 1 w r Dþ D unþ1 þ 12 þ 12r ðwj wnj Þ ¼ fjnþ1 s; > j > > j > nþ1 > nþ1 nþ1 1 1
j – k;
j ¼ k;
j
Multiplying the above equations by wnþ1 h; j ¼ 1; 2; . . . ; J 1 and summing them up with respect to j, we have j
J1 J1 J1 X X 1 1 X nþ1 n nþ1 nþ1 n nþ1 þ ðwnþ1 Þ2 h r wnþ1 Dþ D unþ1 hþ ðw wnj Þwnþ1 h rwnþ1 j j j j k ðwkþ1 wkþ1 Þh rwk ðwk1 wk1 Þh 2 12r j¼1 j j¼1 j¼1 ¼
J1 X fjnþ1 wnþ1 sh: j j¼1
ð5:18Þ
6175
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181
From Lemma 5.2, we have J1 J J1 J1 J1 X X 1X 1X 1X wnþ1 Dþ D unþ1 ¼ D unþ1 D wnþ1 ¼ ðDþ unþ1 Þ2 þ ðDþ unj Þ2 ðDþ unþ1 Dþ unj Þ2 : j j j j j j 2 2 2 j¼1 j¼1 j¼0 j¼0 j¼0
Then (5.18) becomes
i r XJ1 XJ1 r hXJ1 1 1 nþ1 2 n 2 nþ1 n 2 þ kwnþ1 k2 þ ð D u Þ h ð D u Þ h þ ð D u D u Þ h þ þ þ þ þ j j j j j¼0 j¼0 j¼0 2 2 2 12r X
X X J1 J1 J1 1 1 1 nþ1 nþ1 ðwnþ1 Þ2 h ðwnj Þ2 h þ ðwnþ1 wnj Þ2 h þ rhwk ðwnþ1 j j kþ1 þ wk1 Þ j¼1 j¼1 j¼1 2 2 2 J1 X nþ1 ¼ rhwk ðwnkþ1 þ wnk1 Þ þ wnþ1 fjnþ1 sh j j¼1 2 2 n ðwnþ1 ðwnþ1 Þ2 þ ðwnk1 Þ2 k Þ þ ðwkþ1 Þ 6 rh þ k 2 2 2 6 rhðwnþ1 k Þ þ rh
! þ
J1 X ðwnþ1 Þ2 þ ðfjnþ1 sÞ2 j
2
j¼1
h
ðwnkþ1 Þ2 þ ðwnk1 Þ2 1 nþ1 2 s2 nþ1 2 þ kw k þ kf k : 2 2 2
ð5:19Þ
Noting that
" # 2 nþ1 2 ðwnþ1 r XJ1 r XJ1 nþ1 2 kþ1 Þ þ ðwk1 Þ nþ1 n 2 nþ1 nþ1 nþ1 2 þ rhðwnþ1 ð D D ð D Þ h þ rh þ uj þ uj Þ h þ rhwk ðwkþ1 þ wk1 Þ ¼ þ wj k Þ ; j¼0 j¼0;j–k;k1 2 2 2 we can simplify (5.19) as follows
i r XJ1 1 nþ1 2 r h 1 1 nþ1 2 kw k þ kDþ unþ1 k2 kDþ un k2 þ þ ðkwnþ1 k2 kwn k2 þ kwnþ1 wn k2 Þ ð D w Þ h þ þ j j¼0;j–k;k1 2 2 2 4 24r " # 2 nþ1 2 ðwnþ1 ðwnkþ1 Þ2 þ ðwnk1 Þ2 s2 kþ1 Þ þ ðwk1 Þ þ rh 6 kf nþ1 k2 : 2 2 2 Since kwnþ1 k2 P 0;
PJ1
nþ1 2 Þ h j¼0;j–k;k1 ðDþ wj
P 0 and kwnþ1 wn k2 P 0, then we have
" # 2 nþ1 2 i 1 ðwnþ1 ðwnkþ1 Þ2 þ ðwnk1 Þ2 rh 1 s2 kþ1 Þ þ ðwk1 Þ nþ1 2 n 2 nþ1 2 n 2 6 kf nþ1 k2 : kDþ u k kDþ u k þ þ ðkw k kw k Þ þ rh 2 4 24r 2 2 2 Summing up with respect to n, we get
2 ðwnþ1 Þ2 þ ðwnþ1 ðw1 Þ2 þ ðw1k1 Þ2 r 1 1 r 1 1 k1 Þ kDþ unþ1 k2 þ þ kwnþ1 k2 þ rh kþ1 þ kw1 k2 þ rh kþ1 6 kDþ u1 k2 þ 2 4 24r 2 4 24r 2 2 Ts maxkf n k2 : þ 2 n Thus
kDþ unþ1 k2 6 kDþ u1 k2 þ Lemma 5.3 is proved.
1 1 2 þ 1 ku1 u0 k2 þ Th maxkf n k2 : þ n 2r 12r 2
h
5.2. Proof of stability Let f ¼ 0 in Lemma 5.3, thus Theorem 3.1 is proved.
h
5.3. Proof of convergence According to the parallel scheme (3.10)–(3.11), the errors enþ1 ð1 6 n < NÞ follows the below equations: j
! nþ1 2 n n1 h ej 2ej þ ej þ ¼ Gnþ1 ; þ j 2 2 2 12 s h ! Dþ D enþ1 2enj þ en1 s h2 enþ1 j j j Ds enþ1 þ þ j 2 2 12 s2 h
Ds enþ1 j
Dþ D enþ1 j
s
2
nþ1 n nþ1 rðDs enjþ1 Ds enþ1 h þ Gnþ1 ; j jþ1 Þ rðDs ej1 Ds ej1 Þ ¼ U
j – k;
j ¼ k;
enþ1 ¼ 0; j
j ¼ 0; J;
6176
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181 4
2
where Unþ1 ¼ 2r 2 @@tU2 ðxk ; tnþ1 Þ; jGnþ1 j ¼ Oðh Þ. In order to get the error estimates for enþ1 , we assume that enþ1 ¼ pnþ1 þ qnþ1 , j j j j j nþ1 where pnþ1 and q are the solutions of the following problems respectively. j j problem I
! nþ1 2 n n1 h pj 2pj þ pj ¼ Gnþ1 ; j – k; j 2 2 2 12 s h ! Dþ D pnþ1 2pnj þ pn1 s h2 pnþ1 j j j nþ1 n nþ1 Ds pnþ1 þ rðDs pnjþ1 Ds pnþ1 ; þ j jþ1 Þ rðDs pj1 Ds pj1 Þ ¼ Gj 2 2 12 s2 h
Ds pnþ1 j
Dþ D pnþ1 j
pnþ1 ¼ 0; j
s
þ
þ
j ¼ k;
ð5:20Þ
j ¼ 0; J;
p0j ¼ e0j ;
p1j ¼ e1j ; 8j:
problem II
Ds qnþ1 j Ds qnþ1 j
Dþ D qnþ1 j 2
h
Dþ D qnþ1 j 2
h
s
þ
2
s
þ
2
þ þ
! nþ1 2 n n1 h qj 2qj þ qj ¼ 0; 2 12 s ! nþ1 2 n n1 h qj 2qj þ qj 12
s2 2
nþ1 n nþ1 rðDs qnjþ1 Ds qnþ1 h ; jþ1 Þ rðDs qj1 Ds qj1 Þ ¼ U
qnþ1 j
¼ 0;
j – k;
j ¼ k;
j ¼ 0; J;
q0j ¼ q1j ¼ 0; 8j: From Lemma 5.3, we can obtain the estimate for pnþ1 , i.e., j
1 1 2 þ kDþ pn k2 6 kDþ p1 k2 þ þ 1 kp1 p0 k2 þ Th maxkGn k2 2 n 2r 12r 1 1 10 6 þ 1 ðkDþ e1 k þ ke1 e0 kÞ2 þ Oðh Þ: þ 2r 12r 2
ð5:21Þ
nþ1 For estimating qnþ1 , we first consider q ð1 6 n 6 N 1Þ which satisfies the following equations j j
nþ1 Dþ D q j 2
h nþ1 Dþ D q k
nþ1 q 0
¼ 0;
j – k; 2
¼ Unþ1 h ; 2 h nþ1 ¼q ¼ 0: J
nþ1 Then the formula of q is j
nþ1 q j
8 0; j ¼ 0; J; > > < j PJ nþ1 4 ¼ J kþ1 U h ; 1 6 j 6 k; > P > : Jj k1 Unþ1 h4 ; k < j < J: i¼1 J
Hence
nþ1 k ¼ Oðh4 Þ and kq nþ1 q n k ¼ Oðsh3 Þ: kDþ q ~nþ1 nþ1 Define q ¼ qnþ1 q , then we have j j j
! nþ1 2 ~ j 2q ~nj þ q ~n1 h q j þ ¼ Rnþ1 ; j – k; þ j 2 2 2 12 s h ! ~nþ1 ~nj þ q ~n1 Dþ D q 2q s h2 q~nþ1 j j j nþ1 ~nþ1 ~njþ1 Ds q ~nþ1 ~n ~nþ1 Ds q þ rðDs q ; þ j jþ1 Þ rðDs qj1 Ds qj1 Þ ¼ Rj 2 2 2 12 s h ~nþ1 Ds q j
where
~nþ1 Dþ D q j
s
~nþ1 ¼ 0; q j
j ¼ 0; J;
~0j ¼ q 0j ; q
~1j ¼ q 1j ; 8j; q
j ¼ k;
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181
Rnþ1 j
6177
8 qnþ1 2qn þqn1 j j j h2 > nþ1 > Ds q þ 2s þ 12 ; j – k; > j s2 > < nþ1 2q n þq n1 q 2 ¼ D q j j j nþ1 h njþ1 Ds q nþ1 þ 2s þ 12 rðDs q s j > jþ1 Þ s2 > > > : nþ1 n j1 Ds q j1 Þ; j ¼ k; rðDs q 3
~nþ1 k. and kRnþ1 k ¼ Oðh Þ. From Lemma 5.3, we can obtain the estimate of kDþ q
~nþ1 k2 6 kDþ q ~1 k2 þ kDþ q
1 1 ~1 q ~0 k2 þ Th2 maxkRn k2 ¼ Oðh8 Þ: þ 1 kq þ n 2r 12r 2
Thus, 4
~nþ1 k þ kDþ q nþ1 k ¼ Oðh Þ: kDþ qnþ1 k 6 kDþ q Combining with (5.21), we get 4
kDþ enþ1 k 6 kDþ pnþ1 k þ kDþ qnþ1 k 6 CðkDþ e1 k þ ke1 e0 k þ h Þ; where C is a positive constant independent of h and s. Theorem 3.2 is proved.
h
6. Numerical experiments In this section, some numerical results are presented to show the stability, accuracy, and parallelism of the scheme described above, and the computational costs are also presented. All the experiments are run on a cluster consisting of a manager that uses one core of a Xeon 5410 processor and up to 32 computing nodes, each containing two Xeon 5410 processors running 64-bit Linux, i.e., each node consists of 8 processing cores. 6.1. One dimensional numerical experiments We consider the problem defined in equations (1.1)–(1.3) with U 0 ðxÞ ¼ sinðpxÞ. Obviously the exact solution of the equa2 tions is Uðx; tÞ ¼ ep t sinðpxÞ. First, we verify the stability of the scheme by taking the step size h ¼ 103 ; r ¼ 1; 10; 100; 1000 with two subdomains. Fig. 3 clearly shows that the norm of un doesn’t occur blowing up even if r is large enough. This explains the unconditional stability of the scheme. Second, we examine the numerical errors in the solutions. Table 1 shows that the errors for each case are roughly of the 3 same order of magnitude, and the errors appear to be Oðh Þ in each case. We also compare them with those resulting from
Fig. 3. The infinity norm of un v.s. t.
6178
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181
Table 1 Numerical errors for different grid length. 2 Processors
4 Processors
J
ku U k=h
r ¼ 10; T ¼ 1 1000 2000 4000 8000
J
2.06821 2.06715 2.06316 1.94919
r ¼ 100; T ¼ 1 1000 2000 4000 8000
8 Processors ku U k=h
J
kun U n k=h
1000 2000 4000 8000
4.05448 4.05336 4.04938 3.93590
1000 2000 4000 8000
8.07385 8.07276 8.06881 7.95632
2.06703e2 2.06729e2 2.06690e2 2.06820e2
1000 2000 4000 8000
4.04766e2 4.05299e2 4.05307e2 4.05434e2
1000 2000 4000 8000
8.04137e2 8.06955e2 8.07216e2 8.07356e2
r ¼ 1000; T ¼ 1 1000 2000 4000 8000
2.01883e4 2.03740e4 2.06360e4 2.06634e4
1000 2000 4000 8000
4.04424e4 3.93818e4 4.03957e4 4.05137e4
1000 2000 4000 8000
8.03459e4 8.03843e4 8.01579e4 8.06558e4
r ¼ 10000; T ¼ 1 1000 2000 4000 8000
2.02383e6 2.04768e6 2.5982e6 2.06920e6
1000 2000 4000 8000
4.05348e6 4.10234e6 4.00321e6 4.01203e6
1000 2000 4000 8000
8.00374e6 8.04875e6 8.00829e6 8.06558e6
n
n
3
n
n
3
3
Table 2 Numerical errors for the scheme in [8]. 2 Processors
4 Processors
J
kun U n k=h
r ¼ 10; T ¼ 1 1000 2000 4000 8000
8 Processors
J
kun U n k=h
J
kun U n k=h
1.02839 1.00293 1.12303 1.00290
1000 2000 4000 8000
2.03842 2.20394 2.04095 2.03938
1000 2000 4000 8000
4.09382 4.02910 4.10293 3.98201
r ¼ 100; T ¼ 1 1000 2000 4000 8000
1.09282e1 1.02938e1 1.10293e1 1.00927e1
1000 2000 4000 8000
2.09384e1 2.07261e1 2.01628e1 2.10928e1
1000 2000 4000 8000
4.08271e1 4.08726e1 4.09281e1 4.09817e1
r ¼ 1000; T ¼ 1 1000 2000 4000 8000
1.03928e2 1.09282e2 1.07621e2 1.01822e2
1000 2000 4000 8000
2.00971e2 2.00761e2 2.00127e2 2.02872e2
1000 2000 4000 8000
4.02893e2 4.02934e2 4.00029e2 4.07625e2
r ¼ 10000; T ¼ 1 1000 2000 4000 8000
1.00121e3 1.00273e3 1.02922e3 1.00461e3
1000 2000 4000 8000
2.04759e3 2.07162e3 2.00937e3 2.00183e3
1000 2000 4000 8000
4.01123e3 4.08211e3 4.00182e3 4.09837e3
2
2
Table 3 Comparison of scheme. Processors n
1 2 4 8 16 32 64
The third order scheme
The second order scheme in [8]
Time tn (s)
Speed-up t1 =tn
Time tn (s)
Speed-up t 1 =t n
4862.66 2450.82 1250.29 621.74 325.05 165.85 100.45
– 1.98 3.89 7.82 14.96 29.32 48.40
4564.28 2340.66 1164.36 598.98 300.84 157.82 98.97
– 1.95 3.92 7.62 15.21 28.92 46.12
2
6179
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181 Table 4 Comparison of time for two dimensional case. Processors n
Time tn (s)
Speed-up t1 =tn
1 2 4 8 16 32 64
10 h 23 m 19 s 6 h 01 m 29 s 2 h 52 m 14 s 1 h 23 m 54 s 42 m 47 s 21 m 23 s 12 m 58 s
– 1.72 3.61 7.41 14.53 29.06 47.91
Table 5 Numerical errors for two dimensional cases. 8 Processors
16 Processors
Jx Jy
ku U k=h
r ¼ 10; T ¼ 1 100 100 200 200 400 400
10.07867 10.06592 10.02345
r ¼ 100; T ¼ 1 100 100 200 200 400 400
n
n
1.09837e4 1.05629e4 1.07820e4
3
32 Processors
Jx Jy
ku U k=h
100 100 200 200 400 400
40.09821 40.06710 40.01023
100 100 200 200 400 400
n
n
3
4.09272e4 4.09983e4 4.05730e4
Jx Jy
kun U n k=h
100 100 200 200 400 400
160.09802 160.02893 160.07802
100 100 200 200 400 400
3
1.60029e5 1.60789e5 1.60892e5
Fig. 4. Initial condition for the Lotka–Volterra system.
the second order scheme in [8], which are shown in Table 2. Table 1 and Table 2 show that our domain decomposition algorithm has higher order accuracy. Third, we test the speed-up for the scheme. Here we take h ¼ 105 ; s ¼ 106 , i.e., r ¼ 104 ; T ¼ 0:5 and list the time for computing and the speed-up in Table 3 for our domain decomposition algorithm and the scheme in [8]. Table 3 shows that our domain decomposition algorithm is competitive with the existing method in computing time.
6180
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181 Table 6 Comparison of time for Lotka–Volterra system. Processors n
Time tn (s)
Speed-up t 1 =tn
1 2 4 8 16 32 64 128 256
21 h 51 m 32 s 11 h 45 m 8 s 6 h 12 m 36 s 2 h 58 m 57 s 1 h 33 m 37 s 45 m 28 s 28 m 12 s 14 m 33 s 7 m 28 s
– 1.86 3.52 7.37 14.16 28.85 46.51 90.12 175.61
Fig. 5. Numerical solutions for T ¼ 1.
6.2. Two dimensional numerical experiments We consider the problem defined in equations (4.12)–(4.14) with U 0 ðx; yÞ ¼ sinðpxÞ sinðpyÞ. Obviously the exact solution 2 of the equations is Uðx; y; tÞ ¼ e2p t sinðpxÞ sinðpyÞ. The computing time and the speed-up of the scheme are shown in Table 4 3 2 with h ¼ 10 ; s ¼ 10 and T ¼ 1. Numerical errors are shown in Table 5 and also achieve the third order accuracy.
6.3. Application to Lotka–Volterra system The competitive Lotka–Volterra equations are a simple model of the population dynamics of species competing for some common resource. Given two populations, u and v, with logistic dynamics, the Lotka–Volterra formulation adds an additional term to account for the species’ interactions. Thus the competitive Lotka–Volterra equations are:
ut ¼ Du þ uð1 v Þ;
ðx; yÞ 2 X
v t ¼ Dv v ð1 uÞ; ðx; yÞ 2 X v ¼ u ¼ 0; ðx; yÞ 2 @ X: We choose random values for uðx; yÞ and v ðx; yÞ as initial conditions shown in Fig. 4. We take h ¼ 0:001, s ¼ 0:01; T ¼ 1 and show the computing time in Table 6. Fig. 5 shows the numerical solution for T ¼ 1.
W. Hao, S. Zhu / Applied Mathematics and Computation 219 (2013) 6170–6181
6181
7. Conclusion We presented domain decomposition finite difference schemes with unconditional stability and third-order accuracy for the parabolic system. Error estimate and stability of the numerical solutions have been derived for one dimensional case. The scheme is easy to implement the parallelism and is extended in two and three dimensional case. The numerical results demonstrate the good performance of the parallel scheme, namely, unconditional stability, the third order accuracy and high degree of parallelism. References [1] C.N. Dawson, Q. Du, T.F. Dupont, A finite difference domain decomposition algorithm for numerical solution of the Heat equation, Math. Comp. 57 (1991) 63–71. [2] G. Yuan, S. Zhu, L. Shen, Domain decomposition algorithm based on the group explicit formula for the heat equation, Int. J. Comput. Math 82 (2005) 1295–1306. [3] S. Zhu, Z. Yu, J. Zhao, A high-order parallel finite difference algorithm, Appl. Math. Comp. 183 (2006) 365–372. [4] C.N. Dawson, T.F. Dupont, Explicit/implicit conservative domain decomposition procedures for parabolic problems based on block-centered finite differences, SIAM J. Numer. Anal. 31 (1994) 1045–1061. [5] G. Yuan, L. Shen, Stability and convergence of the explicit–implicit conservative domain decomposition procedure for parabolic problems, Comput. Math. Appl. 47 (2004) 793–801. [6] H. Shi, H. Liao, Unconditional stability of corrected explicit–implicit domain decomposition algorithms for parallel approximation of heat equations, SIAM J. Numer. Anal. 44 (2006) 1584–1611. [7] H. Liao, H. Shi, Z. Sun, Corrected explicit–implicit domain decomposition algorithms for two-dimensional semilinear parabolic equations, Sci. China Series A 52 (2009) 2362–2388. [8] Z. Sheng, G. Yuan, X. Hang, Unconditional stability of parallel difference schemes with second order accuracy for parabolic equation, Appl. Math. Compt. 184 (2007) 1015–1031. [9] S. Zhu, Conservative domain decomposition procedure with unconditional stability and second-order accuracy, Appl. Math. Compt. 216 (2010) 3275– 3282. [10] W. Hao, S. Zhu, Parallel iterative methods for parabolic equations, Int. J. Comput. Math 86 (2009) 431–440. [11] Y. Zhou, Applications of discrete functional analysis to the finite difference method, International Academic publishers (1991).