# Chebyshev Wavelets Method for Solution of Nonlinear Fractional Integrodifferential Equations  in a Large Interval

Authors: M. H. Heydari; M. R. Hooshmandasl; F. M. Maalek Ghaini; Ming Li
Journal: Advances in Mathematical Physics
 (2013)

---

## Abstract

An efficient Chebyshev wavelets method for solving a class of nonlinear fractional integrodifferential equations in a large interval is developed, and a new technique for computing nonlinear terms in such equations is proposed. Existence of a unique solution for such equations is proved. Convergence and error analysis of the proposed method are investigated. Moreover in order to show efficiency of the proposed method, the new approach is compared with some numerical methods.


---

## Body

## 1. Introduction

Fractional integrodifferential equations (FIDEs) arise in modelling processes in applied sciences such as physics, engineering, and biology. The nonlinear fractional integro-differential equation (NFIDE) of the type(1)D*αy(x)=g(x,y(x))+λ∫0x‍K(x,t,y(t))dt,0<α≤1,x∈[0,b],y(0)=c,
where D*α is Caputo fractional derivative, α is a parameter describing the order of the fractional derivative, λ, c, and b are fixed constants, and g is a nonlinear continuous function, arise in the mathematical modelling of various physical phenomena, such as heat conduction in materials with memory. Moreover, these equations are encountered in combined conduction, convection, and radiation problems [1–3]. Therefore in recent years, numerous works have been focusing on the solution of these problems. Some of these methods are Adomian decomposition method (ADM) [4], fractional differential transform method (FDTM) [5], and the collocation method [6]. Most of these methods have been utilized in linear problems, and a few number of works have considered nonlinear problems. In [7] Rawashdeh applied Legendre wavelet method for solving fractional Voltera integro-differential equations in the form
(2)D*αy(x)=f(x)+p(x)y(x)+λ∫0x‍K(x,t)y(t)dt,0<α≤1,x∈[0,1].
Also in [8] Awawdeh et al. applied homotopy analysis method (HAM) for solution of (2). Mittal and Nigam [9] applied the ADM for (1) in the form
(3)D*αy(x)=f(x)+p(x)y(x)+λ∫0x‍K(x,t)[y(t)]ndt,0<α≤1,x∈[0,1],n≥1.
In view of successful application of wavelets in approximation theory [10–17], we will use the Chebyshev wavelets for solving a generalized form of the previous described equations of the form
(4)D*αy(x)=f(x)+p(x)F(y(x))+λ∫0x‍K(x,t)G(y(t))dt,0<α≤1,x∈[0,b],b>1,
with high nonlinearity in a large interval under the initial condition y(0)=c.Here, for simplicity, we assume thatf and p are continuous functions on [0,b], K(x,t) is continuous on 0≤t≤x≤b, and also F and G are analytic functions.The method is based on reducing the equation to a system of nonlinear algebraic equations by expanding the solution as Chebyshev wavelet bases with unknown coefficients. We note that, forα=1, (4) is an ordinary integro-differential equation and the method can be easily applied to it. Also the method is fast and mathematically simple and guarantees the necessary accuracy for a small number of basic functions. Moreover in order to show the efficiency of the proposed method the new approach is compared with the ADM, VIM [18–21] and ITM which has been proposed by Daftardar-Gejji and Jafari in 2006 [22, 23].The outline of this paper is as follows: in Section2, some basic definitions of the fractional calculus, Chebyshev wavelets, functional approximation, and operational matrix of fractional integration are described. In Section 3 the existence of a unique solution for the problem is established. In Section 4, Chebyshev wavelets method and error investigation are described. In Section 5, some numerical examples are solved by applying the methods of this article. Finally a conclusion is drawn in Section 6.


## 2. Basic Definitions

### 2.1. Fractional Calculus

Here, we give some basic definitions and properties of the fractional calculus theory which will be used further in this paper.Definition 1.
A real functionf(t),t>0, is said to be in the space Cμ,μ∈ℝ, if there exists a real number p(>μ) such that f(t)=tpf1(t), where f1(t)∈C[0,∞] and it is said to be in the space Cμn if f(n)∈Cμ,n∈ℕ.Definition 2.
The Riemann-Liouville fractional integration operator of orderα≥0 of a function f∈Cμ,μ≥-1, is defined as
(5)(Iαf)(t)={1Γ(α)∫0t‍(t-τ)α-1f(τ)dτ,α>0,f(t),α=0.
It has the following properties:
(6)(i)IαIβ=Iα+β,(ii)IαIβ=IβIα,(iii)(IαIβf)(t)=(IβIαf)(t),(iv)Iαtϑ=Γ(ϑ+1)Γ(α+ϑ+1)tα+ϑ,
where α,β≥0,t>0, and ϑ>-1.
Riemann-Liouville fractional derivative operator of orderα>0 is defined as
(7)Dαf(t)=(ddt)nIn-αf(t),(n-1<α≤n),
where n is an integer and f∈C1n.
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 operatorD*α proposed by Caputo [24].Definition 3.
The fractional derivative operator of orderα>0 in the Caputo sense is defined as
(8)D*αf(t)=1Γ(n-α)∫0t‍(t-τ)n-α-1f(n)(τ)dτ,(n-1<α≤n),
where n is an integer, t>0, and f∈C1n.
The Riemann-Liouville and Caputo operators have a useful property as(9)IαD*αf(t)=f(t)-∑k=0n-1‍f(k)(0+)tkk!,(n-1<α≤n),
where n is an integer, t>0, and f∈C1n.
For more details on the mathematical properties of fractional derivatives and integrals see [24].


### 2.2. Chebyshev Wavelets

Chebyshev waveletsψpq(t)=ψ(k,p^,q,t) have four arguments; k∈ℕ,p=1,2,…,2k-1, and p^=2p-1; moreover q is the degree of the Chebyshev polynomial of the first kind and t is the normalized time, that is, t∈[0,1]. They are defined on the interval [0,1] as [25]
(10)ψpq(t)={2k/2T~q(2kt-p^),p^-12k≤t≤p^+12k,0,otherwise,
where
(11)T~q(t)={1π,q=02πTq(t),q>0,q=0,1,…,M-1, and M is a fixed positive integer. The coefficients in (11) are used for orthonormality of the system. Here, {Tq(t),q∈ℕ∪{0}} is the set of well-known Chebyshev polynomials of degree m which are orthogonal with respect to the weight function w(t)=1/1-t2 on the interval [-1,1] and satisfy the following recursive formula:
(12)T0(t)=1,T1(t)=t,Tq+1(t)=2tTq(t)-Tq-1(t),q=1,2,3….
We should note that in dealing with Chebyshev polynomials the weight function w~(t)=w(2t-1) has to be dilated and translated as wp(t)=w(2kt-p^), to get orthogonal wavelets.


### 2.3. Function Approximation

An arbitrary functionf(t)∈L2(R) defined over [0,1] may be expanded into Chebyshev wavelet basis as
(13)f(t)=∑p=1∞‍∑q=0∞‍cpqψpq(t),
where cpq=(f(t),ψpq(t)) in which (,) denotes the inner product.If the infinite series in (13) is truncated, then (13) can be written as
(14)f(t)≃∑p=12k-1‍∑q=0M-1‍cpqψpq(t)=CTΨ(t),
where C and Ψ(t) are 2k-1M×1 matrices given by
(15)C=[c10,c11,…,c1M-1,c20,…,c2M-1,…,c2k-10,…,c2k-1M-1]T,Ψ(t)=[ψ10(t),ψ11(t),…,ψ1M-1(t),ψ20(t),…,ψ2M-1(t),…,ψ2k-10(t),…,ψ2k-1M-1(t)]T.
Taking the collocation points
(16)ti=(2i-1)2kM,i=1,2,…,m,
where m=2k-1M, we define the wavelet matrix Φm×m as:
(17)Φm×m=[Ψ(12m),Ψ(32m),…,Ψ(2m-12m)].
Indeed Φm×m has a diagonal form (see [15]).


### 2.4. The Operational Matrix of Fractional Integration

The fractional integration of orderα of the vector function Ψ(t) defined in (15) can be expressed as
(18)(IαΨ)(t)=PαΨ(t),
where Pα is the m×m operational matrix of fractional integration of order α. In [25] it is shown that the matrix Pα can be approximated as
(19)Pα≃Pm×mα=Φm×mPBαΦm×m-1,
where PBα is the operational matrix of fractional integration of order α of the Block Pulse functions (BPFs) that has the following form [26]:
(20)PBα=1mα1Γ(α+2)(1ξ1ξ2…ξm-101ξ1…ξm-2001…ξm-3000⋱⋮00001),
where ξi=(i+1)α+1-2iα+1+(i-1)α+1.In AppendixA, some important properties of BPFs are listed. In [15], it is shown that Pm×mα is an upper trigonometric matrix. Also, from (14) and (19) it is concluded that for a function f∈Cμ,μ≥-1, we have
(21)(Iαf)(t)≃CTΦm×mPBαBm(t).


## 2.1. Fractional Calculus

Here, we give some basic definitions and properties of the fractional calculus theory which will be used further in this paper.Definition 1.
A real functionf(t),t>0, is said to be in the space Cμ,μ∈ℝ, if there exists a real number p(>μ) such that f(t)=tpf1(t), where f1(t)∈C[0,∞] and it is said to be in the space Cμn if f(n)∈Cμ,n∈ℕ.Definition 2.
The Riemann-Liouville fractional integration operator of orderα≥0 of a function f∈Cμ,μ≥-1, is defined as
(5)(Iαf)(t)={1Γ(α)∫0t‍(t-τ)α-1f(τ)dτ,α>0,f(t),α=0.
It has the following properties:
(6)(i)IαIβ=Iα+β,(ii)IαIβ=IβIα,(iii)(IαIβf)(t)=(IβIαf)(t),(iv)Iαtϑ=Γ(ϑ+1)Γ(α+ϑ+1)tα+ϑ,
where α,β≥0,t>0, and ϑ>-1.
Riemann-Liouville fractional derivative operator of orderα>0 is defined as
(7)Dαf(t)=(ddt)nIn-αf(t),(n-1<α≤n),
where n is an integer and f∈C1n.
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 operatorD*α proposed by Caputo [24].Definition 3.
The fractional derivative operator of orderα>0 in the Caputo sense is defined as
(8)D*αf(t)=1Γ(n-α)∫0t‍(t-τ)n-α-1f(n)(τ)dτ,(n-1<α≤n),
where n is an integer, t>0, and f∈C1n.
The Riemann-Liouville and Caputo operators have a useful property as(9)IαD*αf(t)=f(t)-∑k=0n-1‍f(k)(0+)tkk!,(n-1<α≤n),
where n is an integer, t>0, and f∈C1n.
For more details on the mathematical properties of fractional derivatives and integrals see [24].


## 2.2. Chebyshev Wavelets

Chebyshev waveletsψpq(t)=ψ(k,p^,q,t) have four arguments; k∈ℕ,p=1,2,…,2k-1, and p^=2p-1; moreover q is the degree of the Chebyshev polynomial of the first kind and t is the normalized time, that is, t∈[0,1]. They are defined on the interval [0,1] as [25]
(10)ψpq(t)={2k/2T~q(2kt-p^),p^-12k≤t≤p^+12k,0,otherwise,
where
(11)T~q(t)={1π,q=02πTq(t),q>0,q=0,1,…,M-1, and M is a fixed positive integer. The coefficients in (11) are used for orthonormality of the system. Here, {Tq(t),q∈ℕ∪{0}} is the set of well-known Chebyshev polynomials of degree m which are orthogonal with respect to the weight function w(t)=1/1-t2 on the interval [-1,1] and satisfy the following recursive formula:
(12)T0(t)=1,T1(t)=t,Tq+1(t)=2tTq(t)-Tq-1(t),q=1,2,3….
We should note that in dealing with Chebyshev polynomials the weight function w~(t)=w(2t-1) has to be dilated and translated as wp(t)=w(2kt-p^), to get orthogonal wavelets.


## 2.3. Function Approximation

An arbitrary functionf(t)∈L2(R) defined over [0,1] may be expanded into Chebyshev wavelet basis as
(13)f(t)=∑p=1∞‍∑q=0∞‍cpqψpq(t),
where cpq=(f(t),ψpq(t)) in which (,) denotes the inner product.If the infinite series in (13) is truncated, then (13) can be written as
(14)f(t)≃∑p=12k-1‍∑q=0M-1‍cpqψpq(t)=CTΨ(t),
where C and Ψ(t) are 2k-1M×1 matrices given by
(15)C=[c10,c11,…,c1M-1,c20,…,c2M-1,…,c2k-10,…,c2k-1M-1]T,Ψ(t)=[ψ10(t),ψ11(t),…,ψ1M-1(t),ψ20(t),…,ψ2M-1(t),…,ψ2k-10(t),…,ψ2k-1M-1(t)]T.
Taking the collocation points
(16)ti=(2i-1)2kM,i=1,2,…,m,
where m=2k-1M, we define the wavelet matrix Φm×m as:
(17)Φm×m=[Ψ(12m),Ψ(32m),…,Ψ(2m-12m)].
Indeed Φm×m has a diagonal form (see [15]).


## 2.4. The Operational Matrix of Fractional Integration

The fractional integration of orderα of the vector function Ψ(t) defined in (15) can be expressed as
(18)(IαΨ)(t)=PαΨ(t),
where Pα is the m×m operational matrix of fractional integration of order α. In [25] it is shown that the matrix Pα can be approximated as
(19)Pα≃Pm×mα=Φm×mPBαΦm×m-1,
where PBα is the operational matrix of fractional integration of order α of the Block Pulse functions (BPFs) that has the following form [26]:
(20)PBα=1mα1Γ(α+2)(1ξ1ξ2…ξm-101ξ1…ξm-2001…ξm-3000⋱⋮00001),
where ξi=(i+1)α+1-2iα+1+(i-1)α+1.In AppendixA, some important properties of BPFs are listed. In [15], it is shown that Pm×mα is an upper trigonometric matrix. Also, from (14) and (19) it is concluded that for a function f∈Cμ,μ≥-1, we have
(21)(Iαf)(t)≃CTΦm×mPBαBm(t).


## 3. Existence and Uniqueness

For investigating existence of a unique solution for initial value problem (4), we reconsider this equation in the following operator form:
(22)𝒟*αy(x)=f(x)+p(x)ℱy(x)+𝒦𝒢y(x),
where we define
(23)𝒦y(x)=λ∫0x‍K(x,t)y(t)dt.
Applying operator ℐα, the inverse of 𝒟*α, on both sides of (22), and considering initial condition yield
(24)y(x)=c+ℐα(f(x)+p(x)ℱy(x)+𝒦𝒢y(x)).
Equation (24), can be written as a fixed point equation 𝒜y=y, where the map 𝒜 is defined as
(25)𝒜y(x)=c+ℐα(f(x)+p(x)ℱy(x)+𝒦𝒢y(x)),
Now let (C[J],∥·∥∞) be the Banach space of all continues functions on J=[0,b],b>1, with the norm ∥f∥∞=maxt∈J|f(t)|. Moreover, suppose that the nonlinear terms ℱ and 𝒢 satisfy Lipschitz conditions on [0,b] as
(26)|ℱy^m(x)-ℱy(x)|≤L1|y^m(x)-y(x)|,|𝒢y^m(x)-𝒢y(x)|≤L2|y^m(x)-y(x)|,
where L1 and L2 are fixed positive constants. Then we have the following.Theorem 4.
IfL1∥p∥∞+L2∥𝒦∥∞<Γ(α+1)/bα, then the initial value problem (4) has a unique solution y∈C[J].Proof.
There are different approaches in applying the Banach fixed point theorem to prove the existence of a unique solution for (4). We give a sketch of two approaches below. For this purpose we define the nonlinear integral operator
(27)𝒜:C[J]⟶C[J],𝒜y(x)≡c+1Γ(α)×∫0x‍(x-t)α-1(f(t)+p(t)ℱy(t)+𝒦𝒢y(t))dt.Approach 1.
Lety~,y∈C[J], then we have
(28)𝒜y~(x)-𝒜y(x)=1Γ(α)×∫0x‍(x-t)α-1×(p(t)[ℱy~(t)-ℱy(t)]+𝒦[𝒢y~(t)-𝒢y(t)])dt.
Then for x>0 we have
(29)|𝒜y~(x)-𝒜y(x)|≤1Γ(α)×∫0x‍(x-t)α-1(L1|p(t)||y~(t)-y(t)|+L2|𝒦||y~(t)-y(t)|)dt,|𝒜y~(x)-𝒜y(x)|≤(L1∥p∥∞+L2∥𝒦∥∞)×∥y~-y∥∞1Γ(α)×∫0x‍(x-t)α-1dt≤(L1∥p∥∞+L2∥𝒦∥∞)×∥y~-y∥∞bαΓ(α+1).
Therefore,
(30)∥𝒜y~-𝒜y∥∞≤ΩL1,L2,p,𝒦,α∥y~-y∥∞,
where ΩL1,L2,p,𝒦,α=(L1∥p∥∞+L2∥𝒦∥∞)(bα/Γ(α+1)).
SinceΩL1,L2,p,𝒦,α<1, by contraction mapping theorem, the initial value problem (4) has a unique solution in C[J].Approach 2.
Let us introduce the norm(31)∥|f(t)|∥=maxt∈Je-βt|f(t)|,
on the space C[J], which is equivalent to the standard norm ∥f∥∞ on C[J]. The parameter β is chosen to satisfy βα>Γ(α+1)/bα. Then for x>0 we have
(32)e-βx|𝒜y~(x)-𝒜y(x)|≤e-βxγ(α)∫0x‍(x-t)α-1eβte-βt×(|p(t)||ℱy~-ℱy|+|𝒦||𝒢y~-𝒢y|)dt≤(L1∥|p|∥+L2∥|𝒦|∥)∥|y~-y|∥×e-βxΓ(α)∫0x‍(x-t)α-1eβtdt=(L1∥|p|∥+L2∥|𝒦|∥)∥|y~-y|∥×e-βxeβxΓ(α)∫0x‍tα-1e-βtdt=(L1∥|p|∥+L2∥|𝒦|∥)∥|y~-y|∥×1Γ(α)1βα∫0βx‍tα-1e-tdt≤(L1∥|p|∥+L2∥|𝒦|∥)∥|y~-y|∥×1Γ(α)1βα∫0∞‍tα-1e-tdt=(L1∥|p|∥+L2∥|𝒦|∥)∥|y~-y|∥1βα.
Then,
(33)∥|𝒜y~-𝒜y|∥≤ΩL1,L2,p,𝒦,α∥|y~-y|∥,
where ΩL1,L2,p,𝒦,α=(L1∥|p|∥+L2∥|𝒦|∥)(bα/Γ(α+1)).
Therefore, the operator𝒜 is a contraction on the Banach space (C[J],∥|·|∥), and so the initial value problem (4) has a unique solution in C[J].


## 4. Chebyshev Wavelets Method

Here, the Chebyshev wavelets expansion together with their operational matrices of fractional order integration are used for numerical solution of nonlinear NFIDE (4). For solving this problem we assume that
(34)D*αy(x)=KmTΨm(x),
where KmT is an unknown vector and Ψm(x) is the expansion of (15) by BPFs, that is, Ψm=Φm×mBm(x), and Bm(x) is defined in Appendix A. By using initial conditions and (9), we have
(35)y(x)=KmTPm×mαΨm(x)+c.
From (35) we have
(36)y(x)=[KmTPm×mαΦm×m+[c,c,…,c]]Bm(x).
This equation can also be written as
(37)y(x)=[KmTPm×mαΦm×m+[c,c,…,c]]×Bm(x)=[a1,a2,…,am]Bm(x),
and from Appendix B, we have
(38)F(y(x))=[F(a1),F(a2),…,F(am)]Bm(x),∫0x‍K(x,t)G(y(t))dt=[G(a1),G(a2),…,G(am)]×∫0x‍K(x,t)Bm(t)dt,
where
(39)∫0x‍K(x,t)Bm(t)dt=[∫0x‍K(x,t)b0(t)dt,∫0x‍K(x,t)b1(t)dt,…,∫0x‍K(x,t)bm-1(t)dt]T.
By substituting (34) and (38) into (4), we obtain
(40)KmTΦm×mBm(x)=f(x)+p(x)[F(a1),F(a2),…,F(am)]Bm(x)+[G(a1),…,G(am)]∫0x‍K(x,t)Bm(t)dt.
Now, from (37) and (40) we have
(41)[a1-c,a2-c,…,am-c](PBα)-1Bm(x)=f(x)+p(x)[F(a1),F(a2),…,F(am)]Bm(x)+[G(a1),…,G(am)]∫0x‍K(x,t)Bm(t)dt.
This is a nonlinear algebraic equation. Now, by taking collocation points, expressed in (16), this equation is transformed into a nonlinear system of algebraic equations with m unknowns ai(i=1,2,…,m). By solving this system and determining ai, we get the numerical solution of problem (4).


### 4.1. Error Investigation of the Chebyshev Wavelets Method

In this section we investigate error of the Chebyshev wavelets method. The representation error can be obtained when a functionf(t) is presented on interval [0,1].Theorem 5 (see [27]).
A functionf(t) defined on [0,1) with bounded second derivative, say |f′′(t)|≤M^, can be expanded as an infinite sum of Chebyshev wavelets, and the series converges uniformly to f(t), that is
(42)f(t)=∑p=1∞‍∑q=0∞‍cpqψpq(t),
where cpq’s are defined in (13).Theorem 6 (see [27]).
Letf(t) be a continuous function defined on [0,1), with bounded second derivative |f′′(t)|≤M^, and let ∑p=12k-1‍∑q=0M-1‍cpqψpq(t) be the approximate solution using Chebyshev wavelets method, respectively. Then for the error bound one has:
(43)σk,M<πM^8(∑p=2k-1+1∞‍1p5∑q=M∞‍1(q-1)4)1/2,
where
(44)σk,M=(∫01‍(f(t)-∑p=12k-1‍∑q=0M-1‍cpqψpq(t))2wp(t)dt)1/2.Lemma 7.
Suppose thatf(t) is approximated on [0,1) by
(45)fm(t)=∑p=12k-1‍∑q=0M-1‍cpqψpq(t),
and moreover suppose that by solving some problems one has found c^pq as an approximation of cpq and put:
(46)f^m(t)=∑p=12k-1‍∑q=0M-1‍c^pqψpq(t).
Then for each t∈[0,1) one has
(47)∥f^m(t)-fm(t)∥≤m∥f^m(t)-f(t)∥∞.Proof .
We have(48)∥f^m(t)-fm(t)∥=(∫01‍[∑p=12k-1‍∑q=0M-1‍c^pqψpq(t)-∑p=12k-1‍∑q=0M-1‍cpqψpq(t)]2×wq(t)dt∫01‍[∑p=12k-1‍∑q=0M-1‍c^pqψpq(t)-∑p=12k-1‍∑q=0M-1‍cpqψpq(t)]2)1/2=(∫01‍∑p=12k-1‍∑q=0M-1‍[c^pq-cpq]2[ψpq(t)]2wq(t)dt)1/2=∑p=12k-1‍∑q=0M-1‍|c^pq-cpq|(∫01‍[ψpq(t)]2wq(t)dt)1/2=∑p=12k-1‍∑q=0M-1‍|c^pq-cpq|≤∑p=12k-1‍∑q=0M-1‍∥f^m(t)-f(t)∥∞=m∥f^m(t)-f(t)∥∞.
This completes the proof.Corollary 8.
Suppose that by solving some problems one obtainsf^m(t) as the approximation of f(t). Then one has:
(49)∥f^m(t)-f(t)∥≤σk,M+m∥f^m(t)-f(t)∥∞.Proof .
For everyt∈[0,1) we have
(50)∥f^m(t)-f(t)∥=∥f^m(t)-f(t)-fm(t)+fm(t)∥≤∥fm(t)-f(t)∥+∥f^m(t)-fm(t)∥.
Now from Theorem 6 and Lemma 7 we have
(51)∥f^m(t)-f(t)∥≤σk,M+m∥f^m(t)-f(t)∥∞.
This completes the proof.In real problems, we often tend to solve equations with unknown exact solutions. Hence, when we apply our method to these problems, we cannot say that this approximate solution is good or bad unless we are able to calculate the error functione^m(x)=f^m(x)-f(x). Then it is necessary to introduce a process for estimating the error function when the exact solution is unknown [28]. Here, we introduce a method to estimate the error function. Let y^m(x) be an approximate solution of (4). Then from (24), it is concluded that y^m(x) satisfies the following equation:
(52)y^m(x)=c+1Γ(α)∫0x‍(x-t)α-1f(t)dt+1Γ(α)∫0x‍(x-t)α-1p(t)F(y^m(t))dt+λΓ(α+1)∫0x‍(x-t)αK(x,t)G(y^m(t))dt+rm(x),
where rm(x) is the perturbation function that depends only on y^m(x) and can be obtained by substituting the computed solution y^m(x) into the equation
(53)rm(x)=y^m(x)-c-1Γ(α)∫0x‍(x-t)α-1f(t)dt-1Γ(α)∫0x‍(x-t)α-1p(t)F(y^m(t))dt-λΓ(α+1)∫0x‍(x-t)α-1K(x,t)G(y^m(t))dt.
Subtracting (53) from (24), we obtain
(54)e^m(x)=1Γ(α)∫0x‍(x-t)α-1p(t)[F(y^m(t))-F(y(t))]+λΓ(α+1)∫0x‍(x-t)αK(x,t)×[G(y^m(t))-G(y(t))]dt+rm(x),
in which we define e^m(x)=y^m(x)-y(x) as the error function. Since F and G are analytic functions, then we can write
(55)F(y^m(x))-F(y(x))=F(y^m(x))-F(y^m(x)-e^m(x))=F′(y^m(x))e^m(x)-12F′′(wm(x))[e^m(x)]2,
and similarly
(56)G(y^m(t))-G(y(t))=G′(y^m(t))e^m(t)-12G′′vm(t)[e^m(t)]2,
where, by Taylor’s Theorem, wm(t)=y^m(t)-θe^m(t),θ∈(0,1), and vm(t)=y^m(t)-δe^m(t),δ∈(0,1). Thus, by substituting (55), and (56) in the error equation (54), we get a nonlinear fractional integral equation in which the error function e^m(t) is unknown. Obviously, we can apply the proposed method for this equation to find an approximation of the error function e^m(x).


## 4.1. Error Investigation of the Chebyshev Wavelets Method

In this section we investigate error of the Chebyshev wavelets method. The representation error can be obtained when a functionf(t) is presented on interval [0,1].Theorem 5 (see [27]).
A functionf(t) defined on [0,1) with bounded second derivative, say |f′′(t)|≤M^, can be expanded as an infinite sum of Chebyshev wavelets, and the series converges uniformly to f(t), that is
(42)f(t)=∑p=1∞‍∑q=0∞‍cpqψpq(t),
where cpq’s are defined in (13).Theorem 6 (see [27]).
Letf(t) be a continuous function defined on [0,1), with bounded second derivative |f′′(t)|≤M^, and let ∑p=12k-1‍∑q=0M-1‍cpqψpq(t) be the approximate solution using Chebyshev wavelets method, respectively. Then for the error bound one has:
(43)σk,M<πM^8(∑p=2k-1+1∞‍1p5∑q=M∞‍1(q-1)4)1/2,
where
(44)σk,M=(∫01‍(f(t)-∑p=12k-1‍∑q=0M-1‍cpqψpq(t))2wp(t)dt)1/2.Lemma 7.
Suppose thatf(t) is approximated on [0,1) by
(45)fm(t)=∑p=12k-1‍∑q=0M-1‍cpqψpq(t),
and moreover suppose that by solving some problems one has found c^pq as an approximation of cpq and put:
(46)f^m(t)=∑p=12k-1‍∑q=0M-1‍c^pqψpq(t).
Then for each t∈[0,1) one has
(47)∥f^m(t)-fm(t)∥≤m∥f^m(t)-f(t)∥∞.Proof .
We have(48)∥f^m(t)-fm(t)∥=(∫01‍[∑p=12k-1‍∑q=0M-1‍c^pqψpq(t)-∑p=12k-1‍∑q=0M-1‍cpqψpq(t)]2×wq(t)dt∫01‍[∑p=12k-1‍∑q=0M-1‍c^pqψpq(t)-∑p=12k-1‍∑q=0M-1‍cpqψpq(t)]2)1/2=(∫01‍∑p=12k-1‍∑q=0M-1‍[c^pq-cpq]2[ψpq(t)]2wq(t)dt)1/2=∑p=12k-1‍∑q=0M-1‍|c^pq-cpq|(∫01‍[ψpq(t)]2wq(t)dt)1/2=∑p=12k-1‍∑q=0M-1‍|c^pq-cpq|≤∑p=12k-1‍∑q=0M-1‍∥f^m(t)-f(t)∥∞=m∥f^m(t)-f(t)∥∞.
This completes the proof.Corollary 8.
Suppose that by solving some problems one obtainsf^m(t) as the approximation of f(t). Then one has:
(49)∥f^m(t)-f(t)∥≤σk,M+m∥f^m(t)-f(t)∥∞.Proof .
For everyt∈[0,1) we have
(50)∥f^m(t)-f(t)∥=∥f^m(t)-f(t)-fm(t)+fm(t)∥≤∥fm(t)-f(t)∥+∥f^m(t)-fm(t)∥.
Now from Theorem 6 and Lemma 7 we have
(51)∥f^m(t)-f(t)∥≤σk,M+m∥f^m(t)-f(t)∥∞.
This completes the proof.In real problems, we often tend to solve equations with unknown exact solutions. Hence, when we apply our method to these problems, we cannot say that this approximate solution is good or bad unless we are able to calculate the error functione^m(x)=f^m(x)-f(x). Then it is necessary to introduce a process for estimating the error function when the exact solution is unknown [28]. Here, we introduce a method to estimate the error function. Let y^m(x) be an approximate solution of (4). Then from (24), it is concluded that y^m(x) satisfies the following equation:
(52)y^m(x)=c+1Γ(α)∫0x‍(x-t)α-1f(t)dt+1Γ(α)∫0x‍(x-t)α-1p(t)F(y^m(t))dt+λΓ(α+1)∫0x‍(x-t)αK(x,t)G(y^m(t))dt+rm(x),
where rm(x) is the perturbation function that depends only on y^m(x) and can be obtained by substituting the computed solution y^m(x) into the equation
(53)rm(x)=y^m(x)-c-1Γ(α)∫0x‍(x-t)α-1f(t)dt-1Γ(α)∫0x‍(x-t)α-1p(t)F(y^m(t))dt-λΓ(α+1)∫0x‍(x-t)α-1K(x,t)G(y^m(t))dt.
Subtracting (53) from (24), we obtain
(54)e^m(x)=1Γ(α)∫0x‍(x-t)α-1p(t)[F(y^m(t))-F(y(t))]+λΓ(α+1)∫0x‍(x-t)αK(x,t)×[G(y^m(t))-G(y(t))]dt+rm(x),
in which we define e^m(x)=y^m(x)-y(x) as the error function. Since F and G are analytic functions, then we can write
(55)F(y^m(x))-F(y(x))=F(y^m(x))-F(y^m(x)-e^m(x))=F′(y^m(x))e^m(x)-12F′′(wm(x))[e^m(x)]2,
and similarly
(56)G(y^m(t))-G(y(t))=G′(y^m(t))e^m(t)-12G′′vm(t)[e^m(t)]2,
where, by Taylor’s Theorem, wm(t)=y^m(t)-θe^m(t),θ∈(0,1), and vm(t)=y^m(t)-δe^m(t),δ∈(0,1). Thus, by substituting (55), and (56) in the error equation (54), we get a nonlinear fractional integral equation in which the error function e^m(t) is unknown. Obviously, we can apply the proposed method for this equation to find an approximation of the error function e^m(x).


## 5. Numerical Examples

In this section, for the sake of comparison, we have selected some examples that their exact solutions exist, which will ultimately show the simplicity, effectiveness, and exactness of the proposed method. All programs are performed by Maple 15 and digits 20. We consider the following test problems.Example 1.
Let us first consider the following NFIDE:(57)D*αy(x)=f(x)+y(x)+∫0x‍xtarcsin(y(t))dty(0)=0x∈[0,2],
where
(58)f(x)=cos(x)-sin(x)-13x4.
The exact solution of this problem for α=1 is y(x)=sin(x). Figure 1 shows the behavior of the numerical solution for m=128(M=4,k=6). From Figure 1, it is seen that the numerical solution is in a very good agreement with the exact solution for α=1. Therefore, we hold that the solutions for α=0.90 and α=0.95 are also credible. Figure 1 also shows that as α→1, the approximate solution tends to y(x)=sin(x), which is the exact solution of this equation in the case α=1.Figure 1
Numerical solution (N.S.) and exact solution (Ex.S).Example 2.
Consider the following NFIDE:(59)D*αy(x)=f(x)+ey(x)+∫0x‍t[y(t)]2dt,y(0)=0,x∈[0,5],
where
(60)f(x)=6x3-αΓ(4-α)-15πx3/2-α2Γ(5/2-α)-20x5+4013x13/2-18x8+ex3-10x3/2.
The exact solution of this problem is y(x)=x3-10x3/2.This NFIDE is now solved by the Chebyshev wavelet method forα=0.75. Numerical solutions of this problem for some values of m are shown in Figure 2. It is obviously seen that computing y0(x) via the standard ADM and ITM is very difficult because analytic computation of ∫0x‍et3-10t3/2dt is impossible. Also computing other components yi(x)(i≥1), via the standard ADM, VIM, and ITM, is very difficult in each iteration and also requires a large amount of computational work. We can avoid such integrations by taking the truncated Maclaurin expansion for the exponential term in (59), that is, ex=1+x+(1/2!)x2+⋯+(1/5!)x5. In order to show efficiency of the proposed method, a comparison between the numerical solution given by the proposed method for m=128, y4(x) of the ADM, ITM, and VIM, and the exact solution is performed in Figure 3. From Figure 3, it can be concluded that the Chebyshev wavelet method for the numerical solution of this problem is very efficient and accurate in comparison with other presented numerical methods.Numerical solutions of Example1, for some values of m.
(a)
(b)
(c)Figure 3
Numerical solutions of Example1 for α=0.75.Example 3.
Consider the following NFIDE:(61)D*αy(x)=f(x)-[y(x)]2+x∫0x‍[y(t)]3dt,y(0)=1,x∈[0,4.5],
where
(62)f(x)=-2x2-αΓ(3-α)+3x1-αΓ(2-α)-17x8+32x7-245x6+94x5+9x4+8x2+6x+1.
The exact solution of this problem is y(x)=1+3x-x2.This problem is now solved by the Chebyshev wavelet method forα=0.85. Numerical solutions of this problem for some values of m are shown in Figure 4. A comparison between the numerical solution given by the proposed method for m=128, y4(x) of the ADM and ITM and y3(x) of the VIM and the exact solution is performed in Figure 5. From Figure 5, it can be concluded that the Chebyshev wavelet method for the numerical solution of this problem is very efficient and accurate in comparison with other presented numerical methods. Moreover, the solution has been derived in a large domain. It can be mentioned that computing the components yi(x) via the ADM, VIM and ITM require a large amount of computational work and memory, such that computing other components for i>4 is very difficult.Numerical solutions of Example2, for some values of m.
(a)
(b)
(c)Figure 5
Numerical solutions of Example2 for α=0.85.Example 4.
Finally, consider the following NFIDE:(63)D*αy(x)=f(x)+∫0x‍e-t[y(t)]2dt,y(0)=0,x∈[0,5],
where
(64)f(x)=6x3-αΓ(4-α)-12x2-αΓ(3-α)+8x1-αΓ(2-α)+e-x×(80+80x+40x2-8x3+22x4-6x5+x6).
The exact solution of this problem is y(x)=x3-6x2+8x.This NFIDE is now solved by the proposed method forα=0.80. Numerical solutiones of this problem for some selected values of m are shown in Figure 6. It is clearly seen that computing the components yi(x) via the standard ADM, VIM, and ITM is very difficult and also requires a large amount of computational work. We can avoid such fractional integrations by taking the truncated Maclaurin expansion for the exponential term in (63). A comparison between the numerical solution given by the proposed method for m=128, y4(x) of the ADM, ITM, and VIM for α=0.80 and α=0.90 and the exact solution is performed in Figure 7. From Figure 7, it can be concluded that the Chebyshev wavelet method for the numerical solution of this problem is very efficient and accurate in comparison with the ADM, VIM, and ITM.Numerical solutions of Example3, for some values of m.
(a)
(b)
(c)Numerical solutions of Example3, for α=0.80 (a) and α=0.90 (b).
(a)
(b)


## 6. Discussion and Conclusions

In this paper, we proposed the Chebyshev wavelets method for the numerical solution of nonlinear fractional integro-differential equations in a large interval. Existence of a unique solution for such equations is proved. Convergence and error analysis of the Chebyshev wavelets expansion is discussed. Efficient approximate solutions have derived for NFIDEs, and the results have been shown remarkable performance. Moreover the new approach has been compared with the ADM, VIM, and ITM in solving various nonlinear fractional integro-differential equations. There are five important points to be noted here.(1)
The proposed method uses a new technique for computation of nonlinear terms in such equations.(2)
The proposed approach can provide the suitable approximate solution in a large interval by using only a few number of collocation points, as shown in the solved examples.(3)
The proposed method is very simple and requires less computational work in comparison with ADM, VIM, and ITM.(4)
The proposed method overcomes the probable difficulties arising in calculating integrals.(5)
The proposed method does not need to use Adomian polynomials and also has no need for the Lagrange multiplier, correction functional, stationary conditions, the variational theory, and so forth, which eliminates the complications that exist in VIM.So in comparison with the ADM, VIM, and ITM, the proposed method is proved to be simpler in principle and more convenient for computer algorithms.



---

*Source: 482083-2013-10-31.xml*