# Novel Second-Order Accurate Implicit Numerical Methods for the Riesz Space Distributed-Order Advection-Dispersion Equations

Authors: X. Wang; F. Liu; X. Chen
Journal: Advances in Mathematical Physics
 (2015)

---

## Abstract

We derive and analyze second-order accurate implicit numerical methods for the Riesz space distributed-order advection-dispersion equations (RSDO-ADE) in one-dimensional (1D) and two-dimensional (2D) cases, respectively. Firstly, we discretize the Riesz space distributed-order advection-dispersion equations into multiterm Riesz space fractional advection-dispersion equations (MT-RSDO-ADE) by using the midpoint quadrature rule. Secondly, we propose a second-order accurate implicit numerical method for the MT-RSDO-ADE. Thirdly, stability and convergence are discussed. We investigate the numerical solution and analysis of the RSDO-ADE in 1D case. Then we discuss the RSDO-ADE in 2D case.  For 2D case, we propose a new second-order accurate implicit alternating direction method, and the stability and convergence of this method are proved. Finally, numerical results are presented to support our theoretical analysis.


---

## Body

## 1. Introduction

Fractional differential equations play a significant role in modeling the so-called anomalous transport phenomena and in the theory of complex systems. In recent years there has been a growing interest in the field of fractional calculus. The books [1–4] are completely devoted to different applications of fractional differential equations in many areas, such as engineering, physics, chemistry, astrophysics, and other sciences and historical summaries of the development of fractional calculus. Fractional kinetics systems are widely applied to describe anomalous diffusion or advection-dispersion processes [5, 6]. For processes lacking such scaling the corresponding description may be given by distributed-order fractional partial differential equations [7]. It has been reported that the dynamical systems describing and solving the real world properties have been undergoing two stages. One is from integer-order dynamic systems to fractional-order dynamic systems, and the other is from fractional-order dynamic systems to distributed-order dynamic systems [8].Furthermore, distributed-order differential equations have recently been investigated for complex dynamical systems, namely, distributed-order dynamic systems, which have been explored to describe some important physical phenomena. Distributed-order differential models are more powerful tools to describe complex dynamical systems than classical and fractional-order models because of their nonlocal properties. Chechkin et al. proposed a distributed-order fractional diffusion equation as a generalization of fractional kinetic equations to describe the random process possessing nonunique diffusion exponent and, hence, nonunique Hurst exponent [9]. An important application of distributed-order equations is to model ultraslow diffusion where a plume of particles spreads at a logarithmic rate (see [10, 11]). Kochubei [12] considered the time distributed-order equation and developed a mathematical theory of this equation and studied the derivatives and integrals of distributed order. This equation is applied in physical literature for modeling diffusion with a logarithmic growth of the mean square displacement. As the order of the fractional derivative is distributed over the unit interval, it is useful for modeling a mixture of delay sources (see [13]). Moreover, distributed-order equations may be viewed as consisting of viscoelastic and viscoinertial elements when the order of the fractional derivative varies from zero to two (see [14, 15]). With the motivation of these applications, some attentions have been paid to the fractional partial differential equations (FPDEs) with distributed order [7, 16, 17]. In [18], a series of distributed-order PI controller design methods are derived and applied to the robust control of wheeled service robots, which can tolerate more structural and parametric uncertain ties than the corresponding fractional-order PI control. In the book [8], two initial applications including distributed-order signal processing and optimal distributed damping are provided as motivating examples to further the investigation in the distributed-order dynamic systems.Caputo [19] pointed out that it is very important to investigate the diffusion in porous media for science and engineering field and for social needs especially in the case of water and pollutants. Nowadays the studies of the dissipative and dispersive properties in diffusion equation with fractional order in time and/or space domain in anelastic and dielectric media have been spread to many phenomena from nonlinearity to statistical mechanics and memory formalisms, to represent the diversified forms of deviations from the classic constitutive laws and several complex mathematical methods. Distributed-order equations were first introduced in time domain [20, 21]. Caputo solved the classic problems of anelastic and dielectric media and of diffusion with distributed order in time domain [22]. Then in [19], Caputo considered an extension of the constitutive relation of diffusion to the case when a space memory mechanism operating the medium is represented by a fractional-order differential equations whose order covers a continuum in a given range and introduced distributed order in the constitutive equation in space domain. During the study, Caputo pointed out that the major difference when using space distributed order in the constitutive equation with the case using the single space fractional derivative is that the solutions found in the distributed-order case are potentially more flexible to represent more complex media and more nonlocal phenomena. The difference between the space memory medium and that in time memory, that is, distributed-order equation in space domain and distributed-order equation in time domain, is that the former is more flexible to represent local phenomena while the latter is more reflexible to represent variations in space.There are some numerical methods for distributed-order partial differential equations which have been proposed. Diethelm and Ford [23] introduced and analyzed a numerical method for the solution of a distributed-order differential equations. Meerschaert et al. [13] provided explicit strong solutions and stochastic analogues for distributed-order time-fractional diffusion equations on bounded domains with Dirichlet boundary conditions. Atanackovic et al. [24] studied waves in a viscoelastic rod of finite length. Viscoelastic material is described by a constitutive equation of fractional distributed-order type with the special choice of weight functions. Prescribing boundary conditions on displacement, they obtained displacement and stress in a stress relaxation test. Morgado and Rebelo [25] took into account an implicit scheme for the numerical approximation of the distributed-order time-fractional reaction-diffusion equation with a nonlinear source term.Ye et al. [26] proposed numerical methods for the time distributed-order and Riesz space fractional diffusions and a distributed-order time-fractional diffusion-wave equation, respectively. Hu et al. [27] also considered a new time distributed-order and two-side space-fractional advection-dispersion equation and a time distributed-order diffusion model, respectively. They discretized the distributed-order equation into a multiterm fractional partial differential equation. Some numerical methods for the multiterm fractional partial differential equation have been investigated. Liu et al. [28] considered the multiterm time-fractional wave-diffusion equations and proposed some computationally effective numerical methods for simulating the multiterm time-fractional wave-diffusion equations. Jiang et al. [29] considered the multiterm modified power law wave equations in a finite domain and derived the fundamental solutions of the multiterm modified power law wave equations with the methods and techniques based on Luchko’s Theorem, a spectral representation of the Laplacian operator, a method of separating variables and fractional derivative techniques. Taking the use of the similar methods they derived the analytical solutions of the three types of the space Caputo-Riesz fractional advection-diffusion equations with Dirichlet nonhomogeneous boundary conditions in [30]. Ye et al. [31] derived series expansion based on a spectral representation of the Laplacian operator in a bounded region and gave some applications for the two- and three-dimensional telegraph equation, power law wave equation, and Szabo wave equation. However, published papers on numerical methods of the fractional partial differential equations (FPDEs) with distributed-order especially the space distributed-order FPDES are sparse. This motivates us to consider effective numerical methods for space distributed-order advection-diffusion equations.In this paper, we consider the Riesz space distributed-order advection-diffusion equation (RSDO-ADE) in one-dimensional (1D) and two-dimensional (2D) cases, respectively. The rest of the paper is organized as follows. We discuss numerical method and analysis in 1D case in Sections2 and 3, respectively. We investigate RSDO-ADE in 2D case and propose a new second-order accurate implicit alternating direction method; the stability and convergence of this method are proved in Section 4. Finally, we give two examples to illustrate the behavior of our numerical methods and demonstrate the effectiveness of our theoretical analysis.


## 2. A Second-Order Accurate Implicit Numerical Method for RSDO-ADE in 1D Case

### 2.1. Discretization of the Integral Term

Consider the following Riesz space distributed-order advection-diffusion equation (RSDO-ADE):(1)
∂
u
x
,
t
∂
t
=
∫
0
1
P
α
∂
α
u
x
,
t
∂
x
α
d
α
+
∫
1
2
Q
β
∂
β
u
x
,
t
∂
x
β
d
β
+
f
x
,
t
,where P
(
α
) and Q
(
β
) are nonnegative weight functions which satisfy the conditions(2)
0
≤
P
α
,
P
α
≢
0
,
α
∈
0,1
,
0
<
∫
0
1
P
α
d
α
<
∞
,
0
≤
Q
β
,
Q
β
≢
0
,
β
∈
1,2
,
0
<
∫
1
2
Q
β
d
β
<
∞
,and the Riesz space fractional derivative operators ∂
α
u
(
x
,
t
)
/
∂
x
α and ∂
β
u
(
x
,
t
)
/
∂
x
β on a finite domain [
0
,
L
] are defined as follows:(3)
∂
α
∂
x
α
u
x
,
t
=
-
c
α
D
x
α
0
u
x
,
t
+
D
L
α
x
u
x
,
t
,
∂
β
∂
x
β
u
x
,
t
=
-
c
β
D
x
β
0
u
x
,
t
+
D
L
β
x
u
x
,
t
,where(4)
c
α
=
1
2
cos
⁡
π
α
/
2
,
0
<
α
<
1
,
c
β
=
1
2
cos
⁡
π
β
/
2
,
1
<
β
<
2
,
D
x
α
0
u
x
,
t
=
1
Γ
1
-
α
∂
∂
x
∫
0
x
x
-
ξ
-
α
u
ξ
,
t
d
ξ
,
D
L
α
x
u
x
,
t
=
-
1
Γ
1
-
α
∂
∂
x
∫
x
L
ξ
-
x
-
α
u
ξ
,
t
d
ξ
,
D
x
β
0
u
x
,
t
=
1
Γ
2
-
β
∂
2
∂
x
2
∫
0
x
x
-
ξ
1
-
β
u
ξ
,
t
d
ξ
,
D
L
β
x
u
x
,
t
=
1
Γ
2
-
β
∂
2
∂
x
2
∫
x
L
ξ
-
x
1
-
β
u
ξ
,
t
d
ξ
,where Γ
(
·
) represents the Euler gamma function.Now we consider (1) in the finite domain [
0
,
L
] with the following initial and boundary conditions:(5)
u
0
,
t
=
0
,
u
L
,
t
=
0
,
t
∈
0
,
T
,
(6)
u
x
,
0
=
ϕ
0
x
,
x
∈
0
,
L
.Firstly, we discretize the integral intervals[
0,1
] of α and [
1,2
] of β by the grid 0
=
ξ
0
<
ξ
1
<
⋯
<
ξ
S
=
1, 1
=
η
0
<
η
1
<
⋯
<
η
S
¯
=
2, and denote Δ
ξ
j
=
ξ
j
-
ξ
j
-
1
=
1
/
S
=
σ and α
j
=
ξ
j
+
ξ
j
-
1
/
2
=
2
j
-
1
/
2
S
,
j
=
1,2
,
…
,
S
,
S
∈
N. Consider Δ
η
j
=
η
j
-
η
j
-
1
=
1
/
S
¯
=
ρ and β
j
=
η
j
+
η
j
-
1
/
2
=
1
+
2
j
-
1
/
2
S
¯
,
j
=
1,2
,
…
,
S
¯
,
S
¯
∈
N
.Then by using the midpoint quadrature rule, we obtain(7)
∫
0
1
P
α
∂
α
u
x
,
t
∂
x
α
d
α
=
∑
j
=
1
S
P
α
j
∂
α
j
u
x
,
t
∂
x
α
j
Δ
ξ
j
+
O
σ
2
,
∫
1
2
Q
β
∂
β
u
x
,
t
∂
x
β
d
β
=
∑
j
=
1
S
¯
Q
β
j
∂
β
j
u
x
,
t
∂
x
β
j
Δ
η
j
+
O
ρ
2
.Thus, Riesz space distributed-order advection-diffusion equation (1) in 1D case is now transformed into the following multiterm Riesz space fractional advection-diffusion equation:(8)
∂
u
x
,
t
∂
t
=
∑
j
=
1
S
P
α
j
∂
α
j
u
x
,
t
∂
x
α
j
1
S
+
∑
j
=
1
S
¯
Q
β
j
∂
β
j
u
x
,
t
∂
x
β
j
1
S
¯
+
f
x
,
t
+
O
σ
2
+
ρ
2
.


### 2.2. A Second-Order Accurate Implicit Numerical Method for RSDO-ADE

Then, we discretize the computing domain[
0
,
L
]
×
[
0
,
T
] by x
i
=
i
h
,
i
=
0,1
,
…
,
M, and t
k
=
k
τ
,
k
=
0,1
,
…
,
N, where h
=
L
/
M and τ
=
T
/
N are the space and time steps, respectively, and M and N are two positive integers. Assume that u
(
x
,
t
)
∈
C
5
(
[
a
,
b
]
×
[
0
,
T
]
) and denote(9)
∂
u
x
i
,
t
n
+
1
/
2
∂
t
=
u
x
i
,
t
n
+
1
-
u
x
i
,
t
n
τ
+
O
τ
2
.The Riesz space fractional derivative operators are discretized as follows [32, 33]:(10)
∂
α
j
u
x
i
,
t
n
∂
x
α
j
=
-
1
h
α
j
∑
k
=
1
M
-
1
g
i
-
k
α
j
u
x
k
,
t
n
+
O
h
2
,
0
<
α
j
<
1
,
j
=
1,2
,
…
,
S
,
∂
β
j
u
x
i
,
t
n
∂
x
β
j
=
-
1
h
β
j
∑
k
=
1
M
-
1
g
i
-
k
β
j
u
x
k
,
t
n
+
O
h
2
,
1
<
β
j
≤
2
,
j
=
1,2
,
…
,
S
¯
,where(11)
g
k
α
j
=
-
1
k
Γ
1
+
α
j
Γ
α
j
/
2
-
k
+
1
Γ
α
j
/
2
+
k
+
1
=
1
-
1
+
α
j
α
j
/
2
+
k
g
k
-
1
α
j
,
j
=
1,2
,
…
,
S
,
(12)
g
k
β
j
=
-
1
k
Γ
1
+
β
j
Γ
β
j
/
2
-
k
+
1
Γ
β
j
/
2
+
k
+
1
=
1
-
1
+
β
j
β
j
/
2
+
k
g
k
-
1
β
j
,
j
=
1,2
,
…
,
S
¯
.Lemma 1.
The coefficientsg
k
α
j
(
0
<
α
j
<
1
) of (11) and g
k
β
j
(
1
<
β
j
≤
2
) of (12) satisfy(1)
g
0
α
j
≥
0 and g
0
β
j
≥
0;
(2)
g
-
k
α
j
=
g
k
α
j
≤
0, g
-
k
β
j
=
g
k
β
j
≤
0, for k
=
±
1
,
±
2
,
…;
(3)
∑
k
=
-
∞
∞
g
k
α
j
=
0 and ∑
k
=
-
∞
∞
g
k
β
j
=
0;
(4)
g
0
α
j
=
∑
k
=
-
∞
,
k
≠
0
∞
g
k
α
j and g
0
β
j
=
∑
k
=
-
∞
,
k
≠
0
∞
g
k
β
j.Proof.
(
1
) According to the property of Γ function, we have(13)
g
0
α
j
=
Γ
1
+
α
j
Γ
α
j
/
2
+
1
2
>
0
,
for
α
j
>
0
;
g
0
β
j
=
Γ
1
+
β
j
Γ
β
j
/
2
+
1
2
>
0
,
for
β
j
>
0
.
(
2
) From (11) and (12), we have(14)
g
-
k
α
j
=
-
1
-
k
Γ
1
+
α
j
Γ
α
j
/
2
-
k
+
1
Γ
α
j
/
2
+
k
+
1
=
-
1
k
Γ
1
+
α
j
Γ
α
j
/
2
-
k
+
1
Γ
α
j
/
2
+
k
+
1
=
g
k
α
j
,
g
-
k
β
j
=
-
1
-
k
Γ
1
+
β
j
Γ
β
j
/
2
-
k
+
1
Γ
β
j
/
2
+
k
+
1
=
-
1
k
Γ
1
+
β
j
Γ
β
j
/
2
-
k
+
1
Γ
β
j
/
2
+
k
+
1
=
g
k
β
j
.For k
=
1, we have(15)
g
1
α
j
=
1
-
1
+
α
j
α
j
/
2
+
1
g
0
α
j
≤
0
,
g
1
β
j
=
1
-
1
+
β
j
β
j
/
2
+
1
g
0
β
j
≤
0
.Assuming that g
k
α
j
≤
0, g
k
β
j
≤
0. Since 0
<
α
j
<
1 and 1
<
β
j
≤
2, we have(16)
1
-
1
+
α
j
α
j
/
2
+
k
+
1
>
0
,
1
-
1
+
β
j
β
j
/
2
+
k
+
1
>
0
.Therefore,(17)
g
k
+
1
α
j
=
1
-
1
+
α
j
α
j
/
2
+
k
+
1
g
k
α
j
≤
0
,
j
=
1,2
,
…
,
s
,
g
k
+
1
β
j
=
1
-
1
+
β
j
β
j
/
2
+
k
+
1
g
k
β
j
≤
0
,
j
=
1,2
,
…
,
s
¯
.
By mathematical induction method, conclusion(
2
) is proved.
(
3
) From the following formula [32],(18)
2
sin
⁡
z
2
α
j
=
∑
k
=
-
∞
+
∞
-
1
k
Γ
α
j
+
1
Γ
α
j
/
2
-
k
+
1
Γ
α
j
/
2
+
k
+
1
e
i
k
z
,we have ∑
k
=
-
∞
+
∞
(
-
1
)
k
Γ
(
α
j
+
1
)
/
Γ
(
α
j
/
2
-
k
+
1
)
Γ
(
α
j
/
2
+
k
+
1
)
=
0; namely, ∑
k
=
-
∞
∞
g
k
α
j
=
0
.
The conclusion that∑
k
=
-
∞
∞
g
k
β
j
=
0 is also obtained by the same method.
(
4
) According to (
3
), we have g
0
α
j
+
∑
k
=
-
∞
,
k
≠
0
+
∞
g
k
α
j
=
0, and using the conclusion g
k
α
j
=
g
-
k
α
j
≤
0, we get ∑
k
=
-
∞
,
k
≠
0
+
∞
g
k
α
j
=
g
0
α
j.
The conclusion that∑
k
=
-
∞
,
k
≠
0
+
∞
g
k
β
j
=
g
0
β
j is also obtained by the same method.Using Crank-Nicholson method and second-order accurate implicit finite difference scheme, we obtain the following discrete form for RSDO-ADE in 1D case:(19)
u
x
i
,
t
n
+
1
-
u
x
i
,
t
n
τ
=
-
1
2
∑
j
=
1
S
P
α
j
S
1
h
α
j
∑
k
=
1
M
-
1
g
i
-
k
α
j
u
x
k
,
t
n
-
1
2
∑
j
=
1
S
P
α
j
S
1
h
α
j
∑
k
=
1
M
-
1
g
i
-
k
α
j
u
x
i
,
t
n
+
1
-
1
2
∑
j
=
1
S
¯
Q
β
j
S
¯
1
h
β
j
∑
k
=
1
M
-
1
g
i
-
k
β
j
u
x
k
,
t
n
-
1
2
∑
j
=
1
S
¯
Q
β
j
S
¯
1
h
β
j
∑
k
=
1
M
-
1
g
i
-
k
β
j
u
x
i
,
t
n
+
1
+
f
x
i
,
t
n
+
1
/
2
+
R
i
n
+
1
/
2
,where the local truncation error R
i
n
+
1
/
2
=
O
(
σ
2
+
ρ
2
+
τ
2
+
h
2
) and f
(
x
i
,
t
n
+
1
/
2
)
=
1
/
2
(
f
(
x
k
,
t
n
+
1
)
+
f
(
x
k
,
t
n
)
).By omitting the local truncation error termR
i
n
+
1
/
2 in (19). We obtain the following second-order accurate implicit numerical method for RSDO-ADE in 1D case:(20)
u
i
n
+
1
-
u
i
n
τ
=
-
1
2
∑
j
=
1
S
P
α
j
S
1
h
α
j
∑
k
=
1
M
-
1
g
i
-
k
α
j
u
k
n
-
1
2
∑
j
=
1
S
P
α
j
S
1
h
α
j
∑
k
=
1
M
-
1
g
i
-
k
α
j
u
k
n
+
1
-
1
2
∑
j
=
1
S
¯
Q
β
j
S
¯
1
h
β
j
∑
k
=
1
M
-
1
g
i
-
k
β
j
u
k
n
-
1
2
∑
j
=
1
S
¯
Q
β
j
S
¯
1
h
β
j
∑
k
=
1
M
-
1
g
i
-
k
β
j
u
k
n
+
1
+
f
i
n
+
1
/
2
,
i
=
1,2
,
…
,
M
-
1
;
n
=
0,1
,
…
,
N
-
1
;
(21)
u
0
k
=
0
,
u
M
k
=
0
,
k
=
1,2
,
…
,
N
,
(22)
u
i
0
=
ϕ
0
x
i
,
i
=
0,1
,
…
,
M
.We define the function space as follows:Λ
(
Ω
)
=
{
u
(
x
,
t
)
∣
∂
5
u
x
,
t
/
∂
x
5
,
∂
4
u
x
,
t
/
∂
x
2
∂
t
2
∈
C
(
Ω
)
}, where Ω
=
[
0
,
L
]
×
[
0
,
T
]. In this paper, we suppose that the problem: (1) satisfies conditions (5) and (6) has a smooth solution u
x
,
t
∈
Λ
(
Ω
), and f
(
x
,
t
), ϕ
0
(
x
) are sufficiently smooth functions.


## 2.1. Discretization of the Integral Term

Consider the following Riesz space distributed-order advection-diffusion equation (RSDO-ADE):(1)
∂
u
x
,
t
∂
t
=
∫
0
1
P
α
∂
α
u
x
,
t
∂
x
α
d
α
+
∫
1
2
Q
β
∂
β
u
x
,
t
∂
x
β
d
β
+
f
x
,
t
,where P
(
α
) and Q
(
β
) are nonnegative weight functions which satisfy the conditions(2)
0
≤
P
α
,
P
α
≢
0
,
α
∈
0,1
,
0
<
∫
0
1
P
α
d
α
<
∞
,
0
≤
Q
β
,
Q
β
≢
0
,
β
∈
1,2
,
0
<
∫
1
2
Q
β
d
β
<
∞
,and the Riesz space fractional derivative operators ∂
α
u
(
x
,
t
)
/
∂
x
α and ∂
β
u
(
x
,
t
)
/
∂
x
β on a finite domain [
0
,
L
] are defined as follows:(3)
∂
α
∂
x
α
u
x
,
t
=
-
c
α
D
x
α
0
u
x
,
t
+
D
L
α
x
u
x
,
t
,
∂
β
∂
x
β
u
x
,
t
=
-
c
β
D
x
β
0
u
x
,
t
+
D
L
β
x
u
x
,
t
,where(4)
c
α
=
1
2
cos
⁡
π
α
/
2
,
0
<
α
<
1
,
c
β
=
1
2
cos
⁡
π
β
/
2
,
1
<
β
<
2
,
D
x
α
0
u
x
,
t
=
1
Γ
1
-
α
∂
∂
x
∫
0
x
x
-
ξ
-
α
u
ξ
,
t
d
ξ
,
D
L
α
x
u
x
,
t
=
-
1
Γ
1
-
α
∂
∂
x
∫
x
L
ξ
-
x
-
α
u
ξ
,
t
d
ξ
,
D
x
β
0
u
x
,
t
=
1
Γ
2
-
β
∂
2
∂
x
2
∫
0
x
x
-
ξ
1
-
β
u
ξ
,
t
d
ξ
,
D
L
β
x
u
x
,
t
=
1
Γ
2
-
β
∂
2
∂
x
2
∫
x
L
ξ
-
x
1
-
β
u
ξ
,
t
d
ξ
,where Γ
(
·
) represents the Euler gamma function.Now we consider (1) in the finite domain [
0
,
L
] with the following initial and boundary conditions:(5)
u
0
,
t
=
0
,
u
L
,
t
=
0
,
t
∈
0
,
T
,
(6)
u
x
,
0
=
ϕ
0
x
,
x
∈
0
,
L
.Firstly, we discretize the integral intervals[
0,1
] of α and [
1,2
] of β by the grid 0
=
ξ
0
<
ξ
1
<
⋯
<
ξ
S
=
1, 1
=
η
0
<
η
1
<
⋯
<
η
S
¯
=
2, and denote Δ
ξ
j
=
ξ
j
-
ξ
j
-
1
=
1
/
S
=
σ and α
j
=
ξ
j
+
ξ
j
-
1
/
2
=
2
j
-
1
/
2
S
,
j
=
1,2
,
…
,
S
,
S
∈
N. Consider Δ
η
j
=
η
j
-
η
j
-
1
=
1
/
S
¯
=
ρ and β
j
=
η
j
+
η
j
-
1
/
2
=
1
+
2
j
-
1
/
2
S
¯
,
j
=
1,2
,
…
,
S
¯
,
S
¯
∈
N
.Then by using the midpoint quadrature rule, we obtain(7)
∫
0
1
P
α
∂
α
u
x
,
t
∂
x
α
d
α
=
∑
j
=
1
S
P
α
j
∂
α
j
u
x
,
t
∂
x
α
j
Δ
ξ
j
+
O
σ
2
,
∫
1
2
Q
β
∂
β
u
x
,
t
∂
x
β
d
β
=
∑
j
=
1
S
¯
Q
β
j
∂
β
j
u
x
,
t
∂
x
β
j
Δ
η
j
+
O
ρ
2
.Thus, Riesz space distributed-order advection-diffusion equation (1) in 1D case is now transformed into the following multiterm Riesz space fractional advection-diffusion equation:(8)
∂
u
x
,
t
∂
t
=
∑
j
=
1
S
P
α
j
∂
α
j
u
x
,
t
∂
x
α
j
1
S
+
∑
j
=
1
S
¯
Q
β
j
∂
β
j
u
x
,
t
∂
x
β
j
1
S
¯
+
f
x
,
t
+
O
σ
2
+
ρ
2
.


## 2.2. A Second-Order Accurate Implicit Numerical Method for RSDO-ADE

Then, we discretize the computing domain[
0
,
L
]
×
[
0
,
T
] by x
i
=
i
h
,
i
=
0,1
,
…
,
M, and t
k
=
k
τ
,
k
=
0,1
,
…
,
N, where h
=
L
/
M and τ
=
T
/
N are the space and time steps, respectively, and M and N are two positive integers. Assume that u
(
x
,
t
)
∈
C
5
(
[
a
,
b
]
×
[
0
,
T
]
) and denote(9)
∂
u
x
i
,
t
n
+
1
/
2
∂
t
=
u
x
i
,
t
n
+
1
-
u
x
i
,
t
n
τ
+
O
τ
2
.The Riesz space fractional derivative operators are discretized as follows [32, 33]:(10)
∂
α
j
u
x
i
,
t
n
∂
x
α
j
=
-
1
h
α
j
∑
k
=
1
M
-
1
g
i
-
k
α
j
u
x
k
,
t
n
+
O
h
2
,
0
<
α
j
<
1
,
j
=
1,2
,
…
,
S
,
∂
β
j
u
x
i
,
t
n
∂
x
β
j
=
-
1
h
β
j
∑
k
=
1
M
-
1
g
i
-
k
β
j
u
x
k
,
t
n
+
O
h
2
,
1
<
β
j
≤
2
,
j
=
1,2
,
…
,
S
¯
,where(11)
g
k
α
j
=
-
1
k
Γ
1
+
α
j
Γ
α
j
/
2
-
k
+
1
Γ
α
j
/
2
+
k
+
1
=
1
-
1
+
α
j
α
j
/
2
+
k
g
k
-
1
α
j
,
j
=
1,2
,
…
,
S
,
(12)
g
k
β
j
=
-
1
k
Γ
1
+
β
j
Γ
β
j
/
2
-
k
+
1
Γ
β
j
/
2
+
k
+
1
=
1
-
1
+
β
j
β
j
/
2
+
k
g
k
-
1
β
j
,
j
=
1,2
,
…
,
S
¯
.Lemma 1.
The coefficientsg
k
α
j
(
0
<
α
j
<
1
) of (11) and g
k
β
j
(
1
<
β
j
≤
2
) of (12) satisfy(1)
g
0
α
j
≥
0 and g
0
β
j
≥
0;
(2)
g
-
k
α
j
=
g
k
α
j
≤
0, g
-
k
β
j
=
g
k
β
j
≤
0, for k
=
±
1
,
±
2
,
…;
(3)
∑
k
=
-
∞
∞
g
k
α
j
=
0 and ∑
k
=
-
∞
∞
g
k
β
j
=
0;
(4)
g
0
α
j
=
∑
k
=
-
∞
,
k
≠
0
∞
g
k
α
j and g
0
β
j
=
∑
k
=
-
∞
,
k
≠
0
∞
g
k
β
j.Proof.
(
1
) According to the property of Γ function, we have(13)
g
0
α
j
=
Γ
1
+
α
j
Γ
α
j
/
2
+
1
2
>
0
,
for
α
j
>
0
;
g
0
β
j
=
Γ
1
+
β
j
Γ
β
j
/
2
+
1
2
>
0
,
for
β
j
>
0
.
(
2
) From (11) and (12), we have(14)
g
-
k
α
j
=
-
1
-
k
Γ
1
+
α
j
Γ
α
j
/
2
-
k
+
1
Γ
α
j
/
2
+
k
+
1
=
-
1
k
Γ
1
+
α
j
Γ
α
j
/
2
-
k
+
1
Γ
α
j
/
2
+
k
+
1
=
g
k
α
j
,
g
-
k
β
j
=
-
1
-
k
Γ
1
+
β
j
Γ
β
j
/
2
-
k
+
1
Γ
β
j
/
2
+
k
+
1
=
-
1
k
Γ
1
+
β
j
Γ
β
j
/
2
-
k
+
1
Γ
β
j
/
2
+
k
+
1
=
g
k
β
j
.For k
=
1, we have(15)
g
1
α
j
=
1
-
1
+
α
j
α
j
/
2
+
1
g
0
α
j
≤
0
,
g
1
β
j
=
1
-
1
+
β
j
β
j
/
2
+
1
g
0
β
j
≤
0
.Assuming that g
k
α
j
≤
0, g
k
β
j
≤
0. Since 0
<
α
j
<
1 and 1
<
β
j
≤
2, we have(16)
1
-
1
+
α
j
α
j
/
2
+
k
+
1
>
0
,
1
-
1
+
β
j
β
j
/
2
+
k
+
1
>
0
.Therefore,(17)
g
k
+
1
α
j
=
1
-
1
+
α
j
α
j
/
2
+
k
+
1
g
k
α
j
≤
0
,
j
=
1,2
,
…
,
s
,
g
k
+
1
β
j
=
1
-
1
+
β
j
β
j
/
2
+
k
+
1
g
k
β
j
≤
0
,
j
=
1,2
,
…
,
s
¯
.
By mathematical induction method, conclusion(
2
) is proved.
(
3
) From the following formula [32],(18)
2
sin
⁡
z
2
α
j
=
∑
k
=
-
∞
+
∞
-
1
k
Γ
α
j
+
1
Γ
α
j
/
2
-
k
+
1
Γ
α
j
/
2
+
k
+
1
e
i
k
z
,we have ∑
k
=
-
∞
+
∞
(
-
1
)
k
Γ
(
α
j
+
1
)
/
Γ
(
α
j
/
2
-
k
+
1
)
Γ
(
α
j
/
2
+
k
+
1
)
=
0; namely, ∑
k
=
-
∞
∞
g
k
α
j
=
0
.
The conclusion that∑
k
=
-
∞
∞
g
k
β
j
=
0 is also obtained by the same method.
(
4
) According to (
3
), we have g
0
α
j
+
∑
k
=
-
∞
,
k
≠
0
+
∞
g
k
α
j
=
0, and using the conclusion g
k
α
j
=
g
-
k
α
j
≤
0, we get ∑
k
=
-
∞
,
k
≠
0
+
∞
g
k
α
j
=
g
0
α
j.
The conclusion that∑
k
=
-
∞
,
k
≠
0
+
∞
g
k
β
j
=
g
0
β
j is also obtained by the same method.Using Crank-Nicholson method and second-order accurate implicit finite difference scheme, we obtain the following discrete form for RSDO-ADE in 1D case:(19)
u
x
i
,
t
n
+
1
-
u
x
i
,
t
n
τ
=
-
1
2
∑
j
=
1
S
P
α
j
S
1
h
α
j
∑
k
=
1
M
-
1
g
i
-
k
α
j
u
x
k
,
t
n
-
1
2
∑
j
=
1
S
P
α
j
S
1
h
α
j
∑
k
=
1
M
-
1
g
i
-
k
α
j
u
x
i
,
t
n
+
1
-
1
2
∑
j
=
1
S
¯
Q
β
j
S
¯
1
h
β
j
∑
k
=
1
M
-
1
g
i
-
k
β
j
u
x
k
,
t
n
-
1
2
∑
j
=
1
S
¯
Q
β
j
S
¯
1
h
β
j
∑
k
=
1
M
-
1
g
i
-
k
β
j
u
x
i
,
t
n
+
1
+
f
x
i
,
t
n
+
1
/
2
+
R
i
n
+
1
/
2
,where the local truncation error R
i
n
+
1
/
2
=
O
(
σ
2
+
ρ
2
+
τ
2
+
h
2
) and f
(
x
i
,
t
n
+
1
/
2
)
=
1
/
2
(
f
(
x
k
,
t
n
+
1
)
+
f
(
x
k
,
t
n
)
).By omitting the local truncation error termR
i
n
+
1
/
2 in (19). We obtain the following second-order accurate implicit numerical method for RSDO-ADE in 1D case:(20)
u
i
n
+
1
-
u
i
n
τ
=
-
1
2
∑
j
=
1
S
P
α
j
S
1
h
α
j
∑
k
=
1
M
-
1
g
i
-
k
α
j
u
k
n
-
1
2
∑
j
=
1
S
P
α
j
S
1
h
α
j
∑
k
=
1
M
-
1
g
i
-
k
α
j
u
k
n
+
1
-
1
2
∑
j
=
1
S
¯
Q
β
j
S
¯
1
h
β
j
∑
k
=
1
M
-
1
g
i
-
k
β
j
u
k
n
-
1
2
∑
j
=
1
S
¯
Q
β
j
S
¯
1
h
β
j
∑
k
=
1
M
-
1
g
i
-
k
β
j
u
k
n
+
1
+
f
i
n
+
1
/
2
,
i
=
1,2
,
…
,
M
-
1
;
n
=
0,1
,
…
,
N
-
1
;
(21)
u
0
k
=
0
,
u
M
k
=
0
,
k
=
1,2
,
…
,
N
,
(22)
u
i
0
=
ϕ
0
x
i
,
i
=
0,1
,
…
,
M
.We define the function space as follows:Λ
(
Ω
)
=
{
u
(
x
,
t
)
∣
∂
5
u
x
,
t
/
∂
x
5
,
∂
4
u
x
,
t
/
∂
x
2
∂
t
2
∈
C
(
Ω
)
}, where Ω
=
[
0
,
L
]
×
[
0
,
T
]. In this paper, we suppose that the problem: (1) satisfies conditions (5) and (6) has a smooth solution u
x
,
t
∈
Λ
(
Ω
), and f
(
x
,
t
), ϕ
0
(
x
) are sufficiently smooth functions.


## 3. Numerical Analysis of the Second-Order Accurate Implicit Numerical Method in 1D Case

In this subsection, we discuss the stability and convergence of the second-order accurate implicit numerical method (20)–(22).Equation (20) can be rewritten as(23)
u
i
n
+
1
+
τ
2
∑
j
=
1
S
P
α
j
S
1
h
α
j
∑
k
=
1
M
-
1
g
i
-
k
α
j
u
k
n
+
1
+
τ
2
∑
j
=
1
S
¯
Q
β
j
S
¯
1
h
β
j
∑
k
=
1
M
-
1
g
i
-
k
β
j
u
k
n
+
1
=
u
i
n
-
τ
2
∑
j
=
1
S
P
α
j
S
1
h
α
j
∑
k
=
1
M
-
1
g
i
-
k
α
j
u
k
n
-
τ
2
∑
j
=
1
S
¯
Q
β
j
S
¯
1
h
β
j
∑
k
=
1
M
-
1
g
i
-
k
β
j
u
k
n
+
τ
f
i
n
+
1
/
2
,
i
=
1,2
,
…
,
M
-
1
;
n
=
0,1
,
…
,
N
-
1
.Further, (23) can be written into the following matrix form:(24)
I
+
A
U
n
+
1
=
I
-
A
U
n
+
τ
F
n
+
1
/
2
,where A
=
A
1
+
A
2
=
[
a
i
,
j
]
(
M
-
1
)
×
(
M
-
1
). Consider(25)
A
1
=
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
0
α
j
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
-
1
α
j
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
-
2
α
j
⋯
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
-
M
+
2
α
j
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
1
α
j
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
0
α
j
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
-
1
α
j
⋯
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
-
M
+
3
α
j
⋮
⋮
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
M
-
2
α
j
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
M
-
3
α
j
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
M
-
4
α
j
⋯
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
0
α
j
,
A
2
=
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
0
β
j
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
-
1
β
j
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
-
2
β
j
⋯
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
-
M
+
2
β
j
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
1
β
j
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
0
β
j
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
-
1
β
j
⋯
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
-
M
+
3
β
j
⋮
⋮
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
M
-
2
β
j
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
M
-
3
β
j
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
M
-
4
β
j
⋯
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
0
β
j
.Consider U
n
=
(
u
1
n
,
u
2
n
,
…
,
u
M
-
1
n
)
T and F
n
+
1
/
2
=
(
f
1
n
+
1
/
2
,
f
2
n
+
1
/
2
,
…
,
f
M
-
1
n
+
1
/
2
)
T and I is a (
M
-
1
)
×
(
M
-
1
) identity matrix.According to Lemma1, we have the following Lemma.Lemma 2.
MatrixA is symmetric and strictly diagonally dominant.


### 3.1. Stability of the Second-Order Accurate Implicit Numerical Method

Lemma 3 (see [33]).
If matrixI
+
A is invertible, then matrices (
I
+
A
)
-
1 and I
-
A commute.Lemma 4 (see [33]).
If bothA and B are symmetric matrices of order n and matrices A and B commute, then A
B is symmetric.Lemma 5 (see [33]).
If matrixA is real and symmetric, then A
2
=
[
ρ
(
A
2
)
]
1
/
2
=
[
ρ
2
(
A
)
]
1
/
2
=
ρ
(
A
)
=
max
⁡
λ
i.Theorem 6.
The second-order accurate implicit numerical method (20)–(22) for RSDO-ADE (1), (5), and (6) is unconditionally stable.Proof.
Assuming thatu
i
n and u
~
i
n are numerical solution and approximation solution of the second-order accurate implicit numerical method (20)–(22), let ε
i
k
=
u
~
i
k
-
u
i
k and Y
n
=
ε
1
n
,
ε
2
n
,
…
,
ε
M
-
1
n
T. Then the error Y
n satisfies the following equation:(26)
Y
n
+
1
=
I
+
A
-
1
I
-
A
Y
n
.
According to the Gershgorin theorem the eigenvalues at each diagonal entry are(27)
a
i
i
=
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
0
α
j
+
τ
2
S
¯
∑
j
=
1
S
¯
P
β
j
h
β
j
g
0
β
j
,with radius r
i
=
∑
k
=
1
,
k
≠
i
M
-
1
τ
/
2
S
∑
j
=
1
S
P
α
j
/
h
α
j
g
i
-
k
α
j
+
τ
/
2
S
¯
∑
j
=
1
S
¯
P
β
j
/
h
β
j
g
i
-
k
β
j
. Using Lemma 1, we obtain(28)
λ
i
-
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
0
α
j
-
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
0
β
j
≤
∑
k
=
1
,
k
≠
i
M
-
1
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
i
-
k
α
j
+
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
i
-
k
β
j
≤
∑
k
=
1
,
k
≠
i
M
-
1
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
i
-
k
α
j
+
∑
k
=
1
,
k
≠
i
M
-
1
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
i
-
k
β
j
≤
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
0
α
j
+
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
0
β
j
.
Further, according to (28), we have(29)
0
<
λ
<
2
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
0
α
j
+
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
0
β
j
.Next, λ is an eigenvalue of A, then 1
+
λ is an eigenvalue of matrix I
+
A, and 1
-
λ
/
1
+
λ is an eigenvalue of matrix (
I
+
A
)
-
1
(
I
+
A
)
. According to (29), we have λ
>
0, and it is easy to check that(30)
1
-
λ
1
+
λ
<
1
.Thus, we obtain the following conclusion:(31)
ρ
I
+
A
-
1
I
-
A
<
1
.
It is obviously thatI and A are both symmetric. Therefore matrices I
-
A, I
+
A, and (
I
+
A
)
-
1 are also symmetric. Since I
+
A is invertible, according to Lemma 3,  (
I
+
A
)
-
1 and I
-
A commute. Using Lemma 4, W
=
(
I
+
A
)
-
1
(
I
-
A
) is symmetric. Based on Lemma 5, we have W
2
=
ρ
(
W
)
<
1.
Therefore, we obtain(32)
Y
n
+
1
2
≤
W
2
Y
n
2
≤
⋯
≤
Y
0
2
.This completes the proof.


### 3.2. Convergence of the Second-Order Accurate Implicit Numerical Method

Now let us consider the convergence of second-order accurate implicit numerical method (20)–(22).Theorem 7.
Assuming that RSDO-ADE (1), (5), and (6) have smooth solution u
(
x
i
,
t
n
)
∈
Λ
(
Ω
), {
u
i
n
} is the solution of second-order accurate implicit numerical method (20)–(22). Let error e
i
n
=
u
(
x
i
,
t
n
)
-
u
i
n and E
n
=
(
e
1
n
,
e
2
n
,
…
,
e
M
-
1
n
)
T. Then there exists a constant C such that(33)
E
n
2
≤
C
σ
2
+
ρ
2
+
τ
2
+
h
2
,
1
≤
n
≤
N
.Proof.
From (19) and (23), we obtain the following error equations:(34)
e
i
n
+
1
-
e
i
n
τ
=
-
1
2
∑
j
=
1
S
P
α
j
S
1
h
α
j
∑
k
=
1
M
-
1
g
i
-
k
α
j
e
k
n
-
1
2
∑
j
=
1
S
P
α
j
S
1
h
α
j
∑
k
=
1
M
-
1
g
i
-
k
α
j
e
k
n
+
1
-
1
2
∑
j
=
1
S
¯
Q
β
j
S
¯
1
h
β
j
∑
k
=
1
M
-
1
g
i
-
k
β
j
e
k
n
-
1
2
∑
j
=
1
S
¯
Q
β
j
S
¯
1
h
β
j
∑
k
=
1
M
-
1
g
i
-
k
β
j
e
k
n
+
1
+
R
i
n
+
1
/
2
,
(35)
e
0
k
=
0
,
e
M
k
=
0
,
k
=
1,2
,
…
,
N
,
(36)
e
i
0
=
0
,
i
=
0,1
,
…
,
M
.We rewrite (34) into the following matrix form:(37)
I
+
A
E
n
+
1
=
I
-
A
E
n
+
τ
C
1
ρ
2
+
σ
2
+
τ
2
+
h
2
I
,where I
=
I
(
M
-
1
)
×
(
M
-
1
) is unix matrix.
Then we have(38)
E
n
+
1
=
I
+
A
-
1
I
-
A
E
n
+
C
1
τ
I
+
A
-
1
ρ
2
+
σ
2
+
τ
2
+
h
2
.Further, we obtain(39)
E
n
+
1
2
=
I
+
A
-
1
I
-
A
E
n
+
C
1
τ
I
+
A
-
1
σ
2
+
ρ
2
+
τ
2
+
h
2
2
≤
I
+
A
-
1
I
-
A
E
n
2
+
C
1
τ
I
+
A
-
1
2
σ
2
+
ρ
2
+
τ
2
+
h
2
≤
E
n
2
+
C
1
τ
σ
2
+
ρ
2
+
τ
2
+
h
2
≤
E
0
2
+
C
1
n
τ
σ
2
+
ρ
2
+
τ
2
+
h
2
≤
C
σ
2
+
ρ
2
+
τ
2
+
h
2
.
This completes the proof.


## 3.1. Stability of the Second-Order Accurate Implicit Numerical Method

Lemma 3 (see [33]).
If matrixI
+
A is invertible, then matrices (
I
+
A
)
-
1 and I
-
A commute.Lemma 4 (see [33]).
If bothA and B are symmetric matrices of order n and matrices A and B commute, then A
B is symmetric.Lemma 5 (see [33]).
If matrixA is real and symmetric, then A
2
=
[
ρ
(
A
2
)
]
1
/
2
=
[
ρ
2
(
A
)
]
1
/
2
=
ρ
(
A
)
=
max
⁡
λ
i.Theorem 6.
The second-order accurate implicit numerical method (20)–(22) for RSDO-ADE (1), (5), and (6) is unconditionally stable.Proof.
Assuming thatu
i
n and u
~
i
n are numerical solution and approximation solution of the second-order accurate implicit numerical method (20)–(22), let ε
i
k
=
u
~
i
k
-
u
i
k and Y
n
=
ε
1
n
,
ε
2
n
,
…
,
ε
M
-
1
n
T. Then the error Y
n satisfies the following equation:(26)
Y
n
+
1
=
I
+
A
-
1
I
-
A
Y
n
.
According to the Gershgorin theorem the eigenvalues at each diagonal entry are(27)
a
i
i
=
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
0
α
j
+
τ
2
S
¯
∑
j
=
1
S
¯
P
β
j
h
β
j
g
0
β
j
,with radius r
i
=
∑
k
=
1
,
k
≠
i
M
-
1
τ
/
2
S
∑
j
=
1
S
P
α
j
/
h
α
j
g
i
-
k
α
j
+
τ
/
2
S
¯
∑
j
=
1
S
¯
P
β
j
/
h
β
j
g
i
-
k
β
j
. Using Lemma 1, we obtain(28)
λ
i
-
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
0
α
j
-
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
0
β
j
≤
∑
k
=
1
,
k
≠
i
M
-
1
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
i
-
k
α
j
+
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
i
-
k
β
j
≤
∑
k
=
1
,
k
≠
i
M
-
1
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
i
-
k
α
j
+
∑
k
=
1
,
k
≠
i
M
-
1
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
i
-
k
β
j
≤
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
0
α
j
+
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
0
β
j
.
Further, according to (28), we have(29)
0
<
λ
<
2
τ
2
S
∑
j
=
1
S
P
α
j
h
α
j
g
0
α
j
+
τ
2
S
¯
∑
j
=
1
S
¯
Q
β
j
h
β
j
g
0
β
j
.Next, λ is an eigenvalue of A, then 1
+
λ is an eigenvalue of matrix I
+
A, and 1
-
λ
/
1
+
λ is an eigenvalue of matrix (
I
+
A
)
-
1
(
I
+
A
)
. According to (29), we have λ
>
0, and it is easy to check that(30)
1
-
λ
1
+
λ
<
1
.Thus, we obtain the following conclusion:(31)
ρ
I
+
A
-
1
I
-
A
<
1
.
It is obviously thatI and A are both symmetric. Therefore matrices I
-
A, I
+
A, and (
I
+
A
)
-
1 are also symmetric. Since I
+
A is invertible, according to Lemma 3,  (
I
+
A
)
-
1 and I
-
A commute. Using Lemma 4, W
=
(
I
+
A
)
-
1
(
I
-
A
) is symmetric. Based on Lemma 5, we have W
2
=
ρ
(
W
)
<
1.
Therefore, we obtain(32)
Y
n
+
1
2
≤
W
2
Y
n
2
≤
⋯
≤
Y
0
2
.This completes the proof.


## 3.2. Convergence of the Second-Order Accurate Implicit Numerical Method

Now let us consider the convergence of second-order accurate implicit numerical method (20)–(22).Theorem 7.
Assuming that RSDO-ADE (1), (5), and (6) have smooth solution u
(
x
i
,
t
n
)
∈
Λ
(
Ω
), {
u
i
n
} is the solution of second-order accurate implicit numerical method (20)–(22). Let error e
i
n
=
u
(
x
i
,
t
n
)
-
u
i
n and E
n
=
(
e
1
n
,
e
2
n
,
…
,
e
M
-
1
n
)
T. Then there exists a constant C such that(33)
E
n
2
≤
C
σ
2
+
ρ
2
+
τ
2
+
h
2
,
1
≤
n
≤
N
.Proof.
From (19) and (23), we obtain the following error equations:(34)
e
i
n
+
1
-
e
i
n
τ
=
-
1
2
∑
j
=
1
S
P
α
j
S
1
h
α
j
∑
k
=
1
M
-
1
g
i
-
k
α
j
e
k
n
-
1
2
∑
j
=
1
S
P
α
j
S
1
h
α
j
∑
k
=
1
M
-
1
g
i
-
k
α
j
e
k
n
+
1
-
1
2
∑
j
=
1
S
¯
Q
β
j
S
¯
1
h
β
j
∑
k
=
1
M
-
1
g
i
-
k
β
j
e
k
n
-
1
2
∑
j
=
1
S
¯
Q
β
j
S
¯
1
h
β
j
∑
k
=
1
M
-
1
g
i
-
k
β
j
e
k
n
+
1
+
R
i
n
+
1
/
2
,
(35)
e
0
k
=
0
,
e
M
k
=
0
,
k
=
1,2
,
…
,
N
,
(36)
e
i
0
=
0
,
i
=
0,1
,
…
,
M
.We rewrite (34) into the following matrix form:(37)
I
+
A
E
n
+
1
=
I
-
A
E
n
+
τ
C
1
ρ
2
+
σ
2
+
τ
2
+
h
2
I
,where I
=
I
(
M
-
1
)
×
(
M
-
1
) is unix matrix.
Then we have(38)
E
n
+
1
=
I
+
A
-
1
I
-
A
E
n
+
C
1
τ
I
+
A
-
1
ρ
2
+
σ
2
+
τ
2
+
h
2
.Further, we obtain(39)
E
n
+
1
2
=
I
+
A
-
1
I
-
A
E
n
+
C
1
τ
I
+
A
-
1
σ
2
+
ρ
2
+
τ
2
+
h
2
2
≤
I
+
A
-
1
I
-
A
E
n
2
+
C
1
τ
I
+
A
-
1
2
σ
2
+
ρ
2
+
τ
2
+
h
2
≤
E
n
2
+
C
1
τ
σ
2
+
ρ
2
+
τ
2
+
h
2
≤
E
0
2
+
C
1
n
τ
σ
2
+
ρ
2
+
τ
2
+
h
2
≤
C
σ
2
+
ρ
2
+
τ
2
+
h
2
.
This completes the proof.


## 4. The Second-Order Accurate Implicit Numerical Method for RSDO-ADE in 2D Case

In this section, a spatially second-order accurate alternating direction difference method for the SRDO-ADE in 2D case is proposed. The stability and convergence of this method are discussed.


### 4.1. The Second-Order Accurate Implicit Alternating Direction Method for the RSDO-ADE in 2D Case

We consider the following RSDO-ADE in 2D case:(40)
∂
u
x
,
y
,
t
∂
t
=
∫
0
1
P
α
∂
α
u
x
,
y
,
t
∂
x
α
d
α
+
∫
0
1
P
α
∂
α
u
x
,
y
,
t
∂
y
α
d
α
+
∫
1
2
Q
β
∂
β
u
x
,
y
,
t
∂
x
β
d
β
+
∫
1
2
Q
β
∂
β
u
x
,
y
,
t
∂
y
β
d
β
+
f
x
,
y
,
t
,where P
(
α
) and Q
(
β
) are nonnegative weight functions which satisfy conditions(41)
0
≤
P
α
,
P
α
≢
0
,
α
∈
0,1
,
0
<
∫
0
1
P
α
d
α
<
∞
,
0
≤
Q
β
,
Q
β
≢
0
,
β
∈
1,2
,
0
<
∫
1
2
Q
β
d
β
<
∞
,and ∂
α
u
(
x
,
y
,
t
)
/
∂
x
α, ∂
α
u
(
x
,
y
,
t
)
/
∂
x
β, ∂
α
u
(
x
,
y
,
t
)
/
∂
y
α, and ∂
α
u
(
x
,
y
,
t
)
/
∂
y
β are the Riesz space fractional derivative operators defined the same as Section 2.Now we consider (40) in a square domain Ω
=
[
0
,
L
]
×
[
0
,
L
] with the following initial and boundary conditions:(42)
u
0
,
y
,
t
=
u
L
,
y
,
t
=
0
,
u
x
,
0
,
t
=
u
x
,
L
,
t
=
0
,
t
∈
0
,
T
,
(43)
u
x
,
y
,
0
=
ψ
x
,
y
,
x
,
y
∈
Ω
.Firstly, we discretize the integral intervals[
0,1
] of α and [
1,2
] of β by the grid 0
=
ξ
0
<
ξ
1
<
⋯
<
ξ
S
=
1 and 1
=
η
0
<
η
1
<
⋯
<
η
S
¯
=
2 and denote Δ
ξ
l
=
ξ
l
-
ξ
l
-
1
=
1
/
S
=
σ and α
l
=
ξ
l
+
ξ
l
-
1
/
2
=
2
l
-
1
/
2
S
,
l
=
1,2
,
…
,
S
,
S
∈
N. Consider Δ
η
l
=
η
l
-
η
l
-
1
=
1
/
S
¯
=
ρ and β
l
=
η
l
+
η
l
-
1
/
2
=
1
+
2
l
-
1
/
2
S
¯
,
l
=
1,2
,
…
,
S
¯
,
S
¯
∈
N.Then by using the similar method in Section2, we can obtain(44)
∫
0
1
P
α
∂
α
u
x
,
y
,
t
∂
x
α
d
α
=
∑
l
=
1
S
P
α
l
∂
α
l
u
x
,
y
,
t
∂
x
α
l
Δ
ξ
l
+
O
σ
2
,
∫
0
1
P
α
∂
α
u
x
,
y
,
t
∂
y
α
d
α
=
∑
l
=
1
S
P
α
l
∂
α
l
u
x
,
y
,
t
∂
y
α
l
Δ
ξ
l
+
O
σ
2
,
∫
1
2
Q
β
∂
β
u
x
,
y
,
t
∂
x
β
d
β
=
∑
l
=
1
S
¯
Q
β
l
∂
β
l
u
x
,
y
,
t
∂
x
β
l
Δ
η
l
+
O
ρ
2
,
∫
1
2
Q
β
∂
β
u
x
,
y
,
t
∂
y
β
d
β
=
∑
l
=
1
S
¯
Q
β
l
∂
β
l
u
x
,
y
,
t
∂
y
β
l
Δ
η
l
+
O
ρ
2
.Thus RSDO-ADE in 2D case is now transformed into the following multiterm fractional equation:(45)
∂
u
x
,
y
,
t
∂
t
=
∑
l
=
1
S
P
α
l
∂
α
l
u
x
,
y
,
t
∂
x
α
l
1
S
+
∑
l
=
1
S
P
α
l
∂
α
l
u
x
,
y
,
t
∂
y
α
l
1
S
+
∑
l
=
1
S
¯
Q
β
l
∂
β
l
u
x
,
y
,
t
∂
x
β
l
1
S
¯
+
∑
l
=
1
S
¯
Q
β
l
∂
β
l
u
x
,
y
,
t
∂
y
β
l
1
S
¯
+
f
x
,
y
,
t
+
O
σ
2
+
ρ
2
.For the numerical simulation of (45), Let h
=
L
/
M be the spatial grid size in the x-direction and in the y-direction; let τ
=
T
/
N be the time step; let x
i
=
i
h
,
i
=
0,1
,
…
,
M; let y
j
=
j
h
,
j
=
0,1
,
…
,
M; let t
n
=
n
τ
,
n
=
0,1
,
…
,
N. Define u
i
,
j
n as the numerical solution to u
(
x
i
,
y
j
,
t
n
).The Riesz space fractional derivative operator is discretized as follows [32, 33]:(46)
∂
α
l
u
x
i
,
y
j
,
t
n
∂
x
α
l
=
-
1
h
α
l
∑
p
=
1
M
-
1
g
i
-
p
α
l
u
x
p
,
y
j
,
t
n
+
O
h
2
,
0
<
α
l
<
1
,
i
=
1,2
,
…
,
M
-
1
,
∂
α
l
u
x
i
,
y
j
,
t
n
∂
y
α
l
=
-
1
h
α
l
∑
q
=
1
M
-
1
g
j
-
q
α
l
u
x
i
,
y
q
,
t
n
+
O
h
2
,
0
<
α
l
<
1
,
j
=
1,2
,
…
,
M
-
1
,
∂
β
l
u
x
i
,
y
j
,
t
n
∂
x
β
l
=
-
1
h
β
l
∑
p
=
1
M
-
1
g
i
-
p
β
l
u
x
p
,
y
j
,
t
n
+
O
h
2
,
1
<
β
l
≤
2
,
i
=
1,2
,
…
,
M
-
1
,
∂
β
l
u
x
i
,
y
j
,
t
n
∂
y
β
l
=
-
1
h
β
l
∑
q
=
1
M
-
1
g
j
-
q
β
l
u
x
i
,
y
q
,
t
n
+
O
h
2
,
1
<
β
l
≤
2
,
j
=
1,2
,
…
,
M
-
1
,where g
l
α
l and g
l
β
l are defined in Section 2.2.Using Crank-Nicholson method and second-order accurate implicit finite difference scheme, we obtain the following discrete form for RSDO-ADE in 2D case:(47)
u
x
i
,
y
j
,
t
n
+
1
-
u
x
i
,
y
j
,
t
n
τ
=
-
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
∑
p
=
1
M
-
1
g
i
-
p
α
l
u
x
p
,
y
j
,
t
n
-
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
∑
p
=
1
M
-
1
g
i
-
p
α
l
u
x
p
,
y
j
,
t
n
+
1
-
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
∑
q
=
1
M
-
1
g
j
-
q
α
l
u
x
i
,
y
q
,
t
n
-
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
∑
q
=
1
M
-
1
g
j
-
q
α
l
u
x
i
,
y
q
,
t
n
+
1
-
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
∑
p
=
1
M
-
1
g
i
-
p
β
l
u
x
p
,
y
j
,
t
n
-
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
∑
p
=
1
M
-
1
g
i
-
p
β
l
u
x
p
,
y
j
,
t
n
+
1
-
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
∑
q
=
1
M
-
1
g
j
-
q
β
l
u
x
i
,
y
q
,
t
n
-
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
∑
q
=
1
M
-
1
g
j
-
q
β
l
u
x
i
,
y
q
,
t
n
+
1
+
f
x
i
,
y
j
,
t
n
+
1
/
2
+
R
i
,
j
n
+
1
/
2
,
i
,
j
=
1,2
,
…
,
M
-
1
;
n
=
0,1
,
…
,
N
-
1
,where R
i
,
j
n
+
1
/
2
=
O
(
σ
2
+
ρ
2
+
h
2
+
τ
2
) and f
(
x
i
,
y
j
,
t
n
+
1
/
2
)
=
1
/
2
(
f
(
x
i
,
y
j
,
t
n
+
1
)
+
f
(
x
i
,
y
j
,
t
n
)
).By omitting the local truncation error termR
i
,
j
n
+
1
/
2 in (47), we obtain the following second-order accurate implicit numerical method for RSDO-ADE in 2D case:(48)
u
i
,
j
n
+
1
-
u
i
,
j
n
τ
=
-
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
∑
p
=
1
M
-
1
g
i
-
p
α
l
u
p
,
j
n
-
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
∑
p
=
1
M
-
1
g
i
-
p
α
l
u
p
,
j
n
+
1
-
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
∑
q
=
1
M
-
1
g
j
-
q
α
l
u
i
,
q
n
-
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
∑
q
=
1
M
-
1
g
j
-
q
α
l
u
i
,
q
n
+
1
-
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
∑
p
=
1
M
-
1
g
i
-
p
β
l
u
p
,
j
n
-
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
∑
p
=
1
M
-
1
g
i
-
p
β
l
u
p
,
j
n
+
1
-
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
∑
q
=
1
M
-
1
g
j
-
q
β
l
u
i
,
q
n
-
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
∑
q
=
1
M
-
1
g
j
-
q
β
l
u
i
,
q
n
+
1
+
f
x
i
,
y
j
,
t
n
+
1
/
2
,
i
,
j
=
1,2
,
…
,
M
-
1
;
n
=
0,1
,
…
,
N
-
1
,
u
0
,
j
n
=
u
M
,
j
n
=
0
,
u
i
,
0
n
=
u
0
,
M
n
=
0
,
n
=
1,2
,
…
,
N
,
u
i
,
j
0
=
ψ
x
i
,
y
j
,
i
,
j
=
0,1
,
…
,
M
.Define the following fractional partial difference operators:(49)
δ
x
u
i
,
j
n
=
-
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
∑
p
=
1
M
-
1
g
i
-
p
α
l
u
p
,
j
n
-
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
∑
p
=
1
M
-
1
g
i
-
p
β
l
u
p
,
j
n
,
δ
y
u
i
,
j
n
=
-
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
∑
q
=
1
M
-
1
g
j
-
q
α
l
u
i
,
q
n
-
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
∑
q
=
1
M
-
1
g
j
-
q
β
l
u
i
,
q
n
.The second-order accurate implicit numerical method for RSDO-ADE in 2D case may be written in the following operator form:(50)
1
-
τ
δ
x
-
τ
δ
y
u
i
,
j
n
+
1
=
1
+
τ
δ
x
+
τ
δ
y
u
i
,
j
n
+
τ
f
x
i
,
y
j
,
t
n
+
1
/
2
.We introduce two additional perturbation errors equal toτ
2
δ
x
δ
y
u
i
,
j
n
+
1 and τ
2
δ
x
δ
y
u
i
,
j
n. Equation (50) is then written in the following directional separation product form:(51)
1
-
τ
δ
x
1
-
τ
δ
y
u
i
,
j
n
+
1
=
1
+
τ
δ
x
1
+
τ
δ
y
u
i
,
j
n
+
τ
f
x
i
,
y
j
,
t
n
+
1
/
2
.The additional perturbation errors are not large compared to the approximation errors for the other terms in (50), and hence (51), which is called the implicit alternation direction method, is consistent with order O
(
σ
2
+
ρ
2
+
h
2
+
τ
2
).Computationally, the implicit alternation direction method defined by (51) can now be solved by the following iterative scheme, at time t
=
t
n
+
1.Step 1.
Solve the problem in thex-direction for each fixed y
j to obtain an intermediate solution u
i
,
j
∗ in the form(52)
1
-
τ
δ
x
u
i
,
j
∗
=
1
+
τ
δ
x
1
+
τ
δ
y
u
i
,
j
n
+
τ
f
x
i
,
y
j
,
t
n
+
1
/
2
,
1
≤
i
≤
M
-
1
.Step 2.
Then solve it in they-direction for each fixed x
i:(53)
1
-
τ
δ
y
u
i
,
j
n
+
1
=
u
i
,
j
∗
,
1
≤
j
≤
M
-
1
.The initial and boundary conditions for the numerical solution u
i
,
j
n
+
1 and u
i
,
j
n are defined from the given initial and boundary conditions (see [34]).The initial and boundary conditions for numerical solutionsu
i
,
j
n
+
1 and u
i
,
j
n are defined from the given initial and boundary conditions. Prior to carrying out step one of solving (52), the boundary conditions for the intermediate solution u
i
,
j
∗ should be set from (53) (which incorporates the values of u
i
,
j
n
+
1 at the boundary); otherwise the order of convergence will be adversely affected. Specifically, for homogeneous Dirichlet boundary conditions (42), we have(54)
u
0
,
j
n
+
1
=
u
0
,
y
j
,
t
n
+
1
=
0
;
u
M
,
j
n
+
1
=
u
L
,
y
j
,
t
n
+
1
=
0
;
u
i
,
0
n
+
1
=
u
x
i
,
0
,
t
n
+
1
=
0
;
u
i
,
M
n
+
1
=
u
x
i
,
L
,
t
n
+
1
=
0
.Thus, we compute the boundary values for u
∗ from(55)
u
0
,
j
∗
=
1
-
τ
δ
y
u
0
,
j
n
+
1
,
u
M
,
j
∗
=
1
-
τ
δ
y
u
M
,
j
n
+
1
.


### 4.2. Numerical Analysis of the Implicit Alternating Direction Method for the RSDO-ADE in 2D Case

In this section, we discuss the stability and convergence of alternating direction method (51) for the SRDO-ADE in 2D case. We need to rewrite (51), (52), and (53) in matrix form [34].Then (51) may be written as(56)
R
H
U
n
+
1
=
R
¯
H
¯
U
n
+
F
n
+
1
/
2
,where matrices R, R
¯, H, and H
¯ represent operators 1
-
τ
δ
x, 1
+
τ
δ
x, 1
-
τ
δ
y, and 1
+
τ
δ
y, respectively, where(57)
U
n
=
u
1,1
n
,
…
,
u
M
-
1,1
n
,
u
1,2
n
,
…
,
u
M
-
1,2
n
,
…
,
u
1
,
M
-
1
n
,
…
,
u
M
-
1
,
M
-
1
k
T
,
F
n
+
1
/
2
=
f
1,1
n
+
1
/
2
,
f
1,2
n
+
1
/
2
,
…
,
f
M
-
1
,
M
-
1
n
+
1
/
2
T
.MatrixR is a block diagonal matrix of (
M
-
1
)
×
(
M
-
1
) blocks of (
M
-
1
)
×
(
M
-
1
) square matrices A
j resulting from (52). We may write R
=
diag
⁡
(
A
1
,
A
2
,
…
,
A
M
-
1
). Similarly, matrix H is a block matrix of (
M
-
1
)
×
(
M
-
1
) blocks of (
M
-
1
)
×
(
M
-
1
) square diagonal matrices resulting from (52). That is, we may write H
=
[
H
i
,
j
], where each H
i
,
j is (
M
-
1
)
×
(
M
-
1
) matrix, such that H
i
,
j is a diagonal matrix H
i
,
j
=
diag
⁡
(
b
i
,
j
,
b
i
,
j
,
…
,
b
i
,
j
), and where the notation b
i
,
j refers to the (
i
,
j
)th entry of matrix B
i defined. We note that matrices A
j
(
j
=
1,2
,
…
,
M
-
1
) and B
i
(
i
=
1,2
,
…
,
M
-
1
) are strictly diagonally dominant. Because their diagonal elements are all positive, these matrices are symmetric and positive definite. Here A
j
=
[
a
i
,
k
(
j
)
]
(
M
-
1
)
×
(
M
-
1
),(58)
a
i
,
k
j
=
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
g
i
-
k
α
l
+
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
g
i
-
k
β
l
,
for
k
<
i
;
1
+
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
g
0
α
l
+
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
g
0
β
l
,
for
k
=
i
;
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
g
k
-
i
α
l
+
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
g
k
-
i
β
l
,
for
k
>
i
.Consider B
i
=
[
b
k
,
j
(
i
)
]
(
M
-
1
)
×
(
M
-
1
),(59)
b
k
,
j
i
=
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
g
k
-
j
α
l
+
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
g
k
-
j
β
l
,
for
j
<
k
;
1
+
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
g
0
α
l
+
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
g
0
β
l
,
for
j
=
k
;
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
g
j
-
k
α
l
+
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
g
j
-
k
β
l
,
for
j
>
k
.Similarly, we can define matrices R
¯ and H
¯.To prove the stability and convergence of the implicit alternating direction method, we need the following lemma in [34].Lemma 8.
LetX
=
[
x
1
,
x
2
,
…
,
x
m
]
T, X
∞
=
max
1
≤
i
≤
m
⁡
x
i. If matrix D
=
(
d
i
,
j
)
m
×
m satisfies conditions(60)
∑
j
=
1
,
j
≠
i
m
d
i
,
j
≤
d
i
,
i
-
1
,
i
=
1,2
,
…
,
m
,then(61)
X
∞
≤
D
X
∞
.Letu
i
,
j
n and u
~
i
,
j
n be the numerical and approximate solutions of implicit alternating direction method (51), respectively, and set(62)
Y
n
=
ε
1,1
n
,
…
,
ε
M
-
1,1
n
,
…
,
ε
M
-
1
,
M
-
1
n
T
,where ε
i
,
j
n
=
u
~
i
,
j
n
-
u
i
,
j
n.Theorem 9.
The second-order accurate implicit alternating direction method (51) for the RSDO-ADE in 2D case is unconditionally stable and there is a positive constant C
1
∗ such that(63)
Y
n
+
1
∞
≤
C
1
∗
Y
0
∞
,
n
=
0,1
,
2
,
…
.Proof.
The errorY
n satisfies the following equation:(64)
R
H
Y
n
+
1
=
R
¯
H
¯
Y
n
.
SinceR, R
¯, H, and H
¯ satisfy conditions of Lemma 8, we obtain(65)
Y
n
+
1
∞
≤
H
Y
n
+
1
∞
≤
R
H
Y
n
+
1
∞
=
R
¯
H
¯
Y
n
∞
≤
R
¯
∞
H
¯
∞
Y
n
∞
≤
C
1
∗
Y
0
∞
.
Now let us consider the convergence of second-order accurate implicit alternating direction method (51). Let u
(
x
i
,
y
j
,
t
n
) be the exact solution of the RSDO-ADE in 2D case, let u
i
,
j
n be the numerical solution of second-order accurate implicit alternating direction method (51). Let e
i
,
j
n
=
u
(
x
i
,
y
j
,
t
n
)
-
u
i
,
j
n and(66)
E
n
=
e
1,1
n
,
…
,
e
M
-
1,1
n
,
…
,
e
M
-
1
,
M
-
1
n
T
.Theorem 10.
Second-order accurate implicit alternating direction method (51) is convergent and there is a positive constant C
2
∗ such that(67)
E
n
+
1
∞
≤
C
2
∗
σ
2
+
ρ
2
+
τ
2
+
h
2
,
n
=
0,1
,
2
,
…
,that is, when u
i
,
j
n tends to u
(
x
i
,
y
j
,
t
n
) at any fixed point at τ and h both tend to zero.Proof.
The errorE
n satisfies the following equation:(68)
R
H
E
n
+
1
=
E
n
+
τ
R
n
+
1
/
2
,where R
n
+
1
/
2
=
(
R
1,1
n
+
1
/
2
,
R
1,2
n
+
1
/
2
,
…
,
R
M
-
1
,
M
-
1
n
+
1
/
2
)
T.
SinceR, R
¯, H, and H
¯ satisfy conditions of Lemma 8, we obtain(69)
E
n
+
1
∞
≤
H
E
n
+
1
∞
≤
R
H
E
n
+
1
∞
=
R
¯
H
¯
E
n
+
τ
R
n
+
1
/
2
∞
≤
R
¯
∞
H
¯
∞
E
n
∞
+
τ
R
n
+
1
/
2
∞
≤
n
+
1
τ
C
2
σ
2
+
ρ
2
+
τ
2
+
h
2
≤
C
2
∗
σ
2
+
ρ
2
+
τ
2
+
h
2
.Therefore second-order accurate implicit alternating direction method (51) is convergent.
This completes the proof.


## 4.1. The Second-Order Accurate Implicit Alternating Direction Method for the RSDO-ADE in 2D Case

We consider the following RSDO-ADE in 2D case:(40)
∂
u
x
,
y
,
t
∂
t
=
∫
0
1
P
α
∂
α
u
x
,
y
,
t
∂
x
α
d
α
+
∫
0
1
P
α
∂
α
u
x
,
y
,
t
∂
y
α
d
α
+
∫
1
2
Q
β
∂
β
u
x
,
y
,
t
∂
x
β
d
β
+
∫
1
2
Q
β
∂
β
u
x
,
y
,
t
∂
y
β
d
β
+
f
x
,
y
,
t
,where P
(
α
) and Q
(
β
) are nonnegative weight functions which satisfy conditions(41)
0
≤
P
α
,
P
α
≢
0
,
α
∈
0,1
,
0
<
∫
0
1
P
α
d
α
<
∞
,
0
≤
Q
β
,
Q
β
≢
0
,
β
∈
1,2
,
0
<
∫
1
2
Q
β
d
β
<
∞
,and ∂
α
u
(
x
,
y
,
t
)
/
∂
x
α, ∂
α
u
(
x
,
y
,
t
)
/
∂
x
β, ∂
α
u
(
x
,
y
,
t
)
/
∂
y
α, and ∂
α
u
(
x
,
y
,
t
)
/
∂
y
β are the Riesz space fractional derivative operators defined the same as Section 2.Now we consider (40) in a square domain Ω
=
[
0
,
L
]
×
[
0
,
L
] with the following initial and boundary conditions:(42)
u
0
,
y
,
t
=
u
L
,
y
,
t
=
0
,
u
x
,
0
,
t
=
u
x
,
L
,
t
=
0
,
t
∈
0
,
T
,
(43)
u
x
,
y
,
0
=
ψ
x
,
y
,
x
,
y
∈
Ω
.Firstly, we discretize the integral intervals[
0,1
] of α and [
1,2
] of β by the grid 0
=
ξ
0
<
ξ
1
<
⋯
<
ξ
S
=
1 and 1
=
η
0
<
η
1
<
⋯
<
η
S
¯
=
2 and denote Δ
ξ
l
=
ξ
l
-
ξ
l
-
1
=
1
/
S
=
σ and α
l
=
ξ
l
+
ξ
l
-
1
/
2
=
2
l
-
1
/
2
S
,
l
=
1,2
,
…
,
S
,
S
∈
N. Consider Δ
η
l
=
η
l
-
η
l
-
1
=
1
/
S
¯
=
ρ and β
l
=
η
l
+
η
l
-
1
/
2
=
1
+
2
l
-
1
/
2
S
¯
,
l
=
1,2
,
…
,
S
¯
,
S
¯
∈
N.Then by using the similar method in Section2, we can obtain(44)
∫
0
1
P
α
∂
α
u
x
,
y
,
t
∂
x
α
d
α
=
∑
l
=
1
S
P
α
l
∂
α
l
u
x
,
y
,
t
∂
x
α
l
Δ
ξ
l
+
O
σ
2
,
∫
0
1
P
α
∂
α
u
x
,
y
,
t
∂
y
α
d
α
=
∑
l
=
1
S
P
α
l
∂
α
l
u
x
,
y
,
t
∂
y
α
l
Δ
ξ
l
+
O
σ
2
,
∫
1
2
Q
β
∂
β
u
x
,
y
,
t
∂
x
β
d
β
=
∑
l
=
1
S
¯
Q
β
l
∂
β
l
u
x
,
y
,
t
∂
x
β
l
Δ
η
l
+
O
ρ
2
,
∫
1
2
Q
β
∂
β
u
x
,
y
,
t
∂
y
β
d
β
=
∑
l
=
1
S
¯
Q
β
l
∂
β
l
u
x
,
y
,
t
∂
y
β
l
Δ
η
l
+
O
ρ
2
.Thus RSDO-ADE in 2D case is now transformed into the following multiterm fractional equation:(45)
∂
u
x
,
y
,
t
∂
t
=
∑
l
=
1
S
P
α
l
∂
α
l
u
x
,
y
,
t
∂
x
α
l
1
S
+
∑
l
=
1
S
P
α
l
∂
α
l
u
x
,
y
,
t
∂
y
α
l
1
S
+
∑
l
=
1
S
¯
Q
β
l
∂
β
l
u
x
,
y
,
t
∂
x
β
l
1
S
¯
+
∑
l
=
1
S
¯
Q
β
l
∂
β
l
u
x
,
y
,
t
∂
y
β
l
1
S
¯
+
f
x
,
y
,
t
+
O
σ
2
+
ρ
2
.For the numerical simulation of (45), Let h
=
L
/
M be the spatial grid size in the x-direction and in the y-direction; let τ
=
T
/
N be the time step; let x
i
=
i
h
,
i
=
0,1
,
…
,
M; let y
j
=
j
h
,
j
=
0,1
,
…
,
M; let t
n
=
n
τ
,
n
=
0,1
,
…
,
N. Define u
i
,
j
n as the numerical solution to u
(
x
i
,
y
j
,
t
n
).The Riesz space fractional derivative operator is discretized as follows [32, 33]:(46)
∂
α
l
u
x
i
,
y
j
,
t
n
∂
x
α
l
=
-
1
h
α
l
∑
p
=
1
M
-
1
g
i
-
p
α
l
u
x
p
,
y
j
,
t
n
+
O
h
2
,
0
<
α
l
<
1
,
i
=
1,2
,
…
,
M
-
1
,
∂
α
l
u
x
i
,
y
j
,
t
n
∂
y
α
l
=
-
1
h
α
l
∑
q
=
1
M
-
1
g
j
-
q
α
l
u
x
i
,
y
q
,
t
n
+
O
h
2
,
0
<
α
l
<
1
,
j
=
1,2
,
…
,
M
-
1
,
∂
β
l
u
x
i
,
y
j
,
t
n
∂
x
β
l
=
-
1
h
β
l
∑
p
=
1
M
-
1
g
i
-
p
β
l
u
x
p
,
y
j
,
t
n
+
O
h
2
,
1
<
β
l
≤
2
,
i
=
1,2
,
…
,
M
-
1
,
∂
β
l
u
x
i
,
y
j
,
t
n
∂
y
β
l
=
-
1
h
β
l
∑
q
=
1
M
-
1
g
j
-
q
β
l
u
x
i
,
y
q
,
t
n
+
O
h
2
,
1
<
β
l
≤
2
,
j
=
1,2
,
…
,
M
-
1
,where g
l
α
l and g
l
β
l are defined in Section 2.2.Using Crank-Nicholson method and second-order accurate implicit finite difference scheme, we obtain the following discrete form for RSDO-ADE in 2D case:(47)
u
x
i
,
y
j
,
t
n
+
1
-
u
x
i
,
y
j
,
t
n
τ
=
-
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
∑
p
=
1
M
-
1
g
i
-
p
α
l
u
x
p
,
y
j
,
t
n
-
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
∑
p
=
1
M
-
1
g
i
-
p
α
l
u
x
p
,
y
j
,
t
n
+
1
-
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
∑
q
=
1
M
-
1
g
j
-
q
α
l
u
x
i
,
y
q
,
t
n
-
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
∑
q
=
1
M
-
1
g
j
-
q
α
l
u
x
i
,
y
q
,
t
n
+
1
-
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
∑
p
=
1
M
-
1
g
i
-
p
β
l
u
x
p
,
y
j
,
t
n
-
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
∑
p
=
1
M
-
1
g
i
-
p
β
l
u
x
p
,
y
j
,
t
n
+
1
-
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
∑
q
=
1
M
-
1
g
j
-
q
β
l
u
x
i
,
y
q
,
t
n
-
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
∑
q
=
1
M
-
1
g
j
-
q
β
l
u
x
i
,
y
q
,
t
n
+
1
+
f
x
i
,
y
j
,
t
n
+
1
/
2
+
R
i
,
j
n
+
1
/
2
,
i
,
j
=
1,2
,
…
,
M
-
1
;
n
=
0,1
,
…
,
N
-
1
,where R
i
,
j
n
+
1
/
2
=
O
(
σ
2
+
ρ
2
+
h
2
+
τ
2
) and f
(
x
i
,
y
j
,
t
n
+
1
/
2
)
=
1
/
2
(
f
(
x
i
,
y
j
,
t
n
+
1
)
+
f
(
x
i
,
y
j
,
t
n
)
).By omitting the local truncation error termR
i
,
j
n
+
1
/
2 in (47), we obtain the following second-order accurate implicit numerical method for RSDO-ADE in 2D case:(48)
u
i
,
j
n
+
1
-
u
i
,
j
n
τ
=
-
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
∑
p
=
1
M
-
1
g
i
-
p
α
l
u
p
,
j
n
-
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
∑
p
=
1
M
-
1
g
i
-
p
α
l
u
p
,
j
n
+
1
-
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
∑
q
=
1
M
-
1
g
j
-
q
α
l
u
i
,
q
n
-
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
∑
q
=
1
M
-
1
g
j
-
q
α
l
u
i
,
q
n
+
1
-
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
∑
p
=
1
M
-
1
g
i
-
p
β
l
u
p
,
j
n
-
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
∑
p
=
1
M
-
1
g
i
-
p
β
l
u
p
,
j
n
+
1
-
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
∑
q
=
1
M
-
1
g
j
-
q
β
l
u
i
,
q
n
-
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
∑
q
=
1
M
-
1
g
j
-
q
β
l
u
i
,
q
n
+
1
+
f
x
i
,
y
j
,
t
n
+
1
/
2
,
i
,
j
=
1,2
,
…
,
M
-
1
;
n
=
0,1
,
…
,
N
-
1
,
u
0
,
j
n
=
u
M
,
j
n
=
0
,
u
i
,
0
n
=
u
0
,
M
n
=
0
,
n
=
1,2
,
…
,
N
,
u
i
,
j
0
=
ψ
x
i
,
y
j
,
i
,
j
=
0,1
,
…
,
M
.Define the following fractional partial difference operators:(49)
δ
x
u
i
,
j
n
=
-
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
∑
p
=
1
M
-
1
g
i
-
p
α
l
u
p
,
j
n
-
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
∑
p
=
1
M
-
1
g
i
-
p
β
l
u
p
,
j
n
,
δ
y
u
i
,
j
n
=
-
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
∑
q
=
1
M
-
1
g
j
-
q
α
l
u
i
,
q
n
-
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
∑
q
=
1
M
-
1
g
j
-
q
β
l
u
i
,
q
n
.The second-order accurate implicit numerical method for RSDO-ADE in 2D case may be written in the following operator form:(50)
1
-
τ
δ
x
-
τ
δ
y
u
i
,
j
n
+
1
=
1
+
τ
δ
x
+
τ
δ
y
u
i
,
j
n
+
τ
f
x
i
,
y
j
,
t
n
+
1
/
2
.We introduce two additional perturbation errors equal toτ
2
δ
x
δ
y
u
i
,
j
n
+
1 and τ
2
δ
x
δ
y
u
i
,
j
n. Equation (50) is then written in the following directional separation product form:(51)
1
-
τ
δ
x
1
-
τ
δ
y
u
i
,
j
n
+
1
=
1
+
τ
δ
x
1
+
τ
δ
y
u
i
,
j
n
+
τ
f
x
i
,
y
j
,
t
n
+
1
/
2
.The additional perturbation errors are not large compared to the approximation errors for the other terms in (50), and hence (51), which is called the implicit alternation direction method, is consistent with order O
(
σ
2
+
ρ
2
+
h
2
+
τ
2
).Computationally, the implicit alternation direction method defined by (51) can now be solved by the following iterative scheme, at time t
=
t
n
+
1.Step 1.
Solve the problem in thex-direction for each fixed y
j to obtain an intermediate solution u
i
,
j
∗ in the form(52)
1
-
τ
δ
x
u
i
,
j
∗
=
1
+
τ
δ
x
1
+
τ
δ
y
u
i
,
j
n
+
τ
f
x
i
,
y
j
,
t
n
+
1
/
2
,
1
≤
i
≤
M
-
1
.Step 2.
Then solve it in they-direction for each fixed x
i:(53)
1
-
τ
δ
y
u
i
,
j
n
+
1
=
u
i
,
j
∗
,
1
≤
j
≤
M
-
1
.The initial and boundary conditions for the numerical solution u
i
,
j
n
+
1 and u
i
,
j
n are defined from the given initial and boundary conditions (see [34]).The initial and boundary conditions for numerical solutionsu
i
,
j
n
+
1 and u
i
,
j
n are defined from the given initial and boundary conditions. Prior to carrying out step one of solving (52), the boundary conditions for the intermediate solution u
i
,
j
∗ should be set from (53) (which incorporates the values of u
i
,
j
n
+
1 at the boundary); otherwise the order of convergence will be adversely affected. Specifically, for homogeneous Dirichlet boundary conditions (42), we have(54)
u
0
,
j
n
+
1
=
u
0
,
y
j
,
t
n
+
1
=
0
;
u
M
,
j
n
+
1
=
u
L
,
y
j
,
t
n
+
1
=
0
;
u
i
,
0
n
+
1
=
u
x
i
,
0
,
t
n
+
1
=
0
;
u
i
,
M
n
+
1
=
u
x
i
,
L
,
t
n
+
1
=
0
.Thus, we compute the boundary values for u
∗ from(55)
u
0
,
j
∗
=
1
-
τ
δ
y
u
0
,
j
n
+
1
,
u
M
,
j
∗
=
1
-
τ
δ
y
u
M
,
j
n
+
1
.


## 4.2. Numerical Analysis of the Implicit Alternating Direction Method for the RSDO-ADE in 2D Case

In this section, we discuss the stability and convergence of alternating direction method (51) for the SRDO-ADE in 2D case. We need to rewrite (51), (52), and (53) in matrix form [34].Then (51) may be written as(56)
R
H
U
n
+
1
=
R
¯
H
¯
U
n
+
F
n
+
1
/
2
,where matrices R, R
¯, H, and H
¯ represent operators 1
-
τ
δ
x, 1
+
τ
δ
x, 1
-
τ
δ
y, and 1
+
τ
δ
y, respectively, where(57)
U
n
=
u
1,1
n
,
…
,
u
M
-
1,1
n
,
u
1,2
n
,
…
,
u
M
-
1,2
n
,
…
,
u
1
,
M
-
1
n
,
…
,
u
M
-
1
,
M
-
1
k
T
,
F
n
+
1
/
2
=
f
1,1
n
+
1
/
2
,
f
1,2
n
+
1
/
2
,
…
,
f
M
-
1
,
M
-
1
n
+
1
/
2
T
.MatrixR is a block diagonal matrix of (
M
-
1
)
×
(
M
-
1
) blocks of (
M
-
1
)
×
(
M
-
1
) square matrices A
j resulting from (52). We may write R
=
diag
⁡
(
A
1
,
A
2
,
…
,
A
M
-
1
). Similarly, matrix H is a block matrix of (
M
-
1
)
×
(
M
-
1
) blocks of (
M
-
1
)
×
(
M
-
1
) square diagonal matrices resulting from (52). That is, we may write H
=
[
H
i
,
j
], where each H
i
,
j is (
M
-
1
)
×
(
M
-
1
) matrix, such that H
i
,
j is a diagonal matrix H
i
,
j
=
diag
⁡
(
b
i
,
j
,
b
i
,
j
,
…
,
b
i
,
j
), and where the notation b
i
,
j refers to the (
i
,
j
)th entry of matrix B
i defined. We note that matrices A
j
(
j
=
1,2
,
…
,
M
-
1
) and B
i
(
i
=
1,2
,
…
,
M
-
1
) are strictly diagonally dominant. Because their diagonal elements are all positive, these matrices are symmetric and positive definite. Here A
j
=
[
a
i
,
k
(
j
)
]
(
M
-
1
)
×
(
M
-
1
),(58)
a
i
,
k
j
=
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
g
i
-
k
α
l
+
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
g
i
-
k
β
l
,
for
k
<
i
;
1
+
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
g
0
α
l
+
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
g
0
β
l
,
for
k
=
i
;
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
g
k
-
i
α
l
+
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
g
k
-
i
β
l
,
for
k
>
i
.Consider B
i
=
[
b
k
,
j
(
i
)
]
(
M
-
1
)
×
(
M
-
1
),(59)
b
k
,
j
i
=
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
g
k
-
j
α
l
+
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
g
k
-
j
β
l
,
for
j
<
k
;
1
+
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
g
0
α
l
+
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
g
0
β
l
,
for
j
=
k
;
1
2
∑
l
=
1
S
P
α
l
S
1
h
α
l
g
j
-
k
α
l
+
1
2
∑
l
=
1
S
¯
Q
β
l
S
¯
1
h
β
l
g
j
-
k
β
l
,
for
j
>
k
.Similarly, we can define matrices R
¯ and H
¯.To prove the stability and convergence of the implicit alternating direction method, we need the following lemma in [34].Lemma 8.
LetX
=
[
x
1
,
x
2
,
…
,
x
m
]
T, X
∞
=
max
1
≤
i
≤
m
⁡
x
i. If matrix D
=
(
d
i
,
j
)
m
×
m satisfies conditions(60)
∑
j
=
1
,
j
≠
i
m
d
i
,
j
≤
d
i
,
i
-
1
,
i
=
1,2
,
…
,
m
,then(61)
X
∞
≤
D
X
∞
.Letu
i
,
j
n and u
~
i
,
j
n be the numerical and approximate solutions of implicit alternating direction method (51), respectively, and set(62)
Y
n
=
ε
1,1
n
,
…
,
ε
M
-
1,1
n
,
…
,
ε
M
-
1
,
M
-
1
n
T
,where ε
i
,
j
n
=
u
~
i
,
j
n
-
u
i
,
j
n.Theorem 9.
The second-order accurate implicit alternating direction method (51) for the RSDO-ADE in 2D case is unconditionally stable and there is a positive constant C
1
∗ such that(63)
Y
n
+
1
∞
≤
C
1
∗
Y
0
∞
,
n
=
0,1
,
2
,
…
.Proof.
The errorY
n satisfies the following equation:(64)
R
H
Y
n
+
1
=
R
¯
H
¯
Y
n
.
SinceR, R
¯, H, and H
¯ satisfy conditions of Lemma 8, we obtain(65)
Y
n
+
1
∞
≤
H
Y
n
+
1
∞
≤
R
H
Y
n
+
1
∞
=
R
¯
H
¯
Y
n
∞
≤
R
¯
∞
H
¯
∞
Y
n
∞
≤
C
1
∗
Y
0
∞
.
Now let us consider the convergence of second-order accurate implicit alternating direction method (51). Let u
(
x
i
,
y
j
,
t
n
) be the exact solution of the RSDO-ADE in 2D case, let u
i
,
j
n be the numerical solution of second-order accurate implicit alternating direction method (51). Let e
i
,
j
n
=
u
(
x
i
,
y
j
,
t
n
)
-
u
i
,
j
n and(66)
E
n
=
e
1,1
n
,
…
,
e
M
-
1,1
n
,
…
,
e
M
-
1
,
M
-
1
n
T
.Theorem 10.
Second-order accurate implicit alternating direction method (51) is convergent and there is a positive constant C
2
∗ such that(67)
E
n
+
1
∞
≤
C
2
∗
σ
2
+
ρ
2
+
τ
2
+
h
2
,
n
=
0,1
,
2
,
…
,that is, when u
i
,
j
n tends to u
(
x
i
,
y
j
,
t
n
) at any fixed point at τ and h both tend to zero.Proof.
The errorE
n satisfies the following equation:(68)
R
H
E
n
+
1
=
E
n
+
τ
R
n
+
1
/
2
,where R
n
+
1
/
2
=
(
R
1,1
n
+
1
/
2
,
R
1,2
n
+
1
/
2
,
…
,
R
M
-
1
,
M
-
1
n
+
1
/
2
)
T.
SinceR, R
¯, H, and H
¯ satisfy conditions of Lemma 8, we obtain(69)
E
n
+
1
∞
≤
H
E
n
+
1
∞
≤
R
H
E
n
+
1
∞
=
R
¯
H
¯
E
n
+
τ
R
n
+
1
/
2
∞
≤
R
¯
∞
H
¯
∞
E
n
∞
+
τ
R
n
+
1
/
2
∞
≤
n
+
1
τ
C
2
σ
2
+
ρ
2
+
τ
2
+
h
2
≤
C
2
∗
σ
2
+
ρ
2
+
τ
2
+
h
2
.Therefore second-order accurate implicit alternating direction method (51) is convergent.
This completes the proof.


## 5. Numerical Results

In order to illustrate the behaviour of our numerical method and demonstrate the effectiveness of our theoretical analysis, some examples are given.Example 1.
Consider the following Riesz space distributed-order advection-dispersion equation in 1D case:(70)
∂
u
x
,
t
∂
t
=
∫
0
1
P
α
∂
α
u
x
,
t
∂
x
α
d
α
+
∫
1
2
Q
β
∂
β
u
x
,
t
∂
x
β
d
β
+
f
x
,
t
,where(71)
P
α
=
-
2
Γ
5
-
α
cos
⁡
π
α
2
,
Q
β
=
-
2
Γ
5
-
β
cos
⁡
π
β
2
.Now we consider (70) in the finite domain [
0,1
] with the following initial and boundary conditions:(72)
u
0
,
t
=
0
,
u
1
,
t
=
0
,
t
∈
0,1
,
u
x
,
0
=
x
2
1
-
x
2
,
x
∈
0,1
,where(73)
f
x
,
t
=
e
t
x
2
1
-
x
2
-
2
e
t
2
x
2
-
2
x
ln
⁡
1
-
x
-
3
+
6
x
3
-
5
x
2
+
x
-
2
ln
⁡
1
-
x
-
2
+
6
x
2
x
3
-
2
x
2
-
x
+
1
ln
⁡
1
-
x
-
1
+
2
x
x
-
1
ln
-
3
⁡
x
+
x
6
x
2
-
13
x
+
5
ln
-
2
⁡
x
+
6
x
x
-
1
2
x
2
-
4
x
+
1
ln
-
1
⁡
x
-
2
e
t
-
2
x
ln
⁡
1
-
x
-
3
+
6
x
2
-
x
-
2
ln
⁡
1
-
x
-
2
+
2
-
6
x
3
+
3
x
2
+
3
x
-
1
ln
⁡
1
-
x
-
1
+
2
x
-
1
ln
-
3
⁡
x
+
6
x
2
-
11
x
+
3
ln
-
2
⁡
x
+
2
6
x
3
-
15
x
2
+
9
x
-
1
ln
-
1
⁡
x
.The exact solution of the above problem is(74)
u
x
,
t
=
e
t
x
2
1
-
x
2
.Table 1 shows the maximum error between the exact solution and the numerical solution obtained by the second-order accurate implicit numerical method described in Section 2 for Example 1 at time t
=
1. Figure 1 also shows the exact solution and the numerical solution. From Table 1 and Figure 1, it can be seen that the numerical results are in good agreement with the theoretical results.Table 1
The maximum errors for the numerical method and the convergence orders withσ
=
ρ
=
1
/
1000 and τ
=
h at time T
=
1.
h
=
τ
Max. error
Error rate
1/20
5.154
e
-
2
—
1/40
1.247
e
-
2
2.047
1/80
3.021
e
-
3
2.045
1/160
7.307
e
-
4
2.048Figure 1
Exact solutions and numerical solutions withσ
=
ρ
=
1
/
1000 and h
=
τ
=
1
/
80 at T
=
1.0.Example 2.
Consider the following Riesz space distributed-order advection-dispersion equation in 2D case:(75)
∂
u
x
,
y
,
t
∂
t
=
∫
0
1
P
α
∂
α
u
x
,
y
,
t
∂
x
α
d
α
+
∫
0
1
P
α
∂
α
u
x
,
y
,
t
∂
y
α
d
α
+
∫
1
2
Q
β
∂
β
u
x
,
y
,
t
∂
x
β
d
β
+
∫
1
2
Q
β
∂
β
u
x
,
y
,
t
∂
y
β
d
β
+
f
x
,
y
,
t
,where(76)
P
α
=
-
2
Γ
5
-
α
cos
⁡
π
α
2
,
Q
β
=
-
2
Γ
5
-
β
cos
⁡
π
β
2
.Now we consider (75) in the finite domain Ω
=
[
0,1
]
×
[
0,1
] with the following initial and boundary conditions:(77)
u
0
,
y
,
t
=
u
1
,
y
,
t
=
u
x
,
0
,
t
=
u
x
,
1
,
t
=
0
,
t
∈
0,1
,
u
x
,
y
,
0
=
x
2
1
-
x
2
y
2
1
-
y
2
,
x
,
y
∈
0,1
×
0,1
,where(78)
f
x
,
y
,
t
=
e
t
x
2
1
-
x
2
y
2
1
-
y
2
-
2
e
t
y
2
1
-
y
2
2
x
2
-
2
x
ln
-
3
⁡
1
-
x
+
6
x
3
-
5
x
2
+
x
-
2
ln
-
2
⁡
1
-
x
+
6
x
2
x
3
-
2
x
2
-
x
+
1
ln
-
1
⁡
1
-
x
+
2
x
x
-
1
ln
-
3
⁡
x
+
x
6
x
2
-
13
x
+
5
ln
-
2
⁡
x
+
6
x
x
-
1
2
x
2
-
4
x
+
1
ln
-
1
⁡
x
-
2
e
t
x
2
1
-
x
2
2
y
2
-
2
y
ln
-
3
⁡
1
-
y
+
6
y
3
-
5
y
2
+
y
-
2
ln
-
2
⁡
1
-
y
+
6
y
2
y
3
-
2
y
2
-
y
+
1
ln
-
1
⁡
1
-
y
+
2
y
y
-
1
ln
-
3
⁡
y
+
y
6
y
2
-
13
y
+
5
ln
-
2
⁡
y
+
6
y
y
-
1
2
y
2
-
4
y
+
1
ln
-
1
⁡
y
-
2
e
t
y
2
1
-
y
2
-
2
x
ln
-
3
⁡
1
-
x
+
6
x
2
-
x
-
2
ln
-
2
⁡
1
-
x
+
2
-
6
x
3
+
3
x
2
+
3
x
-
1
ln
-
1
⁡
1
-
x
+
2
x
-
1
ln
-
3
⁡
x
+
6
x
2
-
11
x
+
3
ln
-
2
⁡
x
+
2
6
x
3
-
15
x
2
+
9
x
-
1
ln
-
1
⁡
x
-
2
e
t
x
2
1
-
x
2
-
2
y
ln
-
3
⁡
1
-
y
+
6
y
2
-
y
-
2
ln
-
2
⁡
1
-
y
+
2
-
6
y
3
+
3
y
2
+
3
y
-
1
ln
-
1
⁡
1
-
y
+
2
y
-
1
ln
-
3
⁡
y
+
6
y
2
-
11
y
+
3
ln
-
2
⁡
y
+
2
6
y
3
-
15
y
2
+
9
y
-
1
ln
-
1
⁡
y
.The exact solution of the above problem is(79)
u
x
,
y
,
t
=
e
t
x
2
1
-
x
2
y
2
1
-
y
2
.Table2 shows the maximum error between the exact solution and the numerical solution obtained by the second-order accurate implicit alternating direction method described in Section 4 for Example 2 at time t
=
1. Figure 2 also shows the exact solution and the numerical solution. From Table 2 and Figure 2, it can be seen that the numerical results are in good agreement with the theoretical results.Table 2
The maximum errors for the numerical method and the convergence orders withσ
=
ρ
=
1
/
1000 and τ
=
h at time T
=
1.
h
=
τ
Max. error
Error rate
1/20
3.592
e
-
2
—
1/40
8.916
e
-
3
2.010
1/80
2.204
e
-
3
2.016
1/160
5.388
e
-
4
2.032Figure 2
Exact solutions and numerical solutions withσ
=
ρ
=
1
/
1000 and τ
=
h
=
1
/
320 at T
=
1.0.


## 6. Conclusion

In this paper, we considered the Riesz space distributed-order advection-dispersion equations in 1D and 2D cases. For 1D case, we discretized the Riesz space distributed-order advection-dispersion equation (RSDO-ADE) into multiterm Riesz space fractional advection-dispersion equations (MT-RSDO-ADE), and a second-order accurate implicit numerical method is proposed using Crank-Nicholson method and a second-order accurate numerical scheme. The stability and convergence are proved. For 2D case, we proposed a new second-order accurate implicit alternating direction method for MT-RSDO-ADE; the stability and convergence of this method are also proved. Finally, numerical results are presented to support our theoretical analysis. This method may be extended to the high-dimensional time, space, and time-space distributed-order partial differential equations. These numerical methods and techniques presented are accurate and effective and can be used to simulate the corresponding physical process.



---

*Source: 590435-2015-11-25.xml*