- Email: [email protected]
A collocation-shooting method for solving fractional boundary value problems
Commun Nonlinear Sci Numer Simulat 15 (2010) 3814–3822
Contents lists available at ScienceDirect
Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns
A collocation-shooting method for solving fractional boundary value problems Qasem M. Al-Mdallal *, Muhammed I. Syam, M.N. Anwar Department of Mathematical Sciences, College of Science, P.O. Box 17551, UAE University, Al-Ain, United Arab Emirates
a r t i c l e
i n f o
Article history: Received 28 April 2009 Received in revised form 5 October 2009 Accepted 18 January 2010 Available online 1 February 2010 AMS subject classifications: 65
a b s t r a c t In this paper, we discuss the numerical solution of special class of fractional boundary value problems of order 2. The method of solution is based on a conjugating collocation and spline analysis combined with shooting method. A theoretical analysis about the existence and uniqueness of exact solution for the present class is proven. Two examples involving Bagley–Torvik equation subject to boundary conditions are also presented; numerical results illustrate the accuracy of the present scheme. Ó 2010 Elsevier B.V. All rights reserved.
Keywords: Collocation method Shooting method Bagley–Torvik equation
1. Introduction Many mathematical models of numerous engineering and physical phenomena involve either ordinary or partial differential equation of fractional order. Researchers have indicated that derivatives of fractional order have the intrinsic property of replicating the history of underlying functions in a weighted form [9]. Further, such property leads to mathematical models of problems governed by differential equations of fractional order. Among such problems different phenomena in areas such as: electromagnetic, acoustics, viscoelasticity, electroanalytical chemistry, neuron modeling and material sciences. For a more detailed discussion of relevant issues the references [6,8,23,26,29,32] give a comprehensive account just to name a few. Furthermore, partial differential equations of fractional orders have been used to describe interesting phenomena such as diffusion problems with certain anomalous features especially for dispersive transport in amorphous semiconductors, glasses, liquid crystals, proteins, and biosystems [18,23–25], and the references therein. Recently, fractional differential equations have found interesting applications in the area of mathematical biology [1,2] due to the relation of such equations with memory that is inherit in corresponding biological systems. The treatment of the above mentioned models takes different facets. For example, existence and uniqueness of solution have been investigated in [11,19,20,33]. Attempts to obtain analytical, semi-analytical, and numerical solutions received considerable attentions of scientists and practitioners. For example, a generalized differential transform [28], is introduced to treat linear partial differential equations of fractional order, the He’s variational iteration method [14], the Adomian decomposition method [3,10,35] and the references therein. The homotopy perturbation technique is used by [4]. Additionally, algorithm based on quadrature formula is introduced in [12], as well as in [13], along with the implementation of the so-called fractional multistep methods discussed in [21,22]. The use of spline approximation [9] and collocation method [34] has been effectively implemented for the solution of initial value * Corresponding author. E-mail addresses: [email protected] (Q.M. Al-Mdallal), [email protected] (M.I. Syam), [email protected] (M.N. Anwar). 1007-5704/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cnsns.2010.01.020
Q.M. Al-Mdallal et al. / Commun Nonlinear Sci Numer Simulat 15 (2010) 3814–3822
3815
fractional differential equations. Additionally, the cubic spline collocation method to solve two classes of special fractional boundary value problems involve the fractional derivative D1=2 has been used in [17]. In the present paper, we consider the Bagley–Torvik boundary value problems of the form
h
i A0 D2 þ A1 D3=2 þ A2 D0 yðtÞ ¼ f ðtÞ;
t 2 ½0; T;
ð1Þ
subject to
yð0Þ ¼ a0 ;
yðTÞ ¼ a1 ;
ð2Þ a
where A0 ; A1 ; A2 ; a0 and a1 are constants with A0 –0, and y 2 L1 ½0; T. Here, D (a is not integer) denotes the fractional differential operator of order a and given by
Da yðxÞ ¼
1 Cðk aÞ
Z
x
ðx tÞka1 yðkÞ ðtÞdt
ð3Þ
0
where k 2 N and satisfies the relation k 1 < a < k. Note that similar study can be done for a general case when the boundary conditions are
c0 yð0Þ þ b0 y0 ð0Þ ¼ a0 ; c1 yðTÞ þ b1 y0 ðTÞ ¼ a1 :
ð4Þ
where c0 ; c1 ; b0 ; b1 ; a0 and a1 are constants. In this paper a technique based on combination between collocation and spline method along with the shooting method is proposed for the solution of problem (1) subject to the boundary conditions (2). It should be noted that the solution of problem (1) subject to specific initial conditions has been discussed by few researchers, see for example [31] and the references therein, whereas, to the best of our knowledge, no work has been done to discuss this problem subject to boundary conditions. This paper is organized as follows: In Section 2, The existence and uniqueness of the exact solution to the present problem are proved. The numerical methods of solution and problem statement are presented in Section 3. Section 4 gives an account of the numerical experiments to be carried out. Section 5 provides the numerical results obtained, concluded by discussion and conclusion remarks. 2. Existence and uniqueness The existence and uniqueness of the exact solution to problem (1) subject to the boundary conditions (2) are discussed herein. Prior to presenting the main result, we recall the following theorem. Theorem 1. The Laplace transform for Caputo’s fractional derivative [30] is given by
£½yðpÞ ðsÞ ¼ sp YðsÞ
n1 X
spk1 yðkÞ ð0Þ;
n 1 < p 6 n:
ð5Þ
k¼0
Our main claim is introduced in the following theorem: Theorem 2. Problem (1) has a unique solution given by
yðtÞ ¼ yh ðtÞ þ cyp ðtÞ where c ¼ a1 yh ðTÞ=yp ðTÞ;
yh ðtÞ ¼ £1
a1 s1=2 þ A1 a1
s5=2 þ A1 s2 þ A2 s
ðtÞ 1=2
is the solution of the homogenous solution to ð3=2Þ
A0 y00h ðtÞ þ A1 yh
ðtÞ þ A2 yh ðtÞ ¼ 0;
0 < t < T;
y0h ð0Þ ¼ a1 :
yh ð0Þ ¼ 0;
ð6Þ
and
yp ðtÞ ¼ £1
FðsÞ þ A0 a0 s þ A1 s1=2 a0 ðtÞ A0 s2 þ A1 s3=2 þ A2
is the particular solution to
A0 y00p ðtÞ þ A1 ypð3=2Þ ðtÞ þ A2 yp ðtÞ ¼ f ðtÞ; yp ð0Þ ¼ a0 ;
y0p ð0Þ
¼ 0:
0 < t < T;
ð7Þ
3816
Q.M. Al-Mdallal et al. / Commun Nonlinear Sci Numer Simulat 15 (2010) 3814–3822
Proof. The proof of this theorem consists of three major steps following similar technique to that presented in [32]. In the first step, we solve a homogenous initial value problem (6) while in the second one we solve a nonhomogeneous initial value problem (7). In the last step, we combine the two solutions to determine the exact solution. The general solutions to problems (6) and (7) can be determined by taking the Laplace transform of both sides and then applying Theorem 1; those solutions are respectively found to be
Y h ðsÞ ¼
A0 a1 þ A1 s1=2 a1 A0 a1 s1=2 þ A1 a1 ¼ A0 s2 þ A1 s3=2 þ A2 A0 s5=2 þ A1 s2 þ A2 s1=2
ð8Þ
Y p ðsÞ ¼
FðsÞ þ A0 a0 s þ A1 s1=2 a0 A0 s2 þ A1 s3=2 þ A2
ð9Þ
and
where Y h ðsÞ ¼ £½yh ðtÞ and Y p ðsÞ ¼ £½yp ðtÞ. Applying the properties of the inverse Laplace transform described in [27] on (8) and (9), we determine the homogenous and particular solutions yh ðtÞ and yp ðtÞ, respectively. Let the general solution has the form
yðtÞ ¼ d1 yh ðtÞ þ d2 yp ðtÞ
ð10Þ
where d1 and d2 are constants need to be determined. It can be easily verified that yðtÞ in (10) is a solution for problem (1). On the other hand, we need to select the constants d1 and d2 so that the solution satisfy the boundary conditions (2); i.e.
a0 ¼ yð0Þ ¼ d1 yh ð0Þ þ d2 yp ð0Þ ¼ d1 a0 a1 ¼ yð1Þ ¼ d1 yh ðTÞ þ d2 yp ðTÞ: Hence, it follows that:
d1 ¼ 1 and d2 ¼
a1 yh ðTÞ yp ðTÞ
ð11Þ
:
By combining the results of (10) and (11).
h
3. Method of solution The following is a brief derivation of the algorithm used to solve problem (1) subject to (2). The method of solution presented herein is based on conjugating collocation approach combined with multiple shooting method. 3.1. Collocation method For sake of simplicity we discuss the solution of (1) as initial value problem with
yð0Þ ¼ a0 ;
^1; y0 ð0Þ ¼ a
ð12Þ
^ 1 is unknown constant which will be determined later. where a We firstly subdivide the interval ½0; T into N uniform subintervals rn ¼ ½t n ; t nþ1 where n ¼ 0; . . . ; N 1. Let Z N ¼ ftn ¼ nh : n ¼ 1; . . . ; N 1g with h ¼ T=N. Assume that the exact solution of Eq. (1) subject to condition (12) can be approximated by an element u in the space of piecewise polynomials of degree m þ d which are d-times continuously differentiable on I represented by
Smþd ðZ N Þ ¼ fu 2 C d ðIÞ : ujrn is a polynomial of degree m þ dg; ðdÞ
for known integers m P 1 and d P 0. It should be noted that the integer m represents the number of collocation points in each subinterval rn ; ðn ¼ 0; . . . ; N 1Þ; those points are defined as
X n ¼ ft n;i ¼ t n þ ci h : i ¼ 1; . . . ; mg;
with 0 6 c1 6 6 cm 6 1;
whereas, the integer d represents the number of the initial conditions minus one (i.e. d ¼ 1 in the present work). ConseðdÞ quently, the exact solution, y, of (1) and (12) needs to be approximated on I by an element u 2 Smþd ðZ N Þ satisfying
h i A0 D2 þ A1 D3=2 þ A2 D0 uðtÞ ¼ f ðtÞ;
t2
N1 [
Xn;
ð13Þ
n¼0
subject to
uð0Þ ¼ a0 ;
^1; u0 ð0Þ ¼ a
On each subinterval,
rn , the spline u can be presented as a piecewise polynomials of degree m þ d of the form
ð14Þ
Q.M. Al-Mdallal et al. / Commun Nonlinear Sci Numer Simulat 15 (2010) 3814–3822
uðtÞ ¼ un ðt n þ shÞ ¼
d X
s aðnÞ s s þ
s¼0
m X
ðnÞ dþr
br
s ; t 2 rn ;
3817
ð15Þ
r¼1
where s 2 ½0; 1. The fractional differential operator of order q for the collocation solution (15) at t ¼ tn þ ci h is found by Blank [9] and given by:
" # q n X d n X m X X h ðnj;qÞ ðjÞ ðnj;qÞ ðjÞ D ðun ðtn þ ci hÞÞ ¼ w as þ wi;rþd br ; Cð1 qÞ j¼0 s¼0 i;s j¼0 r¼1 q
ð16Þ
where q 2 Rþ ; d 2 N and, in general, ðj;qÞ wi;k
8 q q k¼0 > < ðj þ ci Þ dj;0 ðj þ ci 1Þ ; k k ¼ P qþk Q p > dj;0 ðj þ ci 1Þmq Sqm;k ; k P 1; : ðj þ ci Þ pq m¼0
p¼1
where
Sqm;k
¼
Qm
kmþp p¼1 pq
and
dj;0
¼ 0 if j ¼ 0, and 1 otherwise. Applying Blank’s result (16) on (13); one obtains
" # d m X X 2 s2 ðnÞ dþr2 ðnÞ h sðs 1Þci as þ ðd þ rÞðd þ r 1Þci br þ A1 s¼2
þ A2
" d X
r¼1
csi aðnÞ s þ
s¼0
¼ f ðtn;i Þ;
m X
#
" # 3=2 n X d n X m X X h ðnj;3=2Þ ðjÞ ðnj;3=2Þ ðjÞ w as þ wi;rþd br Cð1=2Þ j¼0 s¼0 i;s j¼0 r¼1
ðnÞ
cdþr br i
r¼1
ðfor i ¼ 1; . . . ; mÞ:
It can be shown that the above equations can be expressed in the following m m matrix form ðnÞ
Vb
¼ UaðnÞ þ F
where 1=2
ðVÞi;r ¼ ðd þ rÞðd þ r 1Þcdþr2 þ A1 i
h 2 ð0;3=2Þ þ A2 h cdþr ; w i Cð1=2Þ i;rþd
ði; r ¼ 1; . . . ; mÞ;
1=2
ðUÞi;s ¼ sðs 1Þcis2 þ A1
2
ðFÞi;1 ¼ h f ðt n;i Þ
h 2 ð0;3=2Þ þ A2 h csi ; w Cð1=2Þ i;s
" #! 1=2 n1 X d n1 X m X X h ðnj;3=2Þ ðjÞ ðnj;3=2Þ ðjÞ ; A1 w as þ wi;rþd br Cð1=2Þ j¼0 s¼0 i;s j¼0 r¼1
dn;0
for i ¼ 1; . . . ; m. Here aðnÞ ¼ ½a0 ; . . . ; ad t and b
ðnÞ
ðnÞ
ðnÞ
¼ ½b1 ; . . . ; bm t , where ½t means the transpose of the vector. Note that hs s it is known from the initial conditions, i.e. að0Þ ¼ hs! ddtys ð0Þ ; s ¼ 0; . . . ; d. On the other hand, when ðnÞ
when n ¼ 0, the vector að0Þ
ði ¼ 1; . . . ; m; s ¼ 0; . . . dÞ:
ðnÞ
n P 1, we have to use the smoothness conditions at t ¼ tn which leads to a relation between the unknown vector aðnþ1Þ and ðnÞ
the known vectors aðnÞ and b . In the present work, this relation is given by
aðnþ1Þ ¼ HaðnÞ þ Gb
ðnÞ
where
H¼
1 1 0 1
;
and G ¼
1 1 1 2 3 4
:
In the next subsection, we discuss the implementation of the multiple shooting method in the numerical solution of problem (1) and (2). 3.2. Shooting method For simplicity, we can rewrite problem (1) and (2) in the form
L½y ¼ f ðtÞ;
t 2 ½0; T
ð17Þ
yðTÞ ¼ a1 ;
ð18Þ
subject to
yð0Þ ¼ a0 ;
3818
Q.M. Al-Mdallal et al. / Commun Nonlinear Sci Numer Simulat 15 (2010) 3814–3822
where
L ¼ A0 D2 þ A1 D3=2 þ A2 D0 Firstly we partition the domain ½0; T as
0 ¼ T 0 < T 1 < < T L ¼ T: Hence, theoretically, the solution of problem (17) and (18) on the time interval ½0; T can be determined by solving problem (17) and (18) on the subintervals ½T l ; T lþ1 , for ðl ¼ 0; . . . ; L 1Þ. For the purpose of implementing shooting method, we added extra conditions at the grid points T l , ðl ¼ 0; . . . ; L 1Þ, so that we have the following set of initial value problems
L½xl ¼ f ðtÞ;
t 2 ½T l ; T lþ1
ð19Þ
subject to
xl ðT l Þ ¼ a^ 2l ; x0l ðT l Þ ¼ a^ 2lþ1
ð20Þ
^ j ; ðj ¼ 1; . . . ; 2L 1Þ are unknown real parameters and a ^ 0 ¼ a0 . The solution of initial value problems (19) and (20) where a ^ l ; ðl ¼ 1; . . . ; 2L 1Þ, will be will be obtained by the method described in the previous subsection where the parameters a determined by solving the following system of algebraic equations
x0 ðT 1 ; a^ 1 Þ ¼ x1 ðT 1 ; a^ 2 ; a^ 3 Þ; x00 ðT 1 ; a^ 1 Þ ¼ x01 ðT 1 ; a^ 2 ; a^ 3 Þ; x1 ðT 2 ; a^ 2 ; a^ 3 Þ ¼ x2 ðT 2 ; a^ 4 ; a^ 5 Þ; x01 ðT 2 ; a^ 2 ; a^ 3 Þ ¼ x02 ðT 2 ; a^ 4 ; a^ 5 Þ; .. .
xk2 ðT k1 ; a^ 2k4 ; a^ 2k3 Þ ¼ xk1 ðT k1 ; a^ 2k2 ; a^ 2k1 Þ; x0k2 ðT k1 ; a^ 2k4 ; a^ 2k3 Þ ¼ x0k1 ðT k1 ; a^ 2k2 ; a^ 2k1 Þ; xk1 ðT; a^ 2k2 ; a^ 2k1 Þ ¼ a1 which can be solved numerically using several numerical techniques such as Newton’s method. In summary, this technique is called multiple shooting technique of order L. Note that, we use shooting method of low order ðL ¼ 4Þ in our computations which are done using Matlab. For more details about the multiple shooting method see [5,15,16].
4. Numerical results In this section we will consider the Bagley–Torvik equation which arises in the modeling of the motion of a rigid plate immersed in Newtonian fluid, see Bagley and Torvik [7]. Example 1. Consider the Bagley–Torvik equation
1 1 D2 yðtÞ þ D3=2 yðtÞ þ yðtÞ ¼ 8v ðtÞ 8v ðt 1Þ; 2 2
t 2 ½0; 20;
subject to
yð0Þ ¼ 0;
yð20Þ ¼ 1:48433;
where v ðtÞ is the Heaviside function. Applying the multiple shooting method of order four requires dividing the domain ½0; 20 into four subintervals as follows
0 ¼ T 0 < T 1 ¼ 5 < T 2 ¼ 10 < T 3 ¼ 15 < T 4 ¼ 20: Consequently, our problem is divided into the following initial value problems
L½x1 ¼ v ðtÞ; L½x2 ¼ v ðtÞ; L½x3 ¼ v ðtÞ;
^1; 0 < t < 5; x1 ð0; :Þ ¼ 0; x01 ð0; :Þ ¼ a ^ 2 ; x02 ð5; :Þ ¼ a ^3; 5 < t < 10; x2 ð5; :Þ ¼ a ^ 4 ; x03 ð10; :Þ ¼ a ^5; 10 < t < 15; x3 ð10; :Þ ¼ a
L½x4 ¼ v ðtÞ;
15 < t < 20;
x4 ð15; :Þ ¼ a^ 6 ; x04 ð15; :Þ ¼ a^ 7 :
where
1 1 L½xl ¼ D2 xl ðtÞ þ D3=2 xl ðtÞ þ xl ðtÞ: 2 2 and
ð21Þ
Q.M. Al-Mdallal et al. / Commun Nonlinear Sci Numer Simulat 15 (2010) 3814–3822
3819
vðtÞ ¼ v ðtÞ 8v ðt 1Þ: ^ 1 Þ; x2 ðt; :Þ ¼ x2 ðt; a ^2 ; a ^ 3 Þ; x3 ðt; :Þ ¼ x3 ðt; a ^4; a ^ 5 Þ; and x4 ðt; :Þ ¼ x4 ðt; a ^6; a ^ 7 Þ. Each initial value probHere x1 ðt; :Þ ¼ x1 ðt; a lem of (21) will be solved using the collocation method where we use m ¼ 3 with the collocation parameters
2
c1
2x
3
1 þ1
6 7 4 c2 5 ¼ c3
3
6 x 2þ1 7 6 2 7 4 2 5 x3 þ1 2
where xk ¼ cosð2kþ1 pÞ; ðk ¼ 1; 2; 3Þ are the roots of the Tchebychev polynomial of degree 3. The choice of m ¼ 3 is due to the 6 fact that although the more collocation points are used the more accurate the resulting numerical solution is to be expected, however, this will cause round off errors to have serious effects on the computed results. Needless to say that it is difficult to ^ i ; i ¼ 1; 2; . . . ; 7 are determined by solving the find a specific optimum choice for m. Finally, the values for the parameters a following algebraic system
x1 ð5; a^ 1 Þ ¼ x2 ð5; a^ 2 ; a^ 3 Þ; x01 ð5; a^ 1 Þ ¼ x02 ð5; a^ 2 ; a^ 3 Þ; x2 ð10; a^ 2 ; a^ 3 Þ ¼ x3 ð10; a^ 4 ; a^ 5 Þ; x02 ð10; a^ 2 ; a^ 3 Þ ¼ x03 ð10; a^ 4 ; a^ 5 Þ; x3 ð15; a^ 4 ; a^ 5 Þ ¼ x4 ð15; a^ 6 ; a^ 7 Þ; x03 ð15; a^ 4 ; a^ 5 Þ ¼ x04 ð15; a^ 6 ; a^ 7 Þ; x4 ð20; a^ 6 ; a^ 7 Þ ¼ 1:48433: Using Matlab package we obtain
a^ 1 ¼ 0; a^ 2 ¼ 2:9439356637512617; a^ 3 ¼ 3:50426674344817; a^ 4 ¼ 2:848380869392841; a^ 5 ¼ 1:858241791552473; a^ 6 ¼ 2:168701764974635; a^ 7 ¼ 0:8360760617165397: The graph of the approximated solution is given in Fig. (1). Note that this example was solved theoretically by Podlubny [30] where the compact form of the solution is given by
yðtÞ ¼
Z
t
8Gðt sÞðv ðsÞ v ðs 1ÞÞds
0
where G is the fractional Green’s function defined as
GðtÞ ¼
1 1 X ð1=2Þk X ðj þ kÞ!ðt1=2 Þj : j 1 1 3 k! j¼0 j!2 Cð2 j þ 2 k þ 2 k þ 2Þ k¼0
Here C is the Gamma function. The error between the exact solution and computed one is shown in Fig. 2. It is clearly seen that the two solutions are in excellent agreement. The computed L2 error norm is given by
kuðtÞ yðtÞk ¼
Z
20
ðyðtÞ uðtÞÞ2 dt ¼ 1:7 1010 :
0
Fig. 1. The approximate solution for Example 1.
3820
Q.M. Al-Mdallal et al. / Commun Nonlinear Sci Numer Simulat 15 (2010) 3814–3822
Fig. 2. Computed absolute error between the present numerical solution and the theoretical one by Podlubny [30] for Example 1.
Example 2. Consider the Bagley–Torvik equation
D2 yðtÞ þ D3=2 yðtÞ þ yðtÞ ¼ f ðtÞ;
t 2 ½0; 5;
ð22Þ
subject to
yð0Þ ¼ 0;
yð5Þ ¼ 25
2
where f ðtÞ ¼ t þ 2 þ 4
qffiffiffi t
p
and the exact solution yðtÞ ¼ t 2 .
Using the multiple shooting method of order five then following the same steps in the previous example, the parameters
a^ i ; i ¼ 1; 2; . . . ; 9 are found to be
a^ 1 ¼ 0; a^ 2 ¼ 0:99999999999999; a^ 3 ¼ 1:99999999999999; a^ 4 ¼ 3:99999999999999; a^ 5 ¼ 4; a^ 6 ¼ 8:99999999999999; a^ 7 ¼ 6:0000000000000001; a^ 8 ¼ 16; a^ 9 ¼ 7:99999999999999: The graph of the approximated solution is given in Fig. 3. In addition, Fig. 4 shows the error between the computed and the exact solutions which indicates the efficiency of the present numerical technique. The computed L2 error norm is given by
kuðtÞ yðtÞk ¼
Z
20
ðyðtÞ uðtÞÞ2 dt ¼ 3:78 1012 :
0
Fig. 3. The approximate solution for Example 2.
3821
Q.M. Al-Mdallal et al. / Commun Nonlinear Sci Numer Simulat 15 (2010) 3814–3822
Fig. 4. Computed absolute error between the present numerical solution and the theoretical one for Example 2.
Table 1 Effect of the initial condition, yð0Þ, on the computed L2 error norm. yð0Þ
kuðtÞ yðtÞk (Example 1)
kuðtÞ yðtÞk (Example 2)
0.000
10
1:71 10
3:78 1012
0.001
2:10 1010
3:61 1012
0.001
3:84 109
2:94 1012
0.100
7:14 1010
7:22 1012
0.500
2:92 1010
1:10 1012
1.000
8:73 1010
2:55 1012
In order to demonstrate the effect of the initial condition on the numerical results, Table 1 has been obtained for the above examples. It can be concluded that the effect of the initial condition is not significant on the error. 5. Conclusion In this paper, we have solved special class of fractional linear differential equation of order 2 subject to boundary conditions. A well-known example involving this class is the Bagley–Torvik equation. The method of solution is based on conjugating collocation and spline technique with shooting method. The numerical results for the examples demonstrate the efficiency and accuracy of the present method. The proposed method has the advantage of being both effective and flexible for solving boundary value problems governed by ordinary differential equations of fractional orders. The method combines the efficacy of multiple shooting method along with the flexibility of the collocation approach. To the best of the authors’ knowledge such combination can be considered as a new contribution. Acknowledgement The authors would like to express their appreciation for the valuable comments of the reviewers. The authors also would like to express their sincere appreciation to the United Arab Emirates University Research Affairs for the financial support of Grant nos. 04 13 02 11=07; 01 13 02 11=07. References [1] Ahmed E, Elgazzar AS. On fractional order differential equations model for nonlocal epidemics. Phys A 2007;379(2):607–14. [2] Ahmed E, El-Sayed AMA, Elsaka HAA. Equilibrium points, stability, and numerical solutions of fractional-order predator-prey and rabies models. J Math Anal Appl 2007;325(1):542–53. [3] Al-Mdallal QM. An efficient method for solving fractional Sturm-Liouville problems. Chaos, Solitons Fractals 2009;40(1):183–9. [4] Al-Mdallal QM. On the numerical solution of fractional Sturm-Liouville problems. Int J Comput Math 2009;99999:1. [5] Attili B, Syam M. Efficient shooting method for solving two point boundary value problems. Chaos, Solitons Fractals 2008;35(5):895–903. [6] Babenko Y. Non integer differential equation in engineering, Conference Bordeaux 3–8 July 1994. [7] Bagley RL, Torvik PJ. On the appearance of the fractional derivative in the behavior of real materials. ASME J Appl Mech 1984;51(2):294–8. [8] Beyer H, Kempfle S. Definition of physically consistent damping laws with fractional derivatives. Zeitschrift fr Angewandte Mathematik und Mechanik 1995;75(8):623–35.
3822
Q.M. Al-Mdallal et al. / Commun Nonlinear Sci Numer Simulat 15 (2010) 3814–3822
[9] Blank L. Numerical treatment of differential equations of fractional order. Numerical Analysis Report 287, Manchester Center for Numerical Computational Mathematics; 1996. [10] Daftardar-Geiji V, Jafari H. Adomian decomposition: a tool for solving a system of fractional differential equations. J Math Anal Appl 2005;301(2):508–18. [11] Daftardar-Geiji V, Jafari H. Analysis of a system of nonautonomous fractional differential equations involving Capato derivatives. J Math Anal Appl 2007;328:1026–33. [12] Diethelm K. An algorithm for the numerical solution of differential equations of fractional order. Elect Translat Numer Anal 1997;5:1–6. [13] Diethelm K, Ford NJ, Freed AD. A predictor–corrector approach for the numerical solution of fractional differential equations. Nonlinear Dynam 2002;29:3–22. [14] Ghorbani A, Alavi A. Applications of He’s variational iteration method to solve semidifferential equations of nth order. Math Prob Eng 2008;2008:1–9. [15] Ha SN. A nonlinear shooting method for two point boundary value problems. Comput Math Appl 2001;42:1411–20. [16] Keller HB. Numerical methods for two point boundary value problems. New York: Dover; 1992. [17] Lia M, Jimenezc S, Niea N, Tanga Y,Vazqueze L. Solving Two-point boundary value problems of fractional differential equations by spline collocation methods. Available from
Contents lists available at ScienceDirect
Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns
A collocation-shooting method for solving fractional boundary value problems Qasem M. Al-Mdallal *, Muhammed I. Syam, M.N. Anwar Department of Mathematical Sciences, College of Science, P.O. Box 17551, UAE University, Al-Ain, United Arab Emirates
a r t i c l e
i n f o
Article history: Received 28 April 2009 Received in revised form 5 October 2009 Accepted 18 January 2010 Available online 1 February 2010 AMS subject classifications: 65
a b s t r a c t In this paper, we discuss the numerical solution of special class of fractional boundary value problems of order 2. The method of solution is based on a conjugating collocation and spline analysis combined with shooting method. A theoretical analysis about the existence and uniqueness of exact solution for the present class is proven. Two examples involving Bagley–Torvik equation subject to boundary conditions are also presented; numerical results illustrate the accuracy of the present scheme. Ó 2010 Elsevier B.V. All rights reserved.
Keywords: Collocation method Shooting method Bagley–Torvik equation
1. Introduction Many mathematical models of numerous engineering and physical phenomena involve either ordinary or partial differential equation of fractional order. Researchers have indicated that derivatives of fractional order have the intrinsic property of replicating the history of underlying functions in a weighted form [9]. Further, such property leads to mathematical models of problems governed by differential equations of fractional order. Among such problems different phenomena in areas such as: electromagnetic, acoustics, viscoelasticity, electroanalytical chemistry, neuron modeling and material sciences. For a more detailed discussion of relevant issues the references [6,8,23,26,29,32] give a comprehensive account just to name a few. Furthermore, partial differential equations of fractional orders have been used to describe interesting phenomena such as diffusion problems with certain anomalous features especially for dispersive transport in amorphous semiconductors, glasses, liquid crystals, proteins, and biosystems [18,23–25], and the references therein. Recently, fractional differential equations have found interesting applications in the area of mathematical biology [1,2] due to the relation of such equations with memory that is inherit in corresponding biological systems. The treatment of the above mentioned models takes different facets. For example, existence and uniqueness of solution have been investigated in [11,19,20,33]. Attempts to obtain analytical, semi-analytical, and numerical solutions received considerable attentions of scientists and practitioners. For example, a generalized differential transform [28], is introduced to treat linear partial differential equations of fractional order, the He’s variational iteration method [14], the Adomian decomposition method [3,10,35] and the references therein. The homotopy perturbation technique is used by [4]. Additionally, algorithm based on quadrature formula is introduced in [12], as well as in [13], along with the implementation of the so-called fractional multistep methods discussed in [21,22]. The use of spline approximation [9] and collocation method [34] has been effectively implemented for the solution of initial value * Corresponding author. E-mail addresses: [email protected] (Q.M. Al-Mdallal), [email protected] (M.I. Syam), [email protected] (M.N. Anwar). 1007-5704/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cnsns.2010.01.020
Q.M. Al-Mdallal et al. / Commun Nonlinear Sci Numer Simulat 15 (2010) 3814–3822
3815
fractional differential equations. Additionally, the cubic spline collocation method to solve two classes of special fractional boundary value problems involve the fractional derivative D1=2 has been used in [17]. In the present paper, we consider the Bagley–Torvik boundary value problems of the form
h
i A0 D2 þ A1 D3=2 þ A2 D0 yðtÞ ¼ f ðtÞ;
t 2 ½0; T;
ð1Þ
subject to
yð0Þ ¼ a0 ;
yðTÞ ¼ a1 ;
ð2Þ a
where A0 ; A1 ; A2 ; a0 and a1 are constants with A0 –0, and y 2 L1 ½0; T. Here, D (a is not integer) denotes the fractional differential operator of order a and given by
Da yðxÞ ¼
1 Cðk aÞ
Z
x
ðx tÞka1 yðkÞ ðtÞdt
ð3Þ
0
where k 2 N and satisfies the relation k 1 < a < k. Note that similar study can be done for a general case when the boundary conditions are
c0 yð0Þ þ b0 y0 ð0Þ ¼ a0 ; c1 yðTÞ þ b1 y0 ðTÞ ¼ a1 :
ð4Þ
where c0 ; c1 ; b0 ; b1 ; a0 and a1 are constants. In this paper a technique based on combination between collocation and spline method along with the shooting method is proposed for the solution of problem (1) subject to the boundary conditions (2). It should be noted that the solution of problem (1) subject to specific initial conditions has been discussed by few researchers, see for example [31] and the references therein, whereas, to the best of our knowledge, no work has been done to discuss this problem subject to boundary conditions. This paper is organized as follows: In Section 2, The existence and uniqueness of the exact solution to the present problem are proved. The numerical methods of solution and problem statement are presented in Section 3. Section 4 gives an account of the numerical experiments to be carried out. Section 5 provides the numerical results obtained, concluded by discussion and conclusion remarks. 2. Existence and uniqueness The existence and uniqueness of the exact solution to problem (1) subject to the boundary conditions (2) are discussed herein. Prior to presenting the main result, we recall the following theorem. Theorem 1. The Laplace transform for Caputo’s fractional derivative [30] is given by
£½yðpÞ ðsÞ ¼ sp YðsÞ
n1 X
spk1 yðkÞ ð0Þ;
n 1 < p 6 n:
ð5Þ
k¼0
Our main claim is introduced in the following theorem: Theorem 2. Problem (1) has a unique solution given by
yðtÞ ¼ yh ðtÞ þ cyp ðtÞ where c ¼ a1 yh ðTÞ=yp ðTÞ;
yh ðtÞ ¼ £1
a1 s1=2 þ A1 a1
s5=2 þ A1 s2 þ A2 s
ðtÞ 1=2
is the solution of the homogenous solution to ð3=2Þ
A0 y00h ðtÞ þ A1 yh
ðtÞ þ A2 yh ðtÞ ¼ 0;
0 < t < T;
y0h ð0Þ ¼ a1 :
yh ð0Þ ¼ 0;
ð6Þ
and
yp ðtÞ ¼ £1
FðsÞ þ A0 a0 s þ A1 s1=2 a0 ðtÞ A0 s2 þ A1 s3=2 þ A2
is the particular solution to
A0 y00p ðtÞ þ A1 ypð3=2Þ ðtÞ þ A2 yp ðtÞ ¼ f ðtÞ; yp ð0Þ ¼ a0 ;
y0p ð0Þ
¼ 0:
0 < t < T;
ð7Þ
3816
Q.M. Al-Mdallal et al. / Commun Nonlinear Sci Numer Simulat 15 (2010) 3814–3822
Proof. The proof of this theorem consists of three major steps following similar technique to that presented in [32]. In the first step, we solve a homogenous initial value problem (6) while in the second one we solve a nonhomogeneous initial value problem (7). In the last step, we combine the two solutions to determine the exact solution. The general solutions to problems (6) and (7) can be determined by taking the Laplace transform of both sides and then applying Theorem 1; those solutions are respectively found to be
Y h ðsÞ ¼
A0 a1 þ A1 s1=2 a1 A0 a1 s1=2 þ A1 a1 ¼ A0 s2 þ A1 s3=2 þ A2 A0 s5=2 þ A1 s2 þ A2 s1=2
ð8Þ
Y p ðsÞ ¼
FðsÞ þ A0 a0 s þ A1 s1=2 a0 A0 s2 þ A1 s3=2 þ A2
ð9Þ
and
where Y h ðsÞ ¼ £½yh ðtÞ and Y p ðsÞ ¼ £½yp ðtÞ. Applying the properties of the inverse Laplace transform described in [27] on (8) and (9), we determine the homogenous and particular solutions yh ðtÞ and yp ðtÞ, respectively. Let the general solution has the form
yðtÞ ¼ d1 yh ðtÞ þ d2 yp ðtÞ
ð10Þ
where d1 and d2 are constants need to be determined. It can be easily verified that yðtÞ in (10) is a solution for problem (1). On the other hand, we need to select the constants d1 and d2 so that the solution satisfy the boundary conditions (2); i.e.
a0 ¼ yð0Þ ¼ d1 yh ð0Þ þ d2 yp ð0Þ ¼ d1 a0 a1 ¼ yð1Þ ¼ d1 yh ðTÞ þ d2 yp ðTÞ: Hence, it follows that:
d1 ¼ 1 and d2 ¼
a1 yh ðTÞ yp ðTÞ
ð11Þ
:
By combining the results of (10) and (11).
h
3. Method of solution The following is a brief derivation of the algorithm used to solve problem (1) subject to (2). The method of solution presented herein is based on conjugating collocation approach combined with multiple shooting method. 3.1. Collocation method For sake of simplicity we discuss the solution of (1) as initial value problem with
yð0Þ ¼ a0 ;
^1; y0 ð0Þ ¼ a
ð12Þ
^ 1 is unknown constant which will be determined later. where a We firstly subdivide the interval ½0; T into N uniform subintervals rn ¼ ½t n ; t nþ1 where n ¼ 0; . . . ; N 1. Let Z N ¼ ftn ¼ nh : n ¼ 1; . . . ; N 1g with h ¼ T=N. Assume that the exact solution of Eq. (1) subject to condition (12) can be approximated by an element u in the space of piecewise polynomials of degree m þ d which are d-times continuously differentiable on I represented by
Smþd ðZ N Þ ¼ fu 2 C d ðIÞ : ujrn is a polynomial of degree m þ dg; ðdÞ
for known integers m P 1 and d P 0. It should be noted that the integer m represents the number of collocation points in each subinterval rn ; ðn ¼ 0; . . . ; N 1Þ; those points are defined as
X n ¼ ft n;i ¼ t n þ ci h : i ¼ 1; . . . ; mg;
with 0 6 c1 6 6 cm 6 1;
whereas, the integer d represents the number of the initial conditions minus one (i.e. d ¼ 1 in the present work). ConseðdÞ quently, the exact solution, y, of (1) and (12) needs to be approximated on I by an element u 2 Smþd ðZ N Þ satisfying
h i A0 D2 þ A1 D3=2 þ A2 D0 uðtÞ ¼ f ðtÞ;
t2
N1 [
Xn;
ð13Þ
n¼0
subject to
uð0Þ ¼ a0 ;
^1; u0 ð0Þ ¼ a
On each subinterval,
rn , the spline u can be presented as a piecewise polynomials of degree m þ d of the form
ð14Þ
Q.M. Al-Mdallal et al. / Commun Nonlinear Sci Numer Simulat 15 (2010) 3814–3822
uðtÞ ¼ un ðt n þ shÞ ¼
d X
s aðnÞ s s þ
s¼0
m X
ðnÞ dþr
br
s ; t 2 rn ;
3817
ð15Þ
r¼1
where s 2 ½0; 1. The fractional differential operator of order q for the collocation solution (15) at t ¼ tn þ ci h is found by Blank [9] and given by:
" # q n X d n X m X X h ðnj;qÞ ðjÞ ðnj;qÞ ðjÞ D ðun ðtn þ ci hÞÞ ¼ w as þ wi;rþd br ; Cð1 qÞ j¼0 s¼0 i;s j¼0 r¼1 q
ð16Þ
where q 2 Rþ ; d 2 N and, in general, ðj;qÞ wi;k
8 q q k¼0 > < ðj þ ci Þ dj;0 ðj þ ci 1Þ ; k k ¼ P qþk Q p > dj;0 ðj þ ci 1Þmq Sqm;k ; k P 1; : ðj þ ci Þ pq m¼0
p¼1
where
Sqm;k
¼
Qm
kmþp p¼1 pq
and
dj;0
¼ 0 if j ¼ 0, and 1 otherwise. Applying Blank’s result (16) on (13); one obtains
" # d m X X 2 s2 ðnÞ dþr2 ðnÞ h sðs 1Þci as þ ðd þ rÞðd þ r 1Þci br þ A1 s¼2
þ A2
" d X
r¼1
csi aðnÞ s þ
s¼0
¼ f ðtn;i Þ;
m X
#
" # 3=2 n X d n X m X X h ðnj;3=2Þ ðjÞ ðnj;3=2Þ ðjÞ w as þ wi;rþd br Cð1=2Þ j¼0 s¼0 i;s j¼0 r¼1
ðnÞ
cdþr br i
r¼1
ðfor i ¼ 1; . . . ; mÞ:
It can be shown that the above equations can be expressed in the following m m matrix form ðnÞ
Vb
¼ UaðnÞ þ F
where 1=2
ðVÞi;r ¼ ðd þ rÞðd þ r 1Þcdþr2 þ A1 i
h 2 ð0;3=2Þ þ A2 h cdþr ; w i Cð1=2Þ i;rþd
ði; r ¼ 1; . . . ; mÞ;
1=2
ðUÞi;s ¼ sðs 1Þcis2 þ A1
2
ðFÞi;1 ¼ h f ðt n;i Þ
h 2 ð0;3=2Þ þ A2 h csi ; w Cð1=2Þ i;s
" #! 1=2 n1 X d n1 X m X X h ðnj;3=2Þ ðjÞ ðnj;3=2Þ ðjÞ ; A1 w as þ wi;rþd br Cð1=2Þ j¼0 s¼0 i;s j¼0 r¼1
dn;0
for i ¼ 1; . . . ; m. Here aðnÞ ¼ ½a0 ; . . . ; ad t and b
ðnÞ
ðnÞ
ðnÞ
¼ ½b1 ; . . . ; bm t , where ½t means the transpose of the vector. Note that hs s it is known from the initial conditions, i.e. að0Þ ¼ hs! ddtys ð0Þ ; s ¼ 0; . . . ; d. On the other hand, when ðnÞ
when n ¼ 0, the vector að0Þ
ði ¼ 1; . . . ; m; s ¼ 0; . . . dÞ:
ðnÞ
n P 1, we have to use the smoothness conditions at t ¼ tn which leads to a relation between the unknown vector aðnþ1Þ and ðnÞ
the known vectors aðnÞ and b . In the present work, this relation is given by
aðnþ1Þ ¼ HaðnÞ þ Gb
ðnÞ
where
H¼
1 1 0 1
;
and G ¼
1 1 1 2 3 4
:
In the next subsection, we discuss the implementation of the multiple shooting method in the numerical solution of problem (1) and (2). 3.2. Shooting method For simplicity, we can rewrite problem (1) and (2) in the form
L½y ¼ f ðtÞ;
t 2 ½0; T
ð17Þ
yðTÞ ¼ a1 ;
ð18Þ
subject to
yð0Þ ¼ a0 ;
3818
Q.M. Al-Mdallal et al. / Commun Nonlinear Sci Numer Simulat 15 (2010) 3814–3822
where
L ¼ A0 D2 þ A1 D3=2 þ A2 D0 Firstly we partition the domain ½0; T as
0 ¼ T 0 < T 1 < < T L ¼ T: Hence, theoretically, the solution of problem (17) and (18) on the time interval ½0; T can be determined by solving problem (17) and (18) on the subintervals ½T l ; T lþ1 , for ðl ¼ 0; . . . ; L 1Þ. For the purpose of implementing shooting method, we added extra conditions at the grid points T l , ðl ¼ 0; . . . ; L 1Þ, so that we have the following set of initial value problems
L½xl ¼ f ðtÞ;
t 2 ½T l ; T lþ1
ð19Þ
subject to
xl ðT l Þ ¼ a^ 2l ; x0l ðT l Þ ¼ a^ 2lþ1
ð20Þ
^ j ; ðj ¼ 1; . . . ; 2L 1Þ are unknown real parameters and a ^ 0 ¼ a0 . The solution of initial value problems (19) and (20) where a ^ l ; ðl ¼ 1; . . . ; 2L 1Þ, will be will be obtained by the method described in the previous subsection where the parameters a determined by solving the following system of algebraic equations
x0 ðT 1 ; a^ 1 Þ ¼ x1 ðT 1 ; a^ 2 ; a^ 3 Þ; x00 ðT 1 ; a^ 1 Þ ¼ x01 ðT 1 ; a^ 2 ; a^ 3 Þ; x1 ðT 2 ; a^ 2 ; a^ 3 Þ ¼ x2 ðT 2 ; a^ 4 ; a^ 5 Þ; x01 ðT 2 ; a^ 2 ; a^ 3 Þ ¼ x02 ðT 2 ; a^ 4 ; a^ 5 Þ; .. .
xk2 ðT k1 ; a^ 2k4 ; a^ 2k3 Þ ¼ xk1 ðT k1 ; a^ 2k2 ; a^ 2k1 Þ; x0k2 ðT k1 ; a^ 2k4 ; a^ 2k3 Þ ¼ x0k1 ðT k1 ; a^ 2k2 ; a^ 2k1 Þ; xk1 ðT; a^ 2k2 ; a^ 2k1 Þ ¼ a1 which can be solved numerically using several numerical techniques such as Newton’s method. In summary, this technique is called multiple shooting technique of order L. Note that, we use shooting method of low order ðL ¼ 4Þ in our computations which are done using Matlab. For more details about the multiple shooting method see [5,15,16].
4. Numerical results In this section we will consider the Bagley–Torvik equation which arises in the modeling of the motion of a rigid plate immersed in Newtonian fluid, see Bagley and Torvik [7]. Example 1. Consider the Bagley–Torvik equation
1 1 D2 yðtÞ þ D3=2 yðtÞ þ yðtÞ ¼ 8v ðtÞ 8v ðt 1Þ; 2 2
t 2 ½0; 20;
subject to
yð0Þ ¼ 0;
yð20Þ ¼ 1:48433;
where v ðtÞ is the Heaviside function. Applying the multiple shooting method of order four requires dividing the domain ½0; 20 into four subintervals as follows
0 ¼ T 0 < T 1 ¼ 5 < T 2 ¼ 10 < T 3 ¼ 15 < T 4 ¼ 20: Consequently, our problem is divided into the following initial value problems
L½x1 ¼ v ðtÞ; L½x2 ¼ v ðtÞ; L½x3 ¼ v ðtÞ;
^1; 0 < t < 5; x1 ð0; :Þ ¼ 0; x01 ð0; :Þ ¼ a ^ 2 ; x02 ð5; :Þ ¼ a ^3; 5 < t < 10; x2 ð5; :Þ ¼ a ^ 4 ; x03 ð10; :Þ ¼ a ^5; 10 < t < 15; x3 ð10; :Þ ¼ a
L½x4 ¼ v ðtÞ;
15 < t < 20;
x4 ð15; :Þ ¼ a^ 6 ; x04 ð15; :Þ ¼ a^ 7 :
where
1 1 L½xl ¼ D2 xl ðtÞ þ D3=2 xl ðtÞ þ xl ðtÞ: 2 2 and
ð21Þ
Q.M. Al-Mdallal et al. / Commun Nonlinear Sci Numer Simulat 15 (2010) 3814–3822
3819
vðtÞ ¼ v ðtÞ 8v ðt 1Þ: ^ 1 Þ; x2 ðt; :Þ ¼ x2 ðt; a ^2 ; a ^ 3 Þ; x3 ðt; :Þ ¼ x3 ðt; a ^4; a ^ 5 Þ; and x4 ðt; :Þ ¼ x4 ðt; a ^6; a ^ 7 Þ. Each initial value probHere x1 ðt; :Þ ¼ x1 ðt; a lem of (21) will be solved using the collocation method where we use m ¼ 3 with the collocation parameters
2
c1
2x
3
1 þ1
6 7 4 c2 5 ¼ c3
3
6 x 2þ1 7 6 2 7 4 2 5 x3 þ1 2
where xk ¼ cosð2kþ1 pÞ; ðk ¼ 1; 2; 3Þ are the roots of the Tchebychev polynomial of degree 3. The choice of m ¼ 3 is due to the 6 fact that although the more collocation points are used the more accurate the resulting numerical solution is to be expected, however, this will cause round off errors to have serious effects on the computed results. Needless to say that it is difficult to ^ i ; i ¼ 1; 2; . . . ; 7 are determined by solving the find a specific optimum choice for m. Finally, the values for the parameters a following algebraic system
x1 ð5; a^ 1 Þ ¼ x2 ð5; a^ 2 ; a^ 3 Þ; x01 ð5; a^ 1 Þ ¼ x02 ð5; a^ 2 ; a^ 3 Þ; x2 ð10; a^ 2 ; a^ 3 Þ ¼ x3 ð10; a^ 4 ; a^ 5 Þ; x02 ð10; a^ 2 ; a^ 3 Þ ¼ x03 ð10; a^ 4 ; a^ 5 Þ; x3 ð15; a^ 4 ; a^ 5 Þ ¼ x4 ð15; a^ 6 ; a^ 7 Þ; x03 ð15; a^ 4 ; a^ 5 Þ ¼ x04 ð15; a^ 6 ; a^ 7 Þ; x4 ð20; a^ 6 ; a^ 7 Þ ¼ 1:48433: Using Matlab package we obtain
a^ 1 ¼ 0; a^ 2 ¼ 2:9439356637512617; a^ 3 ¼ 3:50426674344817; a^ 4 ¼ 2:848380869392841; a^ 5 ¼ 1:858241791552473; a^ 6 ¼ 2:168701764974635; a^ 7 ¼ 0:8360760617165397: The graph of the approximated solution is given in Fig. (1). Note that this example was solved theoretically by Podlubny [30] where the compact form of the solution is given by
yðtÞ ¼
Z
t
8Gðt sÞðv ðsÞ v ðs 1ÞÞds
0
where G is the fractional Green’s function defined as
GðtÞ ¼
1 1 X ð1=2Þk X ðj þ kÞ!ðt1=2 Þj : j 1 1 3 k! j¼0 j!2 Cð2 j þ 2 k þ 2 k þ 2Þ k¼0
Here C is the Gamma function. The error between the exact solution and computed one is shown in Fig. 2. It is clearly seen that the two solutions are in excellent agreement. The computed L2 error norm is given by
kuðtÞ yðtÞk ¼
Z
20
ðyðtÞ uðtÞÞ2 dt ¼ 1:7 1010 :
0
Fig. 1. The approximate solution for Example 1.
3820
Q.M. Al-Mdallal et al. / Commun Nonlinear Sci Numer Simulat 15 (2010) 3814–3822
Fig. 2. Computed absolute error between the present numerical solution and the theoretical one by Podlubny [30] for Example 1.
Example 2. Consider the Bagley–Torvik equation
D2 yðtÞ þ D3=2 yðtÞ þ yðtÞ ¼ f ðtÞ;
t 2 ½0; 5;
ð22Þ
subject to
yð0Þ ¼ 0;
yð5Þ ¼ 25
2
where f ðtÞ ¼ t þ 2 þ 4
qffiffiffi t
p
and the exact solution yðtÞ ¼ t 2 .
Using the multiple shooting method of order five then following the same steps in the previous example, the parameters
a^ i ; i ¼ 1; 2; . . . ; 9 are found to be
a^ 1 ¼ 0; a^ 2 ¼ 0:99999999999999; a^ 3 ¼ 1:99999999999999; a^ 4 ¼ 3:99999999999999; a^ 5 ¼ 4; a^ 6 ¼ 8:99999999999999; a^ 7 ¼ 6:0000000000000001; a^ 8 ¼ 16; a^ 9 ¼ 7:99999999999999: The graph of the approximated solution is given in Fig. 3. In addition, Fig. 4 shows the error between the computed and the exact solutions which indicates the efficiency of the present numerical technique. The computed L2 error norm is given by
kuðtÞ yðtÞk ¼
Z
20
ðyðtÞ uðtÞÞ2 dt ¼ 3:78 1012 :
0
Fig. 3. The approximate solution for Example 2.
3821
Q.M. Al-Mdallal et al. / Commun Nonlinear Sci Numer Simulat 15 (2010) 3814–3822
Fig. 4. Computed absolute error between the present numerical solution and the theoretical one for Example 2.
Table 1 Effect of the initial condition, yð0Þ, on the computed L2 error norm. yð0Þ
kuðtÞ yðtÞk (Example 1)
kuðtÞ yðtÞk (Example 2)
0.000
10
1:71 10
3:78 1012
0.001
2:10 1010
3:61 1012
0.001
3:84 109
2:94 1012
0.100
7:14 1010
7:22 1012
0.500
2:92 1010
1:10 1012
1.000
8:73 1010
2:55 1012
In order to demonstrate the effect of the initial condition on the numerical results, Table 1 has been obtained for the above examples. It can be concluded that the effect of the initial condition is not significant on the error. 5. Conclusion In this paper, we have solved special class of fractional linear differential equation of order 2 subject to boundary conditions. A well-known example involving this class is the Bagley–Torvik equation. The method of solution is based on conjugating collocation and spline technique with shooting method. The numerical results for the examples demonstrate the efficiency and accuracy of the present method. The proposed method has the advantage of being both effective and flexible for solving boundary value problems governed by ordinary differential equations of fractional orders. The method combines the efficacy of multiple shooting method along with the flexibility of the collocation approach. To the best of the authors’ knowledge such combination can be considered as a new contribution. Acknowledgement The authors would like to express their appreciation for the valuable comments of the reviewers. The authors also would like to express their sincere appreciation to the United Arab Emirates University Research Affairs for the financial support of Grant nos. 04 13 02 11=07; 01 13 02 11=07. References [1] Ahmed E, Elgazzar AS. On fractional order differential equations model for nonlocal epidemics. Phys A 2007;379(2):607–14. [2] Ahmed E, El-Sayed AMA, Elsaka HAA. Equilibrium points, stability, and numerical solutions of fractional-order predator-prey and rabies models. J Math Anal Appl 2007;325(1):542–53. [3] Al-Mdallal QM. An efficient method for solving fractional Sturm-Liouville problems. Chaos, Solitons Fractals 2009;40(1):183–9. [4] Al-Mdallal QM. On the numerical solution of fractional Sturm-Liouville problems. Int J Comput Math 2009;99999:1. [5] Attili B, Syam M. Efficient shooting method for solving two point boundary value problems. Chaos, Solitons Fractals 2008;35(5):895–903. [6] Babenko Y. Non integer differential equation in engineering, Conference Bordeaux 3–8 July 1994. [7] Bagley RL, Torvik PJ. On the appearance of the fractional derivative in the behavior of real materials. ASME J Appl Mech 1984;51(2):294–8. [8] Beyer H, Kempfle S. Definition of physically consistent damping laws with fractional derivatives. Zeitschrift fr Angewandte Mathematik und Mechanik 1995;75(8):623–35.
3822
Q.M. Al-Mdallal et al. / Commun Nonlinear Sci Numer Simulat 15 (2010) 3814–3822
[9] Blank L. Numerical treatment of differential equations of fractional order. Numerical Analysis Report 287, Manchester Center for Numerical Computational Mathematics; 1996. [10] Daftardar-Geiji V, Jafari H. Adomian decomposition: a tool for solving a system of fractional differential equations. J Math Anal Appl 2005;301(2):508–18. [11] Daftardar-Geiji V, Jafari H. Analysis of a system of nonautonomous fractional differential equations involving Capato derivatives. J Math Anal Appl 2007;328:1026–33. [12] Diethelm K. An algorithm for the numerical solution of differential equations of fractional order. Elect Translat Numer Anal 1997;5:1–6. [13] Diethelm K, Ford NJ, Freed AD. A predictor–corrector approach for the numerical solution of fractional differential equations. Nonlinear Dynam 2002;29:3–22. [14] Ghorbani A, Alavi A. Applications of He’s variational iteration method to solve semidifferential equations of nth order. Math Prob Eng 2008;2008:1–9. [15] Ha SN. A nonlinear shooting method for two point boundary value problems. Comput Math Appl 2001;42:1411–20. [16] Keller HB. Numerical methods for two point boundary value problems. New York: Dover; 1992. [17] Lia M, Jimenezc S, Niea N, Tanga Y,Vazqueze L. Solving Two-point boundary value problems of fractional differential equations by spline collocation methods. Available from