Wavelet Methods for Solving Fractional Order Differential Equations

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.

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 𝜓(𝑥): an oscillatory function 𝜓(𝑥)∈𝐿²(ℜ) with zero mean is a wavelet if it has the following desirable attributes:(1)smoothness: 𝜓(𝑥) is 𝑛 times differentiable and its derivatives are continuous;(2)localization: 𝜓(𝑥) is well localized both in time and frequency domains; that is, 𝜓(𝑥) 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 𝜓(𝑥); that is, ^∞_−∞𝑥^𝑘Ψ(𝑥)𝑑𝑑𝑥=0orequivalently^𝑘𝑑𝜔^𝑘Ψ(𝜔)=0for𝑘=0,1,…,𝑛 (1) 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) suggests that |Ψ(𝜔)| decay at least as |𝜔|⁻¹ or |𝑥|^𝜀−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 {𝑉_𝑗}_𝑗∈𝑍 is said to be MRA of 𝐿²(ℜ) if it possess the following properties: 𝑉_𝑗⊂𝑉_𝑗+1,∀𝑗∈𝑍,_𝑗∈𝑍𝑉_𝑗isdensein𝐿²(𝑅),_𝑗∈𝑍𝑉_𝑗=𝜙,𝑓(𝑥)∈𝑉_𝑗⟺𝑓(2𝑥)∈𝑉_𝑗+1,∀𝑗∈𝑍, (3)(4)(5)(6) where 𝑍 denotes the set of integers. Properties (3)–(5) state that {𝑉_𝑗}_𝑗∈𝑍 is a nested sequence of subspaces that effectively covers 𝐿²(ℜ). 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 𝑉_𝑗. A function 𝜑∈𝐿²(𝑅) is called a scaling function if it generates the nested sequence of subspaces 𝑉_𝑗 and satisfies the dilation equation; namely, 𝜑(𝑥)=_𝑘𝑝_𝑘𝜑(𝑎𝑥−𝑘), (7) with 𝑝_𝑘∈𝑙² and 𝑎 being any rational number.

For each scale 𝑗, since 𝑉_𝑗⊂𝑉_𝑗+1, there exists a unique orthogonal complementary subspace 𝑊_𝑗 of 𝑉_𝑗 in 𝑉_𝑗+1. This subspace 𝑊_𝑗 is called wavelet subspace and is generated by 𝜓_𝑗,𝑘=𝜓(2^𝑗𝑥−𝑘), where 𝜓∈𝐿² is called the wavelet. From the above discussion, these results follow easily: 𝑉_𝑗1𝑉_𝑗2=𝑉_𝑗2,𝑗₁>𝑗₂,𝑊_𝑗1𝑊_𝑗2=0,𝑗₁≠𝑗₂,𝑉_𝑗1𝑊_𝑗2=0,𝑗₁≤𝑗₂. (8)

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 𝐽^𝛼 of a function 𝑓(𝑡) is defined as [3]: 𝐽^𝛼𝑓(𝑡)=1Γ(𝛼)^𝑡₀(𝑡−𝜏)^𝛼−1𝑓(𝜏)𝑑𝜏,𝛼>0,𝛼∈ℜ⁺, (9) where Γ(⋅) is the well-known gamma function, and some properties of the operator 𝐽^𝛼 are as follows: 𝐽^𝛼𝐽^𝛽𝑓(𝑡)=𝐽^𝛼+𝛽𝑓(𝑡),,𝐽𝛼>0,𝛽>0_𝛼𝑡^𝛾=Γ1+𝛾Γ𝑡1+𝛾+𝛼^𝛼+𝛾,𝛾>−1. (10)

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 ₀𝐷^𝛼_𝑡 of a function 𝑓(𝑡) is defined as [3] ₀𝐷^𝛼_𝑡𝑓(𝑡)=1Γ(𝑛−𝛼)^𝑡₀𝑓^𝑛(𝜏)(𝑡−𝜏)^{𝛼−𝑛+1}().𝑑𝜏,𝑛−1<𝛼≤𝑛,𝑛∈𝑁 (11) The following are two basic properties of the Caputo fractional derivative: ₀𝐷^𝛼_𝑡𝑡^𝛽=Γ1+𝛽Γ𝑡1+𝛽−𝛼^𝛽−𝛼𝐽,0<𝛼<𝛽+1,𝛽>−1,_𝛼𝐷^𝛼𝑓(𝑡)=𝑓(𝑡)−^𝑛−1_𝑘=0𝑓^𝑘0⁺𝑡^𝑘,𝑘!𝑛−1<𝛼≤𝑛,𝑛∈𝑁. (12) For Caputo’s derivative we have 𝐷^𝛼𝐶=0,(𝐶isaconstant). (13) Similar to integer order differentiation, Caputo’s derivative is linear: 𝐷^𝛼𝛾𝑓(𝑡)+𝛿𝑔(𝑡)=𝛾𝐷^𝛼𝑓(𝑡)+𝛿𝐷^𝛼𝑔(𝑡), (14) where 𝛾 and 𝛿 are constants, and it satisfies the so-called Leibnitz’s rule 𝐷^𝛼(𝑔(𝑡)𝑓(𝑡))=∑^∞_𝑘=0(^𝛼_𝑘)𝑔^(𝑘)(𝑡)𝐷^𝛼−𝑘𝑓(𝑡), where 𝑓(𝜏) is continuous in [0,𝑡] and 𝑔(𝜏) has continuous derivatives sufficient number of times in [0,𝑡].

2.1.3. Definition of Grunwald-Letnikov Fractional Derivative

The Grunwald-Letnikov fractional derivative is defined as [3, 10] _𝑎𝐷^𝛼_𝑡𝑓(𝑡)=lim_{ℎ→0𝑛ℎ=𝑡−𝑎}ℎ^𝑛−𝑝_𝑟=0(−1)^𝑟𝑝𝑟𝑓(𝑡−𝑟ℎ), (15) where 𝜔^𝑝_𝑟=(−1)^𝑟𝑝_𝑟 and 𝜔^𝑝₀=1,𝜔^𝑝_𝑟=1−𝑝+1𝑟𝜔^𝑝_𝑟−1,𝑟=1,2,…. (16)

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 interval 𝑡∈[0,1) but in general case 𝑡∈[𝐴,𝐵], we divide the interval [𝐴,𝐵] into 𝑚 equal subintervals, each of width Δ𝑡=(𝐵−𝐴)/𝑚. In this case, the orthogonal set of Haar functions are defined in the interval [𝐴,𝐵] by [6] ℎ₀(𝑡)=[],ℎ1𝑡∈𝐴,𝐵0elsewhere,_𝑖(𝑡)=1,𝜁₁(𝑖)≤𝑡<𝜁₂(𝑖)−1,𝜁₂(𝑖)≤𝑡<𝜁₃(𝑖)0,otherwise, (17) where 𝜁₁(𝑖)=𝐴+𝑘−12^𝑗(𝐵−𝐴)=𝐴+𝑘−12^𝑗𝜁𝑚Δ𝑡,₂(𝑖)()=𝐴+𝑘−1/22^𝑗()()𝐵−𝐴=𝐴+𝑘−1/22^𝑗𝜁𝑚Δ𝑡,₃(𝑖)𝑘=𝐴+2^𝑗)𝑘(𝐵−𝐴=𝐴+2^𝑗𝑚Δ𝑡, (18) for 𝑖=1,2,…,𝑚, 𝑚=2^𝐽 and 𝐽 is a positive integer which is called the maximum level of resolution. Here 𝑗 and 𝑘 represent the integer decomposition of the index 𝑖; that is, 𝑖=𝑘+2^𝑗−1, 0≤𝑗<𝑖 and 1≤𝑘<2^𝑗+1.

3.1.1. Operational Matrix for General Order Integration [6]

The integration of the 𝐻_𝑚(𝑡)=[ℎ₀(𝑡),ℎ₁(𝑡),…,ℎ_𝑚−1(𝑡)]^𝑇 can be approximated by [16] ^𝑡₀𝐻_𝑚(𝜏)𝑑𝜏≅𝑄𝐻_𝑚(𝑡), (19) where 𝑄 is called the Haar wavelet operational matrix of integration which is a square matrix of 𝑚-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]: 𝐽^𝛼𝑓(𝑡)=1Γ(𝛼)^𝑡₀(𝑡−𝜏)^𝛼−1𝑓(𝜏)𝑑𝜏,𝑡>0,𝛼∈ℜ⁺, (20) where ℜ⁺ is the set of positive real numbers.

The Haar wavelet operational matrix 𝑄^𝛼 for integration of the general order 𝛼 is given by 𝑄^𝛼𝐻_𝑚(𝑡)=𝐽^𝛼𝐻_𝑚(𝑡)=𝐽^𝛼ℎ₀(𝑡),𝐽^𝛼ℎ₁(𝑡),…,𝐽^𝛼ℎ_𝑚−1(𝑡)^𝑇=𝑄ℎ₀(𝑡),𝑄ℎ₁(𝑡),…,𝑄ℎ_𝑚−1(𝑡)^𝑇, (21) where 𝑄ℎ₀(𝑡)=𝑡^𝛼Γ()[],1+𝛼,𝑡∈𝐴,𝐵0,elsewhere,𝑄ℎ_𝑖(𝑡)=0,𝐴≤𝑡<𝜁₁(𝑖),𝜙₁,𝜁₁(𝑖)≤𝑡<𝜁₂(𝑖),𝜙₂,𝜁₂(𝑖)≤𝑡<𝜁₃(𝑖),𝜙₃,𝜁₃(𝑖)≤𝑡<𝐵, (22) where 𝜙₁=𝑡−𝜁₁(𝑖)^𝛼Γ(),𝜙𝛼+1₂=𝑡−𝜁₁(𝑖)^𝛼Γ(𝛼+1)−2𝑡−𝜁₂(𝑖)^𝛼Γ(),𝜙𝛼+1₃=𝑡−𝜁₁(𝑖)^𝛼Γ(𝛼+1)−2𝑡−𝜁₂(𝑖)^𝛼Γ()+𝛼+1𝑡−𝜁₃(𝑖)^𝛼Γ(),𝛼+1 (23) for 𝑖=1,2,…,𝑚, 𝑚=2^𝐽 and 𝐽 is a positive integer, called the maximum level of resolution. Here 𝑗 and 𝑘 represent the integer decomposition of the index 𝑖. That is, 𝑖=𝑘+2^𝑗−1, 0≤𝑗<𝑖 and 1≤𝑘<2^𝑗+1.

3.1.2. Function Approximation

Any function 𝑦(𝑡)∈𝐿²([0,1)) can be expanded into Haar wavelets by [6, 10, 16] 𝑦(𝑡)=𝑐₀ℎ₀(𝑡)+𝑐₁ℎ₁(𝑡)+𝑐₂ℎ₂(𝑡)+⋯,where𝑐_𝑗=¹₀𝑦(𝑡)ℎ_𝑗(𝑡)𝑑𝑡. (24) If 𝑦(𝑡) is approximated as piecewise constant in each subinterval, the sum in (24) may be terminated after 𝑚 terms and consequently we can write discrete version in the matrix form as 𝐘≈^𝑚−1_𝑖=0𝑐_𝑖ℎ_𝑖𝑡_𝑙=𝐂^𝑇_𝑚𝐻_𝑚, (25) where 𝐘 and 𝐂^𝑇_𝑚 are the 𝑚-dimensional row vectors.

Here 𝐻 is the Haar wavelet matrix of order 𝑚 defined by 𝐻=[𝐡₀,𝐡₁,…,𝐡_𝑚−1]^𝑇; that is, 𝐡𝐻=⁰𝐡₁⋮𝐡_𝑚−1=ℎ_0,0ℎ_0,1⋯ℎ_0,𝑚−1ℎ_1,0ℎ_1,1⋯ℎ_1,𝑚−1⋮ℎ_𝑚−1,0ℎ_𝑚−1,1⋯ℎ_{𝑚−1,𝑚−1}, (26) where 𝐡₀,𝐡₁,…,𝐡_𝑚−1 are the discrete form of the Haar wavelet bases.

The collocation points are given by 𝑡_𝑙=𝐴+(𝑙−0.5)Δ𝑡,𝑙=1,2,…,𝑚. (27) Similarly, a function of two variables 𝑦(𝑥,𝑡)∈𝐿²([0,1]×[0,1]) can be approximated by discrete version in the matrix form of Haar wavelets as 𝑌(𝑥,𝑡)=𝐻^𝑇(𝑥)𝐶𝐻(𝑡), (28) where 𝑌(𝑥,𝑡) is the discrete form of 𝑦(𝑥,𝑡) and 𝐶 is the coefficient matrix of 𝑌, and it can be obtained by the following formula: 𝐶=𝐻^𝑇−1⋅𝑌⋅𝐻⁻¹. (29) Since 𝐻 is orthogonal, (29) becomes 𝐶=𝐻⋅𝑌⋅𝐻⁻¹. (30)

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) where 0≤𝑥<1, 0<𝑡≤1, 𝑎(𝑥)>0, 𝑏(𝑥)>0, 0<𝛼≤1, 1<𝛽≤2, 𝑞(𝑥,𝑡)∈𝐶(𝐷), and 𝐷=[0,1]×[0,1].

Since 𝑢(𝑥,𝑡)∈𝐿²(ℜ), we suppose 𝑢(𝑥,𝑡)≈^𝑚−1_𝑖=0^𝑚−1_𝑗=0𝑐_𝑖𝑗ℎ_𝑖(𝑥)ℎ_𝑗(𝑡). (32) Then we can obtain the discrete form of (32) by taking step Δ=1/𝑚 of 𝑥 and 𝑡: 𝑈=𝐻^𝑇(𝑥)𝐶𝐻(𝑡) (33) and then 𝜕^𝛼𝑢𝜕𝑡^𝛼≈𝜕^𝛼𝑈𝜕𝑡^𝛼=𝐻^𝑇(𝑥)𝐶𝜕^𝛼𝐻(𝑡)𝜕𝑡^𝛼=𝐻^𝑇(𝑥)𝐶𝑄^−𝛼𝐻(𝑡),𝜕^𝛽𝑢𝜕𝑥^𝛽≈𝜕^𝛽𝑈𝜕𝑥^𝛽=𝜕^𝛽𝜕𝑥^𝛽𝐻^𝑇(𝑥)(𝑡)=𝜕𝐶𝐻_𝛽𝜕𝑥^𝛽𝐻(𝑥)^𝑇𝐶𝐻(𝑡)=𝑄^−𝛽𝐻(𝑥)^𝑇𝐶𝐻(𝑡). (34) Here, 𝑞(𝑥,𝑡) is known function; discretizing it, we have 𝐷=𝑞𝑥_𝑖,𝑡_𝑗,𝑖,𝑗=0,1,2,…,𝑚−1. (35) Coefficients 𝑎(𝑥),𝑏(𝑥) can be dispersed into 𝑎(𝑥_𝑖),𝑏(𝑥_𝑖), 𝑖=0,1,2,…,𝑚−1.

Let 𝑎𝑥𝐴=⁰𝑥0⋯00𝑎₁𝑥⋯0⋮⋱⋮0⋯0𝑎_𝑚−1,𝑏𝑥𝐵=⁰𝑥0⋯00𝑏₁𝑥⋯0⋮⋱⋮0⋯0𝑏_𝑚−1. (36) Substituting (33)–(36) in (31), we have 𝐻^𝑇(𝑥)𝐶𝑄^−𝛼𝐻(𝑡)𝜕=−𝐵𝐻(𝑥)𝜕𝑥^𝑇(𝑡)𝜕𝐶𝐻+𝐴^𝛽𝜕𝑥^𝛽𝐻(𝑥)^𝑇𝐶𝐻(𝑡)+𝐷. (37) Equation (37) is a system of algebraic equations in 𝑎(𝑥_𝑖) and 𝑏(𝑥_𝑖). By solving this system of equations using mathematical software, the Haar wavelet coefficients 𝑎(𝑥_𝑖) and 𝑏(𝑥_𝑖) can be obtained.

3.2. Cubic B-Spline Basis Functions on 𝐻²(𝐼)

Let 𝐼=[0,𝐿] be an interval with 𝐿>4, and let 𝐻²(𝐼) be a Sobolev space which contains functions with square integrable second derivatives, and the homogenous Sobolev space 𝐻²₀(𝐼) can be defined by 𝐻²₀(𝐼)=𝑓(𝑡)∈𝐻²(𝐼),𝑓(0)=𝑓^′(0)=𝑓(𝐿)=𝑓^′(𝐿)=0, (38) which is a Hilbert space equipped with inner product: 𝑓,𝑔=_𝐼𝑓^′′(𝑡)𝑔^′′(𝑡)𝑑𝑡. (39) Cai and Wang [18] present a multiresolution analysis (MRA) and a wavelet decomposition for 𝐻²₀(𝐼) by constructing scaling spline functions 𝜑(𝑡)=16⁴_𝑙=04𝑙(−1)^𝑙(𝑡−𝑙)³₊,𝜑_𝑏(𝑡)=32𝑡²₊−11𝑡12³⁺+32(𝑡−1)³₊−34(𝑡−2)³₊+16(𝑡−3)³₊, (40) and wavelet functions 𝜓(𝑡)3=−7𝜑(2𝑡)+127𝜑()−32𝑡−17𝜑(),𝜓2𝑡−2_𝑏(𝑡)=24𝜑13^𝑏()−62𝑡𝜑(),132𝑡 (41) which have compact supports, where 𝑡^𝑛₊=𝑡^𝑛,𝑡>00,𝑡≤0. (42) Then the scaling spaces 𝑉₀⊂𝑉₁⊂𝑉₂⋯⊂𝑉_∞=𝐻²₀(𝐼) of MRA and a wavelet decomposition 𝐻²₀(𝐼)=𝑉₀⊕𝑊₀⊕𝑊₁⊕⋯⊕𝑊_𝑗−1⊕⋯ corresponding to the MRA are obtained by defining 𝑉_𝑗=span𝜑_𝑗,𝑘(𝑡),𝜑_𝑏𝑗(𝑡),𝜑_𝑏𝑗(𝐿−𝑡),0≤𝑘≤𝑛_𝑗,𝑊−4_𝑗=span𝜓_𝑗,𝑘(𝑡),𝑘=−1,0,…,𝑛_𝑗−2, (43) in which 𝑛_𝑗=2^𝑗𝐿, 𝜑_𝑗,𝑘(𝑡)=𝜑2^𝑗𝑡−𝑘,𝜑_𝑏𝑗(𝑡)=𝜑_𝑏2^𝑗𝑡,0≤𝑘≤𝑛_𝑗𝜓−4,_𝑗,𝑘(𝑡)=𝜓2^𝑗𝑡−𝑘,𝑗≥0,𝑘=0,1,…,𝑛_𝑗𝜓−3,_𝑗,−1(𝑡)=𝜓_𝑏2^𝑗𝑡,𝜓_{𝑗,𝑛𝑗−2}(𝑡)=𝜓_𝑏2^𝑗(𝐿−𝑡). (44) For convenience, we set 𝜓_−1,𝑘(𝑡)=𝜑_0,𝑘(𝑡), 0≤𝑘≤𝐿−4, 𝜓_−1,−1(𝑡)=𝜑_𝑏(𝑡), 𝜓_{−1,𝐿−3}(𝑡)=𝜑_𝑏(𝐿−𝑡), and 𝑛₋₁=𝐿−1. Let 𝐵_𝑗={2^−3𝑗/2𝜓_𝑗,𝑘(𝑡)}^𝑛𝑗−2_𝑘=−1, −1≤𝑗<∞. Wang [19] derived that 𝐵=⋃^∞_𝑗=−1𝐵_𝑗 is an unconditional basis of 𝐻²₀(𝐼), which turns out to be a basis of continuous space 𝐶₀(𝐼). For nonhomogeneous Sobolev space 𝐻²(𝐼), Cai and Wang [18] introduced boundary spline functions 𝜂₁(𝑡)=(1−𝑡)³₊,𝜂₂(𝑡)=2𝑡₊−3𝑡²₊+76𝑡³₊−43(𝑡−1)³₊+16(𝑡−2)³₊ (45) to deal with the values of functions at boundary points.

3.2.1. Function Approximation

Any function 𝑓(𝑡)∈𝐻²₀(𝐼) can be uniquely expanded into cubic spline wavelet series by 𝑓(𝑥)=^∞𝑛_𝑗=−1^𝑗_𝑘=1𝑑_𝑗,𝑘𝜓_{𝑗,𝑘−2}(𝑡), (46) with 𝑑_𝑗,𝑘=∫_𝐼𝑓^′′(𝑡)(𝜓^∗_{𝑗,𝑘−2})(𝑡)𝑑𝑡, where 𝜓^∗_𝑗,𝑘(𝑡) are dual functions of 𝜓_𝑗,𝑘. Truncating the infinite series in (46) at 𝐽−1, we get 𝑓_𝐽(𝑡)=^𝐽−1𝑛_𝑗=−1^𝑗_𝑘=1𝑑_𝑗,𝑘𝜓_{𝑗,𝑘−2}(𝑡). (47) From Wang [19] we have 𝑓(𝑡)−𝑓_𝐽(𝑡)²_𝐻20⟶0as𝐽⟶+∞. (48) Hence any function 𝑓(𝑡)∈𝐻²₀(𝐼) can be approximated by 𝑓_𝐽(𝑡) defined in (47). And any function 𝑓(𝑡)∈𝐻²(𝐼) can be approximated by 𝑓_𝐽(𝑡)=𝐼_𝑏,𝐽𝑓(𝑡)+^𝐽−1𝑛_𝑗=−1^𝑗_𝑘=1𝑑_𝑗,𝑘𝜓_{𝑗,𝑘−2}(𝑡), (49) and the approximation order is 𝑂(2^−4𝐽) if 𝑓(𝑡) is sufficiently smooth, where 𝐼_𝑏,𝐽𝑓(𝑡)=𝑎₁𝜂₁2^𝐽𝑡+𝑎₂𝜂₂2^𝐽𝑡+𝑎₃𝜂₂2^𝐽(𝐿−𝑡)+𝑎₄𝜂₁2^𝐽(𝐿−𝑡). (50) Suppose 𝑁=2^𝐽𝐿+3, and Ω_𝐽(𝑡) is a 1×𝑁 vector given by Ω_𝐽(𝑡)=𝜂₁2^𝐽𝑡,𝜂₂2^𝐽𝑡,𝜂₂2^𝐽(),𝜂𝐿−𝑡₁2^𝐽(𝐿−𝑡),𝜓_−1,−1(𝑡),𝜓_−1,0(𝑡)𝜓,…,_{−1,𝑛−1−2}(𝑡),𝜓_0,−1(𝑡),𝜓_0,0(𝑡),𝜓_0,1(𝑡),…,𝜓_0,𝐿−3(𝑡),𝜓_0,𝑛0−2(𝑡),𝜓_1,−1(𝑡),𝜓_1,0(𝑡),𝜓_1,1(𝑡)𝜓,…,_1,2𝐿−3(𝑡),𝜓_1,𝑛1−2(𝑡),…,𝜓_{𝐽−1,−1}(𝑡),𝜓_𝐽−1,0(𝑡),𝜓_𝐽−1,1(𝑡),…,𝜓_{𝐽−1,𝑛𝐽−3}(𝑡),𝜓_{𝐽−1,𝑛𝐽−2}(𝑡≡𝜔₁)(𝑡),𝜔₂(𝑡),…,𝜔_𝑁(𝑡). (51)𝑓_𝐽(𝑡) defined in (49) can be rewritten as 𝑓_𝐽(𝑡)=^𝑁_𝑘=1𝑓_𝑘𝜔_𝑘(𝑡)=Ω_𝐽(𝑡)𝑓, (52) where 𝑓=[𝑓₁,𝑓₂,…,𝑓_𝑁]^𝑇are the wavelet expansion coefficients which can be determined by interpolating at collocation points 𝑡⁽⁻¹⁾₋₁=12^𝐽+1,𝑡⁽⁻¹⁾_𝐿+11=𝐿−2^𝐽+1,𝑡_𝑘⁽⁻¹⁾𝑡=𝑘,𝑘=0,1,2,…,𝐿,_𝑘(𝑗)=𝑘+1.52^𝑗,−1≤𝑘≤2^𝑗𝐿−2,0≤𝑗≤𝐽−1. (53)

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

Consider the fractional differential equation of the form 𝐷^𝛼𝑦(𝑡)=−𝑦(𝑡)+𝑓(𝑡),0<𝑡≤𝐿, (54) with 𝑦^(𝑛)(0)=𝑦₀^(𝑛),𝑛=0,1,2,…,𝑚,𝑚<𝛼≤𝑚+1,𝑚∈𝑁. (55) To solve problem (54), we approximate 𝑦(𝑡) by 𝑦(𝑡)≈^𝑁_𝑘=1𝑐_𝑘𝜔_𝑘(𝑡)=𝑌_𝑗(𝑡), (56) where 𝑐_𝑘=[𝑐₁,𝑐₂,…,𝑐_𝑁]^𝑇 is unknown. The 𝛼 order fractional derivative of 𝑦(𝑡) is approximated by 𝐷^𝛼𝑦(𝑡)≈^𝑁_𝑘=1𝑐_𝑘𝐷^𝛼𝜔_𝑘(𝑡)≡𝐷^𝛼𝑌_𝑗(𝑡). (57) At collocation points, 𝑌_𝑗𝑡_𝑘=𝑌𝑡_𝑘𝐷,𝑘=1,2,…,𝑁,_𝛼𝑌_𝑗𝑡_𝑘=𝐷^𝛼𝑌𝑡_𝑘,𝑘=1,2,…,𝑁. (58) Interpolating the fractional differential equation (54) by 𝑌(𝑡_𝑘) at all collocation points, we obtain 𝐷^𝛼𝑌_𝑗𝑡_𝑘=−𝑌_𝑗𝑡_𝑘+𝑓𝑡_𝑘𝑦,2≤𝑘≤𝑁,_𝑗(𝑛)(0)=𝑦^(𝑛)(0),𝑛=0,1,2,…,𝑚. (59) Letting 𝐵₁=𝜔₁𝑡₂𝜔₂𝑡₂⋯𝜔_𝑁𝑡₂𝜔⋮⋱⋮₁𝑡_𝑁𝜔₂𝑡_𝑁⋯𝜔_𝑁𝑡_𝑁,𝐵₂=𝐷^𝛼𝜔₁𝑡₂𝐷^𝛼𝜔₂𝑡₂⋯𝐷^𝛼𝜔_𝑁𝑡₂𝐷⋮⋱⋮_𝛼𝜔₁𝑡_𝑁𝐷^𝛼𝜔₂𝑡_𝑁⋯𝐷^𝛼𝜔_𝑁𝑡_𝑁, (60) where 𝐵₁ and 𝐵₂ are obtained analytically.

Now (59) can be represented by 𝐴𝑐𝑦=𝑏,_𝑛𝑦=𝑦₀^(𝑛), (61) where 𝐴=𝐵₁+𝐵₂ and 𝑏=𝑓(𝑡₂),…,𝑓(𝑡_𝑁)^𝑇,𝑐_𝑛=𝜔₁^(𝑛)(0),𝜔₂^(𝑛)(0),…,𝜔_𝑁^(𝑛)(0),𝑛=0,1,2,…,𝑚. (62) Let 𝐴^𝑇_𝑚=(𝐴^𝑇,𝑐^𝑇₀,𝑐^𝑇₁,…,𝑐^𝑇_𝑚), 𝑏^𝑇_𝑚=(𝑏^𝑇,𝑦₀⁽⁰⁾,𝑦₀⁽¹⁾,…,𝑦₀^(𝑚)); then (61) can be written as 𝐴_𝑚𝑦=𝑏_𝑚. (63) Consequently the wavelet coefficient 𝑦 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 𝜓(𝑡) called the mother wavelet. When the dilation parameter 𝑎 and the translation parameter 𝑏 vary continuously, we have the following family of continuous wavelets as [13] 𝜓_𝑎,𝑏(𝑡)=|𝑎|^−1/2𝜓𝑡−𝑏𝑎,𝑎,𝑏∈ℜ,𝑎≠0. (64) If we restrict the parameters 𝑎 and 𝑏 to discrete values as 𝑎=𝑎₀^−𝑘, 𝑏=𝑛𝑏₀𝑎₀^−𝑘; 𝑎₀>1, 𝑏₀>0, where 𝑛 and 𝑘 are positive integer, the following family of discrete wavelets are defined as 𝜓_𝑘,𝑛(𝑡)=𝑎₀^𝑘/2𝜓𝑎^𝑘₀𝑡−𝑛𝑏₀, (65) where 𝜓_𝑘,𝑛(𝑡) forms a wavelet basis for 𝐿²(ℜ). In particular, when 𝑎₀=2 and 𝑏₀=1, then 𝜓_𝑘,𝑛(𝑡) forms an orthonormal basis; that is, 𝜓_𝑗,𝑘(𝑡),𝜓_𝑙,𝑚(𝑡)=𝛿_𝑗𝑙𝛿_𝑘𝑚. (66)

3.3.1. Chebyshev Wavelet and Function Approximation

Chebyshev wavelets 𝜓_𝑛,𝑚(𝑡), on the interval [0,1), are defined as [20] 𝜓_𝑛,𝑚(𝑡)=2^𝑘/2𝑇_𝑚2^𝑘𝑡−2𝑛+1,𝑛−12^𝑘−1𝑛≤𝑡<2^𝑘−10,otherwise, (67) where 𝑇_𝑚=1√𝜋,𝑚=0,2𝜋𝑇_𝑚(𝑡),𝑚>0, (68) and 𝑚=0,1,…,𝑀−1, 𝑛=1,2,…,2^𝑘−1,𝑘 is any positive integer, and 𝑇_𝑚(𝑡) are Chebyshev polynomials of the first kind of degree 𝑚 which are orthogonal with respect to the weight function 𝜔(𝑡)=(1/√1−𝑡²), on the interval [−1,1], and 𝑇_𝑚(𝑡) can be determined by the following recurrence formula: 𝑇₀(𝑡)=1,𝑇₁(𝑡)=𝑡,𝑠𝑇_𝑚+1(𝑡)=2𝑡𝑇_𝑚(𝑡)−𝑇_𝑚−1(𝑡),𝑚=1,2,3,…. (69) A function 𝑓(𝑡) defined over [0, 1) may be expanded as 𝑓(𝑡)=^∞∞_𝑛=1𝑚=0𝑐_𝑛,𝑚𝜓_𝑛,𝑚(𝑡), (70) where 𝑐_𝑛𝑚=⟨𝑓(𝑡),𝜓_𝑛,𝑚(𝑡)⟩.

If the infinite series in (70) is truncated, then (70) can be written as 𝑓(𝑡)≈^2𝑘−1_𝑛=1^𝑀−1_𝑚=0𝑐_𝑛,𝑚𝜓_𝑛,𝑚(𝑡)=𝐶^𝑇Ψ(𝑡), (71) where 𝐶 and Ψ(𝑡) are 2^𝑘−1𝑀×1 matrices given by 𝐶≡𝑐_1,0,𝑐_1,1,…,𝑐_1,𝑀−1,𝑐_2,0𝑐,…,_2,𝑀−1,…,𝑐₂^𝑘−1_,0,…,𝑐₂^𝑘−1_,𝑀−1^𝑇,Ψ(𝑡)≡𝜓_1,0,𝜓_1,1,…,𝜓_1,𝑀−1,𝜓_2,0𝜓,…,_2,𝑀−1,…,𝜓_2𝑘−1,0,…,𝜓_{2𝑘−1,𝑀−1}^𝑇. (72)(73) Taking the collocation points as follows: 𝑡_𝑖=2𝑖−12^𝑘𝑀,𝑖=1,2,…,2^𝑘−1𝑀, (74) we define the Chebyshev wavelet matrix Φ_𝑚×𝑚 as Φ_𝑚×𝑚≡Ψ1Ψ32𝑚2𝑚⋯Ψ2𝑚−1.2𝑚 (75)

3.3.2. Operational Matrix of the Fractional Integration

The integration of the vector Ψ(𝑡) defined in (73) can be obtained as ^𝑡₀Ψ(𝜏)𝑑𝜏≈𝑃Ψ(𝑡), (76) where 𝑃 is the 2^𝑘−1𝑀×2^𝑘−1𝑀 operational matrix for integration [13].

We know the Riemann-Liouville fractional integral operator 𝐽^𝛼, of a function 𝑓(𝑡), is defined as [3] 𝐽^𝛼𝑓(𝑡)=1Γ(𝛼)^𝑡₀(𝑡−𝜏)^𝛼−1𝑓(𝜏)1𝑑𝜏=Γ(𝛼)𝑡^𝛼−1∗𝑓(𝑡), (77) where 𝛼∈ℜ⁺(𝛼>0) is the order of integration, Γ(⋅) is the well-known gamma function, and 𝑡^𝛼−1∗𝑓(𝑡) denotes the convolution product of 𝑡^𝛼−1 and 𝑓(𝑡). Now if 𝑓(𝑡) is expanded in Chebyshev wavelets, the Riemann-Liouville fractional integration becomes 𝐽^𝛼𝑓(𝑡)=1Γ(𝛼)𝑡^𝛼−1∗𝑓(𝑡)≅𝐶^𝑇1Γ(𝛼)𝑡^𝛼−1∗Ψ(𝑡). (78) Thus if 𝑡^𝛼−1∗𝑓(𝑡) can be integrated and then expanded in Chebyshev wavelets, the Riemann-Liouville fractional integration can be solved via Chebyshev wavelets.

Also we define a 𝑚-set of Block-Pulse functions (BPF) as [13] 𝑏_𝑖(𝑡)=𝑖1,𝑚≤𝑡<𝑖+1𝑚0,otherwise, (79) where 𝑖=0,1,2,…,(𝑚−1).

The functions 𝑏_𝑖(𝑡) are disjoint and orthogonal; that is, 𝑏_𝑖(𝑡)𝑏_𝑙(𝑡)=𝑏0,𝑖≠𝑙_𝑖(𝑡),𝑖=𝑙,₁₀𝑏_𝑖(𝜏)𝑏_𝑙(𝜏)𝑑𝜏=10,𝑖≠𝑙𝑚,𝑖=𝑙. (80) Similarly, Chebyshev wavelets may be expanded into an 𝑚-term block pulse functions (BPF) as Ψ_𝑚′(𝑡)=Φ_{𝑚′×𝑚′}𝐵_𝑚′(𝑡), (81) where 𝑚^′=2^𝑘−1𝑀 and 𝐵_𝑚′(𝑡)≡𝑏₀(𝑡)𝑏₁(𝑡)⋯𝑏_𝑖(𝑡)⋯𝑏_𝑚−1(𝑡)^𝑇.

In [21], Kilicman and Al Zhour have proposed the block pulse operational matrix of the fractional integration 𝐹^𝛼 as follows: 𝐽^𝛼𝐵_𝑚′(𝑡)≈𝐹^𝛼𝐵_𝑚′(𝑡), (82) where 𝐹^𝛼=1𝑚^𝛼1Γ(𝛼+2)1𝜉₁𝜉₂⋅⋅⋅𝜉_𝑚−101𝜉₁⋅⋅⋅𝜉_𝑚−2001⋅⋅⋅𝜉_𝑚−3⋅⋅⋅⋅⋅⋅000⋅⋅00001 (83) with 𝜉_𝑘=(𝑘+1)^𝛼+1−2𝑘^𝛼+1+(𝑘−1)^𝛼+1.

Next we derive the Chebyshev wavelet operational matrix of the fractional integration. Let 𝐽^𝛼Ψ_𝑚(𝑡)≈𝑃^𝛼_{𝑚′×𝑚′}Ψ_𝑚(𝑡), (84) where matrix 𝑃^𝛼_{𝑚′×𝑚′} is called the Chebyshev wavelet operational matrix of the fractional integration.

Using (81) and (82), we have 𝐽^𝛼Ψ_𝑚′(𝑡)≈𝐽^𝛼Φ_{𝑚′×𝑚′}𝐵_𝑚′(𝑡)=Φ_{𝑚′×𝑚′}𝐽^𝛼𝐵_𝑚′(𝑡)≈Φ_{𝑚′×𝑚′}𝐹^𝛼𝐵_𝑚′(𝑡). (85) From (84) and (85) we get 𝑃^𝛼_{𝑚′×𝑚′}Ψ_𝑚′(𝑡)=Φ_{𝑚′×𝑚′}𝐹^𝛼𝐵_𝑚′(𝑡). (86) Then, the Chebyshev wavelet operational matrix of the fractional integration 𝑃^𝛼_{𝑚′×𝑚′} is given by 𝑃^𝛼_{𝑚′×𝑚′}=Φ_{𝑚′×𝑚′}𝐹^𝛼Φ_𝑚⁻¹_′×𝑚′. (87) 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 𝐷^𝛼𝑦(𝑡)+𝑦(𝑡)²=1,0<𝛼≤1, (88) subject to initial state 𝑦(0)=0.

Let 𝐷^𝛼𝑦(𝑡)=𝐶^𝑇Ψ(𝑡), (89) and then 𝑦(𝑡)=𝐶^𝑇𝑃^𝛼_{𝑚′×𝑚′}Ψ(𝑡). (90) Since Ψ(𝑡)=Φ_{𝑚′×𝑚′}𝐵_𝑚′(𝑡), we have 𝑦(𝑡)=𝐶^𝑇𝑃^𝛼_{𝑚′×𝑚′}Φ_{𝑚′×𝑚′}𝐵_𝑚′(𝑡). (91) Let 𝐶^𝑇𝑃^𝛼_{𝑚′×𝑚′}Φ_{𝑚′×𝑚′}=[𝑎₁,𝑎₂,…,𝑎_𝑚′], and, using (79), we get 𝑦(𝑡)²=𝑎₁𝑏₁(𝑡)+𝑎₂𝑏₂(𝑡)+⋯+𝑎_𝑚′𝑏_𝑚′(𝑡)²=𝑎²₁,𝑎²₂,…,𝑎²_𝑚′𝐵_𝑚′(𝑡). (92) Substituting (89) and (92) in (88), we have 𝐶^𝑇Φ_{𝑚′×𝑚′}𝐵_𝑚′(𝑡)+𝑎²₁,𝑎²₂,…,𝑎²_𝑚′𝐵_𝑚′(𝑡)−[1,1,…,1]𝐵_𝑚′(𝑡)=0. (93) Solving the above nonlinear system of algebraic equations, we can find the vector 𝐶 and consequently the solution for 𝑦(𝑡).

3.4. Legendre Wavelets

Legendre wavelet 𝜓_𝑛,𝑚=𝜓(𝑘,𝑛,𝑚,𝑡) have four arguments: 𝑛=1,…,2^𝑘−1,𝑘 can be assumed as any positive integer, 𝑚 is order of Legendre polynomials, and 𝑡 denotes the time.

The Legendre polynomials of order 𝑚, denoted by 𝐿_𝑚(𝑡), are defined on the interval [−1, 1] and can be determined with the help of the following recurrence formulae [11]: 𝐿₀(𝑡)𝐿=1,₁(𝑡)=𝑡,𝐿_𝑚+1(𝑡)=2𝑚+1𝑚+1𝑡𝐿_𝑚(𝑡)−𝑚𝐿𝑚+1^𝑚−1(𝑡),𝑚=1,2,3,…. (94)(95) The Legendre wavelets are defined on interval [0, 1) by 𝜓_𝑛,𝑚(𝑡)=1𝑚+2^1/22^𝑘/2𝐿_𝑚2^𝑘𝑡−𝑛,𝑛−12^𝑘≤𝑡<𝑛+12^𝑘,0,elsewhere (96) for 𝑘=2,3,…𝑛=2𝑛−1, 𝑛=1,2,3,…,2^𝑘−1, 𝑚=0,1,…,𝑀−1 is the order of the Legendre polynomials, and 𝑀 is a fixed positive integer. The set of Legendre wavelets forms an orthogonal basis of 𝐿²(ℜ).

3.4.1. Function Approximation

A function 𝑓(𝑡)∈𝐿²[0,1] can be expanded in terms of Legendre wavelet as [23] 𝑓(𝑡)=^∞∞_𝑛=1𝑚=0𝑐_𝑛,𝑚𝜓_𝑛,𝑚(𝑡), (97) where 𝑐_𝑛,𝑚=⟨𝑓(𝑡),𝜓_𝑛,𝑚(𝑡)⟩.

If the infinite series in (97) is truncated, then (97) can be written as 𝑓(𝑡)≈^2𝑘−1_𝑛=1^𝑀−1_𝑚=0𝑐_𝑛,𝑚𝜓_𝑛,𝑚(𝑡)=𝐶^𝑇Ψ(𝑡), (98) where 𝐶 and Ψ(𝑥) are 2^𝑘−1𝑀×1 matrices given by 𝐶≡𝑐_1,0,𝑐_1,1,…,𝑐_1,𝑀−1,𝑐_2,0,…,𝑐_2,𝑀−1𝑐,…,_2𝑘−1,0,…,𝑐₂^𝑘−1_,𝑀−1^𝑇,Ψ(𝑡)≡𝜓_1,0,𝜓_1,1,…,𝜓_1,𝑀−1,𝜓_2,0,…,𝜓_2,𝑀−1𝜓,…,_2𝑘−1,0,…,𝜓_{2𝑘−1,𝑀−1}^𝑇. (99) Similarly, a function of two variables 𝑘(𝑥,𝑡)∈𝐿²([0,1]×[0,1]) can be approximated by Legendre wavelet series as 𝑘(𝑥,𝑡)≈Ψ^𝑇(𝑡)𝐾Ψ(𝑥), (100) where 𝐾=[𝑘_𝑖𝑗] is (2^𝑘−1𝑀×2^𝑘−1𝑀) matrix, with 𝑘_𝑖,𝑗=𝜓_𝑖(𝑡),𝑘(𝑥,𝑡),𝜓_𝑗(𝑥). (101)

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]: 𝐷^𝛼𝑦(𝑡)+𝑦(𝑡)=0,0<𝛼≤2, (102) such that 𝑦(0)=1,𝑦^′(0)=0. (103) The corresponding integral representation for (102) and (103) is 𝑦=−𝐽^𝛼𝑦(𝑡)+𝑦(0)+𝑦^′(0)𝑡. (104) Let 𝑦 be approximated using Legendre wavelet as 𝑦=𝐶^𝑇Ψ(𝑡), (105) and then 𝐽^𝛼𝑦=𝐶^𝑇𝐽^𝛼Ψ(𝑡)=𝐶^𝑇𝑃^𝛼_𝑚𝑚×Ψ(𝑡),where𝑚=2^𝑘−1𝑀. (106) Substituting (105) and (106) in (104), we have the following system of algebraic equations: 𝐶^𝑇Ψ(𝑡)+𝐶^𝑇𝑃^𝛼_𝑚𝑚×Ψ(𝑡)=𝑌₀,where𝑌₀=𝑦(0)+𝑦^′(0)𝑡. (107) By solving the linear system (107), we can find the vector 𝐶 and hence consequently the solution 𝑦(𝑡).

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 𝜓_𝑛,𝑚(𝑡)=𝜓(𝑘,𝑛,𝑚,𝑡) have four arguments: 𝑛=0,1,…,2^𝑘−1, 𝑘 can be assumed as any positive integer, 𝑚 is order of Legendre polynomials, and 𝑡 denotes the normalized time. The Legendre multiwavelets are defined on the interval [0, 1) by 𝜓_𝑛,𝑚(𝑡)=(2𝑚+1)^1/22^𝑘/2𝑃_𝑚2^𝑘,𝑛𝑡−𝑛2^𝑘≤𝑡<𝑛+12^𝑘,0,elsewhere, (108) where 𝑚=0,1,…,𝑀−1 and 𝑛=1,2,3,…,2^𝑘−1. The coefficient (2𝑚+1)^1/2 is for orthonormality and 𝑃_𝑚(𝑡) are the well-known shifted Legendre polynomials of order 𝑚 which are defined on the interval [0,1] and can be determined with the help of the following recurrence formula: 𝑃₀(𝑡)𝑃=1,₁(𝑡)=2𝑡−1,𝑃_𝑚+1(𝑡)=2𝑚+1()𝑃𝑚+12𝑡−1_𝑚(𝑡)−𝑚𝑃𝑚+1^𝑚−1(𝑡),𝑚=1,2,3,…. (109)(110) The two-dimensional Legendre multiwavelets are defined as 𝜓_{𝑛1,𝑚1,𝑛2,𝑚2}()=𝑥,𝑡𝐴𝑃_𝑚12^𝑘1𝑥−𝑛₁𝑃_𝑚22^𝑘2𝑡−𝑛₂,𝑛₁2^𝑘1𝑛≤𝑥<¹+12^𝑘1,𝑛₂2^𝑘2𝑛≤𝑡<²+12^𝑘20,elsewhere, (111) where 𝐴=2𝑚₁+12𝑚₂+12^{(𝑘1+𝑘2)/2}. (112)𝑛₁ and 𝑛₂ are defined similarly to 𝑛, 𝑘₁ and 𝑘₂ are any positive integer, 𝑚₁ and 𝑚₂ are the order for Legendre polynomials, and 𝜓_{𝑛1,𝑚1,𝑛2,𝑚2}(𝑥,𝑡) forms a basis for 𝐿²([0,1]×[0,1]).

3.5.1. Function Approximation

A function 𝑓(𝑥,𝑡) defined over [0,1]×[0,1] can be expanded in terms of Legendre multiwavelet as [12] 𝑓(𝑥,𝑡)=^∞∞_𝑛=0^∞_𝑖=0^∞_𝑙=0𝑗=0𝑐_{𝑛,𝑖,𝑙,𝑗}𝜓_𝑛,𝑖(𝑥)𝜓_𝑙,𝑗(𝑡). (113) If the infinite series in (113) is truncated, then (113) can be written as 𝑓(𝑥,𝑡)=^{2𝑘1−1𝑁}_𝑛=0²_𝑖=0^{𝑘2−1𝑀}_𝑙=0𝑗=0𝑐_{𝑛,𝑖,𝑙,𝑗}𝜓_𝑛,𝑖(𝑥)𝜓_𝑙,𝑗(𝑡)=Ψ^𝑇(𝑥)𝐶Ψ(𝑡), (114) where Ψ(𝑥) and Ψ(𝑡) are 2^𝑘1(𝑀₁+1)×1 and 2^𝑘2(𝑀₂+1)×1 matrices, respectively, given by Ψ(𝑥)≡𝜓_0,0(𝑥),𝜓_0,1(𝑥),…,𝜓_0,𝑀1(𝑥),𝜓_1,0(𝑥)𝜓,…,_1,𝑀1(𝑥),…,𝜓_2,0(𝑥),…,𝜓_2,𝑀1(𝑥)𝜓,…,_2𝑘1−1,0(𝑥),…,𝜓_{2𝑘1−1,𝑀1}(𝑥,Ψ(𝑡)≡𝜓_0,0)(𝑡),𝜓_0,1(𝑡),…,𝜓_0,𝑀2(𝑡),𝜓_1,0(𝑡)𝜓,…,_1,𝑀2(𝑡),…,𝜓_2,0(𝑡),…,𝜓_2,𝑀2(𝑡)𝜓,…,_2𝑘2−1,0(𝑡),…,𝜓_{2𝑘2−1,𝑀2}(𝑡)^𝑇. (115) Also, 𝐶 is a 2^𝑘1(𝑀₁+1)×2^𝑘2(𝑀₂+1) matrix whose elements can be calculated from the following formula: 𝐶=¹₀¹₀𝜓_𝑛,𝑖(𝑥)𝜓_𝑙,𝑗(𝑡)𝑓(𝑥,𝑡)𝑑𝑡𝑑𝑥 (116) with 𝑛=0,1,…,2^𝑘1−1, 𝑖=0,…,𝑀₁, 𝑙=0,1,…,2^𝑘2−1, and 𝑗=0,…,𝑀₂.

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

Consider the nonlinear fractional partial differential equation 𝐷^𝛼_𝑡𝑢=𝑁(𝑢)+𝑔(𝑥,𝑡),𝑚<𝛼<𝑚+1,𝑚≥0, (117) with initial condition 𝑢(𝑥,0)=𝑓(𝑥).

Let 𝐹(𝑢)=𝐷^𝛼_𝑡𝑢−𝑁(𝑢)−𝑔(𝑥,𝑡). (118) A Galerkin approximation to (118) is constructed as follows. The approximation 𝑢_𝑁𝑀 are sought in the form of the truncated series 𝑢_𝑁𝑀()=𝑥,𝑡_{2𝑘1−1𝑁𝑛=02𝑖=0𝑘2−1𝑀𝑙=0𝑗=0}𝑡𝑐_{𝑛,𝑖,𝑙,𝑗}×𝜓_𝑛,𝑖(𝑥)𝜓_𝑙,𝑗(𝑡)+𝑢(𝑥,0),for𝑚=0_{2𝑘1−1𝑁𝑛=02𝑖=0𝑘2−1𝑀𝑙=0𝑗=0}𝑡²𝑐_{𝑛,𝑖,𝑙,𝑗}×𝜓_𝑛,𝑖(𝑥)𝜓_𝑙,𝑗(𝑡)+𝑢(𝑥,0)+𝑡𝑢_𝑡(𝑥,0),for𝑚=1, (119) where 𝜓_𝑖𝑗 are Legendre multiwavelet basis.

The expansion coefficients 𝑐_{𝑛,𝑖,𝑙,𝑗} are determined by Galerkin equations 𝐹𝑢_𝑁𝑀,𝜓_𝑛,𝑖𝜓_𝑙,𝑗=0, (120) where 𝐹𝑢_𝑁𝑀,𝜓_𝑛,𝑖𝜓_𝑙,𝑗=¹₀¹₀𝐹𝑢_𝑁𝑀(𝑥,𝑡)𝜓_𝑛,𝑖(𝑥)𝜓_𝑙,𝑗(𝑡)𝑑𝑡𝑑𝑥. (121) Galerkin equations (120) gives a system of 2^𝑘1(𝑁+1)×2^𝑘2(𝑀+1) equations which can be solved for the elements of 𝑐_{𝑛,𝑖,𝑙,𝑗}, 𝑖=0,1,2,…,𝑁, 𝑗=0,1,2,…,𝑀,  𝑛=0,1,2,…,2^𝑘1−1, and 𝑙=0,1,2,…,2^𝑘2−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 𝑘 and 𝑀. 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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This research work was financially supported by SERB, DST, Government of India, under Grant no. SR/S4/MS.:722/11.