An Adaptive Pseudospectral Method for Fractional Order Boundary Value Problems

An adaptive pseudospectral method is presented for solving a class of multiterm fractional boundary value problems (cid:2) FBVP (cid:3) which involve Caputo-type fractional derivatives. The multiterm FBVP is ﬁrst converted into a singular Volterra integrodi ﬀ erential equation (cid:2) SVIDE (cid:3) . By dividing the interval of the problem to subintervals, the unknown function is approximated using a piecewise interpolation polynomial with unknown coe ﬃ cients which is based on shifted Legendre-Gauss (cid:2) ShLG (cid:3) 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 e ﬃ cient and accurate. Abstract the solutions of initial value and boundary value problems for linear and nonlinear fractional di ﬀ erential equations. These methods include ﬁnite di ﬀ erence approximation method (cid:6) 4 (cid:7) , collocation method (cid:6) 5, 6 (cid:7) , the Adomian decomposition method (cid:6) 7, 8 (cid:7) , variational iteration method (cid:6) 9–12 (cid:7) , operational matrix methods (cid:6) 13–16 (cid:7) , and homotopy methods (cid:6) 17, 18 (cid:7) . In (cid:6) 19 (cid:7) suitable spline functions of polynomial form are derived and used to solve linear and nonlinear fractional di ﬀ erential equations. The authors of (cid:6) 20 (cid:7) have investigated the existence and multiplicity of positive solutions of a nonlinear fractional di ﬀ erential equation initial value problem. Furthermore, some physical and geometrical interpretations of fractional operators and fractional di ﬀ erential equations have been of concern to many authors (cid:6) 12, 21, 22 (cid:7) e cient adaptive pseudospectral for (cid:2) of the form


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 where F can be nonlinear in general, 0 < α 1 < α 2 < · · · < α m , l < α m ≤ l 1, L ∈ R, H r 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 : D α q y x 1 Γ n q − α q x 0 y n q t x − t α q 1−n q dt, n q α q 1, q 1, 2, . . . , m, 1.3 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 first l 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 3 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.

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 .

Review of Legendre and Shifted Legendre Polynomials
The Legendre polynomials, P i z , i 0, 1, 2, . . ., are the eigenfunctions of the singular Sturm-Liouville problem Also, they are orthogonal with respect to L 2 inner product on the interval −1, 1 with the weight function w z 1, that is, where δ ij is the Kronecker delta. The Legendre polynomials satisfy the recursion relation where P 0 z 1 and P 1 z z. If P i z is normalized so that P i 1 1, then, for any i, the Legendre polynomials in terms of power of z are where i/2 denotes the integer part of i/2.

Abstract and Applied Analysis
The Legendre-Gauss LG collocation points −1 < z 1 < z 2 < · · · < z N−1 < 1 are the roots of P N−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 The shifted Legendre polynomials on the interval x ∈ a, b are defined by which are obtained by an affine transformation from the Legendre polynomials. The set of shifted Legendre polynomials is a complete L 2 a, b -orthogonal system with the weight function w x 1. Thus, any function f ∈ L 2 a, b can be expanded in terms of shifted Legendre polynomials.
The ShLG collocation points a < x 1 < x 2 < · · · < x N−1 < b on the interval a, b are obtained by shifting the LG points, z j , using the transformation Thanks to the property of the standard LG quadrature, it follows that for any polynomial p of degree at most 2N − 3 on where w j b − a /2 w j , 1 j N − 1 are ShLG quadrature weights.

Function Approximation
Suppose that the interval 0, L is divided into K subintervals I k k − 1 h, kh , k 1, 2, . . . , K, where h L/K. Let y k x be the solution of the problem in 1.1 -1.2 in the subinterval I k . Consider now the ShLG collocation points k − 1 h < x k1 < · · · < x k,N−1 < kh on the kth subinterval I k , k 1, 2, . . . , K, obtained using 2.8 . Obviously, Also, consider two additional noncollocated points x k0 k − 1 h and x kN kh. Let us define is a basis of Lagrange interpolating polynomials on the subinterval I k that satisfy L ki x kj δ ij , where δ ij is the Kronecker delta function. The L 2 I k -orthogonal projection I N : kN is a mapping in a way that for any y k ∈ L 2 I k or equivalently Here, it can be easily seen that for i 0, 1, . . . , N and k 1, 2, . . . , K, we have Thus, by utilizing 2.15 for 2.14 , the approximation of y k x within each subinterval I k can be restated as where Y k and L k x are N 1 × 1 matrices given by Y k y k0 , . . . , y kN T and L k x It is important to observe that the series 2.16 includes the Lagrange polynomials associated with the noncollocated points x k0 k − 1 h and x kN 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 L 1i x at the first subinterval.

Abstract and Applied Analysis
By n m 1 times n m is defined in 1.3 differentiating of 2.16 , we obtain

Convergence Rate
For N 1 we introduce the piecewise polynomials space which is the space of the continuous functions over 0, L whose restrictions on each subinterval I k are polynomials of degree N. Then, for any continuous function y in 0, L , the piecewise interpolation polynomial ψ N y coincides on each subinterval I k with the interpolating polynomial I N y of y k y| I k 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.
and, for 1 u v, if h 1, then

2.20
and if h > 1, then Note that c 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 Abstract and Applied Analysis 7 and, for u 1, if h 1, then and if h > 1, then

2.25
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.

Problem Replacement and the Solution Technique
Consider the multi-term FBVP in 1.1 -1.2 . With substituting the definition of the Caputoderivative 1.3 into 1.1 , we can convert 1.1 into an equivalent SVIDE as

2.26
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 I k , k 1, 2, . . . , K as follows:

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

2.32
Abstract and Applied Analysis 9 In order to use the Gaussian integration formula in the subinterval I k for 2.32 , we transfer the t-interval x k0 , x kj into the τ-interval I k by means of the transformation τ h/x 1j t − x k0 x k0 . Using this transformation, the Gaussian integration formula and 2.10 , we have : G q x 1j .

2.33
By 2.33 , 2.32 may be approximated as

2.34
In addition, substituting 2. 16  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,

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.  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. which lead to the exact solution y x x 2 . This case was solved in 6 using a collocationshooting method. Their computed maximum absolute error and L 2 error norm were 2.00 × 10 −14 and 3.78 × 10 −12 , respectively, which show that our method is more efficient.     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.
The exact solution to this problem is y x x 3 . 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 2 J 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.

3.6
The exact solution to this problem is y x x 2 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.

3.7
For comparison, we choose n 2, a e −2π , f x 105 √ π/16Γ 9/2 − α x 7/2 −α e −2π x 7 , c 0 0, and c 1 1. It is readily verified that the exact solution is y x x 7/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.
Example 3.5. In this example, to show the applicability of the present method for larger interval, we consider the nonlinear FBVP described by The exact solution of this problem is given by y x E α −x α , where E α z ∞ k 0 z k /Γ αk 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 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 integerorder differential equations. where a, b, c, e ∈ R, 0 < α 1 1, 1 < α 2 2 and f x 2ax 2b/Γ 4 − α 2 x 3−α 2 c 2/Γ 4 − α 1 x 3−α 1 2 e 1/3 x 3 3 . The exact solution to this problem is y x For a 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.

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 l 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.
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 z can be used: Now, the following two transforms can be identified in A.9 : sL v z, t L ∂v z, t ∂t , A.10 With substituting A.10 into A.9 , one obtains L σ z, t μρL 1 Γ 1/2 √ t · L v z, t .
A.11  The product of two transforms in A.11 corresponds to the following convolution when evaluating the inverse transformation: 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 m immersed into an infinite Newtonian fluid. The plate is held at a fixed point by means of a spring of stiffness k. It is assumed that the motions of the spring do not influence the motion of the fluid, and that the surface A 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 With v z 0, t y t , a fractional differential equation of degree α 3/2 follows for the displacement of a rigid plate immersed into an infinite Newtonian fluid, as follows: my t ky t 2A μρD 3/2 t y t 0. A.15