International Journal of Mathematics and Mathematical Sciences

Volume 2012 (2012), Article ID 638026, 16 pages

doi:10.1155/2012/638026

Research Article

The Use of Cubic Splines in the Numerical Solution of Fractional Differential Equations

W. K. Zahra¹ and S. M. Elkholy²

¹Department of Physics and Engineering Mathematics, Faculty of Engineering, Tanta University, 31521 Tanta, Egypt

²Department of Engineering Physics and Mathematics, Faculty of Engineering, Kafr El Sheikh University, Kafr El Sheikh, Egypt

Received 31 March 2012; Revised 8 May 2012; Accepted 23 May 2012

Academic Editor: Manfred Moller

Copyright © 2012 W. K. Zahra and S. M. Elkholy. 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.

Fractional calculus became a vital tool in describing many phenomena appeared in physics, chemistry as well as engineering fields. Analytical solution of many applications, where the fractional differential equations appear, cannot be established. Therefore, cubic polynomial spline-function-based method combined with shooting method is considered to find approximate solution for a class of fractional boundary value problems (FBVPs). Convergence analysis of the method is considered. Some illustrative examples are presented.

1. Introduction

Fractional calculus attracted the attention of many researchers because it has recently gained popularity in the investigation of dynamical systems. There are many applications of fractional derivative and fractional integration in several complex systems such as physics, chemistry, fluid mechanics, viscoelasticity, signal processing, mathematical biology, and bioengineering, and various applications in many branches of science and engineering could be found [1–16].

One of the applications where fractional differential equation appears is the equation describing the motion of fluids, which are encountered downhole during the process of oil well logging, through a device that has been designed to measure fluids viscosity. In the oil exploration industry, the fluid viscosity can indicate the permeability of the reservoir formation, its flow characteristics, and the commercial value of the reservoir fluid. So, viscometers are required to measure the thermophysical properties of these fluids. It is hard to simulate reservoir conditions in a laboratory because a reservoir can exhibit temperatures of 20–200°C and pressures of 5–200 MPa. Therefore, a microelectromechanical system (MEMS) instrument has been designed to measure the viscosity of fluids which contains only a single moving part (all others being electrical). This device can operate at high ambient pressures, and the behavior of the device may be analyzed in a manner that allows its design to be optimised see Figure 1, [12, 17].

Figure 1: MEMS instrument.

The fluid flow is governed by the Navier-Stokes equations:

where denotes the fluid velocity, denotes pressure, denotes time, and and are the fluid density and kinematic viscosity, respectively. Then, it was found that the equation governing the motion of the fluid through the instrument is

The above fractional differential equation is well known as Bagley-Torvik equation when which appears in modeling the motion of a rigid plate immersed in a Newtonian fluid [12, 17].

Several methods have been proposed to obtain the analytical solution of fractional differential equations (FDEs) such as Laplace and Fourier transforms, eigenvector expansion, method based on Laguerre integral formula, direct solution based on Grunewald Letnikov approximation, truncated Taylor series expansion, and power series method [9, 18–23]. There are also several methods have recently been proposed to solve FDEs numerically such as fractional Adams-Moulton methods, explicit Adams multistep methods, fractional difference method, decomposition method, variational iteration method, least squares finite element solution, extrapolation method, and the Kansa method which is meshless, easy-to-use, and has been used to handle a broad range of partial differential equation models [24–31]. Also, the authors considered the numerical solution of the fractional boundary value problem (FBVP) , , with Dirichlet boundary conditions using quadratic polynomial spline, [32].

The existence of at least one solution of fractional problems can be seen in [3, 11, 14, 16, 31].

We consider the numerical solution of the following fractional boundary value problem (FBVPs):

Subject to boundary conditions:

where the function is continuous on the interval and the operator represents the Caputo fractional derivative. Where, the Caputo fractional derivative is [22]

when , (1.3) is reduced to the classical second order boundary value problem.

2. Method of Solution

The following is a brief derivation of the algorithm used to solve problem (1.3)-(1.4). The method of solution presented in the following section is based on cubic spline approach combined with shooting method.

2.1. Cubic Spline Solution for FDEs

In order to develop cubic spline approximation for the fractional differential equation (1.3)-(1.4), we would discuss the solution of (1.3) as initial value problem of the form:

Let

be a partition of which divides the interval into -equal parts.

Cubic spline approximation will be built in each subinterval to approximate the solution of (2.1)-(2.2). Starting with the first interval , consider that the cubic polynomial spline segment has the form:

where , and are constants to be determined. It is straightforward to check:

By construction, (2.4) satisfies (2.1) for . Then, for complete determination of the spline in the first interval, we have to find . From (2.4), we have

We will impose that the spline be a solution of the problem (2.1) at the point . Hence, we obtain

From (2.6), (2.7) and using (2.4) we obtain:

Then the spline is fully determined in the first subinterval. In the next subinterval the cubic spline segment has the form:

From which we get

Taking into consideration that this cubic spline is of class , and again all of the coefficients of are determined with exception of . It is easy to check that the spline be a solution of the problem (2.1) at the point , then for determining we will impose that the spline be a solution of the problem (2.1) at the point . Hence, by repeating the previous procedure we obtain

Substituting by into (2.10) and equating the result by (2.11), we get

By this way the spline is totally determined in the subinterval . Iterating this process, let us consider that the cubic spline is constructed until the subinterval then we can define it in the next the subinterval as:

where

Then the cubic spline and easy to check that (2.13) verifies the differential equation (2.1) at the point . The constant can be determined by imposing that the spline be a solution of the problem (2.1) at the point . Hence, we obtain

From (2.14)-(2.15), the spline approximation for the solutions of (1.3) and (2.1) at can be written in the following form:

Lemma 2.1. Let then the error bound associated with (2.16) is .

Proof. For each subinterval , the error terms are

Using, Taylor expansion for , we get

Then (2.16) and (2.18) led to

For the subinterval :

Then, for in (2.19), we get:

In general, it can be written as . Then, it can be proved that .

Cubic spline method presented above can be extended to solve fractional boundary value problems by implementing the shooting method. FBVPs (1.3) with boundary conditions (1.4) will be solved as initial value problem with two guesses for . Using linear interpolation between in the two cases gives the next guess . Then problem (1.3) is resolved again with this new guess and so on.

2.2. Numerical Approximation of Fractional Term

The algorithm used for solving fractional differential equation is based on transforming the fractional derivative into a system of ordinary differential equation. Firstly, the Caputo fractional derivative for can be written as:

We now use the binomial formula [9]:

With (2.23) the expression for can be written as follows with :

The integrals:

are solutions to the following system of differential equations:

According to (2.24)–(2.26) the expression for can be rewritten as:

with satisfying (2.26), (2.27) will represent the fundamental relation used in numerical representation of the fractional term in fractional differential equations. In application, we will use finite number of terms suitably chosen, so (2.27) will be

3. Convergence Analysis

Let be the space of cubic splines with respect to and with smoothness . Also, let us denote by the cubic spline approximation to . This implies that which can be written as .

Without loss of generality, we will consider problem (1.3) with homogeneous Dirichlet boundary conditions [33]:

It will be assumed that satisfy these boundary conditions.

If we assume that the BVP along with boundary conditions (3.1) has a unique solution then there is a Green’s function for the problems

where

where

is a compact operator, since is continuous in , [33].

Lemma 3.1. Consider the following:

Proof. From the Caputo fractional derivative , we get

Using the principle of differentiation under the integral sign, for the function with the form:

We have that

where the functions and are both continuous in both and in some region of the plane, including and , then we can deduce that

Then we have

Changing the order of integration leads to

and this proves the lemma.

Substituting from (3.2)–(3.4) and (3.6) into (1.3) leads to

We will introduce the operator defined by:

which maps to . We also introduce a linear projection that maps to piecewise linear interpolation at the grid points . Then (3.13) can be rewritten as:

and we have also:

By the definition of [33], converges to zero as approaches zero for continuous function . This in turn implies that converges to zero as approaches zero.

Theorem 3.2 (see [34]). If there is large enough, then exists and consists of a sequence of bounded linear operators. Which means, for a constant independent of and , if , then .

Theorem 3.3. Assuming that

(H1) the BVP (1.3) along with boundary conditions (3.1) has a unique solution in , (H2) the BVP along with boundary conditions (3.1) has a unique solution,  then, for some one has

where is a constant and independent of , and .

Proof. Let be a solution for (3.15) and be the solution of (1.3)-(1.4). Then, operating on both sides of (3.15) by the linear projection operator gives

Adding to both sides of (3.18) and subtracting (3.16) from the results lead to

Operating on both sides of (3.19) by leads to

Operating on both sides of (3.20) by the operator and using (3.2)–(3.4), we get

Since the operator is bounded and from Theorem 3.2 the operator is also bounded, then

From [33], we have that

where, .

Substituting from Theorem 3.2 and (3.23) into (3.22) completes the proof.

4. Numerical Examples

We will consider some numerical examples illustrating the solution using cubic spline methods. All calculations are implemented with MATLAB 7, and we used implicit Adams-Bashforth three-step method in approximating the fractional term.

Example 4.1. Consider the initial value problem:

The analytical solution of (4.1), as found in [17], has the following form:

As , we may verify that the solution reduces to .

This example occurs in the mathematical model of (MEMS) instrument [17] and had been solved for various values of , and the solutions are represented in Figures 2–4.

Figures 3 and 4 represent a comparison between our approximate solutions and the analytical solutions for and respectively. The results are tabulated also in Table 1.

The obtained results have good agreement with the exact solution as in Figures (3 and 4) and Table 1 and those published in [17].

Table 1: Numerical results of Example 4.1.

= 1
= 1/5
= 0.005

Analytical
Approx
Analytical
Approx
Analytical
Approx

solution
solution
solution
solution
solution
solution

0
1
1
1
1
1
1

0.125
0.99437
0.993126
0.992747
0.992391
0.992212
0.992212

0.250
0.979919
0.974802
0.971922
0.970148
0.968995
0.968983

0.375
0.958424
0.944545
0.938558
0.933609
0.930733
0.930674

0.500
0.930957
0.904813
0.893615
0.883958
0.878038
0.877899

0.625
0.898335
0.857938
0.838087
0.822499
0.811743
0.811497

0.750
0.861241
0.805442
0.773025
0.750552
0.732892
0.732514

0.875
0.820277
0.748795
0.699540
0.669584
0.642719
0.642193

1
0.775989
0.688838
0.618798
0.580978
0.542633
0.541945

Figure 2: Numerical solutions of Example 4.1 for (a) , (b) , and (c) .

Figure 3: A comparison between the analytical solution and our approximated solution for .

Figure 4: A comparison between the analytical solution and our approximated solution for .

Example 4.2. Consider the initial value problem:

This example had been solved for many methods. Table 2 shows a comparison between the solution of (4.3) by our method, decomposition method and fractional difference method.

Table 2: Numerical results of Example 4.2.

Decomp. Method [29]
Fractional diff. method [29]
Our method

0
0
0
0

0.1
0.039874
0.039473
0.039933

0.2
0.158512
0.157703
0.158981

0.3
0.353625
0.352402
0.353996

0.4
0.622083
0.622083
0.619900

0.5
0.960047
0.957963
0.950455

0.6
1.363093
1.360551
1.348551

0.7
1.826257
1.823267
1.796370

0.8
2.344224
2.340749
2.295551

0.9
2.911278
2.907324
2.899808

1
3.521462
3.517013
3.499200

Example 4.3. Consider the boundary value problem:

The exact of solution (4.4) is

The numerical solutions using shooting method with and led to next guess of initial condition , and the results are represented in Table 3.

Table 3: Numerical results of Example 4.3.

Exact solution
Approx. with
Error

0
0
0
0

0.125
0.015625
0.012842
−3

0.250
0.062500
0.053552
−3

0.375
0.140630
0.138569
−3

0.500
0.250000
0.246192
−3

0.625
0.390630
0.386836
−3

0.750
0.562500
0.553883
−3

0.875
0.765630
0.758442
−3

1
1
0.999971
−5

Example 4.4. Consider the boundary value problem:

The exact of solution (4.6) is

The numerical solutions use shooting method for and with initial guesses and , ; these lead to as the fourth guess for the initial condition and the results are represented in Table 4.

Table 4: Numerical results of Example 4.4.

Exact solution
Approx. solution
Error

0.125
−0.00021
−0.00039
−4

0.250
−0.00292
−0.00346
−4

0.375
−0.01235
−0.01316
−4

0.500
−0.03125
−0.03192
−4

0.625
−0.05722
−0.05713
−5

0.750
−0.07910
−0.07732
−3

0.875
−0.07327
−0.06985
−3

1
0
−4
−4

5. Conclusion

New scheme for solving class of fractional boundary value problem is presented using cubic spline method combined with shooting method. Transforming the fractional derivative into a system of ordinary differential equations is used for approximating the fractional term. Implicit Adams-Bashforth three-step method has been used for approximating this system of ordinary differential equations. Convergence analysis of the method is considered and is shown to be second order. Numerical comparisons between the solution using this new method and the methods introduced in [17, 29] are presented. The obtained numerical results show that the proposed method maintains a remarkable high accuracy which makes it encouraging for dealing with the solution of two-point boundary value problem of fractional order.

Acknowledgment

The authors are grateful to the referees for their suggestions and valuable comments.

References

O. P. Agrawal and P. Kumar, “Comparison of five schemes for fractional differential equations,” in Advances in Fractional Calculus: Theoretical Developments and Applications in Physics and Engineering, J. Sabatier, O. P. Agrawal, and J. A. Tenreiro Machado, Eds., pp. 43–60, 2007.

J. H. Ahlberg, E. N. Nilson, and J. L. Walsh, The Theory of Splines and Their Applications, Academic Press, New York, NY, USA, 1967.

D. Baleanu and S. I. Muslih, “On Fractional Variational Principles,” in Advances in Fractional Calculus: Theoretical Developments and Applications in Physics and Engineering, J. Sabatier, O. P. Agrawal, and J. A. Tenreiro Machado, Eds., pp. 115–126, 2007.

M. M. Benghorbal, Power series solutions of fractional differential equations and symbolic derivatives and integrals [Ph.D. thesis], Faculty of Graduate Studies, The University of Western Ontario, Ontario, Canada, 2004.

B. Bonilla, M. Rivero, and J. J. Trujillo, “Linear differential equations of fractional order,” in Advances in Fractional Calculus: Theoretical Developments and Applications in Physics and Engineering, J. Sabatier, O. P. Agrawal, and J. A. Tenreiro Machado, Eds., pp. 77–91, 2007.

C. X. Jiang, J. E. Carletta, and T. T. Hartley, “Implementation of fractional-order operators on field programmable gate arrays,” in Advances in Fractional Calculus: Theoretical Developments and Applications in Physics and Engineering, J. Sabatier, O. P. Agrawal, and J. A. Tenreiro Machado, Eds., pp. 333–346, 2007.

N. Kosmatov, “Integral equations and initial value problems for nonlinear differential,” Nonlinear Analysis: Theory, Methods & Applications, vol. 70, no. 7, pp. 2521–2529, 2009.

V. Lakshmikantham and A. S. Vatsala, “Basic theory of fractional differential equations,” Nonlinear Analysis:Theory, Methods & Applications, vol. 69, no. 8, pp. 2677–2682, 2008.

K. S. Miller and B. Ross, An Introduction to the Fractional Calculus and Differential Equations, John Wiley & Sons, New York, NY, USA, 1993.

H. Nasuno, N. Shimizu, and M. Fukunaga, “Fractional derivative consideration on nonlinear viscoelastic statical and dynamical behavior under large pre-displacement,” in Advances in Fractional Calculus: Theoretical Developments and Applications in Physics and Engineering, J. Sabatier, O. P. Agrawal, and J. A. Tenreiro Machado, Eds., pp. 363–376, 2007.

A. Ouahab, “Some results for fractional boundary value problem of differential inclusions,” Nonlinear Analysis: Theory, Methods & Applications, vol. 69, no. 11, pp. 3877–3896, 2008.

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

I. Podlubny, I. Petráš, B. M. Vinagre, P. O'Leary, and L'. Dorčák, “Analogue realizations of fractional-order controllers,” Nonlinear Dynamics, vol. 29, no. 1–4, pp. 281–296, 2002.

X. Su and S. Zhang, “Solutions to boundary-value problems for nonlinear differential equations of fractional order,” Electronic Journal of Differential Equations, vol. 2009, no. 26, pp. 1–15, 2009.

M. S. Tavazoei and M. Haeri, “A note on the stability of fractional order systems,” Mathematics and Computers in Simulation, vol. 79, no. 5, pp. 1566–1576, 2009.

X. Su, “Boundary value problem for a coupled system of nonlinear fractional differential equations,” Applied Mathematics Letters, vol. 22, no. 1, pp. 64–69, 2009.

A. D. Fitt, A. R. H. Goodwin, K. A. Ronaldson, and W. A. Wakeham, “A fractional differential equation for a MEMS viscometer used in the oil industry,” Journal of Computational and Applied Mathematics, vol. 229, no. 2, pp. 373–381, 2009.

J. Duan, J. An, and M. Xu, “Solution of system of fractional differential equations by Adomian decomposition method,” Applied Mathematics, A Journal of Chinese Universities Series B, vol. 22, no. 1, pp. 7–12, 2007.

R. Garrappa, “On some explicit Adams multistep methods for fractional differential equations,” Journal of Computational and Applied Mathematics, vol. 229, no. 2, pp. 392–399, 2009.

A. Ghorbani, “Toward a new analytical method for solving nonlinear fractional differential equations,” Computer Methods in Applied Mechanics and Engineering, vol. 197, no. 49-50, pp. 4173–4179, 2008.

E. R. Kaufmann and E. Mboumi, “Positive solutions of a boundary value problem for a nonlinear fractional differential equation,” Electronic Journal of Qualitative Theory of Differential Equations, no. 3, pp. 1–11, 2008.

A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo, Theory and Applications of Fractional Differential Equations, Minsk, Belarus, 1st edition, 2006.

S. Momani and Z. Odibat, “A novel method for nonlinear fractional partial differential equations: combination of DTM and generalized Taylor's formula,” Journal of Computational and Applied Mathematics, vol. 220, no. 1-2, pp. 85–95, 2008.

W. Chen, H. Sun, X. Zhang, and D. Korošak, “Anomalous diffusion modeling by fractal and fractional derivatives,” Computers & Mathematics with Applications, vol. 59, no. 5, pp. 1754–1758, 2010.

W. Chen, L. Ye, and H. Sun, “Fractional diffusion equations by the Kansa method,” Computers & Mathematics with Applications, vol. 59, no. 5, pp. 1614–1620, 2010.

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.

G. J. Fix and J. P. Roop, “Least squares finite-element solution of a fractional order two-point boundary value problem,” Computers & Mathematics with Applications, vol. 48, no. 7-8, pp. 1017–1033, 2004.

L. Galeone and R. Garrappa, “Fractional Adams-Moulton methods,” Mathematics and Computers in Simulation, vol. 79, no. 4, pp. 1358–1367, 2008.

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.

J. P. Roop, Variational solution of the fractional advection dispersion equation [Ph.D. thesis], Clemson University, Clemson, SC, USA, 2004.

F. I. Taukenova and M. Kh. Shkhanukov-Lafishev, “Difference methods for solving boundary value problems for fractional-order differential equations,” Computational Mathematics and Mathematical Physics, vol. 46, no. 10, pp. 1871–1795, 2006.

W. K. Zahra and S. M. Elkholy, “Quadratic spline solution for boundary value problem of fractional order,” Numerical Algorithms, vol. 59, pp. 373–391, 2012.

P. M. Prenter, Splines and Variational Methods, John Wiley & Sons, New York, NY, USA, 1975.

R. D. Russell and L. F. Shampine, “A collocation method for boundary value problems,” Numerische Mathematik, vol. 19, pp. 1–28, 1972.