# Sinc-Chebyshev Collocation Method for a Class of Fractional Diffusion-Wave Equations

Authors: Zhi Mao; Aiguo Xiao; Zuguo Yu; Long Shi
Journal: The Scientific World Journal
 (2014)

---

## Abstract

This paper is devoted to investigating the numerical solution for a class of fractional diffusion-wave equations
with a variable coefficient where the fractional derivatives are described in the Caputo sense.
The approach is based on the collocation technique where the shifted Chebyshev polynomials in time
and the sinc functions in space are utilized, respectively. The problem is reduced to the solution of a system of linear algebraic equations.
Through the numerical example, the procedure is tested and the efficiency of the proposed method is confirmed.


---

## Body

## 1. Introduction

Fractional models have been increasingly shown by many scientists to describe adequately the problems with memory and nonlocal properties in fluid mechanics, viscoelasticity, physics, biology, chemistry, finance, and other areas of applications [1–6]. In particular, the fractional diffusion-wave equation has been used to model many important physical phenomena ranging from amorphous, colloid, glassy, and porous materials through fractals, percolation clusters, and random and disordered media to comb structures, dielectrics and semiconductors, polymers, and biological systems [7–10]. It is a generalization of the classical diffusion-wave equation by replacing the integer-order time derivative with a fractional derivative of order α (1<α<2). This equation can be derived from the anomalous superdiffusion in continuous time random walk which is generally non-Markovian processes [11].Although the considerable work on the numerical solution of fractional diffusion equations has been done [12–15], there are very limited numerical methods for solving the fractional diffusion-wave equations [16–18]. However, all the above mentioned papers dealt with the fractional diffusion-wave equations by finite difference methods. It is well known that any algorithm based on the finite difference discretization of a fractional derivative has to take into account its memory or nonlocal structure; thus this means a high storage requirement [19].In the present paper, we consider the following differential equation with the Caputo fractional derivative and a variable coefficient:(1)∂αu(x,t)∂tα=a(x,t)∂2u(x,t)∂x2+f(x,t),kka<x<b,0<t≤τ,
with the initial conditions,
(2)u(x,0)=ϕ(x),∂u(x,0)∂t=ψ(x),a<x<b,
and the boundary conditions,
(3)u(a,t)=0,u(b,t)=0,0<t≤τ,
where x∈[a,b] and t∈(0,τ] are space and time variables, respectively, a(x,t) is a continuous function, and f(x,t) denotes the field variable. For 1<α<2, the fractional equation (1) is known as the fractional diffusion-wave equation which fills the gaps between the diffusion equation and wave equation [16, 20].We develop a sinc-Chebyshev collocation method to solve numerically problem (1) with (2) and (3). Since a fractional derivative is a nonlocal operator, it is natural to consider a global scheme such as the collocation method for its numerical solution [19, 21]. The required approximate solution is expanded as a series with the elements of shifted Chebyshev polynomials in time and sinc functions in space with unknown coefficients. By utilizing the collocation technique and some properties of the shifted Chebyshev polynomials and sinc functions, the problem is reduced to the solution to a system of linear algebraic equations. And a matrix representation of the system is obtained to calculate the solution. The presented method is effective and convenient.The remainder of this paper is organized as follows: in the next section, we introduce some necessary definitions and relevant results for developing this method. Section3 is devoted to constructing and analyzing the numerical algorithm. As a result, a system of linear algebraic equations is formed and the solution of the considered problem is obtained. In Section 4, the numerical example is given to demonstrate the effectiveness and convergence of the proposed method. A brief conclusion is given in the final section.


## 2. Notations and Some Preliminary Results

In this section, we introduce some basic definitions and derive several preliminary results for developing the presented method.


### 2.1. The Caputo Fractional Derivative

Definition 1 (see [22]).
Letα∈ℝ+. The operator Jaα defined on L1[a,b] by
(4)Jaαf(t)=1Γ(α)∫at(t-s)α-1f(s)ds‍
for a≤t≤b is called the Riemann-Liouville fractional integral operator of order α. For, α=0, we set Ja0:=I, that is, the identity operator.Definition 2 (see [22]).
Letα∈ℝ+ and n=⌈α⌉. The Caputo fractional differential operator DaαC for a≤t≤b is defined as
(5)DaαCf(t)=Jan-αDnf(t)=1Γ(n-α)∫at(t-s)n-α-1f(n)(s)ds‍.


### 2.2. The Composite Translated Sinc Functions

The sinc functions and their properties are discussed in [23, 24]. For any h>0, the translated sinc functions with equidistant space nodes are given as
(6)S(k,h)(z)=sinc(z-khh),k=0,±1,±2,…,
where the sinc functions are defined on the whole real line by
(7)sinc(x)={sin(πx)πx,x≠0,1,x=0.
If f is defined on ℝ, then for any h>0 the series
(8)C(f,h)(z)=∑k=-∞∞f(kh)S(k,h)(z)‍
is called the Whittaker cardinal expansion of f whenever this series converges. f can be approximated by truncating (8).To construct our needed approximations on the interval[a,b], we choose
(9)φ(x)=ln⁡(x-ab-x)
which maps the finite interval [a,b] onto ℝ. The basic functions on [a,b] are taken to be the composite translated sinc functions:
(10)Sφ(k,h)(x)=S(k,h)(φ(x))=sinc(φ(x)-khh).
Thus we may define the inverse image of the equidistant space node {ih} as
(11)xi=φ-1(ih)=a+beih1+eih,i=0,±1,±2,….
The class of functions such that the known exponential convergence rate exists for the sinc interpolation is denoted by B(DE) and defined in the following text.Definition 3 (see [21]).
LetB(DE) be the class of functions f which are analytic in DE and satisfy
(12)∫φ-1(x+L)|f(z)dz|‍⟶0,x⟶±∞,
where L={iυ:|υ|<d≤π/2}, and
(13)∫∂DE|f(z)dz|‍<∞
on the boundary of DE (denoted ∂DE).Theorem 4 (see [21, 23]).
Ifφ′f∈B(DE), then, for all x∈[a,b],
(14)|f(x)-∑k=-∞∞f(xk)Sφ(k,h)(x)‍|≤2N(φ′f)πde-πd/h.
Further, one assumes that there are positive constants C and β so that |f(x)|≤Cexp⁡(-β|φ(x)|). And if one selects h=πd/βN≤2πd/ln⁡2, then,
(15)|dmf(x)dxm-∑k=-NNf(xk)dmdxmSφ(k,h)(x)‍|≤KN(m+1)/2exp⁡(-πdβN)
for all m=0,1,…,n.The above expressions show that the sinc interpolation onB(DE) converges exponentially. We also require the following derivatives of the composite translated sinc functions evaluated at the nodes. Consider
(16)δk,i(0)=[Sφ(k,h)(x)]|x=xi={1,k=i,0,k≠i,(17)δk,i(1)=ddφ[Sφ(k,h)(x)]|x=xi={0,k=i,(-1)i-k(i-k)h,k≠i,δk,i(2)=d2dφ2[Sφ(k,h)(x)]|x=xi={-π23h2,k=i,-2(-1)i-k(i-k)2h2,k≠i.


### 2.3. The Shifted Chebyshev Polynomials

The Chebyshev polynomials{Ti(σ);i=0,1,…} are a well-known family of orthogonal polynomials defined on the interval [-1,1] and can be determined with the aid of the recurrence formulae [25, 26]:
(18)Tn+1(σ)=2xTn(σ)-Tn-1(σ),n=1,2,…,T0(σ)=1,T1(σ)=σ.In order to use these polynomials on the intervalt∈[0,τ], it is necessary to define the so-called shifted Chebyshev polynomials by the variable substitution: σ=(2t/τ)-1. Let the shifted Chebyshev polynomials Ti((2t/τ)-1) be denoted by Tτ,i(t). The analytic form of the shifted Chebyshev polynomials Tτ,i(t) is given by
(19)Tτ,i(t)=i∑k=0i(-1)i-k(i+k-1)!22k(i-k)!(2k)!τktk‍,i=1,2,…,Tτ,0(t)=1.
Specially, Tτ,i(0)=(-1)i and Tτ,i(τ)=1.Caputo’s fractional derivative of orderα>0 for the shifted Chebyshev polynomials Tτ,i(t) is given by
(20)D0αCTτ,i(t)=∑k=⌈α⌉ibi,ktk-α‍,i=⌈α⌉,⌈α⌉+1,…,D0αCTτ,i(t)=0,i=0,1,…,⌈α⌉-1,
where
(21)bi,k=i(-1)i-k(i+k-1)!22kk!(i-k)!(2k)!τkΓ(k-α+1).


## 2.1. The Caputo Fractional Derivative

Definition 1 (see [22]).
Letα∈ℝ+. The operator Jaα defined on L1[a,b] by
(4)Jaαf(t)=1Γ(α)∫at(t-s)α-1f(s)ds‍
for a≤t≤b is called the Riemann-Liouville fractional integral operator of order α. For, α=0, we set Ja0:=I, that is, the identity operator.Definition 2 (see [22]).
Letα∈ℝ+ and n=⌈α⌉. The Caputo fractional differential operator DaαC for a≤t≤b is defined as
(5)DaαCf(t)=Jan-αDnf(t)=1Γ(n-α)∫at(t-s)n-α-1f(n)(s)ds‍.


## 2.2. The Composite Translated Sinc Functions

The sinc functions and their properties are discussed in [23, 24]. For any h>0, the translated sinc functions with equidistant space nodes are given as
(6)S(k,h)(z)=sinc(z-khh),k=0,±1,±2,…,
where the sinc functions are defined on the whole real line by
(7)sinc(x)={sin(πx)πx,x≠0,1,x=0.
If f is defined on ℝ, then for any h>0 the series
(8)C(f,h)(z)=∑k=-∞∞f(kh)S(k,h)(z)‍
is called the Whittaker cardinal expansion of f whenever this series converges. f can be approximated by truncating (8).To construct our needed approximations on the interval[a,b], we choose
(9)φ(x)=ln⁡(x-ab-x)
which maps the finite interval [a,b] onto ℝ. The basic functions on [a,b] are taken to be the composite translated sinc functions:
(10)Sφ(k,h)(x)=S(k,h)(φ(x))=sinc(φ(x)-khh).
Thus we may define the inverse image of the equidistant space node {ih} as
(11)xi=φ-1(ih)=a+beih1+eih,i=0,±1,±2,….
The class of functions such that the known exponential convergence rate exists for the sinc interpolation is denoted by B(DE) and defined in the following text.Definition 3 (see [21]).
LetB(DE) be the class of functions f which are analytic in DE and satisfy
(12)∫φ-1(x+L)|f(z)dz|‍⟶0,x⟶±∞,
where L={iυ:|υ|<d≤π/2}, and
(13)∫∂DE|f(z)dz|‍<∞
on the boundary of DE (denoted ∂DE).Theorem 4 (see [21, 23]).
Ifφ′f∈B(DE), then, for all x∈[a,b],
(14)|f(x)-∑k=-∞∞f(xk)Sφ(k,h)(x)‍|≤2N(φ′f)πde-πd/h.
Further, one assumes that there are positive constants C and β so that |f(x)|≤Cexp⁡(-β|φ(x)|). And if one selects h=πd/βN≤2πd/ln⁡2, then,
(15)|dmf(x)dxm-∑k=-NNf(xk)dmdxmSφ(k,h)(x)‍|≤KN(m+1)/2exp⁡(-πdβN)
for all m=0,1,…,n.The above expressions show that the sinc interpolation onB(DE) converges exponentially. We also require the following derivatives of the composite translated sinc functions evaluated at the nodes. Consider
(16)δk,i(0)=[Sφ(k,h)(x)]|x=xi={1,k=i,0,k≠i,(17)δk,i(1)=ddφ[Sφ(k,h)(x)]|x=xi={0,k=i,(-1)i-k(i-k)h,k≠i,δk,i(2)=d2dφ2[Sφ(k,h)(x)]|x=xi={-π23h2,k=i,-2(-1)i-k(i-k)2h2,k≠i.


## 2.3. The Shifted Chebyshev Polynomials

The Chebyshev polynomials{Ti(σ);i=0,1,…} are a well-known family of orthogonal polynomials defined on the interval [-1,1] and can be determined with the aid of the recurrence formulae [25, 26]:
(18)Tn+1(σ)=2xTn(σ)-Tn-1(σ),n=1,2,…,T0(σ)=1,T1(σ)=σ.In order to use these polynomials on the intervalt∈[0,τ], it is necessary to define the so-called shifted Chebyshev polynomials by the variable substitution: σ=(2t/τ)-1. Let the shifted Chebyshev polynomials Ti((2t/τ)-1) be denoted by Tτ,i(t). The analytic form of the shifted Chebyshev polynomials Tτ,i(t) is given by
(19)Tτ,i(t)=i∑k=0i(-1)i-k(i+k-1)!22k(i-k)!(2k)!τktk‍,i=1,2,…,Tτ,0(t)=1.
Specially, Tτ,i(0)=(-1)i and Tτ,i(τ)=1.Caputo’s fractional derivative of orderα>0 for the shifted Chebyshev polynomials Tτ,i(t) is given by
(20)D0αCTτ,i(t)=∑k=⌈α⌉ibi,ktk-α‍,i=⌈α⌉,⌈α⌉+1,…,D0αCTτ,i(t)=0,i=0,1,…,⌈α⌉-1,
where
(21)bi,k=i(-1)i-k(i+k-1)!22kk!(i-k)!(2k)!τkΓ(k-α+1).


## 3. The Derivation of the Sinc-Chebyshev Collocation Method

In order to solve problem (1) with (2) and (3), first of all, we approximate u(x,t) by the 2m+1 composite translated sinc functions and n+1 shifted Chebyshev polynomials as
(22)um,n(x,t)=∑i=-mm∑j=0ncijSφ(i,h)(x)Tτ,j(t)‍‍.
It is noted that the approximate solution um,n(x,t) satisfies the boundary conditions in (3) since Sφ(i,h)(x),i=-m,-m+1,…,m, tend to zeros when x tends to a and b. For discretizing (1) with (2), the lemma is given as follows.Lemma 5.
Let1<α<2 and xk be spatial collocation points given in (11). Then the following relations hold:
(23)∂αum,n(xk,t)∂tα=∑j=2n∑r=2jckjbj,rtr-α‍‍,∂2um,n(xk,t)∂x2=∑i=-mm∑j=0ncijqi,kTτ,j(t)‍‍,∂um,n(xk,t)∂t=∑j=1n∑r=1jckjdj,rtr-1‍‍,
where qi,k=φ′′(xk)δi,k(1)+[φ′(xk)]2δi,k(2) and dj,r=j(-1)j-r(r(j+r-1)!22r/(j-r)!(2r)!τr).Proof.
By (16), (20), and (22), it follows that
(24)∂αum,n(xk,t)∂tα=∑i=-mm‍∑j=0n‍cijSφ(i,h)(xk)D0αCTτ,j(t)=∑i=-mm∑j=2ncijδi,k(0)∑r=2jbj,rtr-α‍‍‍=∑j=2n∑r=2jckjbj,rtr-α‍‍.
Taking into account (17), we obtain
(25)∂2um,n(xk,t)∂x2=∑i=-mm∑j=0ncijd2dx2[Sφ(i,h)(x)]|x=xkTτ,j(t)‍‍=∑i=-mm∑j=0ncij(φ′′(xk)δi,k(1)+[φ′(xk)]2δi,k(2))Tτ,j(t)‍‍=∑i=-mm∑j=0ncijqi,kTτ,j(t).‍‍
Using (16) and (19), one has
(26)∂um,n(xk,t)∂t=∑i=-mm∑j=0ncijSφ(i,h)(xk)ddtTτ,j(t)‍‍=∑i=-mm∑j=1ncijδi,k(0)j∑r=1j(-1)j-rr(j+r-1)!22r(j-r)!(2r)!τrtr-1‍‍‍=∑j=1n∑r=1jckjdj,rtr-1‍‍.
The proof is completed.We are now ready to solve problem (1) with (2) and (3). A collocation scheme is constructed by substituting (22) for u(x,t) into (1) and evaluating the result at the points xk in (11) and tl. For suitable temporal collocation points, we use the roots tl(l=1,2,…,n-1) of the shifted Chebyshev polynomials Tτ,n-1(t). Therefore, using Lemma 5, we have
(27)∑j=2n∑r=2jckjbj,rtlr-α‍‍-a(xk,tl)∑i=-mm∑j=0ncijqikTτ,j(tl)‍‍=f(xk,tl),k=-m,-m+1,…,m,l=1,2,…,n-1.
Also by applying (22) to the initial conditions (2) and collocating in 2m+1 points xk, we obtain
(28)∑j=0n(-1)jckj=ϕ(xk)‍,k=-m,-m+1,…,m,∑j=1n(-1)j-12j2τckj=ψ(xk)‍,k=-m,-m+1,…,m.To obtain a matrix representation of the above equations, we let(29)A=[A1A2A3],B=[B1B2B3],C=[c11,c12,…,c1,n+1,c21,c22,…,c2,n+1,…,c2m+1,,n+1hc2m+1,1,c2m+1,2,…,c2m+1,n+1]T,p=⌊i(2m+1)⌋,s=⌊j(n+1)⌋,v=j-(s-1)(n+1)-1,Δ=a(xi-p(2m+1)-m-1,tp)qs,i-p(2m+1)-m-1Tτ,v(tp),
where
(30)A1=(aij(1))(2m+1)×[(2m+1)(n+1)],A2=(aij(2))(2m+1)×[(2m+1)(n+1)],A3=(aij(3))[(2m+1)(n-1)]×[(2m+1)(n+1)],aij(1)={(-1)j-(i-1)(n+1)-1,1≤j-(i-1)(n+1)≤n+1,0,else,aij(2)={(-1)j-(i-1)(n+1)2[j-(i-1)(n+1)-1]2τ,(-1)j-(i-1)2≤j-(i-1)(n+1)≤n+1,0,hhhhhhhhhhhhhhhhhhhhhhhhhelse,aij(3)={Δ+∑l=2vbv,ltpl-α‍,s=i-p(2m+1),v≥2,Δ,else,B1=[ϕ(x-m),ϕ(x-m+1),…,ϕ(xm)]T,B2=[ψ(x-m),ψ(x-m+1),…,ψ(xm)]T,B3=[f(x-m,t1),f(x-m+1,t1),…,f(xm,t1),hf(x-m,t2),f(x-m+1,t2),…,hf(xm,t2),…,f(x-m,tn-1),hf(x-m+1,tn-1),…,f(xm,tn-1)]T.
So we get a system of (2m+1)(2n+1) linear equations with (2m+1)(2n+1) unknown parameters cij,i=1,2,…,2m+1,j=1,2,…,n+1. And this system can be expressed in a matrix form
(31)AC=B.
Equation (31) can be solved easily for the unknown coefficients cij. Consequently um,n(x,t) given in (22) can be calculated.


## 4. Numerical Examples

To validate the effectiveness of the proposed method for problem (1) with (2) and (3), we consider the example given in [16]. 						(32)∂αu(x,t)∂tα=∂2u(x,t)∂x2+sin(πx),0<x<1,0<t≤1,u(x,0)=0,∂u(x,0)∂t=0,0<x<1,u(0,t)=0,u(1,t)=0,0<t≤1.
The exact solution of the above problem is [27]
(33)u(x,t)=1π2[1-Eα(-π2tα)]sin(πx),
where Eα(z)=∑k=0∞‍zk/Γ(αk+1) is the one-parameter Mittag-Leffler function.To solve the above problem withα=1.7 by using the method described in Section 3, we choose β=1 and d=π/2, and this leads to h=π/2m. We will report the accuracy and efficiency of the method based on the L2-errors and L∞-errors. Figure 1 gives the 3D diagrams of the numerical and exact solutions on the whole computational domain [0,1]×[0,1] with m=15,n=8. A good agreement of the numerical solution with the exact one is achieved. In Table 1, we list the numerical and exact solutions at some points for different numbers of collocation points with n=m. Furthermore, Figure 2 shows the absolute error function |u(x,t)-um,n(x,t)| obtained by the presented method with m=15 and n=8. In Figure 3, we plot the curves of the absolute errors at t=1 for different numbers of collocation points. From Figures 2 and 3, we see that the proposed method can provide accurate results only using a small number of collocation points.Table 1
Some numerical and exact solutions att=1.
m
∖
(
x
,
t
)
(
0.1,1
)
(
0.2,1
)
(
0.3,1
)
(
0.4,1
)
(
0.5,1
)
8
0.042229
0.080693
0.111213
0.130835
0.137592
12
0.042615
0.081166
0.111755
0.131400
0.138168
16
0.042720
0.081293
0.111903
0.131557
0.138330
20
0.042752
0.081332
0.111948
0.131606
0.138379
Exact solution
0.042779
0.081371
0.111997
0.131661
0.138436Figure 1
Comparison of the numerical and exact solution in the domain[0,1]×[0,1].Figure 2
Plot of the absolute error.Figure 3
The error curves att=1.To explore the dependence of errors on the parametersm,n, we represent the L∞-error and L2-error in semi-log scale. Firstly, the computational investigation is concerned with the spatial error. To this end, we fix the polynomial degree n=20, a value large enough such that the error stemming from the temporal approximation is negligible. In Figure 4, we plot the error as functions of m, where a logarithmic scale is used for the spatial-error-axis. As expected, the error shows an exponential decay, since in this semi-log representation one observes that the error variations are approximately linear versus m [19].Figure 4
L
2- and L∞-errors versus m.Now we check the temporal error, which is more interesting because of the fractional derivative in time. For a similar reason mentioned above, we fix a large enough valuem=20 to avoid contamination of the spatial error. We present the error as a function of the shifted Chebyshev polynomial degree n in Figure 5, where a logarithmic scale is now used for the temporal-error-axis. From Figure 5, it is clearly observed that the temporal error depends on the discretization parameters n.Figure 5
L
2- and L∞-errors versus n.


## 5. Conclusion

In this paper, we develop and analyze the efficient numerical methods for the fractional diffusion-wave equation. Based on the collocation technique, the sinc functions and shifted Chebyshev polynomials are used to reduce the problem to the solution of a system of linear algebraic equations. And a matrix representation of the above equations is obtained. In the numerical example, the solution obtained by this method is in excellent agreement with the exact one. The effectiveness and convergence of the presented method are confirmed through the numerical experimentation. One issue of future work is to develop the theory analysis of the method for the proposed fractional differential equation.



---

*Source: 143983-2014-04-01.xml*