#### 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 [13].

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 [912], operational matrix methods [1316], 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 subject to where can be nonlinear in general, , , , are linear functions, the points , lie in , and denotes the Caputo-fractional derivative of order , defined as follows [23]: 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, , , are the eigenfunctions of the singular Sturm-Liouville problem Also, they are orthogonal with respect to inner product on the interval with the weight function , that is, where is the Kronecker delta. The Legendre polynomials satisfy the recursion relation where and . If is normalized so that , then, for any , the Legendre polynomials in terms of power of are where denotes the integer part of .

The Legendre-Gauss (LG) collocation points are the roots of . No explicit formulas are known for the LG points; however, they are computed numerically using existing subroutines. The LG points have the property that is exact for polynomials of degree at most , where are LG quadrature weights. For more details about Legendre polynomials, see [24].

The shifted Legendre polynomials on the interval are defined by 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 . Thus, any function can be expanded in terms of shifted Legendre polynomials.

The ShLG collocation points on the interval are obtained by shifting the LG points, , using the transformation Thanks to the property of the standard LG quadrature, it follows that for any polynomial of degree at most on , where , are ShLG quadrature weights.

##### 2.2. Function Approximation

Suppose that the interval is divided into subintervals , , where . Let be the solution of the problem in (1.1)-(1.2) in the subinterval . Consider now the ShLG collocation points on the th subinterval , , obtained using (2.8). Obviously, Also, consider two additional noncollocated points and . Let us define where 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 or equivalently where .

Here, it can be easily seen that for and , we have Thus, by utilizing (2.15) for (2.14), the approximation of within each subinterval can be restated as where and are matrices given by and . It is important to observe that the series (2.16) includes the Lagrange polynomials associated with the noncollocated points 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 at the first subinterval.

By times ( is defined in (1.3)) differentiating of (2.16), we obtain where .

##### 2.3. Convergence Rate

For we introduce the piecewise polynomials space which is the space of the continuous functions over whose restrictions on each subinterval are polynomials of degree . Then, for any continuous function in , 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 with . Then and, for , if , then and if , then

Note that denotes a positive constant that depends on , but which is independent of the function and integer . Moreover, we introduce the seminorm of , , , , as

Remark 2.2. Whenever , using (2.19)–(2.22), we get and, for , if , then and if , then Equations (2.23)–(2.25) show that if is infinitely smooth on and , the convergence rate of to is faster than to the power of and any power of , which is superior to that for the global collocation method over . 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 The problem is to find , , satisfying (2.26) and (1.2).

The generally nonlinear SVIDE in (2.26) is given in subinterval , as follows: where and It is important to note that, at the first subinterval, the summations in (2.27) are automatically discarded. For approximating the functions , with the aid of (2.17) and the Gaussian integration formula in the subinterval , given by (2.9), we obtain where the ShLG collocation points on the subinterval are defined by (2.10).

From (2.27), (2.16), and (2.17), for we have We now collocate (2.30) at collocation points , and as 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 In order to use the Gaussian integration formula in the subinterval for (2.32), we transfer the -interval into the -interval by means of the transformation . Using this transformation, the Gaussian integration formula and (2.10), we have By (2.33), (2.32) may be approximated as In addition, substituting (2.16) and (2.17) into the boundary conditions (1.2) yields where . Besides, it is required that the approximate solution and its first derivatives be continuous at the interface of subintervals, that is,

Equation (2.34) for , together with (2.35)-(2.36) gives a system of equations with set of algebraic equations, which can be solved to find the unknowns of the vectors . 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] with boundary conditions where and . Bagley-Torvik equation involving fractional derivative of order or 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 gave the analytical solution with homogeneous initial conditions by using Green’s function, as follows: 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 , , , , , and . It is readily verified that the exact solution of this case is . Using the adaptive pseudospectral method in Section 2 with and , the unknowns , , in (2.16) are found to be 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 and , respectively, which show that our method is more efficient.

Case 2. Set , , , , , and . The exact solution of this case, which was considered in [5], is . 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.

Case 3. For comparison, the same coefficients as considered in [27] have been used here. Set , , , , , , and . Table 2 shows the comparison of solutions of this case by the present method (with , ), GTC method [27] and the analytical solution [2], and the good agreement of our adaptive pseudospectral solution with analytical solution.

Case 4. Set , , , , , and . The numerical solutions obtained by the present method (with , ), 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.

Example 3.2. As a multi-term equation, consider the linear multi-term FBVP described by 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 , while in the present method the obtained algebraic system is of order . 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.

Example 3.3. Consider the nonlinear multi-term FBVP described by The exact solution to this problem is . 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.

Example 3.4. Consider the following nonlinear FBVP: For comparison, we choose , , , , and . It is readily verified that the exact solution is . In Table 6, the absolute errors obtained using the present adaptive pseudospectral method for and 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 , where is the Mittag-Leffler function.
In Table 7, the maximum absolute errors in the interval for 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 , and , 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 , the exact solution is given as . 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.

Example 3.6. Finally consider the nonlinear multi-term FBVP described by where , , and . The exact solution to this problem is .
For , , and the maximum absolute errors obtained using the adaptive pseudospectral method are given in Table 8. Also, for , , , and 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.

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

#### 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: 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 Torvik and Bagley assumed the initial velocity profile in the fluid to be zero and thus (A.2) reduces to 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: thus With insertion of (A.5) in (A.3) the following algebraic equation for the unknown parameter is obtained:

Next, the shear stress relationship of the Newtonian fluid given as can be transformed into the Laplace domain using the above results: Equation (A.8) can be restated as Now, the following two transforms can be identified in (A.9): With substituting (A.10) into (A.9), one obtains 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 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 By substituting (A.12) one obtains With , a fractional differential equation of degree follows for the displacement of a rigid plate immersed into an infinite Newtonian fluid, as follows:

#### Acknowledgment

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