# Wavelet Methods for Solving Fractional Order Differential Equations

Authors: A. K. Gupta; S. Saha Ray
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Fractional calculus is a field of applied mathematics which deals with derivatives and integrals of arbitrary orders. The fractional calculus has gained considerable importance during the past decades mainly due to its application in diverse fields of science and engineering such as viscoelasticity, diffusion of biological population, signal processing, electromagnetism, fluid mechanics, electrochemistry, and many more. In this paper, we review different wavelet methods for solving both linear and nonlinear fractional differential equations. Our goal is to analyze the selected wavelet methods and assess their accuracy and efficiency with regard to solving fractional differential equations. We discuss challenges faced by researchers in this field, and we emphasize the importance of interdisciplinary effort for advancing the study on various wavelets in order to solve differential equations of arbitrary order.


---

## Body

## 1. Introduction to Wavelets

Wavelet theory is a relatively new and emerging area in mathematical research. Wavelet methods have been used to develop accurate and fast algorithms for solving numerically integral and differential equations, especially those whose solutions are highly localized in position and scale. The concept of “wavelets” originated from the study of time-frequency signal analysis, wave propagation, and sampling theory. The main reason behind the discovery of wavelets and wavelet transforms is that the Fourier transform analysis does not contain the local information of signals. So the Fourier transform cannot be used for analyzing signals in joint time and frequency domain. In 1982, Jean Morlet, in collaboration with a group of French engineers, discovered the idea of wavelets transform for the analysis of nonstationary signals (signals containing transients and fractal structures).Wavelet method is an exciting method for solving difficult problems in mathematics, physics, and engineering, with modern applications in diverse fields such as wave propagation, data compression, image processing, pattern recognition, computer graphics, the detection of aircraft, and submarines and improvement in CAT scans and other medical technology. Morlet (1982) [1] first introduced the idea of wavelets as a family of functions constructed from translation and dilation of a single function called the “mother wavelet.”Waveletψ(x): an oscillatory function ψ(x)∈L2(R) with zero mean is a wavelet if it has the following desirable attributes:(1)
smoothness:ψ(x) is n times differentiable and its derivatives are continuous;(2)
localization:ψ(x) is well localized both in time and frequency domains; that is, ψ(x) and its derivatives must decay rapidly. For frequency localization Ψ^(ω) must decay sufficiently fast as ω→∞ and that Ψ^(ω) becomes flat in the neighborhood of ω=0. The flatness is associated with number of vanishing moments of ψ(x); that is,
(1)∫-∞∞xkΨ(x)dx=0orequivalentlydkdωkΨ^(ω)=0∫-∞∞xkΨ(x)dx=orequivlentlyfork=0,1,…,n
in the sense that the larger is the number of vanishing moments the more is the flatness when ω is small;(3)
the admissibility condition(2)∫-∞∞|Ψ^(ω)||ω|dω<∞
suggests that |Ψ^(ω)| decay at least as |ω|-1 or |x|ε-1 for ε>0.Although most of the numerical methods have been successfully applied for many linear and nonlinear differential equations, they have also some drawbacks in regions where singularities or sharp transitions occur. In those cases the solutions may be oscillating and for accurate representation of the results adaptive numerical schemes must be used which complicates the solution. To overcome the above difficulty wavelet transform methods are quite useful.The wavelet algorithms for solving differential equations usually are based on the collocation method. From the beginning of 1980s wavelets have been used invariably for the solution of differential equations. Most of the wavelet algorithms can handle easily periodic boundary conditions. In the present paper, the main idea is to apply wavelet methods, namely, Haar wavelet method, Legendre wavelet methods, Chebyshev wavelet method, B-spline wavelet, and so forth, for solving fractional differential equations (FDE).


### 1.1. Multiresolution Analysis (MRA) [2]

A set of subspaces{Vj}j∈Z is said to be MRA of L2(R) if it possess the following properties:
(3)Vj⊂Vj+1,∀j∈Z,(4)⋃j∈ZVjis dense inL2(R),(5)⋂j∈ZVj=ϕ,(6)f(x)∈Vj⟺f(2x)∈Vj+1,∀j∈Z,
where Z denotes the set of integers. Properties (3)–(5) state that {Vj}j∈Z is a nested sequence of subspaces that effectively covers L2(R). That is, every square integrable function can be approximated as closely as desired by a function that belongs to at least one of the subspaces Vj. A function φ∈L2(R) is called a scaling function if it generates the nested sequence of subspaces Vj and satisfies the dilation equation; namely,
(7)φ(x)=∑kpkφ(ax-k),
with pk∈l2 and a being any rational number.For each scalej, since Vj⊂Vj+1, there exists a unique orthogonal complementary subspace Wj of Vj in Vj+1. This subspace Wj is called wavelet subspace and is generated by ψj,k=ψ(2jx-k), where ψ∈L2 is called the wavelet. From the above discussion, these results follow easily:
(8)Vj1⋂Vj2=Vj2,j1>j2,Wj1⋂Wj2=0,j1≠j2,Vj1⋂Wj2=0,j1≤j2.


## 1.1. Multiresolution Analysis (MRA) [2]

A set of subspaces{Vj}j∈Z is said to be MRA of L2(R) if it possess the following properties:
(3)Vj⊂Vj+1,∀j∈Z,(4)⋃j∈ZVjis dense inL2(R),(5)⋂j∈ZVj=ϕ,(6)f(x)∈Vj⟺f(2x)∈Vj+1,∀j∈Z,
where Z denotes the set of integers. Properties (3)–(5) state that {Vj}j∈Z is a nested sequence of subspaces that effectively covers L2(R). That is, every square integrable function can be approximated as closely as desired by a function that belongs to at least one of the subspaces Vj. A function φ∈L2(R) is called a scaling function if it generates the nested sequence of subspaces Vj and satisfies the dilation equation; namely,
(7)φ(x)=∑kpkφ(ax-k),
with pk∈l2 and a being any rational number.For each scalej, since Vj⊂Vj+1, there exists a unique orthogonal complementary subspace Wj of Vj in Vj+1. This subspace Wj is called wavelet subspace and is generated by ψj,k=ψ(2jx-k), where ψ∈L2 is called the wavelet. From the above discussion, these results follow easily:
(8)Vj1⋂Vj2=Vj2,j1>j2,Wj1⋂Wj2=0,j1≠j2,Vj1⋂Wj2=0,j1≤j2.


## 2. Fractional Calculus

In the last few decades many authors pointed out that derivatives and integrals of noninteger order are very suitable for the description of properties of various real phenomena. Fractional derivatives [3] provide an excellent instrument for the description of memory and hereditary properties of various materials and processes. This is the main advantage of fractional derivatives in comparison with classical integer order models in which such effects are neglected.The mathematical modelling and simulation of systems and processes, based on the description of their properties in terms of fractional derivatives, naturally lead to differential equations of fractional order and to the necessity of solving such equations. However, effective general methods for solving them cannot be found even in the most useful works on fractional derivatives and integrals.Recently, orthogonal wavelets bases are becoming more popular for numerical solutions of differential equations due to their excellent properties such as ability to detect singularities, orthogonality, flexibility to represent a function at different levels of resolution, and compact support. In recent years, there has been a growing interest in developing wavelet based numerical algorithms for solution of fractional differential equations. Wavelets have been successfully applied for the solutions of ordinary and partial differential equations, integral equations, and integrodifferential equations of arbitrary order. Therefore, the main focus of the present paper is the application of different wavelet techniques for solving differential equations of arbitrary order.A computational approach to solve fractional differential equation is an essential work in scientific research. Some methods [4–11] for solving fractional differential equations are available in open literature. The B-spline wavelet collocation method [4] and variational iteration method (VIM) [5] have been applied to solve fractional differential equations. The learned researcher Saha Ray proposed some numerical methods for solving nonlinear fractional differential equations using generalized Haar wavelet method [6, 10] and Adomian decomposition method [7–9]. Haar wavelet method with operational matrices of integration [6] has been applied to solve fractional differential equations (FDE). Legendre wavelet method [11], Legendre multiwavelet Galerkin method [12], and also Chebyshev wavelets [13] can be applied to solve nonlinear fractional differential equations.


### 2.1. Fractional Derivative and Integration

The fractional calculus was first anticipated by Leibnitz, who was one of the founders of standard calculus, in a letter written in 1695. This fractional calculus involves different definitions of the fractional operators such as the Riemann-Liouville fractional derivative, Caputo derivative, Riesz derivative, and Grunwald-Letnikov fractional derivative [3]. Riemann-Liouville fractional derivative is not suitable for real world physical problems since it requires the definition of fractional order initial conditions, which have no physically meaningful explanation yet. Caputo introduced an alternative definition, which has the advantage of defining integer order initial conditions for fractional order differential equations.


#### 2.1.1. Definition of Riemann-Liouville Integral

The most frequently encountered definition of an integral of fractional order is the Riemann-Liouville integral [3], in which the fractional integral operator Jα of a function f(t) is defined as [3]:
(9)Jαf(t)=1Γ(α)∫0t(t-τ)α-1f(τ)dτ,α>0,α∈R+,
where Γ(·) is the well-known gamma function, and some properties of the operator Jα are as follows:
(10)JαJβf(t)=Jα+βf(t),(α>0,β>0),Jαtγ=Γ(1+γ)Γ(1+γ+α)tα+γ,(γ>-1).


#### 2.1.2. Definition of Caputo Fractional Derivative

The fractional derivative introduced by Caputo [3], in the late sixties, is called Caputo fractional derivative. The Caputo fractional derivative Dtα0 of a function f(t) is defined as [3]
(11)Dtα0f(t)=1Γ(n-α)∫0tfn(τ)(t-τ)α-n+1dτ,hhlhh(n-1<α≤n,n∈N).
The following are two basic properties of the Caputo fractional derivative:
(12)Dtα0tβ=Γ(1+β)Γ(1+β-α)tβ-α,0<α<β+1,β>-1,JαDαf(t)=f(t)-∑k=0n-1fk(0+)tkk!,hln-1<α≤n,n∈N.
For Caputo’s derivative we have
(13)DαC=0,(Cisaconstant).
Similar to integer order differentiation, Caputo’s derivative is linear:
(14)Dα(γf(t)+δg(t))=γDαf(t)+δDαg(t),
where γ and δ are constants, and it satisfies the so-called Leibnitz’s rule Dα(g(t)f(t))=∑k=0∞(αk)g(k)(t)Dα-kf(t), where f(τ) is continuous in [0,t] and g(τ) has continuous derivatives sufficient number of times in [0,t].


#### 2.1.3. Definition of Grunwald-Letnikov Fractional Derivative

The Grunwald-Letnikov fractional derivative is defined as [3, 10]
(15)Dtαaf(t)=lim⁡h→0nh=t-ah-p∑r=0n(-1)r(pr)f(t-rh),
where ωrp=(-1)r(pr) and
(16)ω0p=1,ωrp=(1-p+1r)ωr-1p,r=1,2,….


## 2.1. Fractional Derivative and Integration

The fractional calculus was first anticipated by Leibnitz, who was one of the founders of standard calculus, in a letter written in 1695. This fractional calculus involves different definitions of the fractional operators such as the Riemann-Liouville fractional derivative, Caputo derivative, Riesz derivative, and Grunwald-Letnikov fractional derivative [3]. Riemann-Liouville fractional derivative is not suitable for real world physical problems since it requires the definition of fractional order initial conditions, which have no physically meaningful explanation yet. Caputo introduced an alternative definition, which has the advantage of defining integer order initial conditions for fractional order differential equations.


### 2.1.1. Definition of Riemann-Liouville Integral

The most frequently encountered definition of an integral of fractional order is the Riemann-Liouville integral [3], in which the fractional integral operator Jα of a function f(t) is defined as [3]:
(9)Jαf(t)=1Γ(α)∫0t(t-τ)α-1f(τ)dτ,α>0,α∈R+,
where Γ(·) is the well-known gamma function, and some properties of the operator Jα are as follows:
(10)JαJβf(t)=Jα+βf(t),(α>0,β>0),Jαtγ=Γ(1+γ)Γ(1+γ+α)tα+γ,(γ>-1).


### 2.1.2. Definition of Caputo Fractional Derivative

The fractional derivative introduced by Caputo [3], in the late sixties, is called Caputo fractional derivative. The Caputo fractional derivative Dtα0 of a function f(t) is defined as [3]
(11)Dtα0f(t)=1Γ(n-α)∫0tfn(τ)(t-τ)α-n+1dτ,hhlhh(n-1<α≤n,n∈N).
The following are two basic properties of the Caputo fractional derivative:
(12)Dtα0tβ=Γ(1+β)Γ(1+β-α)tβ-α,0<α<β+1,β>-1,JαDαf(t)=f(t)-∑k=0n-1fk(0+)tkk!,hln-1<α≤n,n∈N.
For Caputo’s derivative we have
(13)DαC=0,(Cisaconstant).
Similar to integer order differentiation, Caputo’s derivative is linear:
(14)Dα(γf(t)+δg(t))=γDαf(t)+δDαg(t),
where γ and δ are constants, and it satisfies the so-called Leibnitz’s rule Dα(g(t)f(t))=∑k=0∞(αk)g(k)(t)Dα-kf(t), where f(τ) is continuous in [0,t] and g(τ) has continuous derivatives sufficient number of times in [0,t].


### 2.1.3. Definition of Grunwald-Letnikov Fractional Derivative

The Grunwald-Letnikov fractional derivative is defined as [3, 10]
(15)Dtαaf(t)=lim⁡h→0nh=t-ah-p∑r=0n(-1)r(pr)f(t-rh),
where ωrp=(-1)r(pr) and
(16)ω0p=1,ωrp=(1-p+1r)ωr-1p,r=1,2,….


## 2.1.1. Definition of Riemann-Liouville Integral

The most frequently encountered definition of an integral of fractional order is the Riemann-Liouville integral [3], in which the fractional integral operator Jα of a function f(t) is defined as [3]:
(9)Jαf(t)=1Γ(α)∫0t(t-τ)α-1f(τ)dτ,α>0,α∈R+,
where Γ(·) is the well-known gamma function, and some properties of the operator Jα are as follows:
(10)JαJβf(t)=Jα+βf(t),(α>0,β>0),Jαtγ=Γ(1+γ)Γ(1+γ+α)tα+γ,(γ>-1).


## 2.1.2. Definition of Caputo Fractional Derivative

The fractional derivative introduced by Caputo [3], in the late sixties, is called Caputo fractional derivative. The Caputo fractional derivative Dtα0 of a function f(t) is defined as [3]
(11)Dtα0f(t)=1Γ(n-α)∫0tfn(τ)(t-τ)α-n+1dτ,hhlhh(n-1<α≤n,n∈N).
The following are two basic properties of the Caputo fractional derivative:
(12)Dtα0tβ=Γ(1+β)Γ(1+β-α)tβ-α,0<α<β+1,β>-1,JαDαf(t)=f(t)-∑k=0n-1fk(0+)tkk!,hln-1<α≤n,n∈N.
For Caputo’s derivative we have
(13)DαC=0,(Cisaconstant).
Similar to integer order differentiation, Caputo’s derivative is linear:
(14)Dα(γf(t)+δg(t))=γDαf(t)+δDαg(t),
where γ and δ are constants, and it satisfies the so-called Leibnitz’s rule Dα(g(t)f(t))=∑k=0∞(αk)g(k)(t)Dα-kf(t), where f(τ) is continuous in [0,t] and g(τ) has continuous derivatives sufficient number of times in [0,t].


## 2.1.3. Definition of Grunwald-Letnikov Fractional Derivative

The Grunwald-Letnikov fractional derivative is defined as [3, 10]
(15)Dtαaf(t)=lim⁡h→0nh=t-ah-p∑r=0n(-1)r(pr)f(t-rh),
where ωrp=(-1)r(pr) and
(16)ω0p=1,ωrp=(1-p+1r)ωr-1p,r=1,2,….


## 3. Wavelet Methods for Fractional Differential Equations

Wavelet analysis is a numerical concept which allows representing a function in terms of a set of basis functions, called wavelets, which are localized both in location and scale. Wavelets used in this method are mostly compact support introduced by Daubechies [14]. The wavelet based approximations of ordinary and partial differential equations have been attracting the attention, since the contribution of orthonormal bases of compactly supported wavelet by Daubechies and multiresolution analysis based Fast Wavelet transform algorithm by Beylkin et al. [15] gained momentum to make wavelet approximations attractive.


### 3.1. Haar Wavelets

Haar wavelet functions have been used from 1910 and were introduced by the Hungarian mathematician Alfred Haar. Haar wavelets (which are Daubechies wavelets of order 1) consist of piecewise constant functions on the real line that can take only three values, that is, 0, 1, and −1, and are therefore the simplest orthonormal wavelets with a compact support. Haar wavelet method is to be used due to the following features: being simpler and fast, flexible, and convenient, having small computational costs, and being computationally attractive. The Haar functions are a family of switched rectangular wave forms where amplitudes can differ from one function to another. These properties of Haar wavelets are utilized to reduce the computation of integral equations to some algebraic equations.Usually the Haar wavelets are defined for the intervalt∈[0,1) but in general case t∈[A,B], we divide the interval [A,B] into m equal subintervals, each of width Δt=(B-A)/m. In this case, the orthogonal set of Haar functions are defined in the interval [A,B] by [6]
(17)h0(t)={1t∈[A,B],0elsewhere,hi(t)={1,ζ1(i)≤t<ζ2(i)-1,ζ2(i)≤t<ζ3(i)0,otherwise,
where
(18)ζ1(i)=A+(k-12j)(B-A)=A+(k-12j)mΔt,ζ2(i)=A+(k-(1/2)2j)(B-A)=A+(k-(1/2)2j)mΔt,ζ3(i)=A+(k2j)(B-A)=A+(k2j)mΔt,
for i=1,2,…,m, m=2J and J is a positive integer which is called the maximum level of resolution. Here j and k represent the integer decomposition of the index i; that is, i=k+2j-1, 0≤j<i and 1≤k<2j+1.


#### 3.1.1. Operational Matrix for General Order Integration [6]

The integration of theHm(t)=[h0(t),h1(t),…,hm-1(t)]T can be approximated by [16]
(19)∫0tHm(τ)dτ≅QHm(t),
where Q is called the Haar wavelet operational matrix of integration which is a square matrix of m-dimension. To derive the Haar wavelet operational matrix of the general order of integration, we recall the fractional integral of order α(>0) which is defined by Podlubny [3]:
(20)Jαf(t)=1Γ(α)∫0t(t-τ)α-1f(τ)dτ,t>0,α∈R+,
where R+ is the set of positive real numbers.The Haar wavelet operational matrixQα for integration of the general order α is given by
(21)QαHm(t)=JαHm(t)=[Jαh0(t),Jαh1(t),…,Jαhm-1(t)]T=[Qh0(t),Qh1(t),…,Qhm-1(t)]T,
where
(22)Qh0(t)={tαΓ(1+α),t∈[A,B],0,elsewhere,Qhi(t)={0,A≤t<ζ1(i),ϕ1,ζ1(i)≤t<ζ2(i),ϕ2,ζ2(i)≤t<ζ3(i),ϕ3,ζ3(i)≤t<B,
where
(23)ϕ1=(t-ζ1(i))αΓ(α+1),ϕ2=(t-ζ1(i))αΓ(α+1)-2(t-ζ2(i))αΓ(α+1),ϕ3=(t-ζ1(i))αΓ(α+1)-2(t-ζ2(i))αΓ(α+1)+(t-ζ3(i))αΓ(α+1),
for i=1,2,…,m, m=2J and J is a positive integer, called the maximum level of resolution. Here j and k represent the integer decomposition of the index i. That is, i=k+2j-1, 0≤j<i and 1≤k<2j+1.


#### 3.1.2. Function Approximation

Any functiony(t)∈L2([0,1)) can be expanded into Haar wavelets by [6, 10, 16]
(24)y(t)=c0h0(t)+c1h1(t)+c2h2(t)+⋯,hhlhhhhwherecj=∫01y(t)hj(t)dt.
If y(t) is approximated as piecewise constant in each subinterval, the sum in (24) may be terminated after m terms and consequently we can write discrete version in the matrix form as
(25)Y≈∑i=0m-1cihi(tl)=CmTHm,
where Y and CmT are the m-dimensional row vectors.HereH is the Haar wavelet matrix of order m defined by H=[h0,h1,…,hm-1]T; that is,
(26)H=[h0h1⋮hm-1]=[h0,0h0,1⋯h0,m-1h1,0h1,1⋯h1,m-1⋮hm-1,0hm-1,1⋯hm-1,m-1],
where h0,h1,…,hm-1 are the discrete form of the Haar wavelet bases.The collocation points are given by(27)tl=A+(l-0.5)Δt,l=1,2,…,m.
Similarly, a function of two variables y(x,t)∈L2([0,1]×[0,1]) can be approximated by discrete version in the matrix form of Haar wavelets as
(28)Y(x,t)=HT(x)CH(t),
where Y(x,t) is the discrete form of y(x,t) and C is the coefficient matrix of Y, and it can be obtained by the following formula:
(29)C=(HT)-1·Y·H-1.
Since H is orthogonal, (29) becomes
(30)C=H·Y·H-1.


#### 3.1.3. Application of Haar Wavelet on FDE [17]

In order to show the effectiveness of Haar wavelet method for solving fractional differential equations, we consider the following numerical example of variable coefficient fractional convection diffusion equation:(31)∂αu∂tα=-b(x)∂u∂x+a(x)∂βu∂xβ+q(x,t),
where 0≤x<1, 0<t≤1, a(x)>0, b(x)>0, 0<α≤1, 1<β≤2, q(x,t)∈C(D), and D=[0,1]×[0,1].Sinceu(x,t)∈L2(R), we suppose
(32)u(x,t)≈∑i=0m-1‍∑j=0m-1cijhi(x)hj(t).
Then we can obtain the discrete form of (32) by taking step Δ=1/m of x and t:
(33)U=HT(x)CH(t)
and then
(34)∂αu∂tα≈∂αU∂tα∂αu∂tα=HT(x)C∂αH(t)∂tα=HT(x)CQ-αH(t),∂βu∂xβ≈∂βU∂xβ∂αu∂tα=(∂β∂xβHT(x))CH(t)∂αu∂tα=(∂β∂xβH(x))TCH(t)=(Q-βH(x))TCH(t).
Here, q(x,t) is known function; discretizing it, we have
(35)D=(q(xi,tj)),i,j=0,1,2,…,m-1.
Coefficients a(x),b(x) can be dispersed into a(xi),b(xi), i=0,1,2,…,m-1.Let(36)A=[a(x0)0⋯00a(x1)⋯0⋮⋱⋮0⋯0a(xm-1)],B=[b(x0)0⋯00b(x1)⋯0⋮⋱⋮0⋯0b(xm-1)].
Substituting (33)–(36) in (31), we have
(37)HT(x)CQ-αH(t)=-B(∂∂xH(x))TCH(t)+A(∂β∂xβH(x))TCH(t)+D.
Equation (37) is a system of algebraic equations in a(xi) and b(xi). By solving this system of equations using mathematical software, the Haar wavelet coefficients a(xi) and b(xi) can be obtained.


### 3.2. Cubic B-Spline Basis Functions onH2(I)

LetI=[0,L] be an interval with L>4, and let H2(I) be a Sobolev space which contains functions with square integrable second derivatives, and the homogenous Sobolev space H02(I) can be defined by
(38)H02(I)={f′(0)f′(L)f(t)∈H2(I),f(0)=f′(0)=f(L)=f′(L)=0f′(0)f′(L)},
which is a Hilbert space equipped with inner product:
(39)〈f,g〉=∫If′′(t)g′′(t)dt.
Cai and Wang [18] present a multiresolution analysis (MRA) and a wavelet decomposition for H02(I) by constructing scaling spline functions
(40)φ(t)=16∑l=04(4l)(-1)l(t-l)+3,φb(t)=32t+2-1112t+3+32(t-1)+3-34(t-2)+3+16(t-3)+3,
and wavelet functions
(41)ψ(t)=-37φ(2t)+127φ(2t-1)-37φ(2t-2),ψb(t)=2413φb(2t)-613φ(2t),
which have compact supports, where
(42)t+n={tn,t>00,t≤0.
Then the scaling spaces V0⊂V1⊂V2⋯⊂V∞=H02(I) of MRA and a wavelet decomposition H02(I)=V0⊕W0⊕W1⊕⋯⊕Wj-1⊕⋯ corresponding to the MRA are obtained by defining
(43)Vj=span⁡{φj,k(t),φbj(t),φbj(L-t),0≤k≤nj-4},Wj=span⁡{ψj,k(t),k=-1,0,…,nj-2},
in which nj=2jL,
(44)φj,k(t)=φ(2jt-k),φbj(t)=φb(2jt),hhhhhhhhhhhhhhhhhhhlhh0≤k≤nj-4,ψj,k(t)=ψ(2jt-k),j≥0,k=0,1,…,nj-3,ψj,-1(t)=ψb(2jt),ψj,nj-2(t)=ψb(2j(L-t)).
For convenience, we set ψ-1,k(t)=φ0,k(t), 0≤k≤L-4, ψ-1,-1(t)=φb(t), ψ-1,L-3(t)=φb(L-t), and n-1=L-1. Let Bj={2-3j/2ψj,k(t)}k=-1nj-2, -1≤j<∞. Wang [19] derived that B=⋃j=-1∞Bj is an unconditional basis of H02(I), which turns out to be a basis of continuous space C0(I). For nonhomogeneous Sobolev space H2(I), Cai and Wang [18] introduced boundary spline functions
(45)η1(t)=(1-t)+3,η2(t)=2t+-3t+2+76t+3-43(t-1)+3+16(t-2)+3
to deal with the values of functions at boundary points.


#### 3.2.1. Function Approximation

Any functionf(t)∈H02(I) can be uniquely expanded into cubic spline wavelet series by
(46)f(x)=∑j=-1∞∑k=1njdj,kψj,k-2(t),
with dj,k=∫If′′(t)(ψj,k-2*)(t)dt, where ψj,k*(t) are dual functions of ψj,k. Truncating the infinite series in (46) at J-1, we get
(47)fJ(t)=∑j=-1J-1∑k=1njdj,kψj,k-2(t).
From Wang [19] we have
(48)∥f(t)-fJ(t)∥H022⟶0asJ⟶+∞.
Hence any function f(t)∈H02(I) can be approximated by fJ(t) defined in (47). And any function f(t)∈H2(I) can be approximated by
(49)fJ(t)=Ib,Jf(t)+∑j=-1J-1∑k=1njdj,kψj,k-2(t),
and the approximation order is O(2-4J) if f(t) is sufficiently smooth, where
(50)Ib,Jf(t)=a1η1(2Jt)+a2η2(2Jt)+a3η2(2J(L-t))+a4η1(2J(L-t)).
Suppose N=2JL+3, and ΩJ(t) is a 1×N vector given by
(51)ΩJ(t)=[ψJ-1,nJ-2(t)η1(2Jt),η2(2Jt),η2(2J(L-t)),hhhη1(2J(L-t)),ψ-1,-1(t),ψ-1,0(t),…,hhhψ-1,n-1-2(t),ψ0,-1(t),ψ0,0(t),ψ0,1(t),…,ψ0,L-3(t),hhhψ0,n0-2(t),ψ1,-1(t),ψ1,0(t),ψ1,1(t),…,hhhψ1,2L-3(t),ψ1,n1-2(t),…,ψJ-1,-1(t),hhhψJ-1,0(t),ψJ-1,1(t),…,ψJ-1,nJ-3(t),ψJ-1,nJ-2(t)]≡[ω1(t),ω2(t),…,ωN(t)].fJ(t) defined in (49) can be rewritten as
(52)fJ(t)=∑k=1Nf^kωk(t)=ΩJ(t)f^,
where f^=[f^1,f^2,…,f^N]Tare the wavelet expansion coefficients which can be determined by interpolating at collocation points
(53)t-1(-1)=12J+1,tL+1(-1)=L-12J+1,tk(-1)=k,k=0,1,2,…,L,tk(j)=k+1.52j,-1≤k≤2jL-2,0≤j≤J-1.


#### 3.2.2. Application of Cubic Spline Wavelets on FDE [4]

Consider the fractional differential equation of the form(54)Dαy(t)=-y(t)+f(t),0<t≤L,
with
(55)y(n)(0)=y0(n),n=0,1,2,…,m,hhhhhlm<α≤m+1,m∈N.
To solve problem (54), we approximate y(t) by
(56)y(t)≈∑k=1Nckωk(t)=Yj(t),
where ck=[c1,c2,…,cN]T is unknown. The α order fractional derivative of y(t) is approximated by
(57)Dαy(t)≈∑k=1NckDαωk(t)≡DαYj(t).
At collocation points,
(58)Yj(tk)=Y(tk),k=1,2,…,N,DαYj(tk)=DαY(tk),k=1,2,…,N.
Interpolating the fractional differential equation (54) by Y(tk) at all collocation points, we obtain
(59)DαYj(tk)=-Yj(tk)+f(tk),2≤k≤N,yj(n)(0)=y(n)(0),n=0,1,2,…,m.
Letting
(60)B1=[ω1(t2)ω2(t2)⋯ωN(t2)⋮⋱⋮ω1(tN)ω2(tN)⋯ωN(tN)],B2=[Dαω1(t2)Dαω2(t2)⋯DαωN(t2)⋮⋱⋮Dαω1(tN)Dαω2(tN)⋯DαωN(tN)],
where B1 and B2 are obtained analytically.Now (59) can be represented by
(61)Ay^=b,cny^=y0(n),
where A=B1+B2 and
(62)b=(f(t2),…,f(tN))T,cn=(ω1(n)(0),ω2(n)(0),…,ωN(n)(0)),n=0,1,2,…,m.
Let AmT=(AT,c0T,c1T,…,cmT), bmT=(bT,y0(0),y0(1),…,y0(m)); then (61) can be written as
(63)Amy^=bm.
Consequently the wavelet coefficient y^ can be obtained by solving linear equations (63).


### 3.3. Chebyshev Wavelets

The Legendre wavelet and Chebyshev wavelet are constructed from their corresponding polynomials. These wavelets are useful tools in the numerical computations. The second kind Chebyshev polynomials have many good properties and are widely applied in different disciplines.Wavelets constitute a family of functions constructed from dilation and translation of a single functionψ(t) called the mother wavelet. When the dilation parameter a and the translation parameter b vary continuously, we have the following family of continuous wavelets as [13]
(64)ψa,b(t)=|a|-1/2ψ(t-ba),a,b∈R,a≠0.
If we restrict the parameters a and b to discrete values as a=a0-k, b=nb0a0-k; a0>1, b0>0, where n and k are positive integer, the following family of discrete wavelets are defined as
(65)ψk,n(t)=|a0|k/2ψ(a0kt-nb0),
where ψk,n(t) forms a wavelet basis for L2(R). In particular, when a0=2 and b0=1, then ψk,n(t) forms an orthonormal basis; that is,
(66)〈ψj,k(t),ψl,m(t)〉=δjlδkm.


#### 3.3.1. Chebyshev Wavelet and Function Approximation

Chebyshev waveletsψn,m(t), on the interval [0,1), are defined as [20]
(67)ψn,m(t)={2k/2T~m(2kt-2n+1),n-12k-1≤t<n2k-10,otherwise,
where
(68)T~m={1π,m=0,2πTm(t),m>0,
and m=0,1,…,M-1, n=1,2,…,2k-1,k is any positive integer, and Tm(t) are Chebyshev polynomials of the first kind of degree m which are orthogonal with respect to the weight function ω(t)=(1/1-t2), on the interval [-1,1], and Tm(t) can be determined by the following recurrence formula:
(69)T0(t)=1,T1(t)=t,sTm+1(t)=2tTm(t)-Tm-1(t),m=1,2,3,….
A function f(t) defined over [0, 1) may be expanded as
(70)f(t)=∑n=1∞‍∑m=0∞cn,mψn,m(t),
where cnm=〈f(t),ψn,m(t)〉.If the infinite series in (70) is truncated, then (70) can be written as
(71)f(t)≈∑n=12k-1‍∑m=0M-1cn,mψn,m(t)=CTΨ(t),
where C and Ψ(t) are 2k-1M×1 matrices given by
(72)C≡[c1,0,c1,1,…,c1,M-1,c2,0,…,c2,M-1,…,c2k-1,0,…,c2k-1,M-1]T,(73)Ψ(t)≡[ψ1,0,ψ1,1,…,ψ1,M-1,ψ2,0,…,ψ2,M-1,…,ψ2k-1,0,…,ψ2k-1,M-1]T.
Taking the collocation points as follows:
(74)ti=2i-12kM,i=1,2,…,2k-1M,
we define the Chebyshev wavelet matrix Φm×m as
(75)Φm×m≡[Ψ(12m)Ψ(32m)⋯Ψ(2m-12m)].


#### 3.3.2. Operational Matrix of the Fractional Integration

The integration of the vectorΨ(t) defined in (73) can be obtained as
(76)∫0tΨ(τ)dτ≈PΨ(t),
where P is the 2k-1M×2k-1M operational matrix for integration [13].We know the Riemann-Liouville fractional integral operatorJα, of a function f(t), is defined as [3]
(77)Jαf(t)=1Γ(α)∫0t(t-τ)α-1f(τ)dτ=1Γ(α)tα-1*f(t),
where α∈R+(α>0) is the order of integration, Γ(·) is the well-known gamma function, and tα-1*f(t) denotes the convolution product of tα-1 and f(t). Now if f(t) is expanded in Chebyshev wavelets, the Riemann-Liouville fractional integration becomes
(78)Jαf(t)=1Γ(α)tα-1*f(t)≅CT1Γ(α){tα-1*Ψ(t)}.
Thus if tα-1*f(t) can be integrated and then expanded in Chebyshev wavelets, the Riemann-Liouville fractional integration can be solved via Chebyshev wavelets.Also we define am-set of Block-Pulse functions (BPF) as [13]
(79)bi(t)={1,im≤t<i+1m0,otherwise,
where i=0,1,2,…,(m-1).The functionsbi(t) are disjoint and orthogonal; that is,
(80)bi(t)bl(t)={0,i≠lbi(t),i=l,∫01bi(τ)bl(τ)dτ={0,i≠l1m,i=l.
Similarly, Chebyshev wavelets may be expanded into an m-term block pulse functions (BPF) as
(81)Ψm′(t)=Φm′×m′Bm′(t),
where m′=2k-1M and Bm′(t)≡[b0(t)b1(t)⋯bi(t)⋯bm-1(t)]T.In [21], Kilicman and Al Zhour have proposed the block pulse operational matrix of the fractional integration Fα as follows:
(82)Jα(Bm′(t))≈FαBm′(t),
where
(83)Fα=1mα1Γ(α+2)[1ξ1ξ2···ξm-101ξ1···ξm-2001···ξm-3······000··00001]
with ξk=(k+1)α+1-2kα+1+(k-1)α+1.Next we derive the Chebyshev wavelet operational matrix of the fractional integration. Let(84)Jα(Ψm(t))≈Pm′×m′αΨm(t),
where matrix Pm′×m′α is called the Chebyshev wavelet operational matrix of the fractional integration.Using (81) and (82), we have
(85)JαΨm′(t)≈(JαΦm′×m′Bm′)(t)=Φm′×m′(JαBm′)(t)≈Φm′×m′FαBm′(t).
From (84) and (85) we get
(86)Pm′×m′αΨm′(t)=Φm′×m′FαBm′(t).
Then, the Chebyshev wavelet operational matrix of the fractional integration Pm′×m′α is given by
(87)Pm′×m′α=Φm′×m′FαΦm′×m′-1.
Thus, we derive the operational matrix of fractional integration for the second kind Chebyshev wavelet.


#### 3.3.3. Application of Chebyshev Wavelet on FDE [22]

In this section, we use Chebyshev wavelet operational matrices of the fractional integration to solve nonlinear fractional differential equation. We consider fractional Riccati equation(88)Dαy(t)+[y(t)]2=1,0<α≤1,
subject to initial state y(0)=0.Let(89)Dαy(t)=CTΨ(t),
and then
(90)y(t)=CTPm′×m′αΨ(t).
Since Ψ(t)=Φm′×m′Bm′(t), we have
(91)y(t)=CTPm′×m′αΦm′×m′Bm′(t).
Let CTPm′×m′αΦm′×m′=[a1,a2,…,am′], and, using (79), we get
(92)[y(t)]2=[a1b1(t)+a2b2(t)+⋯+am′bm′(t)]2=[a12,a22,…,am′2]Bm′(t).
Substituting (89) and (92) in (88), we have
(93)CTΦm′×m′Bm′(t)+[a12,a22,…,am′2]Bm′(t)-[1,1,…,1]Bm′(t)=0.
Solving the above nonlinear system of algebraic equations, we can find the vector C and consequently the solution for y(t).


### 3.4. Legendre Wavelets

Legendre waveletψn,m=ψ(k,n,m,t) have four arguments: n=1,…,2k-1,k can be assumed as any positive integer, m is order of Legendre polynomials, and t denotes the time.The Legendre polynomials of orderm, denoted by Lm(t), are defined on the interval [−1, 1] and can be determined with the help of the following recurrence formulae [11]:
(94)L0(t)=1,L1(t)=t,(95)Lm+1(t)=2m+1m+1tLm(t)-mm+1Lm-1(t),hhhhhhhhhhhhhhhhhhhhhhm=1,2,3,….
The Legendre wavelets are defined on interval [0, 1) by
(96)ψn,m(t)={(m+12)1/22k/2Lm(2kt-n^),n^-12k≤t<n^+12k,0,elsewhere
for k=2,3,…n^=2n-1, n=1,2,3,…,2k-1, m=0,1,…,M-1 is the order of the Legendre polynomials, and M is a fixed positive integer. The set of Legendre wavelets forms an orthogonal basis of L2(R).


#### 3.4.1. Function Approximation

A functionf(t)∈L2[0,1] can be expanded in terms of Legendre wavelet as [23]
(97)f(t)=∑n=1∞‍∑m=0∞cn,mψn,m(t),
where cn,m=〈f(t),ψn,m(t)〉.If the infinite series in (97) is truncated, then (97) can be written as
(98)f(t)≈∑n=12k-1‍∑m=0M-1cn,mψn,m(t)=CTΨ(t),
where C and Ψ(x) are 2k-1M×1 matrices given by
(99)C≡[c1,0,c1,1,…,c1,M-1,c2,0,…,c2,M-1,…,c2k-1,0,…,c2k-1,M-1]T,Ψ(t)≡[ψ1,0,ψ1,1,…,ψ1,M-1,ψ2,0,…,ψ2,M-1,…,ψ2k-1,0,…,ψ2k-1,M-1]T.
Similarly, a function of two variables k(x,t)∈L2([0,1]×[0,1]) can be approximated by Legendre wavelet series as
(100)k(x,t)≈ΨT(t)KΨ(x),
where K=[kij] is (2k-1M×2k-1M) matrix, with
(101)ki,j=〈ψi(t),〈k(x,t),ψj(x)〉〉.


#### 3.4.2. Application of Legendre Wavelet on FDE

In order to show the effectiveness of Legendre wavelet method for solving fractional differential equations, we consider the following numerical example [11]:
(102)Dαy(t)+y(t)=0,0<α≤2,
such that
(103)y(0)=1,y′(0)=0.
The corresponding integral representation for (102) and (103) is
(104)y=-Jαy(t)+y(0)+y′(0)t.
Let y be approximated using Legendre wavelet as
(105)y=CTΨ(t),
and then
(106)Jαy=CTJαΨ(t)=CTPm^×m^αΨ(t),hhhhhhhlhhhhwherem^=2k-1M.
Substituting (105) and (106) in (104), we have the following system of algebraic equations:
(107)CTΨ(t)+CTPm^×m^αΨ(t)=Y0,whereY0=y(0)+y′(0)t.
By solving the linear system (107), we can find the vector C and hence consequently the solution y(t).


### 3.5. Legendre Multiwavelets Galerkin Method

In this section, the application of Legendre multiwavelet Galerkin method for providing approximate solutions for initial value problems of fractional differential equations has been discussed. Legendre multiwaveletψn,m(t)=ψ(k,n,m,t) have four arguments: n=0,1,…,2k-1, k can be assumed as any positive integer, m is order of Legendre polynomials, and t denotes the normalized time. The Legendre multiwavelets are defined on the interval [0, 1) by
(108)ψn,m(t)={(2m+1)1/22k/2Pm(2kt-n),n2k≤t<n+12k,0,elsewhere,
where m=0,1,…,M-1 and n=1,2,3,…,2k-1. The coefficient (2m+1)1/2 is for orthonormality and Pm(t) are the well-known shifted Legendre polynomials of order m which are defined on the interval [0,1] and can be determined with the help of the following recurrence formula:
(109)P0(t)=1,P1(t)=2t-1,(110)Pm+1(t)=2m+1m+1(2t-1)Pm(t)-mm+1Pm-1(t),hhhhhhhhhhhhhhhhhlhhhm=1,2,3,….
The two-dimensional Legendre multiwavelets are defined as
(111)ψn1,m1,n2,m2(x,t)={APm1(2k1x-n1)Pm2(2k2t-n2),n12k1≤x<n1+12k1,n22k2≤t<n2+12k20,elsewhere,
where
(112)A=(2m1+1)(2m2+1)2(k1+k2)/2.n1 and n2 are defined similarly to n, k1 and k2 are any positive integer, m1 and m2 are the order for Legendre polynomials, and ψn1,m1,n2,m2(x,t) forms a basis for L2([0,1]×[0,1]).


#### 3.5.1. Function Approximation

A functionf(x,t) defined over [0,1]×[0,1] can be expanded in terms of Legendre multiwavelet as [12]
(113)f(x,t)=∑n=0∞∑i=0∞‍∑l=0∞‍∑j=0∞cn,i,l,jψn,i(x)ψl,j(t).
If the infinite series in (113) is truncated, then (113) can be written as
(114)f(x,t)=∑n=02k1-1‍∑i=0N∑l=02k2-1∑j=0Mcn,i,l,jψn,i(x)ψl,j(t)=ΨT(x)CΨ(t),
where Ψ(x) and Ψ(t) are 2k1(M1+1)×1 and 2k2(M2+1)×1 matrices, respectively, given by
(115)Ψ(x)≡[ψ0,0(x),ψ0,1(x),…,ψ0,M1(x),ψ1,0(x),…,ψ1,M1(x),…,ψ2,0(x),…,ψ2,M1(x),…,ψ2k1-1,0(x),…,ψ2k1-1,M1(x)],Ψ(t)≡[ψ0,0(t),ψ0,1(t),…,ψ0,M2(t),ψ1,0(t),…,ψ1,M2(t),…,ψ2,0(t),…,ψ2,M2(t),…,ψ2k2-1,0(t),…,ψ2k2-1,M2(t)]T.
Also, C is a 2k1(M1+1)×2k2(M2+1) matrix whose elements can be calculated from the following formula:
(116)C=∫01∫01ψn,i(x)ψl,j(t)f(x,t)dtdx
with n=0,1,…,2k1-1, i=0,…,M1, l=0,1,…,2k2-1, and j=0,…,M2.


#### 3.5.2. Application of Legendre Multiwavelet Galerkin Method on Nonlinear Fractional Partial Differential Equations [12]

Consider the nonlinear fractional partial differential equation(117)Dtαu=N(u)+g(x,t),m<α<m+1,m≥0,
with initial condition u(x,0)=f(x).Let(118)F(u)=Dtαu-N(u)-g(x,t).
A Galerkin approximation to (118) is constructed as follows. The approximation uNM are sought in the form of the truncated series
(119)uNM(x,t)={∑n=02k1-1∑i=0N∑l=02k2-1∑j=0Mtcn,i,l,j×ψn,i(x)ψl,j(t)+u(x,0),form=0∑n=02k1-1∑i=0N∑l=02k2-1∑j=0Mt2cn,i,l,j×ψn,i(x)ψl,j(t)+u(x,0)+tut(x,0),form=1,
where ψij are Legendre multiwavelet basis.The expansion coefficientscn,i,l,j are determined by Galerkin equations
(120)〈F(uNM),ψn,iψl,j〉=0,
where
(121)〈F(uNM),ψn,iψl,j〉=∫01∫01F(uNM)(x,t)ψn,i(x)ψl,j(t)dtdx.
Galerkin equations (120) gives a system of 2k1(N+1)×2k2(M+1) equations which can be solved for the elements of cn,i,l,j, i=0,1,2,…,N, j=0,1,2,…,M,  n=0,1,2,…,2k1-1, and l=0,1,2,…,2k2-1.


## 3.1. Haar Wavelets

Haar wavelet functions have been used from 1910 and were introduced by the Hungarian mathematician Alfred Haar. Haar wavelets (which are Daubechies wavelets of order 1) consist of piecewise constant functions on the real line that can take only three values, that is, 0, 1, and −1, and are therefore the simplest orthonormal wavelets with a compact support. Haar wavelet method is to be used due to the following features: being simpler and fast, flexible, and convenient, having small computational costs, and being computationally attractive. The Haar functions are a family of switched rectangular wave forms where amplitudes can differ from one function to another. These properties of Haar wavelets are utilized to reduce the computation of integral equations to some algebraic equations.Usually the Haar wavelets are defined for the intervalt∈[0,1) but in general case t∈[A,B], we divide the interval [A,B] into m equal subintervals, each of width Δt=(B-A)/m. In this case, the orthogonal set of Haar functions are defined in the interval [A,B] by [6]
(17)h0(t)={1t∈[A,B],0elsewhere,hi(t)={1,ζ1(i)≤t<ζ2(i)-1,ζ2(i)≤t<ζ3(i)0,otherwise,
where
(18)ζ1(i)=A+(k-12j)(B-A)=A+(k-12j)mΔt,ζ2(i)=A+(k-(1/2)2j)(B-A)=A+(k-(1/2)2j)mΔt,ζ3(i)=A+(k2j)(B-A)=A+(k2j)mΔt,
for i=1,2,…,m, m=2J and J is a positive integer which is called the maximum level of resolution. Here j and k represent the integer decomposition of the index i; that is, i=k+2j-1, 0≤j<i and 1≤k<2j+1.


### 3.1.1. Operational Matrix for General Order Integration [6]

The integration of theHm(t)=[h0(t),h1(t),…,hm-1(t)]T can be approximated by [16]
(19)∫0tHm(τ)dτ≅QHm(t),
where Q is called the Haar wavelet operational matrix of integration which is a square matrix of m-dimension. To derive the Haar wavelet operational matrix of the general order of integration, we recall the fractional integral of order α(>0) which is defined by Podlubny [3]:
(20)Jαf(t)=1Γ(α)∫0t(t-τ)α-1f(τ)dτ,t>0,α∈R+,
where R+ is the set of positive real numbers.The Haar wavelet operational matrixQα for integration of the general order α is given by
(21)QαHm(t)=JαHm(t)=[Jαh0(t),Jαh1(t),…,Jαhm-1(t)]T=[Qh0(t),Qh1(t),…,Qhm-1(t)]T,
where
(22)Qh0(t)={tαΓ(1+α),t∈[A,B],0,elsewhere,Qhi(t)={0,A≤t<ζ1(i),ϕ1,ζ1(i)≤t<ζ2(i),ϕ2,ζ2(i)≤t<ζ3(i),ϕ3,ζ3(i)≤t<B,
where
(23)ϕ1=(t-ζ1(i))αΓ(α+1),ϕ2=(t-ζ1(i))αΓ(α+1)-2(t-ζ2(i))αΓ(α+1),ϕ3=(t-ζ1(i))αΓ(α+1)-2(t-ζ2(i))αΓ(α+1)+(t-ζ3(i))αΓ(α+1),
for i=1,2,…,m, m=2J and J is a positive integer, called the maximum level of resolution. Here j and k represent the integer decomposition of the index i. That is, i=k+2j-1, 0≤j<i and 1≤k<2j+1.


### 3.1.2. Function Approximation

Any functiony(t)∈L2([0,1)) can be expanded into Haar wavelets by [6, 10, 16]
(24)y(t)=c0h0(t)+c1h1(t)+c2h2(t)+⋯,hhlhhhhwherecj=∫01y(t)hj(t)dt.
If y(t) is approximated as piecewise constant in each subinterval, the sum in (24) may be terminated after m terms and consequently we can write discrete version in the matrix form as
(25)Y≈∑i=0m-1cihi(tl)=CmTHm,
where Y and CmT are the m-dimensional row vectors.HereH is the Haar wavelet matrix of order m defined by H=[h0,h1,…,hm-1]T; that is,
(26)H=[h0h1⋮hm-1]=[h0,0h0,1⋯h0,m-1h1,0h1,1⋯h1,m-1⋮hm-1,0hm-1,1⋯hm-1,m-1],
where h0,h1,…,hm-1 are the discrete form of the Haar wavelet bases.The collocation points are given by(27)tl=A+(l-0.5)Δt,l=1,2,…,m.
Similarly, a function of two variables y(x,t)∈L2([0,1]×[0,1]) can be approximated by discrete version in the matrix form of Haar wavelets as
(28)Y(x,t)=HT(x)CH(t),
where Y(x,t) is the discrete form of y(x,t) and C is the coefficient matrix of Y, and it can be obtained by the following formula:
(29)C=(HT)-1·Y·H-1.
Since H is orthogonal, (29) becomes
(30)C=H·Y·H-1.


### 3.1.3. Application of Haar Wavelet on FDE [17]

In order to show the effectiveness of Haar wavelet method for solving fractional differential equations, we consider the following numerical example of variable coefficient fractional convection diffusion equation:(31)∂αu∂tα=-b(x)∂u∂x+a(x)∂βu∂xβ+q(x,t),
where 0≤x<1, 0<t≤1, a(x)>0, b(x)>0, 0<α≤1, 1<β≤2, q(x,t)∈C(D), and D=[0,1]×[0,1].Sinceu(x,t)∈L2(R), we suppose
(32)u(x,t)≈∑i=0m-1‍∑j=0m-1cijhi(x)hj(t).
Then we can obtain the discrete form of (32) by taking step Δ=1/m of x and t:
(33)U=HT(x)CH(t)
and then
(34)∂αu∂tα≈∂αU∂tα∂αu∂tα=HT(x)C∂αH(t)∂tα=HT(x)CQ-αH(t),∂βu∂xβ≈∂βU∂xβ∂αu∂tα=(∂β∂xβHT(x))CH(t)∂αu∂tα=(∂β∂xβH(x))TCH(t)=(Q-βH(x))TCH(t).
Here, q(x,t) is known function; discretizing it, we have
(35)D=(q(xi,tj)),i,j=0,1,2,…,m-1.
Coefficients a(x),b(x) can be dispersed into a(xi),b(xi), i=0,1,2,…,m-1.Let(36)A=[a(x0)0⋯00a(x1)⋯0⋮⋱⋮0⋯0a(xm-1)],B=[b(x0)0⋯00b(x1)⋯0⋮⋱⋮0⋯0b(xm-1)].
Substituting (33)–(36) in (31), we have
(37)HT(x)CQ-αH(t)=-B(∂∂xH(x))TCH(t)+A(∂β∂xβH(x))TCH(t)+D.
Equation (37) is a system of algebraic equations in a(xi) and b(xi). By solving this system of equations using mathematical software, the Haar wavelet coefficients a(xi) and b(xi) can be obtained.


## 3.1.1. Operational Matrix for General Order Integration [6]

The integration of theHm(t)=[h0(t),h1(t),…,hm-1(t)]T can be approximated by [16]
(19)∫0tHm(τ)dτ≅QHm(t),
where Q is called the Haar wavelet operational matrix of integration which is a square matrix of m-dimension. To derive the Haar wavelet operational matrix of the general order of integration, we recall the fractional integral of order α(>0) which is defined by Podlubny [3]:
(20)Jαf(t)=1Γ(α)∫0t(t-τ)α-1f(τ)dτ,t>0,α∈R+,
where R+ is the set of positive real numbers.The Haar wavelet operational matrixQα for integration of the general order α is given by
(21)QαHm(t)=JαHm(t)=[Jαh0(t),Jαh1(t),…,Jαhm-1(t)]T=[Qh0(t),Qh1(t),…,Qhm-1(t)]T,
where
(22)Qh0(t)={tαΓ(1+α),t∈[A,B],0,elsewhere,Qhi(t)={0,A≤t<ζ1(i),ϕ1,ζ1(i)≤t<ζ2(i),ϕ2,ζ2(i)≤t<ζ3(i),ϕ3,ζ3(i)≤t<B,
where
(23)ϕ1=(t-ζ1(i))αΓ(α+1),ϕ2=(t-ζ1(i))αΓ(α+1)-2(t-ζ2(i))αΓ(α+1),ϕ3=(t-ζ1(i))αΓ(α+1)-2(t-ζ2(i))αΓ(α+1)+(t-ζ3(i))αΓ(α+1),
for i=1,2,…,m, m=2J and J is a positive integer, called the maximum level of resolution. Here j and k represent the integer decomposition of the index i. That is, i=k+2j-1, 0≤j<i and 1≤k<2j+1.


## 3.1.2. Function Approximation

Any functiony(t)∈L2([0,1)) can be expanded into Haar wavelets by [6, 10, 16]
(24)y(t)=c0h0(t)+c1h1(t)+c2h2(t)+⋯,hhlhhhhwherecj=∫01y(t)hj(t)dt.
If y(t) is approximated as piecewise constant in each subinterval, the sum in (24) may be terminated after m terms and consequently we can write discrete version in the matrix form as
(25)Y≈∑i=0m-1cihi(tl)=CmTHm,
where Y and CmT are the m-dimensional row vectors.HereH is the Haar wavelet matrix of order m defined by H=[h0,h1,…,hm-1]T; that is,
(26)H=[h0h1⋮hm-1]=[h0,0h0,1⋯h0,m-1h1,0h1,1⋯h1,m-1⋮hm-1,0hm-1,1⋯hm-1,m-1],
where h0,h1,…,hm-1 are the discrete form of the Haar wavelet bases.The collocation points are given by(27)tl=A+(l-0.5)Δt,l=1,2,…,m.
Similarly, a function of two variables y(x,t)∈L2([0,1]×[0,1]) can be approximated by discrete version in the matrix form of Haar wavelets as
(28)Y(x,t)=HT(x)CH(t),
where Y(x,t) is the discrete form of y(x,t) and C is the coefficient matrix of Y, and it can be obtained by the following formula:
(29)C=(HT)-1·Y·H-1.
Since H is orthogonal, (29) becomes
(30)C=H·Y·H-1.


## 3.1.3. Application of Haar Wavelet on FDE [17]

In order to show the effectiveness of Haar wavelet method for solving fractional differential equations, we consider the following numerical example of variable coefficient fractional convection diffusion equation:(31)∂αu∂tα=-b(x)∂u∂x+a(x)∂βu∂xβ+q(x,t),
where 0≤x<1, 0<t≤1, a(x)>0, b(x)>0, 0<α≤1, 1<β≤2, q(x,t)∈C(D), and D=[0,1]×[0,1].Sinceu(x,t)∈L2(R), we suppose
(32)u(x,t)≈∑i=0m-1‍∑j=0m-1cijhi(x)hj(t).
Then we can obtain the discrete form of (32) by taking step Δ=1/m of x and t:
(33)U=HT(x)CH(t)
and then
(34)∂αu∂tα≈∂αU∂tα∂αu∂tα=HT(x)C∂αH(t)∂tα=HT(x)CQ-αH(t),∂βu∂xβ≈∂βU∂xβ∂αu∂tα=(∂β∂xβHT(x))CH(t)∂αu∂tα=(∂β∂xβH(x))TCH(t)=(Q-βH(x))TCH(t).
Here, q(x,t) is known function; discretizing it, we have
(35)D=(q(xi,tj)),i,j=0,1,2,…,m-1.
Coefficients a(x),b(x) can be dispersed into a(xi),b(xi), i=0,1,2,…,m-1.Let(36)A=[a(x0)0⋯00a(x1)⋯0⋮⋱⋮0⋯0a(xm-1)],B=[b(x0)0⋯00b(x1)⋯0⋮⋱⋮0⋯0b(xm-1)].
Substituting (33)–(36) in (31), we have
(37)HT(x)CQ-αH(t)=-B(∂∂xH(x))TCH(t)+A(∂β∂xβH(x))TCH(t)+D.
Equation (37) is a system of algebraic equations in a(xi) and b(xi). By solving this system of equations using mathematical software, the Haar wavelet coefficients a(xi) and b(xi) can be obtained.


## 3.2. Cubic B-Spline Basis Functions onH2(I)

LetI=[0,L] be an interval with L>4, and let H2(I) be a Sobolev space which contains functions with square integrable second derivatives, and the homogenous Sobolev space H02(I) can be defined by
(38)H02(I)={f′(0)f′(L)f(t)∈H2(I),f(0)=f′(0)=f(L)=f′(L)=0f′(0)f′(L)},
which is a Hilbert space equipped with inner product:
(39)〈f,g〉=∫If′′(t)g′′(t)dt.
Cai and Wang [18] present a multiresolution analysis (MRA) and a wavelet decomposition for H02(I) by constructing scaling spline functions
(40)φ(t)=16∑l=04(4l)(-1)l(t-l)+3,φb(t)=32t+2-1112t+3+32(t-1)+3-34(t-2)+3+16(t-3)+3,
and wavelet functions
(41)ψ(t)=-37φ(2t)+127φ(2t-1)-37φ(2t-2),ψb(t)=2413φb(2t)-613φ(2t),
which have compact supports, where
(42)t+n={tn,t>00,t≤0.
Then the scaling spaces V0⊂V1⊂V2⋯⊂V∞=H02(I) of MRA and a wavelet decomposition H02(I)=V0⊕W0⊕W1⊕⋯⊕Wj-1⊕⋯ corresponding to the MRA are obtained by defining
(43)Vj=span⁡{φj,k(t),φbj(t),φbj(L-t),0≤k≤nj-4},Wj=span⁡{ψj,k(t),k=-1,0,…,nj-2},
in which nj=2jL,
(44)φj,k(t)=φ(2jt-k),φbj(t)=φb(2jt),hhhhhhhhhhhhhhhhhhhlhh0≤k≤nj-4,ψj,k(t)=ψ(2jt-k),j≥0,k=0,1,…,nj-3,ψj,-1(t)=ψb(2jt),ψj,nj-2(t)=ψb(2j(L-t)).
For convenience, we set ψ-1,k(t)=φ0,k(t), 0≤k≤L-4, ψ-1,-1(t)=φb(t), ψ-1,L-3(t)=φb(L-t), and n-1=L-1. Let Bj={2-3j/2ψj,k(t)}k=-1nj-2, -1≤j<∞. Wang [19] derived that B=⋃j=-1∞Bj is an unconditional basis of H02(I), which turns out to be a basis of continuous space C0(I). For nonhomogeneous Sobolev space H2(I), Cai and Wang [18] introduced boundary spline functions
(45)η1(t)=(1-t)+3,η2(t)=2t+-3t+2+76t+3-43(t-1)+3+16(t-2)+3
to deal with the values of functions at boundary points.


### 3.2.1. Function Approximation

Any functionf(t)∈H02(I) can be uniquely expanded into cubic spline wavelet series by
(46)f(x)=∑j=-1∞∑k=1njdj,kψj,k-2(t),
with dj,k=∫If′′(t)(ψj,k-2*)(t)dt, where ψj,k*(t) are dual functions of ψj,k. Truncating the infinite series in (46) at J-1, we get
(47)fJ(t)=∑j=-1J-1∑k=1njdj,kψj,k-2(t).
From Wang [19] we have
(48)∥f(t)-fJ(t)∥H022⟶0asJ⟶+∞.
Hence any function f(t)∈H02(I) can be approximated by fJ(t) defined in (47). And any function f(t)∈H2(I) can be approximated by
(49)fJ(t)=Ib,Jf(t)+∑j=-1J-1∑k=1njdj,kψj,k-2(t),
and the approximation order is O(2-4J) if f(t) is sufficiently smooth, where
(50)Ib,Jf(t)=a1η1(2Jt)+a2η2(2Jt)+a3η2(2J(L-t))+a4η1(2J(L-t)).
Suppose N=2JL+3, and ΩJ(t) is a 1×N vector given by
(51)ΩJ(t)=[ψJ-1,nJ-2(t)η1(2Jt),η2(2Jt),η2(2J(L-t)),hhhη1(2J(L-t)),ψ-1,-1(t),ψ-1,0(t),…,hhhψ-1,n-1-2(t),ψ0,-1(t),ψ0,0(t),ψ0,1(t),…,ψ0,L-3(t),hhhψ0,n0-2(t),ψ1,-1(t),ψ1,0(t),ψ1,1(t),…,hhhψ1,2L-3(t),ψ1,n1-2(t),…,ψJ-1,-1(t),hhhψJ-1,0(t),ψJ-1,1(t),…,ψJ-1,nJ-3(t),ψJ-1,nJ-2(t)]≡[ω1(t),ω2(t),…,ωN(t)].fJ(t) defined in (49) can be rewritten as
(52)fJ(t)=∑k=1Nf^kωk(t)=ΩJ(t)f^,
where f^=[f^1,f^2,…,f^N]Tare the wavelet expansion coefficients which can be determined by interpolating at collocation points
(53)t-1(-1)=12J+1,tL+1(-1)=L-12J+1,tk(-1)=k,k=0,1,2,…,L,tk(j)=k+1.52j,-1≤k≤2jL-2,0≤j≤J-1.


### 3.2.2. Application of Cubic Spline Wavelets on FDE [4]

Consider the fractional differential equation of the form(54)Dαy(t)=-y(t)+f(t),0<t≤L,
with
(55)y(n)(0)=y0(n),n=0,1,2,…,m,hhhhhlm<α≤m+1,m∈N.
To solve problem (54), we approximate y(t) by
(56)y(t)≈∑k=1Nckωk(t)=Yj(t),
where ck=[c1,c2,…,cN]T is unknown. The α order fractional derivative of y(t) is approximated by
(57)Dαy(t)≈∑k=1NckDαωk(t)≡DαYj(t).
At collocation points,
(58)Yj(tk)=Y(tk),k=1,2,…,N,DαYj(tk)=DαY(tk),k=1,2,…,N.
Interpolating the fractional differential equation (54) by Y(tk) at all collocation points, we obtain
(59)DαYj(tk)=-Yj(tk)+f(tk),2≤k≤N,yj(n)(0)=y(n)(0),n=0,1,2,…,m.
Letting
(60)B1=[ω1(t2)ω2(t2)⋯ωN(t2)⋮⋱⋮ω1(tN)ω2(tN)⋯ωN(tN)],B2=[Dαω1(t2)Dαω2(t2)⋯DαωN(t2)⋮⋱⋮Dαω1(tN)Dαω2(tN)⋯DαωN(tN)],
where B1 and B2 are obtained analytically.Now (59) can be represented by
(61)Ay^=b,cny^=y0(n),
where A=B1+B2 and
(62)b=(f(t2),…,f(tN))T,cn=(ω1(n)(0),ω2(n)(0),…,ωN(n)(0)),n=0,1,2,…,m.
Let AmT=(AT,c0T,c1T,…,cmT), bmT=(bT,y0(0),y0(1),…,y0(m)); then (61) can be written as
(63)Amy^=bm.
Consequently the wavelet coefficient y^ can be obtained by solving linear equations (63).


## 3.2.1. Function Approximation

Any functionf(t)∈H02(I) can be uniquely expanded into cubic spline wavelet series by
(46)f(x)=∑j=-1∞∑k=1njdj,kψj,k-2(t),
with dj,k=∫If′′(t)(ψj,k-2*)(t)dt, where ψj,k*(t) are dual functions of ψj,k. Truncating the infinite series in (46) at J-1, we get
(47)fJ(t)=∑j=-1J-1∑k=1njdj,kψj,k-2(t).
From Wang [19] we have
(48)∥f(t)-fJ(t)∥H022⟶0asJ⟶+∞.
Hence any function f(t)∈H02(I) can be approximated by fJ(t) defined in (47). And any function f(t)∈H2(I) can be approximated by
(49)fJ(t)=Ib,Jf(t)+∑j=-1J-1∑k=1njdj,kψj,k-2(t),
and the approximation order is O(2-4J) if f(t) is sufficiently smooth, where
(50)Ib,Jf(t)=a1η1(2Jt)+a2η2(2Jt)+a3η2(2J(L-t))+a4η1(2J(L-t)).
Suppose N=2JL+3, and ΩJ(t) is a 1×N vector given by
(51)ΩJ(t)=[ψJ-1,nJ-2(t)η1(2Jt),η2(2Jt),η2(2J(L-t)),hhhη1(2J(L-t)),ψ-1,-1(t),ψ-1,0(t),…,hhhψ-1,n-1-2(t),ψ0,-1(t),ψ0,0(t),ψ0,1(t),…,ψ0,L-3(t),hhhψ0,n0-2(t),ψ1,-1(t),ψ1,0(t),ψ1,1(t),…,hhhψ1,2L-3(t),ψ1,n1-2(t),…,ψJ-1,-1(t),hhhψJ-1,0(t),ψJ-1,1(t),…,ψJ-1,nJ-3(t),ψJ-1,nJ-2(t)]≡[ω1(t),ω2(t),…,ωN(t)].fJ(t) defined in (49) can be rewritten as
(52)fJ(t)=∑k=1Nf^kωk(t)=ΩJ(t)f^,
where f^=[f^1,f^2,…,f^N]Tare the wavelet expansion coefficients which can be determined by interpolating at collocation points
(53)t-1(-1)=12J+1,tL+1(-1)=L-12J+1,tk(-1)=k,k=0,1,2,…,L,tk(j)=k+1.52j,-1≤k≤2jL-2,0≤j≤J-1.


## 3.2.2. Application of Cubic Spline Wavelets on FDE [4]

Consider the fractional differential equation of the form(54)Dαy(t)=-y(t)+f(t),0<t≤L,
with
(55)y(n)(0)=y0(n),n=0,1,2,…,m,hhhhhlm<α≤m+1,m∈N.
To solve problem (54), we approximate y(t) by
(56)y(t)≈∑k=1Nckωk(t)=Yj(t),
where ck=[c1,c2,…,cN]T is unknown. The α order fractional derivative of y(t) is approximated by
(57)Dαy(t)≈∑k=1NckDαωk(t)≡DαYj(t).
At collocation points,
(58)Yj(tk)=Y(tk),k=1,2,…,N,DαYj(tk)=DαY(tk),k=1,2,…,N.
Interpolating the fractional differential equation (54) by Y(tk) at all collocation points, we obtain
(59)DαYj(tk)=-Yj(tk)+f(tk),2≤k≤N,yj(n)(0)=y(n)(0),n=0,1,2,…,m.
Letting
(60)B1=[ω1(t2)ω2(t2)⋯ωN(t2)⋮⋱⋮ω1(tN)ω2(tN)⋯ωN(tN)],B2=[Dαω1(t2)Dαω2(t2)⋯DαωN(t2)⋮⋱⋮Dαω1(tN)Dαω2(tN)⋯DαωN(tN)],
where B1 and B2 are obtained analytically.Now (59) can be represented by
(61)Ay^=b,cny^=y0(n),
where A=B1+B2 and
(62)b=(f(t2),…,f(tN))T,cn=(ω1(n)(0),ω2(n)(0),…,ωN(n)(0)),n=0,1,2,…,m.
Let AmT=(AT,c0T,c1T,…,cmT), bmT=(bT,y0(0),y0(1),…,y0(m)); then (61) can be written as
(63)Amy^=bm.
Consequently the wavelet coefficient y^ can be obtained by solving linear equations (63).


## 3.3. Chebyshev Wavelets

The Legendre wavelet and Chebyshev wavelet are constructed from their corresponding polynomials. These wavelets are useful tools in the numerical computations. The second kind Chebyshev polynomials have many good properties and are widely applied in different disciplines.Wavelets constitute a family of functions constructed from dilation and translation of a single functionψ(t) called the mother wavelet. When the dilation parameter a and the translation parameter b vary continuously, we have the following family of continuous wavelets as [13]
(64)ψa,b(t)=|a|-1/2ψ(t-ba),a,b∈R,a≠0.
If we restrict the parameters a and b to discrete values as a=a0-k, b=nb0a0-k; a0>1, b0>0, where n and k are positive integer, the following family of discrete wavelets are defined as
(65)ψk,n(t)=|a0|k/2ψ(a0kt-nb0),
where ψk,n(t) forms a wavelet basis for L2(R). In particular, when a0=2 and b0=1, then ψk,n(t) forms an orthonormal basis; that is,
(66)〈ψj,k(t),ψl,m(t)〉=δjlδkm.


### 3.3.1. Chebyshev Wavelet and Function Approximation

Chebyshev waveletsψn,m(t), on the interval [0,1), are defined as [20]
(67)ψn,m(t)={2k/2T~m(2kt-2n+1),n-12k-1≤t<n2k-10,otherwise,
where
(68)T~m={1π,m=0,2πTm(t),m>0,
and m=0,1,…,M-1, n=1,2,…,2k-1,k is any positive integer, and Tm(t) are Chebyshev polynomials of the first kind of degree m which are orthogonal with respect to the weight function ω(t)=(1/1-t2), on the interval [-1,1], and Tm(t) can be determined by the following recurrence formula:
(69)T0(t)=1,T1(t)=t,sTm+1(t)=2tTm(t)-Tm-1(t),m=1,2,3,….
A function f(t) defined over [0, 1) may be expanded as
(70)f(t)=∑n=1∞‍∑m=0∞cn,mψn,m(t),
where cnm=〈f(t),ψn,m(t)〉.If the infinite series in (70) is truncated, then (70) can be written as
(71)f(t)≈∑n=12k-1‍∑m=0M-1cn,mψn,m(t)=CTΨ(t),
where C and Ψ(t) are 2k-1M×1 matrices given by
(72)C≡[c1,0,c1,1,…,c1,M-1,c2,0,…,c2,M-1,…,c2k-1,0,…,c2k-1,M-1]T,(73)Ψ(t)≡[ψ1,0,ψ1,1,…,ψ1,M-1,ψ2,0,…,ψ2,M-1,…,ψ2k-1,0,…,ψ2k-1,M-1]T.
Taking the collocation points as follows:
(74)ti=2i-12kM,i=1,2,…,2k-1M,
we define the Chebyshev wavelet matrix Φm×m as
(75)Φm×m≡[Ψ(12m)Ψ(32m)⋯Ψ(2m-12m)].


### 3.3.2. Operational Matrix of the Fractional Integration

The integration of the vectorΨ(t) defined in (73) can be obtained as
(76)∫0tΨ(τ)dτ≈PΨ(t),
where P is the 2k-1M×2k-1M operational matrix for integration [13].We know the Riemann-Liouville fractional integral operatorJα, of a function f(t), is defined as [3]
(77)Jαf(t)=1Γ(α)∫0t(t-τ)α-1f(τ)dτ=1Γ(α)tα-1*f(t),
where α∈R+(α>0) is the order of integration, Γ(·) is the well-known gamma function, and tα-1*f(t) denotes the convolution product of tα-1 and f(t). Now if f(t) is expanded in Chebyshev wavelets, the Riemann-Liouville fractional integration becomes
(78)Jαf(t)=1Γ(α)tα-1*f(t)≅CT1Γ(α){tα-1*Ψ(t)}.
Thus if tα-1*f(t) can be integrated and then expanded in Chebyshev wavelets, the Riemann-Liouville fractional integration can be solved via Chebyshev wavelets.Also we define am-set of Block-Pulse functions (BPF) as [13]
(79)bi(t)={1,im≤t<i+1m0,otherwise,
where i=0,1,2,…,(m-1).The functionsbi(t) are disjoint and orthogonal; that is,
(80)bi(t)bl(t)={0,i≠lbi(t),i=l,∫01bi(τ)bl(τ)dτ={0,i≠l1m,i=l.
Similarly, Chebyshev wavelets may be expanded into an m-term block pulse functions (BPF) as
(81)Ψm′(t)=Φm′×m′Bm′(t),
where m′=2k-1M and Bm′(t)≡[b0(t)b1(t)⋯bi(t)⋯bm-1(t)]T.In [21], Kilicman and Al Zhour have proposed the block pulse operational matrix of the fractional integration Fα as follows:
(82)Jα(Bm′(t))≈FαBm′(t),
where
(83)Fα=1mα1Γ(α+2)[1ξ1ξ2···ξm-101ξ1···ξm-2001···ξm-3······000··00001]
with ξk=(k+1)α+1-2kα+1+(k-1)α+1.Next we derive the Chebyshev wavelet operational matrix of the fractional integration. Let(84)Jα(Ψm(t))≈Pm′×m′αΨm(t),
where matrix Pm′×m′α is called the Chebyshev wavelet operational matrix of the fractional integration.Using (81) and (82), we have
(85)JαΨm′(t)≈(JαΦm′×m′Bm′)(t)=Φm′×m′(JαBm′)(t)≈Φm′×m′FαBm′(t).
From (84) and (85) we get
(86)Pm′×m′αΨm′(t)=Φm′×m′FαBm′(t).
Then, the Chebyshev wavelet operational matrix of the fractional integration Pm′×m′α is given by
(87)Pm′×m′α=Φm′×m′FαΦm′×m′-1.
Thus, we derive the operational matrix of fractional integration for the second kind Chebyshev wavelet.


### 3.3.3. Application of Chebyshev Wavelet on FDE [22]

In this section, we use Chebyshev wavelet operational matrices of the fractional integration to solve nonlinear fractional differential equation. We consider fractional Riccati equation(88)Dαy(t)+[y(t)]2=1,0<α≤1,
subject to initial state y(0)=0.Let(89)Dαy(t)=CTΨ(t),
and then
(90)y(t)=CTPm′×m′αΨ(t).
Since Ψ(t)=Φm′×m′Bm′(t), we have
(91)y(t)=CTPm′×m′αΦm′×m′Bm′(t).
Let CTPm′×m′αΦm′×m′=[a1,a2,…,am′], and, using (79), we get
(92)[y(t)]2=[a1b1(t)+a2b2(t)+⋯+am′bm′(t)]2=[a12,a22,…,am′2]Bm′(t).
Substituting (89) and (92) in (88), we have
(93)CTΦm′×m′Bm′(t)+[a12,a22,…,am′2]Bm′(t)-[1,1,…,1]Bm′(t)=0.
Solving the above nonlinear system of algebraic equations, we can find the vector C and consequently the solution for y(t).


## 3.3.1. Chebyshev Wavelet and Function Approximation

Chebyshev waveletsψn,m(t), on the interval [0,1), are defined as [20]
(67)ψn,m(t)={2k/2T~m(2kt-2n+1),n-12k-1≤t<n2k-10,otherwise,
where
(68)T~m={1π,m=0,2πTm(t),m>0,
and m=0,1,…,M-1, n=1,2,…,2k-1,k is any positive integer, and Tm(t) are Chebyshev polynomials of the first kind of degree m which are orthogonal with respect to the weight function ω(t)=(1/1-t2), on the interval [-1,1], and Tm(t) can be determined by the following recurrence formula:
(69)T0(t)=1,T1(t)=t,sTm+1(t)=2tTm(t)-Tm-1(t),m=1,2,3,….
A function f(t) defined over [0, 1) may be expanded as
(70)f(t)=∑n=1∞‍∑m=0∞cn,mψn,m(t),
where cnm=〈f(t),ψn,m(t)〉.If the infinite series in (70) is truncated, then (70) can be written as
(71)f(t)≈∑n=12k-1‍∑m=0M-1cn,mψn,m(t)=CTΨ(t),
where C and Ψ(t) are 2k-1M×1 matrices given by
(72)C≡[c1,0,c1,1,…,c1,M-1,c2,0,…,c2,M-1,…,c2k-1,0,…,c2k-1,M-1]T,(73)Ψ(t)≡[ψ1,0,ψ1,1,…,ψ1,M-1,ψ2,0,…,ψ2,M-1,…,ψ2k-1,0,…,ψ2k-1,M-1]T.
Taking the collocation points as follows:
(74)ti=2i-12kM,i=1,2,…,2k-1M,
we define the Chebyshev wavelet matrix Φm×m as
(75)Φm×m≡[Ψ(12m)Ψ(32m)⋯Ψ(2m-12m)].


## 3.3.2. Operational Matrix of the Fractional Integration

The integration of the vectorΨ(t) defined in (73) can be obtained as
(76)∫0tΨ(τ)dτ≈PΨ(t),
where P is the 2k-1M×2k-1M operational matrix for integration [13].We know the Riemann-Liouville fractional integral operatorJα, of a function f(t), is defined as [3]
(77)Jαf(t)=1Γ(α)∫0t(t-τ)α-1f(τ)dτ=1Γ(α)tα-1*f(t),
where α∈R+(α>0) is the order of integration, Γ(·) is the well-known gamma function, and tα-1*f(t) denotes the convolution product of tα-1 and f(t). Now if f(t) is expanded in Chebyshev wavelets, the Riemann-Liouville fractional integration becomes
(78)Jαf(t)=1Γ(α)tα-1*f(t)≅CT1Γ(α){tα-1*Ψ(t)}.
Thus if tα-1*f(t) can be integrated and then expanded in Chebyshev wavelets, the Riemann-Liouville fractional integration can be solved via Chebyshev wavelets.Also we define am-set of Block-Pulse functions (BPF) as [13]
(79)bi(t)={1,im≤t<i+1m0,otherwise,
where i=0,1,2,…,(m-1).The functionsbi(t) are disjoint and orthogonal; that is,
(80)bi(t)bl(t)={0,i≠lbi(t),i=l,∫01bi(τ)bl(τ)dτ={0,i≠l1m,i=l.
Similarly, Chebyshev wavelets may be expanded into an m-term block pulse functions (BPF) as
(81)Ψm′(t)=Φm′×m′Bm′(t),
where m′=2k-1M and Bm′(t)≡[b0(t)b1(t)⋯bi(t)⋯bm-1(t)]T.In [21], Kilicman and Al Zhour have proposed the block pulse operational matrix of the fractional integration Fα as follows:
(82)Jα(Bm′(t))≈FαBm′(t),
where
(83)Fα=1mα1Γ(α+2)[1ξ1ξ2···ξm-101ξ1···ξm-2001···ξm-3······000··00001]
with ξk=(k+1)α+1-2kα+1+(k-1)α+1.Next we derive the Chebyshev wavelet operational matrix of the fractional integration. Let(84)Jα(Ψm(t))≈Pm′×m′αΨm(t),
where matrix Pm′×m′α is called the Chebyshev wavelet operational matrix of the fractional integration.Using (81) and (82), we have
(85)JαΨm′(t)≈(JαΦm′×m′Bm′)(t)=Φm′×m′(JαBm′)(t)≈Φm′×m′FαBm′(t).
From (84) and (85) we get
(86)Pm′×m′αΨm′(t)=Φm′×m′FαBm′(t).
Then, the Chebyshev wavelet operational matrix of the fractional integration Pm′×m′α is given by
(87)Pm′×m′α=Φm′×m′FαΦm′×m′-1.
Thus, we derive the operational matrix of fractional integration for the second kind Chebyshev wavelet.


## 3.3.3. Application of Chebyshev Wavelet on FDE [22]

In this section, we use Chebyshev wavelet operational matrices of the fractional integration to solve nonlinear fractional differential equation. We consider fractional Riccati equation(88)Dαy(t)+[y(t)]2=1,0<α≤1,
subject to initial state y(0)=0.Let(89)Dαy(t)=CTΨ(t),
and then
(90)y(t)=CTPm′×m′αΨ(t).
Since Ψ(t)=Φm′×m′Bm′(t), we have
(91)y(t)=CTPm′×m′αΦm′×m′Bm′(t).
Let CTPm′×m′αΦm′×m′=[a1,a2,…,am′], and, using (79), we get
(92)[y(t)]2=[a1b1(t)+a2b2(t)+⋯+am′bm′(t)]2=[a12,a22,…,am′2]Bm′(t).
Substituting (89) and (92) in (88), we have
(93)CTΦm′×m′Bm′(t)+[a12,a22,…,am′2]Bm′(t)-[1,1,…,1]Bm′(t)=0.
Solving the above nonlinear system of algebraic equations, we can find the vector C and consequently the solution for y(t).


## 3.4. Legendre Wavelets

Legendre waveletψn,m=ψ(k,n,m,t) have four arguments: n=1,…,2k-1,k can be assumed as any positive integer, m is order of Legendre polynomials, and t denotes the time.The Legendre polynomials of orderm, denoted by Lm(t), are defined on the interval [−1, 1] and can be determined with the help of the following recurrence formulae [11]:
(94)L0(t)=1,L1(t)=t,(95)Lm+1(t)=2m+1m+1tLm(t)-mm+1Lm-1(t),hhhhhhhhhhhhhhhhhhhhhhm=1,2,3,….
The Legendre wavelets are defined on interval [0, 1) by
(96)ψn,m(t)={(m+12)1/22k/2Lm(2kt-n^),n^-12k≤t<n^+12k,0,elsewhere
for k=2,3,…n^=2n-1, n=1,2,3,…,2k-1, m=0,1,…,M-1 is the order of the Legendre polynomials, and M is a fixed positive integer. The set of Legendre wavelets forms an orthogonal basis of L2(R).


### 3.4.1. Function Approximation

A functionf(t)∈L2[0,1] can be expanded in terms of Legendre wavelet as [23]
(97)f(t)=∑n=1∞‍∑m=0∞cn,mψn,m(t),
where cn,m=〈f(t),ψn,m(t)〉.If the infinite series in (97) is truncated, then (97) can be written as
(98)f(t)≈∑n=12k-1‍∑m=0M-1cn,mψn,m(t)=CTΨ(t),
where C and Ψ(x) are 2k-1M×1 matrices given by
(99)C≡[c1,0,c1,1,…,c1,M-1,c2,0,…,c2,M-1,…,c2k-1,0,…,c2k-1,M-1]T,Ψ(t)≡[ψ1,0,ψ1,1,…,ψ1,M-1,ψ2,0,…,ψ2,M-1,…,ψ2k-1,0,…,ψ2k-1,M-1]T.
Similarly, a function of two variables k(x,t)∈L2([0,1]×[0,1]) can be approximated by Legendre wavelet series as
(100)k(x,t)≈ΨT(t)KΨ(x),
where K=[kij] is (2k-1M×2k-1M) matrix, with
(101)ki,j=〈ψi(t),〈k(x,t),ψj(x)〉〉.


### 3.4.2. Application of Legendre Wavelet on FDE

In order to show the effectiveness of Legendre wavelet method for solving fractional differential equations, we consider the following numerical example [11]:
(102)Dαy(t)+y(t)=0,0<α≤2,
such that
(103)y(0)=1,y′(0)=0.
The corresponding integral representation for (102) and (103) is
(104)y=-Jαy(t)+y(0)+y′(0)t.
Let y be approximated using Legendre wavelet as
(105)y=CTΨ(t),
and then
(106)Jαy=CTJαΨ(t)=CTPm^×m^αΨ(t),hhhhhhhlhhhhwherem^=2k-1M.
Substituting (105) and (106) in (104), we have the following system of algebraic equations:
(107)CTΨ(t)+CTPm^×m^αΨ(t)=Y0,whereY0=y(0)+y′(0)t.
By solving the linear system (107), we can find the vector C and hence consequently the solution y(t).


## 3.4.1. Function Approximation

A functionf(t)∈L2[0,1] can be expanded in terms of Legendre wavelet as [23]
(97)f(t)=∑n=1∞‍∑m=0∞cn,mψn,m(t),
where cn,m=〈f(t),ψn,m(t)〉.If the infinite series in (97) is truncated, then (97) can be written as
(98)f(t)≈∑n=12k-1‍∑m=0M-1cn,mψn,m(t)=CTΨ(t),
where C and Ψ(x) are 2k-1M×1 matrices given by
(99)C≡[c1,0,c1,1,…,c1,M-1,c2,0,…,c2,M-1,…,c2k-1,0,…,c2k-1,M-1]T,Ψ(t)≡[ψ1,0,ψ1,1,…,ψ1,M-1,ψ2,0,…,ψ2,M-1,…,ψ2k-1,0,…,ψ2k-1,M-1]T.
Similarly, a function of two variables k(x,t)∈L2([0,1]×[0,1]) can be approximated by Legendre wavelet series as
(100)k(x,t)≈ΨT(t)KΨ(x),
where K=[kij] is (2k-1M×2k-1M) matrix, with
(101)ki,j=〈ψi(t),〈k(x,t),ψj(x)〉〉.


## 3.4.2. Application of Legendre Wavelet on FDE

In order to show the effectiveness of Legendre wavelet method for solving fractional differential equations, we consider the following numerical example [11]:
(102)Dαy(t)+y(t)=0,0<α≤2,
such that
(103)y(0)=1,y′(0)=0.
The corresponding integral representation for (102) and (103) is
(104)y=-Jαy(t)+y(0)+y′(0)t.
Let y be approximated using Legendre wavelet as
(105)y=CTΨ(t),
and then
(106)Jαy=CTJαΨ(t)=CTPm^×m^αΨ(t),hhhhhhhlhhhhwherem^=2k-1M.
Substituting (105) and (106) in (104), we have the following system of algebraic equations:
(107)CTΨ(t)+CTPm^×m^αΨ(t)=Y0,whereY0=y(0)+y′(0)t.
By solving the linear system (107), we can find the vector C and hence consequently the solution y(t).


## 3.5. Legendre Multiwavelets Galerkin Method

In this section, the application of Legendre multiwavelet Galerkin method for providing approximate solutions for initial value problems of fractional differential equations has been discussed. Legendre multiwaveletψn,m(t)=ψ(k,n,m,t) have four arguments: n=0,1,…,2k-1, k can be assumed as any positive integer, m is order of Legendre polynomials, and t denotes the normalized time. The Legendre multiwavelets are defined on the interval [0, 1) by
(108)ψn,m(t)={(2m+1)1/22k/2Pm(2kt-n),n2k≤t<n+12k,0,elsewhere,
where m=0,1,…,M-1 and n=1,2,3,…,2k-1. The coefficient (2m+1)1/2 is for orthonormality and Pm(t) are the well-known shifted Legendre polynomials of order m which are defined on the interval [0,1] and can be determined with the help of the following recurrence formula:
(109)P0(t)=1,P1(t)=2t-1,(110)Pm+1(t)=2m+1m+1(2t-1)Pm(t)-mm+1Pm-1(t),hhhhhhhhhhhhhhhhhlhhhm=1,2,3,….
The two-dimensional Legendre multiwavelets are defined as
(111)ψn1,m1,n2,m2(x,t)={APm1(2k1x-n1)Pm2(2k2t-n2),n12k1≤x<n1+12k1,n22k2≤t<n2+12k20,elsewhere,
where
(112)A=(2m1+1)(2m2+1)2(k1+k2)/2.n1 and n2 are defined similarly to n, k1 and k2 are any positive integer, m1 and m2 are the order for Legendre polynomials, and ψn1,m1,n2,m2(x,t) forms a basis for L2([0,1]×[0,1]).


### 3.5.1. Function Approximation

A functionf(x,t) defined over [0,1]×[0,1] can be expanded in terms of Legendre multiwavelet as [12]
(113)f(x,t)=∑n=0∞∑i=0∞‍∑l=0∞‍∑j=0∞cn,i,l,jψn,i(x)ψl,j(t).
If the infinite series in (113) is truncated, then (113) can be written as
(114)f(x,t)=∑n=02k1-1‍∑i=0N∑l=02k2-1∑j=0Mcn,i,l,jψn,i(x)ψl,j(t)=ΨT(x)CΨ(t),
where Ψ(x) and Ψ(t) are 2k1(M1+1)×1 and 2k2(M2+1)×1 matrices, respectively, given by
(115)Ψ(x)≡[ψ0,0(x),ψ0,1(x),…,ψ0,M1(x),ψ1,0(x),…,ψ1,M1(x),…,ψ2,0(x),…,ψ2,M1(x),…,ψ2k1-1,0(x),…,ψ2k1-1,M1(x)],Ψ(t)≡[ψ0,0(t),ψ0,1(t),…,ψ0,M2(t),ψ1,0(t),…,ψ1,M2(t),…,ψ2,0(t),…,ψ2,M2(t),…,ψ2k2-1,0(t),…,ψ2k2-1,M2(t)]T.
Also, C is a 2k1(M1+1)×2k2(M2+1) matrix whose elements can be calculated from the following formula:
(116)C=∫01∫01ψn,i(x)ψl,j(t)f(x,t)dtdx
with n=0,1,…,2k1-1, i=0,…,M1, l=0,1,…,2k2-1, and j=0,…,M2.


### 3.5.2. Application of Legendre Multiwavelet Galerkin Method on Nonlinear Fractional Partial Differential Equations [12]

Consider the nonlinear fractional partial differential equation(117)Dtαu=N(u)+g(x,t),m<α<m+1,m≥0,
with initial condition u(x,0)=f(x).Let(118)F(u)=Dtαu-N(u)-g(x,t).
A Galerkin approximation to (118) is constructed as follows. The approximation uNM are sought in the form of the truncated series
(119)uNM(x,t)={∑n=02k1-1∑i=0N∑l=02k2-1∑j=0Mtcn,i,l,j×ψn,i(x)ψl,j(t)+u(x,0),form=0∑n=02k1-1∑i=0N∑l=02k2-1∑j=0Mt2cn,i,l,j×ψn,i(x)ψl,j(t)+u(x,0)+tut(x,0),form=1,
where ψij are Legendre multiwavelet basis.The expansion coefficientscn,i,l,j are determined by Galerkin equations
(120)〈F(uNM),ψn,iψl,j〉=0,
where
(121)〈F(uNM),ψn,iψl,j〉=∫01∫01F(uNM)(x,t)ψn,i(x)ψl,j(t)dtdx.
Galerkin equations (120) gives a system of 2k1(N+1)×2k2(M+1) equations which can be solved for the elements of cn,i,l,j, i=0,1,2,…,N, j=0,1,2,…,M,  n=0,1,2,…,2k1-1, and l=0,1,2,…,2k2-1.


## 3.5.1. Function Approximation

A functionf(x,t) defined over [0,1]×[0,1] can be expanded in terms of Legendre multiwavelet as [12]
(113)f(x,t)=∑n=0∞∑i=0∞‍∑l=0∞‍∑j=0∞cn,i,l,jψn,i(x)ψl,j(t).
If the infinite series in (113) is truncated, then (113) can be written as
(114)f(x,t)=∑n=02k1-1‍∑i=0N∑l=02k2-1∑j=0Mcn,i,l,jψn,i(x)ψl,j(t)=ΨT(x)CΨ(t),
where Ψ(x) and Ψ(t) are 2k1(M1+1)×1 and 2k2(M2+1)×1 matrices, respectively, given by
(115)Ψ(x)≡[ψ0,0(x),ψ0,1(x),…,ψ0,M1(x),ψ1,0(x),…,ψ1,M1(x),…,ψ2,0(x),…,ψ2,M1(x),…,ψ2k1-1,0(x),…,ψ2k1-1,M1(x)],Ψ(t)≡[ψ0,0(t),ψ0,1(t),…,ψ0,M2(t),ψ1,0(t),…,ψ1,M2(t),…,ψ2,0(t),…,ψ2,M2(t),…,ψ2k2-1,0(t),…,ψ2k2-1,M2(t)]T.
Also, C is a 2k1(M1+1)×2k2(M2+1) matrix whose elements can be calculated from the following formula:
(116)C=∫01∫01ψn,i(x)ψl,j(t)f(x,t)dtdx
with n=0,1,…,2k1-1, i=0,…,M1, l=0,1,…,2k2-1, and j=0,…,M2.


## 3.5.2. Application of Legendre Multiwavelet Galerkin Method on Nonlinear Fractional Partial Differential Equations [12]

Consider the nonlinear fractional partial differential equation(117)Dtαu=N(u)+g(x,t),m<α<m+1,m≥0,
with initial condition u(x,0)=f(x).Let(118)F(u)=Dtαu-N(u)-g(x,t).
A Galerkin approximation to (118) is constructed as follows. The approximation uNM are sought in the form of the truncated series
(119)uNM(x,t)={∑n=02k1-1∑i=0N∑l=02k2-1∑j=0Mtcn,i,l,j×ψn,i(x)ψl,j(t)+u(x,0),form=0∑n=02k1-1∑i=0N∑l=02k2-1∑j=0Mt2cn,i,l,j×ψn,i(x)ψl,j(t)+u(x,0)+tut(x,0),form=1,
where ψij are Legendre multiwavelet basis.The expansion coefficientscn,i,l,j are determined by Galerkin equations
(120)〈F(uNM),ψn,iψl,j〉=0,
where
(121)〈F(uNM),ψn,iψl,j〉=∫01∫01F(uNM)(x,t)ψn,i(x)ψl,j(t)dtdx.
Galerkin equations (120) gives a system of 2k1(N+1)×2k2(M+1) equations which can be solved for the elements of cn,i,l,j, i=0,1,2,…,N, j=0,1,2,…,M,  n=0,1,2,…,2k1-1, and l=0,1,2,…,2k2-1.


## 4. Conclusion and Discussion

In this work, we derive the operational matrices of the fractional integration for Haar wavelet, Legendre wavelet, and Chebyshev wavelets. We have examined several wavelet methods to solve fractional differential equations. In this present paper, we have applied cubic B-spline wavelets for solving fractional differential equations. In order to increase the accuracy of the approximate solution, it is necessary to apply higher order spline wavelet method. The Legendre multiwavelet method can be applied for providing approximate solutions for initial value problems of fractional nonlinear partial differential equations. The main characteristic of this approach is using the properties of Legendre multiwavelet together with the Galerkin method to reduce the nonlinear fractional order partial differential equations (NFPDEs) to the nonlinear system of algebraic equations. Similarly, the Chebyshev wavelets, constructed from the corresponding Chebyshev polynomials, are quite efficient in dealing with the nonlinear FDE like Riccati equation [22]. In this present analysis, a Legendre wavelet operational matrix of fractional order integration is obtained and is used to solve fractional differential equations numerically. It is worth mentioning that results agree well with exact solutions even for small values of k and M. The method is very convenient for solving initial value problems as well as boundary value problems.Using these wavelet methods, the fractional differential equations have been reduced to a system of algebraic equations and this system can be easily solved by any usual methods. Haar wavelet method can also be applied to the fractional differential equations by reducing into a system of algebraic equations. These methods give more accuracy if we increase the order or the level of resolution. The approximate solutions by these aforesaid methods highly agree with exact solutions.



---

*Source: 140453-2014-05-27.xml*