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

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


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 .
The fluid flow is governed by the Navier-Stokes equations: where q denotes the fluid velocity, p denotes pressure, t denotes time, and ρ and v 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 α 0 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 D 2−α y x p x y g x , 0 ≤ α < 1, x ∈ a, b , 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 f x is continuous on the interval a, b and the operator D α represents the Caputo fractional derivative.Where, the Caputo fractional derivative is 22 when α 0, 1.3 is reduced to the classical second order boundary value problem.

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.

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: be a partition of a, b which divides the interval into n-equal parts.

International Journal of Mathematics and Mathematical Sciences
Cubic spline approximation will be built in each subinterval a ih, a i 1 h to approximate the solution of 2.1 -2.2 .Starting with the first interval a, a h , consider that the cubic polynomial spline segment S 0 x has the form: where a 0 , b 0 , c 0 , and d 0 are constants to be determined.It is straightforward to check: By construction, 2.4 satisfies 2.1 for x a.Then, for complete determination of the spline in the first interval, we have to find d 0 .From 2.4 , we have We will impose that the spline be a solution of the problem 2.1 at the point x a h.Hence, we obtain From 2.6 , 2.7 and using 2.4 we obtain:

2.8
Then the spline is fully determined in the first subinterval.In the next subinterval a h, a 2h the cubic spline segment S 1 x has the form: x − a h 3 .

2.9
From which we get

2.10
Taking into consideration that this cubic spline is of class C 2 a, a h ∪ a h, a 2h , and again all of the coefficients of S 1 x are determined with exception of d 1 .It is easy to check that the spline S 1 x be a solution of the problem 2.1 at the point x a h, then for determining d 1 we will impose that the spline be a solution of the problem 2.1 at the point x a 2h.Hence, by repeating the previous procedure we obtain Substituting by x a 2h into 2.10 and equating the result by 2.11 , we get .

2.12
By this way the spline is totally determined in the subinterval a h, a 2h .Iterating this process, let us consider that the cubic spline is constructed until the subinterval a i−1 h, a ih then we can define it in the next the subinterval a ih, a i 1 h as: where Then the cubic spline S x ∈ C 2 i j 0 a j, a j 1 h and easy to check that 2.13 verifies the differential equation 2.1 at the point x a ih.The constant d i can be determined by imposing that the spline be a solution of the problem 2.1 at the point x a i 1 h.Hence, we obtain

2.15
From 2.14 -2.15 , the spline approximation for the solutions of 1.3 and 2.1 at x i a ih, i 1, 2, . . ., n can be written in the following form:

2.19
For the subinterval a, a h :

2.20
Then, for i 1 in 2.19 , we get:

2.21
In general, it can be written as e i 1 e i O h 3 .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 z 1 and z 2 for y a .Using linear interpolation between y b in the two cases gives the next guess z 3 .Then problem 1.3 is resolved again with this new guess and so on.

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 y x can be written as: International Journal of Mathematics and Mathematical Sciences 7 We now use the binomial formula 9 : With 2.23 the expression for D α y x can be written as follows with λ m − α − 1: 2.24 The integrals: s p y m s ds, p 0, 1, 2 . . .

2.25
are solutions to the following system of differential equations: v p x p y m x , v p 0 0, p 0, 1, 2 . . .

2.26
According to 2.24 -2.26 the expression for D α y x can be rewritten as: with v p 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 N suitably chosen, so 2.27 will be

Convergence Analysis
Let S 3 Δ be the space of cubic splines with respect to Δ and with smoothness C 2 a, b .Also, let us denote by y Δ x the cubic spline approximation to y x .This implies that y Δ ∈ S 3  Δ which can be written as y Δ S i x , i 0, 1, 2 . . ., n − 1.
Without loss of generality, we will consider problem 1.3 with homogeneous Dirichlet boundary conditions 33 : y a 0, y b 0.

3.1
It will be assumed that y and y Δ satisfy these boundary conditions.

International Journal of Mathematics and Mathematical Sciences
If we assume that the BVP y x 0 along with boundary conditions 3.1 has a unique solution then there is a Green's function G x, s for the problems where Proof.From the Caputo fractional derivative D α y x , we get s a G t, s z s ds dt.

3.7
Using the principle of differentiation under the integral sign, for the function g x with the form: Φ x, t dt.

3.8
We have that where the functions Φ x, t and ∂/∂x Φ x, t are both continuous in both t and x in some region of the t, x plane, including δ 1 ≤ t ≤ δ 2 and x 0 ≤ x ≤ x 1 , then we can deduce that ∂ m ∂t m G t, s z s ds.

3.10
Then we have Changing the order of integration leads to We also introduce a linear projection P Δ that maps C 2 a, b to S 1 Δ piecewise linear interpolation at the grid points {x i } n 0 .Then 3.13 can be rewritten as: and we have also: By the definition of P Δ 33 , P Δ z − z ∞ converges to zero as h approaches zero for continuous function z x .This in turn implies that P Δ K − K ∞ converges to zero as h approaches zero.
where c k is a constant and independent of y, h and ψ y , h Proof.Let z x be a solution for 3.15 and y x be the solution of 1.3 -1.4 .Then, operating on both sides of 3.15 by the linear projection operator P Δ gives Adding z x to both sides of 3.18 and subtracting 3.16 from the results lead to

3.19
Operating on both sides of 3.19 by I P Δ K −1 leads to

3.20
Operating on both sides of 3.20 by the operator G and using 3.2 -3.4 , we get Since the operator G is bounded and from Theorem 3.2 the operator I P Δ K −1 is also bounded, then From 33 , we have that where, ψ y , h sup{|y τ h − y | : τ, τ h ∈ a, b , h ≤ h}.

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.The analytical solution of 4.1 , as found in 17 , has the following form: As k → 0, we may verify that the solution reduces to y x 1 − ∞ r 0 −1 r x 2 2r / 2 2r cos x.
This example occurs in the mathematical model of MEMS instrument 17 and had been solved for various values of k, 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 k 1/5 and k 0.005 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 .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.
Example 4.3.Consider the boundary value problem: The exact of solution 4.4 is y x x 2 .4.5 The numerical solutions using shooting method with z 1 1 and z 2 0.5 led to next guess of initial condition z 3 −0.03914,and the results are represented in Table 3.
Example 4.4.Consider the boundary value problem: The exact of solution 4.6 is The numerical solutions use shooting method for θ 0.5, β 1 and with initial guesses z 1 1 and z 2 0.5, z 3 −0.01224;these lead to z 4 −5.07E− 4 as the fourth guess for the initial condition and the results are represented in Table 4.

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.
s z Δ s ds Gz Δ x , 3.4

Lemma 3 . 1 .
Consider the following: x, s z s ds D α Gz x .3.6 s dt z s ds, the operator Ky x defined by: maps C 2 a, b to C a, b .

Theorem 3 . 2 Theorem 3 . 3 .
see 34 .If there is N 0 large enough, then { I P Δ K −1 : n ≥ N 0 } exists and consists of a sequence of bounded linear operators.Which means, for a constant δ independent of N 0 and z ∈ C a, b , if n ≥ N 0 , then I P Δ K −1 z ≤ δ z .Assuming that H1 the BVP 1.3 along with boundary conditions 3.1 has a unique solution in C 2 a, b , H2 the BVP y x 0 along with boundary conditions 3.1 has a unique solution, then, for some n ≥ N 0 one has

Example 4 . 1 .
Consider the initial value problem:

Table 1 :
Numerical results of Example 4.1.

Table 2 :
Numerical results of Example 4.2.

Table 3 :
Numerical results of Example 4.3.

Table 4 :
Numerical results of Example 4.4.