# An Adaptive Pseudospectral Method for Fractional Order Boundary Value Problems

Authors: Mohammad Maleki; Ishak Hashim; Majid Tavassoli Kajani; Saeid Abbasbandy
Journal: Abstract and Applied Analysis
 (2012)

---

## Abstract

An adaptive pseudospectral method is presented for solving a class of multiterm fractional boundary value problems (FBVP) which involve Caputo-type fractional derivatives. The multiterm FBVP is first converted into a singular Volterra integrodifferential equation (SVIDE). By dividing the interval of the problem to subintervals, the unknown function is approximated using a piecewise interpolation polynomial with unknown coefficients which is based on shifted Legendre-Gauss (ShLG) collocation points. Then the problem is reduced to a system of algebraic equations, thus greatly simplifying the problem. Further, some additional conditions are considered to maintain the continuity of the approximate solution and its derivatives at the interface of subintervals. In order to convert the singular integrals of SVIDE into nonsingular ones, integration by parts is utilized. In the method developed in this paper, the accuracy can be improved either by increasing the number of subintervals or by increasing the degree of the polynomial on each subinterval. Using several examples including Bagley-Torvik equation the proposed method is shown to be efficient and accurate.


---

## Body

## 1. Introduction

Due to the development of the theory of fractional calculus and its applications, such as in the fields of physics, Bode’s analysis of feedback amplifiers, aerodynamics and polymer rheology, and so forth, many works on the basic theory of fractional calculus and fractional order differential equations have been established [1–3].In general, the analytical solutions for most of the fractional differential equations are not readily attainable, and thus the need for finding efficient computational algorithms for obtaining numerical solutions arises. Recently, there have been many papers dealing with the solutions of initial value and boundary value problems for linear and nonlinear fractional differential equations. These methods include finite difference approximation method [4], collocation method [5, 6], the Adomian decomposition method [7, 8], variational iteration method [9–12], operational matrix methods [13–16], and homotopy methods [17, 18]. In [19] suitable spline functions of polynomial form are derived and used to solve linear and nonlinear fractional differential equations. The authors of [20] have investigated the existence and multiplicity of positive solutions of a nonlinear fractional differential equation initial value problem. Furthermore, some physical and geometrical interpretations of fractional operators and fractional differential equations have been of concern to many authors [12, 21, 22].In the present paper, we intend to introduce an efficient adaptive pseudospectral method for multiterm fractional boundary value problems (FBVP) of the form(1.1)F(x,y(x),Dα1y(x),…,Dαmy(x))=0,x∈[0,L],
subject to
(1.2)Hr(y(ξ0),…,y(l)(ξ0),…,y(ξl),…,y(l)(ξl))=0,r=0,1,…,l,
where F can be nonlinear in general, 0<α1<α2<⋯<αm, l<αm≤l+1, L∈ℝ, Hr are linear functions, the points ξ0,ξ1, …,ξl lie in [0,L], and Dαq denotes the Caputo-fractional derivative of order αq, defined as follows [23]:
(1.3)Dαqy(x)=1Γ(nq-αq)∫0x‍y(nq)(t)(x-t)αq+1-nqdt,nq=[αq]+1,q=1,2,…,m,
where [αq] denotes the integer part of the real number αq. For details about the mathematical properties of fractional derivatives, see [2].In this method, the multi-term FBVP is first converted into a singular Volterra integrodifferential equation (SVIDE). By dividing the interval of the problem to subintervals, the unknown function is approximated using a piecewise interpolation polynomial with unknown coefficients which is based on shifted Legendre-Gauss (ShLG) collocation points. Then the problem is reduced to a system of algebraic equations using collocation. Further, some additional conditions are considered to maintain the continuity of the approximate solution and its firstl derivatives at the interface of subintervals. The singular integrals of SVIDE are converted into nonsingular ones by utilizing integration by parts and thus greatly improve the accuracy and convergence rate of the approximate solution. The main characteristics of the method are that it converts the FBVP into a system of algebraic equations which greatly simplifies it. In addition, in the method developed in this paper, the accuracy can be improved either by increasing the number of subintervals or by increasing the degree of the polynomial on each subinterval. The present adaptive pseudospectral method can be implemented for FBVPs defined in large domains. Moreover, this new algorithm also works well even for some solutions having oscillatory behavior. Numerical examples including Bagley-Torvik equation subject to boundary conditions are also presented to illustrate the accuracy of the present scheme. Finally, in order to have a physical understanding of fractional differential equations, the derivation of Bagley-Torvik equation is given in the appendix.The outline of this paper is as follows. In Section2, some basic properties of Legendre and shifted Legendre polynomials, which are required for our subsequent development, are first presented. Piecewise polynomials interpolation based on ShLG points and its convergence properties are then investigated, and finally the adaptive pseudospectral method for FBVPs is explained. Section 3 is devoted to some numerical examples. In Section 4, a brief conclusion is given. The appendix is given which consists of the derivation of Bagley-Torvik equation.


## 2. The Adaptive Pseudospectral Method for FBVPs

In this section we drive the adaptive pseudospectral method based on ShLG collocation points and apply it to solve the nonlinear multi-term FBVP (1.1)-(1.2).


### 2.1. Review of Legendre and Shifted Legendre Polynomials

The Legendre polynomials,Pi(z), i=0,1,2,…, are the eigenfunctions of the singular Sturm-Liouville problem
(2.1)[(1-z2)Pi′(z)]′+i(i+1)Pi(z)=0.
Also, they are orthogonal with respect to L2 inner product on the interval [-1,1] with the weight function w(z)=1, that is,
(2.2)∫-11‍Pi(z)Pj(z)dz=22i+1δij,
where δij is the Kronecker delta. The Legendre polynomials satisfy the recursion relation
(2.3)Pi+1(z)=2i+1i+1zPi(z)-ii+1Pi-1(z),
where P0(z)=1 and P1(z)=z. If Pi(z) is normalized so that Pi(1)=1, then, for any i, the Legendre polynomials in terms of power of z are
(2.4)Pi(z)=12i∑k=0[i/2]‍(-1)k(ik)(2i-2ki)zi-2k,
where [i/2] denotes the integer part of i/2.The Legendre-Gauss (LG) collocation points-1<z1<z2<⋯<zN-1<1 are the roots of PN-1(z). No explicit formulas are known for the LG points; however, they are computed numerically using existing subroutines. The LG points have the property that
(2.5)∫-11‍p(z)dz=∑j=1N-1‍wjp(zj),
is exact for polynomials of degree at most 2N-3, where
(2.6)wj=2(1-zj2)[PN-1′(zj)]2,j=1,2,…,N-1,
are LG quadrature weights. For more details about Legendre polynomials, see [24].The shifted Legendre polynomials on the intervalx∈[a,b] are defined by
(2.7)P^i(x)=Pi(1b-a(2x-a-b)),i=0,1,2,…,
which are obtained by an affine transformation from the Legendre polynomials. The set of shifted Legendre polynomials is a complete L2[a,b]-orthogonal system with the weight function w(x)=1. Thus, any function f∈L2[a,b] can be expanded in terms of shifted Legendre polynomials.The ShLG collocation pointsa<x1<x2<⋯<xN-1<b on the interval [a,b] are obtained by shifting the LG points, zj, using the transformation
(2.8)xj=12((b-a)zj+a+b),j=1,2,…,N-1.
Thanks to the property of the standard LG quadrature, it follows that for any polynomial p of degree at most 2N-3 on (a,b),
(2.9)∫ab‍p(x)dx=b-a2∫-11‍p(12[(b-a)z+a+b])dz=b-a2∑j=1N-1‍wjp(12[(b-a)zj+a+b])=∑j=1N-1‍w^jp(xj),
where w^j=((b-a)/2)wj, 1⩽j⩽N-1 are ShLG quadrature weights.


### 2.2. Function Approximation

Suppose that the interval[0,L] is divided into K subintervals Ik=[(k-1)h,kh], k=1,2,…,K, where h=L/K. Let yk(x) be the solution of the problem in (1.1)-(1.2) in the subinterval Ik. Consider now the ShLG collocation points (k-1)h<xk1<⋯<xk,N-1<kh on the kth subinterval Ik, k=1,2,…,K, obtained using (2.8). Obviously,
(2.10)xkj=h2(zj+2k-1),j=1,2,…,N-1.
Also, consider two additional noncollocated points xk0=(k-1)h and xkN=kh. Let us define
(2.11)ℜkN=Span{Lk0(x),Lk1(x),…,LkN(x)},x∈Ik,
where
(2.12)Lki(x)=∏l=0,l≠iN‍x-xklxki-xkl,i=0,1,…,N,
is a basis of Lagrange interpolating polynomials on the subinterval Ik that satisfy Lki(xkj)=δij, where δij is the Kronecker delta function. The L2(Ik)-orthogonal projection IN:L2(Ik)→ℜkN is a mapping in a way that for any yk∈L2(Ik)(2.13)〈IN(yk)-yk,ϕk〉=0,∀ϕk∈ℜkN,
or equivalently
(2.14)IN(yk)(x)=∑i=0N‍ykiLki(x),x∈Ik,
where yki=yk(xki).Here, it can be easily seen that fori=0,1,…,N and k=1,2,…,K, we have
(2.15)Lki(x)=L1i(x-xk0),x∈Ik.
Thus, by utilizing (2.15) for (2.14), the approximation of yk(x) within each subinterval Ik can be restated as
(2.16)yk(x)≈IN(yk)(x)=∑i=0N‍ykiL1i(x-xk0)=YkT·Lk(x),x∈Ik,
where Yk and Lk(x) are (N+1)×1 matrices given by Yk=[yk0,…,ykN]T and Lk(x)=[L10(x-xk0),…,L1N(x-xk0)]T. It is important to observe that the series (2.16) includes the Lagrange polynomials associated with the noncollocated points xk0=(k-1)h and xkN=kh. Moreover, it is seen from (2.15)-(2.16) that, in the present adaptive scheme, it is only needed to produce the basis of Lagrange polynomials L1i(x) at the first subinterval.Bynm+1 times (nm is defined in (1.3)) differentiating of (2.16), we obtain
(2.17)yk(r)(x)≈YkT·Lk(r)(x),x∈Ik,
where Lk(r)(x)=(dr/dxr)Lk(x).


### 2.3. Convergence Rate

ForN⩾1 we introduce the piecewise polynomials space
(2.18)ΨIkN={y∈C0([0,L]):yk=y∣Ik∈ℙN(Ik)},
which is the space of the continuous functions over [0,L] whose restrictions on each subinterval Ik are polynomials of degree ⩽N. Then, for any continuous function y in [0,L], the piecewise interpolation polynomial ψN(y) coincides on each subinterval Ik with the interpolating polynomial IN(y) of yk=y|Ik at the ShLG points.In [25], with the aid of the formulas (5.4.33), (5.4.34) of [24], we prove the convergence properties of piecewise interpolation polynomial based on shifted Legendre-Gauss-Radau points in the norms of the Sobolev spaces. Accordingly, the following results for the convergence based on ShLG points hold.Theorem 2.1.
Suppose thaty∈Hv(0,L) with v⩾1. Then
(2.19)∥y-ψN(y)∥L2(0,L)⩽cN-v|y|H0;v;N;h(0,L),
and, for 1⩽u⩽v, if h⩽1, then
(2.20)∥y-ψN(y)∥Hu(0,L)⩽cN2u-(1/2)-v|y|Hu;v;N;h(0,L),
and if h>1, then
(2.21)∥y-ψN(y)∥Hu(0,L)⩽cN2u-(1/2)-v|y|H0;v;N;h(0,L).Note thatc denotes a positive constant that depends on v, but which is independent of the function y and integer N. Moreover, we introduce the seminorm of Hv(0,L), 0⩽u⩽v, N⩾0, h>0, as
(2.22)|y|Hu;v;N;h(0,L)=(∑l=min{v,N+1}v‍h2l-2u∥y(l)∥L2(0,L)2)1/2.Remark 2.2.
WheneverN⩾v-1, using (2.19)–(2.22), we get
(2.23)∥y-ψN(y)∥L2(0,L)⩽cN-vhv∥y(v)∥L2(0,L),
and, for u⩾1, if h⩽1, then
(2.24)∥y-ψN(y)∥Hu(0,L)⩽cN2u-(1/2)-vhv-u∥y(v)∥L2(0,L),
and if h>1, then
(2.25)∥y-ψN(y)∥Hu(0,L)⩽cN2u-(1/2)-vhv∥y(v)∥L2(0,L).
Equations (2.23)–(2.25) show that if y is infinitely smooth on [0,L] and h⩽1, the convergence rate of ψN(y) to y is faster than h to the power of N+1-u and any power of 1/N, which is superior to that for the global collocation method over [0,L]. Thus, the bigger the subinterval length the slower the convergence rate.


### 2.4. Problem Replacement and the Solution Technique

Consider the multi-term FBVP in (1.1)-(1.2). With substituting the definition of the Caputo-derivative (1.3) into (1.1), we can convert (1.1) into an equivalent SVIDE as
(2.26)F(x,y(x),1Γ(n1-α1)∫0x‍y(n1)(t)(x-t)α1+1-n1dt,…,1Γ(nm-αm)∫0x‍y(nm)(t)(x-t)αm+1-nmdt)=0.
The problem is to find y(x), x∈[0,L], satisfying (2.26) and (1.2).The generally nonlinear SVIDE in (2.26) is given in subinterval Ik, k=1,2,…,K as follows:
(2.27)F(x,yk(x),1Γ(n1-α1)[∑s=1k-1‍Λα1,s(x)+∫(k-1)hx‍yk(n1)(t)(x-t)α1+1-n1dt],…,1Γ(nm-αm)[∑s=1k-1‍Λαm,s(x)+∫(k-1)hx‍yk(nm)(t)(x-t)αm+1-nmdt])=0,
where x∈Ik and
(2.28)Λαq,s(x)=∫(s-1)hsh‍ys(nq)(t)(x-t)αq+1-nqdt,q=1,2,…,m;s=1,2,…,k-1.
It is important to note that, at the first subinterval, the summations in (2.27) are automatically discarded. For approximating the functions Λαq,s(x),q=1,…,m;s=1,…,k-1, with the aid of (2.17) and the Gaussian integration formula in the subinterval Is, given by (2.9), we obtain
(2.29)Λαq,s(x)≃h2∑j=1N-1‍wjYsT·Ls(nq)(xsj)(x-xsj)αq+1-nq,
where the ShLG collocation points xsj on the subinterval Is are defined by (2.10).From (2.27), (2.16), and (2.17), for k=1,2,…,K we have
(2.30)F(x,YkT·Lk(x),1Γ(n1-α1)[∑s=1k-1‍Λα1,s(x)+∫xk0x‍YkT·Lk(n1)(t)(x-t)α1+1-n1dt],…,1Γ(nm-αm)[∑s=1k-1‍Λαm,s(x)+∫xk0x‍YkT·Lk(nm)(t)(x-t)αm+1-nmdt])=0,x∈Ik.
We now collocate (2.30) at collocation points xkj, k=1,2,…,K and j=1,…,N-l as
(2.31)F(xkj,YkT·Lk(xkj),1Γ(n1-α1)[∑s=1k-1‍Λα1,s(xkj)+∫xk0xkj‍YkT·Lk(n1)(t)(xkj-t)α1+1-n1dt],…,1Γ(nm-αm)[∑s=1k-1‍Λαm,s(xkj)+∫xk0xkj‍YkT·Lk(nm)(t)(xkj-t)αm+1-nmdt])=0.
The integrals involved in (2.31) are singular. In order to convert them into nonsingular integrals, using integration by parts and with the aid of (2.10) we obtain
(2.32)F([∑s=1k-1‍Λαm,s(xkj)+1nm-αm(x1jnm-αmYkT·Lk(nm)(xk0)+∫xk0xkj‍YkT·Lk(nm+1)(t)(xkj-t)αm-nmdt)]xkj,ykj,1Γ(n1-α1)×[∑s=1k-1‍Λα1,s(xkj)+1n1-α1(x1jn1-α1YkT·Lk(n1)(xk0)+∫xk0xkj‍YkT·Lk(n1+1)(t)(xkj-t)α1-n1dt)],…,1Γ(nm-αm)×[∑s=1k-1‍Λαm,s(xkj)+1nm-αm(x1jnm-αmYkT·Lk(nm)(xk0)+∫xk0xkj‍YkT·Lk(nm+1)(t)(xkj-t)αm-nmdt)])=0.
In order to use the Gaussian integration formula in the subinterval Ik for (2.32), we transfer the t-interval [xk0,xkj] into the τ-interval Ik by means of the transformation τ=(h/x1j)(t-xk0)+xk0. Using this transformation, the Gaussian integration formula and (2.10), we have
(2.33)∫xk0xkj‍YkT·Lk(nq+1)(t)(xkj-t)αq-nqdt=x1jh∫xk0xkN‍YkT·Lk(nq+1)((x1j/h)(τ-xk0)+xk0)(x1j-(x1j/h)(τ-xk0))αq-nqdτ≃x1j2∑p=1N-1‍wpYkT·Lk(nq+1)((x1jx1p/h)+xk0)(x1j-(x1jx1p/h))αq-nq:=Gq(x1j).
By (2.33), (2.32) may be approximated as
(2.34)F(xkj,ykj,1Γ(n1-α1)[∑s=1k-1‍Λα1,s(xkj)+1n1-α1(x1jn1-α1YkT·Lk(n1)(xk0)+G1(x1j))],…,1Γ(nm-αm)[∑s=1k-1‍Λαm,s(xkj)+1nm-αm(x1jnm-αmYkT·Lk(nm)(xk0)+Gm(x1j))])=0.
In addition, substituting (2.16) and (2.17) into the boundary conditions (1.2) yields
(2.35)Hr(Yρ0T·Lρ0(ξ0),…,Yρ0T·Lρ0(l)(ξ0),…,YρlT·Lρl(ξl),…,YρlT·Lρl(l)(ξl))=0,r=0,1,…,l,
where ξr∈Iρr. Besides, it is required that the approximate solution and its first l derivatives be continuous at the interface of subintervals, that is,
(2.36)YkT·Lk(r)(xkN)=Yk+1T·Lk+1(r)(xk+1,0),k=1,2,…,K-1,r=0,1,…,l.Equation (2.34) for k=1,…,K, j=1,…,N-l together with (2.35)-(2.36) gives a system of equations with K(N+1) set of algebraic equations, which can be solved to find the unknowns of the vectors Yk,k=1,2,…,K. Consequently, the unknown functions y(x)|Ik≃yk(x) given in (2.16) can be calculated.


## 2.1. Review of Legendre and Shifted Legendre Polynomials

The Legendre polynomials,Pi(z), i=0,1,2,…, are the eigenfunctions of the singular Sturm-Liouville problem
(2.1)[(1-z2)Pi′(z)]′+i(i+1)Pi(z)=0.
Also, they are orthogonal with respect to L2 inner product on the interval [-1,1] with the weight function w(z)=1, that is,
(2.2)∫-11‍Pi(z)Pj(z)dz=22i+1δij,
where δij is the Kronecker delta. The Legendre polynomials satisfy the recursion relation
(2.3)Pi+1(z)=2i+1i+1zPi(z)-ii+1Pi-1(z),
where P0(z)=1 and P1(z)=z. If Pi(z) is normalized so that Pi(1)=1, then, for any i, the Legendre polynomials in terms of power of z are
(2.4)Pi(z)=12i∑k=0[i/2]‍(-1)k(ik)(2i-2ki)zi-2k,
where [i/2] denotes the integer part of i/2.The Legendre-Gauss (LG) collocation points-1<z1<z2<⋯<zN-1<1 are the roots of PN-1(z). No explicit formulas are known for the LG points; however, they are computed numerically using existing subroutines. The LG points have the property that
(2.5)∫-11‍p(z)dz=∑j=1N-1‍wjp(zj),
is exact for polynomials of degree at most 2N-3, where
(2.6)wj=2(1-zj2)[PN-1′(zj)]2,j=1,2,…,N-1,
are LG quadrature weights. For more details about Legendre polynomials, see [24].The shifted Legendre polynomials on the intervalx∈[a,b] are defined by
(2.7)P^i(x)=Pi(1b-a(2x-a-b)),i=0,1,2,…,
which are obtained by an affine transformation from the Legendre polynomials. The set of shifted Legendre polynomials is a complete L2[a,b]-orthogonal system with the weight function w(x)=1. Thus, any function f∈L2[a,b] can be expanded in terms of shifted Legendre polynomials.The ShLG collocation pointsa<x1<x2<⋯<xN-1<b on the interval [a,b] are obtained by shifting the LG points, zj, using the transformation
(2.8)xj=12((b-a)zj+a+b),j=1,2,…,N-1.
Thanks to the property of the standard LG quadrature, it follows that for any polynomial p of degree at most 2N-3 on (a,b),
(2.9)∫ab‍p(x)dx=b-a2∫-11‍p(12[(b-a)z+a+b])dz=b-a2∑j=1N-1‍wjp(12[(b-a)zj+a+b])=∑j=1N-1‍w^jp(xj),
where w^j=((b-a)/2)wj, 1⩽j⩽N-1 are ShLG quadrature weights.


## 2.2. Function Approximation

Suppose that the interval[0,L] is divided into K subintervals Ik=[(k-1)h,kh], k=1,2,…,K, where h=L/K. Let yk(x) be the solution of the problem in (1.1)-(1.2) in the subinterval Ik. Consider now the ShLG collocation points (k-1)h<xk1<⋯<xk,N-1<kh on the kth subinterval Ik, k=1,2,…,K, obtained using (2.8). Obviously,
(2.10)xkj=h2(zj+2k-1),j=1,2,…,N-1.
Also, consider two additional noncollocated points xk0=(k-1)h and xkN=kh. Let us define
(2.11)ℜkN=Span{Lk0(x),Lk1(x),…,LkN(x)},x∈Ik,
where
(2.12)Lki(x)=∏l=0,l≠iN‍x-xklxki-xkl,i=0,1,…,N,
is a basis of Lagrange interpolating polynomials on the subinterval Ik that satisfy Lki(xkj)=δij, where δij is the Kronecker delta function. The L2(Ik)-orthogonal projection IN:L2(Ik)→ℜkN is a mapping in a way that for any yk∈L2(Ik)(2.13)〈IN(yk)-yk,ϕk〉=0,∀ϕk∈ℜkN,
or equivalently
(2.14)IN(yk)(x)=∑i=0N‍ykiLki(x),x∈Ik,
where yki=yk(xki).Here, it can be easily seen that fori=0,1,…,N and k=1,2,…,K, we have
(2.15)Lki(x)=L1i(x-xk0),x∈Ik.
Thus, by utilizing (2.15) for (2.14), the approximation of yk(x) within each subinterval Ik can be restated as
(2.16)yk(x)≈IN(yk)(x)=∑i=0N‍ykiL1i(x-xk0)=YkT·Lk(x),x∈Ik,
where Yk and Lk(x) are (N+1)×1 matrices given by Yk=[yk0,…,ykN]T and Lk(x)=[L10(x-xk0),…,L1N(x-xk0)]T. It is important to observe that the series (2.16) includes the Lagrange polynomials associated with the noncollocated points xk0=(k-1)h and xkN=kh. Moreover, it is seen from (2.15)-(2.16) that, in the present adaptive scheme, it is only needed to produce the basis of Lagrange polynomials L1i(x) at the first subinterval.Bynm+1 times (nm is defined in (1.3)) differentiating of (2.16), we obtain
(2.17)yk(r)(x)≈YkT·Lk(r)(x),x∈Ik,
where Lk(r)(x)=(dr/dxr)Lk(x).


## 2.3. Convergence Rate

ForN⩾1 we introduce the piecewise polynomials space
(2.18)ΨIkN={y∈C0([0,L]):yk=y∣Ik∈ℙN(Ik)},
which is the space of the continuous functions over [0,L] whose restrictions on each subinterval Ik are polynomials of degree ⩽N. Then, for any continuous function y in [0,L], the piecewise interpolation polynomial ψN(y) coincides on each subinterval Ik with the interpolating polynomial IN(y) of yk=y|Ik at the ShLG points.In [25], with the aid of the formulas (5.4.33), (5.4.34) of [24], we prove the convergence properties of piecewise interpolation polynomial based on shifted Legendre-Gauss-Radau points in the norms of the Sobolev spaces. Accordingly, the following results for the convergence based on ShLG points hold.Theorem 2.1.
Suppose thaty∈Hv(0,L) with v⩾1. Then
(2.19)∥y-ψN(y)∥L2(0,L)⩽cN-v|y|H0;v;N;h(0,L),
and, for 1⩽u⩽v, if h⩽1, then
(2.20)∥y-ψN(y)∥Hu(0,L)⩽cN2u-(1/2)-v|y|Hu;v;N;h(0,L),
and if h>1, then
(2.21)∥y-ψN(y)∥Hu(0,L)⩽cN2u-(1/2)-v|y|H0;v;N;h(0,L).Note thatc denotes a positive constant that depends on v, but which is independent of the function y and integer N. Moreover, we introduce the seminorm of Hv(0,L), 0⩽u⩽v, N⩾0, h>0, as
(2.22)|y|Hu;v;N;h(0,L)=(∑l=min{v,N+1}v‍h2l-2u∥y(l)∥L2(0,L)2)1/2.Remark 2.2.
WheneverN⩾v-1, using (2.19)–(2.22), we get
(2.23)∥y-ψN(y)∥L2(0,L)⩽cN-vhv∥y(v)∥L2(0,L),
and, for u⩾1, if h⩽1, then
(2.24)∥y-ψN(y)∥Hu(0,L)⩽cN2u-(1/2)-vhv-u∥y(v)∥L2(0,L),
and if h>1, then
(2.25)∥y-ψN(y)∥Hu(0,L)⩽cN2u-(1/2)-vhv∥y(v)∥L2(0,L).
Equations (2.23)–(2.25) show that if y is infinitely smooth on [0,L] and h⩽1, the convergence rate of ψN(y) to y is faster than h to the power of N+1-u and any power of 1/N, which is superior to that for the global collocation method over [0,L]. Thus, the bigger the subinterval length the slower the convergence rate.


## 2.4. Problem Replacement and the Solution Technique

Consider the multi-term FBVP in (1.1)-(1.2). With substituting the definition of the Caputo-derivative (1.3) into (1.1), we can convert (1.1) into an equivalent SVIDE as
(2.26)F(x,y(x),1Γ(n1-α1)∫0x‍y(n1)(t)(x-t)α1+1-n1dt,…,1Γ(nm-αm)∫0x‍y(nm)(t)(x-t)αm+1-nmdt)=0.
The problem is to find y(x), x∈[0,L], satisfying (2.26) and (1.2).The generally nonlinear SVIDE in (2.26) is given in subinterval Ik, k=1,2,…,K as follows:
(2.27)F(x,yk(x),1Γ(n1-α1)[∑s=1k-1‍Λα1,s(x)+∫(k-1)hx‍yk(n1)(t)(x-t)α1+1-n1dt],…,1Γ(nm-αm)[∑s=1k-1‍Λαm,s(x)+∫(k-1)hx‍yk(nm)(t)(x-t)αm+1-nmdt])=0,
where x∈Ik and
(2.28)Λαq,s(x)=∫(s-1)hsh‍ys(nq)(t)(x-t)αq+1-nqdt,q=1,2,…,m;s=1,2,…,k-1.
It is important to note that, at the first subinterval, the summations in (2.27) are automatically discarded. For approximating the functions Λαq,s(x),q=1,…,m;s=1,…,k-1, with the aid of (2.17) and the Gaussian integration formula in the subinterval Is, given by (2.9), we obtain
(2.29)Λαq,s(x)≃h2∑j=1N-1‍wjYsT·Ls(nq)(xsj)(x-xsj)αq+1-nq,
where the ShLG collocation points xsj on the subinterval Is are defined by (2.10).From (2.27), (2.16), and (2.17), for k=1,2,…,K we have
(2.30)F(x,YkT·Lk(x),1Γ(n1-α1)[∑s=1k-1‍Λα1,s(x)+∫xk0x‍YkT·Lk(n1)(t)(x-t)α1+1-n1dt],…,1Γ(nm-αm)[∑s=1k-1‍Λαm,s(x)+∫xk0x‍YkT·Lk(nm)(t)(x-t)αm+1-nmdt])=0,x∈Ik.
We now collocate (2.30) at collocation points xkj, k=1,2,…,K and j=1,…,N-l as
(2.31)F(xkj,YkT·Lk(xkj),1Γ(n1-α1)[∑s=1k-1‍Λα1,s(xkj)+∫xk0xkj‍YkT·Lk(n1)(t)(xkj-t)α1+1-n1dt],…,1Γ(nm-αm)[∑s=1k-1‍Λαm,s(xkj)+∫xk0xkj‍YkT·Lk(nm)(t)(xkj-t)αm+1-nmdt])=0.
The integrals involved in (2.31) are singular. In order to convert them into nonsingular integrals, using integration by parts and with the aid of (2.10) we obtain
(2.32)F([∑s=1k-1‍Λαm,s(xkj)+1nm-αm(x1jnm-αmYkT·Lk(nm)(xk0)+∫xk0xkj‍YkT·Lk(nm+1)(t)(xkj-t)αm-nmdt)]xkj,ykj,1Γ(n1-α1)×[∑s=1k-1‍Λα1,s(xkj)+1n1-α1(x1jn1-α1YkT·Lk(n1)(xk0)+∫xk0xkj‍YkT·Lk(n1+1)(t)(xkj-t)α1-n1dt)],…,1Γ(nm-αm)×[∑s=1k-1‍Λαm,s(xkj)+1nm-αm(x1jnm-αmYkT·Lk(nm)(xk0)+∫xk0xkj‍YkT·Lk(nm+1)(t)(xkj-t)αm-nmdt)])=0.
In order to use the Gaussian integration formula in the subinterval Ik for (2.32), we transfer the t-interval [xk0,xkj] into the τ-interval Ik by means of the transformation τ=(h/x1j)(t-xk0)+xk0. Using this transformation, the Gaussian integration formula and (2.10), we have
(2.33)∫xk0xkj‍YkT·Lk(nq+1)(t)(xkj-t)αq-nqdt=x1jh∫xk0xkN‍YkT·Lk(nq+1)((x1j/h)(τ-xk0)+xk0)(x1j-(x1j/h)(τ-xk0))αq-nqdτ≃x1j2∑p=1N-1‍wpYkT·Lk(nq+1)((x1jx1p/h)+xk0)(x1j-(x1jx1p/h))αq-nq:=Gq(x1j).
By (2.33), (2.32) may be approximated as
(2.34)F(xkj,ykj,1Γ(n1-α1)[∑s=1k-1‍Λα1,s(xkj)+1n1-α1(x1jn1-α1YkT·Lk(n1)(xk0)+G1(x1j))],…,1Γ(nm-αm)[∑s=1k-1‍Λαm,s(xkj)+1nm-αm(x1jnm-αmYkT·Lk(nm)(xk0)+Gm(x1j))])=0.
In addition, substituting (2.16) and (2.17) into the boundary conditions (1.2) yields
(2.35)Hr(Yρ0T·Lρ0(ξ0),…,Yρ0T·Lρ0(l)(ξ0),…,YρlT·Lρl(ξl),…,YρlT·Lρl(l)(ξl))=0,r=0,1,…,l,
where ξr∈Iρr. Besides, it is required that the approximate solution and its first l derivatives be continuous at the interface of subintervals, that is,
(2.36)YkT·Lk(r)(xkN)=Yk+1T·Lk+1(r)(xk+1,0),k=1,2,…,K-1,r=0,1,…,l.Equation (2.34) for k=1,…,K, j=1,…,N-l together with (2.35)-(2.36) gives a system of equations with K(N+1) set of algebraic equations, which can be solved to find the unknowns of the vectors Yk,k=1,2,…,K. Consequently, the unknown functions y(x)|Ik≃yk(x) given in (2.16) can be calculated.


## 3. Numerical Examples

In this section we give the computational results of numerical experiments with the method based on preceding sections to support our theoretical discussion.Example 3.1.
In this example, we consider the Bagley-Torvik equation [26]
(3.1)Ay′′(x)+BDαy(x)+Cy(x)=f(x),x∈[0,L],
with boundary conditions
(3.2)y(0)=c0,y(L)=c1,
where A,B,C∈ℝ and A≠0. Bagley-Torvik equation involving fractional derivative of order 1/2 or 3/2 arises in the modeling of the motion of a rigid plate in a Newtonian fluid and a gas in a fluid. Since the Bagley-Torvik equation is a prototype fractional differential equation with two derivatives and represents a general form of the fractional problems, its solution can give many ideas about the solution of similar problems in fractional differential equations. Podlubny [2] has investigated the solution of Bagley-Torvik equation (3.1) and for α=3/2 gave the analytical solution with homogeneous initial conditions by using Green’s function, as follows:
(3.3)y(x)=∫0x‍G3(x-t)f(t)dt,G3(x)=1A∑k=0∞‍(-1)kk!(CA)kx2k+1E1/2,2+(3k/2)(k)(-BAx),Eλ,μ(k)(z)=dkdzkEλ,μ(z)=∑j=0∞‍(j+k)!zjj!Γ(λj+λk+μ),k=0,1,2,…,
where Eλ,μ is the Mittag-Leffler function in two parameters and the G3 three-term Green’s function. However, in practice, these equations can not be evaluated easily for different functions f(x). Several other authors have proposed different techniques for the solution of this equation. A review of the solution techniques for Bagley-Torvik equation can be found in [27].
Here, we solve (3.1) with two-point boundary conditions (3.2) by using the adaptive pseudospectral method. For comparison purposes and in order to demonstrate the efficiency of our method, we investigate the following cases. Further, for completeness, the derivation of Bagley-Torvik equation is given in the appendix.Case 1.
In (3.1)-(3.2) set α=3/2, A=B=C=1, f(x)=x2+2+4x/π, L=5, c0=0, and c1=25. It is readily verified that the exact solution of this case is y(x)=x2. Using the adaptive pseudospectral method in Section 2 with K=1 and N=2, the unknowns yki, k=1,…,K, i=0,…,N in (2.16) are found to be
(3.4)y10=0,y11=6.25,y12=5,
which lead to the exact solution y(x)=x2. This case was solved in [6] using a collocation-shooting method. Their computed maximum absolute error and L2 error norm were 2.00×10-14 and 3.78×10-12, respectively, which show that our method is more efficient.Case 2.
Setα=3/2, A=B=C=1, f(x)=(15/4)x+(15/8)πx+x5/2+1, L=1, c0=1, and c1=2. The exact solution of this case, which was considered in [5], is y(x)=x5/2+1. In Table 1 the maximum absolute errors for different values of K and N are presented. We see from Table 1 that, as stated in Section 2.3, the more rapid convergence rate is obtained with smaller subinterval length.Table 1
Maximum absolute errors for Example3.1, Case 2.
N
=
10
N
=
20
N
=
30
N
=
40
K
=
1
3.5
×
10
-
5
3.8
×
10
-
6
1.5
×
10
-
6
4.4
×
10
-
7
K
=
2
6.9
×
10
-
6
7.1
×
10
-
7
2.0
×
10
-
7
8.4
×
10
-
8Case 3.
For comparison, the same coefficients as considered in [27] have been used here. Set α=3/2, A=1, B=C=1/2, f(x)=8, L=1, c0=0, and c1=2.95258388. Table 2 shows the comparison of solutions of this case by the present method (with K=2, N=40), GTC method [27] and the analytical solution [2], and the good agreement of our adaptive pseudospectral solution with analytical solution.Table 2
Comparison of solutions for Example3.1, Case 3.
x
GTCM [27]
Present method
Analytical [2]
0
0
0
0
0.1
0.03648555
0.03648741
0.03648748
0.2
0.14063472
0.14063951
0.14063962
0.3
0.30747623
0.30748449
0.30748463
0.4
0.53327129
0.53328396
0.53328411
0.5
0.81473561
0.81475679
0.81475695
0.6
1.14880581
1.14883734
1.14883742
0.7
1.53252126
1.53256541
1.53256543
0.8
1.96297499
1.96302931
1.96302925
0.9
2.43745598
2.43733391
2.43733397
1
2.95407000
2.95258388
2.95258388Case 4.
Setα=1/2, A=B=C=1, f(x)=8, L=1, c0=0, and c1=3.10190571. The numerical solutions obtained by the present method (with K=2, N=40), fractional finite difference method (FDM), the Adomian decomposition method (ADM), and the variational iteration method (VIM) from [28] are given in Table 3. The exact solution refers to the closed form series solution given in [28]. Table 3 shows the excellent agreement of our adaptive pseudospectral solution with the exact solution.Table 3
Comparison of solutions for Example3.1, Case 4.
x
FDM [28]
ADM [28]
VIM [28]
Present method
Exact
0
0
0
0
0
0
0.1
0.039473
0.039874
0.039874
0.03975004
0.03975003
0.2
0.157703
0.158512
0.158512
0.15703584
0.15703582
0.3
0.352402
0.353625
0.353625
0.34736999
0.34736998
0.4
0.620435
0.622083
0.622083
0.60469514
0.60469515
0.5
0.957963
0.960047
0.960047
0.92176757
0.92176764
0.6
1.360551
1.363093
1.363093
1.29045651
1.29045656
0.7
1.823267
1.826257
1.826257
1.70200794
1.70200797
0.8
2.340749
2.344224
2.344224
2.14728692
2.14728693
0.9
2.907324
2.911278
2.911278
2.61700100
2.61700101
1
3.517013
3.521462
3.521462
3.10190571
3.10190571Example 3.2.
As a multi-term equation, consider the linear multi-term FBVP described by(3.5)πxD5/2y(x)+10Γ(23)x3D4/3y(x)+Γ(34)x4D1/4y(x)-4077y(x)=12x+54x2+87x3,y(0)=y′(0)=0,y(1)=1.
The exact solution to this problem is y(x)=x3. Since this problem is a third-order equation, it can demonstrate the effect of the continuity conditions (2.36) on the approximate solution. Table 4 compares the maximum absolute errors obtained using the present method for different values of K and N with the errors reported in [13] using operational matrix of fractional derivatives using B-spline functions. Note that in [13], for each value of J, the obtained algebraic system is of order 2J+1, while in the present method the obtained algebraic system is of order K(N+1). It is important to see that our method provides more accurate results with solving lower-order algebraic systems. Further, it is seen that in the present method the accuracy can be improved either by increasing the number of subintervals or by increasing the number of collocation points within each subinterval.Table 4
Comparison of maximum absolute errors for Example3.2.
Present method
Method [13]
N
=
4
N
=
8
N
=
16
N
=
24
J
Error
K
=
1
6.9
×
10
-
4
6.4
×
10
-
5
4.1
×
10
-
6
9.9
×
10
-
7
4
1.5
×
10
-
3
K
=
2
2.2
×
10
-
4
1.2
×
10
-
5
1.3
×
10
-
6
2.5
×
10
-
7
5
6.1
×
10
-
4
K
=
3
8.1
×
10
-
5
4.7
×
10
-
6
4.0
×
10
-
7
9.8
×
10
-
8
6
1.8
×
10
-
4
K
=
4
4.5
×
10
-
5
2.8
×
10
-
6
2.4
×
10
-
7
5.9
×
10
-
8
7
7.2
×
10
-
5Example 3.3.
Consider the nonlinear multi-term FBVP described by(3.6)y′′(x)+Γ(45)x65D6/5y(x)+119Γ(56)x6D1/6y(x)-(y′(x))2=2+110x2,y(0)=1,y(1)=2.
The exact solution to this problem is y(x)=x2+1. In Table 5, we compare the maximum absolute errors obtained using the present adaptive method for different values of K and N with the errors reported in [13] using operational matrix of fractional derivatives using B-spline functions.Table 5
Comparison of maximum absolute errors for Example3.3.
Present method
Method [13]
N
=
10
N
=
20
N
=
30
N
=
40
J
Error
K
=
1
1.3
×
10
-
6
8.9
×
10
-
8
2.0
×
10
-
8
6.7
×
10
-
9
4
1.2
×
10
-
3
K
=
2
7.5
×
10
-
7
5.3
×
10
-
8
1.7
×
10
-
8
4.9
×
10
-
9
5
3.3
×
10
-
4
K
=
3
6.5
×
10
-
7
4.7
×
10
-
8
1.4
×
10
-
8
3.6
×
10
-
9
6
8.1
×
10
-
5
K
=
4
5.0
×
10
-
7
3.6
×
10
-
8
7.9
×
10
-
9
2.8
×
10
-
9
7
2.1
×
10
-
5Example 3.4.
Consider the following nonlinear FBVP:(3.7)Dαy(x)+ayn(x)=f(x),1<α⩽2,y(0)=c0,y(1)=c1.
For comparison, we choose n=2, a=e-2π, f(x)=(105π/16Γ((9/2)-α))x(7/2)-α+e-2πx7, c0=0, and c1=1. It is readily verified that the exact solution is y(x)=x7/2. In Table 6, the absolute errors obtained using the present adaptive pseudospectral method for K=4 and N=40 and different values of α are compared with the errors obtained in [16] using Legendre wavelets, which show that the present method provides more accurate numerical results.Table 6
Comparison of absolute errors for Example3.4.
α
x
=
0.1
x
=
0.3
x
=
0.5
x
=
0.7
x
=
0.9
1.1
Present method
2.7
×
10
-
10
5.0
×
10
-
10
2.1
×
10
-
9
6.4
×
10
-
10
1.3
×
10
-
9
Method [16]
2.9
×
10
-
4
6.0
×
10
-
3
8.4
×
10
-
3
5.8
×
10
-
3
3.4
×
10
-
3
1.3
Present method
2.6
×
10
-
9
9.0
×
10
-
9
2.4
×
10
-
8
1.5
×
10
-
8
9.4
×
10
-
9
Method [16]
2.0
×
10
-
4
3.0
×
10
-
3
4.5
×
10
-
3
3.0
×
10
-
3
4.0
×
10
-
3
1.5
Present method
6.3
×
10
-
9
3.8
×
10
-
8
1.0
×
10
-
7
7.9
×
10
-
8
2.9
×
10
-
8
Method [16]
9.7
×
10
-
5
1.5
×
10
-
3
2.4
×
10
-
3
1.5
×
10
-
3
4.7
×
10
-
3
1.7
Present method
3.1
×
10
-
8
1.6
×
10
-
8
2.0
×
10
-
7
1.8
×
10
-
7
7.9
×
10
-
8
Method [16]
4.8
×
10
-
5
7.5
×
10
-
4
1.3
×
10
-
3
6.0
×
10
-
4
5.4
×
10
-
3
1.9
Present method
2.1
×
10
-
7
3.7
×
10
-
7
2.0
×
10
-
7
1.2
×
10
-
7
2.2
×
10
-
7
Method [16]
2.9
×
10
-
5
3.7
×
10
-
4
6.8
×
10
-
4
9.2
×
10
-
5
6.2
×
10
-
3
2.0
Present method
7.0
×
10
-
14
9.0
×
10
-
14
6.8
×
10
-
14
3.7
×
10
-
13
2.1
×
10
-
13
Method [16]
2.4
×
10
-
5
2.6
×
10
-
4
4.9
×
10
-
4
7.0
×
10
-
5
6.6
×
10
-
3Example 3.5.
In this example, to show the applicability of the present method for larger interval, we consider the nonlinear FBVP described by(3.8)Dαy(x)+y2(x)=Eα2(-xα)-Eα(-xα),1<α⩽2,y(0)=1,y(10)=Eα(-10α).
The exact solution of this problem is given by y(x)=Eα(-xα), where Eα(z)=∑k=0∞‍(zk/Γ(αk+1)) is the Mittag-Leffler function.
In Table7, the maximum absolute errors in the interval [0,10] for α=1.75 and different values of K and N are presented, which shows the efficiency of the present method for FBVPs in large domains. Also, the numerical results for y(x) by adaptive pseudospectral method for K=20, N=30 and α=1.25,1.5,1.75,1.95, and 2 together with the exact solutions are plotted in Figure 1, which indicates that the numerical results are in high agreement with the exact ones. Moreover, Figure 1 demonstrates the efficiency of the present method for solutions having oscillatory behavior. For α=2, the exact solution is given as y(x)=cos(x). Note that as α approaches 2, the numerical solution converges to the analytical solution; that is, in the limit, the solution of the fractional differential equations approaches to that of the integer-order differential equations.Table 7
Maximum absolute errors forα=1.75 for Example 3.5.
N
=
10
N
=
20
N
=
30
K
=
5
3.5
×
10
-
3
8.1
×
10
-
4
3.5
×
10
-
4
K
=
10
1.7
×
10
-
3
3.9
×
10
-
4
1.8
×
10
-
4
K
=
20
9.4
×
10
-
4
2.1
×
10
-
4
9.1
×
10
-
5Figure 1
Comparison ofy(x) for K=20, N=30 with exact solutions for Example 3.5.Example 3.6.
Finally consider the nonlinear multi-term FBVP described by(3.9)ay′′(x)+bDα2y(x)+c(Dα1y(x))2+ey3(x)=f(x),0⩽x⩽2,y(0)=0,y(2)=83,
where a,b,c,e∈ℝ, 0<α1⩽1, 1<α2⩽2 and f(x)=2ax+(2b/Γ(4-α2))x3-α2+c((2/Γ(4-α1))x3-α1)2+e((1/3)x3)3. The exact solution to this problem is y(x)=(1/3)x3.
Fora=b=c=e=1, α1=0.555, and α2=1.455 the maximum absolute errors obtained using the adaptive pseudospectral method are given in Table 8. Also, for a=0.1, b=c=e=0.5, α1=0.219, and α2=1.965 the maximum absolute errors are given in Table 9. Again, it is seen that in the present adaptive pseudospectral method the accuracy is improved either by increasing the number of subintervals or by increasing the number of collocation points within each subinterval.Table 8
Maximum absolute errors fora=b=c=e=1,α1=0.555,α2=1.455 for Example 3.6.
N
=
10
N
=
20
N
=
30
K
=
1
1.4
×
10
-
3
1.7
×
10
-
4
5.0
×
10
-
5
K
=
2
3.3
×
10
-
4
4.1
×
10
-
5
8.3
×
10
-
6
K
=
4
2.2
×
10
-
5
3.0
×
10
-
6
9.1
×
10
-
7Table 9
Maximum absolute errors fora=0.1,b=c=e=0.5,α1=0.219,α2=1.965 for Example 3.6.
N
=
10
N
=
20
N
=
30
K
=
1
1.9
×
10
-
2
1.7
×
10
-
3
6.6
×
10
-
4
K
=
2
2.4
×
10
-
3
4.2
×
10
-
4
1.7
×
10
-
4
K
=
4
3.2
×
10
-
4
5.8
×
10
-
5
2.4
×
10
-
5


## 4. Conclusion

In this work a new adaptive pseudospectral method based on ShLG collocation points has been proposed for solving the multi-term FBVPs. We converted the original FBVP into a SIVDE and then reduced it to a system of algebraic equations using collocation. The difficulty in SIVDE, due to the singularity, is overcome here by utilizing integration by parts. By considering some additional conditions, the continuity of the approximate solution and its firstl derivatives is kept. It was also shown that the accuracy can be improved either by increasing the number of subintervals or by increasing the number of collocation points in subintervals. Moreover, this method is valid for large-domain calculations. The achieved results are compared with exact solutions and with the solutions obtained by some other numerical methods, which demonstrate the convergence, validity, and accuracy of the proposed method.



---

*Source: 381708-2012-12-17.xml*