# An Efficient Iteration Method for Toeplitz-Plus-Band Triangular Systems Generated from Fractional Ordinary Differential Equation

Authors: Chunye Gong; Weimin Bao; Guojian Tang; Changwan Min; Jie Liu
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

It is time consuming to numerically solve fractional differential equations. The fractional ordinary differential equations may produce Toeplitz-plus-band triangular systems. An efficient iteration method for Toeplitz-plus-band triangular systems is presented withO
M
l
o
g
M computational complexity and O
M memory complexity in this paper, compared with the regular solution with O
M
2 computational complexity and O
M
2 memory complexity. M is the discrete grid points. Some methods such as matrix splitting, FFT, compress memory storage and adjustable matrix bandwidth are used in the presented solution. The experimental results show that the presented method compares well with the exact solution and is 4.25 times faster than the regular solution.


---

## Body

## 1. Introduction

Fractional differential equation (FDE) plays an important role in dynamical systems [1] and has more than 300 years of research history [2]. Many analytical solutions and numerical solutions [3–6] have been proposed for FDE, such as finite difference method [7, 8], finite element method [9], and spectral method [10, 11]. In recent times, interest in fractional ordinary differential equations (FODE) has increased [12–15]. The derivatives in the FODE are approximated by linear combinations of function values at the discrete grid points. Compared with integer ordinary differential equations, the FODE has nonlocal effect, which means a grid point may rely on the grid points far away from its position. And a grid point of the classical integer equations may only rely on its several neighboring grid points.For integer order equations, the coefficient matrices are often sparse. Because of the nonlocal property of fractional differential operators, the numerical methods for fractional diffusion equations often generate dense or even full coefficient matrices [16]. This nonlocal property makes the computation of FODE and FDE much heavier than that of the traditional integer equations. The short memory principle [17], parallel computing [18–21], fast Fourier transformation (FFT) [22, 23], multigrid method [24], and preconditioner technologies [25, 26] are used to overcome this heavy computation. Gong et al. presented many parallel algorithms for different FDEs on both traditional and heterogeneous parallel platforms [16, 18]. Diethelm [19] proposed a parallel second-order Adams-Bashforth-Moulton method for a FODE. Wang and Du [26] proposed a superfast-preconditioned iterative method for steady-state two-side space-fractional diffusion equations.The fractional ordinary differential equations may produce Toeplitz-plus-band triangular systems. Toeplitz-plus-band systems were studied by professors Chan and Ng  [27]. They considered the solutions of Hermitian Toeplitz-plus-band systems (
A
n
+
B
n
)
x
=
b, where A
n are n-by-n Toeplitz matrices and B
n are n-by-n band matrices with bandwidth independent of n. A
n and B
n are both Hermitian matrix. The authors proved that if A
n is generated by a nonnegative piecewise continuous function and B
n is positive semidefinite, then there exists a band matrix C
n, with bandwidth independent of n, such that the spectra of C
n
-
1
(
A
n
+
B
n
) are uniformly bounded by a constant independent of n. The band preconditioner was developed for Hermitian Toeplitz systems [28]. The recursive blocked algorithms were proposed for triangular systems and the recursive algorithms lead to an automatic variable blocking that has the potential of matching the memory hierarchies of today’s HPC (high performance computing) systems  [29, 30].This paper focuses on the fractional ordinary differential equation [13]:
(1)
u
′
(
t
)
+
a
(
t
)
0
D
t
α
u
(
t
)
+
b
(
t
)
u
(
t
)
=
f
(
t
)
,
u
(
0
)
=
0
,
where 0
<
α
<
1, 0
<
t
<
T
<
+
∞, a
(
t
)
>
0, and b
(
t
)
>
0. The fractional derivative is in the Caputo form [31].Definet
i
=
i
τ for 0
≤
i
≤
M, where M is a positive integer, and τ
=
T
/
M are step size. Assume u
i to be the numerical approximation to u
(
t
i
) and f
i the numerical approximation to f
(
t
i
). Using the Grünwald approximation, the finite difference scheme for (1) is shown as follows:
(2)
u
i
-
u
i
-
1
τ
+
a
i
τ
-
α
∑
k
=
0
i
‍
w
k
u
i
-
k
n
+
b
i
u
i
=
f
i
,
u
0
=
0
,
where the normalized Grünwald weight w is defined by
(3)
w
0
=
1
,
w
i
=
(
-
1
)
i
α
(
α
-
1
)
⋯
(
α
-
i
+
1
)
i
!
,
g
g
g
g
g
g
g
g
g
g
g
g
g
g
g
g
g
g
g
g
g
i
=
1,2
,
3
,
…
.Equation (2) results in a linear system of equations
(4)
A
U
=
F
,
where U
=
(
u
1
,
u
2
,
…
,
u
M
)
T and F
=
(
f
1
,
f
2
,
…
,
f
M
)
T. If u
0
≠
0, the term should be included in F. A
=
(
a
i
j
)
M
×
M is the coefficient matrix. A is defined by
(5)
a
i
,
j
=
{
0
for
i
<
j
1
τ
+
a
i
w
0
τ
α
+
b
i
for
i
=
j
-
1
τ
+
a
i
w
1
τ
α
for
i
=
j
+
1
a
i
w
i
-
j
τ
α
for
i
>
j
+
1
.


## 2. Method

### 2.1. Analysis

In a more explicit format, matrixA can be represented as
(6)
A
=
(
1
τ
+
a
1
w
0
τ
α
+
b
1
0
⋯
0
-
1
τ
+
a
2
w
1
τ
α
1
τ
+
a
2
w
0
τ
α
+
b
2
⋯
0
⋮
⋮
⋱
⋮
a
M
w
M
-
1
τ
α
a
M
w
M
-
2
τ
α
⋯
1
τ
+
a
M
w
0
τ
α
+
b
n
)
.The linear system (4) can be solved with computational complexity O
(
M
2
), shown in Algorithm 1. The output F equals U.Algorithm 1:Forward substitution for lower triangular matrix.
input: M
,
F
,
A
output: F
(1)for
i
=
1
to
M  by  1
do
(2)
f
i
←
f
i
/
a
i
,
i
(3)
for
j
=
i
+
1
to
M  by  1
do
(4)f
j
←
f
j
-
a
i
,
j
f
jFrom (6), we can see that A has some properties.(1)
Ais a low triangular, diagonal dominant matrix.
(2)
One has|
a
i
,
j
+
1
|
>
|
a
i
,
j
| for 1
≤
i
≤
M, j
<
i. This property is determined by the normalized Grünwald weight w
i and is the mathematical background of short memory principle. This property means that for grid point p, if the distance of grid point p
1 is smaller than that of grid point p
2, p
1 has more impact on p than p
2.
(3)
IfA is split into two matrices B and C, A
=
B
-
C. B is a banded matrix and the bandwidth (number of diagonals) η
>
2. Matrix C can be factorized into a product C
=
D
T. D is a diagonal matrix diag
⁡
{
a
1
,
a
2
,
…
,
a
M
}. T is a Toeplitz matrix with M
-
η nonzero diagonals on its left-bottom part.
(4)
The Toeplitz matrixT can be stored with M
-
η memory space compared with (
M
-
η
)
(
M
-
η
+
1
)
/
2 for a general low triangular matrix with order M
-
η.


### 2.2. Efficient Iteration Method

Equation (4) evolves as follows:
(7)
(
B
-
C
)
U
=
F
(8)
B
U
n
+
1
=
C
U
n
+
F
(9)
B
U
n
+
1
=
D
T
U
n
+
F
(10)
B
U
n
+
1
=
D
(
T
U
n
)
+
F
.So the linear algebra can be solved iteratively, shown in (10). Because D is a diagonal matrix, D keeps associative law and commutative law for matrix-matrix multiplication. The rate of convergence associated with (10) depends on the eigenvalues of the iteration matrix [32]:
(11)
H
=
B
-
1
D
T
=
B
-
1
C
.Assuming errore
n
+
1
=
U
n
+
1
-
U with U satisfies A
U
=
F, then
(12)
e
n
+
1
=
B
-
1
C
e
n
=
H
e
n
=
H
n
e
1
with norm ∥
*
∥ [32]:(13)
∥
e
n
+
1
∥
=
∥
G
n
e
1
∥
=
∥
G
n
∥
×
∥
e
1
∥
.So the spectral radius ofH (ρ
(
H
)) determines the asymptotic behavior of H
n. From Theorem 11.2
.
1 of [32], we can conclude that if and only if ρ
(
H
)
<
1, (10) will converge to A
-
1
F. Generally speaking, the iteration is expected to work well with small ρ
(
H
).Assume the bandwidth of matrixB is η and V
=
D
(
T
U
n
)
+
F. Solving B
U
n
+
1
=
V needs about M
η arithmetical operations. If η is near lo
g
2
M, there are about M
lo
g
2
M arithmetical operations with forward substitution. So the computational complexity of B
U
n
+
1
=
V is O
(
M
lo
g
2
M
).AssumeE
1
=
T
U
n, E
2
=
D
E
1, and V
=
E
2
+
F. The computation of E
2 and V needs M multiplications and M additions, respectively. Because T only has nonzero M
-
η diagonals on its left-bottom part, only the front M
-
η elements of U
n are effective for the multiplication T
U
n. The back η elements of E
1 are zero. So E
1
=
T
U
n can be regarded as a Toeplitz matrix vector multiplication T
1
U
1
n with order M
-
η. It is well known that Toeplitz matrix vector multiplication with order M
-
η can be finished with O
(
M
lo
g
2
M
)
=
O
(
(
M
-
η
)
lo
g
2
(
M
-
η
)
) operations [33]. The Toeplitz matrix vector multiplication T
1
U
1
n can be computed by FFTs by first embedding T
1 into a 2
(
M
-
η
)-by-2
(
M
-
η
) circulant matrix. The cost of circulant matrix vector multiplication is O
(
2
(
M
-
η
)
lo
g
2
(
2
(
M
-
η
)
)
) by using FFTs of length 2
(
M
-
η
).So the cost of each iteration of (10) is O
(
M
lo
g
2
M
). If A is a diagonal dominant matrix, we can expect (10) can converge with not too many iterations. The efficient iteration method is shown in Algorithm 2.Algorithm 2:The efficient iteration method for FODE.
input: M
,
η
,
F
,
A
output:  U
(1)ϵ
←
1
0
-
6
,
δ
=
1.0
(2)while
δ
>
ϵ
do
(3)   set v
b
1
→
M
-
η with u
a
M
-
η
→
1
(4)  V
1
←
FFT
(
v
b
)
(5)  V
2
←
FFT
(
v
a
)
(6)  v
←
V
1
⨀
V
2
(7)  u
1
→
M
←
f
1
→
M
(8)  for
i
=
1
to
M
-
η  by 1
do
(9)u
η
+
i
←
u
η
+
i
+
v
2
(
M
-
η
)
-
i
(10)
for
j
=
1
to
M  by 1
do
(11)u
j
←
u
j
d
j
(12)for
i
=
j
+
1
to
min
⁡
{
j
+
η
-
1
,
M
}  by 1
do
(13)
u
i
←
u
i
-
b
i
-
j
+
1
u
j
(14)
δ
←
max
⁡
(
|
u
1
→
M
-
u
a
1
→
M
|
)
(15)
u
a
1
→
M
←
u
a
1
→
MIn Algorithm2, u
a
1
→
M stands for the value of previous iteration and u
1
→
M stands for the current iteration. V
1
⊙
V
2 stands for v
1
i
v
2
i with 1
≤
i
≤
2
(
M
-
η
)
-
1. d
i equals the reciprocal of a
i
,
i with 1
≤
i
≤
M. b
i stands for a
i
,
1 with 1
≤
i
≤
η. v
a and v
b are 2
(
M
-
η
)
-
1 arrays. v
a
M
-
η
-
i
+
1 equals -
a
η
+
i
,
1.The value ofη can affect the performance of Algorithm 2 shown in Table 1.Table 1
Impact ofη.
Procedure
Bigη
Smallη
Iterations
Less
More
B
U
n
+
1
=
V
Slow
Fast
V
=
D
(
T
U
n
)
+
F
Fast
SlowAlgorithm2 has five features/advantages compared to Algorithm 1.(1)
Split the coefficient matrix and solve the triangular system iteratively.
(2)
Use FFT to compute matrix vector multiplication.
(3)
Precomputed
i.
(4)
Compress storage.
(5)
Adjust parameterη.


## 2.1. Analysis

In a more explicit format, matrixA can be represented as
(6)
A
=
(
1
τ
+
a
1
w
0
τ
α
+
b
1
0
⋯
0
-
1
τ
+
a
2
w
1
τ
α
1
τ
+
a
2
w
0
τ
α
+
b
2
⋯
0
⋮
⋮
⋱
⋮
a
M
w
M
-
1
τ
α
a
M
w
M
-
2
τ
α
⋯
1
τ
+
a
M
w
0
τ
α
+
b
n
)
.The linear system (4) can be solved with computational complexity O
(
M
2
), shown in Algorithm 1. The output F equals U.Algorithm 1:Forward substitution for lower triangular matrix.
input: M
,
F
,
A
output: F
(1)for
i
=
1
to
M  by  1
do
(2)
f
i
←
f
i
/
a
i
,
i
(3)
for
j
=
i
+
1
to
M  by  1
do
(4)f
j
←
f
j
-
a
i
,
j
f
jFrom (6), we can see that A has some properties.(1)
Ais a low triangular, diagonal dominant matrix.
(2)
One has|
a
i
,
j
+
1
|
>
|
a
i
,
j
| for 1
≤
i
≤
M, j
<
i. This property is determined by the normalized Grünwald weight w
i and is the mathematical background of short memory principle. This property means that for grid point p, if the distance of grid point p
1 is smaller than that of grid point p
2, p
1 has more impact on p than p
2.
(3)
IfA is split into two matrices B and C, A
=
B
-
C. B is a banded matrix and the bandwidth (number of diagonals) η
>
2. Matrix C can be factorized into a product C
=
D
T. D is a diagonal matrix diag
⁡
{
a
1
,
a
2
,
…
,
a
M
}. T is a Toeplitz matrix with M
-
η nonzero diagonals on its left-bottom part.
(4)
The Toeplitz matrixT can be stored with M
-
η memory space compared with (
M
-
η
)
(
M
-
η
+
1
)
/
2 for a general low triangular matrix with order M
-
η.


## 2.2. Efficient Iteration Method

Equation (4) evolves as follows:
(7)
(
B
-
C
)
U
=
F
(8)
B
U
n
+
1
=
C
U
n
+
F
(9)
B
U
n
+
1
=
D
T
U
n
+
F
(10)
B
U
n
+
1
=
D
(
T
U
n
)
+
F
.So the linear algebra can be solved iteratively, shown in (10). Because D is a diagonal matrix, D keeps associative law and commutative law for matrix-matrix multiplication. The rate of convergence associated with (10) depends on the eigenvalues of the iteration matrix [32]:
(11)
H
=
B
-
1
D
T
=
B
-
1
C
.Assuming errore
n
+
1
=
U
n
+
1
-
U with U satisfies A
U
=
F, then
(12)
e
n
+
1
=
B
-
1
C
e
n
=
H
e
n
=
H
n
e
1
with norm ∥
*
∥ [32]:(13)
∥
e
n
+
1
∥
=
∥
G
n
e
1
∥
=
∥
G
n
∥
×
∥
e
1
∥
.So the spectral radius ofH (ρ
(
H
)) determines the asymptotic behavior of H
n. From Theorem 11.2
.
1 of [32], we can conclude that if and only if ρ
(
H
)
<
1, (10) will converge to A
-
1
F. Generally speaking, the iteration is expected to work well with small ρ
(
H
).Assume the bandwidth of matrixB is η and V
=
D
(
T
U
n
)
+
F. Solving B
U
n
+
1
=
V needs about M
η arithmetical operations. If η is near lo
g
2
M, there are about M
lo
g
2
M arithmetical operations with forward substitution. So the computational complexity of B
U
n
+
1
=
V is O
(
M
lo
g
2
M
).AssumeE
1
=
T
U
n, E
2
=
D
E
1, and V
=
E
2
+
F. The computation of E
2 and V needs M multiplications and M additions, respectively. Because T only has nonzero M
-
η diagonals on its left-bottom part, only the front M
-
η elements of U
n are effective for the multiplication T
U
n. The back η elements of E
1 are zero. So E
1
=
T
U
n can be regarded as a Toeplitz matrix vector multiplication T
1
U
1
n with order M
-
η. It is well known that Toeplitz matrix vector multiplication with order M
-
η can be finished with O
(
M
lo
g
2
M
)
=
O
(
(
M
-
η
)
lo
g
2
(
M
-
η
)
) operations [33]. The Toeplitz matrix vector multiplication T
1
U
1
n can be computed by FFTs by first embedding T
1 into a 2
(
M
-
η
)-by-2
(
M
-
η
) circulant matrix. The cost of circulant matrix vector multiplication is O
(
2
(
M
-
η
)
lo
g
2
(
2
(
M
-
η
)
)
) by using FFTs of length 2
(
M
-
η
).So the cost of each iteration of (10) is O
(
M
lo
g
2
M
). If A is a diagonal dominant matrix, we can expect (10) can converge with not too many iterations. The efficient iteration method is shown in Algorithm 2.Algorithm 2:The efficient iteration method for FODE.
input: M
,
η
,
F
,
A
output:  U
(1)ϵ
←
1
0
-
6
,
δ
=
1.0
(2)while
δ
>
ϵ
do
(3)   set v
b
1
→
M
-
η with u
a
M
-
η
→
1
(4)  V
1
←
FFT
(
v
b
)
(5)  V
2
←
FFT
(
v
a
)
(6)  v
←
V
1
⨀
V
2
(7)  u
1
→
M
←
f
1
→
M
(8)  for
i
=
1
to
M
-
η  by 1
do
(9)u
η
+
i
←
u
η
+
i
+
v
2
(
M
-
η
)
-
i
(10)
for
j
=
1
to
M  by 1
do
(11)u
j
←
u
j
d
j
(12)for
i
=
j
+
1
to
min
⁡
{
j
+
η
-
1
,
M
}  by 1
do
(13)
u
i
←
u
i
-
b
i
-
j
+
1
u
j
(14)
δ
←
max
⁡
(
|
u
1
→
M
-
u
a
1
→
M
|
)
(15)
u
a
1
→
M
←
u
a
1
→
MIn Algorithm2, u
a
1
→
M stands for the value of previous iteration and u
1
→
M stands for the current iteration. V
1
⊙
V
2 stands for v
1
i
v
2
i with 1
≤
i
≤
2
(
M
-
η
)
-
1. d
i equals the reciprocal of a
i
,
i with 1
≤
i
≤
M. b
i stands for a
i
,
1 with 1
≤
i
≤
η. v
a and v
b are 2
(
M
-
η
)
-
1 arrays. v
a
M
-
η
-
i
+
1 equals -
a
η
+
i
,
1.The value ofη can affect the performance of Algorithm 2 shown in Table 1.Table 1
Impact ofη.
Procedure
Bigη
Smallη
Iterations
Less
More
B
U
n
+
1
=
V
Slow
Fast
V
=
D
(
T
U
n
)
+
F
Fast
SlowAlgorithm2 has five features/advantages compared to Algorithm 1.(1)
Split the coefficient matrix and solve the triangular system iteratively.
(2)
Use FFT to compute matrix vector multiplication.
(3)
Precomputed
i.
(4)
Compress storage.
(5)
Adjust parameterη.


## 3. Numerical Example

The experiment platform is a laptop with Intel(R) Core (TM) i3-3110M CPU, 2 GB main memory, and Windows 7 operating system. The CPU clock frequency is 2.40 GHz. The code is developed with MATLAB R2012a and runs on default double precision floating point operations.The following fractional (α
=
0.8) ordinary differential equation [13] was considered:
(14)
u
′
(
t
)
+
0
D
t
α
u
(
t
)
+
(
1
+
t
)
u
(
t
)
=
f
(
t
)
,
t
>
0
,
u
(
0
)
=
0
,
where f
(
t
)
=
(
14
/
Γ
(
3.8
)
)
t
1.8
+
(
5
/
2
)
t
2
+
(
5
/
Γ
(
3.8
)
)
(
1
+
t
)
t
2.8.The exact solution of (14) is
(15)
u
(
t
)
=
5
Γ
(
3.8
)
t
2.8
.The efficient iteration method of Algorithm2 compares well with the exact solution to the FODE in the test case of (14), shown in Figure 1. The τ is 1.0
/
100. The maximum absolute error is 9.78
×
1
0
-
3. The difference between the efficient iteration method and the forward substitution Algorithm 1 is only 2.37
×
1
0
-
10. The efficient iteration method and the regular forward substitution solution have no noticeable artifacts.Figure 1
Comparison of exact solution to the solution of the fast solution at timet = 1.0.The performance comparison between regular forward substitution solution of Algorithm1 and efficient iteration method of Algorithm 2 is shown in Table 2. Columns 2 and 3 of Table 2 are the runtime and the runtime is recorded in seconds. With M
=
8
×
1
0
4, the maximum speedup is 4.25. Because the speedup increases with M, the bigger M is, the higher the speedup that can be expected is. Because of the 2 GB memory limitation, the compress memory usage is also used in Algorithm 1.Table 2
Performance comparison between regular solution and the presented efficient iteration method.
M
Presented method
Regular solution
Speedup
5
×
1
0
3
0.14
0.27
1.96
1
×
1
0
4
0.41
1.10
2.67
2
×
1
0
4
1.37
4.39
3.20
4
×
1
0
4
4.51
17.74
3.93
8
×
1
0
4
14.30
60.71
4.25The impact ofη on the performance of Algorithm 2 is shown in Table 3. The runtime of the presented method varies with η. So η is a key parameter for the performance of Algorithm 2. In real fractional ordinary applications, the proper η should be chosen.Table 3
Impact ofη for M
=
4
×
1
0
4.
η
Number of iterations
Runtime
1(
⌈
log
2
⁡
M
⌉
+
1
)
93
7.95
2(
⌈
log
2
⁡
M
⌉
+
1
)
57
4.62
3(
⌈
log
2
⁡
M
⌉
+
1
)
43
5.61
4(
⌈
log
2
⁡
M
⌉
+
1
)
36
4.51
5(
⌈
log
2
⁡
M
⌉
+
1
)
31
4.93
6(
⌈
log
2
⁡
M
⌉
+
1
)
28
5.01The presented iteration method should be regarded as an iteration method to solve not only the system generated from FODE but also the more general Toeplitz-plus-band triangular systems. The technology of parallel computing is very useful, but with less mathematical background. Parallel computing is attractive for fractional differential equations [34]. As a part of future work, first, we would like to parallelize the presented solution on shared memory or distributed memory systems. Second, accelerating the presented efficient iteration method on heterogeneous architecture [35–38] should also be interesting.



---

*Source: 194249-2014-06-12.xml*