Journal of Applied Mathematics

Volume 2013 (2013), Article ID 864025, 20 pages

http://dx.doi.org/10.1155/2013/864025

Research Article

Cubic Spline Collocation Method for Fractional Differential Equations

Shui-Ping Yang¹ and Ai-Guo Xiao²

¹Department of Mathematics, Huizhou University, Guangdong 516007, China

²School of Mathematics and Computational Science, Hunan Key Laboratory for Computation and Simulation in Science and Engineering, Xiangtan University, Hunan 411105, China

Received 11 February 2013; Revised 17 April 2013; Accepted 7 May 2013

Academic Editor: Alain Miranville

Copyright © 2013 Shui-Ping Yang and Ai-Guo Xiao. 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.

We discuss the cubic spline collocation method with two parameters for solving the initial value problems (IVPs) of fractional differential equations (FDEs). Some results of the local truncation error, the convergence, and the stability of this method for IVPs of FDEs are obtained. Some numerical examples verify our theoretical results.

1. Introduction

In the recent past years, the use of fractional differential equations (FDEs) has gained considerable popularity in some fields (e.g., nonlinear oscillation of earthquake [1], fluid-dynamic traffic model [2], material viscoelastic theory and physics [3–6], etc.). Based on these requirements, the numerical approach has become very important to solve FDEs and analyze the experimental data which are described in a fractional way.

The numerical approaches to FDEs have been recently studied by numerous authors [6–15]. However, the state of art is far less advanced for general fractional order differential equations. In the recent years, spline collocation methods including wavelet methods have been successfully applied to some initial value problems (IVPs) and initial-boundary value problems of FDEs [16–23]. Pedas and Tamme [16, 18] discussed spline collocation methods for some classes of IVPs of linear multiterm fractional differential equations and obtained the corresponding convergence results. Li et al. [19] used the higher-order piecewise interpolation for the fractional integral and fractional derivatives, proposed a higher-order algorithm based on the Simpson method for FDEs, and gave the error and stability results. Li [21] applied the cubic B-spline wavelet collocation method to FDEs by the approach that such problems are converted into a system of algebraic equations which is suitable for computer programming.

It is worth notice that the spline collocation methods for the IVPs of FDEs are often achieved by their solving of the equivalent IVPs of integral equations with weakly singular kernels. In this paper, the cubic spline collocation method is designed to solve directly the IVPs of general FDEs. And the corresponding theoretical results of the local truncation error, the convergence, and the stability of the cubic spline collocation method for the IVPs of general FDEs are given.

This paper is organized as follows. In Section 2, we propose the spline cubic collocation method for solving the IVPs of FDEs. In Section 3, the corresponding theoretical results of the convergence and the stability are also given. In Section 4, the theoretical results are identified by some numerical examples.

In the following text, we first recall the basic definitions of fractional calculus [6]. Usually, denotes the Caputo fractional derivative of order as

where is the classical differential operator of order , is times continuously differentiable, and denotes the integral operator of order as

where is the identity operator.

As is well known, there are some different definitions of fractional operator except the Caputo fractional derivative. From a theoretical point of view, the most natural approach is the Riemann-Liouville definition as

The relationship between the Caputo definition and the Riemann-Liouville definition can be given by the following Lemma.

Lemma 1 (see [5, 6]). Let and . Assume that is times continuously differentiable, and exist; then

2. Cubic Spline Collocation Method for FDEs

Consider the initial value problem

where is a given continuous mapping and satisfies the Lipschitz condition

It is well known that the IVP (5) is equivalent to the problem [24]

Let be the grid points of the uniform partition of into the subintervals . On each subinterval , the cubic -spline function can be represented as

where

and , .

Let , . We can write (9) as

where . Obviously, , and , for all .

Remark 2. In this algorithm, the initial value must be provided, but it does not need to be given in problem. That is, if we approximate , then the approximation of is given by using the low-order methods. In this paper, we obtained approximation of by using one-order BDF method. That is,

Then

where is a tiny stepsize.

By using the collocation conditions in each subinterval and (5), we have

where are the collocation points, and , is a constant. And according to (10),   we have

where , .

Using the definition of the Caputo fractional derivative in the form (1) and using (13),  we have

Through calculation, we can get

where

In order to obtain the corresponding iteration formula, we denote

Obviously, it follows from the inequalities that the matrix is invertible, and

where

Let

Then, it follows that

For (22), we can obtain their numerical solutions by using the Newton iterative method.

In addition, applying the collocation conditions in each subinterval and using (7), we have

where are the collocation points, , and

where ,,,, .

In this paper, for convenience, the cubic spline collocation method based on the direct discretization with (13) is called direct spline collocation method (DSCM for short); and the cubic spline collocation method based on the fractional order integral equations (7) and (27) is called indirect spline collocation method (ISCM for short).

3. Theoretical Results

Lemma 3 (see [25]). Suppose that , and satisfy

where is independent of . Then

where is a Mittag-Leffler function.

Lemma 4. Assume that the function . Then, for all , one has

where , is a constant and satisfies or , , is infinitesimal of higher order.

Proof. Using the Taylor formula, we have

Applying the integral mean value theorem, we can obtain

where . Since , is bounded on . We have

Consider the following.

(i) If , then

(ii) If , then

Thus,

By means of (30), Lemma 4 is established obviously.

Theorem 5 (the local truncation error). If the analytical solution of the problem (5) , then the local truncation error of DSCM is .

Proof. Applying the Taylor expansion and the definition of , we have

From the definition of the Caputo fractional derivative, we obtain

That is,

where . Moreover,

Supposing , for all , and substituting into (13), we have

where is the local truncation error.

Combining (38), (37) with the Lipschitz condition yields

where is an appropriate positive constant. And the norm is the 1-norm of the matrix in this paper.

To prove the convergence of DSCM for the problem (5), we first give the following lemmas.

Lemma 6. If the analytical solution of the problem (5) and the matrixes and in (14) satisfy that exists, and

then the numerical solutions of ISCM for problem (5) satisfy that .

Proof. For (23), we have

It follows from (7) that the analytical solution of the problem (5) satisfies

Moreover, (44) and (45) yield

By using the Lipschitz condition, we have

It follows from the definition of that

where , . Obviously, for (14) and (46), the matrixes and are continuous on ; thus, we have

Therefore, there exist , such that

Obviously,

Applying the integral mean theorem, we have

It follows from the inequalities and that

Combining (45), (48), and (49) with (51) and defining

where , , we can obtain

where are some positive constants. If the matrix is invertible, then it follows from (47) that we have

Hence,

where is an appropriate positive constant. Take , . Then , , and

Since , we have

where , is an appropriate positive constant. According to (53), we obtain

where is an appropriate constant. Obviously, there exist and such that for . Thus,

Let ; then

where is an appropriate constant. Applying Lemma 3, we get

By using the convergence of the Mittag-Leffler function [6], we have

By means of (24), (47), we obtain

where is a positive constant. Consequently, we have

Remark 7. If , then

is invertible, and

Consider

When , the range of values can be obtained by using Algorithm 8, which can be shown in Figures 1 and 2, respectively.

(a)

(b)

Figure 1: The range of values when (a) or (b).

Figure 2: The range of values when .

Algorithm 8.

 For  For               if                           record ;                     end         end end plot.

Lemma 9. If the analytical solution of the problem (5) , the spline functions and defined as (9) and (24) on the same grids are the numerical solutions of the problem (5) obtained by DSCM and ISCM, respectively; then we have the following conclusions.

(i)If , then , where is the analytical solution of the problem (5).(ii)Assume that and are the numerical solutions obtained by DSCM and ISCM for the problem

and the perturbed problem

where satisfies the Lipschitz condition (6), respectively. Then there exists the constant such that, for all , when , there exists the constant such that ; here is a constant, .

Proof. Firstly, we give the proof of the conclusion (i). Applying DSCM, we have

Since , we obtain

It follows from Lemma 4 that

where   as  .

By means of (70) and (72), we obtain

Applying Lemma 4, we get

where there exists a constant , such that . Based on the definitions of , it follows that

Based on the similar proof to that of Lemma 6, it follows from Lemma 3 and the inequality that

where

and is defined in Lemma 6, is a constant. By means of the convergence of the Mittag-Leffler function [6], we get

From Lemma 6, we have

Now, we give the proof of the conclusion (ii) of Lemma 6. According to the conclusion (i), we get

where , .

According to the proof of the conclusion (i), there exists such that, for all ,

Then,

Obviously, take ; hence, the conclusion (ii) is right.

Theorem 10 (convergence of DSCM). If the analytical solution of the problem (5) and the matrixes and of (14) satisfy that exists, and

then DSCM is convergent, and .

Proof. This result follows directly from Lemmas 6 and 9.

Now, we consider the stability of DSCM.

Theorem 11 (stability of DSCM). If the analytical solution of the problem (5) and the matrixes and in (14) satisfy that exists and

then ISCM is stable; that is if ISCM is, applied to solve the problem (68) and the perturbed problem (69), respectively, then, for all ( is a constant), one has

(i)

where

is a constant, are the numerical solutions of the problem (68), and are the numerical solutions of the perturbed problem (69), and(ii)

where is a constant, is the numerical solution of the problem (68), and is the numerical solution of the problem (69).

Moreover, DSCM is stable.

Proof. According to Lemma 9, if ISCM is stable, then DSCM is stable. In the following text, we give the proof of stability of ISCM.

By using ISCM to solve the problem (68) and the problem (69), we have

The previous equations can be written as

Thus,

Applying the Lipschitz condition, we get

where . Note that the definitions of . Let

where .

It follows from (92), (93), (49), and (51) that

Take ; hence, . According to the fact that is invertible and , we have

Moreover,

where . By means of (94), we have

Obviously, there exist and such that as . Then

Let . We obtain

Applying Lemma 3 yields

Using the convergence of Mittag-Leffler function and denoting

we obtain

Moreover, from the definitions of and , we have

By means of (103), we get

where ; that is, ISCM is stable. According to Lemma 9, DSCM is stable.

4. Illustrative Examples

In order to demonstrate our theoretical results, we apply DSCM to the problem (5) and present some numerical examples in this section. Let

where is the error function which depends on the stepsize .

Example 1. Consider the initial value problem [10, 21]

The exact solution is .

For different values of , the numerical solutions for problem (107) are obtained by using DSCM, fast wavelet collocation method (FWCM for short) reported in [21], and the method reported in [10]. When , the time stepsize ; the absolute errors of DSCM, FWCM, and the method reported in [10] are shown in Table 1. When , ; the absolute errors of these methods are shown in Table 2. The numerical solutions and the exact solution are shown in Figure 3. From the numerical results, the results obtained by DSCM are better than by FWCM and the method reported in [10] in terms of accuracy if the exact solution is sufficiently smooth. Therefore, DSCM is a valid method in solving fractional differential equation.

Table 1: The absolute errors of DSCM () and the methods reported in [10, 21]; .

Error by DSCM
Error by [21]
Error by [10]

Table 2: The absolute errors of DSCM () and the methods reported in [10, 21]; .

Error by DSCM
Error by [21]
Error by [10]

Figure 3: The error between the numerical solutions and the exact solution of Example 1; .

Example 2. Consider the initial value problem [12]

The exact solution is . Obviously, . In order to verify our theoretical results, we modify problem (108) as follows:

The exact solution is .

Tables 3 and 4 list the errors and the error orders of DSCM with and and show that the errors between the numerical solutions and the exact solution are very small, respectively. Figures 4 and 5 illustrate high accuracy of DSCM with , and show that the error is also very small. All of the numerical results show that DSCM for solving nonlinear FDEs is convergent and the method is robust.

Table 3: The error between the numerical solutions and the exact solution of the problem (109); , .

3.457022411

3.419757934

3.447443874

3.442757379

3.450810893

3.452123534

3.462413611

3.460568256

Table 4: The error between the numerical solutions and the exact solution of the problem (109); , .

3.342972530

3.264483517

3.261416159

3.252663415

3.252962303

3.247058188

3.249539075

3.246338153

(a)

(b)

Figure 4: The error between the numerical solutions and the exact solution of the problem (109); , .

(a)

(b)

Figure 5: The error between the numerical solutions and the exact solution of the problem (109); , .

Example 3. Consider the initial value problem

and the perturbed problem

We use the cubic spline collocation method to solve the problems (110) and (111), respectively. Selecting , , , , we obtain the numerical results given in Figure 6. When , , , , the numerical results are given in Figure 7. When , , , , the numerical results are shown in Figure 8.

(a)

(b)

Figure 6: The numerical results of the problems (110) and (111); , , , , , , .

(a)

(b)

Figure 7: The numerical results of the problems (110) and (111); , , , , , .

(a)

(b)

Figure 8: The numerical results of the problems (110) and (111); , , , , , .

From these figures, we can see that the absolute errors of the numerical solutions of the problems (110) and (111) decrease and finally tend to 0 as increases. Thus, we can draw the conclusion that the cubic spline collocation method for nonlinear FDEs is stable. The numerical results verify our theoretical results.

5. Conclusion

In this paper, the cubic spline collocation method with two parameters is successfully applied to the IVPs of general FDEs. The result of the local truncation error of this method is given. And the convergence and stability results of the cubic spline collocation method for the fractional order integral equations which is equivalent to the IVPs of general FDEs are obtained. By using the relationship between the numerical solutions from the cubic spline method for the IVPs of general FDEs and the numerical solutions obtained from the cubic spline method for the corresponding equivalent IVPs of fractional order integral equations, we also obtain some results of the convergence and the stability of the method for the IVPs of general FDEs. Some numerical examples successfully verify our theoretical results and show that the given method is efficient.

Acknowledgments

This work is supported by projects from NSF of China (no. 11271311, no. 11226320, and no. 61144004), Program for Changjiang Scholars and Innovative Research Team in University of China (no. IRT1179), the Aid Program for Science and Technology Innovative Research Team in Higher Educational Institutions of Hunan Province of China, Huizhou Science and Technology Program (no. 20110103), and NSF of Huizhou University (no. 2012YB15).

References

D. Delbosco and L. Rodino, “Existence and uniqueness for a nonlinear fractional differential equation,” Journal of Mathematical Analysis and Applications, vol. 204, no. 2, pp. 609–625, 1996.

J. H. He, “Some applications of nonlinear fractional differential equations and theirapproximations,” Bulletin of Science, Technology, vol. 15, no. 2, pp. 86–90, 1999.

V. V. Uchaikin and R. T. Sibatov, “Fractional theory for transport in disordered semiconductors,” Communications in Nonlinear Science and Numerical Simulation, vol. 13, no. 4, pp. 715–727, 2008.

V. E. Tarasov and G. M. Zaslavsky, “Fractional dynamics of systems with long-range interaction,” Communications in Nonlinear Science and Numerical Simulation, vol. 11, no. 8, pp. 885–898, 2006.

K. Diethelm, The Analysis of Fractional Differential Equations, vol. 2004 of Lecture Notes in Mathematics, Springer, Berlin, Germany, 2010.

I. Podlubny, Fractional Differential Equations, vol. 198 of Mathematics in Science and Engineering, Academic Press, New York, NY, USA, 1999.

Ch. Lubich, “Runge-Kutta theory for Volterra and Abel integral equations of the second kind,” Mathematics of Computation, vol. 41, no. 163, pp. 87–102, 1983.

Ch. Lubich, “Fractional linear multistep methods for Abel-Volterra integral equations of the second kind,” Mathematics of Computation, vol. 45, no. 172, pp. 463–469, 1985.

Ch. Lubich, “Discretized fractional calculus,” SIAM Journal on Mathematical Analysis, vol. 17, no. 3, pp. 704–719, 1986.

K. Diethelm, “An algorithm for the numerical solution of differential equations of fractional order,” Electronic Transactions on Numerical Analysis, vol. 5, pp. 1–6, 1997.

K. Diethelm and G. Walz, “Numerical solution of fractional order differential equations by extrapolation,” Numerical Algorithms, vol. 16, no. 3-4, pp. 231–253, 1997.

K. Diethelm, N. J. Ford, and A. D. Freed, “A predictor-corrector approach for the numerical solution of fractional differential equations,” Nonlinear Dynamics, vol. 29, no. 1–4, pp. 3–22, 2002.

K. Diethelm and N. J. Ford, “Numerical solution of the Bagley-Torvik equation,” BIT, vol. 42, no. 3, pp. 490–507, 2002.

F. Liu, V. Anh, and I. Turner, “Numerical solution of the space fractional Fokker-Planck equation,” Journal of Computational and Applied Mathematics, vol. 166, pp. 209–219, 2004.

F. Liu, P. Zhuang, V. Anh, I. Turner, and K. Burrage, “Stability and convergence of the difference methods for the space-time fractional advection-diffusion equation,” Applied Mathematics and Computation, vol. 191, no. 1, pp. 12–20, 2007.

A. Pedas and E. Tamme, “Spline collocation methods for linear multi-term fractional differential equations,” Journal of Computational and Applied Mathematics, vol. 236, no. 2, pp. 167–176, 2011.

A. Pedas and E. Tamme, “Spline collocation method for integro-differential equations with weakly singular kernels,” Journal of Computational and Applied Mathematics, vol. 197, no. 1, pp. 253–269, 2006.

A. Pedas and E. Tamme, “On the convergence of spline collocation methods for solving fractional differential equations,” Journal of Computational and Applied Mathematics, vol. 235, no. 12, pp. 3502–3514, 2011.

C. Li, A. Chen, and J. Ye, “Numerical approaches to fractional calculus and fractional ordinary differential equation,” Journal of Computational Physics, vol. 230, no. 9, pp. 3352–3368, 2011.

J. T. Kemppainen and K. M. Ruotsalainen, “On the spline collocation method for the single layer equation related to time-fractional diffusion,” Numerical Algorithms, vol. 57, no. 3, pp. 313–327, 2011.

X. Li, “Numerical solution of fractional differential equations using cubic B-spline wavelet collocation method,” Communications in Nonlinear Science and Numerical Simulation, vol. 17, no. 10, pp. 3934–3946, 2012.

F. Dubois and S. Mengué, “Mixed collocation for fractional differential equations,” Numerical Algorithms, vol. 34, no. 2–4, pp. 303–311, 2003.

E. A. Rawashdeh, “Numerical solution of semidifferential equations by collocation method,” Applied Mathematics and Computation, vol. 174, no. 2, pp. 869–876, 2006.

M. Weilbeer, Efficient Numerical Methods for Fractional Differential Equations and their Analytical Background, Papierflieger, 2006.

J. Dixon and S. McKee, “Weakly singular discrete Gronwall inequalities,” Zeitschrift für Angewandte Mathematik und Mechanik, vol. 66, no. 11, pp. 535–544, 1986.