An Adaptive Pseudospectral Method for Fractional Order Boundary Value Problems

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.

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𝑦(𝑥),…,𝐷^𝛼𝑚𝑦(𝑥)=0,𝑥∈[0,𝐿], (1.1) subject to 𝐻_𝑟𝑦𝜉₀,…,𝑦^(𝑙)𝜉₀,…,𝑦𝜉_𝑙,…,𝑦^(𝑙)𝜉_𝑙=0,𝑟=0,1,…,𝑙, (1.2) where 𝐹 can be nonlinear in general, 0<𝛼₁<𝛼₂<⋯<𝛼_𝑚, 𝑙<𝛼_𝑚≤𝑙+1, 𝐿∈ℝ, 𝐻_𝑟 are linear functions, the points 𝜉₀,𝜉₁, …,𝜉_𝑙 lie in [0,𝐿], and 𝐷^𝛼𝑞 denotes the Caputo-fractional derivative of order 𝛼_𝑞, defined as follows [23]: 𝐷^𝛼𝑞𝑦(𝑥)=1Γ𝑛_𝑞−𝛼_𝑞^𝑥0𝑦^(𝑛𝑞)(𝑡)(𝑥−𝑡)^{𝛼𝑞+1−𝑛𝑞}d𝑡,𝑛_𝑞=𝛼_𝑞+1,𝑞=1,2,…,𝑚, (1.3) where [𝛼_𝑞] denotes the integer part of the real number 𝛼_𝑞. 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 first 𝑙 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 Section 2, 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, 𝑃_𝑖(𝑧), 𝑖=0,1,2,…, are the eigenfunctions of the singular Sturm-Liouville problem 1−𝑧²𝑃^′_𝑖(𝑧)^′+𝑖(𝑖+1)𝑃_𝑖(𝑧)=0. (2.1) Also, they are orthogonal with respect to 𝐿² inner product on the interval [−1,1] with the weight function 𝑤(𝑧)=1, that is, ¹₋₁𝑃_𝑖(𝑧)𝑃_𝑗(𝑧)2d𝑧=𝛿2𝑖+1^𝑖𝑗, (2.2) where 𝛿_𝑖𝑗 is the Kronecker delta. The Legendre polynomials satisfy the recursion relation 𝑃_𝑖+1(𝑧)=2𝑖+1𝑖+1𝑧𝑃_𝑖(𝑧)−𝑖𝑃𝑖+1^𝑖−1(𝑧), (2.3) where 𝑃₀(𝑧)=1 and 𝑃₁(𝑧)=𝑧. If 𝑃_𝑖(𝑧) is normalized so that 𝑃_𝑖(1)=1, then, for any 𝑖, the Legendre polynomials in terms of power of 𝑧 are 𝑃_𝑖(𝑧)=12^𝑖[𝑖/2]_𝑘=0(−1)^𝑘𝑖𝑘𝑖𝑧2𝑖−2𝑘_𝑖−2𝑘, (2.4) where [𝑖/2] denotes the integer part of 𝑖/2.

The Legendre-Gauss (LG) collocation points −1<𝑧₁<𝑧₂<⋯<𝑧_𝑁−1<1 are the roots of 𝑃_𝑁−1(𝑧). No explicit formulas are known for the LG points; however, they are computed numerically using existing subroutines. The LG points have the property that ¹₋₁𝑝(𝑧)d𝑧=^𝑁−1_𝑗=1𝑤_𝑗𝑝𝑧_𝑗, (2.5) is exact for polynomials of degree at most 2𝑁−3, where 𝑤_𝑗=21−𝑧²_𝑗𝑃^′_𝑁−1𝑧_𝑗²,𝑗=1,2,…,𝑁−1, (2.6) are LG quadrature weights. For more details about Legendre polynomials, see [24].

The shifted Legendre polynomials on the interval 𝑥∈[𝑎,𝑏] are defined by 𝑃_𝑖(𝑥)=𝑃_𝑖1(𝑏−𝑎2𝑥−𝑎−𝑏),𝑖=0,1,2,…, (2.7) which are obtained by an affine transformation from the Legendre polynomials. The set of shifted Legendre polynomials is a complete 𝐿²[𝑎,𝑏]-orthogonal system with the weight function 𝑤(𝑥)=1. Thus, any function 𝑓∈𝐿²[𝑎,𝑏] can be expanded in terms of shifted Legendre polynomials.

The ShLG collocation points 𝑎<𝑥₁<𝑥₂<⋯<𝑥_𝑁−1<𝑏 on the interval [𝑎,𝑏] are obtained by shifting the LG points, 𝑧_𝑗, using the transformation 𝑥_𝑗=12(𝑏−𝑎)𝑧_𝑗+𝑎+𝑏,𝑗=1,2,…,𝑁−1. (2.8) Thanks to the property of the standard LG quadrature, it follows that for any polynomial 𝑝 of degree at most 2𝑁−3 on (𝑎,𝑏), ^𝑏_𝑎𝑝(𝑥)d𝑥=𝑏−𝑎2¹₋₁𝑝12)]=[(𝑏−𝑎𝑧+𝑎+𝑏d𝑧𝑏−𝑎2^𝑁−1_𝑗=1𝑤_𝑗𝑝12(𝑏−𝑎)𝑧_𝑗+𝑎+𝑏=^𝑁−1_𝑗=1𝑤_𝑗𝑝𝑥_𝑗, (2.9) where 𝑤_𝑗=((𝑏−𝑎)/2)𝑤_𝑗, 1⩽𝑗⩽𝑁−1 are ShLG quadrature weights.

2.2. Function Approximation

Suppose that the interval [0,𝐿] is divided into 𝐾 subintervals 𝐼_𝑘=[(𝑘−1)ℎ,𝑘ℎ], 𝑘=1,2,…,𝐾, where ℎ=𝐿/𝐾. Let 𝑦_𝑘(𝑥) be the solution of the problem in (1.1)-(1.2) in the subinterval 𝐼_𝑘. Consider now the ShLG collocation points (𝑘−1)ℎ<𝑥_𝑘1<⋯<𝑥_{𝑘,𝑁−1}<𝑘ℎ on the 𝑘th subinterval 𝐼_𝑘, 𝑘=1,2,…,𝐾, obtained using (2.8). Obviously, 𝑥_𝑘𝑗=ℎ2𝑧_𝑗+2𝑘−1,𝑗=1,2,…,𝑁−1. (2.10) Also, consider two additional noncollocated points 𝑥_𝑘0=(𝑘−1)ℎ and 𝑥_𝑘𝑁=𝑘ℎ. Let us define ℜ_𝑘𝑁=Span𝐿_𝑘0(𝑥),𝐿_𝑘1(𝑥),…,𝐿_𝑘𝑁(𝑥),𝑥∈𝐼_𝑘, (2.11) where 𝐿_𝑘𝑖(𝑥)=^𝑁_{𝑙=0,𝑙≠𝑖}𝑥−𝑥_𝑘𝑙𝑥_𝑘𝑖−𝑥_𝑘𝑙,𝑖=0,1,…,𝑁, (2.12) is a basis of Lagrange interpolating polynomials on the subinterval 𝐼_𝑘 that satisfy 𝐿_𝑘𝑖(𝑥_𝑘𝑗)=𝛿_𝑖𝑗, where 𝛿_𝑖𝑗 is the Kronecker delta function. The 𝐿²(𝐼_𝑘)-orthogonal projection 𝐼_𝑁∶𝐿²(𝐼_𝑘)→ℜ_𝑘𝑁 is a mapping in a way that for any 𝑦_𝑘∈𝐿²(𝐼_𝑘)𝐼_𝑁𝑦_𝑘−𝑦_𝑘,𝜙_𝑘=0,∀𝜙_𝑘∈ℜ_𝑘𝑁, (2.13) or equivalently 𝐼_𝑁𝑦_𝑘(𝑥)=^𝑁_𝑖=0𝑦_𝑘𝑖𝐿_𝑘𝑖(𝑥),𝑥∈𝐼_𝑘, (2.14) where 𝑦_𝑘𝑖=𝑦_𝑘(𝑥_𝑘𝑖).

Here, it can be easily seen that for 𝑖=0,1,…,𝑁 and 𝑘=1,2,…,𝐾, we have 𝐿_𝑘𝑖(𝑥)=𝐿_1𝑖𝑥−𝑥_𝑘0,𝑥∈𝐼_𝑘. (2.15) Thus, by utilizing (2.15) for (2.14), the approximation of 𝑦_𝑘(𝑥) within each subinterval 𝐼_𝑘 can be restated as 𝑦_𝑘(𝑥)≈𝐼_𝑁𝑦_𝑘(𝑥)=^𝑁_𝑖=0𝑦_𝑘𝑖𝐿_1𝑖𝑥−𝑥_𝑘0=𝑌^𝑇_𝑘⋅𝐿_𝑘(𝑥),𝑥∈𝐼_𝑘, (2.16) where 𝑌_𝑘 and 𝐿_𝑘(𝑥) are (𝑁+1)×1 matrices given by 𝑌_𝑘=[𝑦_𝑘0,…,𝑦_𝑘𝑁]^𝑇 and 𝐿_𝑘(𝑥)=[𝐿₁₀(𝑥−𝑥_𝑘0),…,𝐿_1𝑁(𝑥−𝑥_𝑘0)]^𝑇. It is important to observe that the series (2.16) includes the Lagrange polynomials associated with the noncollocated points 𝑥_𝑘0=(𝑘−1)ℎ and 𝑥_𝑘𝑁=𝑘ℎ. 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 𝐿_1𝑖(𝑥) at the first subinterval.

By 𝑛_𝑚+1 times (𝑛_𝑚 is defined in (1.3)) differentiating of (2.16), we obtain 𝑦_𝑘^(𝑟)(𝑥)≈𝑌^𝑇_𝑘⋅𝐿_𝑘^(𝑟)(𝑥),𝑥∈𝐼_𝑘, (2.17) where 𝐿_𝑘^(𝑟)(𝑥)=(𝑑^𝑟/𝑑𝑥^𝑟)𝐿_𝑘(𝑥).

2.3. Convergence Rate

For 𝑁⩾1 we introduce the piecewise polynomials space Ψ^𝑁_𝐼𝑘=𝑦∈𝐶⁰([0,𝐿])∶𝑦_𝑘=𝑦∣_𝐼𝑘∈ℙ_𝑁𝐼_𝑘, (2.18) which is the space of the continuous functions over [0,𝐿] whose restrictions on each subinterval 𝐼_𝑘 are polynomials of degree ⩽𝑁. Then, for any continuous function 𝑦 in [0,𝐿], the piecewise interpolation polynomial 𝜓_𝑁(𝑦) coincides on each subinterval 𝐼_𝑘 with the interpolating polynomial 𝐼_𝑁(𝑦) of 𝑦_𝑘=𝑦|_𝐼𝑘 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 that 𝑦∈𝐻^𝑣(0,𝐿) with 𝑣⩾1. Then 𝑦−𝜓_𝑁𝑦_{𝐿2(0,𝐿)}⩽𝑐𝑁^−𝑣𝑦_{𝐻0;𝑣;𝑁;ℎ(0,𝐿)}, (2.19) and, for 1⩽𝑢⩽𝑣, if ℎ⩽1, then 𝑦−𝜓_𝑁(𝑦)_{𝐻𝑢(0,𝐿)}⩽𝑐𝑁^{2𝑢−(1/2)−𝑣}𝑦_{𝐻𝑢;𝑣;𝑁;ℎ(0,𝐿)}, (2.20) and if ℎ>1, then 𝑦−𝜓_𝑁(𝑦)_{𝐻𝑢(0,𝐿)}⩽𝑐𝑁^{2𝑢−(1/2)−𝑣}𝑦_{𝐻0;𝑣;𝑁;ℎ(0,𝐿)}. (2.21)

Note that 𝑐 denotes a positive constant that depends on 𝑣, but which is independent of the function 𝑦 and integer 𝑁. Moreover, we introduce the seminorm of 𝐻^𝑣(0,𝐿), 0⩽𝑢⩽𝑣, 𝑁⩾0, ℎ>0, as 𝑦_{𝐻𝑢;𝑣;𝑁;ℎ(0,𝐿)}=^𝑣_{𝑙=min{𝑣,𝑁+1}}ℎ^{2𝑙−2𝑢}𝑦^(𝑙)2_{𝐿2()0,𝐿}^1/2. (2.22)

Remark 2.2. Whenever 𝑁⩾𝑣−1, using (2.19)–(2.22), we get 𝑦−𝜓_𝑁(𝑦)_{𝐿2(0,𝐿)}⩽𝑐𝑁^−𝑣ℎ^𝑣𝑦^(𝑣)_{𝐿2(0,𝐿)}, (2.23) and, for 𝑢⩾1, if ℎ⩽1, then 𝑦−𝜓_𝑁(𝑦)_{𝐻𝑢(0,𝐿)}⩽𝑐𝑁^{2𝑢−(1/2)−𝑣}ℎ^𝑣−𝑢𝑦^(𝑣)_{𝐿2(0,𝐿)}, (2.24) and if ℎ>1, then 𝑦−𝜓_𝑁(𝑦)_{𝐻𝑢(0,𝐿)}⩽𝑐𝑁^{2𝑢−(1/2)−𝑣}ℎ^𝑣𝑦^(𝑣)_{𝐿2(0,𝐿)}. (2.25) Equations (2.23)–(2.25) show that if 𝑦 is infinitely smooth on [0,𝐿] and ℎ⩽1, the convergence rate of 𝜓_𝑁(𝑦) to 𝑦 is faster than ℎ to the power of 𝑁+1−𝑢 and any power of 1/𝑁, which is superior to that for the global collocation method over [0,𝐿]. 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 𝐹(𝑥),1𝑥,𝑦Γ𝑛₁−𝛼₁^𝑥₀𝑦^(𝑛1)(𝑡)(𝑥−𝑡)^{𝛼1+1−𝑛1}1d𝑡,…,Γ𝑛_𝑚−𝛼_𝑚^𝑥₀𝑦^(𝑛𝑚)(𝑡)(𝑥−𝑡)^{𝛼𝑚+1−𝑛𝑚}d𝑡=0. (2.26) The problem is to find 𝑦(𝑥), 𝑥∈[0,𝐿], satisfying (2.26) and (1.2).

The generally nonlinear SVIDE in (2.26) is given in subinterval 𝐼_𝑘, 𝑘=1,2,…,𝐾 as follows: 𝐹𝑥,𝑦_𝑘(𝑥),1Γ𝑛₁−𝛼₁^𝑘−1_𝑠=1Λ_𝛼1,𝑠(𝑥)+^𝑥_{()ℎ𝑘−1}𝑦^(𝑛1)_𝑘(𝑡)(𝑥−𝑡)^{𝛼1+1−𝑛1}1d𝑡,…,Γ𝑛_𝑚−𝛼_𝑚^𝑘−1_𝑠=1Λ_{𝛼𝑚,𝑠}(𝑥)+^𝑥_{()ℎ𝑘−1}𝑦^(𝑛𝑚)_𝑘(𝑡)(𝑥−𝑡)^{𝛼𝑚+1−𝑛𝑚}d𝑡=0, (2.27) where 𝑥∈𝐼_𝑘 and Λ_{𝛼𝑞,𝑠}(𝑥)=_()ℎ^𝑠ℎ_𝑠−1𝑦^(𝑛𝑞)_𝑠(𝑡)(𝑥−𝑡)^{𝛼𝑞+1−𝑛𝑞}d𝑡,𝑞=1,2,…,𝑚;𝑠=1,2,…,𝑘−1. (2.28) It is important to note that, at the first subinterval, the summations in (2.27) are automatically discarded. For approximating the functions Λ_{𝛼𝑞,𝑠}(𝑥),𝑞=1,…,𝑚;𝑠=1,…,𝑘−1, with the aid of (2.17) and the Gaussian integration formula in the subinterval 𝐼_𝑠, given by (2.9), we obtain Λ_{𝛼𝑞,𝑠}(𝑥)≃ℎ2^𝑁−1_𝑗=1𝑤_𝑗𝑌^𝑇_𝑠⋅𝐿^(𝑛𝑞)_𝑠𝑥_𝑠𝑗𝑥−𝑥_𝑠𝑗^{𝛼𝑞+1−𝑛𝑞}, (2.29) where the ShLG collocation points 𝑥_𝑠𝑗 on the subinterval 𝐼_𝑠 are defined by (2.10).

From (2.27), (2.16), and (2.17), for 𝑘=1,2,…,𝐾 we have 𝐹𝑥,𝑌^𝑇_𝑘⋅𝐿_𝑘(𝑥),1Γ𝑛₁−𝛼₁^𝑘−1_𝑠=1Λ_𝛼1,𝑠(𝑥)+^𝑥_𝑥𝑘0𝑌^𝑇_𝑘⋅𝐿^(𝑛1)_𝑘(𝑡)(𝑥−𝑡)^{𝛼1+1−𝑛1}1d𝑡,…,Γ𝑛_𝑚−𝛼_𝑚^𝑘−1_𝑠=1Λ_{𝛼𝑚,𝑠}(𝑥)+^𝑥_𝑥𝑘0𝑌^𝑇_𝑘⋅𝐿^(𝑛𝑚)_𝑘(𝑡)(𝑥−𝑡)^{𝛼𝑚+1−𝑛𝑚}d𝑡=0,𝑥∈𝐼_𝑘. (2.30) We now collocate (2.30) at collocation points 𝑥_𝑘𝑗, 𝑘=1,2,…,𝐾 and 𝑗=1,…,𝑁−𝑙 as 𝐹𝑥_𝑘𝑗,𝑌^𝑇_𝑘⋅𝐿_𝑘𝑥_𝑘𝑗,1Γ𝑛₁−𝛼₁^𝑘−1_𝑠=1Λ_𝛼1,𝑠𝑥_𝑘𝑗+^𝑥𝑘𝑗_𝑥𝑘0𝑌^𝑇_𝑘⋅𝐿^(𝑛1)_𝑘(𝑡)𝑥_𝑘𝑗−𝑡^{𝛼1+1−𝑛1}1d𝑡,…,Γ𝑛_𝑚−𝛼_𝑚^𝑘−1_𝑠=1Λ_{𝛼𝑚,𝑠}𝑥_𝑘𝑗+^𝑥𝑘𝑗_𝑥𝑘0𝑌^𝑇_𝑘⋅𝐿^(𝑛𝑚)_𝑘(𝑡)𝑥_𝑘𝑗−𝑡^{𝛼𝑚+1−𝑛𝑚}d𝑡=0. (2.31) 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 𝐹𝑥_𝑘𝑗,𝑦_𝑘𝑗,1Γ𝑛₁−𝛼₁×^𝑘−1_𝑠=1Λ_𝛼1,𝑠𝑥_𝑘𝑗+1𝑛₁−𝛼₁𝑥^{𝑛1−𝛼1}_1𝑗𝑌^𝑇_𝑘⋅𝐿^(𝑛1)_𝑘𝑥_𝑘0+^𝑥𝑘𝑗_𝑥𝑘0𝑌^𝑇_𝑘⋅𝐿^(𝑛1_𝑘⁺¹⁾(𝑡)𝑥_𝑘𝑗−𝑡^{𝛼1−𝑛1},1d𝑡…,Γ𝑛_𝑚−𝛼_𝑚×^𝑘−1_𝑠=1Λ_{𝛼𝑚,𝑠}𝑥_𝑘𝑗+1𝑛_𝑚−𝛼_𝑚𝑥^{𝑛𝑚−𝛼𝑚}_1𝑗𝑌^𝑇_𝑘⋅𝐿^(𝑛𝑚)_𝑘𝑥_𝑘0+^𝑥𝑘𝑗_𝑥𝑘0𝑌^𝑇_𝑘⋅𝐿^(𝑛𝑚_𝑘⁺¹⁾(𝑡)𝑥_𝑘𝑗−𝑡^{𝛼𝑚−𝑛𝑚}d𝑡=0. (2.32) In order to use the Gaussian integration formula in the subinterval 𝐼_𝑘 for (2.32), we transfer the 𝑡-interval [𝑥_𝑘0,𝑥_𝑘𝑗] into the 𝜏-interval 𝐼_𝑘 by means of the transformation 𝜏=(ℎ/𝑥_1𝑗)(𝑡−𝑥_𝑘0)+𝑥_𝑘0. Using this transformation, the Gaussian integration formula and (2.10), we have ^𝑥𝑘𝑗_𝑥𝑘0𝑌^𝑇_𝑘⋅𝐿^(𝑛𝑞_𝑘⁺¹⁾(𝑡)𝑥_𝑘𝑗−𝑡^{𝛼𝑞−𝑛𝑞}𝑥d𝑡=^1𝑗ℎ^𝑥𝑘𝑁_𝑥𝑘0𝑌^𝑇_𝑘⋅𝐿^(𝑛𝑞_𝑘⁺¹⁾𝑥_1𝑗/ℎ𝜏−𝑥_𝑘0+𝑥_𝑘0𝑥_1𝑗−(𝑥_1𝑗/ℎ)𝜏−𝑥_𝑘0^{𝛼𝑞−𝑛𝑞}≃𝑥d𝜏_1𝑗2^𝑁−1_𝑝=1𝑤_𝑝𝑌^𝑇_𝑘⋅𝐿^(𝑛𝑞_𝑘⁺¹⁾𝑥_1𝑗𝑥_1𝑝/ℎ+𝑥_𝑘0𝑥_1𝑗−(𝑥_1𝑗𝑥_1𝑝/ℎ)^{𝛼𝑞−𝑛𝑞}∶=𝐺_𝑞𝑥_1𝑗. (2.33) By (2.33), (2.32) may be approximated as 𝐹𝑥_𝑘𝑗,𝑦_𝑘𝑗,1Γ𝑛₁−𝛼₁^𝑘−1_𝑠=1Λ_𝛼1,𝑠𝑥_𝑘𝑗+1𝑛₁−𝛼₁𝑥^{𝑛1−𝛼1}_1𝑗𝑌^𝑇_𝑘⋅𝐿^(𝑛1)_𝑘𝑥_𝑘0+𝐺₁𝑥_1𝑗1,…,Γ𝑛_𝑚−𝛼_𝑚^𝑘−1_𝑠=1Λ_{𝛼𝑚,𝑠}𝑥_𝑘𝑗+1𝑛_𝑚−𝛼_𝑚𝑥^{𝑛𝑚−𝛼𝑚}_1𝑗𝑌^𝑇_𝑘⋅𝐿^(𝑛𝑚)_𝑘𝑥_𝑘0+𝐺_𝑚𝑥_1𝑗=0. (2.34) In addition, substituting (2.16) and (2.17) into the boundary conditions (1.2) yields 𝐻_𝑟𝑌^𝑇_𝜌0⋅𝐿_𝜌0𝜉₀,…,𝑌^𝑇_𝜌0⋅𝐿_𝜌^(𝑙)₀𝜉₀,…,𝑌^𝑇_𝜌𝑙⋅𝐿_𝜌𝑙𝜉_𝑙,…,𝑌^𝑇_𝜌𝑙⋅𝐿_𝜌^(𝑙)_𝑙𝜉_𝑙=0,𝑟=0,1,…,𝑙, (2.35) where 𝜉_𝑟∈𝐼_𝜌𝑟. Besides, it is required that the approximate solution and its first 𝑙 derivatives be continuous at the interface of subintervals, that is, 𝑌^𝑇_𝑘⋅𝐿_𝑘^(𝑟)𝑥_𝑘𝑁=𝑌^𝑇_𝑘+1⋅𝐿^(𝑟)_𝑘+1𝑥_𝑘+1,0,𝑘=1,2,…,𝐾−1,𝑟=0,1,…,𝑙. (2.36)

Equation (2.34) for 𝑘=1,…,𝐾, 𝑗=1,…,𝑁−𝑙 together with (2.35)-(2.36) gives a system of equations with 𝐾(𝑁+1) set of algebraic equations, which can be solved to find the unknowns of the vectors 𝑌_𝑘,𝑘=1,2,…,𝐾. Consequently, the unknown functions 𝑦(𝑥)|_𝐼𝑘≃𝑦_𝑘(𝑥) 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] 𝐴𝑦^′′(𝑥)+𝐵𝐷^𝛼𝑦(𝑥)+𝐶𝑦(𝑥)=𝑓(𝑥),𝑥∈[0,𝐿], (3.1) with boundary conditions 𝑦(0)=𝑐₀,𝑦(𝐿)=𝑐₁, (3.2) where 𝐴,𝐵,𝐶∈ℝ and 𝐴≠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: 𝑦(𝑥)=^𝑥₀𝐺₃()𝑓(𝑡)𝐺𝑥−𝑡d𝑡,₃(𝑥)=1𝐴^∞_𝑘=0(−1)^𝑘𝐶𝑘!𝐴^𝑘𝑥^2𝑘+1𝐸^(𝑘)_{1/2,2+(3𝑘/2)}−𝐵𝐴√𝑥,𝐸^(𝑘)_𝜆,𝜇(𝑧)=d^𝑘d𝑧^𝑘𝐸_𝜆,𝜇(𝑧)=^∞_𝑗=0𝑗+𝑘!𝑧^𝑗𝑗!Γ𝜆𝑗+𝜆𝑘+𝜇,𝑘=0,1,2,…, (3.3) where 𝐸_𝜆,𝜇 is the Mittag-Leffler function in two parameters and the 𝐺₃ three-term Green’s function. However, in practice, these equations can not be evaluated easily for different functions 𝑓(𝑥). 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, 𝐴=𝐵=𝐶=1, 𝑓(𝑥)=𝑥²+2+4𝑥/𝜋, 𝐿=5, 𝑐₀=0, and 𝑐₁=25. It is readily verified that the exact solution of this case is 𝑦(𝑥)=𝑥². Using the adaptive pseudospectral method in Section 2 with 𝐾=1 and 𝑁=2, the unknowns 𝑦_𝑘𝑖, 𝑘=1,…,𝐾, 𝑖=0,…,𝑁 in (2.16) are found to be 𝑦₁₀=0,𝑦₁₁=6.25,𝑦₁₂=5, (3.4) which lead to the exact solution 𝑦(𝑥)=𝑥². This case was solved in [6] using a collocation-shooting method. Their computed maximum absolute error and 𝐿₂ error norm were 2.00×10⁻¹⁴ and 3.78×10⁻¹², respectively, which show that our method is more efficient.

Case 2. Set 𝛼=3/2, 𝐴=𝐵=𝐶=1, 𝑓(𝑥)=(15/4)√𝑥+(15/8)√𝜋𝑥+𝑥^5/2+1, 𝐿=1, 𝑐₀=1, and 𝑐₁=2. The exact solution of this case, which was considered in [5], is 𝑦(𝑥)=𝑥^5/2+1. In Table 1 the maximum absolute errors for different values of 𝐾 and 𝑁 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 Example 3.1, Case 2.

Case 3. For comparison, the same coefficients as considered in [27] have been used here. Set 𝛼=3/2, 𝐴=1, 𝐵=𝐶=1/2, 𝑓(𝑥)=8, 𝐿=1, 𝑐₀=0, and 𝑐₁=2.95258388. Table 2 shows the comparison of solutions of this case by the present method (with 𝐾=2, 𝑁=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 Example 3.1, Case 3.

Case 4. Set 𝛼=1/2, 𝐴=𝐵=𝐶=1, 𝑓(𝑥)=8, 𝐿=1, 𝑐₀=0, and 𝑐₁=3.10190571. The numerical solutions obtained by the present method (with 𝐾=2, 𝑁=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 Example 3.1, Case 4.

Example 3.2. As a multi-term equation, consider the linear multi-term FBVP described by √𝜋𝑥𝐷^5/2𝑦(𝑥)2+10Γ3³√𝑥𝐷^4/3𝑦(𝑥)3+Γ4⁴√𝑥𝐷^1/4𝑦(𝑥)−40𝑦(𝑥)77=12𝑥+54𝑥²+87𝑥³,𝑦(0)=𝑦^′(0)=0,𝑦(1)=1. (3.5) The exact solution to this problem is 𝑦(𝑥)=𝑥³. 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 𝐾 and 𝑁 with the errors reported in [13] using operational matrix of fractional derivatives using B-spline functions. Note that in [13], for each value of 𝐽, the obtained algebraic system is of order 2^𝐽+1, while in the present method the obtained algebraic system is of order 𝐾(𝑁+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 Example 3.2.

Example 3.3. Consider the nonlinear multi-term FBVP described by 𝑦^′′(𝑥)4+Γ5⁵√𝑥⁶𝐷^6/5𝑦(𝑥)+119Γ56⁶√𝑥𝐷^1/6𝑦(𝑥)−𝑦^′(𝑥)²1=2+𝑥10²,𝑦(0)=1,𝑦(1)=2. (3.6) The exact solution to this problem is 𝑦(𝑥)=𝑥²+1. In Table 5, we compare the maximum absolute errors obtained using the present adaptive method for different values of 𝐾 and 𝑁 with the errors reported in [13] using operational matrix of fractional derivatives using B-spline functions.

Table 5 

Comparison of maximum absolute errors for Example 3.3.

Example 3.4. Consider the following nonlinear FBVP: 𝐷^𝛼𝑦(𝑥)+𝑎𝑦^𝑛(𝑥)(𝑥)𝑦(0)=𝑓,1<𝛼⩽2,=𝑐₀,𝑦(1)=𝑐₁. (3.7) For comparison, we choose 𝑛=2, 𝑎=𝑒^−2𝜋, 𝑓(𝑥)=(105√𝜋/16Γ((9/2)−𝛼))𝑥^(7/2)−𝛼+𝑒^−2𝜋𝑥⁷, 𝑐₀=0, and 𝑐₁=1. It is readily verified that the exact solution is 𝑦(𝑥)=𝑥^7/2. In Table 6, the absolute errors obtained using the present adaptive pseudospectral method for 𝐾=4 and 𝑁=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 Example 3.4.

Example 3.5. In this example, to show the applicability of the present method for larger interval, we consider the nonlinear FBVP described by 𝐷^𝛼𝑦(𝑥)+𝑦²(𝑥)=𝐸²_𝛼−𝑥^𝛼−𝐸_𝛼−𝑥^𝛼𝑦(0)(),1<𝛼⩽2,=1,𝑦10=𝐸_𝛼−10^𝛼. (3.8) The exact solution of this problem is given by 𝑦(𝑥)=𝐸_𝛼(−𝑥^𝛼), where 𝐸_𝛼(𝑧)=∑^∞_𝑘=0(𝑧^𝑘/Γ(𝛼𝑘+1)) is the Mittag-Leffler function.
In Table 7, the maximum absolute errors in the interval [0,10] for 𝛼=1.75 and different values of 𝐾 and 𝑁 are presented, which shows the efficiency of the present method for FBVPs in large domains. Also, the numerical results for 𝑦(𝑥) by adaptive pseudospectral method for 𝐾=20, 𝑁=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 𝑦(𝑥)=cos(𝑥). 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.

Figure 1 

Comparison of 𝑦(𝑥) for 𝐾=20, 𝑁=30 with exact solutions for Example 3.5.

Example 3.6. Finally consider the nonlinear multi-term FBVP described by 𝑎𝑦^′′(𝑥)+𝑏𝐷^𝛼2𝑦(𝑥)+𝑐𝐷^𝛼1𝑦(𝑥)²+𝑒𝑦³(𝑥)(𝑥)𝑦(0)(2)=8=𝑓,0⩽𝑥⩽2,=0,𝑦3, (3.9) where 𝑎,𝑏,𝑐,𝑒∈ℝ, 0<𝛼₁⩽1, 1<𝛼₂⩽2 and 𝑓(𝑥)=2𝑎𝑥+(2𝑏/Γ(4−𝛼₂))𝑥^3−𝛼2+𝑐((2/Γ(4−𝛼₁))𝑥^3−𝛼1)²+𝑒((1/3)𝑥³)³. The exact solution to this problem is 𝑦(𝑥)=(1/3)𝑥³.
For 𝑎=𝑏=𝑐=𝑒=1, 𝛼₁=0.555, and 𝛼₂=1.455 the maximum absolute errors obtained using the adaptive pseudospectral method are given in Table 8. Also, for 𝑎=0.1, 𝑏=𝑐=𝑒=0.5, 𝛼₁=0.219, and 𝛼₂=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 for 𝑎=𝑏=𝑐=𝑒=1,𝛼₁=0.555,𝛼₂=1.455 for Example 3.6.

Table 9 

Maximum absolute errors for 𝑎=0.1,𝑏=𝑐=𝑒=0.5,𝛼₁=0.219,𝛼₂=1.965 for Example 3.6.

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 first 𝑙 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.

Appendix

The Derivation of Bagley-Torvik Equation

Here, in order to give a physical understanding of fractional differential equations, the derivation of Bagley-Torvik equation, which describes the modeling of the motion of a rigid plate in a Newtonian fluid, is given.

Consider a half-space Newtonian viscous fluid in which certain motions are induced by the general transverse motion of an infinite plate. The equation of motion of the half-space fluid is the diffusion equation: 𝜌𝜕𝑣(𝑧,𝑡)𝜕𝜕𝑡=𝜇²𝑣(𝑧,𝑡)𝜕𝑧², (A.1) where 𝜌 is the fluid density, 𝜇 is the viscosity, and 𝑣(𝑧,𝑡) describes the transverse fluid velocity as a function of 𝑧 and 𝑡. Taking the Laplace transform of (A.1) and using the properties of the Laplace transform, one obtains [𝑣(()𝜕𝜌𝑠𝐿𝑧,𝑡)]−𝜌𝑣𝑧,𝑡=0=𝜇²𝜕𝑧²𝐿[𝑣(𝑧,𝑡)]. (A.2) Torvik and Bagley assumed the initial velocity profile in the fluid to be zero and thus (A.2) reduces to [𝑣(𝜕𝜌𝑠𝐿𝑧,𝑡)]=𝜇²𝜕𝑧²𝐿[𝑣(𝑧,𝑡)]. (A.3) Since the Laplace transformation is evaluated with respect to the time variable, only the following representation for the velocity profile with respect to the depth 𝑧 can be used: 𝑣(𝑧,𝑡)=𝑣(𝑡)𝑒^𝜆𝑧, (A.4) thus 𝐿[𝑣(𝑧,𝑡)]=𝑒^𝜆𝑧𝐿[𝑣(𝑡,𝜕)]₂𝜕𝑧²𝐿[𝑣(𝑧,𝑡)]=𝜆²𝑒^𝜆𝑧𝐿[𝑣(𝑡)]. (A.5) With insertion of (A.5) in (A.3) the following algebraic equation for the unknown parameter 𝜆 is obtained: 𝜆=𝑠𝜌𝜇. (A.6)

Next, the shear stress relationship of the Newtonian fluid given as 𝜎()()𝑧,𝑡=𝜇𝜕𝑣𝑧,𝑡𝜕𝑧 (A.7) can be transformed into the Laplace domain using the above results: 𝐿[𝜎(𝑧,𝑡)]=𝜇𝑠𝜌𝜇𝑒^{(𝑠𝜌/𝜇)𝑧}𝐿[𝑣(𝑡)]=√𝜇𝜌√𝑠𝐿[𝑣(𝑧,𝑡)]. (A.8) Equation (A.8) can be restated as 𝐿[𝜎(𝑧,𝑡)]=√𝑠𝜇𝜌√𝑠𝐿[𝑣(𝑧,𝑡)]. (A.9) Now, the following two transforms can be identified in (A.9): [𝑣(()𝑠𝐿𝑧,𝑡)]=𝐿𝜕𝑣𝑧,𝑡,1𝜕𝑡√𝑠1=𝐿Γ(1/2)√𝑡. (A.10) With substituting (A.10) into (A.9), one obtains 𝐿[𝜎(𝑧,𝑡)]=√1𝜇𝜌𝐿Γ(1/2)√𝑡⋅𝐿[𝑣̇(𝑧,𝑡)]. (A.11) The product of two transforms in (A.11) corresponds to the following convolution when evaluating the inverse transformation: 𝜎(𝑧,𝑡)=√1𝜇𝜌Γ(1/2)^𝑡₀𝑣̇(𝑧,𝜏)(𝑡−𝜏)^1/2d𝜏=√𝜇𝜌𝐷_𝑡^1/2𝑣(𝑧,𝑡), (A.12) which introduces a fractional derivative of degree 𝛼=1/2 within the shear stress-velocity relationship of a half-space Newtonian fluid.

Finally, consider a rigid plate of mass 𝑚 immersed into an infinite Newtonian fluid. The plate is held at a fixed point by means of a spring of stiffness 𝑘. It is assumed that the motions of the spring do not influence the motion of the fluid, and that the surface 𝐴 of the plate is very large, such that the stress-velocity relationship in (A.12) is valid on both sides of the plate. Equilibrium of all forces acting on the plate gives 𝑚𝑦^′′(𝑡)+𝑘𝑦(𝑡)+2𝐴𝜎(𝑧=0,𝑡)=0. (A.13) By substituting (A.12) one obtains 𝑚𝑦^′′(𝑡)+𝑘𝑦(𝑡)+2𝐴√𝜇𝜌𝐷_𝑡^1/2𝑣(𝑧=0,𝑡)=0. (A.14) With 𝑣(𝑧=0,𝑡)=𝑦^′(𝑡), a fractional differential equation of degree 𝛼=3/2 follows for the displacement of a rigid plate immersed into an infinite Newtonian fluid, as follows: 𝑚𝑦^′′(𝑡)+𝑘𝑦(𝑡)+2𝐴√𝜇𝜌𝐷_𝑡^3/2𝑦(𝑡)=0. (A.15)

Acknowledgment

The authors are grateful to the reviewers for the constructive comments which led to an improvement in the paper.