# Numerical Algorithm to Solve a Class of Variable Order Fractional Integral-Differential Equation Based on Chebyshev Polynomials

Authors: Kangwen Sun; Ming Zhu
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The purpose of this paper is to study the Chebyshev polynomials for the solution of a class of variable order fractional integral-differential equation. The properties of Chebyshev polynomials together with the four kinds of operational matrixes of Chebyshev polynomials are used to reduce the problem to the solution of a system of algebraic equations. By solving the algebraic equations, the numerical solutions are acquired. Further some numerical examples are shown to illustrate the accuracy and reliability of the proposed approach and the results have been compared with the exact solution.


---

## Body

## 1. Introduction

Fractional calculus has attracted increasing attention for decades since it plays a vital role in different disciplines of science and engineering [1–3]. Compared with integer order differential equation, fractional differential equation has the advantage that it can better describe some natural physics processes and dynamic system processes, because the fractional order differential operators are nonlocal operators. Many physics, chemistry, and engineering systems can be elegantly modeled with the help of the FDEs, such as dielectric polarization, viscoelastic systems, control theory, chaotic behavior, and electrolyte-electrolyte polarization [4–6]. Since its tremendous applications in several disciplines, considerable attention has been given to the exact and the numerical solutions of fractional differential equations and fractional integral equations. Even numerical approximation of fractional differentiation of rough functions is not easy as it is an ill-posed problem.Other than modeling aspects of these differential equations, the solution techniques and their reliability are rather more significant. In order to obtain the goal of highly accurate and reliable solutions, several methods have been proposed to solve the fractional order differential and fractional order integral equations. The most commonly used methods are Variational Iteration Method [7], Adomian Decomposition Method [8, 9], Generalized Differential Transform Method [10, 11], and Wavelet Method [12, 13].Recently, more and more physicists and mathematicians are finding that numerous important dynamical problems exhibit fractional order behavior which can vary with space and time. This fact indicates that variable order calculus provides an effective mathematical framework for the description of complex dynamical problems. The concept of a variable order operator is a much more recent development, which is a new orientation in engineering. Many researchers have proposed different definitions of variable order differential operators, each of these with a specific meaning to get desired goals. The variable order operator definitions recently proposed in the engineering include the Riemann-Liouville definition, Marchaud definition, Grünwald definition, Caputo definition, and Coimbra definition [14, 15].In this paper, the main objective is to introduce the Chebyshev polynomials method to solve the variable order fractional integral-differential equation. The method is based on reducing the equation to a system of algebraic equations by expanding the solution as Chebyshev polynomials with unknown coefficients. The main characteristic of an operational method is to convert the integral-differential equation into an algebraic one. It not only simplifies the problem but also speeds up the computation.Our study focuses on a class of variable order fractional integral-differential equation as follows:(1)
D
α
t
u
x
,
t
g
x
,
t
+
∂
u
x
,
t
∂
t
+
∫
0
t
u
x
,
T
d
T
+
∫
0
t
u
x
,
T
k
x
,
T
d
T
=
f
x
,
t subject to the initial conditions(2)
u
x
,
0
=
g
x
x
∈
0,1
,
u
0
,
t
=
h
t
t
∈
0,1
,where D
α
x
u
x
,
t
g
x
,
t is fractional derivative of Caputo sense; when g
x
,
t
=
u
x
,
t, the initial problem is changed to nonlinear equation. f
x
,
t, g
x
,
t, u
x
,
t, and k
x
,
t are assumed to be casual functions of time and space on the section x
,
t
∈
0,1
×
0,1, where f
x
,
t, g
x
,
t, and k
x
,
t are known and u
x
,
t is the unknown, 0
<
α
x
≤
1.


## 2. Chebyshev Polynomials and Their Properties

The well-known Chebyshev polynomials are defined on the interval[
-
1,1
] and can be determined with the recurrence formula [16](3)
T
^
n
+
1
x
=
2
x
T
^
n
x
-
T
^
n
-
1
x
,
T
^
0
x
=
1
,
T
^
1
x
=
x
,
n
=
1,2
,
…
.The analytic form of the Chebyshev polynomials T
^
n
(
x
) of degree n is given by(4)
T
^
n
x
=
n
∑
i
=
0
n
/
2
-
1
i
2
n
-
2
i
-
1
n
-
i
-
1
!
i
!
n
-
2
i
!
x
n
-
2
i
,where n
/
2 denotes the integer part of n
/
2 and  n denotes positive integer. The orthogonality condition is(5)
∫
-
1
1
T
^
i
x
T
^
j
x
1
-
x
2
d
x
=
π
,
i
=
j
=
0
π
2
,
i
=
j
≠
0
0
,
i
≠
j
.In order to use these polynomials on the interval [
0,1
], we define the shifted Chebyshev polynomials by introducing the change of variable x
=
2
t
-
1. Therefore, the shifted Chebyshev polynomials are defined as T
n
∗
(
t
)
=
T
^
n
(
2
t
-
1
). The analytic form of the shifted Chebyshev polynomials T
n
∗
(
t
) of degree n is given by(6)
T
n
∗
t
=
n
∑
k
=
0
n
-
1
n
-
k
2
2
k
n
+
k
-
1
!
2
k
!
n
-
k
!
t
k
,
n
=
1,2
,
…
.Let (7)
Φ
t
=
T
0
∗
t
,
T
1
∗
t
,
…
,
T
n
∗
t
T
.The Chebyshev polynomials given by (6) can be expressed in the matrix form(8)
Φ
t
=
A
T
n
t
,where(9)
A
=
1
0
0
⋯
0
-
1
2
0
⋯
0
2
-
1
2
1
!
2
!
2
-
1
1
2
2
2
!
2
!
2
-
1
0
2
4
3
!
4
!
⋯
0
⋮
⋮
⋮
⋱
⋮
n
-
1
n
n
-
1
!
n
!
n
-
1
n
-
1
2
2
n
!
2
!
n
-
1
!
n
-
1
n
-
2
2
4
n
+
1
!
4
!
n
-
2
!
⋯
n
-
1
0
2
2
n
2
n
-
1
!
2
n
!
,
Τ
n
t
=
1
t
⋮
t
n
.Obviously(10)
Τ
n
t
=
A
-
1
Φ
t
.A function u
(
t
)
∈
L
2
(
0,1
) can be expressed in terms of the Chebyshev basis. In practice, only the first (
n
+
1
) term of Chebyshev polynomials is considered. Hence(11)
u
t
≅
∑
i
=
0
n
c
i
T
i
∗
t
=
c
T
Φ
t
,where c
=
[
c
0
,
c
1
,
…
,
c
n
]
T, c
i
(
i
=
0,1
,
2
,
…
,
n
) are called Chebyshev coefficients and c
=
Q
-
1
u
,
Φ
t. The dimension of Q is n
+
1
×
n
+
1; it is called the inner product matrix which is given by(12)
Q
=
∫
0
1
Φ
t
Φ
T
t
d
x
=
∫
0
1
A
T
n
t
A
T
n
t
T
d
t
=
A
∫
0
1
T
n
t
T
n
T
t
d
t
A
T
=
A
H
A
T
,where(13)
H
=
1
1
2
⋯
1
n
+
1
1
2
1
3
⋯
1
n
+
2
⋮
⋮
⋱
⋮
1
n
+
1
1
n
+
2
⋯
1
2
n
+
1
.For the function u
x
,
t
∈
L
2
0,1
×
0,1, we can also obtain its approximation by using Chebyshev polynomials(14)
u
x
,
t
≅
∑
i
=
0
n
∑
j
=
0
n
u
i
j
T
i
∗
t
T
j
∗
t
=
Φ
T
x
U
Φ
t
,where(15)
U
=
u
00
u
01
⋯
u
0
n
u
10
u
11
⋯
u
1
n
⋮
⋮
⋱
⋮
u
n
0
u
n
1
⋯
u
n
n
.Theorem 1 (see [16]).
The error in approximatingu
t by the sum of its first n terms is bounded by the sum of the absolute values of all the neglected coefficients. If(16)
u
n
t
=
∑
i
=
0
n
c
i
T
i
∗
tthen(17)
E
T
m
=
u
t
-
u
n
t
≤
∑
k
=
n
+
1
∞
c
kfor all u
t, all n, and all t
∈
[
-
1,1
].Theorem 2 (see [17]).
The Caputo fractional derivative of orderα for the shifted Chebyshev polynomials can be expressed in terms of the shifted Chebyshev polynomials themselves in the following form:(18)
D
α
T
i
∗
t
=
∑
k
=
α
i
∑
j
=
0
k
-
α
θ
i
,
j
,
k
T
j
∗
t
,where(19)
θ
i
,
j
,
k
=
-
1
i
-
k
2
i
i
+
k
-
1
!
Γ
k
-
α
+
1
/
2
Γ
k
+
1
/
2
i
-
k
!
Γ
k
-
α
-
j
+
1
Γ
k
-
α
+
j
+
1
,
j
=
0,1
,
…
.Theorem 3.
The errorE
T
∗
(
m
)
=
D
α
u
(
t
)
-
D
α
u
n
(
t
) in approximating D
α
u
(
t
) by D
α
u
n
(
t
) is bounded by E
T
∗
(
m
)
≤
∑
i
=
n
+
1
∞
c
i
∑
k
=
α
i
∑
j
=
0
k
-
α
θ
i
,
j
,
k.Proof.
A combination of (17) and (25) leads to (20)
E
T
∗
m
=
D
α
u
t
-
D
α
u
n
t
=
∑
i
=
n
+
1
∞
c
i
∑
k
=
α
i
∑
j
=
0
k
-
α
θ
i
,
j
,
k
T
j
∗
t
,but T
j
∗
≤
1, so we can get(21)
E
T
∗
m
≤
∑
i
=
n
+
1
∞
c
i
∑
k
=
α
i
∑
j
=
0
k
-
α
θ
i
,
j
,
k
,subtracting the truncated series from the infinite series, bounding each term in the difference, summing the bounds, and hence completing the proof of the theorem.


## 3. Operational Matrix of the Chebyshev Polynomials

### 3.1. Fractional Calculus

Before we introduce the Chebyshev polynomials operational matrix of the fractional integration, we first review some basic definitions of fractional calculus, which have been given in [18].Definition 4.
The Riemann-Liouville fractional integral operator of orderα
t:(22)
I
a
+
α
t
u
t
=
1
Γ
α
t
∫
a
+
t
t
-
T
α
t
-
1
u
T
d
T
,
t
>
0
Re
⁡
α
t
>
0
.Definition 5.
Riemann-Liouville fractional derivate with orderα
t:(23)
D
a
+
α
t
u
t
=
1
Γ
m
-
α
t
d
m
d
t
m
∫
a
t
u
τ
t
-
τ
α
t
-
m
+
1
d
τ
m
-
1
≤
α
t
<
m
.Definition 6.
Caputo’s fractional derivate with orderα
t,  0
<
α
t
≤
1: (24)
D
α
t
u
t
=
1
Γ
1
-
α
t
∫
0
+
t
t
-
τ
-
α
t
u
′
τ
d
τ
+
u
0
+
-
u
0
-
t
-
α
t
Γ
1
-
α
t
.If we assume the starting time in a perfect situation, we can obtain the definition as follows:(25)
D
α
t
u
t
=
1
Γ
1
-
α
t
∫
0
+
t
t
-
τ
-
α
t
u
′
τ
d
τ
0
<
α
t
<
1
.Generally, we adopt (25) as the definition of fractional derivate in Caputo sense. With the definition above, we can obtain the following formula 0
<
α
t
≤
1:(26)
D
α
t
c
=
0
D
α
t
x
β
=
0
β
=
0
Γ
β
+
1
Γ
β
+
1
-
α
t
x
β
-
a
t
β
=
1,2
,
3
,
…
.


### 3.2. The Operational Matrix of the Section as∂
u
x
,
t
/
∂
t in terms of Chebyshev Polynomials

The differentiation of vectorΦ
t in (7) can be given by(27)
Φ
′
t
=
D
Φ
t
,where D is the n
+
1
×
n
+
1 operational matrix of derivatives for Chebyshev polynomials. From (8) we have(28)
Φ
′
t
=
A
0
1
⋮
n
t
n
-
1
.Define the n
+
1
×
n matrix V
n
+
1
×
n and vector T
n
∗
t as(29)
V
n
+
1
×
n
=
0
0
⋯
0
1
0
⋯
0
0
2
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
n
,
T
n
∗
t
=
1
t
⋮
t
n
-
1
n
×
1
.Equation (28) may then be restated as(30)
Φ
′
t
=
A
V
n
+
1
×
n
T
n
∗
t
.Now we expand vector T
n
∗
t in terms of Φ
t. From (10), we get (31)
T
n
∗
t
=
B
∗
Φ
t
,where(32)
B
∗
=
A
1
-
1
,
A
2
-
1
,
…
,
A
n
-
1
T
.
A
k
-
1 is kth row of A
-
1,  k
=
1,2
,
…
,
n.Then we have(33)
Φ
′
t
=
A
V
n
+
1
×
n
B
∗
Φ
t
.Therefore we obtain the operational matrix of the section as ∂
u
x
,
t
/
∂
t as follows:(34)
∂
u
x
,
t
∂
t
=
Φ
T
x
U
Φ
′
t
=
Φ
T
x
U
A
V
n
+
1
×
n
B
∗
Φ
t
.


### 3.3. The Operational Matrix of the Section asD
α
t
u
x
,
t
g
x
,
t in terms of Chebyshev Polynomials

If we approximate the functionsu
x
,
t, g
x
,
t with Chebyshev polynomials, they can be written as u
x
,
t
=
Φ
T
x
U
Φ
t and g
x
,
t
=
Φ
T
x
G
Φ
t, where U is unknown and G is known. Then we have(35)
D
α
t
u
x
,
t
g
x
,
t
=
D
α
t
Φ
T
x
U
Φ
t
Φ
T
t
G
Φ
x
=
Φ
T
x
U
D
α
t
Φ
t
Φ
T
t
G
Φ
x
=
Φ
T
x
U
D
α
t
A
T
n
∗
t
A
T
n
∗
t
T
A
T
G
Φ
x
=
Φ
T
x
U
A
D
α
t
T
n
∗
t
T
n
∗
t
T
A
T
G
Φ
x
=
Φ
T
x
U
A
D
α
t
1
t
⋮
t
n
1
t
⋯
t
n
A
T
G
Φ
x
=
Φ
T
x
U
A
D
α
t
1
t
⋯
t
n
t
t
2
⋯
t
n
+
1
⋮
⋮
⋱
⋮
t
n
t
2
n
⋯
t
2
n
A
T
G
Φ
x
=
Φ
T
x
U
A
0
Γ
2
Γ
2
-
α
t
t
1
-
α
t
⋯
Γ
n
+
1
Γ
n
+
1
-
α
t
t
n
-
α
t
Γ
2
Γ
2
-
α
t
t
1
-
α
t
Γ
3
Γ
3
-
α
t
t
2
-
α
t
⋯
Γ
n
+
2
Γ
n
+
2
-
α
t
t
n
+
1
-
α
t
⋮
⋮
⋱
⋮
Γ
n
+
1
Γ
n
+
1
-
α
t
t
n
-
α
t
Γ
n
+
2
Γ
n
+
2
-
α
t
t
n
+
1
-
α
t
⋯
Γ
2
n
+
1
Γ
2
n
+
1
-
α
t
t
2
n
-
α
t
A
T
G
Φ
x
=
Φ
T
x
U
A
M
A
T
G
Φ
x
.Now we define(36)
M
=
0
Γ
2
Γ
2
-
α
t
t
1
-
α
t
⋯
Γ
n
+
1
Γ
n
+
1
-
α
t
t
n
-
α
t
Γ
2
Γ
2
-
α
t
t
1
-
α
t
Γ
3
Γ
3
-
α
t
t
2
-
α
t
⋯
Γ
n
+
2
Γ
n
+
2
-
α
t
t
n
+
1
-
α
t
⋮
⋮
⋱
⋮
Γ
n
+
1
Γ
n
+
1
-
α
t
t
n
-
α
t
Γ
n
+
2
Γ
n
+
2
-
α
t
t
n
+
1
-
α
t
⋯
Γ
2
n
+
1
Γ
2
n
+
1
-
α
t
t
2
n
-
α
t
.
M is called the operational matrix of the section as D
α
t
u
x
,
t
g
x
,
t with Chebyshev polynomials.So we have(37)
D
α
t
u
x
,
t
g
x
,
t
=
Φ
T
x
U
A
M
A
T
G
Φ
x
.


### 3.4. The Operational Matrix of the Section as∫
0
t
u
x
,
T
d
T in terms of Chebyshev Polynomials

The integration of the vectorΦ
t in (7) can be expressed as(38)
∫
0
t
Φ
t
d
T
=
P
Φ
t
,where P is the n
+
1
×
n
+
1 operational matrix of integration for Chebyshev polynomials. So we have(39)
∫
0
t
Φ
T
d
T
=
∫
0
t
A
T
n
T
d
T
=
A
∫
0
t
T
n
T
d
T
=
A
∫
0
t
1
T
⋮
T
n
d
T
=
A
t
1
2
t
2
⋮
1
n
+
1
t
n
+
1
=
A
p
T
p
,where A
p is an n
+
1
×
n
+
1 matrix:(40)
A
p
=
A
1
0
⋯
0
0
1
2
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
1
n
+
1
,
T
p
=
t
t
2
⋮
t
n
+
1
.Now we approximate the elements of vector T
p in terms of Φ
t. By (10), then we have(41)
t
k
=
A
k
+
1
-
1
Φ
t
,
k
=
1
,
…
,
n
,where A
k
+
1
-
1 is the k
+
1th row of A
-
1 for k
=
1
,
…
,
n. We just need to approximate t
n
+
1
=
c
n
+
1
T
Φ
t. By using c
=
Q
-
1
u
,
Φ
t, we have(42)
c
n
+
1
=
Q
-
1
∫
0
1
t
n
+
1
Φ
t
d
t
=
Q
-
1
∫
0
1
t
n
+
1
B
0
,
n
t
d
t
∫
0
1
t
n
+
1
B
1
,
n
t
d
t
⋮
∫
0
1
t
n
+
1
B
n
,
n
t
d
t
=
Q
-
1
2
n
+
2
n
0
2
n
+
1
n
+
1
n
1
2
n
+
1
n
+
2
⋯
n
n
2
n
+
1
2
n
+
1
T
. We define(43)
P
1
=
A
2
-
1
A
3
-
1
⋮
A
n
+
1
-
1
c
n
+
1
.Then we can get T
p
=
P
1
Φ
t. Therefore we have the operational matrix of integration as follows:(44)
P
=
A
p
P
1
.So we have(45)
∫
0
t
u
x
,
T
d
T
=
∫
0
t
Φ
T
x
U
Φ
T
d
T
=
Φ
T
x
U
∫
0
t
Φ
T
d
T
=
Φ
T
x
U
P
Φ
t
=
Φ
T
x
U
A
p
P
1
Φ
t
.


### 3.5. The Operational Matrix of the Section as∫
0
t
u
x
,
T
k
x
,
T
d
T in terms of Chebyshev Polynomials

Firstly, we approximate the functionk
x
,
t with Chebyshev polynomials; it can be written as k
x
,
t
=
Φ
T
x
K
Φ
t, and K is known. So we have(46)
∫
0
t
u
x
,
T
k
x
,
t
d
T
=
∫
0
t
Φ
T
x
U
Φ
T
Φ
T
T
K
Φ
x
d
T
=
Φ
T
x
U
∫
0
t
Φ
T
Φ
T
T
d
T
K
Φ
x
=
Φ
T
x
U
A
∫
0
t
1
T
⋯
T
n
T
T
2
⋯
T
n
+
1
⋮
⋮
⋱
⋮
T
n
T
n
+
1
⋯
T
2
n
d
T
A
T
K
Φ
x
=
Φ
T
x
U
A
t
1
2
t
2
⋯
1
n
+
1
t
n
+
1
1
2
t
2
1
3
t
3
⋯
1
n
+
2
t
n
+
2
⋮
⋮
⋱
⋮
1
n
+
1
t
n
+
1
1
n
+
2
t
n
+
2
⋯
1
2
n
+
1
t
2
n
+
1
A
T
K
Φ
x
=
Φ
T
x
U
A
R
A
T
K
Φ
x
.We define (47)
R
=
t
1
2
t
2
⋯
1
n
+
1
t
n
+
1
1
2
t
2
1
3
t
3
⋯
1
n
+
2
t
n
+
2
⋮
⋮
⋱
⋮
1
n
+
1
t
n
+
1
1
n
+
2
t
n
+
2
⋯
1
2
n
+
1
t
2
n
+
1
.
R is called the operational matrix of the section as ∫
0
t
u
x
,
T
k
x
,
T
d
T in terms of Chebyshev polynomials.Therefore the initial equation is transformed into the products of several dependent matrixes as follows:(48)
Φ
T
x
U
A
M
A
T
G
Φ
x
+
Φ
T
x
U
P
Φ
t
+
Φ
T
x
U
A
V
n
+
1
×
n
B
∗
Φ
t
+
Φ
T
x
U
A
R
A
T
K
Φ
x
=
f
x
,
t
.Dispersing (48) with x
i
,
t
j  i
=
1,2
,
…
,
n
;
j
=
1,2
,
…
,
n, by using a symbolic software such as “Mathematica,” we can get U. So the numerical solution of the original problem is obtained ultimately.


## 3.1. Fractional Calculus

Before we introduce the Chebyshev polynomials operational matrix of the fractional integration, we first review some basic definitions of fractional calculus, which have been given in [18].Definition 4.
The Riemann-Liouville fractional integral operator of orderα
t:(22)
I
a
+
α
t
u
t
=
1
Γ
α
t
∫
a
+
t
t
-
T
α
t
-
1
u
T
d
T
,
t
>
0
Re
⁡
α
t
>
0
.Definition 5.
Riemann-Liouville fractional derivate with orderα
t:(23)
D
a
+
α
t
u
t
=
1
Γ
m
-
α
t
d
m
d
t
m
∫
a
t
u
τ
t
-
τ
α
t
-
m
+
1
d
τ
m
-
1
≤
α
t
<
m
.Definition 6.
Caputo’s fractional derivate with orderα
t,  0
<
α
t
≤
1: (24)
D
α
t
u
t
=
1
Γ
1
-
α
t
∫
0
+
t
t
-
τ
-
α
t
u
′
τ
d
τ
+
u
0
+
-
u
0
-
t
-
α
t
Γ
1
-
α
t
.If we assume the starting time in a perfect situation, we can obtain the definition as follows:(25)
D
α
t
u
t
=
1
Γ
1
-
α
t
∫
0
+
t
t
-
τ
-
α
t
u
′
τ
d
τ
0
<
α
t
<
1
.Generally, we adopt (25) as the definition of fractional derivate in Caputo sense. With the definition above, we can obtain the following formula 0
<
α
t
≤
1:(26)
D
α
t
c
=
0
D
α
t
x
β
=
0
β
=
0
Γ
β
+
1
Γ
β
+
1
-
α
t
x
β
-
a
t
β
=
1,2
,
3
,
…
.


## 3.2. The Operational Matrix of the Section as∂
u
x
,
t
/
∂
t in terms of Chebyshev Polynomials

The differentiation of vectorΦ
t in (7) can be given by(27)
Φ
′
t
=
D
Φ
t
,where D is the n
+
1
×
n
+
1 operational matrix of derivatives for Chebyshev polynomials. From (8) we have(28)
Φ
′
t
=
A
0
1
⋮
n
t
n
-
1
.Define the n
+
1
×
n matrix V
n
+
1
×
n and vector T
n
∗
t as(29)
V
n
+
1
×
n
=
0
0
⋯
0
1
0
⋯
0
0
2
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
n
,
T
n
∗
t
=
1
t
⋮
t
n
-
1
n
×
1
.Equation (28) may then be restated as(30)
Φ
′
t
=
A
V
n
+
1
×
n
T
n
∗
t
.Now we expand vector T
n
∗
t in terms of Φ
t. From (10), we get (31)
T
n
∗
t
=
B
∗
Φ
t
,where(32)
B
∗
=
A
1
-
1
,
A
2
-
1
,
…
,
A
n
-
1
T
.
A
k
-
1 is kth row of A
-
1,  k
=
1,2
,
…
,
n.Then we have(33)
Φ
′
t
=
A
V
n
+
1
×
n
B
∗
Φ
t
.Therefore we obtain the operational matrix of the section as ∂
u
x
,
t
/
∂
t as follows:(34)
∂
u
x
,
t
∂
t
=
Φ
T
x
U
Φ
′
t
=
Φ
T
x
U
A
V
n
+
1
×
n
B
∗
Φ
t
.


## 3.3. The Operational Matrix of the Section asD
α
t
u
x
,
t
g
x
,
t in terms of Chebyshev Polynomials

If we approximate the functionsu
x
,
t, g
x
,
t with Chebyshev polynomials, they can be written as u
x
,
t
=
Φ
T
x
U
Φ
t and g
x
,
t
=
Φ
T
x
G
Φ
t, where U is unknown and G is known. Then we have(35)
D
α
t
u
x
,
t
g
x
,
t
=
D
α
t
Φ
T
x
U
Φ
t
Φ
T
t
G
Φ
x
=
Φ
T
x
U
D
α
t
Φ
t
Φ
T
t
G
Φ
x
=
Φ
T
x
U
D
α
t
A
T
n
∗
t
A
T
n
∗
t
T
A
T
G
Φ
x
=
Φ
T
x
U
A
D
α
t
T
n
∗
t
T
n
∗
t
T
A
T
G
Φ
x
=
Φ
T
x
U
A
D
α
t
1
t
⋮
t
n
1
t
⋯
t
n
A
T
G
Φ
x
=
Φ
T
x
U
A
D
α
t
1
t
⋯
t
n
t
t
2
⋯
t
n
+
1
⋮
⋮
⋱
⋮
t
n
t
2
n
⋯
t
2
n
A
T
G
Φ
x
=
Φ
T
x
U
A
0
Γ
2
Γ
2
-
α
t
t
1
-
α
t
⋯
Γ
n
+
1
Γ
n
+
1
-
α
t
t
n
-
α
t
Γ
2
Γ
2
-
α
t
t
1
-
α
t
Γ
3
Γ
3
-
α
t
t
2
-
α
t
⋯
Γ
n
+
2
Γ
n
+
2
-
α
t
t
n
+
1
-
α
t
⋮
⋮
⋱
⋮
Γ
n
+
1
Γ
n
+
1
-
α
t
t
n
-
α
t
Γ
n
+
2
Γ
n
+
2
-
α
t
t
n
+
1
-
α
t
⋯
Γ
2
n
+
1
Γ
2
n
+
1
-
α
t
t
2
n
-
α
t
A
T
G
Φ
x
=
Φ
T
x
U
A
M
A
T
G
Φ
x
.Now we define(36)
M
=
0
Γ
2
Γ
2
-
α
t
t
1
-
α
t
⋯
Γ
n
+
1
Γ
n
+
1
-
α
t
t
n
-
α
t
Γ
2
Γ
2
-
α
t
t
1
-
α
t
Γ
3
Γ
3
-
α
t
t
2
-
α
t
⋯
Γ
n
+
2
Γ
n
+
2
-
α
t
t
n
+
1
-
α
t
⋮
⋮
⋱
⋮
Γ
n
+
1
Γ
n
+
1
-
α
t
t
n
-
α
t
Γ
n
+
2
Γ
n
+
2
-
α
t
t
n
+
1
-
α
t
⋯
Γ
2
n
+
1
Γ
2
n
+
1
-
α
t
t
2
n
-
α
t
.
M is called the operational matrix of the section as D
α
t
u
x
,
t
g
x
,
t with Chebyshev polynomials.So we have(37)
D
α
t
u
x
,
t
g
x
,
t
=
Φ
T
x
U
A
M
A
T
G
Φ
x
.


## 3.4. The Operational Matrix of the Section as∫
0
t
u
x
,
T
d
T in terms of Chebyshev Polynomials

The integration of the vectorΦ
t in (7) can be expressed as(38)
∫
0
t
Φ
t
d
T
=
P
Φ
t
,where P is the n
+
1
×
n
+
1 operational matrix of integration for Chebyshev polynomials. So we have(39)
∫
0
t
Φ
T
d
T
=
∫
0
t
A
T
n
T
d
T
=
A
∫
0
t
T
n
T
d
T
=
A
∫
0
t
1
T
⋮
T
n
d
T
=
A
t
1
2
t
2
⋮
1
n
+
1
t
n
+
1
=
A
p
T
p
,where A
p is an n
+
1
×
n
+
1 matrix:(40)
A
p
=
A
1
0
⋯
0
0
1
2
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
1
n
+
1
,
T
p
=
t
t
2
⋮
t
n
+
1
.Now we approximate the elements of vector T
p in terms of Φ
t. By (10), then we have(41)
t
k
=
A
k
+
1
-
1
Φ
t
,
k
=
1
,
…
,
n
,where A
k
+
1
-
1 is the k
+
1th row of A
-
1 for k
=
1
,
…
,
n. We just need to approximate t
n
+
1
=
c
n
+
1
T
Φ
t. By using c
=
Q
-
1
u
,
Φ
t, we have(42)
c
n
+
1
=
Q
-
1
∫
0
1
t
n
+
1
Φ
t
d
t
=
Q
-
1
∫
0
1
t
n
+
1
B
0
,
n
t
d
t
∫
0
1
t
n
+
1
B
1
,
n
t
d
t
⋮
∫
0
1
t
n
+
1
B
n
,
n
t
d
t
=
Q
-
1
2
n
+
2
n
0
2
n
+
1
n
+
1
n
1
2
n
+
1
n
+
2
⋯
n
n
2
n
+
1
2
n
+
1
T
. We define(43)
P
1
=
A
2
-
1
A
3
-
1
⋮
A
n
+
1
-
1
c
n
+
1
.Then we can get T
p
=
P
1
Φ
t. Therefore we have the operational matrix of integration as follows:(44)
P
=
A
p
P
1
.So we have(45)
∫
0
t
u
x
,
T
d
T
=
∫
0
t
Φ
T
x
U
Φ
T
d
T
=
Φ
T
x
U
∫
0
t
Φ
T
d
T
=
Φ
T
x
U
P
Φ
t
=
Φ
T
x
U
A
p
P
1
Φ
t
.


## 3.5. The Operational Matrix of the Section as∫
0
t
u
x
,
T
k
x
,
T
d
T in terms of Chebyshev Polynomials

Firstly, we approximate the functionk
x
,
t with Chebyshev polynomials; it can be written as k
x
,
t
=
Φ
T
x
K
Φ
t, and K is known. So we have(46)
∫
0
t
u
x
,
T
k
x
,
t
d
T
=
∫
0
t
Φ
T
x
U
Φ
T
Φ
T
T
K
Φ
x
d
T
=
Φ
T
x
U
∫
0
t
Φ
T
Φ
T
T
d
T
K
Φ
x
=
Φ
T
x
U
A
∫
0
t
1
T
⋯
T
n
T
T
2
⋯
T
n
+
1
⋮
⋮
⋱
⋮
T
n
T
n
+
1
⋯
T
2
n
d
T
A
T
K
Φ
x
=
Φ
T
x
U
A
t
1
2
t
2
⋯
1
n
+
1
t
n
+
1
1
2
t
2
1
3
t
3
⋯
1
n
+
2
t
n
+
2
⋮
⋮
⋱
⋮
1
n
+
1
t
n
+
1
1
n
+
2
t
n
+
2
⋯
1
2
n
+
1
t
2
n
+
1
A
T
K
Φ
x
=
Φ
T
x
U
A
R
A
T
K
Φ
x
.We define (47)
R
=
t
1
2
t
2
⋯
1
n
+
1
t
n
+
1
1
2
t
2
1
3
t
3
⋯
1
n
+
2
t
n
+
2
⋮
⋮
⋱
⋮
1
n
+
1
t
n
+
1
1
n
+
2
t
n
+
2
⋯
1
2
n
+
1
t
2
n
+
1
.
R is called the operational matrix of the section as ∫
0
t
u
x
,
T
k
x
,
T
d
T in terms of Chebyshev polynomials.Therefore the initial equation is transformed into the products of several dependent matrixes as follows:(48)
Φ
T
x
U
A
M
A
T
G
Φ
x
+
Φ
T
x
U
P
Φ
t
+
Φ
T
x
U
A
V
n
+
1
×
n
B
∗
Φ
t
+
Φ
T
x
U
A
R
A
T
K
Φ
x
=
f
x
,
t
.Dispersing (48) with x
i
,
t
j  i
=
1,2
,
…
,
n
;
j
=
1,2
,
…
,
n, by using a symbolic software such as “Mathematica,” we can get U. So the numerical solution of the original problem is obtained ultimately.


## 4. Numerical Examples

To demonstrate the efficiency and the practicability of the proposed method based on Chebyshev polynomials method, we present some examples and find their solution via the method described in the previous section.Example 1.
Consider(49)
D
t
/
3
u
x
,
t
x
+
t
+
1
+
∂
u
x
,
t
∂
t
+
∫
0
t
u
x
,
T
d
T
+
∫
0
t
u
x
,
T
x
+
t
d
T
=
f
x
,
t
u
x
,
0
=
x
2
u
0
,
t
=
t
2
x
,
t
∈
0,1
×
0,1
,where(50)
f
x
,
t
=
2
t
+
t
3
3
+
t
4
4
+
t
3
x
3
+
t
x
2
+
t
2
x
2
2
+
t
x
3
-
3
t
1
-
t
/
3
6
t
9
+
8
t
-
6
-
9
+
t
t
x
-
9
+
t
-
6
+
t
-
3
+
t
Γ
1
-
t
/
3
.The exact solution of the above equation is u
x
,
t
=
x
2
+
t
2.
Takingn
=
2, dispersing x
i
=
k
i
/
3
-
1
/
6,  t
j
=
k
j
/
3
-
1
/
6
k
i
=
1,2
,
3
;
k
j
=
1,2
,
3, we can get the matrix U as follows:(51)
U
=
0
0
1
0
0.00217
0.99905
1
1.00036
1.99974
.The absolute error between the exact solution and the numerical solution is displayed in Figure 1.
Takingn
=
3, dispersing x
i
=
k
i
/
4
-
1
/
8,  t
j
=
k
j
/
4
-
1
/
8  k
i
=
1,2
,
3,4
;
k
j
=
1,2
,
3,4, we can get the matrix U as follows:(52)
U
=
0
0
1
3
1
0
0
1
3
1
1
3
1
3
2
3
4
3
1
1
4
3
2
.The absolute error between the exact solution and the numerical solution is displayed in Figure 2.Figure 1
The absolute error between the numerical solution and the exact solution whenn
=
2.Figure 2
The absolute error between the numerical solution and the exact solution whenn
=
3.Example 2.
Consider(53)
D
sin
⁡
t
/
3
u
x
,
t
x
t
+
∂
u
x
,
t
∂
t
+
∫
0
t
u
x
,
T
d
T
+
∫
0
t
u
x
,
T
x
+
t
d
T
=
f
x
,
t
u
x
,
0
=
1
+
x
2
u
0
,
t
=
1
+
t
2
x
,
t
∈
0,1
×
0,1
,where(54)
f
x
,
t
=
2
1
+
x
+
t
+
t
+
3
t
2
2
+
t
3
+
t
4
4
+
3
t
x
+
3
t
2
x
+
t
3
x
+
3
t
x
2
+
3
t
2
x
2
2
+
t
x
3
-
3
t
1
-
sin
⁡
t
/
3
x
54
1
+
t
+
x
2
+
1
+
x
sin
⁡
t
-
3
5
+
4
t
+
5
x
+
1
+
x
sin
⁡
t
Γ
1
-
sin
⁡
t
/
3
-
9
+
sin
⁡
t
-
6
+
sin
⁡
t
-
3
+
sin
⁡
t
.The exact solution of the above problem is u
x
,
t
=
1
+
x
+
t
2.
Takingn
=
2, dispersing x
i
=
k
i
/
3
-
1
/
6,  t
j
=
k
j
/
3
-
1
/
6
i
=
1,2
,
3
;
j
=
1,2
,
3, we can obtain the matrix U as follows:(55)
U
=
1
2
4
2
3.50385
5.99903
4
6.00062
8.99959
.The absolute error between the exact solution and the numerical solution is displayed in Table 1.
Takingn
=
3, dispersing x
i
=
k
i
/
3
-
1
/
6,  t
j
=
k
j
/
3
-
1
/
6  i
=
1,2
,
3
;
j
=
1,2
,
3, the matrix U is displayed as follows:(56)
U
=
1
5
3
8
3
4
5
3
23
9
34
9
16
3
8
3
34
9
47
9
7
4
16
3
7
9
.The absolute errorbetween the exact solution and the numerical solution is displayed in Table 2.
Takingn
=
4, dispersing x
i
=
k
i
/
5
-
1
/
10,  t
j
=
k
j
/
5
-
1
/
10  i
=
1,2
,
…
,
5
;
j
=
1,2
,
…
,
5, the matrix U is displayed as follows:(57)
U
=
1
3
2
13
6
3
4
3
2
2.125
2.91667
31
8
5
13
6
2.91667
3.83333
4.91667
6.1667
3
3.875
4.91667
6.125
7.5
4
5
37
6
7.5
9
.The absolute error between the exact solution and the numerical solution is displayed in Table 3.
Wheng
x
,
t
=
u
x
,
t, the initial equation becomes nonlinear equation. Example 3 describes the situation.Table 1
The absolute error between the numerical solution and the exact solution whenn
=
2.
t
=
0.1
t
=
0.3
t
=
0.5
t
=
0.7
t
=
0.9
x
=
0.0
0
0
0
0
0
x
=
0.1
2.2157
e
-
004
4.8610
e
-
004
5.1546
e
-
004
2.8785
e
-
004
1.5781
e
-
004
x
=
0.2
3.9341
e
-
004
8.7893
e
-
004
9.8541
e
-
004
5.3056
e
-
004
2.9527
e
-
004
x
=
0.3
5.5481
e
-
004
1.3762
e
-
003
1.7659
e
-
003
7.9111
e
-
004
3.5739
e
-
004
x
=
0.4
6.0429
e
-
004
1.6415
e
-
003
1.7751
e
-
003
8.4552
e
-
004
4.1018
e
-
004
x
=
0.5
6.5633
e
-
004
1.6016
e
-
003
1.6325
e
-
003
8.8767
e
-
004
4.3163
e
-
004
x
=
0.6
6.3392
e
-
004
1.4087
e
-
003
1.6103
e
-
003
9.1437
e
-
004
3.8174
e
-
004
x
=
0.7
5.6511
e
-
004
1.2886
e
-
003
1.5435
e
-
003
8.6002
e
-
004
2.9051
e
-
004
x
=
0.8
4.9215
e
-
004
1.2345
e
-
003
1.2001
e
-
003
7.5921
e
-
004
1.7393
e
-
004
x
=
0.9
3.8106
e
-
004
7.8650
e
-
004
8.754
e
-
004
5.5526
e
-
004
1.3122
e
-
005
x
=
1.0
1.8274
e
-
004
3.8541
e
-
004
4.3935
e
-
003
3.8035
e
-
004
1.9023
e
-
004Table 2
The absolute error between the numerical solution and the exact solution whenn
=
3.
t
=
0.1
t
=
0.3
t
=
0.5
t
=
0.7
t
=
0.9
x
=
0.0
0
0
0
0
0
x
=
0.1
1.87636
e
-
015
2.43249
e
-
015
4.37659
e
-
016
3.19872
e
-
015
4.65976
e
-
015
x
=
0.2
2.01345
e
-
015
2.58254
e
-
015
1.39627
e
-
015
5.14303
e
-
015
7.08487
e
-
015
x
=
0.3
2.67658
e
-
015
3.92180
e
-
015
4.47129
e
-
016
4.56189
e
-
015
5.18520
e
-
015
x
=
0.4
1.79236
e
-
015
2.62554
e
-
015
4.72489
e
-
016
5.64307
e
-
015
5.34316
e
-
015
x
=
0.5
1.67836
e
-
015
2.82354
e
-
015
4.93589
e
-
016
1.37827
e
-
015
3.34771
e
-
015
x
=
0.6
1.56727
e
-
015
3.10922
e
-
015
8.60978
e
-
016
1.93427
e
-
015
8.28178
e
-
015
x
=
0.7
2.17645
e
-
015
1.77326
e
-
015
1.39936
e
-
016
3.18642
e
-
015
2.26454
e
-
015
x
=
0.8
8.73278
e
-
016
2.82354
e
-
015
0
2.63545
e
-
015
4.24089
e
-
016
x
=
0.9
4.56389
e
-
016
1.93236
e
-
015
8.65478
e
-
016
4.62389
e
-
015
5.20703
e
-
015
x
=
1.0
1.80036
e
-
015
2.79454
e
-
015
8.44378
e
-
016
2.76054
e
-
015
5.70912
e
-
015Table 3
The absolute error between the numerical solution and the exact solution whenn
=
4.
t
=
0.1
t
=
0.3
t
=
0.5
t
=
0.7
t
=
0.9
x
=
0.0
0
0
0
0
0
x
=
0.1
4.23485
e
-
015
2.53645
e
-
016
3.67471
e
-
015
1.34255
e
-
014
1.84674
e
-
014
x
=
0.2
2.65749
e
-
015
4.43989
e
-
016
2.30435
e
-
015
1.25053
e
-
014
5.96979
e
-
015
x
=
0.3
2.34298
e
-
016
1.83436
e
-
015
3.34257
e
-
015
1.80056
e
-
015
1.34590
e
-
014
x
=
0.4
3.37542
e
-
015
1.65736
e
-
015
1.43257
e
-
015
1.56836
e
-
015
1.57838
e
-
014
x
=
0.5
2.39005
e
-
015
3.12362
e
-
015
1.56436
e
-
015
3.45782
e
-
015
1.10141
e
-
014
x
=
0.6
1.29456
e
-
015
2.64564
e
-
015
1.45637
e
-
015
4.87698
e
-
015
1.43013
e
-
014
x
=
0.7
1.69936
e
-
015
4.56239
e
-
016
3.56342
e
-
015
4.59329
e
-
015
4.56750
e
-
014
x
=
0.8
5.42507
e
-
015
2.76554
e
-
015
2.34525
e
-
015
1.58931
e
-
014
5.97889
e
-
014
x
=
0.9
7.86761
e
-
015
7.89961
e
-
015
3.27652
e
-
015
3.31246
e
-
014
4.14391
e
-
014
x
=
1.0
7.18553
e
-
015
1.54690
e
-
014
0
5.45820
e
-
014
5.56484
e
-
014Example 3.
Consider(58)
D
t
/
3
u
2
x
,
t
+
D
t
/
4
u
x
,
t
+
∂
2
u
x
,
t
∂
t
2
=
f
x
,
t
u
x
,
0
=
x
2
,
u
0
,
t
=
t
2
x
,
t
∈
0,1
×
0,1
,where(59)
f
x
,
t
=
2
+
36
t
2
-
t
/
3
54
t
2
+
-
12
+
t
-
9
+
t
x
2
-
12
+
t
-
9
+
t
-
6
+
t
-
3
+
t
Γ
1
-
t
/
3
+
32
t
2
-
t
/
4
32
-
12
t
+
t
2
Γ
1
-
t
/
4
.The exact solution of the above equation is u
x
,
t
=
x
2
+
t
2.
This is a nonlinear variable order fractional differential equation; the numerical solution can also be gained with the method proposed in Section3 when n
≤
2.
Takingn
=
2, dispersing x
i
=
k
i
/
2
-
1
/
4,  t
j
=
k
j
/
2
-
1
/
4  k
i
=
1,2
;
k
j
=
1,2, we can obtain the matrix U as follows:(60)
U
=
0
0
1
0
0
1
1
1
2
.The numerical solution obtained by our method and the exact solution are shown in Figures 3 and 4. The absolute error between the exact solution and the numerical solution is displayed in Figure 5.
Whenn
≥
3, the computation is very large and getting the numerical solution is a very difficult thing.Figure 3
The numerical solution for Example3 of n
=
2.Figure 4
The exact solution for Example3.Figure 5
The absolute error for Example3 of n
=
2.From Figures1–5, Tables 1–3, we can see that the absolute errors are very small and only a small number of Chebyshev polynomials are needed. Compared with the other methods proposed in [19, 20], the method in this paper has significant advantages. The calculating results also show that combined with Chebyshev polynomials the method in this paper can be effectively used in the numerical solution of the fractional equation. From the above results, the numerical solutions are in good agreement with the exact solution.


## 5. Conclusion

In the present paper, the application and scope of the Chebyshev polynomials have been extended to a class of variable order fractional integral-differential equation successfully. Actually we derive four kinds of operational matrixes using Chebyshev polynomials and use these to solve the variable order fractional integral-differential equation numerically. By solving the system of algebraic equations, numerical solutions are obtained. Numerical examples illustrate the powerfulness of the proposed method. The solutions obtained using the suggested method show that numerical solutions are in very good coincidence with the exact solution. The method can be applied by developing for the other fractional problem.



---

*Source: 902161-2015-10-01.xml*