# High-Order Compact Difference Scheme for the Numerical Solution of Time Fractional Heat Equations

Authors: Ibrahim Karatay; Serife R. Bayramoglu
Journal: The Scientific World Journal
 (2014)

---

## Abstract

A high-order finite difference scheme is proposed for solving time fractional heat equations. The time fractional derivative is described in the Riemann-Liouville sense. In the proposed scheme a new second-order discretization, which is based on Crank-Nicholson method, is applied for the time fractional part and fourth-order accuracy compact approximation is applied for the second-order space derivative. The spectral stability and the Fourier stability analysis of the difference scheme are shown. Finally a detailed numerical analysis, including tables, figures, and error comparison, is given to demonstrate the theoretical results and high accuracy of the proposed scheme.


---

## Body

## 1. Introduction

In the last decades, more and more attention has been placed on the development and research of fractional differential equations, because they can describe many phenomena, physical and chemical processes more accurately than classical integer order differential equations [1–4]. And the finite difference method is an efficient tool for solving fractional partial differential equations.There are many different discretizations in time variable equipped with the compact difference scheme in spatial variable. The approximations given in [4–9] are of the order O
(
τ
μ
+
h
4
), where 1
≤
μ
<
2. Here, we propose a method for the time fractional differential heat equations with the accuracy of order O
(
τ
2
+
h
4
).In this work, we consider the following time fractional heat equation:(1)
∂
M
α
u
(
t
,
x
)
∂
t
α
=
∂
2
u
(
t
,
x
)
∂
x
2
+
f
(
t
,
x
)
,
0
<
x
<
1
,
0
<
t
<
1
,
u
(
0
,
x
)
=
r
(
x
)
,
0
≤
x
≤
1
,
u
(
t
,
0
)
=
0
,
u
(
t
,
1
)
=
0
,
0
≤
t
≤
1
,
where the term ∂
M
α
u
(
t
,
x
)
/
∂
t
α denotes α-order modifying Riemann-Liouville fractional derivative [10] given with the following formula:
(2)
∂
M
α
u
(
t
,
x
)
∂
t
α
=
{
1
Γ
(
1
-
α
)
∂
∂
t
∫
0
t
‍
u
(
s
,
x
)
-
u
(
0
,
x
)
(
t
-
s
)
α
d
s
,
if
0
<
α
<
1
,
∂
∂
t
u
(
t
,
x
)
,
if
α
=
1
,
where Γ
(
·
) is the gamma function.Remark 1.
Ifr
(
x
)
=
0, then the Riemann-Liouville and the modified Riemann-Liouville fractional derivatives are identical, since the Riemann-Liouville derivative is given by the following formula [11]:
(3)
∂
α
u
(
t
,
x
)
∂
t
α
=
{
1
Γ
(
1
-
α
)
∂
∂
t
∫
0
t
‍
u
(
s
,
x
)
(
t
-
s
)
α
d
s
,
if
0
<
α
<
1
,
∂
∂
t
u
(
t
,
x
)
,
if
α
=
1
.
If r
(
x
) is nonzero, then there are some problems about the existence of the solutions for the heat equation (1). To rectify the situation two main approaches can be used; the modified Riemann-Liouville fractional derivative can be used [10] or the initial condition should be modified [12]. We chose the first approach in our work.


## 2. Discretization of the Problem

In this section we introduce the basic ideas for the numerical solution of the time fractional heat equation (1) by compact finite difference scheme.For some positive integersM and N, the grid sizes in space and time for the finite difference algorithm are defined by h
=
1
/
M and τ
=
1
/
N, respectively. The grid points in the space interval [
0,1
] are the numbers x
j
=
j
h, j
=
0,1
,
2
,
…
,
M, and the grid points in the time interval [
0,1
] are labeled t
k
=
k
τ, k
=
0,1
,
2
,
…
,
N. The values of the functions U and f at the grid points are denoted by U
j
k
=
U
(
t
k
,
x
j
) and f
j
k
=
f
(
t
k
,
x
j
), respectively.As in the classical Crank-Nicholson difference scheme, we use the approximation [13] to the fractional derivative ∂
α
U
(
t
,
x
)
/
∂
t
α at (
t
k
+
1
/
2
,
x
j
), and then
(4)
∂
α
U
(
t
k
+
1
/
2
,
x
j
)
∂
t
α
=
∂
∂
t
H
(
t
k
+
1
/
2
,
x
j
)
=
H
(
t
k
+
1
,
x
j
)
-
H
(
t
k
,
x
j
)
τ
+
O
(
τ
2
)
=
w
0
U
j
k
+
1
+
∑
r
=
1
k
(
w
r
-
w
r
-
1
)
U
j
k
+
1
-
r
-
w
k
U
j
0
+
O
(
τ
2
)
,
where H
(
t
,
x
)
=
(
1
/
Γ
(
1
-
α
)
)
∫
0
t
‍
(
(
u
(
s
,
x
)
-
u
(
0
,
x
)
)
/
(
t
-
s
)
α
)
d
s, w
0
=
b
0
-
a
0, and w
k
=
a
k
-
1
-
a
k
+
(
k
+
1
)
b
k
-
(
k
-
1
)
b
k
-
1, for 1
≤
k
≤
N.Definition 2.
Define the average operatorκ
:
ϑ
→
ϑ as follows:
(5)
(
κ
g
)
j
=
{
g
0
,
j
=
0
,
1
12
(
g
j
-
1
+
10
g
j
+
g
j
+
1
)
,
1
≤
j
≤
M
-
1
,
g
M
,
j
=
M
,
where g
=
(
g
0
,
g
1
,
…
,
g
M
) is a grid function and ϑ is the space of the grid functions.Lemma 3.
Supposep
(
x
)
∈
C
6
[
0,1
], then
(6)
κ
∂
2
p
(
x
i
)
∂
x
2
=
1
h
2
[
p
(
x
i
-
1
)
-
2
p
(
x
i
)
+
p
(
x
i
+
1
)
]
+
O
(
h
4
)
,
1
≤
i
≤
M
-
1
.Proof (see [14]).
We use the Taylor expansion of each term aboutx
i, and then we obtain the following truncation error for any i, where 1
≤
i
≤
M
-
1:
(7)
κ
∂
2
p
(
x
i
)
∂
x
2
-
1
h
2
[
p
(
x
i
-
1
)
-
2
p
(
x
i
)
+
p
(
x
i
+
1
)
]
=
1
12
[
p
′′
(
x
i
+
1
)
+
10
p
′′
(
x
i
)
+
p
′′
(
x
i
-
1
)
]
-
1
h
2
[
p
(
x
i
+
1
)
-
2
p
(
x
i
)
+
p
(
x
i
-
1
)
]
=
h
4
360
∫
0
1
[
p
(
6
)
(
x
i
+
s
h
)
+
p
(
6
)
(
x
i
-
s
h
)
]
×
(
1
-
s
)
3
[
5
-
3
(
1
-
s
)
2
]
d
s
=
h
4
240
p
(
6
)
(
ξ
i
)
,
ξ
i
∈
(
x
i
-
1
,
x
i
+
1
)
.


## 3. Compact Finite Difference Scheme

If we apply the operatorκ to both side of (1),
(8)
κ
∂
M
α
u
(
t
,
x
)
∂
t
α

=

κ
∂
2
u
(
t
,
x
)
∂
x
2

+

κ
f
(
t
,
x
)
,
0
<
x
<
1
,
0
<
t
<
1
,
u
(
0
,
x
)
=
r
(
x
)
,
0
≤
x
≤
1
,
u
(
t
,
0
)
=
0
,
u
(
t
,
1
)
=
0
,
0
≤
t
≤
1
,
and use the approximation (4), then we obtain the following difference scheme which is accurate of order O
(
τ
2
+
h
4
);
(9)
1
12
[
w
0
U
j
-
1
k
+
1
+
∑
r
=
1
k
(
w
r
-
w
r
-
1
)
U
j
-
1
k
+
1
-
r
-
w
k
U
j
-
1
0
]
+
10
12
[
w
0
U
j
k
+
1
+
∑
r
=
1
k
(
w
r
-
w
r
-
1
)
U
j
k
+
1
-
r
-
w
k
U
j
0
]
+
1
12
[
w
0
U
j
+
1
k
+
1
+
∑
r
=
1
k
(
w
r
-
w
r
-
1
)
U
j
+
1
k
+
1
-
r
-
w
k
U
j
+
1
0
]
=
[
U
j
+
1
k
+
1
-
2
U
j
k
+
1
+
U
j
-
1
k
+
1
2
h
2
+
U
j
+
1
k
-
2
U
j
k
+
U
j
-
1
k
2
h
2
]
+
1
12
[
f
(
t
k
+
τ
2
,
x
j
-
1
)
+
10
f
(
t
k
+
τ
2
,
x
j
)
+
f
(
t
k
+
τ
2
,
x
j
+
1
)
]
,
0
≤
k
≤
N
-
1
,
1
≤
j
≤
M
-
1
,
U
j
0
=
r
(
x
j
)
,
1
≤
j
≤
M
-
1
,
U
0
k
=
0
,
U
M
k
=
0
,
0
≤
k
≤
N
.
The difference scheme above can be written in matrix form as follows:
(10)
A
U
j
+
1
+
B
U
j
+
A
U
j
-
1
=
φ
j
,
where φ
j
=
[
φ
j
0
,
φ
j
1
,
φ
j
2
,
…
,
φ
j
N
]
T, φ
j
0
=
r
(
x
j
), φ
j
k
=
f
(
t
k
+
1
/
2
,
x
j
), 1
≤
k
≤
N, 1
≤
j
≤
M, and U
j
=
[
U
j
0
,
U
j
1
,
U
j
2
,
…
,
U
j
N
]
T.HereA
(
N
+
1
)
×
(
N
+
1
) and B
(
N
+
1
)
×
(
N
+
1
) are the matrices of the form:(11)
A
=
[
0
-
w
0
12
-
1
2
h
2
w
0
12
-
1
2
h
2
-
w
1
12
w
1
-
w
0
12
-
1
2
h
2
w
0
12
-
1
2
h
2
-
w
2
12
w
2
-
w
1
12
w
1
-
w
0
12
-
1
2
h
2
w
0
12
-
1
2
h
2
⋱
⋱
⋱
-
w
N
-
1
12
w
N
-
1
-
w
N
-
2
12
⋯
w
2
-
w
1
12
w
1
-
w
0
12
-
1
2
h
2
w
0
12
-
1
2
h
2
]
,
B
=
[
1
-
10
w
0
12
+
1
h
2
10
w
0
12
+
1
h
2
-
10
w
1
12
10
(
w
1
-
w
0
)
12
+
1
h
2
10
w
0
12
+
1
h
2
-
10
w
2
12
10
(
w
2
-
w
1
)
12
10
(
w
1
-
w
0
)
12
+
1
h
2
10
w
0
12
+
1
h
2
⋱
⋱
⋱
-
10
w
N
-
1
12
10
(
w
N
-
1
-
w
N
-
2
)
12
⋯
10
(
w
2
-
w
1
)
12
10
(
w
1
-
w
0
)
12
+
1
h
2
10
w
0
12
+
1
h
2
]
.We note that the unspecified entries are zeros at the matrices above.Using the idea of the modified Gauss-Elimination method, we can convert (10) into the following form:
(12)
U
j
=
α
j
+
1
U
j
+
1
+
β
j
+
1
,
j
=
M
-
1
,
…
,
2,1
,
0
.This way, the two-step form of difference schemes in (10) is transformed to one-step method as in (12).Now, we need to determine the matricesα
j
+
1 and β
j
+
1 satisfying the last equality. Since U
0
=
α
1
U
1
+
β
1
=
0, we can select α
1
=
O
(
N
+
1
)
×
(
N
+
1
) and β
1
=
O
(
N
+
1
)
×
1. Combining the equalities U
j
=
α
j
+
1
U
j
+
1
+
β
j
+
1, and U
j
-
1
=
α
j
U
j
+
β
j and the matrix equation (10), we have
(13)
(
A
+
B
α
j
+
1
+
A
α
j
α
j
+
1
)
U
j
+
1
+
(
B
β
j
+
1
+
A
α
j
β
j
+
1
+
A
β
j
)
=
φ
j
.
Then, we write
(14)
A
+
B
α
j
+
1
+
A
α
j
α
j
+
1
=
0
,
B
β
j
+
1
+
A
α
j
β
j
+
1
+
A
β
j
=
φ
j
,
where 1
≤
j
≤
M
-
1.So, we obtain the following pair of formulas [15]:
(15)
α
j
+
1
=
-
(
B
+
A
α
j
)
-
1
A
,
β
j
+
1
=
(
B
+
A
α
j
)
-
1
(
φ
j
-
A
β
j
)
,
where 1
≤
j
≤
M
-
1.


## 4. The Spectral Stability of the Method

The spectral stability analysis is done by using the analysis of the eigenvalues of the iteration matrixα
j (1
≤
j
≤
M) of the scheme (12).Letρ
(
A
) denote the spectral radius of a matrix A, that is, the maximum of the absolute value of the eigenvalues of the matrix A.We will prove thatρ
(
α
j
)
<
1, (
1
≤
j
≤
M
), by inductionsinceα
1 is a zero matrix ρ
(
α
1
)
=
0
<
1.Moreover,α
2
=
-
B
-
1
A, and α
2 is a lower triangular matrix of the following form:(16)
α
2
=
[
0
*
-
w
0
/
12
+
1
/
2
h
2
10
w
0
/
12
+
1
/
h
2
⋱
*
*
-
w
0
/
12
+
1
/
2
h
2
10
w
0
/
12
+
1
/
h
2
]
(
N
+
1
)
×
(
N
+
1
)
.Thereforeρ
(
α
2
)
=
ρ
(
-
B
-
1
A
)
=
|
(
1
/
2
h
2
-
w
0
/
12
)
/
(
1
/
h
2
+
10
w
0
/
12
)
| < |
(
1
/
2
h
2
-
w
0
/
12
)
/
(
1
/
2
h
2
+
9
w
0
/
12
+
1
/
2
h
2
+
w
0
/
12
)
|
<
|
(
1
/
2
h
2
-
w
0
/
12
)
/
(
1
/
2
h
2
+
w
0
/
12
)
|.Since(17)
w
0
=
b
0
-
a
0
=
τ
-
α
(
1
-
α
)
Γ
(
1
-
α
)
-
τ
-
α
(
2
-
α
)
Γ
(
1
-
α
)
=
τ
-
α
Γ
(
3
-
α
)
>
0
and 1
/
h
2
>
0, we can write ρ
(
α
2
)
<
1.Now, assumeρ
(
α
j
)
<
1. After some calculations we find that(18)
α
j
+
1
=
[
-
(
B
+
A
α
j
)
-
1
A
]
(
N
+
1
)
×
(
N
+
1
)
=
[
0
*
-
(
w
0
/
12
-
1
/
2
h
2
)
B
2,2
+
A
2,2
α
j
2,2
*
*
⋱
*
*
-
(
w
0
/
12
-
1
/
2
h
2
)
B
N
+
1
,
N
+
1
+
A
N
+
1
,
N
+
1
α
j
N
+
1
,
N
+
1
]
,and we already know that A
i
,
i
=
w
0
/
12
-
1
/
2
h
2, B
i
,
i
=
10
w
0
/
12
+
1
/
h
2 and α
j
i
,
i
=
λ
j, nonzero eigenvalue of α
j for 2
≤
i
≤
N
+
1. 						(19)
ρ
(
α
j
+
1
)
=
|
1
/
2
h
2
-
w
0
/
12
1
/
h
2
+
10
w
0
/
12
+
(
w
0
/
12
-
1
/
2
h
2
)
λ
i
|
.
If 1
/
2
h
2
-
w
0
/
12
≥
0, then we have two subcases. (a)
If0
<
λ
j
<
1,
(20)
ρ
(
α
j
)
=
|
1
/
2
h
2
-
w
0
/
12
1
/
h
2

+

10
w
0
/
12

+

(
w
0
/
12

-

1
/
2
h
2
)
λ
j
|

=
1
/
2
h
2
-
w
0
/
12
1
/
h
2
+
10
w
0
/
12
-
(
1
/
2
h
2
-
w
0
/
12
)
λ
j
≤
1
/
2
h
2
-
w
0
/
12
1
/
h
2
+
10
w
0
/
12
-
(
1
/
2
h
2
-
w
0
/
12
)
=
1
/
2
h
2
-
w
0
/
12
1
/
2
h
2
+
11
w
0
/
12
≤
1
.
(b)
If-
1
<
λ
j
<
0,
(21)
ρ
(
α
j
+
1
)
=
|
1
/
2
h
2

-

w
0
/
12
1
/
h
2

+

10
w
0
/
12

+

(
w
0
/
12

-

1
/
2
h
2
)
λ
j
|


=
1
/
2
h
2

-

w
0
/
12
1
/
h
2

+

10
w
0
/
12

-

(
1
/
2
h
2

-

w
0
/
12
)
λ
j


≤
1
/
2
h
2
-
w
0
/
12
1
/
h
2
+
10
w
0
/
12
=
1
/
2
h
2
-
w
0
/
12
1
/
2
h
2
+
w
0
/
12
+
1
/
2
h
2
+
9
w
0
/
12
≤
1
/
2
h
2
-
w
0
/
12
1
/
2
h
2
+
w
0
/
12
<
1
.If 1
/
2
h
2
-
w
0
/
12
<
0, then we have two subcases.(a)
If0
<
λ
j
<
1,
(22)
ρ
(
α
j
+
1
)
=
|

1
/
2
h
2

-

w
0
/
12

1
/
h
2

+

10
w
0
/
12

+

(
w
0
/
12

-

1
/
2
h
2
)
λ
j
|


=

-
1
/
2
h
2

+

w
0
/
12

1
/
h
2

+

10
w
0
/
12

+

(
w
0
/
12

-

1
/
2
h
2
)
λ
j


≤
w
0
/
12
-
1
/
2
h
2
w
0
/
12
+
1
/
2
h
2
+
1
/
2
h
2
+
9
w
0
/
12
=
w
0
/
12
-
1
/
2
h
2
w
0
/
12
+
1
/
2
h
2
≤
1
.
(b)
If-
1
<
λ
j
<
0,
(23)
ρ
(
α
j
+
1
)
=
|

1
/
2
h
2

-

w
0
/
12
1
/
h
2

+

10
w
0
/
12

+

(
w
0
/
12

-

1
/
2
h
2
)
λ
j
|



=

-
1
/
2
h
2

+

w
0
/
12
1
/
h
2

+

10
w
0
/
12

+

(
w
0
/
12

-

1
/
2
h
2
)
λ
j



≤
w
0
/
12
-
1
/
2
h
2
1
/
h
2

+

10
w
0
/
12

-

(
w
0
/
12

-

1
/
2
h
2
)

=
w
0
/
12
-
1
/
2
h
2
9
w
0
/
12
+
3
/
2
h
2

=
w
0
/
12

-

1
/
2
h
2
w
0
/
12

+

1
/
2
h
2

+

(
8
w
0
/
12

+

1
/
h
2
)

≤
w
0
/
12
-
1
/
2
h
2
w
0
/
12
+
1
/
2
h
2
<
1
.So, we have proved that wheneverρ
(
α
j
)
<
1 then it follows that ρ
(
α
j
+
1
)
<
1. So ρ
(
α
j
)
<
1 for any j, where 1
≤
j
≤
M.Remark 4 (see [16]).
It is well known that for anyA
∈
ℝ
N
×
N, A
m
→
0 as m
→
∞ if and only if ρ
(
A
)
≤
1. We note that if A is normal, then ∥
A
∥
=
ρ
(
A
) but when the matrix A is not normal the spectral radius gives no indication of the magnitude of the roundoff error for finite M. In this case a condition of the form ρ
(
A
)
≤
1 guarantees eventual decay of the errors, but does not control the intermediate growth of the errors. Then, it is easy to understand that ρ
(
A
)
≤
1 is a necessary condition for stability but not always sufficient.


## 5. The Fourier Stability of the Method

We analyze the stability of the difference scheme by a Fourier analysis. LetU
j
k be the approximate solution and define ρ
j
k
=
u
j
k
-
U
j
k, k
=
0,1
,
…
,
N
-
1, j
=
1
,
…
,
M
-
1. Then, we write ρ
j
k
=
d
k
e
i
j
h
β and obtain the following roundoff error equation for (9):
(24)
1
12
[
w
0
ρ
j
-
1
k
+
1
+
∑
r
=
1
k
(
w
r
-
w
r
-
1
)
ρ
j
-
1
k
+
1
-
r
-
w
k
ρ
j
-
1
0
]

+
10
12
[
w
0
ρ
j
k
+
1

+

∑
r
=
1
k
(
w
r

-

w
r
-
1
)
ρ
j
k
+
1
-
r

-

w
k
ρ
j
0
]

+
1
12
[
∑
r
=
1
k
(
w
r
-
w
r
-
1
)
w
0
ρ
j
+
1
k
+
1
+
∑
r
=
1
k
(
w
r
-
w
r
-
1
)
ρ
j
+
1
k
+
1
-
r
-
w
k
ρ
j
+
1
0
]
=
[
ρ
j
+
1
k
+
1
-
2
ρ
j
k
+
1
+
ρ
j
-
1
k
+
1
2
h
2
+
ρ
j
+
1
k
-
2
ρ
j
k
+
ρ
j
-
1
k
2
h
2
]
,
0
≤
k
≤
N
-
1
,
1
≤
j
≤
M
-
1
,
ρ
0
k
=
ρ
M
k
=
0
.We now define the grid functions:(25)
ρ
k
(
x
)
=
{
ρ
j
k
,
when
x
j
-
h
/
2
<
x
<
x
j
+
h
/
2
0
,
when
0
≤
x
≤
h
2
or
L
-
h
2
<
x
≤
L
,
and then ρ
k
(
x
) can be expanded in a Fourier series as follows:
(26)
ρ
k
(
x
)
=
∑
l
=
-
∞
∞
d
k
(
l
)
e
i
2
π
l
x
/
L
,
k
=
1,2
,
…
,
N
,
where d
k
(
l
)
=
1
/
L
∫
0
L
ρ
k
(
x
)
e
-
i
2
π
l
x
/
L and we introduce the following norm:
(27)
∥
ρ
k
∥
2
=
(
∑
j
=
1
M
-
1
h
|
ρ
j
k
|
2
)
1
/
2
=
[
∫
0
L
|
ρ
k
(
x
)
|
2
d
x
]
1
/
2
,
and applying the Parseval equality ∫
0
L
|
ρ
k
(
x
)
|
2
d
x
=
∑
l
=
-
∞
∞
|
d
k
(
l
)
|
2, we obtain
(28)
∥
ρ
k
∥
2
2
=
∑
l
=
-
∞
∞
|
d
k
(
l
)
|
2
.Based on the above analysis we can suppose that the solution of (24) has the following form ρ
j
k
=
d
k
e
i
j
h
β, where β
=
2
π
l
/
L and L
=
1. Substituting the above expression into (24), we obtain
(29)
1
12
[
∑
r
=
1
k
(
w
r
-
w
r
-
1
)
w
0
d
k
+
1
e
i
(
j
-
1
)
h
β
+
∑
r
=
1
k
(
w
r
-
w
r
-
1
)
d
k
+
1
-
r
e
i
(
j
-
1
)
h
β
-
w
k
d
0
e
i
(
j
-
1
)
h
β
]
+
10
12
[
∑
r
=
1
k
(
w
r
-
w
r
-
1
)
w
0
d
k
+
1
e
i
j
h
β
+
∑
r
=
1
k
(
w
r
-
w
r
-
1
)
d
k
+
1
-
r
e
i
j
h
β
-
w
k
d
0
e
i
j
h
β
]
+
1
12
[
∑
r
=
1
k
(
w
r
-
w
r
-
1
)
w
0
d
k
+
1
e
i
(
j
+
1
)
h
β
+
∑
r
=
1
k
(
w
r
-
w
r
-
1
)
d
k
+
1
-
r
e
i
(
j
+
1
)
h
β
-
w
k
d
0
e
i
(
j
+
1
)
h
β
∑
r
=
1
k
(
w
r
-
w
r
-
1
)
]
=
[
d
k
+
1
e
i
(
j
+
1
)
h
β
-
2
d
k
+
1
e
i
(
j
)
h
β
+
d
k
+
1
e
i
(
j
-
1
)
h
β
2
h
2
+
d
k
e
i
(
j
+
1
)
h
β
-
2
d
k
e
i
(
j
)
h
β
+
d
k
e
i
(
j
-
1
)
h
β
2
h
2
]
,
0
≤
k
≤
N
-
1
,
1
≤
j
≤
M
-
1
,
ρ
0
k
=
ρ
M
k
=
0
.After simplifications, we write(30)
d
k
+
1
(
w
0
(
2
cos
⁡
(
β
h
)
+
10
)
12
+
1
h
2
(
1
-
cos
⁡
(
β
h
)
)
)
=
d
k
(
1
h
2
(
cos
⁡
(
β
h
)
-
1
)
)
+
(
2
cos
⁡
(
β
h
)
+
10
)
12
×
[
∑
r
=
1
k
(
w
r
-
1
-
w
r
)
d
k
+
1
-
r
]
+
(
w
k
(
2
cos
⁡
(
β
h
)
+
10
)
12
)
d
0
.Theorem 5.
|
d
k
|
≤
|
d
0
| for k
=
1,2
,
…
,
N, where d
k is the solution of (30).Proof.
We will use mathematical induction for the proof.
We start withk
=
0, and then
(31)
d
1
(
w
0
(
2
cos
⁡
(
β
h
)
+
10
)
12
+
1
h
2
(
1
-
cos
⁡
(
β
h
)
)
)
=
d
0
(
1
h
2
(
cos
⁡
(
β
h
)
-
1
)
)
+
(
w
0
(
2
cos
⁡
(
β
h
)
+
10
)
12
)
d
0
.
Then
(32)
|
d
1
|
=
|
(
w
0
(
2
cos
⁡
(
β
h
)

+

10
)
/
12
)

-

(
1
/
h
2
)

(
1
-
cos
⁡
(
β
h
)
)

(
w
0
(
2
cos
⁡
⁡
(
β
h
)
+
10
)
/
12
)
+
(
1
/
h
2
)
(
1
-
cos
⁡
(
β
h
)
)
|
×
|
d
0
|
≤
|
d
0
|
,
and therefore |
d
1
|
≤
|
d
0
|.
Now, assume that|
d
n
|
≤
|
d
0
|, n
=
1,2
,
…
,
k. We need to prove that n
=
k
+
1. Indeed,
(33)
|
d
k
+
1
|
(
w
0
(
2
cos
⁡
(
β
h
)
+
10
)
12
+
1
h
2
(
1
-
cos
⁡
(
β
h
)
)
)
≤
|
d
k
|
|
1
h
2
(
cos
⁡
(
β
h
)
-
1
)
|
+
|
(
2
cos
⁡
⁡
(
β
h
)
+
10
)
12
|
∑
r
=
1
k
|
w
r
-
1
-
w
r
|
·
|
d
k
+
1
-
r
|
+
|
w
k
(
2
cos
⁡
⁡
(
β
h
)
+
10
)
12
|
|
d
0
|
,
|
d
k
+
1
|
|
w
0
(
2
cos
⁡
(
β
h
)
+
10
)
12
+
1
h
2
(
1
-
cos
⁡
(
β
h
)
)
|
≤
|
d
0
|
|
1
h
2
(
cos
⁡
(
β
h
)
-
1
)
|
+
|
(
2
cos
⁡
(
β
h
)
+
10
)
12
|
∑
r
=
1
k
|
w
r
-
1
-
w
r
|
·
|
d
0
|
+
|
w
k
(
2
cos
⁡
⁡
(
β
h
)
+
10
)
12
|
|
d
0
|
,
|
d
k
+
1
|
≤
(
|
w
k
(
2
cos
⁡
⁡
(
β
h
)
+
10
)
12
|
|
1
h
2
(
cos
⁡
⁡
(
β
h
)
-
1
)
|
+
|
(
2
cos
⁡
⁡
(
β
h
)
+
10
)
12
(
w
0
-
w
k
)
|
+
|
w
k
(
2
cos
⁡
⁡
(
β
h
)
+
10
)
12
|
)
×
(
|
w
0
12
(
2
cos
⁡
⁡
(
β
h
)
+
10
)
+
1
h
2
(
1
-
cos
⁡
⁡
(
β
h
)
)
|
)
-
1
×
|
d
0
|
,
|
d
k
+
1
|
≤
(
-
1
h
2
(
cos
⁡
(
β
h
)
-
1
)
+
(
2
cos
⁡
⁡
(
β
h
)
+
10
)
12
(
w
0
-
w
k
)
+
w
k
(
2
cos
⁡
⁡
(
β
h
)
+
10
)
12
)
×
(
w
0
12
(
2
cos
⁡
⁡
(
β
h
)
+
10
)
+
1
h
2
(
1
-
cos
⁡
(
β
h
)
)
)
-
1
×
|
d
0
|
,
|
d
k
+
1
|
≤
(
w
0
12
(
2
cos
⁡
⁡
(
β
h
)
+
10
)
+
1
h
2
(
1
-
cos
⁡
(
β
h
)
)
)
×
(
w
0
12
(
2
cos
⁡
(
β
h
)
+
10
)
+
1
h
2
(
1
-
cos
⁡
(
β
h
)
)
)
-
1
×
|
d
0
|
=
|
d
0
|
.Theorem 6.
The finite difference scheme (9) is unconditionally stable.Proof.
Using the last theorem and Parseval equality, we obtain(34)
∥
ρ
k
∥
2
≤
∥
ρ
0
∥
2
,
k
=
1,2
,
…
,
N
,
which means the proposed difference scheme is unconditionally stable.


## 6. Numerical Analysis

Example 7.
One has(35)
∂
α
u
(
t
,
x
)
∂
t
α
=
∂
2
u
(
t
,
x
)
∂
x
2
+
4
!
Γ
(
5
-
α
)
t
7
/
2
x
5
(
1
-
x
)
+
10
(
t
4
+
1
)
x
4
-
20
(
t
4
+
1
)
(
1
-
x
)
x
3
,
(
0
<
x
<
1
,
0
<
t
<
1
)
,
u
(
0
,
x
)
=
0
,
0
≤
x
≤
1
,
u
(
t
,
0
)
=
0
,
u
(
t
,
1
)
=
0
,
0
≤
t
≤
1
.Exact solution of this problem isU
(
t
,
x
)
=
(
t
4
+
1
)
x
5
(
1
-
x
). The solution by the proposed scheme is given in Figure 1. The errors when solving this problem are listed in Table 1 for various values of time and space nodes.Table 1
Error table for Example7.
M
N
α = 0.3
α = 0.5
α = 0.8
Error
Error
Error
8
8
0.0014442
0.0014297
0.0014076
16
32
0.0000910
0.0000908
0.0000898
32
128
0.0000058
0.0000058
0.0000058(a) The approximate solutions of Example7 by the proposed method when N
=
32, M
=
32 and α
=
0.5. (b) The errors for some values of M and N when t
=
1 and α
=
0.5.
(a)
(b)The errors in Table1 are calculated by the following formula:
(36)
max
⁡
0
≤
n
≤
M
0
≤
k
≤
N
|
u
(
t
k
,
x
n
)
-
U
n
k
|
.It can be concluded from Table1 and Figure 1 above that when the time-step size is reduced by a factor of 1/4 and the spatial step size is reduced by a factor of 1/2, then the error decreases by about 1/16. The numerical results support the claim about the order of the convergence.


## 7. Conclusion

In this work, the compact difference scheme was successfully applied to solve the time fractional heat equations. The second order approximation for the Riemann-Liouville fractional derivative is equipped with the higher order compact difference schemes. The Fourier analysis and the spectral stability method are used to show that the proposed scheme is unconditionally stable. Numerical results are in good agreement with the theoretical results.



---

*Source: 642989-2014-02-13.xml*