# A Parallel Spectral Element Method for Fractional Lorenz System

Authors: Yanhui Su
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## Abstract

We provide a parallel spectral element method for the fractional Lorenz system numerically. The detailed construction and implementation of the method are presented. Thanks to the spectral accuracy of the presented method, the storage requirement due to the “global time dependence” can be considerably relaxed. Also, the parareal method combining spectral element method reduces the computing time greatly. Finally, we have tested the chaotic behaviors of fractional Lorenz system. Our numerical results are in excellent agreement with the results from other numerical methods.


---

## Body

## 1. Introduction

In 1963, Lorenz [1] attempted to set up a (much simpler and easier-to-analyze) system of differential equations that would explain some of the unpredictable behavior of the weather. The classical Lorenz system is(1)
x
′
=
σ
y
-
x
,
y
′
=
ρ
x
-
y
-
x
z
,
z
′
=
x
y
-
b
z
,where all three parameters, ρ, σ, b, are assumed to be positive, and moreover, σ
>
b
+
1. Lorenz made an important discovery that all nonequilibrium solutions tend eventually to the same complicated set, the Lorenz attractor. The behavior of the chaotic system is one of the most important subjects in modern mathematics.In [2], I. Grigorenko and E. Grigorenko generalized the Lorenz system to the fractional case by replacing the usual derivatives by the fractional one and found some interesting chaotic behaviors of fractional Lorenz system. This generation is natural. As we know, fractional derivatives play a more and more important role in modern sciences. They appear in the investigation of transport dynamics in complex systems which are characterized by the anomalous diffusion and nonexponential patterns [3]. It has been shown that many important physical and biological systems can be described by fractional differential equations (see [3–6], etc.). This has led to some investigations on the theoretical analysis and scientifical computation of the fractional differential equations in recent years (for details see [7–13] and the references therein). However, the fact that the fractional derivative is a nonlocal operator makes the results of fractional calculus in many respects far from complete.First, let us recall that, for0
<
α
<
1 and f
(
t
) defined on the interval [
a
,
b
], the (Caputo) fractional derivative is defined as(2)
D
a
t
α
f
t
=
1
Γ
1
-
α
∫
a
t
f
′
τ
t
-
τ
α
d
τ
,
t
∈
a
,
b
,where, Γ
(
·
) is the Gamma function. And for simplicity, we denote c
(
α
)
=
1
/
Γ
(
1
-
α
). Thus, the fractional generalization of the Lorenz system is(3)
D
0
t
α
x
=
σ
y
-
x
,
D
0
t
β
y
=
ρ
x
-
y
-
x
z
r
,
D
0
t
γ
z
=
x
y
-
b
z
,where 0
<
α
,
β
,
γ
≤
1, and r
≥
1 is an integer. In [2], the authors pointed out that the system can have an effective noninteger dimension Σ defined as the sum of the orders of all involved derivatives. They found that the system with Σ
<
3 can also exhibit chaotic behavior. Furthermore, a striking finding is that there is a critical value of the effective dimension Σ
c
r, under which the system undergoes a transition from chaotic daynamics to regular one.The numerical scheme of [2] is based on the linearization of system (3), on each step of integration and iterative application of the exact solution of linear fractional differential equations. Also, homotopy methods to solve the (fractional) Lorenz system were proposed in [14, 15]. It has been known that the feature of the fractional derivatives makes the design of accurate and fast methods difficult. Unlike the integer derivatives, which are local in the sense that the derivative of a function at a certain point depends only on the function in the vicinity of this point, presence of the integral in the noninteger order derivatives makes the problem global. This means that the solution at a time t
k depends on the solutions at all previous time levels t
<
t
k. The fact that all previous solutions have to be saved to compute the solution at the current time level would make the storage very expensive if low-order methods are employed for temporal discretization. Thus, it is very desirable to use high-order methods for efficient computations of the numerical solution of this kind of problem. Among the high-order methods, spectral methods are known to be very useful for their exponential rate of convergence, which is also demonstrated for solving fractional differential equations [8, 9]. Thanks to the spectral accuracy, the storage requirement due to the “global time dependence” can be considerably relaxed. However, a drawback of the spectral method is also well-known; that is, the matrix associated with the spectral method is full and the computational cost grows more quickly than that of a low-order method. Thus for long integration it is desirable to combine the spectral method with domain decomposition techniques to avoid using a single high-degree polynomial (see [16] for details). This motivates us to consider the parallel spectral element method, since, as compared to a low-order method, this method needs fewer grid points to produce highly accurate solutions, and, furthermore, the parallel method reduces the computing time greatly.This paper is organized as follows. In Section2, we construct and implement the parallel spectral element schema for the numerical solution of the fractional Lorenz system. In Section 3, we have tested the chaotic behavior of fractional Lorenz system using our numerical method. The numerical results are in excellent agreement with the results from other numerical methods. Finally, Section 4 includes some concluding remarks.


## 2. The Implementation of the Parallel Spectral Element Method

In this section, we implement the parallel spectral element method for fractional Lorenz system (3). First, we divide the interval I
=
[
0
,
T
] into M subintervals I
m
=
[
t
m
-
1
,
t
m
], m
=
1
,
…
,
M. The length of each subinterval is Δ
t
=
T
/
M, and t
m
=
m
Δ
t, m
=
0,1
,
…
,
M. According to this partition, (3) is divided into the following system of M equations:(4)
D
t
m
-
1
t
α
x
m
=
σ
y
m
-
x
m
-
c
α
∑
k
=
1
m
-
1
∫
t
k
-
1
t
k
x
k
′
τ
t
-
τ
α
d
τ
,
D
t
m
-
1
t
β
y
m
=
ρ
x
m
-
y
m
-
x
m
z
m
r
-
c
β
∑
k
=
1
m
-
1
∫
t
k
-
1
t
k
y
k
′
τ
t
-
τ
β
d
τ
,
D
t
m
-
1
t
γ
z
m
=
x
m
y
m
-
b
z
m
-
c
γ
∑
k
=
1
m
-
1
∫
t
k
-
1
t
k
z
k
′
τ
t
-
τ
γ
d
τ
,
t
∈
I
m
,with x
m
=
x
|
I
m, y
m
=
y
|
I
m, and z
m
=
z
|
I
m, m
=
1,2
,
…
,
M. Let X
m
=
{
x
m
,
y
m
,
z
m
}; then the solution of (4) can be expressed by the solvers S
m as follows:(5)
X
m
t
=
S
m
X
1
,
…
,
X
m
-
1
,
t
∈
I
m
.In most cases S
m is not realizable and needs to be approximated.The main idea of the parallel method consists of two steps. First, design two families of discrete solvers, fine solverF
m and coarse solver G
m, to approximate S
m. The fine approximate solver F
m is defined by the spectral collocation method with polynomial degree of N:(6)
U
m
t
=
F
m
U
1
,
…
,
U
m
-
1
,with the polynomial U
m
=
(
U
m
,
V
m
,
W
m
) of degree N. The coarse solver operator G
m is defined in the same way, using polynomials of a lower degree N
~
<
N. Second, set up the iteration:(7)
U
m
k
t
=
G
m
U
1
k
,
…
,
U
m
-
1
k
+
F
m
U
1
k
-
1
,
…
,
U
m
-
1
k
-
1
-
G
m
U
1
k
-
1
,
…
,
U
m
-
1
k
-
1
,
k
≥
1
,with the initial value(8)
U
m
0
t
=
G
m
U
1
0
,
…
,
U
m
-
1
0
,where subscript m refers to the element number, subscript k refers to the iteration number, and the polynomial U
m
k
=
(
U
m
k
,
V
m
k
,
W
m
k
) is the approximate solution of X
m at the time interval I
m. The parallel scheme can be seen as the predictor-corrector method, with the predictor G
m
(
U
1
k
,
…
,
U
m
-
1
k
) and the corrector C
m
=
:
F
m
(
U
1
k
-
1
,
…
,
U
m
-
1
k
-
1
)
-
G
m
(
U
1
k
-
1
,
…
,
U
m
-
1
k
-
1
). The key point of the parallel method is that the corrector C
m can be computed in parallel for 1
≤
m
≤
M. The predictor term G
m
(
U
1
k
,
…
,
U
m
-
1
k
) in (7) must be computed sequentially, but this is relatively cheap if N
~ is small compared to N. Thus, provided the iteration converges rapidly, we can use multiple processors to obtain the high-accuracy solution in a small multiple of the time needed to compute a low-accuracy solution.Now we begin to describe the fine approximate solverF
m. We first define P
N as the polynomial space of degree less than or equal to N. We denote by J
N
α
,
β
(
x
), x
∈
[
-
1,1
], the Jacobi polynomial of degree N with respect to the weight function ω
α
,
β
(
x
)
=
(
1
-
x
)
α
(
1
+
x
)
β, -
1
<
α
,
β
<
1. Let x
i
α
,
β be the points of the Gauss-Labatto-Jacobi (GJ) quadrature formula, defined by(9)
J
N
+
1
α
,
β
x
i
α
,
β
=
0
,
i
=
0
,
…
,
N
,arranged by increasing order: x
0
α
,
β
<
x
1
α
,
β
<
⋯
<
x
N
α
,
β. The associated weights of the GLJ quadrature formula are denoted by ω
i
α
,
β, 0
≤
i
≤
N. Particularly, let ξ
m
i be the Gauss-Labatto-Legendre (GLL) points on the interval I
m.Thus, the spectral collocation method for fractional Lorenz system (3) on each subinterval I
m is to find U
m
(
t
)
,
V
m
(
t
)
,
W
m
(
t
)
∈
P
N
(
I
m
), such that for all 0
≤
i
≤
N
(10)
D
t
m
-
1
t
α
U
m
ξ
m
i
+
σ
U
m
ξ
m
i
-
σ
V
m
ξ
m
i
=
-
c
α
∑
k
=
1
m
-
1
∫
t
k
-
1
t
k
U
k
′
τ
d
τ
ξ
m
i
-
τ
α
,
D
t
m
-
1
t
β
V
m
ξ
m
i
-
ρ
U
m
ξ
m
i
+
V
m
ξ
m
i
+
U
m
ξ
m
i
W
m
ξ
m
i
r
=
-
c
β
∑
k
=
1
m
-
1
∫
t
k
-
1
t
k
V
k
′
τ
d
τ
ξ
m
i
-
τ
β
,
D
t
m
-
1
t
γ
W
m
ξ
n
i
+
b
W
m
ξ
n
i
-
U
m
ξ
n
i
V
m
ξ
n
i
=
-
c
γ
∑
k
=
1
m
-
1
∫
t
k
-
1
t
k
W
k
′
τ
d
τ
ξ
m
i
-
τ
γ
.For the nonlinear problem (10), we choose to use the Lagrangian polynomials as a basis of the approximation spaces. Let h
m
i
t
∣
i
=
0,1
,
…
,
N be the Lagrangian polynomials associated with the points ξ
m
i
∣
i
=
0,1
,
…
,
N. That is, h
m
i
(
t
)
∈
P
N
(
I
m
) such that h
m
i
(
ξ
m
k
)
=
δ
i
k, where δ denotes the Kronecker function. It is seen that the set {
h
m
i
(
t
)
;
i
=
0
,
…
,
N
} forms a basis of P
N
(
I
m
):(11)
P
N
I
m
=
span
⁡
⁡
h
m
i
t
,
i
=
0
,
…
,
N
.By expressing U
m, V
m, and W
m and all nonlinear term (
W
m
(
t
)
)
r in this basis,(12)
U
m
t
=
∑
i
=
0
N
u
m
i
h
m
i
t
,
V
m
t
=
∑
i
=
0
N
v
m
i
h
m
i
t
,
W
m
t
=
∑
i
=
0
N
w
m
i
h
m
i
t
,
W
m
t
r
=
∑
i
=
0
N
w
m
i
r
h
m
i
t
.We arrive at the matrix statement of (10):(13)
A
m
1
u
m
+
σ
u
m
-
σ
v
m
=
∑
k
=
1
m
-
1
B
m
,
k
1
u
k
,
A
m
2
v
m
-
ρ
u
m
+
v
m
+
u
m
∗
w
m
r
=
∑
k
=
1
m
-
1
B
m
,
k
2
v
k
,
A
m
3
w
m
+
b
w
m
-
u
m
v
m
=
∑
k
=
1
m
-
1
B
m
,
k
3
w
k
,where u
m
=
(
u
m
i
)
(
N
+
1
), v
m
=
(
v
m
i
)
(
N
+
1
), and w
m
=
(
w
m
i
)
(
N
+
1
) are the unknown vectors. ∗ is acting on two vectors (
u
∗
v
)
i
=
u
i
∗
v
i. The matrices A
m
1, A
m
2, A
m
3, B
m
,
k
1, B
m
,
k
2, and B
m
,
k
3 are (
N
+
1
)
×
(
N
+
1
) matrices. They are computed as follows:(14)
A
m
1
i
,
j
=
c
α
∫
t
m
-
1
ξ
m
i
h
m
j
′
τ
ξ
m
i
-
τ
α
d
τ
,
B
m
,
k
1
i
,
j
=
-
c
α
∫
t
k
-
1
t
k
h
k
j
′
τ
ξ
m
i
-
τ
α
d
τ
,
A
m
2
i
,
j
=
c
β
∫
t
m
-
1
ξ
m
i
h
m
j
′
τ
ξ
m
i
-
τ
β
d
τ
,
B
m
,
k
2
i
,
j
=
-
c
β
∫
t
k
-
1
t
k
h
k
j
′
τ
ξ
m
i
-
τ
β
d
τ
,
A
m
3
i
,
j
=
c
γ
∫
t
m
-
1
ξ
m
i
h
m
j
′
τ
ξ
m
i
-
τ
γ
d
τ
,
B
m
,
k
3
i
,
j
=
-
c
γ
∫
t
k
-
1
t
k
h
k
j
′
τ
ξ
m
i
-
τ
γ
d
τ
.We are now led to compute the integrals∫
t
m
-
1
ξ
m
i
‍
(
h
m
j
′
(
τ
)
/
(
ξ
m
i
-
τ
)
α
)
d
τ and ∫
t
k
-
1
t
k
‍
h
k
j
′
(
τ
)
/
(
ξ
m
i
-
τ
)
α
d
τ one by one.For the first one, we employ Gauss-type quadrature to calculate it. It is well-known that the numerical quadrature(15)
∫
-
1
1
u
t
1
-
t
-
α
d
t
≃
∑
k
=
0
N
u
x
k
-
α
,
0
ω
k
-
α
,
0is exact for all functions u such that u
∈
P
2
N
-
1
(
-
1,1
). As a result, by the linear variable transformation(16)
τ
θ
=
ξ
m
i
-
t
m
-
1
2
θ
+
ξ
m
i
+
t
m
-
1
2
,
θ
∈
-
1,1
,it holds(17)
∫
t
m
-
1
ξ
m
i
h
m
j
′
τ
ξ
m
i
-
τ
α
d
τ
=
ξ
m
i
-
t
m
-
1
2
1
-
α
∫
-
1
1
1
-
θ
-
α
h
m
j
′
τ
θ
d
θ
=
ξ
m
i
-
t
m
-
1
2
1
-
α
∑
k
=
0
N
h
m
j
′
τ
x
k
-
α
,
0
ω
k
-
α
,
0
,since h
m
j
′
(
τ
(
θ
)
)
∈
P
2
N
-
1
(
I
) for all j
=
0,1
,
…
,
N.For the second one, it is important to point out that the Gauss-type quadrature fails ifτ approach ξ
m
i. Here we design a new technique to calculate this integral precisely. We first rewrite the integral ∫
t
k
-
1
t
k
‍
(
h
k
j
′
(
τ
)
/
(
ξ
m
i
-
τ
)
α
)
d
τ into the form(18)
∫
t
k
-
1
t
k
h
k
j
′
τ
ξ
m
i
-
τ
α
d
τ
=
2
α
t
k
-
t
k
-
1
α
∫
-
1
1
h
j
′
θ
ξ
~
m
i
-
θ
α
d
θ
,by the linear transformation(19)
τ
=
t
k
-
t
k
-
1
2
θ
+
t
k
+
t
k
-
1
2
,where h
j
′
(
θ
)
=
h
k
j
′
(
τ
(
θ
)
) is the Lagrangian polynomials based on the GLL points x
j
0,0, j
=
0,1
,
…
,
N, and ξ
~
m
i
=
(
2
ξ
m
i
-
t
k
-
t
k
-
1
)
/
(
t
k
-
t
k
-
1
). Next, we express h
j by the Legendre basis:(20)
h
j
θ
=
∑
l
=
0
N
c
j
l
L
l
θ
,where c
j
l
=
L
l
θ
,
h
j
θ
/
(
L
l
(
θ
)
,
L
l
(
θ
)
).Now we turn to compute∫
-
1
1
‍
(
L
l
′
(
θ
)
/
(
ξ
~
m
i
-
θ
)
α
)
d
θ. By employing the property(21)
L
l
p
1
=
1
2
p
C
l
p
∏
k
=
1
p
l
+
k
,
L
l
p
-
1
=
-
1
l
+
p
1
2
p
C
l
p
∏
k
=
1
p
l
+
k
,
p
=
0
,
…
,
l
,and integrating ∫
-
1
1
‍
(
L
l
′
(
θ
)
/
(
ξ
~
m
i
-
θ
)
α
)
d
θ by parts, we get(22)
∫
-
1
1
L
l
′
θ
ξ
~
m
i
-
θ
α
d
θ
=
-
∑
p
=
1
l
1
2
p
C
l
p
∏
k
=
1
p
l
+
k
k
-
α
ξ
~
m
i
-
1
p
-
α
-
-
1
l
+
p
1
2
p
C
l
p
∏
k
=
1
p
l
+
k
k
-
α
ξ
~
m
i
+
1
p
-
α
.Thus, applying (18), (20), and (22) yields(23)
B
m
,
k
1
i
,
j
=
-
c
α
2
α
t
k
-
t
k
-
1
α
∑
l
=
0
N
c
j
l
∫
-
1
1
L
l
′
θ
ξ
~
m
i
-
θ
α
d
θ
=
c
α
2
α
t
k
-
t
k
-
1
α
·
∑
l
=
0
N
c
j
l
∑
p
=
1
l
1
2
p
C
l
p
∏
k
=
1
p
l
+
k
k
-
α
ξ
~
m
i
-
1
p
-
α
-
-
1
l
+
p
1
2
p
C
l
p
∏
k
=
1
p
l
+
k
k
-
α
ξ
~
m
i
+
1
p
-
α
.Similarly, we can compute the matrices B
m
,
k
2 and B
m
,
k
3.Equation (13) is a nonlinear nonsymmetric system. We use the Jacobian-free Newton-Krylov methods to solve it. Let U
m
=
(
u
m
,
v
m
,
w
m
)
T; we rewrite (13) in the following form: for 1
≤
m
≤
M,(24)
F
1
U
m
=
A
m
1
u
m
+
σ
u
m
-
σ
v
m
-
∑
k
=
1
m
-
1
B
m
,
k
1
u
k
=
0
,
F
2
U
m
=
A
m
2
v
m
-
ρ
u
m
+
v
m
+
u
m
∗
w
m
r
-
∑
k
=
1
m
-
1
B
m
,
k
2
v
k
=
0
,
F
3
U
m
=
A
m
3
w
m
+
b
w
m
-
u
m
∗
v
m
-
∑
k
=
1
m
-
1
B
m
,
k
1
w
k
=
0
.The iteration can be described as follows:(i)
The initial guessU
m
(
0
) is given.
(ii)
Using BI-Conjugate-Gradient-stab (BICGstab) method solve the following linear system:(25)
J
F
i
U
m
k
δ
U
=
-
F
i
U
m
k
,
whereδ
U
=
(
δ
u
,
δ
v
,
δ
w
)
T, k is the iteration index, and J
F
i
(
U
(
k
)
) is the Jacobian matrix of function F
i at U
m
(
k
). We need not compute the matrix J
F
i
(
U
m
(
k
)
). In fact, only the matrix-vector product is needed in the BICGstab method. Let d
F
i
(
U
m
(
k
)
,
δ
U
) be the Fréchet derivative of the function F
i at U
m
(
k
) with respect to the direction δ
U. Consider(26)
J
F
1
U
m
k
δ
U
=
lim
ε
→
0
⁡
F
1
U
k
+
ε
δ
U
-
F
1
U
k
ε
=
A
m
1
δ
u
+
σ
δ
u
-
σ
δ
v
,
J
F
2
U
m
k
δ
U
=
lim
ε
→
0
⁡
F
2
U
m
k
+
ε
δ
U
-
F
2
U
m
k
ε
=
A
m
2
δ
u
-
ρ
δ
u
+
δ
v
+
w
m
r
∗
δ
u
+
r
u
m
∗
w
m
r
-
1
∗
δ
w
,
J
F
3
U
m
k
δ
w
=
d
F
3
U
m
k
,
δ
U
=
lim
ε
→
0
⁡
F
3
U
m
k
+
ε
δ
U
-
F
3
U
m
k
ε
=
A
m
3
δ
w
+
b
δ
w
-
u
m
∗
δ
v
-
v
m
∗
δ
u
.
(iii)
UpdateU
m
(
k
) by(27)
U
m
k
+
1
=
U
m
k
+
ω
·
δ
U
,
whereω is the correction factor, which is selected such that(28)
F
U
m
k
+
1
<
1
-
λ
ω
F
U
m
k
,
whereλ
=
1
/
10 is a scale parameter; we test ω
=
1,1
/
4,1
/
16
,
…, until (28) is satisfied.
(iv)
Stop the iteration, if(29)
δ
U
U
m
k
≤
10
-
15
.


## 3. Numerical Results

For a comparison with [2, 14], we take the initial values x
(
0
)
=
10, y
(
0
)
=
0, and z
(
0
)
=
10. In our numerical experiments, we take T
=
20, Δ
t
=
0.004, N
~
=
1, and N
=
5.First, we consider the caseρ
=
28, σ
=
10, b
=
8
/
3, and r
=
1. Note that, in this case, the equilibrium points of (3) are the origin, and(30)
Q
±
=
±
6
2
,
±
6
2
,
27
.As pointed in [2], the solution can exhibit chaotic behavior if the effective dimension Σ
=
α
+
β
+
γ is greater than the critical dimension Σ
c
r. In our numerical experiments, we take α
=
β
=
γ
=
0.99 and α
=
β
=
γ
=
0.98 and α
=
β
=
γ
=
0.97. The phase portraits of x
,
y
,
z are shown in Figures 1 and 2. From these figures, it may be concluded that, for α
=
β
=
γ
=
0.98 up to α
=
β
=
γ
=
1, the system is chaotic, while for α
=
β
=
γ
<
0.97, the system turns to be regular and the solution converges to one of the two equilibrium points Q
±. This also coincides with other numerical results.Figure 1
The phase portraits ofx-y-z for fractional Lorenz system α
=
β
=
γ
=
0.98 (left) and α
=
β
=
γ
=
0.99 (right) for the case r
=
1.
(a)
(b)Figure 2
The phase portraits ofx-y-z for fractional Lorenz system α
=
β
=
γ
=
0.97 for the case r
=
1.Secondly, we taker
=
2, while other parameters ρ, σ and b are given above. Furthermore, we take r
=
3, while other parameters ρ, σ and b are given above. The equilibrium points of (3) are the origin, and(31)
Q
±
=
±
2
2
3
,
±
2
2
3
,
3
3
.The results are similar to the case of r
=
1. However, the critical efficient dimension Σ
c
r is much smaller than the case r
=
1. From Figure 3, it may be concluded that, for α
=
β
=
γ
=
0.94 up to α
=
β
=
γ
=
1, the system is chaotic, while for α
=
β
=
γ
<
0.93, the system turns to be regular and the solution converges to one of the two equilibrium points Q
±.Figure 3
The phase portraits ofx-y-z for fractional Lorenz systems α
=
β
=
γ
=
0.93 (left) and α
=
β
=
γ
=
0.94 (right) for the case r
=
2.
(a)
(b)Furthermore, if we taker
=
3, while other parameter as given above. The equilibrium points are the origin, and(32)
Q
±
=
±
2
2
,
±
2
2
,
3
.We can see from Figure 4 that the critical efficient dimension Σ
c
r is much smaller than the case r
=
2.Figure 4
The phase portraits ofx-y-z for fractional Lorenz systems α
=
β
=
γ
=
0.91 (left) and α
=
β
=
γ
=
0.92 (right) for the case r
=
3.
(a)
(b)


## 4. Efficiency of the Method

### 4.1. Comparison with Other Methods

In Table1, we compare the spectral method with classical finite difference method and finite element method with respect to convergence rate, computational complexity, and storage.Table 1
Comparison with FDM and FEM.
Spectral method
FDM
FEM
Convergence rate
O
(
e
-
c
N
)
≥

N
α
-
r
≥

N
-
2
Complexity
O
(
N
3
)
O
(
N
3
)
O
(
N
3
)
Storage
22
O
(
1
0
9
)
O
(
1
0
5
)Spectral method is well-known for its exponential convergency, that is,e
-
c
N, where N is the freedom of the discretized space. The convergence rate is slower than N
-
2 for finite element method [7, 12] and slower than N
α
-
r
,
r
≤
4 for finite difference method [10, 17–20].For partial differential equations of integral order, the final matrix is full for spectral method while sparse for both finite difference method and finite element method. Compared with spectral method, sparse matrix is the main advantage of finite difference method and finite element method for solving integral differential equations. However, for the fractional differential equations, the final matrix for all three methods is full. The presence of the integral in the noninteger order derivatives means that the solution at a timet
k depends on the solutions at all previous time levels t
<
t
k. This fact makes the final matrix full. So the computation complexity is about O
(
N
3
) with O
(
N
2
) for matrix-vector multiple and O
(
N
) for iteration steps.Furthermore, all previous solutions have to be saved to compute the solution at the current time level which would make the storage very expensive if low-order methods are employed for temporal discretization. For example, we consider fractional ordinary differential equations:(33)
D
0
t
α
x
=
σ
y
-
x
+
g
1
t
,
D
0
t
β
y
=
ρ
x
-
y
-
x
z
r
+
g
2
t
,
D
0
t
γ
z
=
x
y
-
b
z
+
g
3
t
,with x
(
t
)
=
y
(
t
)
=
z
(
t
)
=
s
i
n
π
t, t
∈
[
0,1
]. If 1
0
-
10 precision is desired, α
=
β
=
γ
=
0.9, then 22, O
(
1
0
9
), and O
(
1
0
5
) mounts of storage are used to approximate an exact solution by spectral method, finite difference method [10] (r
=
2), and finite element method [12] (r
=
2), respectively. In view of these three advantages, spectral method is the national choice for the fractional differential equation.


### 4.2. Efficiency of the Parallel Method

The parallelism efficiency is demonstrated by a cost comparison between the parallel in time scheme and the classical sequential scheme based on the corresponding fine mesh.First, the classical sequential scheme based on the fine mesh consists of solving the problemsF
m
t
,
U
1
,
U
2
,
…
,
U
m
-
1 consecutively for 1
≤
m
≤
M. The computational complexity is equal to the sum of all the elementary operations:(34)
C
F
=
∑
m
=
1
M
C
F
n
=
O
M
N
3
.If we implement the scheme in a parallel architecture with enough processors, then the total computational time corresponds to the cost to solve a sequential set of M coarse subproblems and a fine subproblem in a single processor. If K is the number of iterations required to achieve the desired convergence of the parallel in time algorithm, then the total computational complexity is(35)
K
O
M
N
~
2
+
N
M
N
~
+
O
N
2
+
O
N
M
N
~
.Comparing (34) with (35), we obtain a speed up (i.e., the percentage with respect to the sequential scheme) close to(36)
O
K
N
~
2
N
3
+
K
N
~
N
2
+
K
N
M
-
1
.


## 4.1. Comparison with Other Methods

In Table1, we compare the spectral method with classical finite difference method and finite element method with respect to convergence rate, computational complexity, and storage.Table 1
Comparison with FDM and FEM.
Spectral method
FDM
FEM
Convergence rate
O
(
e
-
c
N
)
≥

N
α
-
r
≥

N
-
2
Complexity
O
(
N
3
)
O
(
N
3
)
O
(
N
3
)
Storage
22
O
(
1
0
9
)
O
(
1
0
5
)Spectral method is well-known for its exponential convergency, that is,e
-
c
N, where N is the freedom of the discretized space. The convergence rate is slower than N
-
2 for finite element method [7, 12] and slower than N
α
-
r
,
r
≤
4 for finite difference method [10, 17–20].For partial differential equations of integral order, the final matrix is full for spectral method while sparse for both finite difference method and finite element method. Compared with spectral method, sparse matrix is the main advantage of finite difference method and finite element method for solving integral differential equations. However, for the fractional differential equations, the final matrix for all three methods is full. The presence of the integral in the noninteger order derivatives means that the solution at a timet
k depends on the solutions at all previous time levels t
<
t
k. This fact makes the final matrix full. So the computation complexity is about O
(
N
3
) with O
(
N
2
) for matrix-vector multiple and O
(
N
) for iteration steps.Furthermore, all previous solutions have to be saved to compute the solution at the current time level which would make the storage very expensive if low-order methods are employed for temporal discretization. For example, we consider fractional ordinary differential equations:(33)
D
0
t
α
x
=
σ
y
-
x
+
g
1
t
,
D
0
t
β
y
=
ρ
x
-
y
-
x
z
r
+
g
2
t
,
D
0
t
γ
z
=
x
y
-
b
z
+
g
3
t
,with x
(
t
)
=
y
(
t
)
=
z
(
t
)
=
s
i
n
π
t, t
∈
[
0,1
]. If 1
0
-
10 precision is desired, α
=
β
=
γ
=
0.9, then 22, O
(
1
0
9
), and O
(
1
0
5
) mounts of storage are used to approximate an exact solution by spectral method, finite difference method [10] (r
=
2), and finite element method [12] (r
=
2), respectively. In view of these three advantages, spectral method is the national choice for the fractional differential equation.


## 4.2. Efficiency of the Parallel Method

The parallelism efficiency is demonstrated by a cost comparison between the parallel in time scheme and the classical sequential scheme based on the corresponding fine mesh.First, the classical sequential scheme based on the fine mesh consists of solving the problemsF
m
t
,
U
1
,
U
2
,
…
,
U
m
-
1 consecutively for 1
≤
m
≤
M. The computational complexity is equal to the sum of all the elementary operations:(34)
C
F
=
∑
m
=
1
M
C
F
n
=
O
M
N
3
.If we implement the scheme in a parallel architecture with enough processors, then the total computational time corresponds to the cost to solve a sequential set of M coarse subproblems and a fine subproblem in a single processor. If K is the number of iterations required to achieve the desired convergence of the parallel in time algorithm, then the total computational complexity is(35)
K
O
M
N
~
2
+
N
M
N
~
+
O
N
2
+
O
N
M
N
~
.Comparing (34) with (35), we obtain a speed up (i.e., the percentage with respect to the sequential scheme) close to(36)
O
K
N
~
2
N
3
+
K
N
~
N
2
+
K
N
M
-
1
.



---

*Source: 682140-2015-05-20.xml*