# Toward the Approximate Solution for Fractional Order Nonlinear Mixed Derivative and Nonlocal Boundary Value Problems

Authors: Hammad Khalil; Mohammed Al-Smadi; Khaled Moaddy; Rahmat Ali Khan; Ishak Hashim
Journal: Discrete Dynamics in Nature and Society
 (2016)

---

## Abstract

The paper is devoted to the study of operational matrix method for approximating solution for nonlinear coupled system fractional differential equations. The main aim of this paper is to approximate solution for the problem under two different types of boundary conditions,m^-point nonlocal boundary conditions and mixed derivative boundary conditions. We develop some new operational matrices. These matrices are used along with some previously derived results to convert the problem under consideration into a system of easily solvable matrix equations. The convergence of the developed scheme is studied analytically and is conformed by solving some test problems.


---

## Body

## 1. Introduction

Fractional calculus is generalization of integer order integration and differentiation to its noninteger (fractional) order counterpart. Fractional calculus has proved to be a valuable tool in modeling many natural phenomena of physics, chemistry, engineering, aerodynamics, electrodynamics of complex medium, polymer rheology, and so forth [1]. It is well known that fractional order operator is a nonlocal operator. Unlike integer order models the fractional order models have the property that the current state of the system depends on all its historic states. This makes fractional models interesting and remains as an active and hot area of research.Differential equations are used to model different physical and engineering phenomena. The aim is to know about the future state of the system under consideration. It is some time necessary to model differential equations with different kinds of initial and boundary conditions (depends on the nature of the problem). Various types of conditions which are used as boundary conditions are multipoint local boundary conditions, integral type boundary conditions, multipoint nonlocal boundary conditions, and mixed derivative boundary conditions. The multipoint boundary value problems appear in wave propagation and in elastic stability. For example, the vibrations of a guy wire of a uniform cross section composed ofm sections of different densities can be molded as a multipoint boundary value problem (see [2, 3] and the references therein).In this paper, we are interested in finding approximate solutions for the following generalized class of nonlinear coupled systems of fractional order differential equations:(1)uσ1=ft,uγ1,vγ1,u,v,vσ2=gt,uγ1,vγ1,u,v,where u=u(t) and v=v(t) are the unknown solutions to be determined. These solutions are defined on finite time domain; that is, t∈[0,τ]. The superscripts represent the order of derivatives defined in Caputo sense and are defined as 1<σ1, σ2≤2, and 0<γ1≤1. f and g are nonlinear functional of u and v and their fractional derivatives. The method developed in this paper is designed for approximating u and v under any of the following two types of boundary conditions. (i)
M-point nonlocal boundary conditionsare as follows:(2)u0=u0,uτ=∑i=1m-2ζiuξi,v0=v0,vτ=∑i=1m-2ρivξi,whereu0 and v0 are some real constants and ξi are intermediate points of the domain defined as (3)0<ξ1<ξ2⋯<ξm-2<τ.(ii)
Mixed derivative boundary conditionsare as follows:(4)u′0=a1u0+a2uτ,u′τ=b1u0+b2uτ,v′0=c1v0+c2vτ,v′τ=d1v0+d2vτ,whereai, bi, ci, and di are some real constants.We develop new operational matrices based on Jacobi polynomials and use these operational matrices to develop new formulation for finding approximate solutions for the problems. The technique basically lies in the domain of spectral methods. Spectral methods are widely used for solutions to differential equations, functions approximations, and variational problems. These methods demand approximation of solutions for a problem by truncated series of smooth global functions and provide very accurate approximations for a smooth solution with relatively few degrees of freedom. The main reason behind this accuracy is the behavior of spectral coefficientsan, which tend to zero faster than any algebraic power of their index n.Among others, some of the well known mathematicians who successfully applied spectral method are Dehghan [4–9], Bhrawy [10, 11], Doha et al. [12], and Saadatmandi [7, 13, 14]. In these papers the authors solved many scientific problems using spectral methods. Also Dehghan and Shakeri [8] used a seminumerical technique for solving the multipoint boundary value problems. In [4, 6] an efficient way is developed for the construction of operational matrices. In [5] the authors used operational matrices of Bernstein polynomials to solve nonlinear age-structured model. We are motivated by the work of Bhrawy and Al-Shomrani [11], who solved fractional order differential equations (both linear and nonlinear) with m-point boundary conditions (local) via shifted Legendre polynomials and collocation technique.To solve the nonlinear coupled system of boundary value problem we implement operational matrix method combined with quasilinearization method. Quasilinearization method was first introduced by Bellman and Kalaba [15] to solve nonlinear ordinary or partial differential equations as a generalization of the Newton-Raphson method. The origin of this method lies in the theory of dynamic programming. In this method, the nonlinear equations are expressed as a sequence of linear equations and these equations are solved recursively. The main advantage of this method is that it converges monotonically and quadratically to the exact solution to the original equations [16]. Also some other interesting works in which qasilinearization method is used for scientific problems are available in [17–19].In our previous work, we have constructed some efficient methods for the numerical simulations of couple systems of linear fractional differential equations and fractional order partial differential equations [20–22]. Local and nonlocal boundary value problems can be found in [23–25]. For the readers who are new to the field, we recommend studying our previous work in order to get a better understanding.


## 2. Preliminaries

In this section, we recall some basic definitions and known results from fractional calculus. More details can be found in [1].Definition 1.
Given an interval[a,b]⊂R, the Riemann-Liouville fractional order integral of order α∈R+ of a function ϕ∈(L1[a,b],R) is defined by(5)aItαϕt=1Γα∫att-sα-1ϕsds,provided the integral on right hand side exists.Definition 2.
For a given functionϕ(t)∈Cn[a,b], the Caputo fractional order derivative of order α is defined as (6)Dαϕt=1Γn-α∫atϕnst-sα+1-nds,n-1≤α<n,n∈N,provided the right side is pointwise defined on (a,∞), where n=α+1 in case α is not an integer and n=α in case α is an integer.Hence, it follows thatIαtk=(Γ1+k/Γ1+k+α)tk+αforα>0, k≥0, DαC=0, for a constant C, and(7)Dαtk=Γ1+kΓ1+k-αtk-α,fork≥α.


### 2.1. The Shifted Jacobi Polynomials

The well known two parametric Jacobi polynomials defined on[0,τ], with parameters α and β, are given by the following relation:(8)Pη,iα,βt=∑k=0i℧i,kα,βtk,i=0,1,2,3,…,where(9)℧i,kα,β=-1i-kΓi+β+1Γi+k+α+β+1Γk+β+1Γi+α+β+1i-k!k!τk.The Jacobi polynomials are orthogonal and the orthogonality condition is(10)∫0τPτ,iα,βtPτ,jα,βtwτα,βtdt=Rτ,jα,βδi,j,where wτ(α,β)(t)=(τ-t)αtβ is the weight function and Rτ,j(α,β)=τα+β+1Γ(j+α+1)Γ(j+β+1)/(2j+α+β+1)Γ(j+1)Γ(j+α+β+1) is known as normalization constant.Any squared integrable functionv(t) on [0,τ] can be approximated by shifted Jacobi polynomials as follows: v(t)=∑j=0∞cjPη,j(α,β)(t), where the coefficient cj can easily be calculated using (10). In practice we are interested in the truncated series, which can also be written in vector form as(11)vt=HMTΨMt,where M=m+1, HMT is the coefficient vector, and ΨM(t) is M terms function vector.One important property of shifted Jacobi polynomials which will be frequently used in our analysis is given as(12)maxx∈0,τ⁡Pτ,iα,βx≤C^i,ϵ,where C^(i,ϵ)=Γ(i+ϵ+1)/Γ(i+1)Γ(ϵ+1) and ϵ=max⁡(α,β).Lemma 3.
The definite integral of the product of weight function with three Jacobi polynomials over the domain[0,τ] is constant, defined as ϝ(l,m,n)(i,j,k), defined by(13)∫0τwτα,βtPτ,iα,βtPτ,jα,βtPτ,kα,βtdt=ϝl,m,ni,j,k,where(14)ϝl,m,ni,j,k=∑l=0i∑m=0j∑n=0k℧i,lα,β℧j,mα,β℧k,nα,βΥl,m,n,Υl,m,n=Γl+m+n+β+1Γα+1τl+m+n+α+β+1Γl+m+n+α+β+1.Proof.
In view of (8), we obtain(15)∫0τwτα,βtPτ,iα,βtPτ,jα,βtPτ,kα,βtdt=∑l=0i℧i,lα,β∑m=0j℧j,mα,β∑n=0k℧k,nα,β∫0τtl+m+n+βτ-tαdt.Using convolution theorem of Laplace transform, we have (16)∫0τtl+m+n+βτ-tαdt=Γl+m+n+β+1Γα+1τl+m+n+α+β+1Γl+m+n+α+β+1=Υl,m,n.Using the value of Υ(l,m,n) in (15), we obtain (17)∫0τwτα,βtPτ,iα,βtPτ,jα,βtPτ,kα,βtdt=ϝl,m,ni,j,k.


## 2.1. The Shifted Jacobi Polynomials

The well known two parametric Jacobi polynomials defined on[0,τ], with parameters α and β, are given by the following relation:(8)Pη,iα,βt=∑k=0i℧i,kα,βtk,i=0,1,2,3,…,where(9)℧i,kα,β=-1i-kΓi+β+1Γi+k+α+β+1Γk+β+1Γi+α+β+1i-k!k!τk.The Jacobi polynomials are orthogonal and the orthogonality condition is(10)∫0τPτ,iα,βtPτ,jα,βtwτα,βtdt=Rτ,jα,βδi,j,where wτ(α,β)(t)=(τ-t)αtβ is the weight function and Rτ,j(α,β)=τα+β+1Γ(j+α+1)Γ(j+β+1)/(2j+α+β+1)Γ(j+1)Γ(j+α+β+1) is known as normalization constant.Any squared integrable functionv(t) on [0,τ] can be approximated by shifted Jacobi polynomials as follows: v(t)=∑j=0∞cjPη,j(α,β)(t), where the coefficient cj can easily be calculated using (10). In practice we are interested in the truncated series, which can also be written in vector form as(11)vt=HMTΨMt,where M=m+1, HMT is the coefficient vector, and ΨM(t) is M terms function vector.One important property of shifted Jacobi polynomials which will be frequently used in our analysis is given as(12)maxx∈0,τ⁡Pτ,iα,βx≤C^i,ϵ,where C^(i,ϵ)=Γ(i+ϵ+1)/Γ(i+1)Γ(ϵ+1) and ϵ=max⁡(α,β).Lemma 3.
The definite integral of the product of weight function with three Jacobi polynomials over the domain[0,τ] is constant, defined as ϝ(l,m,n)(i,j,k), defined by(13)∫0τwτα,βtPτ,iα,βtPτ,jα,βtPτ,kα,βtdt=ϝl,m,ni,j,k,where(14)ϝl,m,ni,j,k=∑l=0i∑m=0j∑n=0k℧i,lα,β℧j,mα,β℧k,nα,βΥl,m,n,Υl,m,n=Γl+m+n+β+1Γα+1τl+m+n+α+β+1Γl+m+n+α+β+1.Proof.
In view of (8), we obtain(15)∫0τwτα,βtPτ,iα,βtPτ,jα,βtPτ,kα,βtdt=∑l=0i℧i,lα,β∑m=0j℧j,mα,β∑n=0k℧k,nα,β∫0τtl+m+n+βτ-tαdt.Using convolution theorem of Laplace transform, we have (16)∫0τtl+m+n+βτ-tαdt=Γl+m+n+β+1Γα+1τl+m+n+α+β+1Γl+m+n+α+β+1=Υl,m,n.Using the value of Υ(l,m,n) in (15), we obtain (17)∫0τwτα,βtPτ,iα,βtPτ,jα,βtPτ,kα,βtdt=ϝl,m,ni,j,k.


## 3. Operational Matrices

In this section we first recall some previous results and then derive some new operational matrices. The following results are known [22].Lemma 4.
LetΨM(t) be the function vector as defined in (11); then the γ-order integration of ΨM(t) is given by IγΨMt=HM×Mτ,γΨMt, where HM×Mτ,γ is the operational matrix of integration of order γ. The entries Θi,j,k of HM×Mτ,γ are defined as(18)Θi,j,k=∑k=0iΛi,k,γSj,where(19)Λi,k,γ=-1i-kΓi+β+1Γi+k+α+β+1Γ1+kΓk+β+1Γi+α+β+1i-k!k!Γ1+k+γτk,Sj=∑l=0j-1j-l2j+α+β+1Γj+1Γj+l+α+β+1Γk+γ+l+β+1Γα+1τγΓj+α+1Γl+β+1j-l!l!Γk+γ+l+β+α+2.Proof.
The detailed  proof of this lemma is available in [22].Lemma 5.
LetΨM(t) be the function vector as defined in (11); then the γ-order derivative of ΨM(t) is given by DγΨM(t)=GM×Mτ,γΨM(t), where GM×Mτ,γ is the operational matrix of fractional differentiation of order γ. The entries Φi,j,k of GM×Mτ,γ are defined as(20)Φi,j,k=∑k=γiΛi,k,γSj,where(21)Λi,k,γ=-1i-kΓi+β+1Γi+k+α+β+1Γ1+kΓk+β+1Γi+α+β+1i-k!k!Γ1+k-γτk,Sj=∑l=0j-1j-l2j+α+β+1Γj+1Γj+l+α+β+1Γk+γ+l+β+1Γα+1τγΓj+α+1Γl+β+1j-l!l!Γk+γ+l+β+α-2.Proof.
The detailed proof of this lemma is available in [22].In [22], these operational matrices are used for solutions to coupled systems of initial value problems for fractional order differential equations. These matrices do not have the ability to handle m-point nonlocal boundary value problems. Hence the need to develop new operational matrices which have the ability to handle boundary conditions as well as initial conditions has been of great importance. The new operational matrices include the operational matrices studied in [22] as a factor and have ability to solve fractional order boundary value problems.Lemma 6.
Letu(t) and ϕn(t) be any functions defined on [0,τ]. Then(22)ϕntDσut=WMTKϕnσΨMt,where WMT is the Jacobi coefficients vector of u(t) defined by (11) and(23)Kϕnσ=GM×Mτ,σJM×Mτ,ϕn.The matrix GM×Mτ,σ is the operational matrix of derivative as defined in Lemma 5 and(24)JM×Mη,ϕn=Θr+1,s+1,r=0,1,…,m,s=0,1,…,m.The entries are defined by the relation(25)Θr,s=1Rτ,sα,β∑i=0mciϝl,m,ni,r,s,where ci are the Jacobi coefficients of ϕn(t) and ϝ(l,m,n)(i,r,s) is as defined in Lemma 3.Proof.
Using (11) and Lemma 5, we obtain(26)ϕntDσut=ϕntWMTGM×Mτ,σΨMt,where u(t)=WMTΨM(t). With rearranging, the above equation can be written as(27)ϕntDσYt=WMTGM×Mτ,σΨMt︷,where(28)ΨMt︷=ϕntPτ,0α,βtϕntPτ,1α,βt⋯ϕntPτ,rα,βt⋯ϕntPτ,mα,βtT.Approximating ϕn(t) with Jacobi polynomials, we have(29)ϕnt=∑i=0mCiPτ,iα,βt.Hence (28) takes the form(30)ΨMt︷=ℵ0tℵ1t⋯ℵrt⋯ℵmtT,where(31)ℵrt=∑i=0mCiPτ,iα,βtPτ,rα,βt,r=0,1,…,m.Now approximating ℵr(t) with Jacobi polynomials, we have(32)ℵrt=∑s=0mdsrPτ,sα,βt,where(33)dsr=1Rτ,sα,β∫0τℵrtWτα,βtPτ,sα,βtdt.Using (31) in (33) we obtain(34)dsr=1Rτ,sα,β∑i=0mci∫0τWτα,βtPτ,iα,βtPτ,rα,βtPτ,sα,βtdt,which in view of Lemma 3 takes the form(35)dsr=1Rτ,sα,β∑i=0mciϝl,m,ni,r,s=Θr,s.Repeating the procedure for r=0,1,…,m and s=0,1,…,m, we obtain(36)ΨMt︷=JM×Mη,ϕnΨMt.Using (36) in (27), we get(37)ϕntDσut=WMTGM×Mτ,σJM×Mτ,ϕnΨMt=WMTKϕnσΨMt.Lemma 7.
Forϕnc=ctn, where c and n are real constants, u(t)=WMTΨM(t), and, for 0<η≤τ, the following holds:(38)ϕnc0Iησut=WMTQM×Mσ,ϕnc,ηΨMt,where QM×M(σ,ϕnc,η) is M×M operational matrix and is defined by(39)Qσ,ϕnc,η=di+1,j+1,i=0,1,…,m,j=0,1,…,m.
And the entriesd(i,j) are defined by the relation(40)di,j=∑k=0i℧i,kα,β∑l=0j℧j,lα,βcΓ1+kΓl+β+n+1Γα+1Rη,jα,βΓ1+k+σΓl+β+n+α+1τl+β+n+α+1ηk+σ.Proof.
Using the definition of fractional integral (5), on (41)ut=∑i=0mciPτ,iα,βt,we obtain(42)0Iησut=∑i=0mci0IησPτ,iα,βt,which in view of the definition of Jacobi polynomials yields(43)0Iησut=∑i=0mci∑k=0i℧i,kα,β0Iησtk.Evaluating the integral and after simplification, we obtain(44)0Iησut=∑i=0mci∑k=0i℧i,kα,βΓ1+kηk+σΓ1+k+σ.For simplicity, use the following notation:(45)℧i,kα,β,η,σ︷=∑k=0i℧i,kα,βΓ1+kηk+σΓ1+k+σ.From (44) and (45), it follows that(46)ϕnc0Iησut=∑i=0mci℧i,kα,β,η,σ︷ctn.We may also write(47)℧i,kα,β,η,σ︷ctn=∑j=0mdi,jPτ,jα,βt,where d(i,j) are given by(48)di,j=℧i,kα,β,η,σ︷Rτ,jα,β∫0τctnwτα,βtPτ,jα,βtdt.After simplification, we have(49)di,j=c℧i,kα,β,η,σ︷Rτ,jα,β∑l=0j℧j,lα,β∫0τtl+β+nτ-tαdt,=c℧i,kα,β,η,σ︷Rη,jα,β∑l=0j℧j,lα,βΓl+β+n+1Γα+1Γl+β+n+α+1τl+β+n+α+1.Now using relation (45) we can write(50)di,j=∑k=0i℧i,kα,β∑l=0j℧j,lα,βcΓ1+kΓl+β+n+1Γα+1Rτ,jα,βΓ1+k+σΓl+β+n+α+1τl+β+n+α+1ηk+σ.Using (50) in (46), we obtain(51)ϕnc0IησYt=∑i=0m∑j=0mCidi,jPτ,jα,βt,which can be written in matrix form as(52)ϕnc0IησYt=KMTQM×Mσ,ϕnc,ηΨMt, where the entries of the matrix QM×M(σ,ϕnc,η) are as defined in (50).


## 4. Application of the Operational Matrices

In this section we apply the new operational matrices to solve coupled system of nonlinear fractional differential equations. The following lemmas are of basic importance in our further investigation.Lemma 8.
IfDσu(t)=WMTΨM(t), where 1<σ≤2, and(53)u0=u0,uτ=∑i=1m-2ζiuξi,where ζi and ξi are as defined in (2), and if ∑i=0m-2ζiξi-τ≠0, then u(t) can be written as(54)ut=WMTV1ΨMt,and the matrix V1 is defined as(55)V1=HM×Mτ,σ+QM×Mσ,ϕ1ρ,τ-∑i=1m-2QM×Mσ,ϕ1ϱi,ξi+F1MT,where ρ=τ/(∑i=0m-2ζiξi-τ), ϱi=ξi/(∑i=0m-2ζiξi-τ), and F1MTΨM(t)=u0+(t/∑i=0m-2ζiξi-τ)-∑i=1m-2ζiu0+u0.Proof.
Consider(56)Dσut=WMTΨMt.Application of fractional integration of order σ implies that(57)ut=WMTHM×Mτ,σΨMt+c0+c1t,where c0 and c1 are constants of integration. Using the initial condition u(0)=u0 implies that c0=u0. Application of nonlocal boundary condition implies that (58)∑i=1m-2ζiuξi=WMT∑i=1m-2ζiHM×Mτ,σΨMξi+∑i=1m-2ζiu0+c1∑i=1m-2ζiξi,uτ=WMTHM×Mτ,σΨMτ+u0+c1τ.Equating the above two equations we get (59)c1=1∑i=0m-2ζiξi-τWMTHM×Mτ,σΨMτ-WMT∑i=1m-2ζiHM×Mτ,σΨMξi-∑i=1m-2ζiu0+u0.Using the value of c1 in (57) and making use of Lemma 6 we get the following estimates:(60)ut=WMTHM×Mτ,σΨMt+WMTQM×Mσ,ϕ1ρ,τΨMt-WMT∑i=1m-2QM×Mσ,ϕ1ϱi,ξiΨMt+F1MTΨMt,where ρ=1/(∑i=0m-2ζiξi-τ), ϱi=ξi/(∑i=0m-2ζiξi-τ), and (61)F1MTΨMt=u0+t∑i=0m-2ζiξi-τ-∑i=1m-2ζiu0+u0.After simplification we get (62)ut=WMTHM×Mτ,σ+QM×Mσ,ϕ1ρ,τ-∑i=1m-2QM×Mσ,ϕ1ϱi,ξi+F1MTΨMt,which completes the proof.Lemma 9.
IfDσu(t)=WMTΨM(t), where 1<σ≤2 and(63)u′0=a1u0+a2uτ,u′τ=b1u0+b2uτ,then(64)ut=WMTV2ΨMt,where the matrix V2 is defined as(65)V2=HM×Mτ,σ+QM×Mσ,ϕ0d1,τ+QM×Mσ-1,ϕ0d2,τ+QM×Mσ,ϕ1d3,τ+QM×Mσ-1,ϕ1d4,τ,where d1=(λ1λ3(1-a2τ)-a2λ2)/λ2λ3, d2=-λ3(1-a2τ)/λ2, d3=λ1/λ2, and d4=-λ3/λ2. Here one assumes that λ2,λ3≠0 and are defined as λ1=(b1+b2)a2-(a1+a2)b2, λ2=(b1+b2)(1-a2τ)-(a1+a2)(1-b2τ), and λ3=(a1+a2).Proof.
Consider(66)Dσut=WMTΨMt,which implies that(67)ut=WMTHM×Mτ,σΨMt+c0+c1t.From the above equation we can also write(68)u′t=WMTHM×Mτ,σ-1ΨMt+c1.Now using u′(0)=a1u(0)+a2u(τ) and u′(τ)=b1u(0)+b2u(τ) we get(69)1-a2τc1-a1+a2c0=a2WMTHM×Mτ,σΨMτ.1-b2τc1-b1+b2c0=b2WMTHM×Mτ,σΨMτ-WMTHM×Mτ,σ-1ΨMτ,which can be easily solved for c0 and c1, and we can write (70)c1=λ1λ2WMTHM×Mτ,σΨMτ-λ3λ2WMTHM×Mτ,σ-1ΨMτ,c0=λ1λ31-a2τ-a2λ2λ2λ3WMTHM×Mτ,σΨMτ-λ31-a2τλ2WMTHM×Mτ,σ-1ΨMτ,where λ1=(b1+b2)a2-(a1+a2)b2, λ2=(b1+b2)(1-a2τ)-(a1+a2)(1-b2τ), and λ3=(a1+a2).
Now using the value ofc0 and c1 in (67) we get the following estimates:(71)ut=WMTHM×Mτ,σΨMt+d1WMTHM×Mτ,σΨMτ+d2WMTHM×Mτ,σ-1ΨMτ+d3tWMTHM×Mτ,σΨMτ+d4tWMTHM×Mτ,σ-1ΨMτ,where d1=(λ1λ3(1-a2τ)-a2λ2)/λ2λ3, d2=-λ3(1-a2τ)/λ2, d3=λ1/λ2, and d4=-λ3/λ2.
In view of Lemma6 we can write the above relation as(72)ut=WMTHM×Mτ,σΨMt+WMTQM×Mσ,ϕ0d1,τΨMt+WMTQM×Mσ-1,ϕ0d2,τΨMt+WMTQM×Mσ,ϕ1d3,τΨMt+WMTQM×Mσ-1,ϕ1d4,τΨMt,which can be written as(73)ut=WMTHM×Mτ,σ+QM×Mσ,ϕ0d1,τ+QM×Mσ-1,ϕ0d2,τ+QM×Mσ,ϕ1d3,τ+QM×Mσ-1,ϕ1d4,τΨMt,which completes proof of the lemma.


### 4.1. Couple System of FDEs with Variable Coefficients

Consider the following coupled system of FDEs:(74)uσ1t=χ11tut+χ12tu1γt+χ13tvt+χ14tv1γt+ft,vσ1t=χ21tut+χ22tu1γt+χ23tvt+χ24tv1γt+gt,subject to one of the above discussed boundary conditions ((2) or (4)). Assume that u(σ1)(t)=XMTΨM(t) and v(σ1)(t)=YMTΨM(t). Then in view of Lemma 8 or Lemma 9 (depends on the given set of boundary conditions) we can write(75)ut=XMTViuΨMt,vt=YMTVivΨMt.Note that the subscript i is used to distinguish between different types of boundary condition. If we are given nonlocal boundary conditions we will use i=1 and if we are given mix derivative boundary condition we will use i=2. Superscripts u and v are used to distinguish between conditions for u(t) and v(t), respectively. Now in view of Lemma 6, we can write (74) as(76)XMTΨMt=XMTViuKχ110ΨMt+XMTViuKχ12γ1ΨMt+YMTVivKχ130ΨMt+YMTVivKχ14γ1ΨMt+FMTΨMt,YMTΨMt=XMTViuKχ210ΨMt+XMTViuKχ22γ1ΨMt+YMTVivKχ230ΨMt+YMTVivKχ24γ1ΨMt+GMTΨMt.which can be converted into easily solvable matrix equation (by following the same step as in [22] from (51) to (69)). The generalized form of the resulting matrix equation is given as (77)AX+XB=C,where A, B, and C are of compatible size. The vector X is the unknown solution to be determined. The solution method of such type of equations is explained in [26] and a detailed algorithm is presented. MatLab command dlyap uses the same algorithm presented in [26] for solution to such type of matrix equations.


### 4.2. Convergence Analysis

ConsiderΩi,ju and Ωi,jv as entries of Viu, and Viv, respectively. Also assume that Λi,jχ11,γ1 are entries of Kχ24γ1. Now we can write(78)uσ1t=∑i=0∞uiPτ,iα,βt=XMTΨMt+∑i=m+1∞uiPτ,iα,βt,vσ2t=∑i=0∞viPτ,iα,βt=YMTΨMt+∑i=m+1∞viPτ,iα,βt.Using the above estimates and using the same procedure we can write(79)XMTΨMt-XMTViuKχ110ΨMt-XMTViuKχ12γ1ΨMt-YMTVivKχ130ΨMt-YMTVivKχ14γ1ΨMt-FMTΨMt=Rut,YMTΨMt-XMTViuKχ210ΨMt-XMTViuKχ22γ1ΨMt-YMTVivKχ230ΨMt-YMTVivKχ24γ1ΨMt-GMTΨMt=Rvt,where(80)Rut=∑i=m+1∞ui∑h=0mΩi,hu∑j=0mΛh,jχ11,0Pτ,jα,βt+∑i=m+1∞ui∑h=0mΩi,hu∑j=0mΛh,jχ12,γ1Pτ,jα,βt∑i=m+1∞vi∑h=0mΩi,hv∑j=0mΛh,jχ13,0Pτ,jα,βt+∑i=m+1∞vi∑h=0mΩi,hv∑j=0mΛh,jχ14,γ1Pτ,jα,βt+∑i=m+1∞fiPτ,iα,βt-∑i=m+1∞uiPτ,iα,βt,Rvt=∑i=m+1∞ui∑h=0mΩi,hu∑j=0mΛh,jχ21,0Pτ,jα,βt+∑i=m+1∞ui∑h=0mΩi,hu∑j=0mΛh,jχ22,γ1Pτ,jα,βt∑i=m+1∞vi∑h=0mΩi,hv∑j=0mΛh,jχ23,0Pτ,jα,βt+∑i=m+1∞vi∑h=0mΩi,hv∑j=0mΛh,jχ24,γ1Pτ,jα,βt+∑i=m+1∞giPτ,iα,βt-∑i=m+1∞uiPτ,iα,βt,where fi and gi are the approximation coefficients of f(t) and g(t), respectively. In view of maximum norm and (12) we can write(81)Rut≤∑i=m+1∞ui∑h=0m∑j=0mΩi,huΛh,jχ11,0C^j,ϵ+∑i=m+1∞ui∑h=0m∑j=0mΩi,huΛh,jχ12,γ1C^j,ϵ∑i=m+1∞vi∑h=0m∑j=0mΩi,hvΛh,jχ13,0C^j,ϵ+∑i=m+1∞vi∑h=0m∑j=0mΩi,hvΛh,jχ14,γ1C^j,ϵ+∑i=m+1∞fiC^i,ϵ-∑i=m+1∞uiC^i,ϵ,Rvt=∑i=m+1∞ui∑h=0m∑j=0mΩi,huΛh,jχ21,0C^j,ϵ+∑i=m+1∞ui∑h=0m∑j=0mΩi,huΛh,jχ22,γ1C^j,ϵ∑i=m+1∞vi∑h=0m∑j=0mΩi,hvΛh,jχ23,0C^j,ϵ+∑i=m+1∞vi∑h=0m∑j=0mΩi,hvΛh,jχ24,γ1C^j,ϵ+∑i=m+1∞giC^i,ϵ-∑i=m+1∞uiC^i,ϵ.In view of Lemma 2.1 of [22] the expansion coefficients ui and vi decay faster than any of the algebraic order of index i. We can write ui, vi, fi, and gi approach to zero as i→∞ and as a result we conclude that(82)Rut⟶0asm⟶∞,Rvt⟶0asm⟶∞. We see that the convergence of the scheme depends on the decay of expansion coefficients, and the decay of expansion coefficients depends solely on the smoothness of solution for the problem. Therefore if the solution for the problem is smooth we get high accuracy at small scale level, and the accuracy will increase if we increase the scale level.


### 4.3. Nonlinear FDEs withm-Point Boundary Conditions

To solve nonlinear FDEs we will implement operational matrix method combined with qasilinearization method. Consider the following nonlinear FDE:(83)uσ1=ft,uγ1,vγ1,u,v,vσ2=gt,uγ1,vγ1,u,v.In qasilinearization method the nonlinear differential equation is converted into linear differential equation with variable coefficients. The nonlinear part of the differential equation is linearized about some function. In most cases the initial function is selected according to the physical problem under consideration, although in some cases it is selected as the solution for the linear part [27]. We will first solve the linear part of (83):(84)uσ1=L1u,v,u1γ,v1γ,vσ2=L2u,v,u1γ,v1γ,under given boundary conditions using the method developed in the previous section. L1 and L2 are the linear part of f and g, respectively. The solution for this problem will be labeled as u0(t) and v0(t). The next step is to linearize the nonlinear part of f and g about u0(t) and v0(t) using Taylor series expansion. As a result the nonlinear differential equation will be converted into linear differential equation with variable coefficients, which can be solved by the proposed. The whole process can be seen as a recurrence relation as(85)ur+1σ1=fr+∂fr∂uur+1-ur+∂fr∂vvr+1-vrt∂fr∂uγ1ur+1γ1-urγ1+∂fr∂vγ1vr+1γ1-vrγ1,vr+1σ1=gr+∂gr∂uur+1-ur+∂gr∂vvr+1-vrt∂gr∂uγ1ur+1γ1-urγ1+∂gr∂vγ1vr+1γ1-vrγ1,where fr=f(t,urγ1,vrγ1,ur,vr) and gr=g(t,urγ1,vrγ1,ur,vr). The above equation is linear fractional differential equations with variable coefficients and can be easily solved with the method developed in the previous section.


## 4.1. Couple System of FDEs with Variable Coefficients

Consider the following coupled system of FDEs:(74)uσ1t=χ11tut+χ12tu1γt+χ13tvt+χ14tv1γt+ft,vσ1t=χ21tut+χ22tu1γt+χ23tvt+χ24tv1γt+gt,subject to one of the above discussed boundary conditions ((2) or (4)). Assume that u(σ1)(t)=XMTΨM(t) and v(σ1)(t)=YMTΨM(t). Then in view of Lemma 8 or Lemma 9 (depends on the given set of boundary conditions) we can write(75)ut=XMTViuΨMt,vt=YMTVivΨMt.Note that the subscript i is used to distinguish between different types of boundary condition. If we are given nonlocal boundary conditions we will use i=1 and if we are given mix derivative boundary condition we will use i=2. Superscripts u and v are used to distinguish between conditions for u(t) and v(t), respectively. Now in view of Lemma 6, we can write (74) as(76)XMTΨMt=XMTViuKχ110ΨMt+XMTViuKχ12γ1ΨMt+YMTVivKχ130ΨMt+YMTVivKχ14γ1ΨMt+FMTΨMt,YMTΨMt=XMTViuKχ210ΨMt+XMTViuKχ22γ1ΨMt+YMTVivKχ230ΨMt+YMTVivKχ24γ1ΨMt+GMTΨMt.which can be converted into easily solvable matrix equation (by following the same step as in [22] from (51) to (69)). The generalized form of the resulting matrix equation is given as (77)AX+XB=C,where A, B, and C are of compatible size. The vector X is the unknown solution to be determined. The solution method of such type of equations is explained in [26] and a detailed algorithm is presented. MatLab command dlyap uses the same algorithm presented in [26] for solution to such type of matrix equations.


## 4.2. Convergence Analysis

ConsiderΩi,ju and Ωi,jv as entries of Viu, and Viv, respectively. Also assume that Λi,jχ11,γ1 are entries of Kχ24γ1. Now we can write(78)uσ1t=∑i=0∞uiPτ,iα,βt=XMTΨMt+∑i=m+1∞uiPτ,iα,βt,vσ2t=∑i=0∞viPτ,iα,βt=YMTΨMt+∑i=m+1∞viPτ,iα,βt.Using the above estimates and using the same procedure we can write(79)XMTΨMt-XMTViuKχ110ΨMt-XMTViuKχ12γ1ΨMt-YMTVivKχ130ΨMt-YMTVivKχ14γ1ΨMt-FMTΨMt=Rut,YMTΨMt-XMTViuKχ210ΨMt-XMTViuKχ22γ1ΨMt-YMTVivKχ230ΨMt-YMTVivKχ24γ1ΨMt-GMTΨMt=Rvt,where(80)Rut=∑i=m+1∞ui∑h=0mΩi,hu∑j=0mΛh,jχ11,0Pτ,jα,βt+∑i=m+1∞ui∑h=0mΩi,hu∑j=0mΛh,jχ12,γ1Pτ,jα,βt∑i=m+1∞vi∑h=0mΩi,hv∑j=0mΛh,jχ13,0Pτ,jα,βt+∑i=m+1∞vi∑h=0mΩi,hv∑j=0mΛh,jχ14,γ1Pτ,jα,βt+∑i=m+1∞fiPτ,iα,βt-∑i=m+1∞uiPτ,iα,βt,Rvt=∑i=m+1∞ui∑h=0mΩi,hu∑j=0mΛh,jχ21,0Pτ,jα,βt+∑i=m+1∞ui∑h=0mΩi,hu∑j=0mΛh,jχ22,γ1Pτ,jα,βt∑i=m+1∞vi∑h=0mΩi,hv∑j=0mΛh,jχ23,0Pτ,jα,βt+∑i=m+1∞vi∑h=0mΩi,hv∑j=0mΛh,jχ24,γ1Pτ,jα,βt+∑i=m+1∞giPτ,iα,βt-∑i=m+1∞uiPτ,iα,βt,where fi and gi are the approximation coefficients of f(t) and g(t), respectively. In view of maximum norm and (12) we can write(81)Rut≤∑i=m+1∞ui∑h=0m∑j=0mΩi,huΛh,jχ11,0C^j,ϵ+∑i=m+1∞ui∑h=0m∑j=0mΩi,huΛh,jχ12,γ1C^j,ϵ∑i=m+1∞vi∑h=0m∑j=0mΩi,hvΛh,jχ13,0C^j,ϵ+∑i=m+1∞vi∑h=0m∑j=0mΩi,hvΛh,jχ14,γ1C^j,ϵ+∑i=m+1∞fiC^i,ϵ-∑i=m+1∞uiC^i,ϵ,Rvt=∑i=m+1∞ui∑h=0m∑j=0mΩi,huΛh,jχ21,0C^j,ϵ+∑i=m+1∞ui∑h=0m∑j=0mΩi,huΛh,jχ22,γ1C^j,ϵ∑i=m+1∞vi∑h=0m∑j=0mΩi,hvΛh,jχ23,0C^j,ϵ+∑i=m+1∞vi∑h=0m∑j=0mΩi,hvΛh,jχ24,γ1C^j,ϵ+∑i=m+1∞giC^i,ϵ-∑i=m+1∞uiC^i,ϵ.In view of Lemma 2.1 of [22] the expansion coefficients ui and vi decay faster than any of the algebraic order of index i. We can write ui, vi, fi, and gi approach to zero as i→∞ and as a result we conclude that(82)Rut⟶0asm⟶∞,Rvt⟶0asm⟶∞. We see that the convergence of the scheme depends on the decay of expansion coefficients, and the decay of expansion coefficients depends solely on the smoothness of solution for the problem. Therefore if the solution for the problem is smooth we get high accuracy at small scale level, and the accuracy will increase if we increase the scale level.


## 4.3. Nonlinear FDEs withm-Point Boundary Conditions

To solve nonlinear FDEs we will implement operational matrix method combined with qasilinearization method. Consider the following nonlinear FDE:(83)uσ1=ft,uγ1,vγ1,u,v,vσ2=gt,uγ1,vγ1,u,v.In qasilinearization method the nonlinear differential equation is converted into linear differential equation with variable coefficients. The nonlinear part of the differential equation is linearized about some function. In most cases the initial function is selected according to the physical problem under consideration, although in some cases it is selected as the solution for the linear part [27]. We will first solve the linear part of (83):(84)uσ1=L1u,v,u1γ,v1γ,vσ2=L2u,v,u1γ,v1γ,under given boundary conditions using the method developed in the previous section. L1 and L2 are the linear part of f and g, respectively. The solution for this problem will be labeled as u0(t) and v0(t). The next step is to linearize the nonlinear part of f and g about u0(t) and v0(t) using Taylor series expansion. As a result the nonlinear differential equation will be converted into linear differential equation with variable coefficients, which can be solved by the proposed. The whole process can be seen as a recurrence relation as(85)ur+1σ1=fr+∂fr∂uur+1-ur+∂fr∂vvr+1-vrt∂fr∂uγ1ur+1γ1-urγ1+∂fr∂vγ1vr+1γ1-vrγ1,vr+1σ1=gr+∂gr∂uur+1-ur+∂gr∂vvr+1-vrt∂gr∂uγ1ur+1γ1-urγ1+∂gr∂vγ1vr+1γ1-vrγ1,where fr=f(t,urγ1,vrγ1,ur,vr) and gr=g(t,urγ1,vrγ1,ur,vr). The above equation is linear fractional differential equations with variable coefficients and can be easily solved with the method developed in the previous section.


## 5. Illustrative Examples

We show the applicability of the method by solving some test problems. Where available we compare our results with results obtained with other methods.Example 1.
Consider the following linear fractional order three-point boundary value problem [28]: (86)D3/2Xt+e-3ππXt=gt,t∈0,1,X0=0,X1=-125196X25,where the forcing term g(t) is defined as (87)gt=e-3ππt240t2-74t+33+4e3πt128t2-148t+33.The exact solution for this problem is X(t)=t5-37t3/20+33t2/40. This problem is solved in [28] using Haar wavelets (HW). In [29] this problem is solved using improved reproducing kernel method (RKM). We compare the absolute error obtained with the proposed method with the error obtained with HW and RKM. The results are shown in Table 1.Table 1
Comparison of absolute error of Example1 obtained with the proposed method and its comparison with Haar wavelets [28] and RKM [29].
M = 8 (HW) [28] n = 8 (RKM) [29] M = 8 (PM) M = 32 (HW) [28] M = 25 (PM) t = 0.1 6.97 × 10-5 2.41 × 10-5 4.192 × 10-6 4.42 × 10-7 6.05 × 10-9 t = 0.2 4.59 × 10-5 2.10 × 10-5 2.13 × 10-5 1.78 × 10-7 3.15 × 10-9 t = 0.3 4.59 × 10-5 3.42 × 10-6 1.47 × 10-5 1.80 × 10-7 1.04 × 10-9 t = 0.4 1.33 × 10-4 2.23 × 10-5 4.45 × 10-6 4.42 × 10-7 4.04 × 10-9 t = 0.5 1.37 × 10-4 2.26 × 10-5 2.28 × 10-5 5.35 × 10-7 6.10 × 10-9 t = 0.6 1.13 × 10-4 2.93 × 10-5 9.55 × 10-6 4.41 × 10-7 7.98 × 10-9 t = 0.7 5.80 × 10-5 3.27 × 10-5 1.09 × 10-5 1.78 × 10-7 1.80 × 10-9 t = 0.8 6.68 × 10-5 3.00 × 10-5 1.36 × 10-5 1.78 × 10-7 8.82 × 10-9 t = 0.9 6.27 × 10-4 6.15 × 10-5 1.50 × 10-5 4.40 × 10-7 1.34 × 10-8Example 2.
Consider the following fractional differential equations having variable coefficients and five-point boundary conditions [29, 30]:(88)D1.3Xt+cos⁡tDXt+2Xt=ft,t∈0,1,X0=0,X18+2X12+3149X78=X1,where ft=2t2+2tcos⁡(t)+Γ(3)/Γ1.7t0.7. The exact solution for the problem is X(t)=t2. We approximate solution for this problem with the proposed method and compare its absolute error with the error reported in [29, 30] (in these references reproducing kernel method and improved reproducing kernel method are used to solve this problem). We observe that the solution obtained with the proposed method is in good agreement with the exact solution. In Table 2 absolute error at different value of t is compared with the error reported in [29, 30].Table 2
Comparison of absolute error of Example2 obtained with the proposed method and its comparison with RKM [30] and improved reproducing kernel method [29].
t = 0.1 1.11 × 10-7 3.09 × 10-8 7.35 × 10-8 t = 0.2 2.87 × 10-7 1.69 × 10-8 2.43 × 10-9 t = 0.3 5.51 × 10-7 3.55 × 10-9 9.54 × 10-9 t = 0.4 9.23 × 10-7 1.89 × 10-9 2.17 × 10-10 t = 0.5 1.41 × 10-6 2.94 × 10-10 9.37 × 10-11 t = 0.6 2.03 × 10-6 8.47 × 10-9 7.58 × 10-11 t = 0.7 2.77 × 10-6 2.97 × 10-8 4.72 × 10-11 t = 0.8 3.64 × 10-6 3.94 × 10-8 8.24 × 10-12 t = 0.9 4.62 × 10-6 5.23 × 10-8 3.57 × 10-10Example 3.
Consider the following nonlinear fractional differential equations with 6-point nonlocal boundary value problem:(89)D1.8Xt+X′t+Xt=Xt+X′t2-XtX′t+ft,t∈0,1,X0=1,12X12+23X58+14X34-0.646X78=X1,where the forcing term f(t) is defined as (90)ft=20t4-4t2-2t+11t5+4t6-25t8-5t9-t10+2229536516744740625t16/5144115188075855872-4.The exact and unique solution for the problem is X(t)=t5-2t+1. We approximate the solution for this problem with the proposed iterative method. This example is solved using the parameters α=0, β=2, and M=8. We observe that the method provides good approximation to solution for the problem. The approximate solution approaches the exact solution as we make iteration. Absolute difference of exact and approximate solution at different stages of iteration is shown in Table 3. It can be easily seen that at fourth iteration the absolute difference is less than 10-8.Table 3
Absolute difference of exact and approximate solution of Example3 at different stages of the iteration.
1st iteration 2nd iteration 3rd iteration 4th iteration t = 0.1 4.42 × 10-2 1.87 × 10-3 9.17 × 10-6 1.65 × 10-8 t = 0.2 4.18 × 10-2 3.04 × 10-3 1.12 × 10-6 2.03 × 10-8 t = 0.3 2.70 × 10-2 3.56 × 10-3 2.30 × 10-6 1.80 × 10-8 t = 0.4 1.42 × 10-2 3.93 × 10-3 2.92 × 10-6 2.77 × 10-8 t = 0.5 6.92 × 10-3 4.09 × 10-3 3.25 × 10-6 2.70 × 10-8 t = 0.6 3.83 × 10-3 4.04 × 10-3 3.33 × 10-6 2.04 × 10-8 t = 0.7 2.81 × 10-3 3.96 × 10-3 3.26 × 10-6 2.52 × 10-8 t = 0.8 2.76 × 10-3 3.80 × 10-3 3.15 × 10-6 2.85 × 10-8 t = 0.9 3.33 × 10-3 3.57 × 10-3 2.95 × 10-6 1.94 × 10-8 t = 1.0 4.58 × 10-3 3.30 × 10-3 2.71 × 10-6 2.48 × 10-8Example 4.
Consider the following nonlinear differential equation with 6-point nonlocal boundary conditions:(91)D2Xt+X′t+Xt=X′t2+XtX′t+ft,t∈0,1,X0=0,0.8X12+0.75X58+0.50X34-1.11X78=X1.The forcing term f(t) is defined as (92)ft=-sint/329-2sin⁡t/3sin⁡t/623-16sin⁡t/629+89.
The exact solution for the problem isX(t)=cos⁡(t/3). We carry out the same analysis as in the previous example. The same conclusion is made; the solution converges rapidly to the exact solution for the problem. In Table 4 absolute error at different stages of iteration is given. For this example the absolute error is less than 10-16 at fourth iteration. For this example we set α=1, β=1, and M=10.Table 4
Absolute difference of exact and approximate solution of Example4 at different stages of the iteration.
t = 0.1 1.61 × 10-5 5.27 × 10-8 3.10 × 10-15 1.24 × 10-17 t = 0.2 2.80 × 10-5 8.29 × 10-8 4.88 × 10-15 2.82 × 10-17 t = 0.3 3.62 × 10-5 9.93 × 10-8 5.63 × 10-15 2.23 × 10-17 t = 0.4 4.13 × 10-5 1.07 × 10-7 2.16 × 10-15 1.26 × 10-16 t = 0.5 4.41 × 10-5 1.11 × 10-7 5.09 × 10-15 1.12 × 10-17 t = 0.6 4.52 × 10-5 1.13 × 10-7 6.21 × 10-15 2.22 × 10-17 t = 0.7 4.53 × 10-5 1.13 × 10-7 6.10 × 10-15 2.22 × 10-17 t = 0.8 4.48 × 10-5 1.12 × 10-7 5.88 × 10-15 1.01 × 10-17 t = 0.9 4.40 × 10-5 1.10 × 10-7 5.77 × 10-15 4.42 × 10-17 t = 1.0 4.28 × 10-5 1.08 × 10-7 5.44 × 10-15 2.21 × 10-17Example 5.
Consider the Lane-Emden equations [31]: (93)uσt+ϕ1tu′t-c11u2t-c12utvt=0,vσt+ϕ2tv′t-c21u2t-c22utvt=0.Systems of Lane-Emden equations arise in the modeling of several physical phenomena, such as pattern formation, population evolution, and chemical reactions. The parameters c11, c12, c21, and c22 can be considered from the actual chemical reaction or dynamics under consideration. We solve this problem with the proposed method subject to the following type of boundary conditions:(94)u′0=2u0-u1,u′1=u0+2u1,v′0=-2v0+v1,v′1=v0-2vτ.In [31], Adomian decomposition method is employed to get approximate solution to this problem. We compare our results with the results reported in [31]. We observe that the current method provides a very good approximation to the solutions for the problem. We calculate the error remainder term: (95)ERnu=unσt+ϕ1tun′t-c11un2t-c12untvnt,ERnv=vnσt+ϕ2tvn′t-c21un2t-c22untvnt,where un and vn are the approximate solution for the problem at different stages of iteration n. We calculate the maximal error remainder term:(96)MERnu=maxt∈0,1⁡ERnu,MERnv=maxt∈0,1⁡ERnv.The comparison of results of the proposed method with Adomain decomposition is presented in Table 5.Table 5
Maximal error of Example5 obtained with the proposed method and its comparison with error reported in [31].
M E R n u [31] M E R n v [31] M E R n u (PM) M E R n v (PM) n = 2 0.86666 1.14667 1.12626 3.92536 n = 3 0.44282 0.57063 0.59223 2.63829 n = 4 0.23094 0.29382 0.01280 0.89223 n = 5 0.12671 0.15988 1.12 × 10-3 7.29 × 10-2 n = 6 0.07155 0.08982 9.36 × 10-3 4.96 × 10-2 n = 7 0.04160 0.05202 3.10 × 10-3 9.41 × 10-3 n = 8 0.02467 0.03077 8.72 × 10-3 8.67 × 10-3 n = 9 0.01490 0.01855 6.94 × 10-3 6.62 × 10-4 n = 10 0.00913 0.01134 1.76 × 10-3 1.93 × 10-4 n = 11 0.00566 0.00702 6.82 × 10-4 3.97 × 10-5 n = 12 0.00354 0.00439 7.32 × 10-5 1.35 × 10-5 n = 13 0.00224 0.00277 2.36 × 10-5 6.78 × 10-5 n = 14 0.00142 0.00176 7.91 × 10-6 0.41 × 10-6 n = 15 0.00091 0.00112 9.45 × 10-7 8.93 × 10-6 n = 16 0.00058 0.00072 7.34 × 10-7 8.54 × 10-6


## 6. Conclusion and Future Work

Form the above analysis and observation we conclude that the method works very well and efficiently solve fractional order linear and nonlinear differential equation and coupled system under multipoint nonlocal boundary conditions and mixed derivative boundary conditions. When the absolute error is compared with reproducing kernel method and Adomain decomposition method we observe that the method yields very accurate solution. In this work shifted Jacobi polynomials are used; it may also be possible that the reader may get more accurate solution by using another class of orthogonal polynomials like Bernstein polynomials or Hermite polynomials and so forth. Our future work is related to investigating the best choice of orthogonal polynomials.



---

*Source: 5601821-2016-09-06.xml*