# A New Scheme for Solving Multiorder Fractional Differential Equations Based on Müntz–Legendre Wavelets

Authors: Haifa Bin Jebreen; Fairouz Tchier
Journal: Complexity
 (2021)

---

## Abstract

In this study, we apply the pseudospectral method based on Müntz–Legendre wavelets to solve the multiorder fractional differential equations with Caputo fractional derivative. Using the operational matrix for the Caputo derivative operator and applying the Chebyshev and Legendre zeros, the problem is reduced to a system of linear algebraic equations. We illustrate the reliability, efficiency, and accuracy of the method by some numerical examples. We also compare the proposed method with others and show that the proposed method gives better results.


---

## Body

## 1. Introduction

This paper is dedicated to the numerical solution of the multiorder fractional differential equation with Caputo fractional derivative based on Müntz–Legendre wavelets. Letα∈ℝ+, ℕ∍n≔−α, for α∉ℕ. Further, let αj∈ℝ+j=1,…,σ∈ℕ which satisfies(1)0=α0<α1<⋯<ασ<α.We aim to compute the approximate solution of the linear or nonlinear fractional differential equation(2)Dc0αyx=fx,yx,Dc0α1yx,…,Dc0ασyx,x∈0,1,yη0=qη,qη∈ℝ,η=0,1,…,n−1,where Dc0α is the Caputo fractional differential operator [1] and the function fx,y,y1,…,yσ∈C0,1 with y1,…,yσ∈G (G⊂ℂ is an open set) satisfies the Lipschitz condition(3)fx,y,y1,…,yσx−fx,u,u1,…,uσx≤ℒσ∑j=0σyj−uj,yj,uj∈G,where ℒσ>0 is independent of x. We also assume for simplicity that(4)fα0,…,ασt,y≔fx,yx,Dc0α1yx,…,Dc0ασyx.The existence and uniqueness ofL1 solution of the multiorder fractional differential equation with the Riemann–Liouville and Caputo fractional derivatives under the assumption that fx,t∈L10,1 satisfies the Lipschitz condition with respect to the second variable are investigated by Kilbas et al. [1]. The previous investigation is based on converting equation (2) into the equivalent Volterra integral equation and then solving it. But in this study, we solve the problem directly.There exist some numerical methods that solve the desired problem. Laksetani et al. [2] introduced an operational method using B-spline functions to solve the multiorder fractional differential equation(5)Fux,Dc0β1,…,Dc0βm=gx,βi∈ℝ.In this paper, the operational matrix of the Caputo fractional derivative has been constructed directly, and then using the collocation method, the problem is solved. In [3], the authors applied the collocation method to solve the fractional differential equation(6)Dcαyt=ft,yt,α∈0,1.Their investigation is based on the collocation method using the Chebyshev–Gauss–Lobatto collocation points. The main advantage of this method is its superior accuracy. We can point out its exponential convergence too. But solving this problem is not a big challenge because it lacks multiorder fractional derivatives. Dehestani et al. [4] applied the fractional-Lucas optimization method to solve the multidimensional and multiorder fractional differential equation with Caputo fractional derivative. To this end, they used the operational matrix of fractional derivative for Lucas functions and reduced the problem into a linear or nonlinear system. The result shows the accuracy and efficiency of the method. A special type of equation (2) is considered by Bhrawy et al. [5] as(7)Dcαyt+γyt=ft.To solve this equation, they utilized the Laguerre tau technique. To this end, firstly fractional-order generalized Laguerre functions are introduced and Caputo fractional-order derivative is represented by these bases. In [6], after introducing the fractional-order Legendre functions and the operational matrix of Caputo derivative, the multiorder fractional differential equation is solved. For more details, we refer the readers to [7–10].The fractional differential equations are applied to model various physical phenomena, such as heat conduction, viscoelasticity, dynamical behavior of quantum particles, and laxation and diffusion problems [11–16].This paper is organized as follows. In Section2, we introduce the Müntz–Legendre wavelets, and we construct the operational matrix of fractional integration and Caputo fractional derivative. In Section 3, the pseudospectral method is applied to solve the generalized Cauchy-type problems with Caputo fractional differentiation based on Müntz–Legendre wavelets, and then the error analysis is investigated. Section 4 is devoted to some numerical examples to illustrate the accuracy and efficiency of the proposed method.


## 2. Müntz–Legendre Wavelets

In the last decade, wavelets have been able to get a special place in numerical analysis and especially in the numerical solution of equations [17–22]. As you know, one of the ways to get wavelets is to use multiresolution analysis (MRA). MRA is a family of nested spaces that satisfies certain circumstances [23], namely,(8)0⊂⋯⊂V−1⊂V0⊂V1⊂⋯⊂L2Ω,where Ω is a bounded interval or is equal to ℝ.In this paper, we apply Müntz–Legendre wavelets to solve the multiorder fractional differential equations. To this end, we give a brief introduction to Müntz–Legendre wavelets. Assume that the spaceVJ (J∈ℤ+∪0) is spanned by a set of bases which are called multiscaling functions or mother wavelets, i.e.,(9)VJ=spanψJ,bm:b∈ℬ,m∈ℳ,where ℬ≔0,1,…,2J−1 and ℳ≔0,1,…,r−1, r∈ℕ. The parameter r is called the multiplicity parameter, and J is the refinement level. In the following, we introduce the functions ψJ,bm.Assume thatλk≔kμ where μ is a real constant. Denote by Lmx the Müntz–Legendre polynomials [24] which are defined on Ω≔0,1 as(10)Lmx=∑k=0mlk,mxλk,x∈Ω,where the coefficient lk,m is defined by [24](11)lk,m≔∏i=0m−1λk+λi+1∏i=0,i≠kmλk−λi.Among the properties of these functions, we can mention their orthogonality. These polynomials form an orthogonal system that satisfies the following relation:(12)Lm′x,Lmx=∫01Lm′xLmxdx=δm′,m2λm+1,m≥m′,where δm′,m is used for the Kronecker symbol and is given by(13)δm′,m≔1,m=m′,0,m≥m′.Now we are ready to introduce the functionsψJ,bm. The Müntz–Legendre wavelets on the interval 0,1 are defined as follows [24]:(14)ψJ,bm=2J/22λm+1Lm2Jx−b,b2J≤x≤b+12J,0,otherwise.LetPJ be an operator that projects any function fx∈L20,1 onto the subspace VJ as follows:(15)fx≈PJfx=∑b=02J−1∑m=0r−1fb,mψJ,bmx,where the coefficients fb,m are evaluated by(16)fb,m=f,ψJ,bm=∫01fxψJ,bmxdx.LetF and Ψx be vectors of dimension N=2Jr whose br+m+1-th element is fb,m and ψJ,bmx, respectively. Hence, it follows from equation (15) that(17)fx≈PJfx=FTΨx,where the superscript T is used for the matrix transpose.It follows from [24] that there are some error estimates in the sense of Sobolev norms.Lemma 1 (see [24]).
Letn≥0 and r>n. If f∈Hn0,1, then(18)f−PJfL20,1≤cr−1−n2J−1−nfnL20,1,and for s≥1, we have(19)f−PJfHs0,1≤cr−12s−1/2−n2J−1s−nfnL20,1,where Hn0,1 is the Sobolev space and(20)fHn0,1=∑j=0nfjL20,121/2.


### 2.1. Representation of the Caputo Fractional Derivative Operator in Müntz–Legendre Wavelets

Recall that the Riemann–Liouville fractional integral operatorℐ0α (α∈ℝ+) is determined by(21)ℐaαfx≔1Γα∫axx−sα−1fsds,x∈a,b.Note that iff∈L1a,b, then the function ℐaαf∈L1a,b. We know that there is an operator that satisfies the relation(22)DRLaα=Dnℐan−α,where D≔d/dx and DRLaα is called the Riemann–Liouville fractional derivative operator. There is also another fractional derivative operator that satisfies the relation Dcaαfx≔ℐan−αDnfx and is called the Caputo fractional derivative.In this section, we would like to represent the Caputo fractional derivative operator in Müntz–Legendre wavelets. To this end, we first construct the operational matrix for fractional integral operatorℐ0α. Then, applying the operational matrix of derivative D for Müntz–Legendre wavelets [25] and relation Dcaαfx≔ℐan−αDnfx, we can find the operational matrix of fractional derivative for Müntz–Legendre wavelets.Applying the projection operatorP, the fractional integral operator ℐ0α acting on the vector function Ψx can be approximated by(23)Pℐ0αΨx=IαxΨx,α∈0,1,where Iαx is the operational matrix of integral for the Müntz–Legendre wavelets.To facilitate the evaluation of the operational matrix elements of fractional integration for the Müntz–Legendre wavelets, it is necessary to introduce the piecewise fractional-order Taylor functions. For a fixedJ∈ℤ+∪0, these functions are constructed as(24)ϕJ,bm=tλm,b2J≤x≤b+12J,0,otherwise,b∈ℬ,m∈ℳ.LetΦx be a vector of dimension N whose br+m+1-th element is ϕJ,bmx.To derive matrixIαx, we first introduce a matrix of dimension N×N that is used to transfer the Müntz–Legendre wavelets Ψx to the piecewise fractional-order Taylor functions Φx. Assume that there is a matrix of dimension N×N such that(25)Ψx=T−1Φx,where TJ−1 stands for the inverse of the matrix T. The matrix T is called the transformation matrix whose i,j-th element is evaluated by(26)Ti,j=Φix,Ψjx=∫01ΦixΨjxdx,i,j=1,…,N.LetΛ be a vector of dimension r whose i-th element is xλi, and thus it follows from equation (24) that(27)Φx=Λ,…,ΛT.It is easy to verify that the Riemann–Liouville fractional integration of the power functionsxκ is equal to power functions of the same form, i.e.,(28)ℐ0αxκ=Γκ+1Γκ+α+1xκ+α.This gives rise to find thei-th element of ℐ0αΦx, via(29)ℐ0αΦix=Γλi+1Γλi+α+1xλi+α.Thus, there exists a diagonal matrixIΦ,αx of dimension N×N such that(30)ℐ0αΦx=IΦ,αxΦx.The matrixIΦ,αx elements are obtain as follows:(31)IΦ,αx=diagIT,αx,…,IT,αx,where IT,αx≔xαG (ℐ0αΛx=IT,αxΛx) and G is a diagonal matrix of the form(32)Gi,i=Γλi+1Γλi+α+1.To derive the operational matrix of integral for the Müntz–Legendre wavelets, it follows from equation (25) that(33)Pℐ0αΨx=Pℐ0αT−1Φx=T−1IΦ,αxΦx=T−1IΦ,αxTΨx.This gives rise to(34)Iαx≔T−1IΦ,αxT.Now usingDcaα≔ℐan−αDn, we can introduce the operational matrix for the Caputo fractional derivative(35)Dαx≔DnIn−αx.


## 2.1. Representation of the Caputo Fractional Derivative Operator in Müntz–Legendre Wavelets

Recall that the Riemann–Liouville fractional integral operatorℐ0α (α∈ℝ+) is determined by(21)ℐaαfx≔1Γα∫axx−sα−1fsds,x∈a,b.Note that iff∈L1a,b, then the function ℐaαf∈L1a,b. We know that there is an operator that satisfies the relation(22)DRLaα=Dnℐan−α,where D≔d/dx and DRLaα is called the Riemann–Liouville fractional derivative operator. There is also another fractional derivative operator that satisfies the relation Dcaαfx≔ℐan−αDnfx and is called the Caputo fractional derivative.In this section, we would like to represent the Caputo fractional derivative operator in Müntz–Legendre wavelets. To this end, we first construct the operational matrix for fractional integral operatorℐ0α. Then, applying the operational matrix of derivative D for Müntz–Legendre wavelets [25] and relation Dcaαfx≔ℐan−αDnfx, we can find the operational matrix of fractional derivative for Müntz–Legendre wavelets.Applying the projection operatorP, the fractional integral operator ℐ0α acting on the vector function Ψx can be approximated by(23)Pℐ0αΨx=IαxΨx,α∈0,1,where Iαx is the operational matrix of integral for the Müntz–Legendre wavelets.To facilitate the evaluation of the operational matrix elements of fractional integration for the Müntz–Legendre wavelets, it is necessary to introduce the piecewise fractional-order Taylor functions. For a fixedJ∈ℤ+∪0, these functions are constructed as(24)ϕJ,bm=tλm,b2J≤x≤b+12J,0,otherwise,b∈ℬ,m∈ℳ.LetΦx be a vector of dimension N whose br+m+1-th element is ϕJ,bmx.To derive matrixIαx, we first introduce a matrix of dimension N×N that is used to transfer the Müntz–Legendre wavelets Ψx to the piecewise fractional-order Taylor functions Φx. Assume that there is a matrix of dimension N×N such that(25)Ψx=T−1Φx,where TJ−1 stands for the inverse of the matrix T. The matrix T is called the transformation matrix whose i,j-th element is evaluated by(26)Ti,j=Φix,Ψjx=∫01ΦixΨjxdx,i,j=1,…,N.LetΛ be a vector of dimension r whose i-th element is xλi, and thus it follows from equation (24) that(27)Φx=Λ,…,ΛT.It is easy to verify that the Riemann–Liouville fractional integration of the power functionsxκ is equal to power functions of the same form, i.e.,(28)ℐ0αxκ=Γκ+1Γκ+α+1xκ+α.This gives rise to find thei-th element of ℐ0αΦx, via(29)ℐ0αΦix=Γλi+1Γλi+α+1xλi+α.Thus, there exists a diagonal matrixIΦ,αx of dimension N×N such that(30)ℐ0αΦx=IΦ,αxΦx.The matrixIΦ,αx elements are obtain as follows:(31)IΦ,αx=diagIT,αx,…,IT,αx,where IT,αx≔xαG (ℐ0αΛx=IT,αxΛx) and G is a diagonal matrix of the form(32)Gi,i=Γλi+1Γλi+α+1.To derive the operational matrix of integral for the Müntz–Legendre wavelets, it follows from equation (25) that(33)Pℐ0αΨx=Pℐ0αT−1Φx=T−1IΦ,αxΦx=T−1IΦ,αxTΨx.This gives rise to(34)Iαx≔T−1IΦ,αxT.Now usingDcaα≔ℐan−αDn, we can introduce the operational matrix for the Caputo fractional derivative(35)Dαx≔DnIn−αx.


## 3. Pseudospectral Method

To derive the numerical solution of equation (2), the approximate solution can be approximated by Müntz–Legendre wavelets as follows:(36)y≈Pyx=YTΨx,where Y is a vector of dimension N that should be determined. A similar expression is valid for Fα0,…,ασx,Pyx∈Cγ0,1 where Cγ0,1 is the space of functions that satisfies x−aγfx∈C0,1 for 0≤γ<1.(37)Pfα0,…,ασx,Pyxx=FTΨx.To compute the elements of matrixF, let yJ=Py; then, using Dcaα≔ℐan−αDn,(38)Pfα0,…,ασx,yJxx=fx,yJx,Dcaα1yJx,…,DcaασyJx=fx,yJx,ℐ0n1−α1Dn1yJx,…,ℐ0nσ−ασDnσyJx,ℐ0nj−αjDnjyJx≔YTℐ0nj−αjDnjΨx=YTℐ0nj−αjDnjΨx=YTDnjInj−αjxΨx,j=1,…,σ.It follows from equations (36) and (38) that we can compute the residual in approximating equation (2) as follows:(39)rJx=DαYT−FTxΨx.We aim to reduce the residual to zero. One of the available methods is to use the pseudospectral method such thatrJxi=0 where xi are the collocation points. In this paper, we use the shifted Legendre and Chebyshev polynomial zeros. This gives rise to a system of linear or nonlinear algebraic equations that should be solved to find the unknown coefficients Y. To apply the initial conditions (2), we replace the first n equations of the obtained system of the pseudospectral method with them.(40)rJxiη≔YTDηΨ0−qη,η=0,…,n−1.


### 3.1. Convergence Analysis

It follows from [1] that the fractional integration operators ℐ0α is bounded in Lp0,1 (see [1] Lemma 1(a))(41)ℐ0αyLp0,1≤CαyLp0,1,Cα=1Γα+1.Also there is an optimal error estimate in term of error between the Müntz–Legendre polynomials derivativeDy and the exact derivative PDy, via(42)Dy−PDyL20,1≤CDr−11−nyHn,r−10,1,where CD is a constant and(43)gHn,r−10,12=∑l=minn,r−1r−1gltL20,12is a seminorm.Theorem 1.
Letn=−α and α∈ℝ+. Assume that fα0,…,ασt,y∈Hn0,1 and satisfies the Lipschitz condition (3). Also, assume that y and yJ are the exact and the approximate solutions (39) of equation (2), respectively.
Ify is a sufficiently smooth function, then the overall error(44)ex=Dc0αyx−PDc0αyx−fα0,…,ασt,y+Pfα0,…,ασt,yJsatisfies(45)exL20,1⟶0,asr⟶∞,orJ⟶∞.Proof.
Iffα0,…,ασt,y∈Hn0,1, then by Lemma 1, we can write(46)fα0,…,ασt,yJ−PJfα0,…,ασt,yJL20,1≤c1r−1−n2J−1−nfα0,…,ασnt,yJL20,1.
Using the Lipschitz conditions (3) and (41), we can write the following bound:(47)fα0,…,ασt,y−fα0,…,ασt,yJL20,1≤ℒσ∑j=0σDc0αjy−yJL20,1=ℒσ∑j=0σℐ0nj−αjDnjy−yJL20,1≤ℒσ∑j=0σCnj−αjDnjy−yJL20,1.
Now, applying equation (41) and Theorem 2 [24], we have(48)fα0,…,ασt,y−fα0,…,ασt,yJL20,1≤ℒσ∑j=0σCnj−αjCDnjr−11−nynj−1Hn,r−10,1≤ℒσ∑j=0σCnj−αjCDnjr−11−n2J−1−nyn+nj−1L20,1.
Also, using Lemma1 and equation (41), we have(49)Dc0αyx−PDc0αyxL20,1≤c2r−1−n2J−1−nDc0αyL20,1=≤c2r−1−n2J−1−nℐ0n−αDnyL20,1≤Cn−αc2r−1−n2J−1−nynL20,1.
Subtracting equation (2) from equation (39), one can write(50)ex=Dc0αyx−PDc0αyx−fα0,…,ασt,y+Pfα0,…,ασt,yJ+fα0,…,ασt,yJ−fα0,…,ασt,yJ.
Taking norms from both sides of equation (50) and using equations (46)–(49), we end up with(51)exL20,1≤Dc0αyx−PDc0αyxL20,1+fα0,…,ασt,y−fα0,…,ασt,yJL20,1+fα0,…,ασt,yJ−Pfα0,…,ασt,yJL20,1≤Cn−αc2r−1−n2J−1−nynL20,1+ℒσ∑j=0σCnj−αjCDnjr−11−n2J−1−nyn+nj−1L20,1+c1r−1−n2J−1−nfα0,…,ασnt,yJL20,1.
Suppose thatC1≔maxCnj−αjCDnj (j=0,…,σ) and C≔maxCn−αc2,C1,c1; then, we can obtain the following from equation (51):(52)exL20,1≤Cr−1−n2J−1−nynL20,1+ℒσ∑j=0σyn+nj−1L20,1+fα0,…,ασnt,yJL20,1.
It is easy to show thatexL20,1⟶0 as r⟶∞ or J⟶∞.


## 3.1. Convergence Analysis

It follows from [1] that the fractional integration operators ℐ0α is bounded in Lp0,1 (see [1] Lemma 1(a))(41)ℐ0αyLp0,1≤CαyLp0,1,Cα=1Γα+1.Also there is an optimal error estimate in term of error between the Müntz–Legendre polynomials derivativeDy and the exact derivative PDy, via(42)Dy−PDyL20,1≤CDr−11−nyHn,r−10,1,where CD is a constant and(43)gHn,r−10,12=∑l=minn,r−1r−1gltL20,12is a seminorm.Theorem 1.
Letn=−α and α∈ℝ+. Assume that fα0,…,ασt,y∈Hn0,1 and satisfies the Lipschitz condition (3). Also, assume that y and yJ are the exact and the approximate solutions (39) of equation (2), respectively.
Ify is a sufficiently smooth function, then the overall error(44)ex=Dc0αyx−PDc0αyx−fα0,…,ασt,y+Pfα0,…,ασt,yJsatisfies(45)exL20,1⟶0,asr⟶∞,orJ⟶∞.Proof.
Iffα0,…,ασt,y∈Hn0,1, then by Lemma 1, we can write(46)fα0,…,ασt,yJ−PJfα0,…,ασt,yJL20,1≤c1r−1−n2J−1−nfα0,…,ασnt,yJL20,1.
Using the Lipschitz conditions (3) and (41), we can write the following bound:(47)fα0,…,ασt,y−fα0,…,ασt,yJL20,1≤ℒσ∑j=0σDc0αjy−yJL20,1=ℒσ∑j=0σℐ0nj−αjDnjy−yJL20,1≤ℒσ∑j=0σCnj−αjDnjy−yJL20,1.
Now, applying equation (41) and Theorem 2 [24], we have(48)fα0,…,ασt,y−fα0,…,ασt,yJL20,1≤ℒσ∑j=0σCnj−αjCDnjr−11−nynj−1Hn,r−10,1≤ℒσ∑j=0σCnj−αjCDnjr−11−n2J−1−nyn+nj−1L20,1.
Also, using Lemma1 and equation (41), we have(49)Dc0αyx−PDc0αyxL20,1≤c2r−1−n2J−1−nDc0αyL20,1=≤c2r−1−n2J−1−nℐ0n−αDnyL20,1≤Cn−αc2r−1−n2J−1−nynL20,1.
Subtracting equation (2) from equation (39), one can write(50)ex=Dc0αyx−PDc0αyx−fα0,…,ασt,y+Pfα0,…,ασt,yJ+fα0,…,ασt,yJ−fα0,…,ασt,yJ.
Taking norms from both sides of equation (50) and using equations (46)–(49), we end up with(51)exL20,1≤Dc0αyx−PDc0αyxL20,1+fα0,…,ασt,y−fα0,…,ασt,yJL20,1+fα0,…,ασt,yJ−Pfα0,…,ασt,yJL20,1≤Cn−αc2r−1−n2J−1−nynL20,1+ℒσ∑j=0σCnj−αjCDnjr−11−n2J−1−nyn+nj−1L20,1+c1r−1−n2J−1−nfα0,…,ασnt,yJL20,1.
Suppose thatC1≔maxCnj−αjCDnj (j=0,…,σ) and C≔maxCn−αc2,C1,c1; then, we can obtain the following from equation (51):(52)exL20,1≤Cr−1−n2J−1−nynL20,1+ℒσ∑j=0σyn+nj−1L20,1+fα0,…,ασnt,yJL20,1.
It is easy to show thatexL20,1⟶0 as r⟶∞ or J⟶∞.


## 4. Numerical Implementation

In this section, we reported the numerical results for some examples to show the accuracy and efficiency of the method. To this end, we have performed all numerical computations in Maple and Matlab simultaneously. Wherever collocation nodes have not been reported, we have used Legendre nodes.Example 1.
Let us dedicate the first example to the following one.(53)Dc0αyx+Dc0α1yx+yx2=x2−x2+2x2−αΓ3−α+2x2−α1Γ3−α1−2x1−α1Γ2−α1,y0=0,y′0=−1.
The exact solution is reported in [4] as follows:(54)yx=x2−x.
To show the ability and efficiency of the proposed method, Tables1 and 2 are reported. The absolute error for proposed method is compared with the fractional-Lucas optimization method [4] and Chebyshev wavelet method [26] in Table 1. The results illustrate that the proposed method is flexible against other methods and gives a better approximation. We show L2–error, L∞–error, and CPU time for different values of α and α1 taking μ=0.5 in Table 2.Table 1
Comparison of the absolute error for Example1.
xOur methodFractional-Lucas optimization method [4] (m=2)Chebyshev wavelet method [26] (M=6,k=4)Chebyshev nodes (r=5,J=1,μ=0.5)Legendre nodes (r=5,J=1,μ=0.5)0.11.3620×10−492.0240×10−492.8022×10−178.8658×10−60.26.2780×10−495.2750×10−491.1208×10−148.5359×10−60.31.3344×10−489.4020×10−492.5219×10−148.1318×10−60.42.1909×10−481.4212×10−484.4835×10−147.6897×10−60.53.1020×10−481.9086×10−487.0055×10−147.1843×10−60.64.0646×10−482.4139×10−481.0087×10−136.7665×10−60.74.9381×10−482.8187×10−481.3730×10−136.3058×10−60.85.8985×10−483.3089×10−481.7934×10−135.8497×10−60.96.8068×10−483.7540×10−482.2697×10−135.4018×10−6Table 2
L2–error, L∞–error, and CPU time taking different values for α and α1 for Example 1.
α=1.75α=1.80α1=0.2α1=0.5α1=0.8α1=0.2α1=0.5α1=0.8L2–error5.52×10−492.41×10−482.76×10−481.35×10−471.16×10−471.33×10−47L∞–error8.92×10−493.66×10−484.24×10−482.10×10−471.78×10−472.07×10−47CPU time0.9060.9691.0310.9530.9850.968Example 2.
The second example is dedicated to the following equation:(55)Dc0αyx+yx2=x+xα+1Γα+22,y0=0y′0=−1.
The exact solution is reported in [5] as follows:(56)yx=xα+1Γα+2,0<α≤2.
In Table3, we compare the maximum of absolute value error of our method with fractional-order generalized Laguerre functions (FGLFs) [5]. In this example, we set the value of μ equal to α. L2–errors taking different vales for r and absolute value of error for different values of α are plotted in Figure 1.Table 3
The maximum of absolute value error taking different values forr and α for Example 2.
r=8r=12α=0.5α=0.7α=0.8α=0.5α=0.7α=0.8Our method8.31×10−469.72×10−51.56×10−47.55×10−461.89×10−53.32×10−5FGLFs [5]1.42×10−144.28×10−36.08×10−31.42×10−141.74×10−31.87×10−3Figure 1
Plot ofL2–errors taking different vales for r (a) and absolute value of error for different values of α (b) for Example 2.
(a)(b)Example 3.
We dedicated this example to the Bagley–Torvik equation [27, 28]:(57)D2yx+cD03/2yx+yx=1+x,y0=1,y′0=1.
The exact solution can be found in [27, 28], which is yx=x+1.
In Figure2, we plot the exact solution along with the approximate solution. Also, in this figure, the absolute value of error is reported for μ=0.5, r=3, and J=1. In Table 4, we compare the proposed method with the Bessel collocation method [28]. We observe that our method gives a better result than [28]. In this example, we put μ=0.5.Figure 2
Plot of the exact solution along with the approximate solution (a) and absolute value of error (b) for Example3.
(a)(b)Table 4
Comparison of the absolute error for Example3.
Our methodBessel collocation method [28]r=3N=6N=90.11.0×10−516.1919×10−169.3742×10−160.24.0×10−511.0292×10−153.9634×10−150.37.0×10−511.4779×10−154.2834×10−150.45.3×10−501.9697×10−153.2975×10−150.59.7×10−502.4941×10−152.0455×10−150.61.47×10−493.0365×10−151.0277×10−150.71.68×10−493.5805×10−153.4773×10−160.82.11×10−494.1090×10−156.9289×10−170.92.14×10−494.6047×10−152.3947×10−161.02.10×10−495.0500×10−151.2006×10−16


## 5. Conclusion

In this paper, we apply the pseudospectral method based on Müntz–Legendre wavelets to solve the multiorder fractional differential equations with Caputo fractional derivative. To this end, we represent the Caputo fractional derivative operator in the Müntz–Legendre wavelets. The results illustrate that by selecting the proper value forμ, the proposed method gives better results than others. The most important advantage of this method over other methods is its flexibility and ease of use. In most cases, the approximate solution is very close to the exact solution and we can almost say that the exact solution is obtained.



---

*Source: 9915551-2021-07-31.xml*