- Email: [email protected]
Numerical solutions for systems of fractional differential equations by the decomposition method
Applied Mathematics and Computation 162 (2005) 1351–1365
www.elsevier.com/locate/amc
Numerical solutions for systems of fractional differential equations by the decomposition method Shaher Momani
*,1,
Kamel Al-Khaled
2
Department of Mathematics, United Arab Emirates University, P.O. Box 17551, Al-Ain, United Arab Emirates
Abstract In this paper we use Adomian decomposition method to solve systems of nonlinear fractional differential equations and a linear multi-term fractional differential equation by reducing it to a system of fractional equations each of order at most unity. We begin by showing how the decomposition method applies to a class of nonlinear fractional differential equations and give two examples to illustrate the efficiency of the method. Moreover, we show how the method can be applied to a general linear multi-term equation and solve several applied problems. Ó 2004 Elsevier Inc. All rights reserved. Keywords: Fractional differential equation; Decomposition method; Multi-term equations
1. Introduction Differential equations involving derivatives of noninteger order have shown to be adequate models for various physical phenomena in areas like rheology, damping laws, diffusion processes, etc. A review of some applications of fractional calculus in continuum and statistical mechanics is given by Mainardi [15]. *
Corresponding author. E-mail addresses: [email protected] (S. Momani), [email protected] (K. AlKhaled). 1 On leave from Department of Mathematics, Mutah University, P.O. Box 7, Mutah, Jordan. 2 On leave from Department of Mathematics, Jordan University of Science and Technology, IRBID 22110, Jordan. 0096-3003/$ - see front matter Ó 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2004.03.014
1352
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
The numerical solution of differential equations of an integer order has for a long time been a standard topic in numerical and computational mathematics. However, in spite of a large number of recently formulated applied problems, the state of art is far less advanced for generalized order equations, and only a small number of algorithms for the numerical solution of such equations has been suggested, and most of these numerical schemes deals with linear single term equations of order less than unity (see for example [3,6,8,16]). Recently, there has been some attempt to solve linear problems with multiple fractional derivatives problems [9,10,20]. Not much has been done for the Nonlinear problems and only two numerical schemes have been proposed to solve the nonlinear fractional differential equations [11,21]. In this paper, we will consider a system of nonlinear fractional differential equations. This fractional system of equations is obtained by replacing the first time derivative term by a fractional derivative of order a > 0. The Derivatives are understood in the Caputo sense. The general response expression contains a parameter describing the order of the fractional derivative that can be varied to obtain various responses. In the case of a ¼ 1, the fractional system equations reduces to the standard system of ordinary differential equations. Moreover, we consider a general linear multi-term equation by reducing the equation to a system of ordinary and fractional equations of order at most unity. We introduce several examples to illustrate the main ideas of this work. The Adomian decomposition method [1,2] will be applied for computing solutions to the systems of fractional differential equations considered in this paper. The method provide the solutions in the form of a power series with easily computed terms. The Adomian decomposition method has many advantages over the classical techniques mainly, it avoids discretization and provides an efficient numerical solution with high accuracy, minimal calculations. For more details in using the decomposition method for similar problems see [21,22]. The paper is organized as follows. A brief review of the fractional calculus theory is given in Section 2. In Section 3 we use the decomposition method to construct our numerical solutions for a system of nonlinear fractional equations. In Section 4 we use the decomposition method to construct our numerical solutions for a general linear multi-term fractional equation. We present some examples to show the efficiency and simplicity of the method.
2. Preliminaries and notations In order to proceed, we need the following definitions of fractional derivatives and integrals. We first introduce the Riemann–Liouville definition of fractional derivative operator Jaa .
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
1353
Definition 2.1. Let a 2 Rþ . The operator Jaa , defined on the usual Lebesque space L1 ½a; b by Z x 1 a1 ðx tÞ f ðtÞ dt; Jaa f ðxÞ ¼ CðaÞ a Ja0 f ðxÞ ¼ f ðxÞ; for a 6 x 6 b, is called the Riemann–Liouville fractional integral operator of order a. Properties of the operator J a can be found in [14], we mention the following: For f 2 L1 ½a; b, a; b P 0 and c > 1 1. 2. 3. 4.
Jaa f ðxÞ exists for almost every x 2 ½a; b. Jaa Jab f ðxÞ ¼ Jaaþb f ðxÞ, Jaa Jab f ðxÞ ¼ Jab Jaa f ðxÞ, Cðcþ1Þ Jaa ðx aÞc ¼ Cðaþcþ1Þ ðx aÞaþc ,
The Riemann–Liouville derivative have certain disadvantages when trying to model real-world phenomena with fractional differential equations. Therefore, we shall introduce now a modified fractional differential operator Da proposed by Caputo in his work on the theory of viscoelasticity [7]. For more information about Caputo definition and properties see [7,14,15]. Definition 2.2. The fractional derivative of f ðxÞ in the Caputo sense is defined as Da f ðxÞ ¼ J ma Dm f ðxÞ ¼
1 Cðm aÞ
Z
x
ðx tÞ
ma1 ðmÞ
0
for m 1 < a 6 m, m 2 N, x > 0. Also, we need here two of its basic properties. Lemma 2.1. If m 1 < a 6 m, and f 2 L1 ½a; b, then Daa Jaa f ðxÞ ¼ f ðxÞ; and, Jaa Daa f ðxÞ ¼ f ðxÞ
m1 X k¼0
f ðkÞ ð0þ Þ
ðx aÞk ; k!
x > 0:
f
ðtÞ dt;
ð2:1Þ
1354
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
The Caputo fractional derivative is considered here because it allows traditional initial and boundary conditions to be included in the formulation of the problem. In this paper, we are concerned with providing good quality algorithms for the solution of a system of fractional differential equations of the general form: da yi ¼ fi ðt; y1 ; y2 ; . . . ; yn Þ; dta
i ¼ 1; 2; . . . ; n;
ð2:2Þ
subject to the initial conditions yik ð0Þ ¼ cik ;
k ¼ 0; 1; . . . ; m 1; m 1 < a 6 m; m 2 N;
ð2:3Þ
where cik is a specified constant vector and yi ðtÞ is the solution vector. The fractional derivative in Eq. (2.2) is considered in the Caputo sense. The reason for adopting the Caputo definition is as follows: when we interpret the fractional derivative in Eq. (2.2) as Caputo fractional derivatives then, with suitable conditions on the forcing function fi ðt; y1 ; y2 ; . . . ; yn Þ and with initial values yik ð0Þ ¼ cik , k ¼ 0; 1; . . . ; m 1 specified, the system has a unique solution. If we were to interpret the fractional derivative as Riemann–Liouville fractional derivatives, we would have to specify our initial conditions in terms of fractional integrals and their derivatives. The initial conditions required by the Caputo definition coincide with identifiable physical states, and this leads to the preference, among modellers, for the Caputo definition [11]. For more details on the existence and uniqueness of solutions to the fractional differential equations (2.2) see [9,17–20].
3. Analysis of the numerical method The system of equations (2.2) can be written in terms of operator from as Dat yi ¼ fi ðt; y1 ; y2 ; . . . ; yn Þ;
i ¼ 1; 2; . . . ; n;
ð3:1Þ
where the fractional differential operator Dat is defined as in Eq. (2.1) denoted by Dat ¼
oa : ota
Following similar procedures for system of ordinary differential equation [13] and operating with J a ¼ J0a on both sides of Eq. (3.1) yields J a Dat yi ¼ J a fi ðt; y1 ; y2 ; . . . ; yn Þ: Therefore, it follows that
ð3:2Þ
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
yi ðtÞ ¼
m1 X
ðkÞ
yi ð0þ Þ
k¼0
tk þ J a fi ðt; y1 ðtÞ; y2 ðtÞ; . . . ; yn ðtÞÞ: k!
1355
ð3:3Þ
The Adomain’s decomposition method [1,2] assumes a series solution for yi ðtÞ given by yi ðtÞ ¼
1 X
fi;j ðt; y1 ðtÞ; y2 ðtÞ; . . . ; yn ðtÞÞ:
ð3:4Þ
j¼0
Substituting the initial condition (2.3) into (3.3) and identifying the zeroth component yi;0 by the term arising from the initial condition and from the source term, then we have the following recursive relations yi;0 ð0Þ ¼
m1 X
ðkÞ
yi ð0þ Þ
k¼0
tk ; k!
i ¼ 1; 2; . . . ; n;
yi;jþ1 ðtÞ ¼ J a fi;j ðt; y1 ; y2 ; . . . ; yn Þ;
j ¼ 0; 1; . . .
ð3:5Þ
The values of the natural number m can be 1 and 2 corresponding to 0 < a 6 1 and 1 < a 6 2, respectively. The components of yi;j ðtÞ, j P 1 are determined in the following recursive way yi;1 ¼ J a fi;0 ðt; y1 ; y2 ; . . . ; yn Þ; yi;2 ¼ J a fi;1 ðt; y1 ; y2 ; . . . ; yn Þ; yi;3 ¼ J a fi;2 ðt; y1 ; y2 ; . . . ; yn Þ; .. .
ð3:6Þ
As a result, the series solution is given by yi ðtÞ ¼ yi;0 þ
1 X
½J a fi;j ðt; y1 ðtÞ; y2 ðtÞ; . . . ; yn ðtÞÞ:
ð3:7Þ
j¼0
Define the c-term approximation solution as /i;c ¼
c1 X
yi;j ðtÞ;
ð3:8Þ
j¼0
and the exact solution yi ðtÞ is given by yi ðtÞ ¼ lim /i;c : c!1
ð3:9Þ
To demonstrate the effectiveness of the method for solving nonlinear systems of fractional differential equations, we consider here the following two examples.
1356
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
Example 3.1. Consider the following nonlinear system of fractional differential equations: da y1 ¼ y1 ; dta a d y2 ¼ y1 y22 ; dta da y3 ¼ y22 ; dta
ð3:10Þ
subject to the the initial condition y1 ð0Þ ¼ 1;
y2 ð0Þ ¼ 0;
y3 ð0Þ ¼ 0;
where 0 < a 6 1. For the special case when a ¼ 1, this system represents a nonlinear reaction and was found in [12]. In order to solve the above system and in addition to Eq. (3.5) we define the nonlinear term by P1 Ny ¼ y22 ¼ j¼0 Aj ðtÞ, and where Aj is the appropriate Adomian polynomials generated for the specific nonlinearities in Eq. (3.10), a definition for Aj is given as " 1 dj Aj ðt0 ; . . . ; tj ; y0 ; . . . ; yj Þ ¼ N j! dkj
j X k¼0
ktk
j X k¼0
!# kyk
;
j P 0:
k¼0
ð3:11Þ This formula is easy to compute by using Mathematica software or by setting a computer code to get as many polynomials as we need in the calculation of the numerical as well as explicit solutions. For more details about Adomian decomposition method and a general formula of Adomian polynomials, we refer the readers to [1,2]. In this example, we have Ny ¼ y22 and by using Eq. (3.11), we get 2 ; A0 ¼ y2;0
A1 ¼ 2y2;0 y2;1 ; 2 A2 ¼ y2;1 þ 2y2;0 y2;2 ;
A3 ¼ 2y2;0 y2;3 þ 2y2;1 y2;2 ; .. . Using the alternate algorithm for computing the Adomian polynomials. The Adomian decomposition series (3.5) leads to the following scheme:
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
y1;0 ¼ 1;
y1;jþ1 ¼ J a ½y1;j ;
y2;0 ¼ 0; y3;0 ¼ 0;
y2;jþ1 ¼ J a ½y1;j Aj ; y3;jþ1 ¼ J a ½Aj :
1357
ð3:12Þ
Using the above recursive relationship and Mathematica and upon using (3.9), the first four terms of the decomposition series are given by ta t2a t3a þ ; Cða þ 1Þ Cð2a þ 1Þ Cð3a þ 1Þ " # ta t2a Cð2a þ 1Þ t3a þ 1 ; y2 ¼ 2 Cða þ 1Þ Cð2a þ 1Þ Cða þ 1Þ Cð3a þ 1Þ
y1 ¼ 1
y3 ¼
t3a þ 2 Cða þ 1Þ Cð3a þ 1Þ
Cð2a þ 1Þ
Setting a ¼ 1, we obtain the solution obtained by Kaya [13] which corresponding to a system of ordinary differential equations. Fig. 1 shows two profile solutions for y1 , the solid line represents the solution for the integer derivative ða ¼ 1Þ, while the dotted line represents the fractional derivative when a ¼ 12. It is easy to conclude that the solution for fractional derivative ð0 < a < 1Þ blows slower than the integer derivative solution as time increases. Example 3.2. Consider the following nonlinear system of fractional differential equations:
y1 t 35 30 25 20 15 10 5 t 2
4
6
8
10
Fig. 1. The approximate solution y1 ðtÞ, for a ¼ 1 (solid line) and a ¼ 12 (dotted line).
1358
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
da y1 ¼ y1 þ y2 y3 ; dta a d y2 ¼ y1 y2 y3 2y22 ; dta da y3 ¼ y22 ; dta
ð3:13Þ
subject to the initial conditions y1 ð0Þ ¼ 1;
y2 ð0Þ ¼ 2;
y3 ð0Þ ¼ 0:
ð3:14Þ
Substituting the initial conditions (3.14) into Eq. (3.5) and using Eq. (3.11) to calculate the Adomian polynomials, yields the following recursive relations: y1;0 ¼ 1;
y1;jþ1 ¼ J a ½y2;j y3;j y1;j ;
y2;0 ¼ 2; y3;0 ¼ 0;
y2;jþ1 ¼ J a ½y1;j y2;j y3;j 2Aj ; y3;jþ1 ¼ J a ½Aj :
Using the above recursive relationship and Mathematica, the first few terms of the decomposition series are given by y1 ¼ 1
ta t2a 28Cð2a þ 1Þ t3a þ ; 2 Cða þ 1Þ Cð2a þ 1Þ Cða þ 1Þ Cð3a þ 1Þ
y2 ¼ 2
7ta 55t2a 28Cð2a þ 1Þ t3a þ ; 2 Cða þ 1Þ Cð2a þ 1Þ Cða þ 1Þ Cð3a þ 1Þ
y3 ¼
4ta 28t2a Cða þ 1Þ Cð2a þ 1Þ
4. Linear differential equation Consider the fractional differential equation of order na of the form: dðnaÞ yðtÞ dððn1ÞaÞ yðtÞ dðaÞ yðtÞ þ a þ þ a þ a0 yðtÞ ¼ f ðtÞ; n1 1 dtna dtðn1Þa dta
ð4:1Þ
with only the integer order initial conditions y ði1Þ ð0Þ ¼ ci1 ;
i ¼ 1; 2; . . . ; n;
ð4:2Þ
where ai and ci are arbitrary constants for i ¼ 0; 1; 2; . . . ; n 1. The numerical approximation of the solution of a linear multi-term fractional differential equation (4.1) and (4.2) can be calculated by reducing the
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
1359
problem to a system of ordinary and fractional equations each of order a, 0 < a 6 1. To convert Eq. (4.1) to a system of equations, we assume y1 ¼ y and ðaÞ ðaÞ ðaÞ y1 ¼ y2 , y2 ¼ y3 ; . . . ; yn1 ¼ yn then ynðaÞ ¼ a0 y1 a1 y2 an1 yn þ f ðtÞ; which defines a linear time-invariant system in the canonical form: 2 ðaÞ 3 32 y 3 2 3 2 y1 0 1 0 ... 0 0 0 1 6 ðaÞ 7 7 607 y 6 y2 7 6 0 0 1 ... 0 0 76 2 7 6 7 76 6 ðaÞ 7 6 y3 7 6 y3 7 6 0 0 0 ... 0 0 76 7 607 6 7 7 6 6 .. 7 .. .. .. 6 .. 7 þ 6 6 ... 7f ðtÞ; 7 6 .. 7 ¼ 6 ... . 76 . ... . . . 7 6 7 6 . 7 6 7 6 6 ðaÞ 7 4 0 0 0 ... 0 1 54 yn1 5 4 0 5 4y 5 n1 a0 a1 a2 . . . an2 an1 1 yn ynðaÞ ð4:3Þ The linear system (4.3) can be expressed in a more compact form as: Y ðaÞ ¼ AY þ Bf ðtÞ:
ð4:4Þ
Operating with J a on both sides of Eq. (4.4) yields Y ¼ J a ½AY þ J a ½Bf ðtÞ:
ð4:5Þ
Following Adomian method analysis, the linear system (4.5) is transformed into a set of recursive relations given by Y0 ¼ C þ J a ½Bf ðtÞ; Yjþ1 ¼ J a ½AYj ; j P 0;
ð4:6Þ
where C is the column vector with ci , i ¼ 0; 1; . . . ; n 1. It is essential feature of the decomposition method that the zeroth component Y0 is defined always by all terms that arise from initial data and from integrating the inhomogeneous terms. Having defined the zeroth component Y0 , the remaining components Yj , j P 1, can be completely determined by using Eq. (4.6). 4.1. Applications Example 4.1. Consider the homogeneous differential equation of order 4a of the form: dð4aÞ yðtÞ dð2aÞ yðtÞ 2 þ yðtÞ ¼ 0; dt4a dt2a yð0Þ ¼ 0; y ð0Þ ¼ 1:
ð4:7Þ
1360
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365 1=2
1=2
1=2
If we choose a ¼ 12, and yðtÞ ¼ y1 , y1 ¼ y2 , y2 ¼ y3 , and y3 (4.7) can be reduced to the following matrix form: 2 ð1=2Þ 3 2 32 3 y1 0 1 0 0 y1 6 ð1=2Þ 7 6 6 y2 7 6 y2 7 6 0 0 1 0 7 76 7: 6 ð1=2Þ 7 ¼ 4 0 0 0 1 54 y3 5 4 y3 5 ð1=2Þ 1 0 2 0 y4 y
¼ y4 , then Eq.
ð4:8Þ
4
Operating with J 1=2 on both sides of the system of equations (4.8) and comparing with Eq. (4.6) we get y1;0 ¼ 0;
y1;jþ1 ¼ J ð1=2Þ ½y2;j ;
y2;0 ¼ 0;
y2;jþ1 ¼ J ð1=2Þ ½y3;j ;
y3;0 ¼ 1;
y3;jþ1 ¼ J ð1=2Þ ½y4;j ;
y4;0 ¼ 0;
y4;jþ1 ¼ J ð1=2Þ ½2y3;j y1;j :
The first few terms of the decomposition series are given by t2 t3 y1 ¼ t 1 þ t þ þ þ ; 2 6 t1=2 2t3=2 3t5=2 þ þ þ ; Cð3=2Þ Cð5=2Þ Cð7=2Þ 3t2 þ ; y3 ¼ 1 þ 2t þ 2 2t1=2 3t3=2 5t5=2 þ þ þ : y4 ¼ Cð3=2Þ Cð5=2Þ Cð7=2Þ y2 ¼
It is easily verified that the exact solution of this problem is yðtÞ ¼ tet : Example 4.2. We next consider the following initial value problem for the inhomogeneous Bagley–Torvik equation [4]: d2 yðtÞ dð3=2Þ yðtÞ þ B þ CyðtÞ ¼ f ðtÞ; dt2 dt3=2 yð0Þ ¼ 1; y ð0Þ ¼ 1;
A
0 6 t 6 5;
ð4:9Þ
pffiffiffiffiffiffi where A ¼ M, B ¼ 2S lq, and C ¼ K. This problem describes the motion of a large plate of the surface S and mass M in a Newtonian fluid of viscosity l and density q. The plate is hanging on a massless spring of stiffness K. The function f ðtÞ represents the loading force. In
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
1361
order to make comparison with the numerical solution of [9], we choose A ¼ B ¼ C ¼ 1 and f ðtÞ ¼ Cðt þ 1Þ. As in the previous example, let yðtÞ ¼ y1 , y11=2 ¼ y2 , y21=2 ¼ y3 , and y31=2 ¼ y4 , then Eq. (4.9) can be reduced to the following matrix form: 2 ð1=2Þ 3 2 32 3 2 3 y1 0 1 0 0 y1 0 6 ð1=2Þ 7 6 6 y2 7 6 0 7 6 y2 7 6 0 0 1 0 7 76 7 þ 6 7ðt þ 1Þ: ð4:10Þ 6 ð1=2Þ 7 ¼ 4 0 0 0 1 54 y3 5 4 0 5 4 y3 5 ð1=2Þ 1 0 0 1 1 y4 y 4
Operating with J 1=2 on both sides of the system of equations (4.10) and comparing with Eq. (4.6) we obtain y1;0 ¼ 1;
y1;jþ1 ¼ J ð1=2Þ ½y2;j ;
y2;0 ¼ 0;
y2;jþ1 ¼ J ð1=2Þ ½y3;j ;
y3;0 ¼ 1;
y3;jþ1 ¼ J ð1=2Þ ½y4;j :
Since the computation depends heavily on y4;0 , we will use here a slight modification [23] for defining the components of y4;j by choosing y4;0 to be the simplest term in (4.10) and add the remaining terms to y4;1 and y4;2 . This will ease the computation considerably. Thus y4;0 ¼ 0; t1=2 þ J ð1=2Þ ½y4;0 y1;0 ; Cð3=2Þ t3=2 þ J ð1=2Þ ½y4;1 y1;1 ; ¼ Cð3=2Þ
y4;1 ¼ y4;2
y4;jþ1 ¼ J ð1=2Þ ½y4;j y1;j ;
j P 2:
The first few terms of the decomposition series are given by y1 ¼ 1 þ t; t1=2 ; Cð3=2Þ t2 t2 y3 ¼ 1 þ ; 2 2 t3=2 t3=2 : y4 ¼ Cð3=2Þ Cð5=2Þ y2 ¼
It is easily verified that the exact solution of this problem is yðtÞ ¼ t þ 1:
1362
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
Table 1 Step size
Decomposition method
Error at t ¼ 5 by [9]
0.5 0.25 0.125 0.0625
0 0 0 0
)0.15131473519232 )0.04684102179946 )0.01602947553912 )0.00562770408881
Table 1 shows the resulting error at t ¼ 5 obtained by the numerical method in [9] and compared with the exact solution obtained by the decomposition method. This demonstrates the importance of our algorithm in solving linear multi-term fractional differential equations. Example 4.3. Consider the fractional differential equations of the form: dyðtÞ da yðtÞ þ a1 þ a2 yðtÞ ¼ f ðtÞ; dt dta yð0Þ ¼ c; 0 < a < 1;
t P 0;
ð4:11Þ
where a1 and a2 are positive constants and c is an arbitrary constant. This fractional differential equation corresponds to the Basset Problem [5]. It represents a classical problem in fluid dynamics where the unsteady motion of a particle accelerates in a viscous fluid due to the gravity of force. In this example, we choose f ðtÞ ¼ 1 and a vanishing initial velocity ðc ¼ 0Þ and consider the case for the ordinary Basset problem, a ¼ 1=2. Eq. (4.11) can be converted into the following system of equations y1 ¼ yðtÞ; ð1=2Þ
y1
¼ y2 ;
ð1=2Þ y2
¼ a2 y1 a1 y2 þ 1:
ð4:12Þ
Operating with J 1=2 on both sides of Eq. (4.12) yields y1 ¼ J 1=2 y2 ; y2 ¼ J 1=2 ð1Þ þ J 1=2 ½a2 y1 a1 y2 :
ð4:13Þ
The Adomian procedures leads to the following scheme: y1;0 ¼ 0;
y1;jþ1 ¼ J ð1=2Þ ½y2;j ;
y2;0 ¼ J ð1=2Þ ð1Þ;
y2;jþ1 ¼ J ð1=2Þ ½a2 y1;j a1 y2;j :
ð4:14Þ
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
1363
The first five terms of the decomposition series are given by: y1 ¼ B0 t þ B1 t3=2 þ B2 t2 þ B3 t5=2 þ B4 t3 þ ; y2 ¼ D0 t1=2 þ D1 t þ D1 t3=2 þ D2 t3=2 þ D3 t2 þ D4 t5=2 þ ;
ð4:15Þ
where 1 ; Cð2Þ a1 B1 ¼ ; Cð5=2Þ B0 ¼
B2 ¼
1 ða2 a2 Þ; Cð3Þ 1
B3 ¼
1 ð2a1 a2 a31 Þ; Cð7=2Þ
B4 ¼
1 ða2 3a2 a1 þ a41 Þ; Cð4Þ 2
1 ; Cð3=2Þ a1 D1 ¼ ; Cð2Þ
D0 ¼
D2 ¼
1 ða2 a2 Þ; Cð5=2Þ 1
D3 ¼
1 ð2a1 a2 a31 Þ; Cð3Þ
D4 ¼
1 ða2 3a2 a1 þ a41 Þ: Cð7=2Þ 2
Using the formulation of the problem given by Mainardi [15], we let a2 ¼
9 1 þ 2v
1=2 ;
where v ¼ 10, 100. Fig. 2 shows the profile solution for y1 . It is easy to conclude that our solution is in agreement with the analytical solution of Mainardi [15] and the numerical solution of Edwards et al. [11] for short times. We should mention that only 5 terms of the decomposition series were used for our approximation solution. It is evident that the solution can be improved by adding more terms of the decomposition series.
1364
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
y t 0.5 0.4 0.3 0.2 0.1 t 0.2
0.4
0.6
0.8
1
Fig. 2. The approximate solution yðtÞ, for v ¼ 10 (solid line) and v ¼ 100 (dotted line).
5. Conclusions The fundamental goal of this work has been to propose an efficient algorithm for the solution of linear and nonlinear systems of fractional differential equations. This goal has been achieved by using Adomian decomposition method and an approximation series solution can be obtained to any desired number of terms. For linear systems of equations the algorithm gives exact solution, and for nonlinear equations it provides an approximate solution with good accuracy. Moreover, we have shown that the use of systems of fractional differential equations of different degrees can provide a convergent method for the solution of multi-term equations. References [1] G. Adomian, A review of the decomposition method in applied mathematics, J. Math. Anal. Appl. 135 (1988) 501–544. [2] G. Adomian, Solving Frontier Problems of Physics: The Decomposition Method, Kluwer Academic Publishers, Boston, 1994. [3] C.T.H. Baker, M.S. Derakhshan, Stability barriers to the construction of q,r-reducible and fractional quadrature rules, in: Numerical Integration III, in: H. Brab, B.G. Hammerlin (Eds.), Ser. Numer. Math., vol. 85, Birkhauser, Basel, 1988, pp. 1–15. [4] R.L. Bagley, P.J. Torvik, On the appearance of the fractional derivative in the behavior of real materials, J. Appl. Mech. 51 (1984) 294–298. [5] A.B. Basset, On the descent of a sphere in a viscous liquid, Quart. J. Math. 42 (1910) 369–381. [6] L. Blank, Numerical treatment of differential equations of fractional order, MCCM Numerical Analysis Report No. 287, The University of Manchester, 1996. [7] M. Caputo, Linear models of dissipation whose Q is almost frequency independent. Part II, J. Roy. Astr. Soc. 13 (1967) 529–539. [8] K. Diethelm, An algorithm for the numerical solution of differential equations of fractional order, Electron. Trans. Numer. Anal. 5 (1997) 1–6.
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
1365
[9] K. Diethelm, N.J. Ford, Numerical solution of the Bagley–Torvik equation, BIT 42 (2002) 490–507. [10] K. Diethelm, Y. Luchko, Numerical solution of linear multi-term differential equations of fractional order, J. Comput. Anal. Appl., in press. [11] J.T. Edwards, N.J. Ford, A.C. Simpson, The numerical solution of linear multi-term fractional differential equations: systems of equations, J. Comput. Appl. Math. 148 (2002) 401–418. [12] T.E. Hull, W.H. Enright, B.M. Fellen, A.E. Sedgwick, Comparing numerical methods for ordinary differential equations, SIAM J. Numer. Anal. 9 (1972) 603. [13] D. Kaya, A reliable method for the numerical solution of the kinetics problems, Appl. Math. Comput, in press. [14] AY. Luchko, R. Groreflo, The initial value problem for some fractional differential equations with the Caputo derivative, Preprint series A08–98, Fachbreich Mathematik und Informatik, Freic Universitat Berlin, 1998. [15] F. Mainardi, Fractional calculus: ‘Some basic problems in continuum and statistical mechanics’, in: A. Carpinteri, F. Mainardi (Eds.), Fractals and Fractional Calculus in Continuum Mechanics, Springer-Verlag, New York, 1997, pp. 291–348. [16] R. Metzler;, J. Klafter, Boundary value problems fractional diffusion equations, Physica A 278 (2000) 107–125. [17] K.S. Miller, B. Ross, An Introduction to the Fractional Calculus and Fractional Differential Equations, John Wiley and Sons, Inc., New York, 1993. [18] S.M. Momani, On the existence of solutions of a system of ordinary differential equations of fractional order, Far East J. Math. Sci. 1 (2) (1999) 265–270. [19] K.B. Oldham, J. Spanier, The Fractional Calculus, Academic Press, New York, 1974. [20] I. Podlubny, Fractional Differential Equations, Academic Press, New York, 1999. [21] N.T. Shawagfeh, Analytical approximate solutions for nonlinear fractional differential equations, Appl. Math. Comput. 131 (2–3) (2002) 517–529. [22] A.M. Wazwaz, Blow-up for solutions of some linear wave equations with mixed nonlinear boundary conditions, Appl. Math. Comput. 123 (2001) 133–140. [23] A.M. Wazwaz, A reliable modification of Adomian’s decomposition method, Appl. Math. Comput. 92 (1998) 1–7.
www.elsevier.com/locate/amc
Numerical solutions for systems of fractional differential equations by the decomposition method Shaher Momani
*,1,
Kamel Al-Khaled
2
Department of Mathematics, United Arab Emirates University, P.O. Box 17551, Al-Ain, United Arab Emirates
Abstract In this paper we use Adomian decomposition method to solve systems of nonlinear fractional differential equations and a linear multi-term fractional differential equation by reducing it to a system of fractional equations each of order at most unity. We begin by showing how the decomposition method applies to a class of nonlinear fractional differential equations and give two examples to illustrate the efficiency of the method. Moreover, we show how the method can be applied to a general linear multi-term equation and solve several applied problems. Ó 2004 Elsevier Inc. All rights reserved. Keywords: Fractional differential equation; Decomposition method; Multi-term equations
1. Introduction Differential equations involving derivatives of noninteger order have shown to be adequate models for various physical phenomena in areas like rheology, damping laws, diffusion processes, etc. A review of some applications of fractional calculus in continuum and statistical mechanics is given by Mainardi [15]. *
Corresponding author. E-mail addresses: [email protected] (S. Momani), [email protected] (K. AlKhaled). 1 On leave from Department of Mathematics, Mutah University, P.O. Box 7, Mutah, Jordan. 2 On leave from Department of Mathematics, Jordan University of Science and Technology, IRBID 22110, Jordan. 0096-3003/$ - see front matter Ó 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2004.03.014
1352
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
The numerical solution of differential equations of an integer order has for a long time been a standard topic in numerical and computational mathematics. However, in spite of a large number of recently formulated applied problems, the state of art is far less advanced for generalized order equations, and only a small number of algorithms for the numerical solution of such equations has been suggested, and most of these numerical schemes deals with linear single term equations of order less than unity (see for example [3,6,8,16]). Recently, there has been some attempt to solve linear problems with multiple fractional derivatives problems [9,10,20]. Not much has been done for the Nonlinear problems and only two numerical schemes have been proposed to solve the nonlinear fractional differential equations [11,21]. In this paper, we will consider a system of nonlinear fractional differential equations. This fractional system of equations is obtained by replacing the first time derivative term by a fractional derivative of order a > 0. The Derivatives are understood in the Caputo sense. The general response expression contains a parameter describing the order of the fractional derivative that can be varied to obtain various responses. In the case of a ¼ 1, the fractional system equations reduces to the standard system of ordinary differential equations. Moreover, we consider a general linear multi-term equation by reducing the equation to a system of ordinary and fractional equations of order at most unity. We introduce several examples to illustrate the main ideas of this work. The Adomian decomposition method [1,2] will be applied for computing solutions to the systems of fractional differential equations considered in this paper. The method provide the solutions in the form of a power series with easily computed terms. The Adomian decomposition method has many advantages over the classical techniques mainly, it avoids discretization and provides an efficient numerical solution with high accuracy, minimal calculations. For more details in using the decomposition method for similar problems see [21,22]. The paper is organized as follows. A brief review of the fractional calculus theory is given in Section 2. In Section 3 we use the decomposition method to construct our numerical solutions for a system of nonlinear fractional equations. In Section 4 we use the decomposition method to construct our numerical solutions for a general linear multi-term fractional equation. We present some examples to show the efficiency and simplicity of the method.
2. Preliminaries and notations In order to proceed, we need the following definitions of fractional derivatives and integrals. We first introduce the Riemann–Liouville definition of fractional derivative operator Jaa .
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
1353
Definition 2.1. Let a 2 Rþ . The operator Jaa , defined on the usual Lebesque space L1 ½a; b by Z x 1 a1 ðx tÞ f ðtÞ dt; Jaa f ðxÞ ¼ CðaÞ a Ja0 f ðxÞ ¼ f ðxÞ; for a 6 x 6 b, is called the Riemann–Liouville fractional integral operator of order a. Properties of the operator J a can be found in [14], we mention the following: For f 2 L1 ½a; b, a; b P 0 and c > 1 1. 2. 3. 4.
Jaa f ðxÞ exists for almost every x 2 ½a; b. Jaa Jab f ðxÞ ¼ Jaaþb f ðxÞ, Jaa Jab f ðxÞ ¼ Jab Jaa f ðxÞ, Cðcþ1Þ Jaa ðx aÞc ¼ Cðaþcþ1Þ ðx aÞaþc ,
The Riemann–Liouville derivative have certain disadvantages when trying to model real-world phenomena with fractional differential equations. Therefore, we shall introduce now a modified fractional differential operator Da proposed by Caputo in his work on the theory of viscoelasticity [7]. For more information about Caputo definition and properties see [7,14,15]. Definition 2.2. The fractional derivative of f ðxÞ in the Caputo sense is defined as Da f ðxÞ ¼ J ma Dm f ðxÞ ¼
1 Cðm aÞ
Z
x
ðx tÞ
ma1 ðmÞ
0
for m 1 < a 6 m, m 2 N, x > 0. Also, we need here two of its basic properties. Lemma 2.1. If m 1 < a 6 m, and f 2 L1 ½a; b, then Daa Jaa f ðxÞ ¼ f ðxÞ; and, Jaa Daa f ðxÞ ¼ f ðxÞ
m1 X k¼0
f ðkÞ ð0þ Þ
ðx aÞk ; k!
x > 0:
f
ðtÞ dt;
ð2:1Þ
1354
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
The Caputo fractional derivative is considered here because it allows traditional initial and boundary conditions to be included in the formulation of the problem. In this paper, we are concerned with providing good quality algorithms for the solution of a system of fractional differential equations of the general form: da yi ¼ fi ðt; y1 ; y2 ; . . . ; yn Þ; dta
i ¼ 1; 2; . . . ; n;
ð2:2Þ
subject to the initial conditions yik ð0Þ ¼ cik ;
k ¼ 0; 1; . . . ; m 1; m 1 < a 6 m; m 2 N;
ð2:3Þ
where cik is a specified constant vector and yi ðtÞ is the solution vector. The fractional derivative in Eq. (2.2) is considered in the Caputo sense. The reason for adopting the Caputo definition is as follows: when we interpret the fractional derivative in Eq. (2.2) as Caputo fractional derivatives then, with suitable conditions on the forcing function fi ðt; y1 ; y2 ; . . . ; yn Þ and with initial values yik ð0Þ ¼ cik , k ¼ 0; 1; . . . ; m 1 specified, the system has a unique solution. If we were to interpret the fractional derivative as Riemann–Liouville fractional derivatives, we would have to specify our initial conditions in terms of fractional integrals and their derivatives. The initial conditions required by the Caputo definition coincide with identifiable physical states, and this leads to the preference, among modellers, for the Caputo definition [11]. For more details on the existence and uniqueness of solutions to the fractional differential equations (2.2) see [9,17–20].
3. Analysis of the numerical method The system of equations (2.2) can be written in terms of operator from as Dat yi ¼ fi ðt; y1 ; y2 ; . . . ; yn Þ;
i ¼ 1; 2; . . . ; n;
ð3:1Þ
where the fractional differential operator Dat is defined as in Eq. (2.1) denoted by Dat ¼
oa : ota
Following similar procedures for system of ordinary differential equation [13] and operating with J a ¼ J0a on both sides of Eq. (3.1) yields J a Dat yi ¼ J a fi ðt; y1 ; y2 ; . . . ; yn Þ: Therefore, it follows that
ð3:2Þ
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
yi ðtÞ ¼
m1 X
ðkÞ
yi ð0þ Þ
k¼0
tk þ J a fi ðt; y1 ðtÞ; y2 ðtÞ; . . . ; yn ðtÞÞ: k!
1355
ð3:3Þ
The Adomain’s decomposition method [1,2] assumes a series solution for yi ðtÞ given by yi ðtÞ ¼
1 X
fi;j ðt; y1 ðtÞ; y2 ðtÞ; . . . ; yn ðtÞÞ:
ð3:4Þ
j¼0
Substituting the initial condition (2.3) into (3.3) and identifying the zeroth component yi;0 by the term arising from the initial condition and from the source term, then we have the following recursive relations yi;0 ð0Þ ¼
m1 X
ðkÞ
yi ð0þ Þ
k¼0
tk ; k!
i ¼ 1; 2; . . . ; n;
yi;jþ1 ðtÞ ¼ J a fi;j ðt; y1 ; y2 ; . . . ; yn Þ;
j ¼ 0; 1; . . .
ð3:5Þ
The values of the natural number m can be 1 and 2 corresponding to 0 < a 6 1 and 1 < a 6 2, respectively. The components of yi;j ðtÞ, j P 1 are determined in the following recursive way yi;1 ¼ J a fi;0 ðt; y1 ; y2 ; . . . ; yn Þ; yi;2 ¼ J a fi;1 ðt; y1 ; y2 ; . . . ; yn Þ; yi;3 ¼ J a fi;2 ðt; y1 ; y2 ; . . . ; yn Þ; .. .
ð3:6Þ
As a result, the series solution is given by yi ðtÞ ¼ yi;0 þ
1 X
½J a fi;j ðt; y1 ðtÞ; y2 ðtÞ; . . . ; yn ðtÞÞ:
ð3:7Þ
j¼0
Define the c-term approximation solution as /i;c ¼
c1 X
yi;j ðtÞ;
ð3:8Þ
j¼0
and the exact solution yi ðtÞ is given by yi ðtÞ ¼ lim /i;c : c!1
ð3:9Þ
To demonstrate the effectiveness of the method for solving nonlinear systems of fractional differential equations, we consider here the following two examples.
1356
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
Example 3.1. Consider the following nonlinear system of fractional differential equations: da y1 ¼ y1 ; dta a d y2 ¼ y1 y22 ; dta da y3 ¼ y22 ; dta
ð3:10Þ
subject to the the initial condition y1 ð0Þ ¼ 1;
y2 ð0Þ ¼ 0;
y3 ð0Þ ¼ 0;
where 0 < a 6 1. For the special case when a ¼ 1, this system represents a nonlinear reaction and was found in [12]. In order to solve the above system and in addition to Eq. (3.5) we define the nonlinear term by P1 Ny ¼ y22 ¼ j¼0 Aj ðtÞ, and where Aj is the appropriate Adomian polynomials generated for the specific nonlinearities in Eq. (3.10), a definition for Aj is given as " 1 dj Aj ðt0 ; . . . ; tj ; y0 ; . . . ; yj Þ ¼ N j! dkj
j X k¼0
ktk
j X k¼0
!# kyk
;
j P 0:
k¼0
ð3:11Þ This formula is easy to compute by using Mathematica software or by setting a computer code to get as many polynomials as we need in the calculation of the numerical as well as explicit solutions. For more details about Adomian decomposition method and a general formula of Adomian polynomials, we refer the readers to [1,2]. In this example, we have Ny ¼ y22 and by using Eq. (3.11), we get 2 ; A0 ¼ y2;0
A1 ¼ 2y2;0 y2;1 ; 2 A2 ¼ y2;1 þ 2y2;0 y2;2 ;
A3 ¼ 2y2;0 y2;3 þ 2y2;1 y2;2 ; .. . Using the alternate algorithm for computing the Adomian polynomials. The Adomian decomposition series (3.5) leads to the following scheme:
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
y1;0 ¼ 1;
y1;jþ1 ¼ J a ½y1;j ;
y2;0 ¼ 0; y3;0 ¼ 0;
y2;jþ1 ¼ J a ½y1;j Aj ; y3;jþ1 ¼ J a ½Aj :
1357
ð3:12Þ
Using the above recursive relationship and Mathematica and upon using (3.9), the first four terms of the decomposition series are given by ta t2a t3a þ ; Cða þ 1Þ Cð2a þ 1Þ Cð3a þ 1Þ " # ta t2a Cð2a þ 1Þ t3a þ 1 ; y2 ¼ 2 Cða þ 1Þ Cð2a þ 1Þ Cða þ 1Þ Cð3a þ 1Þ
y1 ¼ 1
y3 ¼
t3a þ 2 Cða þ 1Þ Cð3a þ 1Þ
Cð2a þ 1Þ
Setting a ¼ 1, we obtain the solution obtained by Kaya [13] which corresponding to a system of ordinary differential equations. Fig. 1 shows two profile solutions for y1 , the solid line represents the solution for the integer derivative ða ¼ 1Þ, while the dotted line represents the fractional derivative when a ¼ 12. It is easy to conclude that the solution for fractional derivative ð0 < a < 1Þ blows slower than the integer derivative solution as time increases. Example 3.2. Consider the following nonlinear system of fractional differential equations:
y1 t 35 30 25 20 15 10 5 t 2
4
6
8
10
Fig. 1. The approximate solution y1 ðtÞ, for a ¼ 1 (solid line) and a ¼ 12 (dotted line).
1358
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
da y1 ¼ y1 þ y2 y3 ; dta a d y2 ¼ y1 y2 y3 2y22 ; dta da y3 ¼ y22 ; dta
ð3:13Þ
subject to the initial conditions y1 ð0Þ ¼ 1;
y2 ð0Þ ¼ 2;
y3 ð0Þ ¼ 0:
ð3:14Þ
Substituting the initial conditions (3.14) into Eq. (3.5) and using Eq. (3.11) to calculate the Adomian polynomials, yields the following recursive relations: y1;0 ¼ 1;
y1;jþ1 ¼ J a ½y2;j y3;j y1;j ;
y2;0 ¼ 2; y3;0 ¼ 0;
y2;jþ1 ¼ J a ½y1;j y2;j y3;j 2Aj ; y3;jþ1 ¼ J a ½Aj :
Using the above recursive relationship and Mathematica, the first few terms of the decomposition series are given by y1 ¼ 1
ta t2a 28Cð2a þ 1Þ t3a þ ; 2 Cða þ 1Þ Cð2a þ 1Þ Cða þ 1Þ Cð3a þ 1Þ
y2 ¼ 2
7ta 55t2a 28Cð2a þ 1Þ t3a þ ; 2 Cða þ 1Þ Cð2a þ 1Þ Cða þ 1Þ Cð3a þ 1Þ
y3 ¼
4ta 28t2a Cða þ 1Þ Cð2a þ 1Þ
4. Linear differential equation Consider the fractional differential equation of order na of the form: dðnaÞ yðtÞ dððn1ÞaÞ yðtÞ dðaÞ yðtÞ þ a þ þ a þ a0 yðtÞ ¼ f ðtÞ; n1 1 dtna dtðn1Þa dta
ð4:1Þ
with only the integer order initial conditions y ði1Þ ð0Þ ¼ ci1 ;
i ¼ 1; 2; . . . ; n;
ð4:2Þ
where ai and ci are arbitrary constants for i ¼ 0; 1; 2; . . . ; n 1. The numerical approximation of the solution of a linear multi-term fractional differential equation (4.1) and (4.2) can be calculated by reducing the
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
1359
problem to a system of ordinary and fractional equations each of order a, 0 < a 6 1. To convert Eq. (4.1) to a system of equations, we assume y1 ¼ y and ðaÞ ðaÞ ðaÞ y1 ¼ y2 , y2 ¼ y3 ; . . . ; yn1 ¼ yn then ynðaÞ ¼ a0 y1 a1 y2 an1 yn þ f ðtÞ; which defines a linear time-invariant system in the canonical form: 2 ðaÞ 3 32 y 3 2 3 2 y1 0 1 0 ... 0 0 0 1 6 ðaÞ 7 7 607 y 6 y2 7 6 0 0 1 ... 0 0 76 2 7 6 7 76 6 ðaÞ 7 6 y3 7 6 y3 7 6 0 0 0 ... 0 0 76 7 607 6 7 7 6 6 .. 7 .. .. .. 6 .. 7 þ 6 6 ... 7f ðtÞ; 7 6 .. 7 ¼ 6 ... . 76 . ... . . . 7 6 7 6 . 7 6 7 6 6 ðaÞ 7 4 0 0 0 ... 0 1 54 yn1 5 4 0 5 4y 5 n1 a0 a1 a2 . . . an2 an1 1 yn ynðaÞ ð4:3Þ The linear system (4.3) can be expressed in a more compact form as: Y ðaÞ ¼ AY þ Bf ðtÞ:
ð4:4Þ
Operating with J a on both sides of Eq. (4.4) yields Y ¼ J a ½AY þ J a ½Bf ðtÞ:
ð4:5Þ
Following Adomian method analysis, the linear system (4.5) is transformed into a set of recursive relations given by Y0 ¼ C þ J a ½Bf ðtÞ; Yjþ1 ¼ J a ½AYj ; j P 0;
ð4:6Þ
where C is the column vector with ci , i ¼ 0; 1; . . . ; n 1. It is essential feature of the decomposition method that the zeroth component Y0 is defined always by all terms that arise from initial data and from integrating the inhomogeneous terms. Having defined the zeroth component Y0 , the remaining components Yj , j P 1, can be completely determined by using Eq. (4.6). 4.1. Applications Example 4.1. Consider the homogeneous differential equation of order 4a of the form: dð4aÞ yðtÞ dð2aÞ yðtÞ 2 þ yðtÞ ¼ 0; dt4a dt2a yð0Þ ¼ 0; y ð0Þ ¼ 1:
ð4:7Þ
1360
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365 1=2
1=2
1=2
If we choose a ¼ 12, and yðtÞ ¼ y1 , y1 ¼ y2 , y2 ¼ y3 , and y3 (4.7) can be reduced to the following matrix form: 2 ð1=2Þ 3 2 32 3 y1 0 1 0 0 y1 6 ð1=2Þ 7 6 6 y2 7 6 y2 7 6 0 0 1 0 7 76 7: 6 ð1=2Þ 7 ¼ 4 0 0 0 1 54 y3 5 4 y3 5 ð1=2Þ 1 0 2 0 y4 y
¼ y4 , then Eq.
ð4:8Þ
4
Operating with J 1=2 on both sides of the system of equations (4.8) and comparing with Eq. (4.6) we get y1;0 ¼ 0;
y1;jþ1 ¼ J ð1=2Þ ½y2;j ;
y2;0 ¼ 0;
y2;jþ1 ¼ J ð1=2Þ ½y3;j ;
y3;0 ¼ 1;
y3;jþ1 ¼ J ð1=2Þ ½y4;j ;
y4;0 ¼ 0;
y4;jþ1 ¼ J ð1=2Þ ½2y3;j y1;j :
The first few terms of the decomposition series are given by t2 t3 y1 ¼ t 1 þ t þ þ þ ; 2 6 t1=2 2t3=2 3t5=2 þ þ þ ; Cð3=2Þ Cð5=2Þ Cð7=2Þ 3t2 þ ; y3 ¼ 1 þ 2t þ 2 2t1=2 3t3=2 5t5=2 þ þ þ : y4 ¼ Cð3=2Þ Cð5=2Þ Cð7=2Þ y2 ¼
It is easily verified that the exact solution of this problem is yðtÞ ¼ tet : Example 4.2. We next consider the following initial value problem for the inhomogeneous Bagley–Torvik equation [4]: d2 yðtÞ dð3=2Þ yðtÞ þ B þ CyðtÞ ¼ f ðtÞ; dt2 dt3=2 yð0Þ ¼ 1; y ð0Þ ¼ 1;
A
0 6 t 6 5;
ð4:9Þ
pffiffiffiffiffiffi where A ¼ M, B ¼ 2S lq, and C ¼ K. This problem describes the motion of a large plate of the surface S and mass M in a Newtonian fluid of viscosity l and density q. The plate is hanging on a massless spring of stiffness K. The function f ðtÞ represents the loading force. In
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
1361
order to make comparison with the numerical solution of [9], we choose A ¼ B ¼ C ¼ 1 and f ðtÞ ¼ Cðt þ 1Þ. As in the previous example, let yðtÞ ¼ y1 , y11=2 ¼ y2 , y21=2 ¼ y3 , and y31=2 ¼ y4 , then Eq. (4.9) can be reduced to the following matrix form: 2 ð1=2Þ 3 2 32 3 2 3 y1 0 1 0 0 y1 0 6 ð1=2Þ 7 6 6 y2 7 6 0 7 6 y2 7 6 0 0 1 0 7 76 7 þ 6 7ðt þ 1Þ: ð4:10Þ 6 ð1=2Þ 7 ¼ 4 0 0 0 1 54 y3 5 4 0 5 4 y3 5 ð1=2Þ 1 0 0 1 1 y4 y 4
Operating with J 1=2 on both sides of the system of equations (4.10) and comparing with Eq. (4.6) we obtain y1;0 ¼ 1;
y1;jþ1 ¼ J ð1=2Þ ½y2;j ;
y2;0 ¼ 0;
y2;jþ1 ¼ J ð1=2Þ ½y3;j ;
y3;0 ¼ 1;
y3;jþ1 ¼ J ð1=2Þ ½y4;j :
Since the computation depends heavily on y4;0 , we will use here a slight modification [23] for defining the components of y4;j by choosing y4;0 to be the simplest term in (4.10) and add the remaining terms to y4;1 and y4;2 . This will ease the computation considerably. Thus y4;0 ¼ 0; t1=2 þ J ð1=2Þ ½y4;0 y1;0 ; Cð3=2Þ t3=2 þ J ð1=2Þ ½y4;1 y1;1 ; ¼ Cð3=2Þ
y4;1 ¼ y4;2
y4;jþ1 ¼ J ð1=2Þ ½y4;j y1;j ;
j P 2:
The first few terms of the decomposition series are given by y1 ¼ 1 þ t; t1=2 ; Cð3=2Þ t2 t2 y3 ¼ 1 þ ; 2 2 t3=2 t3=2 : y4 ¼ Cð3=2Þ Cð5=2Þ y2 ¼
It is easily verified that the exact solution of this problem is yðtÞ ¼ t þ 1:
1362
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
Table 1 Step size
Decomposition method
Error at t ¼ 5 by [9]
0.5 0.25 0.125 0.0625
0 0 0 0
)0.15131473519232 )0.04684102179946 )0.01602947553912 )0.00562770408881
Table 1 shows the resulting error at t ¼ 5 obtained by the numerical method in [9] and compared with the exact solution obtained by the decomposition method. This demonstrates the importance of our algorithm in solving linear multi-term fractional differential equations. Example 4.3. Consider the fractional differential equations of the form: dyðtÞ da yðtÞ þ a1 þ a2 yðtÞ ¼ f ðtÞ; dt dta yð0Þ ¼ c; 0 < a < 1;
t P 0;
ð4:11Þ
where a1 and a2 are positive constants and c is an arbitrary constant. This fractional differential equation corresponds to the Basset Problem [5]. It represents a classical problem in fluid dynamics where the unsteady motion of a particle accelerates in a viscous fluid due to the gravity of force. In this example, we choose f ðtÞ ¼ 1 and a vanishing initial velocity ðc ¼ 0Þ and consider the case for the ordinary Basset problem, a ¼ 1=2. Eq. (4.11) can be converted into the following system of equations y1 ¼ yðtÞ; ð1=2Þ
y1
¼ y2 ;
ð1=2Þ y2
¼ a2 y1 a1 y2 þ 1:
ð4:12Þ
Operating with J 1=2 on both sides of Eq. (4.12) yields y1 ¼ J 1=2 y2 ; y2 ¼ J 1=2 ð1Þ þ J 1=2 ½a2 y1 a1 y2 :
ð4:13Þ
The Adomian procedures leads to the following scheme: y1;0 ¼ 0;
y1;jþ1 ¼ J ð1=2Þ ½y2;j ;
y2;0 ¼ J ð1=2Þ ð1Þ;
y2;jþ1 ¼ J ð1=2Þ ½a2 y1;j a1 y2;j :
ð4:14Þ
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
1363
The first five terms of the decomposition series are given by: y1 ¼ B0 t þ B1 t3=2 þ B2 t2 þ B3 t5=2 þ B4 t3 þ ; y2 ¼ D0 t1=2 þ D1 t þ D1 t3=2 þ D2 t3=2 þ D3 t2 þ D4 t5=2 þ ;
ð4:15Þ
where 1 ; Cð2Þ a1 B1 ¼ ; Cð5=2Þ B0 ¼
B2 ¼
1 ða2 a2 Þ; Cð3Þ 1
B3 ¼
1 ð2a1 a2 a31 Þ; Cð7=2Þ
B4 ¼
1 ða2 3a2 a1 þ a41 Þ; Cð4Þ 2
1 ; Cð3=2Þ a1 D1 ¼ ; Cð2Þ
D0 ¼
D2 ¼
1 ða2 a2 Þ; Cð5=2Þ 1
D3 ¼
1 ð2a1 a2 a31 Þ; Cð3Þ
D4 ¼
1 ða2 3a2 a1 þ a41 Þ: Cð7=2Þ 2
Using the formulation of the problem given by Mainardi [15], we let a2 ¼
9 1 þ 2v
1=2 ;
where v ¼ 10, 100. Fig. 2 shows the profile solution for y1 . It is easy to conclude that our solution is in agreement with the analytical solution of Mainardi [15] and the numerical solution of Edwards et al. [11] for short times. We should mention that only 5 terms of the decomposition series were used for our approximation solution. It is evident that the solution can be improved by adding more terms of the decomposition series.
1364
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
y t 0.5 0.4 0.3 0.2 0.1 t 0.2
0.4
0.6
0.8
1
Fig. 2. The approximate solution yðtÞ, for v ¼ 10 (solid line) and v ¼ 100 (dotted line).
5. Conclusions The fundamental goal of this work has been to propose an efficient algorithm for the solution of linear and nonlinear systems of fractional differential equations. This goal has been achieved by using Adomian decomposition method and an approximation series solution can be obtained to any desired number of terms. For linear systems of equations the algorithm gives exact solution, and for nonlinear equations it provides an approximate solution with good accuracy. Moreover, we have shown that the use of systems of fractional differential equations of different degrees can provide a convergent method for the solution of multi-term equations. References [1] G. Adomian, A review of the decomposition method in applied mathematics, J. Math. Anal. Appl. 135 (1988) 501–544. [2] G. Adomian, Solving Frontier Problems of Physics: The Decomposition Method, Kluwer Academic Publishers, Boston, 1994. [3] C.T.H. Baker, M.S. Derakhshan, Stability barriers to the construction of q,r-reducible and fractional quadrature rules, in: Numerical Integration III, in: H. Brab, B.G. Hammerlin (Eds.), Ser. Numer. Math., vol. 85, Birkhauser, Basel, 1988, pp. 1–15. [4] R.L. Bagley, P.J. Torvik, On the appearance of the fractional derivative in the behavior of real materials, J. Appl. Mech. 51 (1984) 294–298. [5] A.B. Basset, On the descent of a sphere in a viscous liquid, Quart. J. Math. 42 (1910) 369–381. [6] L. Blank, Numerical treatment of differential equations of fractional order, MCCM Numerical Analysis Report No. 287, The University of Manchester, 1996. [7] M. Caputo, Linear models of dissipation whose Q is almost frequency independent. Part II, J. Roy. Astr. Soc. 13 (1967) 529–539. [8] K. Diethelm, An algorithm for the numerical solution of differential equations of fractional order, Electron. Trans. Numer. Anal. 5 (1997) 1–6.
S. Momani, K. Al-Khaled / Appl. Math. Comput. 162 (2005) 1351–1365
1365
[9] K. Diethelm, N.J. Ford, Numerical solution of the Bagley–Torvik equation, BIT 42 (2002) 490–507. [10] K. Diethelm, Y. Luchko, Numerical solution of linear multi-term differential equations of fractional order, J. Comput. Anal. Appl., in press. [11] J.T. Edwards, N.J. Ford, A.C. Simpson, The numerical solution of linear multi-term fractional differential equations: systems of equations, J. Comput. Appl. Math. 148 (2002) 401–418. [12] T.E. Hull, W.H. Enright, B.M. Fellen, A.E. Sedgwick, Comparing numerical methods for ordinary differential equations, SIAM J. Numer. Anal. 9 (1972) 603. [13] D. Kaya, A reliable method for the numerical solution of the kinetics problems, Appl. Math. Comput, in press. [14] AY. Luchko, R. Groreflo, The initial value problem for some fractional differential equations with the Caputo derivative, Preprint series A08–98, Fachbreich Mathematik und Informatik, Freic Universitat Berlin, 1998. [15] F. Mainardi, Fractional calculus: ‘Some basic problems in continuum and statistical mechanics’, in: A. Carpinteri, F. Mainardi (Eds.), Fractals and Fractional Calculus in Continuum Mechanics, Springer-Verlag, New York, 1997, pp. 291–348. [16] R. Metzler;, J. Klafter, Boundary value problems fractional diffusion equations, Physica A 278 (2000) 107–125. [17] K.S. Miller, B. Ross, An Introduction to the Fractional Calculus and Fractional Differential Equations, John Wiley and Sons, Inc., New York, 1993. [18] S.M. Momani, On the existence of solutions of a system of ordinary differential equations of fractional order, Far East J. Math. Sci. 1 (2) (1999) 265–270. [19] K.B. Oldham, J. Spanier, The Fractional Calculus, Academic Press, New York, 1974. [20] I. Podlubny, Fractional Differential Equations, Academic Press, New York, 1999. [21] N.T. Shawagfeh, Analytical approximate solutions for nonlinear fractional differential equations, Appl. Math. Comput. 131 (2–3) (2002) 517–529. [22] A.M. Wazwaz, Blow-up for solutions of some linear wave equations with mixed nonlinear boundary conditions, Appl. Math. Comput. 123 (2001) 133–140. [23] A.M. Wazwaz, A reliable modification of Adomian’s decomposition method, Appl. Math. Comput. 92 (1998) 1–7.