Abstract and Applied Analysis

Volume 2012 (2012), Article ID 381708, 19 pages

doi:10.1155/2012/381708

Research Article

An Adaptive Pseudospectral Method for Fractional Order Boundary Value Problems

Mohammad Maleki,¹ Ishak Hashim,¹ Majid Tavassoli Kajani,² and Saeid Abbasbandy³

¹School of Mathematical Sciences, Universiti Kebangsaan Malaysia, 43600 Bangi Selangor, Malaysia

²Department of Mathematics, Islamic Azad University, Khorasgan Branch, Isfahan, Iran

³Department of Mathematics, Imam Khomeini International University, Ghazvin 34149-16818, Iran

Received 3 June 2012; Revised 17 August 2012; Accepted 30 August 2012

Academic Editor: Dumitru Bǎleanu

Copyright © 2012 Mohammad Maleki et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

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.

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

Table 2: Comparison of solutions for Example 3.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.95258388

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.

Table 3: Comparison of solutions for Example 3.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.10190571

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.

Table 4: Comparison of maximum absolute errors for Example 3.2.

Present method
Method [13]

Error

4

5

6

7

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.

Table 5: Comparison of maximum absolute errors for Example 3.3.

Present method
Method [13]

Error

4

5

6

7

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.

Table 6: Comparison of absolute errors for Example 3.4.

α

1.1
Present method

Method [16]

1.3
Present method

Method [16]

1.5
Present method

Method [16]

1.7
Present method

Method [16]

1.9
Present method

Method [16]

2.0
Present method

Method [16]

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.

Table 7: Maximum absolute errors for for Example 3.5.

Figure 1: Comparison of for , with exact solutions for Example 3.5.

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.

Table 8: Maximum absolute errors for for Example 3.6.

Table 9: Maximum absolute errors for 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:

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.

References

D. Baleanu, K. Diethelm, E. Scalas, and J. J. Trujillo, Fractional Calculus: Models and Numerical Methods, vol. 3 of Series on Complexity, Nonlinearity and Chaos, World Scientific Publishing, Hackensack, NJ, USA, 2012.

I. Podlubny, Fractional Differential Equations, vol. 198, Academic Press, San Diego, Calif, USA, 1999.

S. Das, Functional Fractional Calculus for System Identification and Controls, Springer, New York, NY, USA, 2008.

K. Moaddy, I. Hashim, and S. Momani, “Non-standard finite difference schemes for solving fractional-order Rössler chaotic and hyperchaotic systems,” Computers & Mathematics with Applications, vol. 62, no. 3, pp. 1068–1074, 2011.

A. Pedas and E. Tamme, “Piecewise polynomial collocation for linear boundary value problems of fractional differential equations,” Journal of Computational and Applied Mathematics, vol. 236, no. 13, pp. 3349–3359, 2012.

Q. M. Al-Mdallal, M. I. Syam, and M. N. Anwar, “A collocation-shooting method for solving fractional boundary value problems,” Communications in Nonlinear Science and Numerical Simulation, vol. 15, no. 12, pp. 3814–3822, 2010.

N. J. Ford and J. A. Connolly, “Systems-based decomposition schemes for the approximate solution of multi-term fractional differential equations,” Journal of Computational and Applied Mathematics, vol. 229, no. 2, pp. 382–391, 2009.

S. Momani and Z. Odibat, “Analytical solution of a time-fractional Navier-Stokes equation by Adomian decomposition method,” Applied Mathematics and Computation, vol. 177, no. 2, pp. 488–494, 2006.

G. Wu and E. W. M. Lee, “Fractional variational iteration method and its application,” Physics Letters A, vol. 374, no. 25, pp. 2506–2509, 2010.

S. Abbasbandy, “An approximation solution of a nonlinear equation with Riemann-Liouville's fractional derivatives by He's variational iteration method,” Journal of Computational and Applied Mathematics, vol. 207, no. 1, pp. 53–58, 2007.

J.-H. He, “A short remark on fractional variational iteration method,” Physics Letters A, vol. 375, no. 38, pp. 3362–3364, 2011.

J. H. He, G. C. Wu, and F. Austin, “The variational iteration method which should be followed,” Nonlinear Science Letters A, vol. 1, no. 1, pp. 1–30, 2010.

M. Lakestani, M. Dehghan, and S. Irandoust-pakchin, “The construction of operational matrix of fractional derivatives using B-spline functions,” Communications in Nonlinear Science and Numerical Simulation, vol. 17, no. 3, pp. 1149–1162, 2012.

M. U. Rehman and R. A. Khan, “A numerical method for solving boundary value problems for fractional differential equations,” Applied Mathematical Modelling, vol. 36, no. 3, pp. 894–907, 2012.

E. H. Doha, A. H. Bhrawy, and S. S. Ezz-Eldien, “A new Jacobi operational matrix: an application for solving fractional differential equations,” Applied Mathematical Modelling, vol. 36, no. 10, pp. 4931–4943, 2012.

M. U. Rehman and R. Ali Khan, “The Legendre wavelet method for solving fractional differential equations,” Communications in Nonlinear Science and Numerical Simulation, vol. 16, no. 11, pp. 4163–4173, 2011.

O. Abdulaziz, I. Hashim, and S. Momani, “Application of homotopy-perturbation method to fractional IVPs,” Journal of Computational and Applied Mathematics, vol. 216, no. 2, pp. 574–584, 2008.

I. Hashim, O. Abdulaziz, and S. Momani, “Homotopy analysis method for fractional IVPs,” Communications in Nonlinear Science and Numerical Simulation, vol. 14, no. 3, pp. 674–684, 2009.

A. Al-Rabtah, S. Momani, and M. A. Ramadan, “Solving linear and nonlinear fractional differential equations using spline functions,” Abstract and Applied Analysis, vol. 2012, Article ID 426514, 9 pages, 2012.

D. Baleanu, H. Mohammadi, and S. Rezapour, “Positive solutions of an initial value problem for nonlinear fractional differential equations,” Abstract and Applied Analysis, vol. 2012, Article ID 837437, 7 pages, 2012.

J.-H. He, S. K. Elagan, and Z. B. Li, “Geometrical explanation of the fractional complex transform and derivative chain rule for fractional calculus,” Physics Letters A, vol. 376, no. 4, pp. 257–259, 2012.

M. Moshrefi-Torbati and J. K. Hammond, “Physical and geometrical interpretation of fractional operators,” Journal of the Franklin Institute B, vol. 335, no. 6, pp. 1077–1086, 1998.

M. Caputo, “Linear models of dissipation whose Q is almost frequency independent,” Journal of the Royal Australian Society, vol. 13, Part 2, pp. 529–539, 1967.

C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral Methods: Fundamentals in Single Domains, Scientific Computation, Springer, Berlin, Germany, 2006.

M. Maleki, I. Hashim, and S. Abbasbandy, “Solution of time-varying delay systems using an adaptive collocation method,” Applied Mathematics and Computation, vol. 219, no. 4, pp. 1434–1448, 2012.

P. J. Torvik and R. L. Bagley, “On the appearance of the fractional derivative in the behavior of real materials,” Journal of Applied Mechanics, vol. 51, no. 2, pp. 294–298, 1984.

Y. Çenesiz, Y. Keskin, and A. Kurnaz, “The solution of the Bagley-Torvik equation with the generalized Taylor collocation method,” Journal of the Franklin Institute, vol. 347, no. 2, pp. 452–466, 2010.

S. Momani and Z. Odibat, “Numerical comparison of methods for solving linear differential equations of fractional order,” Chaos, Solitons & Fractals, vol. 31, no. 5, pp. 1248–1255, 2007.