# A Legendre tau-Spectral Method for Solving Time-Fractional Heat Equation with Nonlocal Conditions

Authors: A. H. Bhrawy; M. A. Alghamdi
Journal: The Scientific World Journal
 (2014)

---

## Abstract

We develop the tau-spectral method to solve the time-fractional heat equation (T-FHE) with nonlocal condition. In order to achieve highly accurate solution of this problem, the operational matrix of fractional integration (described in the Riemann-Liouville sense) for shifted Legendre polynomials is investigated in conjunction with tau-spectral scheme and the Legendre operational polynomials are used as the base function. The main advantage in using the presented scheme is that it converts the T-FHE with nonlocal condition to a system of algebraic equations that simplifies the problem. For demonstrating the validity and applicability of the developed spectral scheme, two numerical examples are presented. The logarithmic graphs of the maximum absolute errors is presented to achieve the exponential convergence of the proposed method. Comparing between our spectral method and other methods ensures that our method is more accurate than those solved similar problem.


---

## Body

## 1. Introduction

In recent years, many engineering and physical phenomena can be successfully described by models of fractional differential equations (FDEs); see, for instance, [1–7]. Thus many researchers have been interested in studying the properties of fractional calculus and finding stable and robust numerical and analytical schemes for solving FDEs such as spectral tau method [8–10], Crank-Nicolson method [11], compact finite difference approximation [12], Legendre wavelets method [13], Haar wavelet operational matrix method [14], iterative Laplace transform method [15], Lie symmetry analysis method [16], and other methods [17–20].Recently, spectral methods [21–23] have been applied to solve ordinary FDEs (see [24, 25]) while in [26, 27] the authors introduced the operational matrices of fractional derivatives with the help of the spectral methods to solve FDEs. This is not all; the partial FDEs are also investigated by using the spectral methods. In [28–31], the tau and collocation spectral methods are implemented in combination with the operational matrices of fractional integration for approximating the solution of some classes of space-fractional differential equations.The T-FHE is a generalization of the classical heat equation obtained by replacing the first order time derivative by a fractional derivative of orderν, 0<ν≤1. Ali and Jassim [32] used the homotopy perturbation method to solve the T-FHE, while in [33] the authors introduced a general iteration formula of variational iteration method for a solution of the T-FHE. Moreover, in [34] the differential transform method is applied to solve the T-FHE. In addition, Rostamy and Karimi [35] constructed the Bernstein operational matrix for the fractional derivatives and used it together with spectral method to solve the T-FHE.In this paper, we consider the T-FHE with the nonlocal condition [36]:
(1)∂νu(x,t)∂tν-∂2u(x,t)∂x2=q(x,t),0<x≤1,0<t≤1,
subject to
(2)u(x,0)=u(x,1)+f(x),0<x≤1,u(0,t)=g0(t),u(1,t)=g1(t),0<t≤1,
where 0<ν≤1, u(x,t) is the temperature as a function of space x and time t, and q(x,t) is known source term. Our main aim is to achieve highly accurate solution of the T-FHE with nonlocal conditions (1) and (2). The tau-spectral method is applied based on the shifted Legendre polynomials as a basis function with the help of the operational matrix of fractional integration of such polynomials. Two numerical examples are introduced and solved using the presented technique to show its accuracy and validity. Also, we introduce comparisons between our numerical results and those obtained using the implicit difference approximation (IDA).This paper is arranged in the following way: in Section2 we introduce some definitions and notations of fractional calculus with some properties of Legendre polynomials. In Section 3 we apply our algorithm for the solution of the T-FHE with nonlocal condition. In Section 4 two numerical examples and comparisons between our results and those obtained by the IDA are introduced. Also in Section 5, a conclusion is presented.


## 2. Preliminaries and Notations

### 2.1. Fractional Calculus Definitions

Riemann-Liouville and Caputo fractional definitions are the two most used from other definitions of fractional derivatives which have been introduced recently.Definition 1.
The integral of orderγ≥0 (fractional) according to Riemann-Liouville is given by
(3)Iνf(x)=1Γ(ν)∫0x‍(x-t)ν-1f(t)dt,ν>0,x>0,I0f(x)=f(x),
where
(4)Γ(ν)=∫0∞‍xν-1e-xdx
is gamma function.The operatorIν satisfies the following properties:
(5)IνIμf(x)=Iν+μf(x),IνIμf(x)=IμIνf(x),Iνxβ=Γ(β+1)Γ(β+1+ν)xβ+ν.Definition 2.
The Caputo fractional derivative of orderν is defined by
(6)Dνf(x)=1Γ(m-ν)∫0x‍(x-t)m-ν-1dmdtmf(t)dt,fffffffffm-1<ν≤m,x>0,
where m is the ceiling function of ν.The operatorDγ satisfies the following properties:
(7)DνC=0,(Cisconstant)IνDνf(x)=f(x)-∑i=0m-1‍f(i)(0+)xii!,Dνxβ=Γ(β+1)Γ(β+1-ν)xβ-ν,Dν(λf(x)+μg(x))=λDνf(x)+μDνg(x).


### 2.2.  Shifted Legendre Polynomials

Assuming that the Legendre polynomial of degreej is denoted by Pj(z) (defined on the interval (-1,1)), then Pj(z) may be generated by the recurrence formulae
(8)Pj+1(z)=2j+1j+1zPj(z)-jj+1Pj-1(z),j=1,2,…,P0(z)=1,P1(z)=z.Consideringz=2x-1, Legendre polynomials are defined on the interval (0,1) that may be called shifted Legendre polynomials Pj*(x) that were generated using the following recurrence formulae:
(9)Pj+1*(x)=2j+1j+1(2x-1)Pj*(x)-jj+1Pj-1*(x),gggggggggggggggrgggggj=1,2,…,P0*(x)=1,P1*(x)=2x-1.The orthogonality relation is(10)∫01‍Pi*(x)Pj*(x)dx={12j+1,fori=j,0,fori≠j.The explicit analytical form of shifted Legendre polynomialPj*(x) of degree j may be written as
(11)Pj*(x)=∑k=0j‍(-1)j+k(j+k)!xk(j-k)!(k!)2,
and this in turn enables one to get
(12)Pj*(0)=(-1)j,Pj*(1)=1.Any square integrable functionu(x) defined on the interval (0,1) may be expressed in terms of shifted Legendre polynomials as
(13)u(x)=∑j=0∞‍ajPj*(x),
from which the coefficients aj are given by
(14)aj=(2j+1)∫01‍u(x)Pj*(x)dx,j=0,1,….If we approximateu(x) by the first (N+1)-terms, then we can write
(15)uN(x)=∑j=0N‍ajPj*(x),
which alternatively may be written in the matrix form
(16)uN(x)≃ATΨN(x),
with
(17)AT≡[a0,a1,…,aN],ΨN(x)≡[P0*(x),P1*(x),…,PN*(x)]T.
Similarly, let u(x,t) be an infinitely differentiable function defined on 0<x≤1 and 0<t≤1. Then it may be expressed as
(18)uM,N(x,t)≃∑i=0M∑j=0N‍uijPi*(t)Pj*(x)=ΨMT(t)UΨN(x),
with
(19)U=(u00u01⋯u0Nu10u11⋯u1N⋮⋮⋯⋮uM0uM1⋯uMN),(20)uij=(2i+1)(2j+1)∬01u(x,t)Pi*(t)Pj*(x)dxdt,ggggggggggggi=0,1,…,M,j=0,1,…,N.Theorem 3.
The first derivative of the shifted Legendre vectorΨN(x) may be expressed as
(21)dΨN(x)dx=DΨN(x),
where D is the (N+1)×(N+1) operational matrix of derivative given by
(22)D=(dij)={2(2j+1),forj=i-k,k=1,3,…,N,ifNisodd,k=1,3,…,N-1,ifNiseven,0,otherwise.Repeated use of (21) enables one to write
(23)dqΨN(x)dxq=DqΨN(x),
where q is a natural number and Dq means matrix power.Theorem 4.
The Riemann-Liouville fractional integral of orderν of the shifted Legendre polynomial vector ΨM(t) is given by
(24)IνΨM(t)=PνΨM(t),
where Pν is the (M+1)×(M+1) operational matrix of fractional integral of order ν and is defined by
(25)Pν=(∑k=00‍ξ(0,0,k)∑k=00‍ξ(0,1,k)⋯∑k=00‍ξ(0,M,k)∑k=01‍ξ(1,0,k)∑k=01‍ξ(1,1,k)⋯∑k=01‍ξ(1,M,k)⋮⋮⋯⋮∑k=0i‍ξ(i,0,k)∑k=0i‍ξ(i,1,k)⋯∑k=0i‍ξ(i,M,k)⋮⋮⋯⋮∑k=0M‍ξ(M,0,k)∑k=0M‍ξ(M,1,k)⋯∑k=0M‍ξ(M,M,k)),
where
(26)ξ(i,j,k)=(2j+1)∑l=0j‍(×(l!)2(k+l+ν+1))-1(-1)i+j+k+l(i+k)!(l+j)!×((l!)2(i-k)!k!Γ(k+ν+1)(j-l)!×(l!)2(k+l+ν+1))-1).(See [37] for proof.)


## 2.1. Fractional Calculus Definitions

Riemann-Liouville and Caputo fractional definitions are the two most used from other definitions of fractional derivatives which have been introduced recently.Definition 1.
The integral of orderγ≥0 (fractional) according to Riemann-Liouville is given by
(3)Iνf(x)=1Γ(ν)∫0x‍(x-t)ν-1f(t)dt,ν>0,x>0,I0f(x)=f(x),
where
(4)Γ(ν)=∫0∞‍xν-1e-xdx
is gamma function.The operatorIν satisfies the following properties:
(5)IνIμf(x)=Iν+μf(x),IνIμf(x)=IμIνf(x),Iνxβ=Γ(β+1)Γ(β+1+ν)xβ+ν.Definition 2.
The Caputo fractional derivative of orderν is defined by
(6)Dνf(x)=1Γ(m-ν)∫0x‍(x-t)m-ν-1dmdtmf(t)dt,fffffffffm-1<ν≤m,x>0,
where m is the ceiling function of ν.The operatorDγ satisfies the following properties:
(7)DνC=0,(Cisconstant)IνDνf(x)=f(x)-∑i=0m-1‍f(i)(0+)xii!,Dνxβ=Γ(β+1)Γ(β+1-ν)xβ-ν,Dν(λf(x)+μg(x))=λDνf(x)+μDνg(x).


## 2.2.  Shifted Legendre Polynomials

Assuming that the Legendre polynomial of degreej is denoted by Pj(z) (defined on the interval (-1,1)), then Pj(z) may be generated by the recurrence formulae
(8)Pj+1(z)=2j+1j+1zPj(z)-jj+1Pj-1(z),j=1,2,…,P0(z)=1,P1(z)=z.Consideringz=2x-1, Legendre polynomials are defined on the interval (0,1) that may be called shifted Legendre polynomials Pj*(x) that were generated using the following recurrence formulae:
(9)Pj+1*(x)=2j+1j+1(2x-1)Pj*(x)-jj+1Pj-1*(x),gggggggggggggggrgggggj=1,2,…,P0*(x)=1,P1*(x)=2x-1.The orthogonality relation is(10)∫01‍Pi*(x)Pj*(x)dx={12j+1,fori=j,0,fori≠j.The explicit analytical form of shifted Legendre polynomialPj*(x) of degree j may be written as
(11)Pj*(x)=∑k=0j‍(-1)j+k(j+k)!xk(j-k)!(k!)2,
and this in turn enables one to get
(12)Pj*(0)=(-1)j,Pj*(1)=1.Any square integrable functionu(x) defined on the interval (0,1) may be expressed in terms of shifted Legendre polynomials as
(13)u(x)=∑j=0∞‍ajPj*(x),
from which the coefficients aj are given by
(14)aj=(2j+1)∫01‍u(x)Pj*(x)dx,j=0,1,….If we approximateu(x) by the first (N+1)-terms, then we can write
(15)uN(x)=∑j=0N‍ajPj*(x),
which alternatively may be written in the matrix form
(16)uN(x)≃ATΨN(x),
with
(17)AT≡[a0,a1,…,aN],ΨN(x)≡[P0*(x),P1*(x),…,PN*(x)]T.
Similarly, let u(x,t) be an infinitely differentiable function defined on 0<x≤1 and 0<t≤1. Then it may be expressed as
(18)uM,N(x,t)≃∑i=0M∑j=0N‍uijPi*(t)Pj*(x)=ΨMT(t)UΨN(x),
with
(19)U=(u00u01⋯u0Nu10u11⋯u1N⋮⋮⋯⋮uM0uM1⋯uMN),(20)uij=(2i+1)(2j+1)∬01u(x,t)Pi*(t)Pj*(x)dxdt,ggggggggggggi=0,1,…,M,j=0,1,…,N.Theorem 3.
The first derivative of the shifted Legendre vectorΨN(x) may be expressed as
(21)dΨN(x)dx=DΨN(x),
where D is the (N+1)×(N+1) operational matrix of derivative given by
(22)D=(dij)={2(2j+1),forj=i-k,k=1,3,…,N,ifNisodd,k=1,3,…,N-1,ifNiseven,0,otherwise.Repeated use of (21) enables one to write
(23)dqΨN(x)dxq=DqΨN(x),
where q is a natural number and Dq means matrix power.Theorem 4.
The Riemann-Liouville fractional integral of orderν of the shifted Legendre polynomial vector ΨM(t) is given by
(24)IνΨM(t)=PνΨM(t),
where Pν is the (M+1)×(M+1) operational matrix of fractional integral of order ν and is defined by
(25)Pν=(∑k=00‍ξ(0,0,k)∑k=00‍ξ(0,1,k)⋯∑k=00‍ξ(0,M,k)∑k=01‍ξ(1,0,k)∑k=01‍ξ(1,1,k)⋯∑k=01‍ξ(1,M,k)⋮⋮⋯⋮∑k=0i‍ξ(i,0,k)∑k=0i‍ξ(i,1,k)⋯∑k=0i‍ξ(i,M,k)⋮⋮⋯⋮∑k=0M‍ξ(M,0,k)∑k=0M‍ξ(M,1,k)⋯∑k=0M‍ξ(M,M,k)),
where
(26)ξ(i,j,k)=(2j+1)∑l=0j‍(×(l!)2(k+l+ν+1))-1(-1)i+j+k+l(i+k)!(l+j)!×((l!)2(i-k)!k!Γ(k+ν+1)(j-l)!×(l!)2(k+l+ν+1))-1).(See [37] for proof.)


## 3. Legendre tau-Spectral Method

In this section, the Legendre operational matrix of fractional integrals is applied with the help of Legendre tau-spectral method to solve the T-FHE with the nonlocal condition.Consider the T-FHE with the nonlocal condition(27)∂νu(x,t)∂tν-∂2u(x,t)∂x2=q(x,t),0<x≤1,0<t≤1,0<ν≤1,u(x,0)=u(x,1)+f(x),0<x≤1,u(0,t)=g0(t),u(1,t)=g1(t),0<t≤1.
We integrate (27) of order ν and making use of (7), we have
(28)u(x,t)-(u(x,1)+f(x))-Itν(∂2u(x,t)∂x2)=Itνq(x,t),0<ν≤1,u(0,t)=g0(t),u(1,t)=g1(t),0<t≤1.In order to use tau-spectral method based on the shifted Legendre operational matrix for fractional integrals to solve the fully integrated problem (28), we approximate (x,t),f(x), and q(x,t) by the shifted Legendre polynomials as
(29)uM,N(x,t)=ΨMT(t)UΨN(x),fN(x)=ΨMT(t)FΨN(x),qM,N(x,t)=ΨMT(t)QΨN(x),
where U is the unknown coefficients (M+1)×(N+1) matrix and F and Q are known matrices that can be written as
(30)F=(f0f1⋯fN00⋯0⋮⋮⋯⋮00⋯0),Q=(q00q01⋯q0Nq10q11⋯q1N⋮⋮⋯⋮qM0qM1⋯qMN),
where fj and qij are given as in (14) and (20), respectively.Using (29), it is easy to write
(31)uM,N(x,1)=ΨMT(1)UΨN(x),=(∑k=0M‍uk0∑k=0M‍uk1⋯∑k=0M‍ukN)ΨN(x),=ΨMT(t)VΨN(x),
where V is a (M+1)×(N+1) matrix that can be written as
(32)V=(∑k=0M‍uk0∑k=0M‍uk1⋯∑k=0M‍ukN00⋯0⋮⋮⋯⋮00⋯0).
Making use of (23), (24), and (29) enables one to write
(33)Itν(∂2u(x,t)∂x2)≃(ItνΨMT(t))U(D2ΨN(x)),=ΨMT(t)PνTUD2ΨN(x).
In addition, if we use (24) and (29), we obtain
(34)ItνqM,N(x,t)=ΨMT(t)PνTQΨN(x).
Equations (31) and (34) enable one to write the residual RM,N(x,t) for (28) in the form
(35)RM,N(x,t)=ΨMT(t)[U-V-F-PνTUD2-PνTQ]ΨN(x).
As in a typical tau method (see [22, 38, 39]) we generate (M+1)(N-1) linear algebraic equations in the unknown expansion coefficients, uij, i=0,1,…, M;j=0,1,…,N-2, namely;
(36)∬01RM,N(x,t)Pi*(t)Pj*(x)dxdt=0,ggi=0,1,…,M,j=0,1,…,N-2,
and the rest of linear algebraic equations are obtained from the boundary conditions, as
(37)ΨMT(t)UΨN(0)=g0(ti),ΨMT(t)UΨN(1)=g1(ti),fi=0,1,…,M,
where ti,i=0,1,…,M are the roots of PM+1*(t). The number of the unknown coefficients uij is equal to (M+1)(N+1) and can be obtained from (36) and (37). Consequently uM,N(x,t) given in (29) can be calculated.


## 4. Numerical Experiments

In order to highlight the accuracy of the presented scheme, we implement it to solve two numerical examples, and also comparisons between their exact solutions with the approximate solutions achieved using the presented scheme and with those achieved using other methods are made.Example 1.
We consider the following problem [36]:
(38)∂0.5u(x,t)∂t0.5-∂2u(x,t)∂x2=2t1.5sin(2πx)Γ(1.5)+4π2t2sin(2πx),u(x,0)=u(x,1)-sin(2πx),0<x≤1,u(0,t)=0,u(1,t)=0,0<t≤1,
with exact solution u(x,t)=t2sin(2πx).Karatay et al. [36] introduced this problem and applied the IDA method to approximate its solution at various choices of time and space nodes M and N.We apply our numerical scheme for this problem. In order to show that our scheme is more accurate than the IDA method, in Table1, we compare the maximum absolute errors (MAEs) achieved using our scheme with those obtained using the IDA [36] method at different values of N, (N=M). Moreover, Figure 1 plots the absolute error function at M=N=18, while Figure 2 plots the absolute error function for t=0.5 at M=N=18.Table 1
Comparison of our scheme with the IDA [36] at various choices of N, (N=M) for Example 1.
Our scheme
IDA [36]
N
=
M
MAEs
N
=
M
MAEs
6
2.541507
·
1
0
-
2
4
2.297695
·
1
0
-
1
8
8.774320
·
1
0
-
4
8
5.383793
·
1
0
-
2
10
2.118341
·
1
0
-
5
16
1.391800
·
1
0
-
2
12
1.832678
·
1
0
-
6
32
3.843610
·
1
0
-
3
14
1.037991
·
1
0
-
6
64
1.152111
·
1
0
-
3
16
6.590323
·
1
0
-
7
128
3.844224
·
1
0
-
4
18
4.383403
·
1
0
-
7
256
1.447756
·
1
0
-
4Figure 1
Absolute error function atN=M=18 for Example 1.Figure 2
Absolute error function atN=M=18 with t=0.5 for Example 1.Example 2.
Consider the following problem:(39)∂νu(x,t)∂tν-∂2u(x,t)∂x2=2t2-νln⁡(1+x(1-x))Γ(3-ν)+t2(2x2-2x-1)(x2-x+1)2,0<ν≤1,u(x,0)=u(x,1)-ln⁡(1+x(1-x)),0<x≤1,u(0,t)=0,u(1,t)=0,0<t≤1,
with exact solution u(x,t)=t2ln⁡(1+x(1-x)).Karatay et al. [36] introduced this problem and solved it for two choices of ν, ν=0.45,0.95 at different values of N and M. Table 2 lists the MAEs for ν=0.45,0.95 using our scheme at N=M=4,6,8,10,12,14,16,18 and a comparison with those obtained in [36] at M=16,N=2,4,8,16,32,64,128,256 and N=16, M=2,4,8,16,32,64,128,256. Figures 3 and 4 plot the absolute error functions at M=N=20 with ν1=0.45 and ν=0.95, respectively. Finally, in order to demonstrate the convergence of the proposed method, in Figure 5, we plot the logarithmic graphs of the maximum absolute errors (log10Error) at two choices of ν,ν=0.45,0.95 and various choices of N, (N=M), by using the presented algorithm.Table 2
Comparison of our scheme with the IDA [36] at various choices of N and M for Example 2.
Our scheme
IDA [36]
N
ν
=
0.45
ν
=
0.95
M
=
16
N
=
16
N
ν
=
0.45
ν
=
0.95
M
ν
=
0.45
ν
=
0.95
4
1.08
·
1
0
-
3
1.10
·
1
0
-
3
2
1.46
·
1
0
-
2
7.45
·
1
0
-
3
2
3.88
·
1
0
-
2
4.30
·
1
0
-
2
6
5.84
·
1
0
-
5
5.86
·
1
0
-
5
4
7.59
·
1
0
-
3
4.06
·
1
0
-
3
4
9.51
·
1
0
-
3
9.47
·
1
0
-
3
8
6.33
·
1
0
-
6
5.01
·
1
0
-
6
8
4.03
·
1
0
-
3
2.30
·
1
0
-
3
8
3.65
·
1
0
-
3
2.96
·
1
0
-
3
10
1.86
·
1
0
-
6
1.15
·
1
0
-
6
16
2.24
·
1
0
-
3
1.41
·
1
0
-
3
16
2.24
·
1
0
-
3
1.41
·
1
0
-
3
12
9.03
·
1
0
-
7
5.12
·
1
0
-
7
32
1.35
·
1
0
-
3
9.59
·
1
0
-
4
32
1.89
·
1
0
-
3
1.02
·
1
0
-
3
14
5.07
·
1
0
-
7
2.81
·
1
0
-
7
64
9.02
·
1
0
-
4
7.33
·
1
0
-
4
64
1.81
·
1
0
-
3
9.31
·
1
0
-
4
16
3.06
·
1
0
-
7
1.69
·
1
0
-
7
128
6.78
·
1
0
-
4
6.19
·
1
0
-
4
128
1.78
·
1
0
-
3
9.07
·
1
0
-
4
18
1.95
·
1
0
-
7
1.69
·
1
0
-
7
256
5.66
·
1
0
-
4
5.62
·
1
0
-
4
256
1.78
·
1
0
-
3
9.01
·
1
0
-
4Figure 3
Absolute error function atN=M=20 with γ=0.45 for Example 2.Figure 4
Absolute error function atN=M=20 with γ=0.95 for Example 2.Figure 5
Absolute error function atN=M=20 with γ=0.95 for Example 2.From Tables1 and 2 and Figures 1 and 2 introduced above, it is shown that the proposed scheme is more accurate than the IDA method introduced by Karatay et al. [36].


## 5. Conclusion

An effective and accurate numerical scheme was developed to approximate the solution of the T-FHE with the nonlocal condition. The developed approach is based on the Legendre tau-spectral method combined with the operational matrix of fractional integration (described in the Riemann-Liouville sense) for orthogonal polynomials. A good approximation of the exact solution was achieved by using a limited number of the basis function.The logarithmic graphs of the maximum absolute errors were presented to achieve the exponential convergence of the proposed method. Comparisons between our approximate solutions of test problems with their exact solutions and the approximate solutions achieved by the IDA method were introduced to confirm the validity and accuracy of our scheme.



---

*Source: 706296-2014-06-25.xml*