# Chebyshev Wavelet Finite Difference Method:  A New Approach for Solving Initial and Boundary Value  Problems of Fractional Order

Authors: A. Kazemi Nasab; A. Kılıçman; Z. Pashazadeh Atabakan; S. Abbasbandy
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

A new method based on a hybrid of Chebyshev wavelets and finite difference methods is introduced for solving linear and nonlinear fractional differential equations. The useful properties of the Chebyshev wavelets and finite difference method are utilized to reduce the computation of the problem to a set of linear or nonlinear algebraic equations. This method can be considered as a nonuniform finite difference method. Some examples are given to verify and illustrate the efficiency and simplicity of the proposed method.


---

## Body

## 1. Introduction

The study of fractional calculus dates back to 17th century, starting by G. W. Leibnitz (1695, 1697) and L. Euler (1730) [1, 2] and then has been developed by many researchers in different disciplines. In the year 1823, Liouville and Abel introduced the theory of fractional derivatives and integrals; for more details, please refer to [3, 4]. Fractional calculus has received much more attention from scientists and engineers in recent years. Many researchers in various fields found that derivatives of noninteger order are useful for the description of some natural physics phenomena and dynamic system processes such as damping laws and diffusion process [5, 6]. In general, it is difficult to solve fractional differential equations analytically. Therefore, it is necessary to introduce some reliable and efficient numerical algorithms to solve them. During the past decades, an increasing number of numerical methods were being developed. These methods include homotopy analysis method [7], homotopy perturbation method [8–10], variational iteration method [9–13], finite difference method [5, 14–18], Adomian decomposition method [19–23], fractional differential transform method [24, 25], predictor-corrector method [26], fractional linear multistep method [27], extrapolation method [28], integral transform [29], and generalized block pulse operational matrix method [30, 31].In recent years, wavelets have received considerable attention by researchers in different fields of science and engineering. One advantage of wavelet analysis is the ability to perform local analysis [32]. Wavelet analysis is able to reveal signal aspects that other analysis methods miss, such as trends, breakdown points, and discontinuities. In comparison with other orthogonal functions, multiresolution analysis aspect of wavelets permits the accurate representation of a variety of functions and operators. In other words, we can change M and k simultaneously to get more accurate solution. Another benefit of wavelet method for solving equations is that after discreting the coefficients matrix of algebraic equations is sparse. So the use of wavelet methods for solving equations is computationally efficient. In addition, the solution is convergent. The operational matrix of fractional order integration for the Chebyshev wavelet, Legendre wavelet, and Haar wavelet has been introduced in [33–35] to solve the differential equations of fractional order. A CAS wavelet operational matrix of fractional order integration has been developed by Saeedi et al., to solve fractional nonlinear integrodifferential equations [36, 37].The paper is organized as follows. Section2 included some necessary definitions and mathematical preliminaries of fractional calculus, Chebyshev polynomials, and Chebyshev wavelets. In Section 3, we introduce the Chebyshev finite difference method. In Section 4, Chebyshev wavelet finite difference method (CWFD) is presented. Section 5 included convergence analysis of the proposed method. In Section 6, the proposed approach is used to approximate the fractional differential equations. As a result the fractional differential equation is converted to a system of algebraic equations which is simply solved. Some illustrative examples of different types are given to demonstrate the efficiency and accuracy of the method in Section 7. In Section 8, concluding remarks are given.


## 2. Preliminaries and Notations

In this section, we present some notations, definitions, and preliminary facts that will be used further in this work.


### 2.1. Fractional Integral and Derivative

There are several definitions of fractional integral and derivatives [2, 38, 39], including Riemann-Liouville, Caputo, Weyl, Hadamard, Marchaud, Riesz, Grunwald-Letnikov, and Erdelyi-Kober. The most commonly used definition is of Riemann-Liouville and Caputo. One of the drawbacks of Riemann-Liouville method is that it cannot incorporate the nonzero initial condition at lower limit. The Caputo fractional derivative allows the utilization of initial and boundary conditions involving integer order derivatives, which have clear physical interpretations. Therefore, in this study we will use the Caputo fractional derivative Dα proposed by Podlubny [14].Definition 1.
A real functionf(x), x>0, is said to be in the space Cμ, μ∈ℝ if there exists a real number p>μ such that f(x)=xpf1(x), where f1(x)∈C[0,∞), and it is said to be in the space Cμn as if and only if fn(x)∈Cμ, n∈ℕ.Definition 2.
The Riemann-Liouville fractional integration of orderα≥0 of a function f∈Cμ, μ≥-1 is defined as [14]
(1)(Iαf)(x)=1Γ(α)∫0x‍(x-τ)α-1f(τ)dτ,(I0f)(x)=f(x).Definition 3.
The fractional derivative off(x) in the Caputo sense is defined as [14]
(2)(Dαf)(x)=(In-αf(n))(x)=1Γ(n-α)∫0x(x-τ)n-α-1f(n)(τ)dτ,
for n-1<α≤n, n∈ℕ, x>0, f∈C-1n.
Some basic properties of the fractional operator are as follows [14, 40], for f∈Cμn, μ≥-1,α,β≥0, and γ≥-1 we have(1)
IαIβf(x)=Iα+βf(x),(2)
IαIβf(x)=IβIαf(x),(3)
Iαxγ=(Γ(γ+1)/Γ(γ+α+1))xγ+α,(4)
DαIαf(x)=f(x),(5)
IαDαf(x)=f(x)-∑k=0n-1‍f(k)(0+)(xk/k!),x>0,n-1<α≤n,(6)
(3)Dαxk={0,k≤⌈α⌉,k∈ℤ+Γ(k+1)Γ(k+1-α)xk-α,k≥⌈α⌉,k∈ℤ+,where⌈α⌉ denotes the smallest integer greater than or equal to α.


### 2.2.  Chebyshev Polynomials

Chebyshev polynomials of the first kind of degreem can be defined as follows:
(4)Tm(x)=cosmα,α=arccosx,
which are orthogonal with respect to the weight function w(x)=1/1-x2 [41],
(5)∫-11‍Tm(x)Tn(x)w(x)dx=π2cmδnm,
where δnm is the Kronecker delta function and
(6)cm={2,m=01,m≥1.Chebyshev polynomials also satisfy the following recursive formula:(7)T0(x)=1,T1(x)=x,Tm+1(x)=2xTm(x)-Tm-1(x),m=1,2,….The set of Chebyshev polynomials is a complete orthogonal set in the Hilbert space£w2[-1,1]. A function f∈£w2[-1,1] can be written in terms of Chebyshev polynomials as
(8)f(x)=∑m=0∞‍f^mTm(x),f^m=2πcm∫-11‍f(x)Tm(x)w(x)dx.


### 2.3. Wavelets and Chebyshev Wavelets

Wavelets have been very successfully used in many scientific and engineering fields. They constitute a family of functions constructed from dilation and transformation of a single function called the mother waveletψ(x); we have the following family of continuous wavelets as [42, 43]:
(9)ψa,b(x)=|a|-1/2ψ(x-ba),a,b∈R,a≠0.If we seta=a0-k, b=nb0a0-k; a0>1, b0>0, we will get the following family of discrete wavelet which forms a wavelet basis for L2(ℝ):
(10)ψa,b(x)=|a0k/2|ψ(a0kx-nb0),a,b∈ℤ.In particular, whena0=2 and b0=1, then ψa,b(x) forms an orthonormal basis.Chebyshev waveletsψn,m=ψ(k,n,m,t) have four arguments; n=1,…,2k-1,k can assume any positive integer, m is degree of Chebyshev polynomials of the first kind, and t denotes the time. Consider
(11)ψn,m(x)={2k/2pmTm(2kx-2n+1),n-12k-1≤x<n2k-10,otherwise,
where
(12)pm={1π,m=0,2π,m≥1,
and m=0,1,…,M and n=1,…,2k-1. In (11) the coefficients are used for orthonormality. We should note that in dealing with the Chebyshev wavelets, the weight function w~(x)=w(2x-1) has to be dilated and translated to get orthogonal wavelets as follows:
(13)wn(x)=w(2kx-2n+1).In view of (5), Chebyshev wavelets are an orthonormal set with respect to the weight function wn(x) because
(14)〈ψn,m(x),ψr,s(x)〉wn=π2pmpscmδnrδms.Lemma 4.
The family of Chebyshev wavelets{ψn,m(x)∣n=1,2,…,2k-1,m∈ℕ∪{0}} forms an orthonormal basis for L2[0,1] with respect to the weight function wn(x) [44].
Fork=2 and M=3, Chebyshev wavelets are as follows:
(15)ψ1,0=2πψ1,1=22π(4x-1)ψ1,2=22π[2(4x-1)2-1]ψ1,3=22π(256x3-192x2+36x-1)}0≤x<12ψ2,0=2πψ2,1=22π(4x-3)ψ2,2=22π[2(4x-3)2-1]ψ2,3=22π(256x3-576x2+420x-99)}12≤x<1.


### 2.4. Function Approximation

A functionf(t) defined over [0,1) may be expanded as
(16)f(x)=∑n=1∞‍∑m=0∞‍cn,mψn,m(x),
where cn,m=(f(x),ψn,m(x))wn, in which (·,·) denotes the inner product in Lwn2[0,1]. If we consider truncated series in (16), we obtain
(17)f(x)=∑n=12k-1‍∑m=0M‍cn,mψn,m(x)=CTΨ(x),
where C and Ψ(x) are 2k-1(M+1)×1 matrices given by
(18)C=[c1,0,c1,1,…,c1,M,c2,0,c2,1,…,c2,M,…,c2k-1,0,c2k-1,1,…,c2k-1,M]T,Ψ(x)=[ψ1,0,…,ψ1,M,ψ2,0,…,ψ2,M,…,ψ2k-1,0,ψ2k-1,1,…,ψ2k-1,M]T.


## 2.1. Fractional Integral and Derivative

There are several definitions of fractional integral and derivatives [2, 38, 39], including Riemann-Liouville, Caputo, Weyl, Hadamard, Marchaud, Riesz, Grunwald-Letnikov, and Erdelyi-Kober. The most commonly used definition is of Riemann-Liouville and Caputo. One of the drawbacks of Riemann-Liouville method is that it cannot incorporate the nonzero initial condition at lower limit. The Caputo fractional derivative allows the utilization of initial and boundary conditions involving integer order derivatives, which have clear physical interpretations. Therefore, in this study we will use the Caputo fractional derivative Dα proposed by Podlubny [14].Definition 1.
A real functionf(x), x>0, is said to be in the space Cμ, μ∈ℝ if there exists a real number p>μ such that f(x)=xpf1(x), where f1(x)∈C[0,∞), and it is said to be in the space Cμn as if and only if fn(x)∈Cμ, n∈ℕ.Definition 2.
The Riemann-Liouville fractional integration of orderα≥0 of a function f∈Cμ, μ≥-1 is defined as [14]
(1)(Iαf)(x)=1Γ(α)∫0x‍(x-τ)α-1f(τ)dτ,(I0f)(x)=f(x).Definition 3.
The fractional derivative off(x) in the Caputo sense is defined as [14]
(2)(Dαf)(x)=(In-αf(n))(x)=1Γ(n-α)∫0x(x-τ)n-α-1f(n)(τ)dτ,
for n-1<α≤n, n∈ℕ, x>0, f∈C-1n.
Some basic properties of the fractional operator are as follows [14, 40], for f∈Cμn, μ≥-1,α,β≥0, and γ≥-1 we have(1)
IαIβf(x)=Iα+βf(x),(2)
IαIβf(x)=IβIαf(x),(3)
Iαxγ=(Γ(γ+1)/Γ(γ+α+1))xγ+α,(4)
DαIαf(x)=f(x),(5)
IαDαf(x)=f(x)-∑k=0n-1‍f(k)(0+)(xk/k!),x>0,n-1<α≤n,(6)
(3)Dαxk={0,k≤⌈α⌉,k∈ℤ+Γ(k+1)Γ(k+1-α)xk-α,k≥⌈α⌉,k∈ℤ+,where⌈α⌉ denotes the smallest integer greater than or equal to α.


## 2.2.  Chebyshev Polynomials

Chebyshev polynomials of the first kind of degreem can be defined as follows:
(4)Tm(x)=cosmα,α=arccosx,
which are orthogonal with respect to the weight function w(x)=1/1-x2 [41],
(5)∫-11‍Tm(x)Tn(x)w(x)dx=π2cmδnm,
where δnm is the Kronecker delta function and
(6)cm={2,m=01,m≥1.Chebyshev polynomials also satisfy the following recursive formula:(7)T0(x)=1,T1(x)=x,Tm+1(x)=2xTm(x)-Tm-1(x),m=1,2,….The set of Chebyshev polynomials is a complete orthogonal set in the Hilbert space£w2[-1,1]. A function f∈£w2[-1,1] can be written in terms of Chebyshev polynomials as
(8)f(x)=∑m=0∞‍f^mTm(x),f^m=2πcm∫-11‍f(x)Tm(x)w(x)dx.


## 2.3. Wavelets and Chebyshev Wavelets

Wavelets have been very successfully used in many scientific and engineering fields. They constitute a family of functions constructed from dilation and transformation of a single function called the mother waveletψ(x); we have the following family of continuous wavelets as [42, 43]:
(9)ψa,b(x)=|a|-1/2ψ(x-ba),a,b∈R,a≠0.If we seta=a0-k, b=nb0a0-k; a0>1, b0>0, we will get the following family of discrete wavelet which forms a wavelet basis for L2(ℝ):
(10)ψa,b(x)=|a0k/2|ψ(a0kx-nb0),a,b∈ℤ.In particular, whena0=2 and b0=1, then ψa,b(x) forms an orthonormal basis.Chebyshev waveletsψn,m=ψ(k,n,m,t) have four arguments; n=1,…,2k-1,k can assume any positive integer, m is degree of Chebyshev polynomials of the first kind, and t denotes the time. Consider
(11)ψn,m(x)={2k/2pmTm(2kx-2n+1),n-12k-1≤x<n2k-10,otherwise,
where
(12)pm={1π,m=0,2π,m≥1,
and m=0,1,…,M and n=1,…,2k-1. In (11) the coefficients are used for orthonormality. We should note that in dealing with the Chebyshev wavelets, the weight function w~(x)=w(2x-1) has to be dilated and translated to get orthogonal wavelets as follows:
(13)wn(x)=w(2kx-2n+1).In view of (5), Chebyshev wavelets are an orthonormal set with respect to the weight function wn(x) because
(14)〈ψn,m(x),ψr,s(x)〉wn=π2pmpscmδnrδms.Lemma 4.
The family of Chebyshev wavelets{ψn,m(x)∣n=1,2,…,2k-1,m∈ℕ∪{0}} forms an orthonormal basis for L2[0,1] with respect to the weight function wn(x) [44].
Fork=2 and M=3, Chebyshev wavelets are as follows:
(15)ψ1,0=2πψ1,1=22π(4x-1)ψ1,2=22π[2(4x-1)2-1]ψ1,3=22π(256x3-192x2+36x-1)}0≤x<12ψ2,0=2πψ2,1=22π(4x-3)ψ2,2=22π[2(4x-3)2-1]ψ2,3=22π(256x3-576x2+420x-99)}12≤x<1.


## 2.4. Function Approximation

A functionf(t) defined over [0,1) may be expanded as
(16)f(x)=∑n=1∞‍∑m=0∞‍cn,mψn,m(x),
where cn,m=(f(x),ψn,m(x))wn, in which (·,·) denotes the inner product in Lwn2[0,1]. If we consider truncated series in (16), we obtain
(17)f(x)=∑n=12k-1‍∑m=0M‍cn,mψn,m(x)=CTΨ(x),
where C and Ψ(x) are 2k-1(M+1)×1 matrices given by
(18)C=[c1,0,c1,1,…,c1,M,c2,0,c2,1,…,c2,M,…,c2k-1,0,c2k-1,1,…,c2k-1,M]T,Ψ(x)=[ψ1,0,…,ψ1,M,ψ2,0,…,ψ2,M,…,ψ2k-1,0,ψ2k-1,1,…,ψ2k-1,M]T.


## 3. Chebyshev Finite Difference Method

Clenshaw and Curtis [45] introduced the following approximation of the function f(x) denoted by (PM)f(x) as
(19)(PM)f(x)=∑m=0M′′fmTm(x),fm=2M∑k=0M′′f(xk)Tm(xk),
where the summation symbol with double primes denotes a sum with both the first and last terms halved. Moreover, xm are the extrema of the Mth-order Chebyshev polynomial TM(x) and are defined as
(20)xm=cos(mπM),m=0,1,2,…,M.These well-known Chebyshev-Gauss-Lobatto interpolated points,(21)xM=-1<xM-1<⋯<x1<x0=1,
are the zeros of (1-x2)(dTM(x)/dx). Using (4) we have,
(22)Tm(xk)=cos(mkπM),
so fm can be rewritten as
(23)fm=2M∑k=0M∑′′f(xk)cos(mkπM).The first two derivatives of the functionf(x) at the points xm,m=0,1,…,M are given by [46] as
(24)f(n)(xm)=∑j=0M‍dm,j(n)f(xj),n=1,2,
where
(25)dm,j(1)=4θjM∑k=1M‍∑l=0,(k+l)oddk-1‍kθkclTk(xj)Tl(xm)=4θjM∑k=1M‍∑l=0,(k+l)oddk-1‍kθkcl×cos(kjπM)cos(lmπM),dm,j(2)=2θjM∑k=2M‍∑l=0,(k+l)evenk-2‍k(k2-l2)θkclTk(xj)Tl(xm)=2θjM∑k=2M‍∑l=0,(k+l)evenk-2‍k(k2-l2)θkclcos(kjπM)cos(lmπM),
with θ0=θM=1/2, θj=1 for  j=1,2,…,M-1.As can be seen from (24), the first two derivatives of the function f(x) at any point of the Chebyshev-Gauss-Lobatto points are expanded as a linear combination of the values of the function as these points.


## 4. Chebyshev Wavelet Finite Difference Method

In this section, we present the Chebyshev wavelet finite difference (CWFD) method. Considerxnm, n=1,2,…,2k-1, m=0,1,…,M, as the corresponding Chebyshev-Gauss-Lobatto collocation points at the nth subinterval [(n-1)/2(k-1),n/2(k-1)] such that
(26)xnm=12k(xm+2n-1).A functionf(x) can be written in terms of Chebyshev wavelet basis functions as follows:
(27)(PM)f(x)=∑n=12k-1∑m=0M′′cnmψnm(x),
where cnm, n=1,2,…,2k-1, m=0,1,…,M, are the expansion coefficients of the function f(x) at the subinterval [(n-1)/2k-1,n/2k-1] and ψn,m(x), n=1,2,…,2k-1, m=0,1,…,M, are defined in (11).In view of (11) and (19), we can obtain the coefficients cnm as
(28)cnm=1(2k/2pm)2·2M∑p=0M′′∑f(xnp)ψnm(xnp)=12k/2pm·2M∑p=0M′′f(xnp)cos(mpπM).Using (25), the first two derivatives of the function f(x) at the points xnm, n=1,2,…,2k-1, m=0,1,…,M, can be obtained as
(29)f(n)(xnm)=∑j=0M‍dn,m,j(n)f(xnj),n=1,2,
where
(30)dn,m,j(1)=4θjM∑k=1M‍∑l=0,(k+l)oddk-1‍kθkclpkplψnk(xnj)ψnl(xnm)=4θjM∑k=1M‍∑l=0,(k+l)oddk-1‍2kkθkclcos(kjπM)cos(lmπM),dm,j(2)=2θjM∑k=2M‍∑l=0,(k+l)evenk-2‍2kk(k2-l2)θkclpkpl222222222222222222×ψnk(xnj)ψnl(xnm)=2θjM∑k=2M‍∑l=0,(k+l)evenk-2‍4kk(k2-l2)θkcl×cos(kjπM)cos(lmπM).


## 5. Convergence Analysis

Lemma 5.
If the Chebyshev wavelet expansion of a continuous functionf(x) converges uniformly, then the Chebyshev wavelet expansion converges to the function f(x) [47].Theorem 6.
A functionf(x)∈Lwn2[0,1), with bounded second derivative, say |f′′(x)|≤B, can be expanded as an infinite sum of Chebyshev wavelets, and the series converges uniformly to f(x); that is,
(31)f(x)=∑n=1∞‍∑m=0∞‍c^nmψnm(x).Proof.
We have(32)c^n,m=(f(x),ψn,m(x))wn,k=∫01‍f(x)ψnm(x)wn,k(x)dx=∫(n-1)/2(k-1)n/2(k-1)‍2k/2pmf(x)Tm(2kx-2n+1)w×(2kx-2n+1)dx,
if m>1, by substituting 2kx-2n+1=cosα, it yields
(33)c^nm=12k/2∫0π‍f(cosα+2n-12k)2πcosmαdα.
Using integration by part, we get
(34)c^nm=22k/2πf(cosα+2n-12k)(sinmαm)|0π+223k2mπ∫0π‍f′(cosα+2n-12k)sinmαsinαdα;
the first part is zero, therefore,
(35)c^nm=223k/2mπ∫0π‍f′(cosα+2n-12k)sinmαsinαdα.
Using integration by part again, it yields(36)c^nm=123k/2m2πf′(cosα+2n-12k)×(sin(m-1)αm-1-sin(m+1)αm+1)|0π+125k/2m2π∫0π‍f′′(cosα+2n-12k)rm(α)dα,
where
(37)rm(α)=sinα(sin(m-1)αm-1-sin(m+1)αm+1).
Thus, we get(38)|c^nm|=|125k/2m2π∫0π‍f′′(cosα+2n-12k)rm(α)dα|≤(125k/2m2π)∫0π‍|f′′(cosα+2n-12k)rm(α)|dα≤B25k/2m2π∫0π‍|rm(α)|dα.
However(39)∫0π‍|rm(α)|dα=∫0π‍|sinα(sin(m-1)αm-1-sin(m+1)αm+1)|dα≤∫0π‍|sinαsin(m-1)αm-1|+|sinαsin(m+1)αm+1|dα≤2mπm2-1.
Sincen≤2k-1, we obtain
(40)|c^nm|≤2πB(2n)5/2(m2-1).
Now, ifm=1, by using (35), we have
(41)|c^n1|<2π(2n)3/2max0≤x≤1|f′(x)|.
It is mentioned in [42] that {ψn0}n=1∞ form an orthogonal system constructed by Haar scaling function with respect to the weight function w(x), so ∑n=1∞‍(2k/√π)c^n0ψn0(x) is convergent. Hence, we will have
(42)|∑n=1∞‍∑m=0∞‍c^nmψnm(x)|≤|2kπ∑n=1∞‍c^n0ψn0(x)|+∑n=1∞‍∑m=1∞‍|c^nm||ψnm(x)|≤|2kπ∑n=1∞‍c^n0ψn0(x)|+2(k+1)2π∑n=1∞‍∑m=1∞‍|c^nm|<∞.
Therefore, with the aid of Lemma5, the series ∑n=1∞‍∑m=1∞‍c^nmψnm(x) converges to f(x) uniformly [47].Theorem 7.
Supposef(x)∈Lwn2[0,1) with bounded second derivative, say |f′′(t)|≤B; then its Chebyshev wavelet finite difference expansion converges uniformly to f(x); that is,
(43)∑n=12k-1∑m=0∞′cnmψnm(x)=f(x),
where the summation symbol with prime denotes a sum with the first term halved.Proof.
From Theorem6, we have
(44)f(x)=∑n=12k-1‍∑m=0∞‍c^nmψnm(x),
where
(45)c^n,m=(f(x),ψn,m(x))wn.
This series converges tof(x) uniformly. We first show that cnm converges to c^nm. We have
(46)c^n,m=(f(x),ψn,m(x))wn,k=∫01‍f(x)ψnm(x)wn,k(x)dx=∫(n-1)/2(k-1)n/2(k-1)‍2k/2pmf(x)Tm(2kx-2n+1)×w(2kx-2n+1)dx,
by substituting 2kx-2n+1=cosα, it yields
(47)c^nm=pm2k/2∫0π‍f(cosα+2n-12k)cosmαdα.
Using the trapezoidal rule for integration withM equidistant subintervals gives
(48)c^nm≃pm2k/2·πM∑j=0M′′f(cos(αj)+2n-12k)×cos(mαj),αj=jπM.
From (26) and (28), for 1<m<M, we have
(49)c^nm≃12k/2pm·2M∑j=0M′′f(xnj)cos(mjπM)=cnm.
According to approximation error for the trapezoidal rule, we have(50)|c^nm-cnm|≤Rπ312M2,
where
(51)R=max{d2dα2(f(cos(α)+2n-12k)cos(mα)),0≤α≤πd2dα2}.
In view of (35) and triangle inequality, we get
(52)|cnm|≤|cnm-c^nm|+|c^nm|≤Rπ312M2+2πB(2n)5/2(m2-1)≤Rπ312(m2-1)+2πB(2n)5/2(m2-1)≤C(m2-1),
where
(53)C=Rπ312+2πB(2n)5/2.
Because |ψnm(x)|≤2k/2pm, it is understandable that
(54)|∑n=12k-1‍∑m=0∞′cnmψnm(x)|≤|∑n=12k-1‍12cn0ψn0(x)|+|∑n=12k-1‍cn1ψn1(x)|+∑n=12k-1‍∑m=2∞‍|cnm||ψnm(x)|≤2k/22π∑n=12k-1‍|cn0|+2(k+1)/2π∑n=12k-1‍|cn1|+2(k+1)/2π×∑n=12k-1‍∑m=2∞‍|cnm|<∞.
Therefore, in view of Lemma5, series (43) is uniformly convergent to f(x).Theorem 8 (accuracy estimation).
Supposef(x)∈Lwn2[0,1) with bounded second derivative, say |f′′(x)|≤B, then one has the following accuracy estimation:
(55)σk,M≤(∑n=12k-1‍∑m=M+1∞‍C2(m2-1)2)1/2,
where
(56)σk,M=(∫01‍[f(x)-∑n=12k-1‍∑m=0M‍cnmψn,m(x)]2wn(x)dx)1/2,C=Rπ312+2πB(2n)5/2.Proof .
We have(57)σk,M2=∫01‍[f(x)-∑n=12k-1‍∑m=0M‍cnmψn,m(x)]2wn(x)dx=∫01‍[∑n=12k-1‍∑m=0M‍∑n=12k-1‍∑m=0∞‍cnmψn,m(x)-∑n=12k-1‍∑m=0M‍cnmψn,m(x)]2wn(x)dx=∫01‍∑n=12k-1‍∑m=M+1∞‍cnm2(ψn,m(x))2wn(x)dx=∑n=12k-1‍∑m=M+1∞‍cnm2∫01‍(ψn,m(x))2wn(x)dx.
We know that the family{ψn,m(x)∣n=1,2,…,2k-1,m∈ℕ∪{0}} forms orthonormal basis for Lwn2[0,1), so ∫01‍(ψn,m(x))2wn(x)dx=1. Therefore, in view of (52), we will have
(58)σk,M2=∑n=12k-1‍∑m=M+1∞‍cnm2≤∑n=12k-1‍∑m=M+1∞‍C2(m2-1)2,
where
(59)C=Rπ312+2πB(2n)5/2.


## 6. Discretization of Problem

In this section, the Chebyshev wavelet finite difference method (CWFD) is used for solving the following general form:(60)H(t,y(t),Dβ1y(t),…,Dβky(t),…,Dβpy(t))=0,aaaaaaaaaaaaaaaaaaaaaaaaaat∈I=[0,L],
subject to the conditions
(61)Bq(y(ρ0),…,y(l)(ρ0),y(ρ1),…,y(l)(ρ1),…,y(ρl),…,y(l)(ρl)),q=0,1,…,l,
where 0<β1<β2<⋯<βp, l<βp≤l+1, ρj, j=0,1,…,l are located in [0,L] and H can be linear or nonlinear while Bq are linear functions.We suppose the interval[0,1] is divided into 2k-1 subintervals In=[(n-1)/2k-1,n/2k-1], n=1,2,…,2k-1. We also consider the shifted Chebyshev-Gauss-Lobatto collocation points (n-1)/2k-1=tn0<tn1<⋯<tn,M-1<tn,M=n/2k-1 on the nth subinterval In, n=1,2,…,2k-1, where tns is defined as follows:
(62)tns=12k(ts+2n-1),s=1,2,…,M-l.In order to obtain the solutiony(t) in (60), we first approximate y(t) according to (27) and rewrite Dβiy(t), i=1,…,p as fractional derivatives of y(t), in the Caputo sense using (2). Collocating (60) at the shifted Chebyshev-Gauss-Lobatto points tns, n=1,2,…,2k-1, s=0,1,…,M-l, we get
(63)H(tns,y(tns),Dβ1y(tns),…,Dβky(tns),…,Dβpy(tns))=0.
We then continue as follows:
(64)Dβiy(tns)=1Γ(ni-βi)∫0tns‍(tns-τ)ni-βi-1y(ni)(τ)dτ11111111111111ni=[βi]+1,i=1,…,p=1Γ(ni-βi)∫0(n-1)/2(k-1)(tns-τ)ni-βi-1y(ni)(τ)dτ+1Γ(ni-βi)∫(n-1)/2(k-1)tns(tns-τ)ni-βi-1y(ni)(τ)dτ.To calculate the first integral in the above summation, we use the Clenshaw-Curtis quadrature formula [48]:
(65)∫-11‍f(x)dx=∑r=0N‍wrf(xr),
where xr are Chebyshev-Gauss-Lobatto nodes and the weights wr are given by
(66)w0=wN={1N2,Nodd,1N2-1,Neven,wr=2Nγr[1-∑k=1[N2]‍2γ2k(4k2-1)cos2krπN],2222222222222222r=1,2,…,N-1,
where γ0=γN=2 and γr=1, for r=1,2,…,N-1.In this paper we setN=M. We obtain
(67)1Γ(ni-βi)∫0(n-1)/2k-1‍(tns-τ)ni-βi-1y(ni)(τ)dτ=1Γ(ni-βi)×∑r1=1n-1‍∫(r1-1)/2k-1r1/2k-1‍(tns-τ)ni-βi-1y(ni)(τ)dτ=1Γ(ni-βi)·12k∑r1=1n-1‍∑r2=0N‍wr2(tns-tr1r2)ni-βi-1y(ni)(tr1r2).Using integration by part and in view of (62), we convert second singular integral to nonsingular one as follows: get the second integral
(68)1Γ(ni-βi)∫(n-1)/2k-1tns‍(tns-τ)ni-βi-1y(ni)(τ)dτ=(tns-((n-1)/2k-1))ni-βi(ni-βi)Γ(ni-βi)y(ni)(n-12k-1)+(2k-1tns-n+1)2k(ni-βi)Γ(ni-βi),∑s1=0N‍(tns-(tns-n-12k-1)(2k-1tns1-n+1)-n-12k-1)ni-βiy(ni+1)((tns-n-12k-1)(2k-1tns1-n+1)+n-12k-1).Replacing (67) and (68) into (64), we obtain Dβiy(tns), i=1,2,…,p and then replace them into (64) to get 2k-1(M-l) equations.Furthermore, substituting (27) and (29) into (61), we get (l+1) equations.Moreover, we should impose continuity condition on the approximate solution and its firstl derivatives at the interface between subintervals which results in (2k-1-1)(l+1) equations. 						(69)y(r)(tn+1,M)=y(r)(tn,0),r=0,1,…,l,n=1,2,…,2k-1-1.We will totally have a system of2k-1(M+1) algebraic equations, which can be solved for the y(tns). Consequently, we obtain the solution y(t) to the given (60) using (27) and (28).


## 7. Illustrative Examples

In this section, we consider some numerical examples for the fractional equation to demonstrate the validity of the proposed method (CWFD) in solving fractional differential equation. These examples are considered because closed form solutions are available for them, and they have been solved using other numerical methods. This allows one to compare the results obtained using this method with the analytical solution and the solutions obtained using other methods.Example 1.
As the first example, we consider a nonlinear equation defined as follows [49, 50]:
(70)Dαy(t)=40320Γ(9-α)t8-α-3Γ(5+α/2)Γ(5-α/2)t4-(α/2)+94Γ(α+1)+(32t(α/2)-t4)3-[y(t)](3/2),
such that
(71)y(0)=0,y′(0)=0,
the second initial condition is for α>1 only. The exact solution of (70) and (71) is given as [26]
(72)y(t)=t8-3t4+α/2+94tα.
It should be noted that forα<1, the slope of the solution at t=0 goes to infinity. Therefore, one can expect a large numerical error near t=0.
We applied the method introduced in Section6 and solved this problem for different values of α and M, k. Figure 1 shows the analytical and numerical results for M=9,6,3 and k=1 and α=0.75,1.5. It can be seen that increasing the values of M and k results in more accurate solution. In Figure 2 numerical results for M=12, k=4, and α=0.5,0.75,0.95 and α=1 (exact solution) and also α=1.5,1.75,1.95, and α=2 (exact solution) are plotted. We see that as α approaches 1 and 2, the solution of the fractional differential equation approaches to that of the integer-order differential equation. In Table 1, we compare the results obtained by the present method using M=12, k=6 with those in [50]. As can be seen, our results are much more better than those obtained by Saadatmandi and Dehghan [50]. Furthermore, as it is expected, the absolute error is reduced as α approaches an integer value.Table 1
Comparison of absolute error iny(t) in [50] and obtained using our method for Example 1.
α
t
=
0.1
t
=
0.3
t
=
0.7
t
=
0.9
Our result
Ref. [50]
Our result
Ref. [50]
Our result
Ref. [50]
Our result
Ref. [50]
0.2
8.6
×
1
0
-
5
2.2
×
1
0
-
1
2.6
×
1
0
-
5
2.3
×
1
0
-
1
1.1
×
1
0
-
5
5.3
×
1
0
-
1
1.0
×
1
0
-
5
1.7
×
1
0
0
0.4
1.7
×
1
0
-
4
6.3
×
1
0
-
2
5.3
×
1
0
-
5
6.0
×
1
0
-
2
2.0
×
1
0
-
5
1.2
×
1
0
-
1
1.8
×
1
0
-
5
3.0
×
1
0
-
1
0.6
1.9
×
1
0
-
4
1.5
×
1
0
-
2
7.5
×
1
0
-
5
1.3
×
1
0
-
2
3.0
×
1
0
-
5
2.1
×
1
0
-
2
2.7
×
1
0
-
5
3.7
×
1
0
-
2
0.8
1.2
×
1
0
-
4
2.9
×
1
0
-
3
7.1
×
1
0
-
5
2.1
×
1
0
-
3
4.8
×
1
0
-
5
2.5
×
1
0
-
3
5.2
×
1
0
-
5
2.1
×
1
0
-
3
1.2
5.1
×
1
0
-
3
1.9
×
1
0
-
3
1.8
×
1
0
-
2
1.6
×
1
0
-
3
3.5
×
1
0
-
2
2.9
×
1
0
-
3
4.0
×
1
0
-
2
1.6
×
1
0
-
2
1.4
2.8
×
1
0
-
4
2.0
×
1
0
-
4
1.3
×
1
0
-
3
1.6
×
1
0
-
3
2.8
×
1
0
-
3
4.9
×
1
0
-
3
3.3
×
1
0
-
3
3.3
×
1
0
-
2
1.6
4.3
×
1
0
-
5
6.3
×
1
0
-
5
3.2
×
1
0
-
5
7.3
×
1
0
-
4
2.9
×
1
0
-
5
2.3
×
1
0
-
3
6.0
×
1
0
-
5
1.3
×
1
0
-
2
1.8
3.0
×
1
0
-
5
3.8
×
1
0
-
5
7.3
×
1
0
-
5
2.0
×
1
0
-
4
1.5
×
1
0
-
4
5.9
×
1
0
-
4
1.9
×
1
0
-
4
2.8
×
1
0
-
3Comparison ofy(t) for M=9,6,3 and k=2 and (a) α=0.75, (b) α=1.5, with exact solution, for Example 1.
(a)
(b)Comparison ofy(t) for M=12, k=4, and (a) α=0.5,0.75,0.95,1, (b) α=1.5,1.75,1.95,2 for Example 1.
(a)
(b)Example 2.
As the second example, we consider the following initial value problem in the case of the inhomogeneous Bagley-Torvik equation [51]:
(73)D2y(t)+D1.5y(t)+y(t)=g(t),
where g(t)=1+t subject to the following initial states
(74)y(0)=1,y′(0)=1.
The exact solution of this problem isy(t)=1+t.
We solve this fractional initial value problem by applying the method described in Section6 with M=k=3. The absolute error is shown in Figure 3, which shows that the numerical solution is in perfect agreement with the exact solution. We set Digits=50. However, we get the absolute error less or equal to 10-18 when we set Digits=20.Figure 3
Absolute error for Example2.Example 3.
Following El-Mesiry et al. [52] and Li [33], we consider the following nonlinear fractional differential equation:
(75)aD2.0y(t)+bDα2y(t)+cDα1y(t)+e[y(t)]3=f(t),0<α1≤1,1<α2≤2,
where
(76)f(t)=2aΓ(2)t+2bΓ(4-α2)t3-α2+2cΓ(4-α1)t3-α1+e[13t3]3,
subject to
(77)y(0)=y′(0)=0.
The exact solution of this problem isy(t)=(1/3)t3.
Fora=b=c=e=1, α1=0.333, α2=1.234, in Table 2, we compare the absolute errors in the solution with those obtained using Chebyshev wavelet method [33]. It can be seen that our results are much more accurate. We also show absolute errors for different values of M and k in Table 3. It is observed that the absolute errors in solution are reduced by increasing the values of M and k.Table 2
Comparison of absolute error in solution of our method and reported in [33] for Example 3.
t
Our  method
Chebyshev  wavelets  method  [33]
M
=
9, k=3
m
=
384
0.1
4.318636
×
1
0
-
10
3.260279
×
1
0
-
7
0.2
6.606933
×
1
0
-
9
7.929112
×
1
0
-
7
0.3
4.151355
×
1
0
-
8
1.158235
×
1
0
-
6
0.4
9.317364
×
1
0
-
8
1.184173
×
1
0
-
6
0.5
1.561308
×
1
0
-
7
1.981458
×
1
0
-
6
0.6
2.914365
×
1
0
-
7
1.681103
×
1
0
-
6
0.7
4.226536
×
1
0
-
7
2.482836
×
1
0
-
6
0.8
5.374743
×
1
0
-
7
2.783712
×
1
0
-
6
0.9
3.388977
×
1
0
-
6
2.343684
×
1
0
-
6Table 3
Obtained absolute errors for different values ofM and k for Example 3.
t
M
=
9, k=3
M
=
12, k=3
M
=
15, k=3
M
=
20, k=4
0.1
4.318636
×
1
0
-
10
1.665675
×
1
0
-
10
7.328301
×
1
0
-
11
6.046732
×
1
0
-
11
0.2
6.606933
×
1
0
-
9
2.590659
×
1
0
-
9
1.134788
×
1
0
-
9
8.408221
×
1
0
-
10
0.3
4.151355
×
1
0
-
8
1.580848
×
1
0
-
8
7.193868
×
1
0
-
9
2.628054
×
1
0
-
9
0.4
9.317364
×
1
0
-
8
3.453152
×
1
0
-
8
1.573303
×
1
0
-
8
5.658458
×
1
0
-
9
0.5
1.561308
×
1
0
-
7
5.795390
×
1
0
-
8
2.634232
×
1
0
-
8
1.014053
×
1
0
-
8
0.6
2.914365
×
1
0
-
7
1.074635
×
1
0
-
7
4.947062
×
1
0
-
8
1.724892
×
1
0
-
8
0.7
4.226536
×
1
0
-
7
1.550796
×
1
0
-
7
7.144234
×
1
0
-
7
2.685408
×
1
0
-
8
0.8
5.374743
×
1
0
-
7
1.810074
×
1
0
-
7
7.570217
×
1
0
-
8
3.913804
×
1
0
-
8
0.9
3.388977
×
1
0
-
6
2.093040
×
1
0
-
6
1.384330
×
1
0
-
6
5.057264
×
1
0
-
8Example 4.
Consider the following boundary value problem in the case of the inhomogeneous Bagley-Torvik equation [53, 54]
(78)D2y(t)+D3/2y(t)+y(t)=t2+4tπ+2,y(0)=0,y(50)=2500,
where the exact solution is y(t)=t2.
It is worth mentioning that the method presented above only can be employed for solving this problem fort∈[0,1]. That is because the first kind Chebyshev wavelet is defined on interval [0,1]. However, the variable t in this example is defined on interval [0,50], and we should replace Chebyshev wavelet bases ψ(t) with ψ(t/50) in the discrete procedure. We applied the method described in Section 6 and got almost exact solution for t∈[0,50], although other authors had solved the problem for t∈[0,1] and t∈[0,5]. It can be seen from Figure 4 that the approximate solution and exact solution are closely overlapped for any t∈[0,50].Figure 4
Approximate and exact solutions for Example4 in the interval [0,50].Example 5.
Consider the following boundary value problem for nonlinear fractional order differential equation:(79)Dαy(t)+e-2π[y(t)]n=r(t),y(0)=a,y(1)=b,
where 1<α≤2, a,b∈ℝ, n∈ℕ, and r(t) is a given problem. For α=1.5, n=2, a=0, b=1, and r(t)=((15√π)/8Γ((7/2)-α))t(5/2)-α+e-2πt5, it can be easily verified that the exact solution is y(t)=t5/2. We use the introduced scheme in Section 6 to solve this example with M=20, k=5. Absolute error in solution for different values of α is presented in Table 4 which confirm the accuracy and efficiency of the proposed method. As it is expected, the absolute error is reduced as α approaches an integer value. Furthermore, we plot exact and numerical solutions for α=2 in Figure 5.Table 4
Obtained absolute errors forM=20, k=5 and different values of α for Example 5.
t
α
=
1.1
α
=
1.3
α
=
1.5
α
=
1.7
α
=
1.9
α
=
2
0.1
4.3764
×
1
0
-
8
1.4878
×
1
0
-
7
2.0748
×
1
0
-
7
1.7265
×
1
0
-
7
2.5621
×
1
0
-
6
5.3321
×
1
0
-
10
0.2
2.2918
×
1
0
-
8
2.5440
×
1
0
-
7
3.9721
×
1
0
-
7
2.2383
×
1
0
-
7
4.1274
×
1
0
-
6
4.7397
×
1
0
-
10
0.3
5.2510
×
1
0
-
9
3.3687
×
1
0
-
7
5.5097
×
1
0
-
7
2.3383
×
1
0
-
7
5.0983
×
1
0
-
6
4.1473
×
1
0
-
10
0.4
8.8262
×
1
0
-
9
3.9515
×
1
0
-
7
6.6732
×
1
0
-
7
2.0325
×
1
0
-
7
5.5447
×
1
0
-
6
3.5549
×
1
0
-
10
0.5
2.0156
×
1
0
-
8
4.2244
×
1
0
-
7
7.2700
×
1
0
-
7
1.6890
×
1
0
-
7
5.5744
×
1
0
-
6
2.9625
×
1
0
-
10
0.6
2.6409
×
1
0
-
8
4.1667
×
1
0
-
7
7.3063
×
1
0
-
7
1.1386
×
1
0
-
7
5.1809
×
1
0
-
6
2.3700
×
1
0
-
10
0.7
2.7986
×
1
0
-
8
3.7630
×
1
0
-
7
6.6892
×
1
0
-
7
6.2553
×
1
0
-
8
4.4105
×
1
0
-
6
1.7776
×
1
0
-
10
0.8
2.4855
×
1
0
-
8
2.9259
×
1
0
-
7
5.2546
×
1
0
-
7
2.8409
×
1
0
-
8
3.2886
×
1
0
-
6
1.1851
×
1
0
-
10
0.9
1.5051
×
1
0
-
8
1.6962
×
1
0
-
6
3.0942
×
1
0
-
7
4.8299
×
1
0
-
9
1.8027
×
1
0
-
6
5.9254
×
1
0
-
11Figure 5
Approximate and exact solutions forα=2 for Example 5.Example 6.
Consider the following boundary value problem:(80)Dαy(t)=Dβy(t)-g(t),y(0)=0,y(1)=0,
where 1<α≤2, 0<β≤1 and g(t)=et-1+1. The exact solution of the problem is not known generally. It can be easily verified that for α=2, β=1, the exact solution is y(t)=t(1-et-1). The problem (80) is solved numerically for integer order case in [55, 56] using Haar wavelet method and combined homotopy perturbation method, respectively. We solve the problem using M=12, k=3. In Table 5, the absolute errors are presented which confirm that the proposed method is more accurate. The numerical results for β=1 and different values of α plotted in Figure 6 show that as α tends to 2, the solution of fractional differential equation approaches to that of the integer-order differential equation.Table 5
Comparison of absolute errors between the proposed method withM=12, k=3, Haar wavelet method, and fourth order HPM, for different values of α for Example 6.
t
Current  method
Haar  wavelet  [55]
Fourth  order  HPM  [56]
0.1
9.0
×
1
0
-
20
3.4
×
1
0
-
8
5.2
×
1
0
-
6
0.2
2.1
×
1
0
-
19
2.7
×
1
0
-
8
2.7
×
1
0
-
8
0.3
2.1
×
1
0
-
19
2.9
×
1
0
-
8
1.1
×
1
0
-
9
0.4
2.0
×
1
0
-
20
3.6
×
1
0
-
8
7.9
×
1
0
-
6
0.5
9.0
×
1
0
-
20
4.0
×
1
0
-
8
3.6
×
1
0
-
6
0.6
5.0
×
1
0
-
19
5.2
×
1
0
-
8
1.4
×
1
0
-
6
0.7
1.1
×
1
0
-
18
6.6
×
1
0
-
8
5.3
×
1
0
-
6
0.8
1.2
×
1
0
-
18
7.8
×
1
0
-
8
6.5
×
1
0
-
6
0.9
6.2
×
1
0
-
19
9.4
×
1
0
-
8
4.5
×
1
0
-
6Figure 6
Chebyshev wavelet finite difference solutions forβ=1 and different values of α for Example 6.Example 7.
Consider the following linear fractional differential equation [26, 49, 50, 57]:
(81)Dαy(t)+y(t)=0,0<α≤2,
such that
(82)y(0)=1,y′(0)=0.
The conditiony′(0)=0 is only for 1<α≤2. The exact solution of (74) and (75) is given by
(83)y(t)=Ea(-tα),
where Ea(z)=∑k=0∞‍zk/Γ(αk+1) is the Mittag-Leffler function of order α. It can be easily verified that for α=1 and α=2, the exact solutions are y(t)=exp(-t) and y(t)=cos(t), respectively.
We solved the problem using the proposed method for different values ofα.
The absolute errors forα=0.2,0.4,0.6,0.8 (with M=6, k=6) and α=1.2,1.4,1.6,1.8 (with M=6, k=8) are shown in Table 6. It can be seen that our results are more accurate than the ones reported in [50]. The numerical solutions for α=0.5,0.75,0.95, and 1 (Figure 7(a)) and α=1.5,1.75,1.95, and 2 (Figure 7(b)) are plotted in Figure 7. It should be noted that as α approaches 1 and 2, the numerical solutions converge to the analytical solution y(t)=exp(-t) and y(t)=cos(t), respectively.Table 6
Comparison of absolute errors iny(t) in [50] and the ones obtained by our method for Example 7.
α
t
=
0.1
t
=
0.3
t
=
0.7
t
=
0.9
Our result
Ref. [50]
Our result
Ref. [50]
Our result
Ref. [50]
Our result
Ref. [50]
0.2
2.1
×
1
0
-
4
2.9
×
1
0
-
1
7.3
×
1
0
-
5
4.5
×
1
0
-
1
3.1
×
1
0
-
5
3.7
×
1
0
-
1
2.5
×
1
0
-
5
2.0
×
1
0
-
1
0.4
3.6
×
1
0
-
4
3.9
×
1
0
-
1
1.4
×
1
0
-
4
5.1
×
1
0
-
1
6.0
×
1
0
-
5
3.3
×
1
0
-
1
4.7
×
1
0
-
5
2.2
×
1
0
-
1
0.6
3.7
×
1
0
-
4
6.7
×
1
0
-
3
1.7
×
1
0
-
4
2.0
×
1
0
-
5
7.8
×
1
0
-
5
4.4
×
1
0
-
3
6.0
×
1
0
-
5
4.6
×
1
0
-
3
0.8
2.4
×
1
0
-
4
1.1
×
1
0
-
3
1.5
×
1
0
-
4
2.1
×
1
0
-
4
8.0
×
1
0
-
5
8.7
×
1
0
-
4
6.3
×
1
0
-
5
5.8
×
1
0
-
4
1.2
2.5
×
1
0
-
4
3.1
×
1
0
-
3
1.9
×
1
0
-
4
2.8
×
1
0
-
3
7.1
×
1
0
-
5
3.6
×
1
0
-
3
2.2
×
1
0
-
5
1.8
×
1
0
-
3
1.4
4.9
×
1
0
-
5
1.0
×
1
0
-
3
4.5
×
1
0
-
5
7.0
×
1
0
-
4
2.1
×
1
0
-
5
1.1
×
1
0
-
3
6.9
×
1
0
-
6
2.4
×
1
0
-
4
1.6
1.4
×
1
0
-
5
3.0
×
1
0
-
4
1.6
×
1
0
-
5
1.3
×
1
0
-
4
8.8
×
1
0
-
6
3.0
×
1
0
-
4
3.0
×
1
0
-
6
6.2
×
1
0
-
7
1.8
4.0
×
1
0
-
6
6.1
×
1
0
-
5
5.9
×
1
0
-
6
1.4
×
1
0
-
5
3.5
×
1
0
-
6
5.3
×
1
0
-
5
1.8
×
1
0
-
6
8.8
×
1
0
-
6Comparison ofy(t) with (a) α=0.5,0.75,0.95,1, (b) α=1.5,1.75,1.95,2, for Example 7.
(a)
(b)Example 8.
As the last example, consider the following linear multiterm fractional boundary value problem:(84)4(t+1)D2.5y(t)+4D1.5y(t)+1t+1y(t)=t+π
subject to
(85)y(0)=π,y′(0)=π2,y(1)=2π.
The exact solution of this problem isy(t)=π(t+1). This problem was solved in [58] by operational matrix of fractional derivatives using B-spline functions. We compare maximum absolute errors obtained by the introduced method with the ones reported in [58] in Table 7. It should be noted that the algebraic system of equations obtained by the method [58] is of order 2J+1, while the resulting one by the current method is of order 2k-1(M+1). It can be seen from Table 7 that our results are more accurate while we need to solve a system of algebraic equations of lower order.Table 7
Comparison of maximum absolute errors iny(t) in [58] and the ones obtained by current method for Example 8.
Presented method
Ref. [58]
M
=
6
M
=
8
M
=
10
M
=
12
J
Error
k
=
2
3.3
×
1
0
-
7
2.8
×
1
0
-
7
1.4
×
1
0
-
7
8.2
×
1
0
-
8
4
2.6
×
1
0
-
4
k
=
3
1.8
×
1
0
-
7
8.2
×
1
0
-
8
3.9
×
1
0
-
8
2.2
×
1
0
-
8
5
7.1
×
1
0
-
5
k
=
4
6.6
×
1
0
-
8
2.4
×
1
0
-
8
1.1
×
1
0
-
8
6.2
×
1
0
-
9
6
1.5
×
1
0
-
5
k
=
5
2.4
×
1
0
-
8
8.6
×
1
0
-
9
4.0
×
1
0
-
9
2.2
×
1
0
-
9
7
4.2
×
1
0
-
6


## 8. Conclusion

An efficient and accurate method based on hybrid of Chebyshev wavelets and finite difference methods was introduced. The fractional derivative is described in the Caputo sense. The useful properties of Chebyshev wavelets and finite difference method make it a computationally efficient method for solving the problems. The main advantage of the present method is the ability to represent smooth and especially piecewise smooth functions properly. Several examples are given to demonstrate the powerfulness of the proposed method. We note that the accuracy can be enhanced either by increasing the number of subintervals or by increasing the number of collocation points in subintervals properly. The validity and accuracy of the method were investigated for large intervals as well.



---

*Source: 916456-2013-12-12.xml*