Second-order LOD multigrid method for multidimensional Riesz fractional diffusion equation

0.27 MB·English

by  Minghua Chen

Second-order LOD multigrid method for
multidimensional Riesz fractional diffusion equation
Minghua Chen, Yantao Wang, Xiao Cheng, Weihua Deng∗
School of Mathematics and Statistics, Lanzhou University, Lanzhou 730000, P.R. China
3
1
0
2 Abstract

n
a Weproposealocallyonedimensional (LOD)finitedifference methodformultidimensional
J
Riesz fractional diffusion equation with variable coefficients on a finite domain. The
2
1
numerical method is second-order convergent in both space and time directions, and its

] unconditional stability is strictly proved. Comparing with the popular first-order finite
A
N difference method for fractional operator, the form of obtained matrix algebraic equation
h. is changed from (I A)uk+1 = uk + bk+1 to (I A)uk+1 = (I + B)uk + ˜bk+1/2; the
t − −
a three matrices A, A and B are all Toeplitz-like, i.e., they have completely same structure
m
e e
and the computational count for matrix vector multiplication is (NlogN); and the
[
e e O

1 computational costs for solving the two matrix algebraic equations are almost the same.
v
3 The LOD-multigrid method is used to solve the resulting matrix algebraic equation, and
4
the computational count is (NlogN) and the required storage is (N), where N is the
6
O O
2
number of grid points. Finally, the extensive numerical experiments are performed to
.
1
0 show the powerfulness of the second-order scheme and the LOD-multigrid method.
3
1
Keywords: Riesz fractional diffusion equation; Second-order discretization; Toeplitz and
:
v
Circulant matrices; Multigrid method
i
X
r Mathematics Subject Classification (2010): 35R11, 65M06, 65M55
a
1. Introduction
In recent years considerable interests in fractional calculus have been stimulated by
the applications in physical, chemical, biological, and engineering, etc., areas [9]. The def-
initions of fractional calculus are versatile, e.g., Riemann-Liouville derivative, Gru¨nwald-
Letnikov derivative, Caputo derivative, Weyl derivative and Riesz derivative et al [11, 14],
and they are not completely equivalent. Depending on the particular applied field, some-
times one of its definitions is more popular than others. For example, the Riesz fractional
∗Corresponding author. E-mail: [email protected].
Preprint submitted to Elsevier
derivative appears in the continuous limit of lattice models with long-range interactions
[19]. This paper focuses on the multidimensional Riesz fractional diffusion equations.
Nowadays, the finite difference discretization for space fractional derivatives is expe-
riencing rapid development, including the Riesz fractional derivative; such as, Yang et
al numerically study the Riesz space fractional PDEs with two different fractional orders
1 < α 2 and 0 < β < 1 [24]; Zhuang et al consider a variable-order fractional advection-
≤
diffusion equation with a nonlinear source term on a finite domain [26]. In the last two
years, for the space fractional derivatives, we notice that two different second-order dis-
cretization schemes are developed [17, 20]; even the third-order discrezation scheme is
obtained [25] if a compact difference operator is performed on the discrezation scheme
given in [20].
Another topicrelatedto effectively solving theequations involving fractionaloperators
is about how to efficiently solve the resulting matrix algebraic equations. The ‘unlucky’
thing is that the matrix in the matrix algebraic equation is usually full because of the
nonlocal properties of the fractional operators, and the ‘lucky’ thing, as pointed out in
[12, 21, 22], is that the matrix has some special structure, i.e., the matrix is Toeplitz-
like matrix, and the count of its matrix vector multiplication is (NlogN) by using
O
the constructed circulant matrix and fast Fourier transform, and the required storage
is (N), where N is the number of grid points. Pang and Sun [12] successfully use
O
the multigrid method (MGM) to efficiently solve the resulting matrix algebraic equation
of the one dimensional fractional diffusion equation by using first-order discretization
scheme [10]. Here we further extend the MGM to solve the matrix algebraic equation of
the multidimensional Riesz fractional diffusion equation discretized by the second-order
scheme.
We use the LOD strategy to solve the multidimensional Riesz fractional diffusion
equation. The LOD methods include alternating direction (AD) methods and fractional
step procedures [6]. The AD methods were first introduced in three papers [4, 7, 13]
by Douglas, Peaceman, and Rachford. The Peaceman and Rachford (PR-AD) method
works well for two-dimensional problems. However, it can not be extended to higher
dimensional problems. Douglas (D-AD) method [4, 5, 6] are valid for any dimensional
equations. And PR-AD and D-ADareequivalent in two-dimensional problems. BothPR-
AD and D-AD schemes are used in this paper to discretize the multidimensional Riesz
fractional diffusion equation. In eachdimension, the obtainedmatrixalgebraic equationis
solved by MGM. Although the spacial fractional derivative is discretized by second-order
scheme, for any single dimension the form of the obtained matrix algebraic equation is
(I A)uk+1 = (I + B)uk + ˜bk+1/2, in fact the corresponding form for the first-order
−
e e
2
discretization scheme is (I A)uk+1 = uk +bk+1; the three matrices A, A and B are all
−
Toeplitz-like, and they have completely same structure and the computational count for
e e
matrix vector multiplication is (NlogN); and the computational costs for solving the
O
two matrix algebraic equations are almost the same. In other words, the second-order
scheme improves the accuracy but almost without increasing the computational cost.
More concretely, in this paper using the second-order accurate and unconditionally
stable computational scheme and LOD-MGM, we solve the following multidimensional
variable coefficients Riesz fractional diffusion equation with the computational count
(NlogN) and storage (N),
O O
∂u(x,y,z,t) ∂αu(x,y,z,t) ∂βu(x,y,z,t)
= c(x,y,z,t) +d(x,y,z,t)
∂t ∂ x α ∂ y β

| | | |
 ∂γu(x,y,z,t)

 +e(x,y,z,t) ∂ z γ +f(x,y,z,t), (1.1)

 | |


u(x,y,z,0) = u (x,y,z) for (x,y,z) Ω,
0
∈

 u(x,y,z,t) = 0 for (x,y,z,t) ∂Ω (0,T],
 ∈ ×



in the domain Ω = (x ,x ) (y ,y ) (z ,z ), 0 < t T, where the orders of the Riesz
L R L R L R
× × ≤
fractional derivatives are 1 < α,β,γ < 2; f(x,y,z,t) is a source term and the variable
coefficients c(x,y,z,t) 0, d(x,y,z,t) 0, e(x,y,z,t) 0; theRiesz fractional derivative
≥ ≥ ≥
for n N, n 1 < ν n, is defined as [3, 19]
∈ − ≤
∂νu(x,y,z,t)
= κ Dν + Dν u(x,y,z,t), (1.2)
∂ x ν − ν xL x x xR
| |
(cid:0) (cid:1)
where the coefficient κ = 1 , and
ν 2cos(νπ/2)
1 ∂n x
Dνu(x,y,z,t) = (x ξ)n−ν−1u(ξ,y,z,t)dξ, (1.3)
xL x Γ(n ν)∂xn −
− ZxL
( 1)n ∂n xR
Dν u(x,y,z,t) = − (ξ x)n−ν−1u(ξ,y,z,t)dξ, (1.4)
x xR Γ(n ν)∂xn −
− Zx
are the left and right Riemann-Liouville space fractional derivatives, respectively.
The outline of this paper is as follows. In the next section, we introduce the second-
order finite difference discretizations for the Riesz fractional derivatives; and the full
discretization of (1.1) is derived, where the Crank-Nicolson scheme and LOD method
are combined together. We theoretically prove the presented finite difference scheme is
unconditionally stable in Section 3. In Section 4 we propose a V-cycle LOD-MGM for
the resulting system of (1.1). To show the powerfulness of the second-order scheme and
LOD-MGM, the extensive numerical experiments are performed in Section 5. Finally, we
conclude the paper with some remarks in the last section.
3
2. Derivation of the finite difference scheme
In this section, we derive the full discretization schemes of (1.1). The first subsection
introduces the second-order finite difference discretizations for the Riesz fractional deriva-
tives in a finite domain. Then in the second subsection, we present the scheme for the
one dimensional case of (1.1). The third and fourth subsections detailedly provide the
two dimensional case of (1.1) and (1.1) itself, respectively.
2.1. Discretizations for the Riesz fractional derivatives
Take the mesh points x = x +i∆x,i = 0,1,...,N , y = y +j∆y,j = 0,1,...,N ,
i L x j L y
z = z +l∆z,l = 0,1,...,N and t = k∆t,k = 0,1,...,N , where ∆x = (x x )/N ,
l L z k t R L x
−
∆y = (y y )/N , ∆z = (z z )/N , ∆t = T/N , i.e., ∆x, ∆y and ∆z are the uniform
R L y R L z t
− −
space stepsizes in the corresponding directions, ∆t the time stepsize. For ν (1,2), the
∈
left and right Riemann-Liouville space fractional derivatives (1.3) and (1.4) have the
second-order approximation operators δ uk and δ uk , respectively, given in a
ν,+x i,j,l ν,−x i,j,l
finite domain [3, 17], where uk denotes the approximated value of u(x ,y ,z ,t ).
i,j,l i j l k
The approximation operator of (1.3) is defined by [3, 17]
i+1
1
δ uk := gν uk , (2.1)
ν,+x i,j,l Γ(4 ν)(∆x)ν m i−m+1,j,l
− m=0
X
and there exists
Dνu(x,y,z,t) = δ uk + (∆x)2, (2.2)
xL x ν,+x i,j,l O
where
1, m = 0,
 4+23−ν, m = 1,
−

gν =  6 25−ν +33−ν, m = 2, (2.3)
m  −

 (m+1)3−ν 4m3−ν +6(m 1)3−ν
− −
 4(m 2)3−ν +(m 3)3−ν, m 3.

 − − − ≥


Analogously, the approximation operator of (1.4) is described as [3]

1 Nx−i+1
δ uk := gν uk , (2.4)
ν,−x i,j,l Γ(4 ν)(∆x)ν m i+m−1,j,l
− m=0
X
and it holds that
Dν u(x,y,z,t) = δ uk + (∆x)2, (2.5)
x xR ν,−x i,j,l O
where gν is defined by (2.3).
m
Combining (2.2) and (2.5), we obtain the approximation operator of the (Riemann-
Liouville) Riesz fractional derivative
4
∂νu(x,y ,z ,t )
j l k = κ Dν + Dν u(x,y ,z ,t )
∂|x|ν (cid:12)x=xi − ν xL x x xR j l k x=xi
(cid:0) (cid:1) (cid:12)
(cid:12) = κ δ +δ uk + (∆x)(cid:12)2
(cid:12) − ν ν,+x ν,−x i,j,l O
(cid:0) κ i+(cid:1)1 Nx−i+1
= − ν gν uk + gν uk + (∆x)2
Γ(4 ν)∆xν m i−m+1,j,l m i+m−1,j,l O
!
− m=0 m=0
X X
κ i+1 Nx
= − ν gν uk + gν uk + (∆x)2
Γ(4 ν)∆xν i−m+1 m,j,l m−i+1 m,j,l O
!
− m=0 m=i−1
X X
κ Nx
:= − ν gν uk + (∆x)2,
Γ(4 ν)∆xν i,m m,j,l O
− m=0
X
(2.6)
e
where
gν , m < i 1,
i−m+1 −
 gν +gν, m = i 1,
0 2 −

gν =  2gν, m = i, (2.7)
i,m  1


 gν +gν, m = i+1,
0 2
e
 gν , m > i+1,
 m−i+1



with i = 1,...,N 1, together with the Dirichlet boundary conditions that define uk
x −  0,j,l
and uk as appropriate.
Nx,j,l
Taking ν = 2, both Eq. (2.2) and (2.5) reduce to the following form
∂2u(x ,y,z,t) u(x ,y,z,t) 2u(x ,y,z,t)+u(x ,y,z,t)
i = i+1 − i i−1 + (∆x)2.
∂x2 (∆x)2 O
Similarly, it is easy to get the one-dimensional and two-dimensioanl case of (2.1)-(2.7).
2.2. Numerical scheme for 1D
Consider the one-dimensional Riesz fractional diffusion equation
∂u(x,t) ∂αu(x,t)
= c(x,t) +f(x,t). (2.8)
∂t ∂ x α
| |
In the time direction, we use the Crank-Nicolson scheme. Taking the uniform time
step ∆t and space step ∆x, and setting ck = c(x ,t ) and fk+1/2 = f(x ,t ), where
i i k i i k+1/2
t = (t +t )/2, the full discretization of (2.8) has the following form
k+1/2 k k+1
uk+1 uk κ ck+1/2 Nx uk +uk+1
i − i = − α i gα m m +fk+1/2. (2.9)
∆t Γ(4 α)∆xα i,m 2 i
− m=0
X
e
5
Then (2.9) can be expressed as
∆t ∆t
1 δ′′ uk+1 = 1+ δ′′ uk +∆tfk+1/2, (2.10)
− 2 α,x i 2 α,x i i
(cid:18) (cid:19) (cid:18) (cid:19)
where
κ ck+1/2 Nx κ ck+1/2 Nx
δ′′ uk := − α i gα uk ; δ′′ uk+1 := − α i gα uk+1.
α,x i Γ(4 α)∆xα i,m m α,x i Γ(4 α)∆xα i,m m
− m=0 − m=0
X X
e e
Putting ξk+1/2 = −∆tκαcik+1/2, the system of equations given by (2.10) takes the form
i 2Γ(4−α)∆xα
(I Ak+1/2)Uk+1 = (I +Ak+1/2)Uk +∆tFk+1/2, (2.11)
−
where I is the (N 1) (N 1) identity matrix,
x x
− × −
Uk = [uk,uk,...,uk ]T, Fk+1/2 = [fk+1/2,fk+1/2,...,fk+1/2]T,
1 2 Nx−1 1 2 Nx−1
and the discretizations at the interior x-gridpoints define the entries of the matrix Ak+1/2,
k+1/2
A for i = 1,...,N 1 and m = 1,...,N 1 are defined by
i,m x − x −
gα ξk+1/2, m < i 1,
i−m+1 i −
 (gα +gα)ξk+1/2, m = i 1,
0 2 i −

Ak+1/2 =  2gαξk+1/2, m = i, (2.12)
i,m  1 i

 (gα +gα)ξk+1/2, m = i+1,
0 2 i
 gα ξk+1/2, m > i+1.
 m−i+1 i




2.3. LOD scheme for 2D 
Consider the following two-dimensional Riesz fractional diffusion equation
∂u(x,y,t) ∂αu(x,y,t) ∂βu(x,y,t)
= c(x,y,t) +d(x,y,t) +f(x,y,t). (2.13)
∂t ∂ x α ∂ y β
| | | |
Analogously we still use the Crank-Nicolson scheme to do the discretization in time
direction. Taking uk as the approximated value of u(x ,y ,t ), ck = c(x ,y ,t ), dk =
i,j i j k i,j i j k i,j
k+1/2
d(x ,y ,t ), t = (t + t )/2, f = f(x ,y ,t ), ∆x = (x x )/N , and
i j k n+1/2 n n+1 i,j i j k+1/2 R − L x
∆y = (y y )/N , for the uniform space steps ∆x,∆y and the time stepsize ∆t, the
R L y
−
resulting discretization of (2.13) can be written as
uk+1 uk κ ck+1/2 Nx uk +uk+1
i,j − i,j = − α i,j gα m,j m,j
∆t Γ(4 α)∆xα i,m 2
− m=0
X (2.14)
κ dk+1/2 Ney uk +uk+1
+ − β i,j gβ i,m i,m +fk+1/2.
Γ(4 β)∆xβ j,m 2 i,j
− m=0
X
e
6
Similarly, we define
κ ck+1/2 Nx κ ck+1/2 Nx
δ′′ uk = − α i,j gα uk ; δ′′ uk+1 = − α i,j gα uk+1;
α,x i,j Γ(4 α)∆xα i,m m,j α,x i,j Γ(4 α)∆xα i,m m,j
− m=0 − m=0
X X
κ dk+1/2 Ny e κ dk+1/2 Ny e
δ′′ uk = − β i,j gβ uk ; δ′′ uk+1 = − β i,j gβ uk+1;
β,y i,j Γ(4 β)∆yβ j,m i,m β,y i,j Γ(4 β)∆yβ j,m i,m
− m=0 − m=0
X X
then Eq. (2.14) can be rewritten aes e
∆t ∆t ∆t ∆t
1 δ′′ δ′′ uk+1 = 1+ δ′′ + δ′′ uk +∆tfk+1/2. (2.15)
− 2 α,x − 2 β,y i,j 2 α,x 2 β,y i,j i,j
(cid:18) (cid:19) (cid:18) (cid:19)
For the two-dimensional Riesz fractional diffusion equation (2.13), the relevant per-
turbation of (2.15) is of the form
∆t ∆t ∆t ∆t
1 δ′′ 1 δ′′ uk+1 = 1+ δ′′ 1+ δ′′ uk +∆tfk+1/2.
− 2 α,x − 2 β,y i,j 2 α,x 2 β,y i,j i,j
(cid:18) (cid:19)(cid:18) (cid:19) (cid:18) (cid:19)(cid:18) (cid:19)
(2.16)
The scheme (2.16) differs from (2.15) by a perturbation [18]
(∆t)2
δ′′ δ′′ (uk+1 uk ),
4 α,x β,y i,j − i,j
whichmaybededucedbydistributingtheoperatorproductsin(2.16). Since(uk+1 uk )is
i,j − i,j
an (∆t) term, it follows that the perturbation contributes an ((∆t)2) error component
O O
to the truncation error of (2.15). Thus, the scheme (2.16) has a truncation error also
((∆x)2)+ ((∆y)2)+ ((∆t)2).
O O O
For efficiently solving system (2.16), the following techniques can be used:
D-AD scheme [4, 6]:
∆t ∆t
1 δ′′ u∗ = 1+ δ′′ +∆tδ′′ uk +∆tfk+1/2; (2.17)
− 2 α,x i,j 2 α,x β,y i,j i,j
(cid:18) (cid:19) (cid:18) (cid:19)
∆t ∆t
1 δ′′ uk+1 = u∗ δ′′ uk ; (2.18)
− 2 β,y i,j i,j − 2 β,y i,j
(cid:18) (cid:19)
where u∗ is an intermediate solution. Subtracting (2.18) from (2.17), we obtain
i,j
∆t
u∗ = uk+1+ δ′′ uk uk+1 .
i,j i,j 2 β,y i,j − i,j
PR-AD scheme [13]: (cid:0) (cid:1)
∆t ∆t ∆t
1 δ′′ u∗ = 1+ δ′′ uk + fk+1/2; (2.19)
− 2 α,x i,j 2 β,y i,j 2 i,j
(cid:18) (cid:19) (cid:18) (cid:19)
∆t ∆t ∆t
1 δ′′ uk+1 = 1+ δ′′ u∗ + fk+1/2; (2.20)
− 2 β,y i,j 2 α,x i,j 2 i,j
(cid:18) (cid:19) (cid:18) (cid:19)
with intermediate solution u∗ . Subtracting (2.20) from (2.19), we have
i,j
∆t ∆t
2u∗ = 1 δ′′ uk+1 + 1+ δ′′ uk .
i,j − 2 β,y i,j 2 β,y i,j
(cid:18) (cid:19) (cid:18) (cid:19)
7
2.4. D-AD scheme for 3D
Similarly, the resulting discretization of (1.1) can be written as,
∆t ∆t ∆t
1 δ′′ δ′′ δ′′ uk+1
− 2 α,x − 2 β,y − 2 γ,z i,j,l
(cid:18) (cid:19) (2.21)
∆t ∆t ∆t
= 1+ δ′′ + δ′′ + δ′′ uk +fk+1/2∆t.
2 α,x 2 β,y 2 γ,z i,j,l i,j,l
(cid:18) (cid:19)
The perturbation equation of (2.21) is of the form
∆t ∆t ∆t
1 δ′′ 1 δ′′ 1 δ′′ uk+1
− 2 α,x − 2 β,y − 2 γ,z i,j,l
(cid:18) (cid:19)(cid:18) (cid:19)(cid:18) (cid:19) (2.22)
∆t ∆t ∆t
= 1+ δ′′ 1+ δ′′ 1+ δ′′ uk +fk+1/2∆t.
2 α,x 2 β,y 2 γ,z i,j,l i,j,l
(cid:18) (cid:19)(cid:18) (cid:19)(cid:18) (cid:19)
The scheme (2.22) differs from (2.21) by the perturbation term
(∆t)2 (∆t)3
(δ′′ δ′′ +δ′′ δ′′ +δ′′ δ′′ )(uk+1 uk ) δ′′ δ′′ δ′′ (uk+1 uk ).
4 α,x β,y α,x γ,z β,y γ,z i,j,l − i,j,l − 8 α,x β,y γ,z i,j,l − i,j,l
The system of the equations defined by (2.22) can be solved by the D-AD scheme [5, 6]
∆t ∆t
1 δ′′ uk,1 = 1+ δ′′ +∆tδ′′ +∆tδ′′ uk +∆tfk+1/2; (2.23)
− 2 α,x i,j,l 2 α,x β,y γ,z i,j,l i,j,l
(cid:18) (cid:19) (cid:18) (cid:19)
∆t ∆t
1 δ′′ uk,2 = uk,1 δ′′ uk ; (2.24)
− 2 β,y i,j,l i,j,l− 2 β,y i,j,l
(cid:18) (cid:19)
∆t ∆t
1 δ′′ uk+1 = uk,2 δ′′ uk . (2.25)
− 2 γ,z i,j,l i,j,l − 2 γ,z i,j,l
(cid:18) (cid:19)
For maintaining the consistency, we need to carefully specify the boundary conditions
of un,1 and uk,2. According to (2.23)-(2.25), we obtain
i,j,l i,j,l
∆t (∆t)2
uk,1 = uk+1 + δ′′ +δ′′ uk uk+1 + δ′′ δ′′ uk+1;
i,j,l i,j,l 2 β,y γ,z i,j,l− i,j,l 4 β,y γ,z i,j,l
∆t (cid:0) (cid:1)(cid:0) (cid:1)
uk,2 = uk+1 + δ′′ uk uk+1 .
i,j,l i,j,l 2 γ,z i,j,l − i,j,l
(cid:0) (cid:1)
3. Convergence and Stability Analysis
We show the convergence for one-dimensional and multidimensional Riesz fractional
diffusion equationbyproving theconsistency andstability (according toLax’s equivalence
theorem).
Lemma 3.1 ([3]). The coefficients gν , ν (1,2) defined in (2.7) satisfy
i,m ∈
e
(1) gν < 0, gν > 0 (m = i);
i,i i,m 6
Nx Nx
(2) e gν <e0 and gν > gν .
i,m − i,i i,m
m=0 m=0,m6=i
X X
e e e
8
3.1. The stability of the numerical methods in 1D
Theorem 3.2. The Crank-Nicholson scheme (2.11) of the Riesz fractional diffusion
equation (2.9) with 1 < α < 2 is unconditionally stable.
Proof. First, we prove that the eigenvalues of the matrix Ak+1/2 have negative real parts.
Note that Ak+1/2 = gαξk+1/2, and from Lemma 3.1 we obtain
i,i i,i i
Nx Nx
er = Ak+1/2 = ξk+1/2 gα < Ak+1/2. (3.1)
i | i,m | i i,m − i,i
m=0,m6=i m=0,m6=i
X X
e
According to the Greschgorin theorem [8], the eigenvalues of the matrix Ak+1/2 are in the
disks centered at Ak+1/2, with radius r , i.e., the eigenvalues λ of the matrix Ak+1/2 satisfy
i,i i
k+1/2
λ A r , (3.2)
| − i,i | ≤ i
thus, the eigenvalues of the matrix Ak+1/2 have negative real parts. Similarly, we can
prove that the eigenvalues of the matrix I Ak+1/2 have a magnitude greater than 1 and
−
invertible.
Note that λ is an eigenvalue of the matrix Ak+1/2 if and only if 1 λ is an eigenvalue
−
of the matrix I Ak+1/2, if and only if (1 λ)−1(1 + λ) is an eigenvalue of the matrix
− −
(I Ak+1/2)−1(I+Ak+1/2). Since (λ) < 0, it implies that (1 λ)−1(1+λ) < 1. Hence,
− ℜ | − |
the spectral radius of the matrix (I Ak+1/2)−1(I +Ak+1/2) is less than 1.
−
3.2. The stability of the numerical methods in 2D
Under a commutativity assumption for the operators 1 ∆tδ′′ and 1 ∆tδ′′ in
− 2 α,x − 2 β,y
(2.15), the PR-AD scheme and D-AD scheme will be shown to be unconditionally stable.
(cid:0) (cid:1) (cid:0) (cid:1)
The commutativity assumption for these two operators is a common practice in estab-
lishing stability of the classical AD methods for the diffusion [6, 18]. The commutativity
k+1/2 k+1/2
of these operators implies that the matrices A and A given in (3.7) commute.
2D,∆x 2D,∆y
Theorem 3.3. Both the D-AD scheme (2.17)-(2.18) and PR-AD scheme (2.19)-(2.20),
k+1/2
defined by (2.16), are unconditionally stable for α,β (1,2), if the matrices A and
∈ 2D,∆x
k+1/2
A commute.
2D,∆y
Proof. D-AD scheme (2.17)-(2.18) can be expressed in the form
(I Ak+1/2)U∗ = (I +Ak+1/2 +2Ak+1/2)Uk +∆tFk+1/2; (3.3)
− 2D,∆x 2D,∆x 2D,∆y
(I Ak+1/2)Uk+1 = U∗ Ak+1/2Uk; (3.4)
− 2D,∆y − 2D,∆y
and PR-AD scheme (2.19)-(2.20) is of the form
∆t
(I Ak+1/2)U∗ = (I +Ak+1/2)Uk + Fk+1/2; (3.5)
− 2D,∆x 2D,∆y 2
9
∆t
(I Ak+1/2)Uk+1 = (I +Ak+1/2)U∗ + Fk+1/2; (3.6)
− 2D,∆y 2D,∆x 2
where the matrices Ak+1/2 and Ak+1/2 denote the operators ∆tδ′′ and ∆tδ′′ , and
2D,∆x 2D,∆y 2 α,x 2 β,y
Uk = [uk ,uk ,...,uk ,uk ,uk ,...,uk ,...,uk ,uk ,...,uk ]T,
1,1 2,1 Nx−1,1 1,2 2,2 Nx−1,2 1,Ny−1 2,Ny−1 Nx−1,Ny−1
U∗ = [u∗ ,u∗ ,...,u∗ ,u∗ ,u∗ ,...,u∗ ,...,u∗ ,u∗ ,...,u∗ ]T,
1,1 2,1 Nx−1,1 1,2 2,2 Nx−1,2 1,Ny−1 2,Ny−1 Nx−1,Ny−1
and the vector Fk+1/2 absorbs the source terms fk+1/2 and the Dirichlet boundary con-
i,j
k+1/2 k+1/2
ditions at time t = t in the discretized equation. The matrices A and A are
k+1 2D,∆x 2D,∆y
matrices of size (N 1)(N 1) (N 1)(N 1).
x y x y
− − × − −
Let us cancel the intermediate solution U∗, then D-AD scheme and PR-AD scheme
have the same form
(I Ak+1/2)(I Ak+1/2)Uk+1 = (I +Ak+1/2)(I +Ak+1/2)Uk +∆tFk+1/2, (3.7)
− 2D,∆x − 2D,∆y 2D,∆x 2D,∆y
from (3.7) we have the perturbation equation
(I Ak+1/2)(I Ak+1/2)Ek+1 = (I +Ak+1/2)(I +Ak+1/2)Ek,
− 2D,∆x − 2D,∆y 2D,∆x 2D,∆y
where
Ek = [ek ,ek ,...,ek ,ek ,ek ,...,ek ,...,ek ,ek ,...,ek ]T,
1,1 2,1 Nx−1,1 1,2 2,2 Nx−1,2 1,Ny−1 2,Ny−1 Nx−1,Ny−1
and ek = u(x ,y ,t ) uk , consequently
i,j i j k − i,j
Ek = [(I Ak+1/2)−1(I Ak+1/2)−1(I +Ak+1/2)(I +Ak+1/2)]kE0.
− 2D,∆y − 2D,∆x 2D,∆x 2D,∆y
k+1/2 k+1/2
Since the matrices A and A commute, it can be written as
2D,∆x 2D,∆y
Ek = [(I Ak+1/2)−1(I +Ak+1/2)]k[(I Ak+1/2)−1(I +Ak+1/2)]kE0.
− 2D,∆x 2D,∆x − 2D,∆y 2D,∆y
According to Theorem 3.2, similarly, it is easy to check that the eigenvalues of the matrix
k+1/2 k+1/2
A andA havenegativerealparts. Then, boththespectralradiusofthematrixes
2D,∆x 2D,∆y
(I Ak+1/2)−1(I +Ak+1/2) and (I Ak+1/2)−1(I +Ak+1/2) are less than 1, therefore the
− 2D,∆x 2D,∆x − 2D,∆y 2D,∆y
sequence [(I Ak+1/2)−1(I+Ak+1/2)]k and [(I Ak+1/2)−1(I+Ak+1/2)]k converge to zero
− 2D,∆x 2D,∆x − 2D,∆y 2D,∆y
matrix [16]. Hence, the difference scheme (2.16) is unconditionally stable.
3.3. The stability of the numerical methods in 3D
Theorem 3.4. The D-AD scheme (2.23)-(2.25), defined by (2.22), is unconditionally
k+1/2 k+1/2 k+1/2
stable for α,β,γ (1,2), if the matrices A , A and A commute.
∈ 3D,∆x 3D,∆y 3D,∆z
10