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

Authors: W. K. Zahra; S. M. Elkholy
Journal: International Journal of Mathematics and Mathematical Sciences
 (2012)

---

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


---

## Body

## 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 Figure1, [12, 17].Figure 1
MEMS instrument.The fluid flow is governed by the Navier-Stokes equations:(1.1)qt+(q⋅∇)q=-1ρ∇p+v∇2q,∇⋅q=0,
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
(1.2)y′′(x)+kπD1.5y(x)+αy′(x)=0,y(0)=1,y′(0)=0.
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) D2-α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):(1.3)y′′(x)+θDαy+βy=f(x),m-1≤α<m,x∈[a,b].
Subject to boundary conditions:
(1.4)y(a)=ya,y(b)=yb,
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]
(1.5)Dαy(x)=1Γ(m-α)∫0x(x-s)m-α-1y(m)(s)ds,α>0,m-1<α<m,
when α=0, (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:
(2.1)y′′(x)+θDαy+βy=f(x),0≤α<1,x∈[a,b],(2.2)y(a)=ya,y′(a)=ya′.Let(2.3)Δn:xi=a+ih,x0=a,xn=b,h=(b-a)n,i=0,1,2,…,n-1
be a partition of [a,b] which divides the interval into n-equal parts.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 S0(x) has the form:
(2.4)S0(x)=a0+b0(x-a)+c02(x-a)2+d06(x-a)3,
where a0,b0,c0, and d0 are constants to be determined. It is straightforward to check:
(2.5)S0(a)=a0=ya,S0′(a)=b0=ya′,S0′′(a)=c0=y′′(a)=f(a)-βya.
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 d0. From (2.4), we have
(2.6)S0′′(x)=c0+d0(x-a).
We will impose that the spline be a solution of the problem (2.1) at the point x=a+h. Hence, we obtain
(2.7)S0′′(a+h)=y′′(a+h)=f(a+h)-βy(a+h)-θDαy|x=a+h.
From (2.6), (2.7) and using (2.4) we obtain:(2.8)(h+βh36)d0=f(a+h)-β(ya+ya′h+h22y′′(a))-y′′(a)-θDαy|x=a+h.
Then the spline is fully determined in the first subinterval. In the next subinterval [a+h,a+2h] the cubic spline segment S1(x) has the form:
(2.9)S1(x)=S0(a+h)+S0′(a+h)(x-(a+h))+S0′′(a+h)2(x-(a+h))2+d16(x-(a+h))3.
From which we get
(2.10)S1′′(x)=S0′′(a+h)+d1(x-(a+h)).
Taking into consideration that this cubic spline is of class C2([a,a+h]∪[a+h,a+2h]), and again all of the coefficients of S1(x) are determined with exception of d1. It is easy to check that the spline S1(x) be a solution of the problem (2.1) at the point x=a+h, then for determining d1 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
(2.11)S1′′(a+2h)=y′′(a+2h)=f(a+2h)-βy(a+2h)-θDαy|x=a+2h.
Substituting by x=a+2h into (2.10) and equating the result by (2.11), we get
(2.12)(h+βh36)d1=f(a+2h)-β(S0(a+h)+hS0′(a+h)+h22S0′′(a+h))-S0′′(a+h)-θDαy|x=a+2h.
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:
(2.13)Si(x)=ψi+di6(x-(a+ih))3,
where
(2.14)ψi(x)=∑k=021k!Sk-1(k)(a+ih)(x-(a+ih))k.
Then the cubic spline S(x)∈C2(⋃j=0i[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 di 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)(h+βh36)di=f(a+(i+1)h)-βψi(a+(i+1)h)-ψi′′(a+(i+1)h)-θDαy|x=a+(i+1)h.
From (2.14)-(2.15), the spline approximation for the solutions of (1.3) and (2.1) at xi=a+ih,i=1,2,…,n can be written in the following form:
(2.16)Si(xi+1)=∑k=021k!hkSi-1(k)(a+ih)+h36di,di=1h[Si′′(a+(i+1)h)-Si-1′′(a+ih)],i=0,1,2,….Lemma 2.1.
Lety∈C4[a,b] then the error bound associated with (2.16) is |e(x)|=O(h2).Proof.
For each subinterval[a+ih,a+(i+1)h], the error terms are
(2.17)ei+1=y(xi+1)-Si(xi+1),i=0,1,2,…,n-1.
Using, Taylor expansion for y(xi+1), we get
(2.18)y(a+(i+1)h)=y(a+ih)+hy′(a+ih)+h22y′′(a+ih)+h36y′′′(a+ih)+O(h4).
Then (2.16) and (2.18) led to
(2.19)ei+1=y(xi+1)-Si(xi+1)=ei+hei′+(h22)ei′′+(h36)ei′′′+O(h4),i=1,2,…,n-1.
For the subinterval [a,a+h]:
(2.20)e1=y(a+h)-S0(a+h)=h3[y′′′(a)-d0]6+O(h4)=O(h3),e1′=y′(a+h)-S0′(a+h)=O(h2),e1′′=y′′(a+h)-S0′′(a+h)=O(h).
Then, for i=1 in (2.19), we get:
(2.21)e2=y(a+2h)-S1(a+2h)=e1+he1′+(h22)e1′′+(h36)e1′′′+O(h4)=e1+O(h3)=O(h3).
In general, it can be written as ei+1=ei+O(h3). Then, it can be proved that |e(x)|=nO(h3)=O(h2).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 z1andz2 for y′(a). Using linear interpolation between y(b) in the two cases gives the next guess z3. 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 fory(x) can be written as:
(2.22)Dαy(x)=xm-α-1Γ(m-α)∫0x(1-sx)m-α-1y(m)(s)ds,α>0,m-1<α<m.
We now use the binomial formula [9]:
(2.23)(1+z)λ=∑p=0∞(λp)zp=∑p=0∞(-1)pΓ(p-λ)Γ(-λ)p!zp,|z|<1.
With (2.23) the expression for Dαy(x) can be written as follows with λ=m-α-1:
(2.24)Dαy(x)=xλΓ(λ+1)∫0xy(m)(s)[∑p=0∞Γ(p-λ)Γ(-λ)p!(sx)p]ds,α>0,m-1<α<m.
The integrals:
(2.25)vp=∫0xspy(m)(s)ds,p=0,1,2…
are solutions to the following system of differential equations:
(2.26)vp′=xpy(m)(x),vp(0)=0,p=0,1,2…
According to (2.24)–(2.26) the expression for Dαy(x) can be rewritten as:
(2.27)Dαy(x)=xm-α-1Γ(m-α)∑p=0∞(Γ(p-m+α+1)Γ(-m+α+1)p!xpvp),α>0,m-1<α<m,
with vp 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
(2.28)Dαy(x)≅xm-α-1Γ(m-α)∑p=0N(Γ(p-m+α+1)Γ(-m+α+1)p!xpvp),α>0,m-1<α<m.


## 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:
(2.1)y′′(x)+θDαy+βy=f(x),0≤α<1,x∈[a,b],(2.2)y(a)=ya,y′(a)=ya′.Let(2.3)Δn:xi=a+ih,x0=a,xn=b,h=(b-a)n,i=0,1,2,…,n-1
be a partition of [a,b] which divides the interval into n-equal parts.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 S0(x) has the form:
(2.4)S0(x)=a0+b0(x-a)+c02(x-a)2+d06(x-a)3,
where a0,b0,c0, and d0 are constants to be determined. It is straightforward to check:
(2.5)S0(a)=a0=ya,S0′(a)=b0=ya′,S0′′(a)=c0=y′′(a)=f(a)-βya.
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 d0. From (2.4), we have
(2.6)S0′′(x)=c0+d0(x-a).
We will impose that the spline be a solution of the problem (2.1) at the point x=a+h. Hence, we obtain
(2.7)S0′′(a+h)=y′′(a+h)=f(a+h)-βy(a+h)-θDαy|x=a+h.
From (2.6), (2.7) and using (2.4) we obtain:(2.8)(h+βh36)d0=f(a+h)-β(ya+ya′h+h22y′′(a))-y′′(a)-θDαy|x=a+h.
Then the spline is fully determined in the first subinterval. In the next subinterval [a+h,a+2h] the cubic spline segment S1(x) has the form:
(2.9)S1(x)=S0(a+h)+S0′(a+h)(x-(a+h))+S0′′(a+h)2(x-(a+h))2+d16(x-(a+h))3.
From which we get
(2.10)S1′′(x)=S0′′(a+h)+d1(x-(a+h)).
Taking into consideration that this cubic spline is of class C2([a,a+h]∪[a+h,a+2h]), and again all of the coefficients of S1(x) are determined with exception of d1. It is easy to check that the spline S1(x) be a solution of the problem (2.1) at the point x=a+h, then for determining d1 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
(2.11)S1′′(a+2h)=y′′(a+2h)=f(a+2h)-βy(a+2h)-θDαy|x=a+2h.
Substituting by x=a+2h into (2.10) and equating the result by (2.11), we get
(2.12)(h+βh36)d1=f(a+2h)-β(S0(a+h)+hS0′(a+h)+h22S0′′(a+h))-S0′′(a+h)-θDαy|x=a+2h.
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:
(2.13)Si(x)=ψi+di6(x-(a+ih))3,
where
(2.14)ψi(x)=∑k=021k!Sk-1(k)(a+ih)(x-(a+ih))k.
Then the cubic spline S(x)∈C2(⋃j=0i[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 di 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)(h+βh36)di=f(a+(i+1)h)-βψi(a+(i+1)h)-ψi′′(a+(i+1)h)-θDαy|x=a+(i+1)h.
From (2.14)-(2.15), the spline approximation for the solutions of (1.3) and (2.1) at xi=a+ih,i=1,2,…,n can be written in the following form:
(2.16)Si(xi+1)=∑k=021k!hkSi-1(k)(a+ih)+h36di,di=1h[Si′′(a+(i+1)h)-Si-1′′(a+ih)],i=0,1,2,….Lemma 2.1.
Lety∈C4[a,b] then the error bound associated with (2.16) is |e(x)|=O(h2).Proof.
For each subinterval[a+ih,a+(i+1)h], the error terms are
(2.17)ei+1=y(xi+1)-Si(xi+1),i=0,1,2,…,n-1.
Using, Taylor expansion for y(xi+1), we get
(2.18)y(a+(i+1)h)=y(a+ih)+hy′(a+ih)+h22y′′(a+ih)+h36y′′′(a+ih)+O(h4).
Then (2.16) and (2.18) led to
(2.19)ei+1=y(xi+1)-Si(xi+1)=ei+hei′+(h22)ei′′+(h36)ei′′′+O(h4),i=1,2,…,n-1.
For the subinterval [a,a+h]:
(2.20)e1=y(a+h)-S0(a+h)=h3[y′′′(a)-d0]6+O(h4)=O(h3),e1′=y′(a+h)-S0′(a+h)=O(h2),e1′′=y′′(a+h)-S0′′(a+h)=O(h).
Then, for i=1 in (2.19), we get:
(2.21)e2=y(a+2h)-S1(a+2h)=e1+he1′+(h22)e1′′+(h36)e1′′′+O(h4)=e1+O(h3)=O(h3).
In general, it can be written as ei+1=ei+O(h3). Then, it can be proved that |e(x)|=nO(h3)=O(h2).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 z1andz2 for y′(a). Using linear interpolation between y(b) in the two cases gives the next guess z3. 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 fory(x) can be written as:
(2.22)Dαy(x)=xm-α-1Γ(m-α)∫0x(1-sx)m-α-1y(m)(s)ds,α>0,m-1<α<m.
We now use the binomial formula [9]:
(2.23)(1+z)λ=∑p=0∞(λp)zp=∑p=0∞(-1)pΓ(p-λ)Γ(-λ)p!zp,|z|<1.
With (2.23) the expression for Dαy(x) can be written as follows with λ=m-α-1:
(2.24)Dαy(x)=xλΓ(λ+1)∫0xy(m)(s)[∑p=0∞Γ(p-λ)Γ(-λ)p!(sx)p]ds,α>0,m-1<α<m.
The integrals:
(2.25)vp=∫0xspy(m)(s)ds,p=0,1,2…
are solutions to the following system of differential equations:
(2.26)vp′=xpy(m)(x),vp(0)=0,p=0,1,2…
According to (2.24)–(2.26) the expression for Dαy(x) can be rewritten as:
(2.27)Dαy(x)=xm-α-1Γ(m-α)∑p=0∞(Γ(p-m+α+1)Γ(-m+α+1)p!xpvp),α>0,m-1<α<m,
with vp 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
(2.28)Dαy(x)≅xm-α-1Γ(m-α)∑p=0N(Γ(p-m+α+1)Γ(-m+α+1)p!xpvp),α>0,m-1<α<m.


## 3. Convergence Analysis

LetSΔ3 be the space of cubic splines with respect to Δ and with smoothness C2[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Δ=Si(x),i=0,1,2…,n-1.Without loss of generality, we will consider problem (1.3) with homogeneous Dirichlet boundary conditions [33]:
(3.1)y(a)=0,y(b)=0.
It will be assumed that yandyΔ satisfy these boundary conditions.If we assume that the BVPy′′(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
(3.2)z=y′′,zΔ=yΔ′′,
where
(3.3)y(x)=∫abG(x,s)z(s)ds=Gz(x),(3.4)yΔ(x)=∫abG(x,s)zΔ(s)ds=GzΔ(x),
where
(3.5)G(x,s)={(x-s)-(x-a)(b-s)(b-a),a≤s≤x≤b,-(x-a)(b-s)(b-a),a≤x≤s≤b,G is a compact operator, since G(x,s) is continuous in [a,b]×[a,b], [33].Lemma 3.1.
Consider the following:(3.6)Dαy(x)=Dα∫abG(x,s)z(s)ds=∫ab(DαG(x,s))z(s)ds=DαGz(x).Proof.
From the Caputo fractional derivativeDαy(x), we get
(3.7)Dαy(x)=Dα∫s=as=bG(x,s)z(s)ds=1Γ(m-α)∫t=at=x(x-t)m-α-1(dmdtm[∫s=as=bG(t,s)z(s)ds])dt.
Using the principle of differentiation under the integral sign, for the function g(x) with the form:
(3.8)g(x)=∫δ1(x)δ2(x)Φ(x,t)dt.
We have that
(3.9)dg(x)dx=∫δ1(x)δ2(x)∂∂xΦ(x,t)dt+Φ(x,δ2)dδ2dx-Φ(x,δ1)dδ1dx,
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 x0≤x≤x1, then we can deduce that
(3.10)dmdtm[∫abG(t,s)z(s)ds]=∫ab∂m∂tmG(t,s)z(s)ds.
Then we have
(3.11)Dαy(x)=1Γ(m-α)∫t=at=x(x-t)m-α-1[∫s=as=b∂m∂tmG(t,s)z(s)ds]dt.
Changing the order of integration leads to
(3.12)Dαy(x)=1Γ(m-α)∫s=as=b[∫t=at=x(x-t)m-α-1∂m∂tmG(t,s)z(s)dt]ds,Dαy(x)=∫s=as=b[1Γ(m-α)∫t=0t=x(x-t)m-α-1∂m∂tmG(t,s)dt]z(s)ds,Dαy(x)=∫ab(DαG(x,s))z(s)ds=DαGz(x),
and this proves the lemma.Substituting from (3.2)–(3.4) and (3.6) into (1.3) leads to
(3.13)z(x)+θDαGz(x)+βGz(x)=f(x).
We will introduce the operator Ky(x) defined by:
(3.14)Ky(x)=θ∫ab(DαG(x,s))y(s)ds+β∫abG(x,s)y(s)ds,
which maps C2[a,b] to C[a,b]. We also introduce a linear projection PΔ that maps C2[a,b] to SΔ1 piecewise linear interpolation at the grid points {xi}0n. Then (3.13) can be rewritten as:
(3.15)z(x)+Kz(x)=f(x),
and we have also:
(3.16)zΔ(x)+KzΔ(x)=f(x).By the definition ofPΔ [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.Theorem 3.2 (see [34]).
If there isN0 large enough, then {(I+PΔK)-1:n≥N0} exists and consists of a sequence of bounded linear operators. Which means, for a constant δ independent of N0 and z∈C[a,b],  if n≥N0, then ∥(I+PΔK)-1z∥≤δ∥z∥.Theorem 3.3.
Assuming that(H1)
the BVP (1.3) along with boundary conditions (3.1) has a unique solution in C2[a,b],(H2)
the BVPy′′(x)=0 along with boundary conditions (3.1) has a unique solution,then, for somen≥N0 one has
(3.17)‖y-yΔ‖∞≤ck‖y(k+2)‖hk,∀y∈Ck+2[a,b],1≤k≤2,‖y-yΔ‖∞≤c0ψ(y′′,h),∀y∈C2[a,b],where ck is a constant and independent of y, h and ψ(y′′,h)=sup{|y′′(τ+h~)-y′′|:τ,τ+h~∈[a,b],h~≤h}.Proof.
Letz(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
(3.18)PΔz(x)+PΔKz(x)=PΔf(x).
Adding z(x) to both sides of (3.18) and subtracting (3.16) from the results lead to
(3.19)(I+PΔK)(z(x)-zΔ(x))=z(x)-PΔz(x).
Operating on both sides of (3.19) by (I+PΔK)-1 leads to
(3.20)z(x)-zΔ(x)=(I+PΔK)-1(z(x)-PΔz(x)).
Operating on both sides of (3.20) by the operator G and using (3.2)–(3.4), we get
(3.21)y(x)-yΔ(x)=G(I+PΔK)-1(y′′(x)-PΔy′′(x)).
Since the operator G is bounded and from Theorem 3.2 the operator (I+PΔK)-1 is also bounded, then
(3.22)‖y(x)-yΔ(x)‖≤‖G‖‖(I+PΔK)-1‖‖(y′′(x)-PΔy′′(x))‖.
From [33], we have that
(3.23)‖(y′′(x)-PΔy′′(x))‖≤ck‖y(k+2)‖hk,∀y∈Ck+2[a,b],1≤k≤2,(3.24)‖y-yΔ‖∞≤c0ψ(y′′,h),∀y∈C2[a,b],
where, ψ(y′′,h)=sup{|y′′(τ+h~)-y′′|:τ,τ+h~∈[a,b],h~≤h}.
Substituting from Theorem3.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:(4.1)y′′(x)+kπD1.5y(x)+y(x)=0,y(0)=1,y′(0)=0.
The analytical solution of (4.1), as found in [17], has the following form:
(4.2)y(x)=1-∑j=0∞∑r=0∞(-1)r(-kπ)j(j+r)!x2+2r+j/2j!r!(2+2r+j/2)Γ(2+2r+j/2).
As k→0, we may verify that the solution reduces to y(x)=1-∑r=0∞(((-1)rx2+2r)/(2+2r))=cosx.
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.
Figures3 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].Table 1
Numerical results of Example4.1.
x
k= 1
k = 1/5
k = 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.541945Figure 2
Numerical solutions of Example4.1 for (a) k=1, (b) k=0.2, and  (c) k=0.005.Figure 3
A comparison between the analytical solution and our approximated solution fork=1/5.Figure 4
A comparison between the analytical solution and our approximated solution fork=0.005.Example 4.2.
Consider the initial value problem:(4.3)y′′(x)+D0.5y(x)+y(x)=8,y(0)=y′(0)=0.
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 Example4.2.
x
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.499200Example 4.3.
Consider the boundary value problem:(4.4)y′′(x)+0.5D0.5y(x)+y(x)=2+x(1Γ(2.5)x0.5+x),y(0)=0,y(1)=1.
The exact of solution (4.4) is
(4.5)y(x)=x2.
The numerical solutions using shooting method with z1=1 and z2=0.5 led to next guess of initial condition z3=-0.03914, and the results are represented in Table 3.Table 3
Numerical results of Example4.3.
x
Exact solution
Approx. withz3
Error
0
0
0
0
0.125
0.015625
0.012842
2.78
E−3
0.250
0.062500
0.053552
8.95
E−3
0.375
0.140630
0.138569
2.06
E−3
0.500
0.250000
0.246192
3.81
E−3
0.625
0.390630
0.386836
3.79
E−3
0.750
0.562500
0.553883
8.62
E−3
0.875
0.765630
0.758442
7.19
E−3
1
1
0.999971
2.90
E−5Example 4.4.
Consider the boundary value problem:(4.6)y′′(x)+θD0.3y(x)+βy(x)=-12x2+x3[20+θ(120Γ(5.7)x1.7-24Γ(4.7)x0.7)+β(x2-x)],y(0)=y(1)=0.
The exact of solution (4.6) is
(4.7)y(x)=x4(x-1).
The numerical solutions use shooting method for θ=0.5,β=1 and with initial guesses z1=1 and z2=0.5,z3=-0.01224; these lead to z4=-5.07E-4 as the fourth guess for the initial condition and the results are represented in Table 4.Table 4
Numerical results of Example4.4.
x
Exact solution
Approx. solution
Error
0.125
−0.00021
−0.00039
1.73
E−4
0.250
−0.00292
−0.00346
5.35
E−4
0.375
−0.01235
−0.01316
7.98
E−4
0.500
−0.03125
−0.03192
6.74
E−4
0.625
−0.05722
−0.05713
-
9.50
E−5
0.750
−0.07910
−0.07732
-
1.78
E−3
0.875
−0.07327
−0.06985
-
3.42
E−3
1
0
9.44
E−4
-
9.44
E−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.



---

*Source: 638026-2012-08-15.xml*