Abstract

The operational matrices of fractional-order integration for the Legendre and Chebyshev wavelets are derived. Block pulse functions and collocation method are employed to derive a general procedure for forming these matrices for both the Legendre and the Chebyshev wavelets. Then numerical methods based on wavelet expansion and these operational matrices are proposed. In this proposed method, by a change of variables, the multiorder fractional differential equations (MOFDEs) with nonhomogeneous initial conditions are transformed to the MOFDEs with homogeneous initial conditions to obtain suitable numerical solution of these problems. Numerical examples are provided to demonstrate the applicability and simplicity of the numerical scheme based on the Legendre and Chebyshev wavelets.

1. Introduction

Fractional-order differential equations (FODEs), as generalizations of classical integer-order differential equations, are increasingly used to model some problems in fluid flow, mechanics, viscoelasticity, biology, physics, engineering, and other applications. Fractional derivatives provide an excellent instrument for the description of memory and hereditary properties of various materials and processes [1ā€“6]. Fractional differentiation and integration operators are also used for extensions of the diffusion and wave equations [7]. The solutions of FODEs are much involved, because in general, there exists no method that yields an exact solution for FODEs, and only approximate solutions can be derived using linearization or perturbation methods. Several methods have been suggested to solve fractional differential equations (see [8] and references therein). Also there are different methods for solving MOFDEs as a special kind of FODEs [9ā€“17]. However, few papers have reported applications of wavelets in solving fractional differential equations [8, 18ā€“21]. In view of successful application of wavelet operational matrices in numerical solution of integral and differential equations [22ā€“27], together with the characteristics of wavelet functions, we believe that they should be applicable in solving MOFDEs.

In this paper, the operational matrices of fractional-order integrations are derived and a general procedure based on collocation method and block pulse functions for forming these matrices is presented. Then, by application of these matrices, a numerical method for solving MOFDEs with nonhomogeneous conditions is presented. In the proposed method by a change of variables the MOFDEs with nonhomogeneous conditions transforms to the MOFDEs with homogeneous conditions. This way, we are able to obtain the exact solutions for such problems. In the proposed method, the Legendre and Chebyshev wavelet expansions along with operational matrices of fractional-order integrations are employed to reduce the MOFDE to systems of nonlinear algebraic equations. Illustrative examples of nonlinear types are given to demonstrate the efficiency and applicability of the proposed method. As numerical results show, the proposed method is efficient and simple in implementation for both the Legendre and the Chebyshev wavelets. Moreover, for both these kinds of wavelets, numerical results have a good agreement with the exact solutions and the numerical results presented in other works.

The paper is organized as follows. In Section 2, we review some necessary definitions and mathematical preliminaries of fractional calculus and wavelets that are required for establishing our results. In Section 3 the Legendre and Chebyshev operational matrices of integration are derived. In Section 4 an application of the Legendre and Chebyshev operational matrices for solving the MOFDEs is presented. In Section 5 the proposed method is applied to several numerical examples. Finally, a conclusion is given in Section 6.

2. Basic Definitions

2.1. Fractional Calculus

We give some basic definitions and properties of the fractional calculus theory which are used further in this paper.

Definition 2.1. A real function š‘“(š‘”),š‘”>0, is said to be in the space š¶šœ‡, šœ‡āˆˆā„ if there exists a real number š‘(>šœ‡) such that š‘“(š‘”)=š‘”š‘š‘“1(š‘”), where š‘“1(š‘”)āˆˆš¶[0,āˆž], and it is said to be in the space š¶š‘›šœ‡ if š‘“(š‘›)āˆˆš¶šœ‡, š‘›āˆˆā„•.

Definition 2.2. The Riemann-Liouville fractional integration operator of order š›¼ā‰„0, of a function š‘“āˆˆš¶šœ‡, šœ‡ā‰„āˆ’1, is defined as follows [4]: š¼š›¼1š‘“(š‘”)=ī€œĪ“(š›¼)š‘”0(š‘”āˆ’šœ)š›¼āˆ’1š¼š‘“(šœ)š‘‘šœ,0š‘“(š‘”)=š‘“(š‘”),(2.1) and according to [4], we have š¼š›¼š¼š›½š‘“(š‘”)=š¼š›¼+š›½š¼š‘“(š‘”),š›¼š¼š›½š‘“(š‘”)=š¼š›½š¼š›¼š¼š‘“(š‘”),š›¼š‘”šœ—=Ī“(šœ—+1)š‘”Ī“(š›¼+šœ—+1)š›¼+šœ—,(2.2) where š›¼,š›½ā‰„0, š‘”>0 and šœ—>āˆ’1.

Definition 2.3. The fractional derivative of order š›¼>0 in the Riemann-Liouville sense is defined as [4] š·š›¼ī‚€š‘‘š‘“(š‘”)=ī‚š‘‘š‘”š‘›š¼š‘›āˆ’š›¼š‘“(š‘”),(š‘›āˆ’1<š›¼ā‰¤š‘›),(2.3) where š‘› is an integer and š‘“āˆˆš¶š‘›1.
The Riemann-Liouville derivatives have certain disadvantages when trying to model real-world phenomena with fractional differential equations. Therefore, we will now introduce a modified fractional differential operator š·š›¼āˆ— proposed by Caputo [5].

Definition 2.4. The fractional derivative of order š›¼>0 in the Caputo sense is defined as [5] š·š›¼āˆ—1š‘“(š‘”)=ī€œĪ“(š‘›āˆ’š›¼)š‘”0(š‘”āˆ’šœ)š‘›āˆ’š›¼āˆ’1š‘“(š‘›)(šœ)š‘‘šœ,(š‘›āˆ’1<š›¼ā‰¤š‘›),(2.4) where š‘› is an integer š‘”>0 and š‘“āˆˆš¶š‘›1. Caputos integral operator has a useful property: š¼š›¼š·š›¼āˆ—š‘“(š‘”)=š‘“(š‘”)āˆ’š‘›āˆ’1ī“š‘˜=0š‘“(š‘˜)ī€·0+ī€øš‘”š‘˜š‘˜!,(š‘›āˆ’1<š›¼ā‰¤š‘›),(2.5) where š‘› is an integer š‘”>0 and š‘“āˆˆš¶š‘›1.

2.2. Wavelets

Wavelets constitute a family of functions constructed from dilations and translations 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 [24]: šœ“š‘Ž,š‘(š‘”)=|š‘Ž|āˆ’1/2šœ“ī‚€š‘”āˆ’š‘š‘Žī‚,š‘Ž,š‘āˆˆā„,š‘Žā‰ 0.(2.6)

If we restrict the parameters š‘Ž and š‘ to discrete values as š‘Ž=š‘Ž0āˆ’š‘˜, š‘=š‘›š‘0š‘Ž0āˆ’š‘˜, š‘Ž0>1, š‘0>0, for š‘› and š‘˜ positive integers, we have the following family of discrete wavelets: šœ“š‘˜,š‘›||š‘Ž(š‘”)=0||š‘˜/2šœ“ī€·š‘Žš‘˜0š‘”āˆ’š‘›š‘0ī€ø,(2.7) where šœ“š‘˜,š‘›(š‘”) forms a wavelet basis for šæ2(ā„). In particular, when š‘Ž0=2 and š‘0=1, šœ“š‘˜,š‘›(š‘”) forms an orthonormal basis. That is (šœ“š‘˜,š‘›(š‘”),šœ“š‘™,š‘š(š‘”))=š›æš‘˜š‘™š›æš‘›š‘š.

2.2.1. The Legendre Wavelets

The Legendre wavelets šœ“š‘›,š‘š(š‘”)=šœ“(š‘˜,Ģ‚š‘›,š‘š,š‘”) have four arguments; š‘˜āˆˆā„•, š‘›=1,2,ā€¦,2š‘˜āˆ’1, and Ģ‚š‘›=2š‘›āˆ’1; moreover, š‘š is the order of the Legendre polynomials and š‘” is the normalized time, and they are defined on the interval [0,1) as šœ“š‘›,š‘šāŽ§āŽŖāŽØāŽŖāŽ©ī‚™(š‘”)=1š‘š+22š‘˜/2š‘š‘šī€·2š‘˜ī€ø,š‘”āˆ’Ģ‚š‘›Ģ‚š‘›āˆ’12š‘˜ā‰¤š‘”<Ģ‚š‘›2š‘˜,0,otherwise,(2.8) where š‘š=0,1,ā€¦,š‘€āˆ’1 and š‘€ is a fixed positive integer. The coefficient āˆšš‘š+1/2 in (2.8) is for orthonormality, the dilation parameter is š‘Ž=2āˆ’š‘˜, and the translation parameter is š‘=Ģ‚š‘›2āˆ’š‘˜. Here, š‘ƒš‘š(š‘”) are the well-known Legendre polynomials of order š‘š which are orthogonal with respect to the weight function š‘¤(š‘”)=1 on the interval [āˆ’1,1] and satisfy the following recursive formula: š‘0(š‘”)=1,š‘1(š‘”)=š‘”,š‘š‘š+1ī‚€(š‘”)=2š‘š+1ī‚š‘š+1š‘”š‘š‘šī‚€š‘š(š‘”)āˆ’ī‚š‘š‘š+1š‘šāˆ’1(š‘”),š‘š=1,2,3,ā€¦.(2.9)

2.2.2. The Chebyshev Wavelets

The Chebyshev wavelets šœ“š‘›,š‘š(š‘”)=šœ“(š‘˜,Ģ‚š‘›,š‘š,š‘”) have four arguments; š‘˜āˆˆā„•, š‘›=1,2,ā€¦,2š‘˜āˆ’1, and Ģ‚š‘›=2š‘›āˆ’1; moreover, š‘š is the order of the Chebyshev polynomials of the first kind and š‘” is the normalized time, and they are defined on the interval [0,1) as šœ“š‘›,š‘šīƒÆ2(š‘”)=š‘˜/2ī‚š‘‡š‘šī€·2š‘˜ī€ø,š‘”āˆ’Ģ‚š‘›Ģ‚š‘›āˆ’12š‘˜ā‰¤š‘”<Ģ‚š‘›2š‘˜,0,otherwise,(2.10) where ī‚š‘‡š‘šāŽ§āŽŖāŽØāŽŖāŽ©1(š‘”)=āˆššœ‹ī‚™,š‘š=0,2šœ‹š‘‡š‘š(š‘”),š‘š>0,(2.11) where š‘š=0,1,ā€¦,š‘€āˆ’1 and š‘€ is a fixed positive integer. The coefficients in (2.11) are used for orthonormality. Here, š‘‡š‘š(š‘”) are the well-known Chebyshev polynomials of order š‘š which are orthogonal with respect to the weight function āˆšš‘¤(š‘”)=1/1āˆ’š‘”2 on the interval [āˆ’1,1] and satisfy the following recursive formula: š‘‡0(š‘”)=1,š‘‡1(š‘”)=š‘”,š‘‡š‘š+1(š‘”)=2š‘”š‘‡š‘š(š‘”)āˆ’š‘‡š‘šāˆ’1,š‘š=1,2,3ā€¦.(2.12) We should note that in dealing with the Chebyshev polynomials the weight function ī‚š‘¤=š‘¤(2š‘”āˆ’1) has to be dilated and translated as follows: š‘¤š‘›ī€·2(š‘”)=š‘¤š‘˜ī€øš‘”āˆ’Ģ‚š‘›,(2.13) to get orthogonal wavelets.

2.2.3. Function Approximation

A function š‘“(š‘”) defined over [0,1) may be expanded as follows: š‘“(š‘”)=āˆžī“āˆžš‘›=1ī“š‘š=0š‘š‘›,š‘ššœ“š‘›,š‘š(š‘”),(2.14) by the Legendre or Chebyshev wavelets, where š‘š‘›,š‘š=(š‘“(š‘”),šœ“š‘›,š‘š(š‘”)) in which (,) denotes the inner product.

If the infinite series in (2.14) is truncated, then (2.14) can be written as š‘“(š‘”)ā‰ˆ2š‘˜āˆ’1ī“š‘›=1š‘€āˆ’1ī“š‘š=0š‘š‘›,š‘ššœ“š‘›,š‘š(š‘”)=š¶š‘‡ĪØ(š‘”),(2.15) where š¶ and ĪØ(š‘”) are 2š‘˜āˆ’1š‘€Ć—1 matrices given by ī€ŗš‘š¶=10,š‘11,ā€¦,š‘1š‘€āˆ’1,š‘20,ā€¦,š‘2š‘€āˆ’1,ā€¦,š‘2š‘˜āˆ’10,ā€¦,š‘2š‘˜āˆ’1š‘€āˆ’1ī€»š‘‡,ī€ŗšœ“ĪØ=10(š‘”),šœ“11(š‘”),ā€¦,šœ“1š‘€āˆ’1(š‘”),šœ“20(š‘”),ā€¦,šœ“2š‘€āˆ’1(š‘”),ā€¦,šœ“2š‘˜āˆ’10(š‘”),ā€¦,šœ“2š‘˜āˆ’1š‘€āˆ’1ī€»(š‘”)š‘‡.(2.16)

Taking the collocation points š‘”š‘–=(2š‘–āˆ’1)2š‘˜š‘€,š‘–=1,2,ā€¦,š‘š,(2.17) where š‘š=2š‘˜āˆ’1š‘€, we define the wavelet matrix Ī¦š‘šĆ—š‘š as Ī¦š‘šĆ—š‘š=ī‚ƒĪØī‚€1ī‚ī‚€32š‘š,ĪØī‚ī‚€2š‘š,ā€¦,ĪØ2š‘šāˆ’12š‘šī‚ī‚„.(2.18)

Indeed Ī¦š‘šĆ—š‘š has the following form: Ī¦š‘šĆ—š‘š=āŽ”āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ£āŽ¤āŽ„āŽ„āŽ„āŽ„āŽ„āŽ„āŽ„āŽ„āŽ¦š“0ā‹Æ000š“0ā‹Æ000š“ā‹Æ0ā‹®ā‹®ā‹±ā‹±ā‹®00ā‹Æ0š“,(2.19) where š“ is a š‘€Ć—š‘€ matrix given by āŽ”āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ£šœ“š“=10ī‚€1ī‚šœ“2š‘š10ī‚€3ī‚2š‘šā‹Æšœ“10ī‚€2š‘€āˆ’1ī‚šœ“2š‘š11ī‚€1ī‚šœ“2š‘š11ī‚€3ī‚2š‘šā‹Æšœ“11ī‚€2š‘€āˆ’1ī‚šœ“2š‘šā‹®ā‹®ā‹®ā‹®1š‘€āˆ’1ī‚€1ī‚šœ“2š‘š1š‘€āˆ’1ī‚€3ī‚2š‘šā‹Æšœ“1š‘€āˆ’1ī‚€2š‘€āˆ’1ī‚āŽ¤āŽ„āŽ„āŽ„āŽ„āŽ„āŽ„āŽ„āŽ¦2š‘š.(2.20)

For example, for š‘€=4 and š‘˜=2, the Legendre and Chebyshev matrices can be expressed as Ī¦8Ɨ8=āŽ”āŽ¢āŽ¢āŽ£āŽ¤āŽ„āŽ„āŽ¦š“00š“,(2.21) where for the Legendre matrix we have āŽ”āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ£āŽ¤āŽ„āŽ„āŽ„āŽ„āŽ„āŽ„āŽ¦š“=1.414211.414211.414211.41421āˆ’1.83712āˆ’.6123720.6123721.837121.08703āˆ’1.28468āˆ’1.284681.087030.2630851.25697āˆ’1.25697āˆ’.263085,(2.22) and for the Chebyshev matrix we have āŽ”āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ£āŽ¤āŽ„āŽ„āŽ„āŽ„āŽ„āŽ„āŽ¦š“=1.128381.128381.128381.12838āˆ’1.19683āˆ’.398942āˆ’.3989421.196830.199471āˆ’1.39630āˆ’1.396300.1994710.8976211.09709āˆ’1.09709āˆ’.897621.(2.23)

3. Operational Matrix of Fractional Integration

The integration of the vector ĪØ(š‘”) defined in (2.16) can be obtained by ī€œš‘”0ĪØ(šœ)š‘‘šœā‰ˆš‘ƒĪØ(š‘”),(3.1) where š‘ƒ is the š‘šĆ—š‘š operational matrix for integration [24].

Now, we derive the wavelet operational matrix of fractional integration. For this purpose, we rewrite the Riemann-Liouville fractional integration as (š¼š›¼1š‘“)(š‘”)=ī€œĪ“(š›¼)š‘”0(š‘”āˆ’šœ)š›¼āˆ’11š‘“(šœ)š‘‘šœ=š‘”Ī“(š›¼)š›¼āˆ’1āˆ—š‘“(š‘”),(3.2) where š›¼āˆˆā„ is the order of the integration and š‘”š›¼āˆ’1āˆ—š‘“(š‘”) denotes the convolution product of š‘”š›¼āˆ’1 and š‘“(š‘”). Now if š‘“(š‘”) is expanded by the Legendre and Chebyshev wavelets, as shown in (2.15), the Riemann-Liouville fractional integration becomes (š¼š›¼1š‘“)(š‘”)=Ī“š‘”(š›¼)š›¼āˆ’1āˆ—š‘“(š‘”)ā‰ˆš¶š‘‡1Ī“ī€½š‘”(š›¼)š›¼āˆ’1ī€¾āˆ—ĪØ(š‘”).(3.3)

So that, if š‘”š›¼āˆ’1āˆ—š‘“(š‘”) can be integrated, then by expanding š‘“(š‘”) in the Legendre and Chebyshev wavelets, the Riemann-Liouville fractional integration can be solved via the Legendre and Chebyshev wavelets.

Also, we define an š‘š-set of block pulse functions (BPFs) as š‘š‘–īƒÆš‘–(š‘”)=1,š‘šā‰¤š‘”<(š‘–+1)š‘š,0,otherwise,(3.4) where š‘–=0,1,2,ā€¦,(š‘šāˆ’1).

The functions š‘š‘–(š‘”) are disjoint and orthogonal, that is, š‘š‘–(š‘”)š‘š‘™ī‚»š‘(š‘”)=0,š‘–ā‰ š‘™,š‘–(ī€œš‘”),š‘–=š‘™,10š‘š‘–(šœ)š‘š‘™īƒÆ1(šœ)š‘‘šœ=0,š‘–ā‰ š‘™,š‘š,š‘–=š‘™.(3.5)

Similarly, the Legendre and Chebyshev wavelets may be expanded into š‘š-term block pulse functions (BPFs) as follows: ĪØš‘š(š‘”)=Ī¦š‘šĆ—š‘ššµš‘š(š‘”),(3.6) where šµš‘š(š‘”)=[š‘0(š‘”),š‘1(š‘”),ā€¦,š‘š‘–(š‘”),ā€¦,š‘š‘šāˆ’1(š‘”)]š‘‡.

In [28], Kilicman and Al Zhour have given the block pulse operational matrix of the fractional integration š¹š›¼ as follows: ī€·š¼š›¼šµš‘šī€ø(š‘”)ā‰ˆš¹š›¼šµš‘š(š‘”),(3.7) where š¹š›¼=1š‘šš›¼1āŽ”āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ£Ī“(š›¼+2)1šœ‰1šœ‰2ā‹Æšœ‰š‘šāˆ’101šœ‰1ā‹Æšœ‰š‘šāˆ’2001ā‹Æšœ‰š‘šāˆ’3āŽ¤āŽ„āŽ„āŽ„āŽ„āŽ„āŽ„āŽ„āŽ„āŽ¦000ā‹±ā‹®00001,(3.8) and šœ‰š‘–=(š‘–+1)š›¼+1āˆ’2š‘–š›¼+1+(š‘–āˆ’1)š›¼+1. Next, we derive the Legendre and Chebyshev wavelet operational matrices of the fractional integration. Let ī€·š¼š›¼ĪØš‘šī€ø(š‘”)ā‰ˆš‘ƒš›¼š‘šĆ—š‘šĪØš‘š(š‘”),(3.9) where the matrix š‘ƒš›¼š‘šĆ—š‘š is called the wavelet operational matrix of the fractional integration. Using (3.6) and (3.7), we have ī€·š¼š›¼ĪØš‘šī€øī€·š¼(š‘”)ā‰ˆš›¼Ī¦š‘šĆ—š‘ššµš‘šĆ—š‘šī€ø(š‘”)=Ī¦š‘šĆ—š‘šī€·š¼š›¼šµš‘šī€ø(š‘”)ā‰ˆĪ¦š‘šĆ—š‘šš¹š›¼šµš‘š(š‘”),(3.10) and from (3.9) and (3.10), we get š‘ƒš›¼š‘šĆ—š‘šĪØš‘š(š‘”)ā‰ˆš‘ƒš›¼š‘šĆ—š‘šĪ¦š‘šĆ—š‘ššµš‘š(š‘”)ā‰ˆĪ¦š‘šĆ—š‘šš¹š›¼šµš‘š(š‘”).(3.11)

Thus, the Legendre and Chebyshev wavelet operational matrices of the fractional integration š‘ƒš›¼š‘šĆ—š‘š can be approximately expressed by š‘ƒš›¼š‘šĆ—š‘š=Ī¦š‘šĆ—š‘šš¹š›¼Ī¦āˆ’1š‘šĆ—š‘š.(3.12) Also, from (3.3) and (3.12), we obtain (š¼š›¼š‘“)(š‘”)ā‰ˆš¶š‘‡Ī¦š‘šĆ—š‘šš¹š›¼šµš‘š(š‘”).(3.13) Here, for the Legendre and Chebyshev wavelets, we can obtain Ī¦š‘šĆ—š‘šš¹š›¼ as follows: Ī¦š‘šĆ—š‘šš¹š›¼=āŽ”āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ£šµ1šµ2šµ3Ģ‚š‘˜ā‹Æšµ0šµ1šµ2Ģ‚š‘˜ā‹Æšµāˆ’100šµ1Ģ‚š‘˜ā‹Æšµāˆ’2000ā‹±ā‹®0000šµ1āŽ¤āŽ„āŽ„āŽ„āŽ„āŽ„āŽ„āŽ„āŽ„āŽ¦,(3.14) where šµš‘Ÿ, Ģ‚Ģ‚š‘Ÿ=1,2,ā€¦,š‘˜(š‘˜=2š‘˜āˆ’1), are š‘€Ć—š‘€ matrices given by šµš‘Ÿ=ī‚€š‘(š‘Ÿ)š‘–š‘—ī‚š‘€Ć—š‘€,(3.15) such that for š‘Ÿ=1, š‘(1)š‘–š‘—=š‘—ī“š‘™=1š‘Žš‘–š‘™šœ‰ā€²š‘—āˆ’š‘™,šœ‰ī…ž0=1,(3.16) and for š‘Ÿ=2,3,ā€¦,š‘€,šµš‘Ÿ=ī‚€š‘(š‘Ÿ)š‘–š‘—ī‚š‘€Ć—š‘€Ģ‚š‘,š‘™=2,3,ā€¦,š‘˜,(š‘Ÿ)š‘–š‘—=š‘€ī“š‘™=1š‘Žš‘–š‘™šœ‰ī…ž(š‘Ÿāˆ’1)š‘€āˆ’š‘™+š‘—,(3.17) where š‘Žš‘–š‘—=ĪØ1(š‘–āˆ’1)ī‚µ2š‘—āˆ’1ī‚¶2š‘š,(3.18) and šœ‰ī…žš‘—=(1/š‘šš›¼Ī“(š›¼+2))šœ‰š‘—, š‘—=0,1,ā€¦š‘šāˆ’1.

Now from (3.12) and (3.14), we have š‘ƒš›¼š‘šĆ—š‘š=āŽ”āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ¢āŽ£š·1š·2š·3Ģ‚š‘˜ā‹Æš·0š·1š·2Ģ‚š‘˜ā‹Æš·āˆ’100š·1Ģ‚š‘˜ā‹Æš·āˆ’2000ā‹±ā‹®0000š·1āŽ¤āŽ„āŽ„āŽ„āŽ„āŽ„āŽ„āŽ„āŽ„āŽ¦,(3.19) where š·š‘–=šµš‘–š“āˆ’1Ģ‚š‘˜,š‘–=1,2,ā€¦.

4. Applications and Results

In this section, the Legendre and Chebyshev wavelet expansions together with their operational matrices of fractional-order integration are used to obtain numerical solution of MOFDEs.

Consider the following nonlinear MOFDE: š›½0(š‘”)š·š›¼āˆ—š‘¦(š‘”)+š‘ 1ī“š‘–=1š›½š‘–(š‘”)š·š›¼š‘–āˆ—š‘¦(š‘”)+š‘ 2ī“š‘–=1š›¾š‘–(š‘”)š‘¦š‘–(š‘”)=š‘“(š‘”),š‘¦(š‘—)(0)=š‘š‘—,š‘—=0,1,ā€¦š‘›āˆ’1,(4.1) where 0ā‰¤š‘”<1, š‘›āˆ’1<š›¼ā‰¤š‘›, 0<š›¼1<š›¼2<ā‹Æ<š›¼š‘ 1<š›¼, š‘›, š‘ 1 and š‘ 2 are fixed positive integers, š·š›¼āˆ— denotes the Caputo fractional derivative of order š›¼, š‘“ is a known function of š‘”, š‘š‘—, š‘—=0,1,ā€¦,š‘›āˆ’1, are arbitrary constants, and š‘¦ is an unknown function to be determined later. To solve this problem, we apply the following scheme.

Suppose š‘¦(š‘”)=š‘›āˆ’1ī“š‘—=0š‘š‘—š‘”š‘—š‘—!+š‘Œ(š‘”).(4.2)

Under this change of variable, we have š‘¦š‘–īƒ©(š‘”)=š‘›āˆ’1ī“š‘—=0š‘š‘—š‘”š‘—īƒŖš‘—!+š‘Œ(š‘”)š‘–=š‘–ī“š‘ =0āŽ›āŽœāŽœāŽš‘–š‘ āŽžāŽŸāŽŸāŽ īƒ©š‘›āˆ’1ī“š‘—=0š‘š‘—š‘”š‘—īƒŖš‘—!š‘ š‘Œ(š‘”)š‘–āˆ’š‘ ,š·š›¼āˆ—š‘¦(š‘”)=š·š›¼āˆ—š·š‘Œ(š‘”),š›¼š‘–āˆ—š‘¦(š‘”)=š·š›¼š‘–āˆ—š‘Œ(š‘”),š‘›āˆ’1<š›¼š‘–š·<š›¼,š›¼š‘–āˆ—š‘¦(š‘”)=š·š›¼š‘–āˆ—š‘Œ(š‘”)+ā„Žš‘–(š‘”),š›¼š‘–ā‰¤š‘›āˆ’1,(4.3) where ā„Žš‘–(š‘”)=š·š›¼š‘–āˆ—īƒ©š‘›āˆ’1ī“š‘—=0š‘š‘—š‘”š‘—īƒŖš‘—!.(4.4)

Now, by substituting (4.3) into (4.1), we transform the nonlinear MOFDE (4.1) with nonhomogeneous conditions to a nonlinear MOFDE with homogeneous conditions as follows: Ģ‚š›½0(š‘”)š·š›¼āˆ—š‘Œ(š‘”)+š‘ 1ī“š‘–=1Ģ‚š›½š‘–(š‘”)š·š›¼š‘–āˆ—š‘Œ(š‘”)+š‘ 2ī“š‘–š‘–=1ī“š‘ =0Ģ‚š›¾š‘–,š‘ (š‘”)š‘Œš‘–āˆ’š‘ ī(š‘”)=š‘“(š‘”),š‘Œ(š‘—)(0)=0,š‘—=0,1,ā€¦š‘›āˆ’1,(4.5) where Ģ‚š›½0(š‘”)=š›½0Ģ‚š›½(š‘”),š‘–(š‘”)=š›½š‘–(š‘”),š‘›āˆ’1<š›¼š‘–Ģ‚š›½<š›¼,š‘–(š‘”)=š›½š‘–(š‘”)ā„Žš‘–(š‘”),š›¼š‘–ā‰¤š‘›āˆ’1,Ģ‚š›¾š‘–,š‘ (š‘”)=š›¾š‘–āŽ›āŽœāŽœāŽš‘–š‘ āŽžāŽŸāŽŸāŽ īƒ©(š‘”)š‘›āˆ’1ī“š‘—=0š‘š‘—š‘”š‘—īƒŖš‘—!š‘ ,(4.6) and īš‘“(š‘”) is a known function.

We assume that š·š›¼āˆ—š‘Œ(š‘”) is given by š·š›¼āˆ—š‘Œ(š‘”)=š¾š‘‡š‘šĪØš‘š(š‘”),(4.7) where š¾š‘‡š‘š is an unknown vector and ĪØš‘š(š‘”) is the vector which is defined in (3.6). By using initial conditions and (2.4), we have š·š›¼š‘–āˆ—š‘Œ(š‘”)=š¾š‘‡š‘šš‘ƒš›¼āˆ’š›¼š‘–š‘šĆ—š‘šĪØš‘š(š‘”),š‘–=1,2,ā€¦,š‘ 1,(4.8)š‘Œ(š‘”)=š¾š‘‡š‘šš‘ƒš›¼š‘šĆ—š‘šĪØš‘š(š‘”).(4.9)

Since ĪØš‘š=Ī¦š‘šĆ—š‘ššµš‘š(š‘”), (4.9) can be rewritten as follows: š‘Œ(š‘”)=š¾š‘‡š‘šš‘ƒš›¼š‘šĆ—š‘šĪ¦š‘šĆ—š‘ššµš‘šī€ŗš‘Ž(š‘”)=1,š‘Ž2,ā€¦,š‘Žš‘šī€»šµš‘š(š‘”).(4.10)

Now by using (3.4), we obtain [š‘Œ](š‘”)š‘–āˆ’š‘ =ī€ŗš‘Ž1š‘–āˆ’š‘ ,š‘Ž2š‘–āˆ’š‘ ,ā€¦,š‘Žš‘šš‘–āˆ’š‘ ī€»šµš‘š(š‘”),š‘–=1,2,ā€¦,š‘ 2,š‘ =0,1,ā€¦,š‘–.(4.11)

Moreover, we expand functions Ģ‚š›½š‘–(š‘”), Ģ‚š›¾š‘–,š‘ (š‘”), and īš‘“(š‘”) by wavelets as follows: Ģ‚š›½š‘–(š‘”)=ĪØš‘š(š‘”)š‘‡Ģƒš›½š‘–,š‘š,š‘–=0,1,ā€¦,š‘ 1,Ģ‚š›¾š‘–,š‘ (š‘”)=ĪØš‘š(š‘”)š‘‡Ģƒš›¾š‘–,š‘ ,š‘š,š‘–=1,2,ā€¦,š‘ 2ī,š‘ =0,1,ā€¦,š‘–,š‘“(š‘”)=š‘“š‘‡š‘šĪØš‘š(š‘”),(4.12) where Ģƒš›½š‘–,š‘š, Ģƒš›¾š‘–,š‘š, and š‘“š‘‡š‘š are known vectors. By substituting (4.7), (4.8), and (4.11)-(4.12) into (4.5), we obtain ĪØš‘š(š‘”)š‘‡Ģƒš›½0,š‘šš¾š‘‡š‘šĪØš‘š(š‘”)+š‘ 1ī“š‘–=1ĪØš‘š(š‘”)š‘‡Ģƒš›½š‘–,š‘šš¾š‘‡š‘šš‘ƒš›¼āˆ’š›¼š‘–š‘šĆ—š‘šĪØš‘š(š‘”)+š‘ 2ī“š‘–š‘–=1ī“š‘ =0ĪØš‘š(š‘”)š‘‡Ģƒš›¾š‘–,š‘ ,š‘šī€ŗš‘Ž1š‘–āˆ’š‘ ,š‘Ž2š‘–āˆ’š‘ ,ā€¦,š‘Žš‘šš‘–āˆ’š‘ ī€»šµš‘š(š‘”)=šµš‘š(š‘”)š‘‡[]1,1,ā€¦,1š‘‡š‘“š‘‡š‘šĪØš‘š(š‘”).(4.13)

Now, from (3.12), (4.10), and (4.13), we have šµš‘š(š‘”)š‘‡īƒ¬Ī¦š‘‡š‘šĆ—š‘šĢƒš›½0,š‘šī€ŗš‘Ž1,š‘Ž2,ā€¦,š‘Žš‘šī€»š¹š›¼āˆ’1+š‘ 1ī“š‘–=1Ī¦š‘‡š‘šĆ—š‘šĢƒš›½š‘–,š‘šī€ŗš‘Ž1,š‘Ž2,ā€¦,š‘Žš‘šī€»š¹š›¼āˆ’1š¹š›¼š‘–š›¼āˆ’š›¼+š‘ 2ī“š‘–š‘–=1ī“š‘ =0Ī¦š‘‡š‘šĆ—š‘šĢƒš›¾š‘–,š‘ ,š‘šī€ŗš‘Ž1š‘–āˆ’š‘ ,š‘Ž2š‘–āˆ’š‘ ,ā€¦,š‘Žš‘šš‘–āˆ’š‘ ī€»īƒ­šµš‘š(š‘”)=šµš‘š(š‘”)š‘‡[]1,1,ā€¦,1š‘‡š‘“š‘‡š‘šĪ¦š‘šĆ—š‘ššµš‘š(š‘”).(4.14)

This is a nonlinear algebraic equation for unknown vector [š‘Ž1,š‘Ž2,ā€¦,š‘Žš‘š]. Here, by taking collocation points, expressed in (2.17), we transform (4.14) into a nonlinear system of algebraic equations. This nonlinear system can be solved by the Newton iteration method for unknown vector [š‘Ž1,š‘Ž2,ā€¦,š‘Žš‘š]. Therefore, š‘Œ(š‘”) as the solution of (4.5) is š‘Œī€ŗš‘Ž(š‘”)=1,š‘Ž2,ā€¦,š‘Žš‘šī€»šµš‘š(š‘”),(4.15) and finally, š‘¦(š‘”) as the solution of (4.1) will be š‘¦(š‘”)=š‘›āˆ’1ī“š‘—=0š‘š‘—š‘”š‘—+ī€ŗš‘Žš‘—!1,š‘Ž2,ā€¦,š‘Žš‘šī€»šµš‘š(š‘”).(4.16)

In cases where all of coefficients Ģ‚š›½š‘–(š‘”) and Ģ‚š›¾š‘–,š‘ (š‘”), are constants, that is, Ģ‚š›½š‘–Ģ‚š›½(š‘”)=š‘– and Ģ‚š›¾š‘–,š‘ (š‘”)=Ģ‚š›¾š‘–,š‘  we can reduce (4.13) and (4.14), respectively, as follows: Ģ‚š›½0š¾š‘‡š‘šĪØš‘š(š‘”)+š‘ 1ī“š‘–=1Ģ‚š›½š‘–š¾š‘‡š‘šš‘ƒš›¼āˆ’š›¼š‘–š‘šĆ—š‘šĪØš‘š(š‘”)+š‘ 2ī“š‘–š‘–=1ī“š‘ =0Ģ‚š›¾š‘–,š‘ ī€ŗš‘Ž1š‘–āˆ’š‘ ,š‘Ž2š‘–āˆ’š‘ ,ā€¦,š‘Žš‘šš‘–āˆ’š‘ ī€»šµš‘š(š‘”)=š‘“š‘‡š‘šĪØš‘šīƒ¬ī€ŗš‘Ž(š‘”),1,š‘Ž2,ā€¦,š‘Žš‘šī€»īƒ¬Ģ‚š›½0š¹š›¼āˆ’1+š¹š›¼āˆ’1š‘ 1ī“š‘–=1Ģ‚š›½š‘–š¹š›¼āˆ’š›¼š‘–īƒ­+š‘ 2ī“š‘–š‘–=1ī“š‘ =0Ģ‚š›¾š‘–,š‘ ī€ŗš‘Ž1š‘–āˆ’š‘ ,š‘Ž2š‘–āˆ’š‘ ,ā€¦,š‘Žš‘šš‘–āˆ’š‘ ī€»īƒ­šµš‘š(š‘”)=š‘“š‘‡š‘šĪ¦š‘šĆ—š‘ššµš‘š(š‘”).(4.17)

This is a nonlinear system of algebraic equations for unknown vector [š‘Ž1,š‘Ž2,ā€¦,š‘Žš‘š] and can be solved directly without use of collocation points by the Newton iteration method. Here, as before we obtain š‘Œī€ŗš‘Ž(š‘”)=1,š‘Ž2,ā€¦,š‘Žš‘šī€»šµš‘š(š‘”),(4.18) as the solution of (4.5), and finally, š‘¦(š‘”) as the solution of (4.1) will be š‘¦(š‘”)=š‘›āˆ’1ī“š‘—=0š‘š‘—š‘”š‘—+ī€ŗš‘Žš‘—!1,š‘Ž2,ā€¦,š‘Žš‘šī€»šµš‘š(š‘”).(4.19)

5. Numerical Examples

In this section we demonstrate the efficiency of the proposed wavelet collocation method for the numerical solution of MOFDEs. These examples are considered because either closed form solutions are available for them, or they have also been solved using other numerical schemes, by other authors.

Example 5.1. Consider the homogeneous Bagley-Torvik equation [13, 16]: š·2āˆ—š‘¦(š‘”)+š·āˆ—0.5š‘¦(š‘”)+š‘¦(š‘”)=0,(5.1) subject to the following initial conditions: š‘¦(0)=1,š‘¦ī…ž(0)=0.(5.2) Here, we apply the method of Section 4 for solving this problem. By setting š‘¦(š‘”)=1+š‘Œ(š‘”).(5.3) Equation (5.1) can be rewritten as š·2āˆ—š‘Œ(š‘”)+š·āˆ—0.5š‘Œ(š‘”)+š‘Œ(š‘”)+1=0,(5.4) subject to the initial condition š‘Œ(0)=š‘Œī…ž(0)=0.(5.5) By using the proposed method, we get a system of algebraic equations for both the Legendre and the Chebyshev wavelets for (5.4) as follows: ī€ŗš‘Ž1,š‘Ž2,ā€¦,š‘Žš‘šī€»ī‚ƒš¹2āˆ’1+š¹2āˆ’1š¹1.5ī‚„+ī€ŗš‘Ž1,š‘Ž2,ā€¦,š‘Žš‘šī€»+[]1,1,ā€¦,1=0.(5.6) Since the linear system of algebraic equations (5.6) is the same for both kinds of wavelets, we have the same numerical solution for the Legendre and Chebyshev wavelets. By solving this linear system, we get numerical solutions for (5.4) as follows: š‘Œī€ŗš‘Ž(š‘”)=1,š‘Ž2,ā€¦,š‘Žš‘šī€»šµš‘š(š‘”).(5.7) Finally the solution of the original problem (5.1) is š‘¦ī€ŗš‘Ž(š‘”)=1+1,š‘Ž2,ā€¦,š‘Žš‘šī€»šµš‘š(š‘”).(5.8) Figure 1 shows the behavior of the numerical solution for š‘š=160(š‘€=5,š‘˜=6). As Figure 1 shows, this method is very efficient for numerical solution of this problem and the solution can be derived in a large interval [0,40].

Example 5.2. Consider the nonhomogeneous Bagley-Torvik equation [10, 13ā€“15, 17, 20]: š‘Žš·2āˆ—š‘¦(š‘”)+š‘š·āˆ—1.5š‘¦(š‘”)+š‘š‘¦(š‘”)=š‘“(š‘”),(5.9) where š‘“(š‘”)=š‘(1+š‘”),(5.10) subject to the initial condition š‘¦(0)=š‘¦ī…ž(0)=1,(5.11) which has the exact solution š‘¦(š‘”)=1+š‘”.(5.12) Here, we solve this problem by applying the method of Section 4. By setting š‘¦(š‘”)=1+š‘”+š‘Œ(š‘”).(5.13) Equation (5.9) can be rewritten as š‘Žš·2āˆ—š‘Œ(š‘”)+š‘š·āˆ—1.5š‘Œ(š‘”)+š‘š‘Œ(š‘”)=0,(5.14) subject to the following initial conditions: š‘Œ(0)=š‘Œī…ž(0)=0.(5.15) The system of algebraic equations for both the Legendre and the Chebyshev wavelets for (5.14) has the following form: ī€ŗš‘Ž1,š‘Ž2,ā€¦,š‘Žš‘šī€»ī‚ƒš‘Žš¹2āˆ’1+š‘š¹2āˆ’1š¹0.5+š‘š¼š‘šĆ—š‘šī‚„=0.(5.16) Here, the linear system of algebraic equations (5.16) is nonsingular and so has only the trivial solution, that is, š‘Œ(š‘”)=0. Then the solution of the original problem (5.9) is š‘¦(š‘”)=1+š‘”,(5.17) which is the exact solution. It is mentionable that this equation has been solved by Diethelm and Luchko [10] (for š‘Ž=š‘=š‘=1) with error 1.60šøāˆ’4, Diethelm and Ford [17] (for š‘Ž=1, š‘=š‘=0.5) with error 5.62šøāˆ’3, while in [20] Li and Zhao obtained approximation solution by the Haar wavelet.

Example 5.3. Consider the following nonhomogenous MOFDE [13]: š‘Žš·š›¼āˆ—š‘¦(š‘”)+š‘š·š›¼2āˆ—š‘¦(š‘”)+š‘š·š›¼1āˆ—š‘¦(š‘”)+š‘’š‘¦(š‘”)=š‘“(š‘”),0<š›¼1ā‰¤1,1<š›¼2ā‰¤2,3<š›¼<4,(5.18) where š‘“(š‘”)=2š‘Ī“ī€·3āˆ’š›¼2ī€øš‘”2āˆ’š›¼2+2š‘Ī“ī€·3āˆ’š›¼1ī€øš‘”2āˆ’š›¼1ī€·š‘”+š‘’2ī€øāˆ’š‘”,(5.19) subject to the following initial conditions: š‘¦(0)=0,š‘¦ī…ž(0)=āˆ’1,š‘¦ī…žī…ž(0)=2,š‘¦ī…žī…žī…ž(0)=0,(5.20) which has the exact solution š‘¦(š‘”)=š‘”2āˆ’š‘”.(5.21) Here, we apply the method of Section 4 for solving this problem. By setting š‘¦(š‘”)=āˆ’š‘”+š‘”2+š‘Œ(š‘”).(5.22) Equation (5.18)can be rewritten as follows: š‘Žš·š›¼āˆ—š‘Œ(š‘”)+š‘š·š›¼2āˆ—š‘Œ(š‘”)+š‘š·š›¼1āˆ—š‘Œ(š‘”)+š‘’š‘Œ(š‘”)=0,0<š›¼1ā‰¤1,1<š›¼2ā‰¤2,3<š›¼<4,(5.23) subject to the initial conditions š‘Œ(0)=0,š‘Œī…ž(0)=0,š‘Œī…žī…ž(0)=0,š‘Œī…žī…žī…ž(0)=0.(5.24) The system of algebraic equations corresponding to the Legendre and Chebyshev wavelets for (5.23) has the following form: ī€ŗš‘Ž1,š‘Ž2,ā€¦,š‘Žš‘šī€»ī‚ƒš‘Žš¹š›¼āˆ’1+š‘š¹š›¼āˆ’1š¹š›¼āˆ’š›¼2+š‘š¹š›¼āˆ’1š¹š›¼āˆ’š›¼1+š‘’š¼š‘šĆ—š‘šī‚„=0.(5.25) Here, the linear system of algebraic (5.25) is nonsingular and so has only the trivial solution, that is, š‘Œ(š‘”)=0. Then the solution of the original problem (5.18) is š‘¦(š‘”)=š‘”2āˆ’š‘”,(5.26) which is the exact solution. This equation has been solved by El-Sayed et al. [13] (for š‘Ž=š‘=š‘=š‘’=1, š›¼1=0.77, š›¼2=1.44, and š›¼=0.91) with error 9.53šøāˆ’4.

Example 5.4. Consider the following nonlinear MOFDE [15]: š·3āˆ—š‘¦(š‘”)+š·āˆ—2.5[]š‘¦(š‘”)+š‘¦(š‘”)2=š‘”4,(5.27) subject to the initial conditions š‘¦(0)=š‘¦ī…ž(0)=0,š‘¦ī…žī…ž(0)=2,(5.28) which has the exact solution š‘¦(š‘”)=š‘”2.(5.29) Here, we apply the method of Section 4 for solving this problem. By setting š‘¦(š‘”)=š‘”2+š‘Œ(š‘”).(5.30) Equation (5.27) can be rewritten as follows: š·3āˆ—š‘Œ(š‘”)+š·āˆ—2.5[]š‘Œ(š‘”)+š‘Œ(š‘”)2+2š‘”2š‘Œ(š‘”)=0,(5.31) subject to the following initial conditions: š‘Œ(0)=0,š‘Œī…ž(0)=0,š‘Œī…žī…ž(0)=0.(5.32) The nonlinear algebraic equation corresponding to the Legendre and Chebyshev wavelets for (5.31) has the following form: šµš‘š(š‘”)š‘‡ī‚ƒ[]1,1,ā€¦,1š‘‡ī€ŗš‘Ž1,š‘Ž2,ā€¦,š‘Žš‘šī€»ī‚ƒš¹3āˆ’1+š¹3āˆ’1š¹0.5ī‚„+[]1,1,ā€¦,1š‘‡ī€ŗš‘Ž21,š‘Ž22,ā€¦,š‘Ž2š‘šī€»+Ī¦š‘‡š‘šĆ—š‘šĢƒš›¾š‘šī€ŗš‘Ž1,š‘Ž2,ā€¦,š‘Žš‘ššµī€»ī€»š‘š(š‘”)=0,(5.33) where Ģƒš›¾š‘š is a known constant vector corresponding to the kind of wavelet expansions and Ī¦š‘šĆ—š‘š is the wavelet matrix. By taking collocation points expressed in (2.17), we transform (5.33) into a nonlinear system of algebraic equations. Then applying Newton iteration method for solving this nonlinear system, we obtain only trivial solution, that is, š‘Œ(š‘”)=0. Then the solution of the original problem (5.27) is š‘¦(š‘”)=š‘”2,(5.34) which is the exact solution.

Example 5.5. Consider the following nonlinear MOFDE [14]: š‘Žš·3āˆ—š‘¦(š‘”)+š‘š·š›¼2āˆ—š‘¦(š‘”)+š‘š·š›¼1āˆ—[]š‘¦(š‘”)+š‘’š‘¦(š‘”)2=š‘“(š‘”),0<š›¼1ā‰¤1,1<š›¼2<2,(5.35) where š‘“(š‘”)=š‘š‘”1āˆ’š›¼1Ī“ī€·2āˆ’š›¼1ī€ø+š‘’š‘”2,(5.36) subject to the initial conditions š‘¦(0)=0,š‘¦ī…ž(0)=1,š‘¦ī…žī…ž(0)=0,(5.37) which has the exact solution š‘¦(š‘”)=š‘”.(5.38) Here, we apply the method of Section 4 for solving this problem. By setting š‘¦(š‘”)=š‘”+š‘Œ(š‘”).(5.39) Equation (5.35) can be rewritten as follows: š‘Žš·3āˆ—š‘Œ(š‘”)+š‘š·š›¼2āˆ—š‘Œ(š‘”)+š‘š·š›¼1āˆ—[]š‘Œ(š‘”)+š‘’š‘Œ(š‘”)2+2š‘’š‘”š‘Œ(š‘”)=0,(5.40) subject to the initial conditions š‘Œ(0)=0,š‘Œī…ž(0)=0,š‘Œī…žī…ž(0)=0.(5.41) The nonlinear algebraic equation corresponding to the Legendre and Chebyshev wavelets for (5.40) has the following form: šµš‘š(š‘”)š‘‡ī‚ƒ[]1,1,ā€¦,1š‘‡ī€ŗš‘Ž1,š‘Ž2,ā€¦,š‘Žš‘šī€»ī‚ƒš‘Žš¹3āˆ’1+š‘š¹3āˆ’1š¹3āˆ’š›¼2+š‘š¹3āˆ’1š¹3āˆ’š›¼1ī‚„+[]1,1,ā€¦,1š‘‡ī€ŗš‘Ž21,š‘Ž22,ā€¦,š‘Ž2š‘šī€»+Ī¦š‘‡š‘šĆ—š‘šĢƒš›¾š‘šī€ŗš‘Ž1,š‘Ž2,ā€¦,š‘Žš‘ššµī€»ī€»š‘š(š‘”)=0,(5.42) where Ģƒš›¾š‘š is a known constant vector corresponding to the kind of wavelet expansion and Ī¦š‘šĆ—š‘š is the wavelet matrix. By taking collocation points expressed in (2.17), we transform (5.42) into a nonlinear system of algebraic equations. By applying the Newton iteration method for solving this nonlinear system, we obtain only trivial solution, that is, š‘Œ(š‘”)=0. Then the solution of the original problem (5.35) is š‘¦(š‘”)=š‘”,(5.43) which is the exact solution.

6. Conclusion

In this paper a general formulation for the Legendre and Chebyshev wavelet operational matrices of fractional-order integration has been derived. Then a numerical method based on Legendre and Chebyshev wavelets expansions together with these matrices are proposed to obtain the numerical solutions of MOFDEs. In this proposed method, by a change of variables, the MOFDEs with nonhomogeneous conditions are transformed to the MOFDEs with homogeneous conditions. The exact solutions for some MOFDEs are obtained by our method. The proposed method is very simple in implementation for both the Legendre and the Chebyshev wavelets. As the numerical results show, the method is very efficient for the numerical solution of MOFDEs and only a few number of wavelet expansion terms are needed to obtain a good approximate solution for these problems.