- Email: [email protected]
An approximate method for numerical solution of fractional differential equations
ARTICLE IN PRESS
Signal Processing 86 (2006) 2602–2610 www.elsevier.com/locate/sigpro
An approximate method for numerical solution of fractional differential equations Pankaj Kumar, Om Prakash Agrawal Mechanical Engineering, Southern Illinois University, Carbondale, IL 62901, USA Received 25 May 2005; received in revised form 15 October 2005; accepted 6 December 2005 Available online 2 March 2006
Abstract This paper presents a numerical solution scheme for a class of fractional differential equations (FDEs). In this approach, the FDEs are expressed in terms of Caputo type fractional derivative. Properties of the Caputo derivative allow one to reduce the FDE into a Volterra type integral equation. Once this is done, a number of numerical schemes developed for Volterra type integral equation can be applied to find numerical solution of FDEs. In this paper the total time is divided into a set of small intervals, and between two successive intervals the unknown functions are approximated using quadratic polynomials. These approximations are substituted into the transformed Volterra type equation to obtain a set of equations. Solution of these equations provides the solution of the FDE. The method is applied to solve two types of FDEs, linear and nonlinear. Results obtained using the scheme presented here agree well with the analytical solutions and the numerical results presented elsewhere. Results also show that the numerical scheme is stable. r 2006 Published by Elsevier B.V. Keywords: Fractional differential equation; Caputo derivative; Fractional-order system
1. Introduction This paper deals with numerical solution of fractional differential equations (FDEs). Fractional derivatives (FDs) and fractional integrals (FIs) have received considerable interest in recent years. In many applications, FDs and FIs provide more accurate models of systems under consideration. For example, FDs have been used successfully to model frequency dependent damping behavior of many viscoelastic materials. Bagley and Torvik Corresponding author. Tel.: +1 618 453 7090;
fax: +1 618 453 7658. E-mail addresses: [email protected] (P. Kumar), [email protected] (O.P. Agrawal). 0165-1684/$ - see front matter r 2006 Published by Elsevier B.V. doi:10.1016/j.sigpro.2006.02.007
[1–3] provided a review of work done in this area prior to 1980, and showed that half-order FD models describe the frequency dependence of the damping materials very well. Other authors have demonstrated applications of FDs and FIs in the areas of electrochemical processes [4,5], dielectric polarization [6], colored noise [7], viscoelastic materials [8–11] and chaos [12]. Mainardi [13] and Rossikhin and Shitikova [14] presented survey of the application of FDs and FIs, in general to solid mechanics, and in particular to modeling of viscoelastic damping. Magin [15–17] presented a three part critical review of applications of fractional calculus in bioengineering. Applications of FDs and FIs in other fields and related mathematical tools and techniques could be found in [15–20].
ARTICLE IN PRESS P. Kumar, O.P. Agrawal / Signal Processing 86 (2006) 2602–2610
Modeling of systems using FDs lead, in most cases, a set of FDEs. Several methods have been proposed to find, in general, the response of these systems, and in particular, the solution of the resulting FDEs. These methods include Laplace transforms [2,3,18–22], Fourier transforms [23,24], modal synthesis and eigenvector expansion [25], numerical methods [26–33], and a method based on the Laguerre integral formula [34]. However, many of these methods cannot be applied to nonlinear FDEs. Further, as pointed out by Diethelm et al. [35], many of these methods are essentially ad hoc methods applicable to specific type of FDEs, and it is not clear if they can be generalized. Further, to the authors knowledge, a systematic convergence analysis of these methods does not appear in the literature. Recently, there has been a growing interest to develop approximate numerical techniques for FDEs which are numerically stable and which can be applied to both linear and nonlinear FDEs. Several papers presented in this area are cited in [35–37]. Many of these techniques take advantage of the fact that the FDEs can be reduced to Volterra type integral equations. Therefore, the numerical schemes for Volterra type integral equations (see [38,39]) can also be applied for the solution of FDEs. Diethelm et al. [35] presented a PECE type method for numerical solution of FDEs, where P, C, and E stands for prediction, correction, and evaluation. Thus, the authors extend the Adams– Bashforth–Moulton type predictor–corrector scheme for ordinary DEs to FDEs. The method presented takes advantage of the fact that the FDEs considered in the paper can be converted into a Volterra-integral equation. These authors also present some error analysis and an extension of Richardson extrapolation to improve the numerical accuracy. The technique presented in [35] has been further analyzed and extended for multi-term and multi-order systems in [36,40–43]. Ford and Simpson [44] presented a numerical scheme for FDEs of order greater than 1. In this formulation, the FDEs of order greater than 1 are reduced to FDEs of order less than 1, and the resulting systems are solved using a numerical scheme. In all these techniques, the unknown functions between the node points are approximated using a linear function. Kumar and Agrawal [45] present numerical scheme for FDEs of order greater than 1. The scheme allows the solution yðtÞ and its derivatives to be continuous at the time node points.
2603
This paper also presents a numerical scheme to find approximate solution of a FDE. The scheme is based on the fact that a class of FDEs can be converted to a Volterra type integral equations. In particular, we use quadratic approximating functions to develop the algorithm, and show that the method can be used to find numerical solution of a class of FDEs. Two examples, linear and nonlinear, are solved to demonstrate the efficiency and accuracy of the solution. It is demonstrated that the numerical scheme is stable. 2. The numerical algorithm Several definitions of a FD have been proposed. These definitions include Riemann–Liouville, Grunwald–Letnikov, Weyl, Caputo, Marchaud, and Riesz FDs [18,20,22,46]. Here, we formulate the problem in terms of the Caputo FD, which is defined as [20], n Z t 1 d Da yðtÞ ¼ ðt tÞna1 yðtÞ dt, Gðn aÞ 0 dt ðn 1oaonÞ,
ð1Þ
where a is the order of the derivative, and n is the smallest integer greater than a. According to Butzer and Westphil [46] and Diethelm et al. [36], Eq. (1) had already been introduced in a 19th century paper of Liouville, and it was used by Rabotnov a year before Caputos paper was published. However, in the literature, the FD defined by Eq. (1) is known as Caputo FD. In the discussion to follow, we consider the following initial value problem containing Caputo FD: find the solution of Da yðtÞ ¼ f ðt; yðtÞÞ
(2)
subject to the following initial conditions: yðkÞ ð0Þ ¼ yðkÞ 0 ;
k ¼ 0; 1; . . . ; n 1,
(3)
where f is an arbitrary function, yðkÞ ðtÞ is the kth derivative of y, and yðkÞ 0 , k ¼ 0; 1; . . . ; n 1 are the specified initial conditions. It is assumed that this function is continuous with respect to both its arguments in the domain of integration, and that it satisfies Lipschitz condition with respect to its second argument. In pure mathematics, Riemann–Liouville (RL) derivative is more commonly used than Caputo derivative. However, Caputo FD is considered here because Riemann–Liouville based FDEs require
ARTICLE IN PRESS P. Kumar, O.P. Agrawal / Signal Processing 86 (2006) 2602–2610
2604
FDs and integrals of y at t ¼ 0. Generally, the physical meanings of these conditions are not known, and for practical applications, they are unavailable. Lorenzo and Hartley [47] have discussed the problem of finding the correct form of the initial conditions in a more general setting (not necessarily assuming that the entire history of the process can be observed [35]). Ref. [48] does the same for FDEs containing Caputo FDs. It is demonstrated in Diethelm and Ford [42] that the problem described by Eqs. (2) and (3) is equivalent to the following Volterra-integral equation: yðtÞ ¼ gðtÞ þ
1 GðaÞ
Z
t
ðt tÞa1 f ðt; yðtÞÞ dt,
(4)
0
where gðtÞ is given as gðtÞ ¼
n1 X k¼0
yðkÞ 0
tk . k!
(5)
In order to explain the quadratic polynomial based numerical method, let us assume that we want to integrate the FDE defined by Eq. (2) from 0 to T. To accomplish this, divide the time T into N equal parts, and let h ¼ T=N be the time interval of each part. The times at the grid points are given as tj ¼ jh, j ¼ 0; 1; . . . ; N. Also assume that the approximate numerical values for yðtÞ have been determined at the grid points tj , j ¼ 0; 1; . . . ; 2m, tj oT. The basic idea of the proposed method is to numerically obtain the values of the function yðtÞ at the next two time points t2mþ1 and t2mþ2 , and repeat the process to advance the integration until the final time T is reached. For simplicity in the discussion to follow, we will use the following notations: yðtj Þ ¼ yðjhÞ yj , gðtj Þ ¼ gðjhÞ gj , and f ðtj ; yðtj ÞÞ ¼ f ðjh; yðjhÞÞ F j . The method developed here will require determination of the integral in Eq. (4) twice at every step. There are two different methods to accomplish this. First, approximate yðtÞ using some approximating functions and determine the integral in Eq. (4) using a numerical technique. This will require initial guesses for y2mþ1 and y2mþ2 without which the integration cannot be performed. Second, approximate yðtÞ and f ðt; yðtÞÞ both using the approximating functions and determine the integral in Eq. (4) explicitly. Note that in this case y2mþ1 and y2mþ2 will appear as coefficients through function f . In this paper, the second approach will be used.
We now provide further details of the algorithm. Let us first determine the value of yðtÞ at t1 and t2 . Using quadratic interpolation functions, yðtÞ and f ðt; yðtÞÞ can be approximated in the interval ½0; t2 as yðtÞ
2 X
f0;k ðtÞyk
(6)
k¼0
and f ðt; yðtÞÞ
2 X
f0;k ðtÞf k ,
(7)
k¼0
where f k f ðtk ; yðtk ÞÞ is the value of the function at the kth time point, fj;k ðtÞ, k ¼ 0; 1; and 2 are the quadratic interpolating functions, and the subscript ðj; kÞ represents the kth approximating function at the step j þ 1, j ¼ 0; . . . ; m. Let us first determine the values of yðtÞ at t1 and t2 . Substituting Eq. (7) into Eq. (4), and integrating, we obtain yðt2 Þ gðt2 Þ þ
2 X
c0;k f k ,
(8)
k¼0
where c0;k
1 ¼ GðaÞ
Z
t2
ðt2 tÞa1 f0;k ðtÞ dt,
(9)
0
which can be exactly computed. Note that Eq. (8) requires the values of f (or indirectly, the values of y) at t1 and t2 . To determine y1 , approximate f ðt; yðtÞÞ on ½t0 ; t1 as X f ðt; yðtÞÞ c0;k ðtÞf k , (10) k¼0;1=2;1
where t1=2 ¼ h=2, f 1=2 ¼ f ðt1=2 ; yðt1=2 ÞÞ, and c0;k ðtÞ, k ¼ 0; 1=2; and 1 are another set of quadratic interpolating functions. Functions f0;k , k ¼ 1, 2, 3, are given as f0;1 ðtÞ ¼
ðt hÞðt 2hÞ , 2h2
f0;2 ðtÞ ¼
ðtÞðt 2hÞ , h2
f0;3 ðtÞ ¼
ðt hÞðtÞ . 2h2
Functions c0;k , k ¼ 1, 2, 3, are defined in a similar way.
ARTICLE IN PRESS P. Kumar, O.P. Agrawal / Signal Processing 86 (2006) 2602–2610
Substituting Eq. (10) into Eq. (4), and integrating, we obtain X yðt1 Þ gðt1 Þ þ d 0;k f k , (11) k¼0;1=2;1
where d 0;k ¼
1 GðaÞ
Z
t1
ðt1 tÞa1 c0;k ðtÞ dt,
(12)
0
which, like c0;k in Eq. (9), can be exactly computed. The value of f 1=2 is determined from Eq. (7), which leads to f 1=2 ¼ 38f 0 þ 34f 1 18f 2 .
(13)
Here, we have explicitly used the properties of the quadratic polynomials. In case of other set of polynomials, the coefficients of f k would be different. Substituting Eq. (13) into Eq. (11), we obtain yðt1 Þ ¼ gðt1 Þ þ d 0 f 0 þ d 1=2 ð38f 0 þ 34f 1 18f 2 Þ þ d 1 f 1 . (14) Note that Eqs. (8) and (14) represent a set of two equations in terms of two unknowns (y1 and y2 ), which can be determined using Newton–Raphson method, fixed point iteration, or another nonlinear solver. In this study, we use Newton–Raphson method to solve these nonlinear equations. This method requires initial guess for y1 and y2 . For a41, the value of the slope at t ¼ 0 could be used to make a better guess for y1 and y2 . However, in this study, for both ao1 and a41, we take the initial guess for these variables as y0 . Note that in each iteration, the time advances by 2h. Let us now assume that the values of y are known at tk , k ¼ 0; 1; . . . ; 2m, and we want to determine y2mþ1 and y2mþ2 . Following the above approach, y2mþ2 and y2mþ1 can be written as Z t2m 1 yðt2mþ2 Þ ¼ gðt2mþ2 Þ þ ðt2mþ2 tÞa1 GðaÞ 0 2mþ2 X f ðt; yðtÞÞ dt þ cm;k f k ð15Þ k¼2m
and Z t2m 1 ðt2mþ1 tÞa1 GðaÞ 0 f ðt; yðtÞÞ dt þ d 2m f 2m þ d 2mþ1=2
yðt2mþ1 Þ ¼ gðt2mþ1 Þ þ
ð38f 2m þ 34f 2mþ1 18f 2mþ2 Þ þ d 2mþ1 f 2mþ1 ,
ð16Þ
2605
where cm;k , k ¼ 2m; 2m þ 1; 2m þ 2 and d m;k , k ¼ 2m; 2m þ 1=2; 2m þ 1 are coefficients which are determined in the same way as c0;k , k ¼ 0; 1; 2 and d 0;k , k ¼ 0; 1=2; 1. Note that the integral in Eqs. (15) and (16) can be determined numerically since the values of yðtÞ at tk , k ¼ 0; 1; . . . ; 2m are known. These equations contain two unknowns, yðt2mþ1 Þ and yðt2mþ2 Þ, which are determined using Newton–Raphson method. In this study, we take y2m as the initial guess for both yðt2mþ1 Þ and yðt2mþ2 Þ. Thus, one can integrate Eq. (1) over the desired interval. As a special case, consider the following nonlinear system: Da yðtÞ ¼ yðtÞ. For this case, f k ¼ yk , and Eqs. (16) and (15) reduce to " #" # " # y2mþ1 a11 a12 b1 ¼ , (17) y a21 a22 b2 2mþ2 where a11 ¼ 1 þ 34d 2mþ1=2 þ d 2mþ1 ,
(18)
a12 ¼ 18d 2mþ1=2 ,
(19)
a21 ¼ c2m;2mþ1 ,
(20)
a22 ¼ 1 þ c2mþ2 ,
(21)
Z t2m 1 ðt2mþ1 tÞa1 yðtÞ dt GðaÞ 0 ðd 2m þ 38d 2mþ1=2 Þy2m ,
b1 ¼ gðt2mþ1 Þ
ð22Þ
and Z t2m 1 b2 ¼ gðt2mþ2 Þ ðt2mþ2 tÞa1 GðaÞ 0 yðtÞ dt cm;2m y2m .
ð23Þ
Eq. (17) is a set of linear simultaneous equation which can be solved using a linear solver. Note the following two additional remarks. First, Eq. (1) contains only one yðtÞ. In case yðtÞ is a vector function, the same algorithm can be applied. However, in that case all yðtÞ and f ðt; yðtÞÞ must be replaced with appropriate vector functions. Second, the scheme requires storage of all past values of y. This is typical of many algorithms for FDEs. This could pose some problems especially when the dimension of y is large as in case of a finite element model of a fractional system. In case, the system has near short term memory, the distance past values of yðtÞ could be neglected to improve the storage space requirements and the computational efficiency.
ARTICLE IN PRESS P. Kumar, O.P. Agrawal / Signal Processing 86 (2006) 2602–2610
2606
3. Numerical results
such that
To demonstrate the effectiveness of this scheme, we consider two examples, one linear and one nonlinear. These examples are considered because closed form solutions are available for them, and they have also been solved using other numerical schemes. This allows one to compare the results obtained using this scheme with the analytical solution and the solutions obtained using other schemes.
yð0Þ ¼ 1;
3.1. Example 1: linear equation As the first example, we consider a linear equation, which is defined as follows: find the numerical solution of the FDE Da yðtÞ ¼ yðtÞ;
0oao2,
(24)
1
y (t)
0.8
0.6
0.4
0.2
0
0
2
4 Time t (sec.)
6
Fig. 1. Comparison of yðtÞ for a ¼ 0:75 (O: Analytical, X: This scheme).
y0 ð0Þ ¼ 0.
(25)
The second initial condition is for a41 only. The solution of Eqs. (24) and (25) is given as yðtÞ ¼ E a ðta Þ,
(26)
where E a ðzÞ ¼
1 X
zk Gðak þ 1Þ k¼0
(27)
is the Mittag–Leffler function of order a. Results are obtained for various values of a and h, some of which are presented here. In each case we take T ¼ 6:4 s. This interval was considered because it is close to the period of the system for a ¼ 2. Fig. 1 shows the comparison of yðtÞ for a ¼ 0:75 obtained analytically and the quadratic numerical scheme developed here. For this case, we take h ¼ 0:1 s. Observe that the two results practically overlap. To highlight the convergence, the numerical results are presented for a ¼ 0:75 and h ¼ 0.4, 0.2, 0.1, 0.05, and 0.025 in Table 1. Note that as the step size is reduced, the error is also reduced, as expected. The ratio of the error R ¼ Eð2hÞ=EðhÞ is for most part very close to 3.37, which indicates that the error is of the order of 1.75 (or EðhÞ ¼ Oðh1:75 Þ). Fig. 2 shows the numerical results for yðtÞ for h ¼ 0:025 and a ¼ 0:25, 0.5, 0.75, 0.95, and 1. The analytical results are not plotted in this figure because they almost coincide with the numerical results. For a ¼ 1, the exact solution is given as yðtÞ ¼ et . Note that as a approaches close 1, the numerical solution converges to the analytical solution yðtÞ ¼ et , i.e. in the limit, the solution of the FDE approaches to that of the integer order DE.
Table 1 Error in function yðtÞ for a ¼ 0:75 and different values of h t#
h ¼ 0:4
h ¼ 0:2
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
0 0.8 1.6 2.4 3.2 4 4.8 5.6 6.4
0 9:832E 4 3:108E 4 2:813E 4 2:100E 4 1:564E 4 1:191E 4 9:287E 5 7:411E 5
0 3:906E 4 2:113E 4 1:243E 4 8:003E 5 5:513E 5 3:998E 5 3:018E 5 2:352E 5
0 1:509E 4 7:031E 5 3:955E 5 2:486E 5 1:687E 5 1:211E 5 9:069E 6 7:029E 6
0 4:730E 5 2:139E 5 1:190E 5 7:433E 6 5:022E 6 3:592E 6 2:685E 6 2:078E 6
0 1:423E 5 6:382E 6 3:537E 6 2:205E 6 1:487E 6 1:063E 6 7:939E 7 6:139E 7
ARTICLE IN PRESS P. Kumar, O.P. Agrawal / Signal Processing 86 (2006) 2602–2610 1
0.8
y (t)
0.6
0.4
0.2
0
0
2
4 Time t (sec.)
6
Fig. 2. Comparison of yðtÞ for different as (O: a ¼ 0:25, X: a ¼ 0:5, þ: a ¼ 0:75, D: a ¼ 0:95, *: a ¼ 1.).
1
y (t)
0.5
0
-0.4 0
2
4 Time t (sec.)
6
Fig. 3. Comparison of yðtÞ for a ¼ 1:5 (O: Analytical, X: This Scheme).
2607
We further present a set of results for a41. The results for ao1 and a41 are presented separately because there is a jump in the slope of yðtÞ at a ¼ 1. Fig. 3 shows the comparison of yðtÞ for a ¼ 1:5 and h ¼ 0:4 obtained analytically and the numerical scheme developed here. Once again the two results practically overlap. To highlight the convergence, Table 2 presents the numerical results for a ¼ 1:5 and h ¼ 0.4, 0.2, 0.1, 0.05, and 0.025. As observed above, for this case also, as the step size is reduced, the error is also reduced. For this case, the ratio of the error R is approximately 5.5 which indicates that the error is of the order of 2.5 (or EðhÞ ¼ Oðh2:5 Þ). Thus, observing the convergence results for ao1 and a41, it can be suggested that the order of error is 1 þ a (or EðhÞ ¼ Oðh1þa Þ), i.e. the order of the error depends not only on h but also on the order of the derivative a as demonstrated in [35]. Fig. 4 shows the numerical results for yðtÞ for h ¼ 0:025 and a ¼ 1:25, 1.5, 1.75, 1.95, and 2. Once again, the analytical results are not plotted in this figure because they almost coincide with the numerical results. For a ¼ 2, the exact solution is given as yðtÞ ¼ cos ðtÞ. Note that once again as a approaches close 2, the numerical solution converges to that of integer order DE. The convergence results shown in Figs. 2 and 4 are important as they demonstrate that the FDEs and their results approach, in the limit, to integer order DEs and their analytical solutions. Table 3 compares the errors in the solution at t ¼ 1:0 reported in [35] and obtained using this scheme for a ¼ 0:5 and 1.25, and for h ¼ 0:1, 0.05 and 0.025. Observe that this scheme gives lower error. This is because the method presented here uses higher order scheme. Similar trend is observed for other values of a and h.
Table 2 Error in function yðtÞ for a ¼ 1:5 different values of h t#
h ¼ 0:4
h ¼ 0:2
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
0 0.8 1.6 2.4 3.2 4 4.8 5.6 6.4
0 1:253E 3 4:999E 4 1:346E 4 3:702E 4 2:915E 4 1:026E 4 4:302E 5 9:386E 5
0 1:664E 4 8:669E 5 6:038E 6 3:438E 5 3:696E 5 2:108E 5 4:275E 6 5:095E 6
0 2:615E 5 1:598E 5 3:683E 6 3:707E 6 5:525E 6 4:008E 6 1:647E 6 6:904E 8
0 4:415E 6 2:920E 6 8:685E 7 4:766E 7 9:064E 7 7:327E 7 3:572E 7 4:760E 8
0 7:671E 7 5:267E 7 1:716E 7 7:061E 8 1:554E 7 1:317E 7 6:824E 8 1:266E 8
ARTICLE IN PRESS P. Kumar, O.P. Agrawal / Signal Processing 86 (2006) 2602–2610
2608
1.8
1
1.6 1.4
0.5
y (t)
y (t)
1.2
0
1 0.8 0.6
-0.5
0.4 Increasing α
0.2
-1
0
2
4 Time t (sec.)
6
Fig. 4. Comparison of yðtÞ for different as (O: a ¼ 1:25, X: a ¼ 1:5, þ: a ¼ 1:75, D: a ¼ 1:95, *: a ¼ 2.).
Table 3 Comparison of error in yðtÞ reported in [35] and obtained using this scheme a ¼ 0:5 h
Ref. [35]
a ¼ 1:5 This scheme
Ref. [35]
This scheme
0.1 1:30E 03 3:98E 04 5:46E 04 2:45E 05 0.05 3:93E 04 1:43E 04 1:28E 04 4:22E 06 0.025 1:26E 04 5:06E 05 3:04E 05 7:40E 07
3.2. Example 2: nonlinear equation As the second example, we consider a nonlinear equation defined as follows: find the numerical solution of Da yðtÞ ¼
40320 8a Gð5 þ a=2Þ 4a=2 t t 3 Gð9 aÞ Gð5 a=2Þ 9 þ Gða þ 1Þ 4 3 3 a=2 4 þ t t ½yðtÞ3=2 2
ð28Þ
such that yð0Þ ¼ 0;
y0 ð0Þ ¼ 0.
(29)
As before, the second initial condition is for a41 only. The exact solution of Eqs. (28) and (29) is given as [35], yðtÞ ¼ t8 3t4þa=2 þ 94ta .
(30)
0 0
0.4 Time t (sec.)
0.8
1
Fig. 5. Comparison of yðtÞ for as ¼ 0:25, 0.75, 1.25, and 1.75 (O: Analytical, X: This Scheme).
Note that for ao1, the slope of the solution at t ¼ 0 goes to infinity. Therefore, one can expect a large numerical error near t ¼ 0. Numerical results are obtained for various values of a and h, some of which are shown below. Fig. 5 shows the analytical and numerical results for h ¼ 0:05 and a ¼ 0:25, 0.75, 1.25, and 1.75. It can be observed that (1) the analytical and numerical results practically overlap. Same observation is made for other values of a, (2) in spite of a very large slope at t ¼ 0, the scheme gives very accurate results, and (3) as expected, at t ¼ 1 the values of yðtÞ for all values of a converges to 0.25 (see Eq. (30)). Table 4 shows the errors in the numerical solutions for a ¼ 0:75 and 1.5, and h ¼ 0:1, 0.05, and 0:025 s. Note that errors decrease as the step size is decreased. The same trend is observed for other values of a and h. For the tried values of a, the error ratios R ¼ Eð2hÞ=EðhÞ suggests no particular order of convergence. However, for a ¼ 0:75 and 1.5, the average values of error of convergence obtained are close to 12 and 15. Table 5 compares the errors in the solution at t ¼ 1:0 reported in [35] and obtained using this scheme for a ¼ 0:25 and 1.25, and h ¼ 0:05. Here we take the values obtained using Richardson extrapolation reported in [35]. Observe that once again this scheme gives lower errors. For a ¼ 0:25, this scheme gives much lower error than Richardson extrapolation scheme does. This may be because for a ¼ 0:25 the slope of yðtÞ near t ¼ 0 changes very rapidly and the linear scheme is unable to capture it accurately. It should be pointed out that this scheme
ARTICLE IN PRESS P. Kumar, O.P. Agrawal / Signal Processing 86 (2006) 2602–2610
2609
Table 4 Error in function yðtÞ for a ¼ 0:75 and 1.5, and for different values of h a ¼ 0:75
a ¼ 1:5
t
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0 1:967E 4 1:374E 4 3:060E 4 2:519E 4 1:117E 4 1:220E 4 7:829E 4 6:911E 4 3:243E 3 2:867E 3
0 6:138E 6 1:246E 5 1:752E 5 1:988E 5 1:692E 5 4:758E 6 2:171E 5 6:861E 5 1:425E 4 2:488E 4
0 5:584E 7 1:037E 6 1:401E 6 1:531E 6 1:213E 6 1:381E 7 2:093E 6 5:953E 6 1:194E 5 2:043E 5
0 1:447E 4 1:058E 4 2:066E 4 1:594E 4 3:803E 5 7:468E 5 6:058E 4 2:928E 4 2:054E 3 1:161E 3
0 4:013E 6 7:124E 6 9:188E 6 9:576E 6 7:308E 6 1:152E 6 1:033E 5 2:879E 5 5:630E 5 9:566E 5
0 2:719E 7 4:651E 7 5:838E 7 5:837E 7 3:968E 7 6:166E 8 8:899E 7 2:201E 6 4:134E 6 6:872E 6
Table 5 Comparison of error in yðtÞ reported in [35] and obtained using this scheme
h
a ¼ 0:25
a ¼ 1:25
Ref. [35]
This scheme Ref. [35]
0.1 2:50E 01 2:51E 03 0.05 1:50E 01 3:47E 04 0.025 4:09E 02 4:26E 05
the step size is reduced. In the limit, as a approaches and integer value, the scheme provides solution for the integer order system. Results also suggest that the scheme is numerically stable.
This scheme
5:53E 03 1:14E 04 2:80E 04 3:14E 05 1:63E 05 3:12E 06
gives more accurate results also for other values of a and h. 4. Conclusions A scheme for approximate numerical solution of a class of fractional differential equations (FDEs) was presented. The FDEs were expressed in terms of Caputo type fractional derivative and the properties of this derivative were used to reduce the FDE into a Volterra type integral equation. The time interval was divided into a set of grid points, and the unknown and known functions yðtÞ and f ðt; yðtÞÞ were approximated over three successive node points using quadratic polynomials. These polynomials were substituted into the Volterra type equations to obtain a set of algebraic equations, and schemes were discussed to solve these equations and to obtain the solution for yðtÞ. Two examples, one linear and the other nonlinear, were solved to demonstrate the performance of the scheme. Results obtained using this scheme agrees with analytical solutions and the numerical results obtained using other schemes. It is shown that results converge as
References [1] R.L. Bagley, P.J. Torvik, A theoretical basis for the application of fractional calculus to viscoelasticity, J. Rheology 27 (3) (1983) 201–210. [2] R.L. Bagley, P.J. Torvik, Fractional calculus—a different approach to the analysis of viscoelastically damped structures, AIAA J. 21 (5) (1983) 741–748. [3] R.L. Bagley, P.J. Torvik, Fractional calculus in the transient analysis of viscoelastically damped structures, AIAA J. 23 (6) (1985) 918–925. [4] M. Ichise, Y. Nagayanagi, T. Kojima, An analog simulation of non-integer order transfer functions for analysis of electrode processes, J. Electronical. Chem. Interfacial Electrochem. 33 (1971) 253–265. [5] H.H. Sun, B. Onaral, Y. Tsao, Application of positive reality principle to metal electrode linear polarization phenomena, IEEE Trans. Biomed. Eng. BME-31 (10) (October 1984) 664–674. [6] H.H. Sun, A.A. Abdelwahab, B. Onaral, Linear approximation of transfer function with a pole of fractional order, IEEE Trans. Automat. Control AC-29 (5) (1984) 441–444. [7] B. Mandelbrot, Some noises with 1/f spectrum, a bridge between direct current and white noise, IEEE Trans. Inform. Theory 13 (2) (April 1967) 289–298. [8] R.L. Bagley, R.A. Calico, Fractional order state equations for the control of viscoelastic structures, J. Guid. Control Dynam. 14 (2) (March–April 1991) 304–311. [9] R.C. Koeller, Application of fractional calculus to he theory of viscoelasticity, J. Appl. Mech. 51 (June 1984) 299–307. [10] R.C. Koeller, Polynomial operators. Stieltjes convolution and fractional calculus in hereditary mechanics, Acta Mech. 58 (1986) 251–264.
ARTICLE IN PRESS 2610
P. Kumar, O.P. Agrawal / Signal Processing 86 (2006) 2602–2610
[11] S.B. Skaar, A.N. Michel, R.K. Miller, Stability of viscoelastic control systems, IEEE Trans. Automat. Control AC33 (4) (April 1988) 348–357. [12] T.T. Hartley, C.F. Lorenzo, H.K. Qammar, Chaoes in a fractional order chua system, IEEE Trans. Circuits Systems I 42 (8) (August 1995) 485–490. [13] F. Mainardi, Fractional calculus: some basic problem in continuum and statistical mechanics, in: A. Carpinteri, F. Mainardi (Eds.), Fractals and Fractional Calculus in Continuum Mechanics, Springer, Wein, New York, 1997, pp. 291–348. [14] Y.A. Rossikhin, M.V. Shitikova, Applications of fractional calculus to dynamic problems of linear and nonlinear hereditary mechanics of solids, Appl. Mech. Rev. 50 (1997) 15–67. [15] R.L. Magin, Fractional calculus in bioengineering, Crit. Rev. Biomed. Eng. 32 (1) (2004) 1–104. [16] R.L. Magin, Fractional calculus in bioengineering—part 2, Crit. Rev. Biomed. Eng. 32 (2) (2004) 105–193. [17] R.L. Magin, Fractional calculus in bioengineering—part 3, Crit. Rev. Biomed. Eng. 32 (3/4) (2004) 194–377. [18] K.B. Oldham, J. Spanier, The Fractional Calculus, Academic Press, New York, 1974. [19] S.G. Samko, A.A. Kilbas, O.I. Marichev, Fractional Integrals and Derivatives—Theory and Applications, Gordon and Breach Science Publishers, Longhorne, PA, 1993. [20] I. Podlubny, Fractional Differential Equations, Academic Press, New York, 1999. [21] R. Gorenflo, F. Mainardi, Fractional calculus: integral and differential equations of fractional order, in: A. Carpinteri, F. Mainardi (Eds.), Fractals and Fractional Calculus in Continuum Mechanics, Springer, Wein, New York, 1997, pp. 223–276. [22] K.S. Miller, B. Ross, An Introduction to the Fractional Calculus and Fractional Differential Equations, Wiley, New York, 1993. [23] L. Gaul, P. Klein, S. Kempfle, Impulse response function of an oscillator with fractional derivative in damping description, Mech. Res. Comm. 16 (5) (1989) 4447–4472. [24] L. Gaul, P. Klein, S. Kempfle, Damping description involving fractional operators, Mech. Systems Signal Process. 5 (2) (1991) 8–88. [25] L.E. Suarez, A. Shokooh, An eigenvector expansion method for the solution of motion containing fractional derivatives, ASME. J. Appl. Mech. 64 (1997) 629–635. [26] M. Enelund, B.L. Josefson, Time domain finite element analysis of viscoelastic structures with fractional derivatives constitutive relations, AIAA J. 35 (10) (1997) 1630–1637. [27] M. Enelund, G.A. Lesieutre, Damping described by fading memory—analysis and application to fractional derivative models, Internat. J. Solid Structures 36 (29) (1999) 939–970. [28] M. Enelund, P. Olsson, Time domain modeling of damping using anelastic displacement fields and fractional calculus, Internat. J. Solids Structures 36 (6) (1999) 939–970. [29] R. Gorenflo, Fractional calculus: some numerical methods, in: A. Carpinteri, F. Mainardi (Eds.), Fractals and Fractional Calculus in Continuum Mechanics, Springer, Wein, New York, 1997, pp. 277–290. [30] C.G. Koh, J.M. Kelly, Application of fractional derivatives to seismic analysis of base-isolated models, Earthquake Eng. Structural Dynamics 19 (1990) 229–241.
[31] J. Padovan, Computational algorithms and finite element formulation involving fractional operators, Comput. Mech. 2 (1987) 271–287. [32] P. Ruge, N. Wagner, Time-domain solutions for vibration systems with feding memory, European Conference of Computational Mechanics, Munchen, Germany, August 31–September 3, 1999. [33] A. Shokooh, L.E. Suarez, A comparison of numerical methods in applied to a fractional model of damping, J. Vibration Controls 5 (1999) 331–354. [34] L. Yuan, O.P. Agrawal, A numerical scheme for dynamic systems containing fractional derivatives, Trans. ASME. J. Vibration and Acoustics 124 (2002) 321–324. [35] K. Diethelm, N.J. Ford, A.D. Freed, A predictor–corrector approach for the numerical solution of fractional differential equations, Nonlinear Dynamics 29 (1–4) (2002) 3–22. [36] K. Diethelm, N.J. Ford, A.D. Freed, Y. Luchko, Algorithms for the fractional calculus: a selection of numerical methods, Comput. Methods Appl. Mech. Eng. 194 (2005) 743–773. [37] P. Kumar, O.P. Agrawal, A cubic scheme for numerical solution of fractional differential equations, in: Proceedings of the Fifth EUROMECH Nonlinear Dynamics Conference, ENOC-2005, pp. 1479–1487. [38] P. Linz, Analytical and Numerical Methods for Volterra Equations, SIAM, Philadelphia, 1985. [39] S. Shaw, J.R. Whiteman, Applications and numerical analysis of partial differential Volterra equations: a brief survey, Comput. Method Appl. Mech. Eng. 150 (1997) 397–409. [40] K. Diethelm, Efficient solution of multi-term fractional differential equation using PðECÞm E methods, Computing 71 (2003) 305–319. [41] K. Diethelm, N.J. Ford, Multi-order fractional differential equations and their numerical solution, Appl. Math. Comput. 154 (3) (2004) 621–640. [42] K. Diethelm, N.J. Ford, Analysis of fractional differential equations, J. Math. Anal. Appl. 265 (2002) 229–248. [43] J.T. Edwards, N.J. Ford, C.A. Simpson, The numerical solution of linear multi-term fractional differential equations: systems of equations, J. Comput. Appl. Math. 148 (2) (2002) 401–418. [44] N.J. Ford, C.A. Simpson, The approximate solution of fractional differential equations of order greater than 1, Internet download, URL hhttp://www.chester.ac.uk/maths/ nevillepub.htmli May 2003. [45] P. Kumar, O.P. Agrawal, Numerical scheme for the solution of fractional differential equations, in: Proceedings of the 2005 ASME Design Engineering Technical Conferences and Computer and Information Engineering Conference, Long Beach, California, September 24–28, 2005. [46] P.L. Butzer, U. Westphil, An introduction to fractional calculus, in: R. Hilfer (Ed.), Applications of Fractional Calculus in Physics, World Scientific, New Jersey, 2000, pp. 1–85. [47] C.F. Lorenzo, T.T. Hartley, Initialized fractional calculus, Internat. J. Appl. Math. 3 (2000) 249–265. [48] B.N. Achar, C.F. Lorenzo, T.T. Hartley, Initialization issue of the Caputo fractional derivative, in: Proceedings of the 2005 ASME Design Engineering Technical Conferences and Computer and Information Engineering Conference, Long Beach, California, September 24–28, 2005.
Signal Processing 86 (2006) 2602–2610 www.elsevier.com/locate/sigpro
An approximate method for numerical solution of fractional differential equations Pankaj Kumar, Om Prakash Agrawal Mechanical Engineering, Southern Illinois University, Carbondale, IL 62901, USA Received 25 May 2005; received in revised form 15 October 2005; accepted 6 December 2005 Available online 2 March 2006
Abstract This paper presents a numerical solution scheme for a class of fractional differential equations (FDEs). In this approach, the FDEs are expressed in terms of Caputo type fractional derivative. Properties of the Caputo derivative allow one to reduce the FDE into a Volterra type integral equation. Once this is done, a number of numerical schemes developed for Volterra type integral equation can be applied to find numerical solution of FDEs. In this paper the total time is divided into a set of small intervals, and between two successive intervals the unknown functions are approximated using quadratic polynomials. These approximations are substituted into the transformed Volterra type equation to obtain a set of equations. Solution of these equations provides the solution of the FDE. The method is applied to solve two types of FDEs, linear and nonlinear. Results obtained using the scheme presented here agree well with the analytical solutions and the numerical results presented elsewhere. Results also show that the numerical scheme is stable. r 2006 Published by Elsevier B.V. Keywords: Fractional differential equation; Caputo derivative; Fractional-order system
1. Introduction This paper deals with numerical solution of fractional differential equations (FDEs). Fractional derivatives (FDs) and fractional integrals (FIs) have received considerable interest in recent years. In many applications, FDs and FIs provide more accurate models of systems under consideration. For example, FDs have been used successfully to model frequency dependent damping behavior of many viscoelastic materials. Bagley and Torvik Corresponding author. Tel.: +1 618 453 7090;
fax: +1 618 453 7658. E-mail addresses: [email protected] (P. Kumar), [email protected] (O.P. Agrawal). 0165-1684/$ - see front matter r 2006 Published by Elsevier B.V. doi:10.1016/j.sigpro.2006.02.007
[1–3] provided a review of work done in this area prior to 1980, and showed that half-order FD models describe the frequency dependence of the damping materials very well. Other authors have demonstrated applications of FDs and FIs in the areas of electrochemical processes [4,5], dielectric polarization [6], colored noise [7], viscoelastic materials [8–11] and chaos [12]. Mainardi [13] and Rossikhin and Shitikova [14] presented survey of the application of FDs and FIs, in general to solid mechanics, and in particular to modeling of viscoelastic damping. Magin [15–17] presented a three part critical review of applications of fractional calculus in bioengineering. Applications of FDs and FIs in other fields and related mathematical tools and techniques could be found in [15–20].
ARTICLE IN PRESS P. Kumar, O.P. Agrawal / Signal Processing 86 (2006) 2602–2610
Modeling of systems using FDs lead, in most cases, a set of FDEs. Several methods have been proposed to find, in general, the response of these systems, and in particular, the solution of the resulting FDEs. These methods include Laplace transforms [2,3,18–22], Fourier transforms [23,24], modal synthesis and eigenvector expansion [25], numerical methods [26–33], and a method based on the Laguerre integral formula [34]. However, many of these methods cannot be applied to nonlinear FDEs. Further, as pointed out by Diethelm et al. [35], many of these methods are essentially ad hoc methods applicable to specific type of FDEs, and it is not clear if they can be generalized. Further, to the authors knowledge, a systematic convergence analysis of these methods does not appear in the literature. Recently, there has been a growing interest to develop approximate numerical techniques for FDEs which are numerically stable and which can be applied to both linear and nonlinear FDEs. Several papers presented in this area are cited in [35–37]. Many of these techniques take advantage of the fact that the FDEs can be reduced to Volterra type integral equations. Therefore, the numerical schemes for Volterra type integral equations (see [38,39]) can also be applied for the solution of FDEs. Diethelm et al. [35] presented a PECE type method for numerical solution of FDEs, where P, C, and E stands for prediction, correction, and evaluation. Thus, the authors extend the Adams– Bashforth–Moulton type predictor–corrector scheme for ordinary DEs to FDEs. The method presented takes advantage of the fact that the FDEs considered in the paper can be converted into a Volterra-integral equation. These authors also present some error analysis and an extension of Richardson extrapolation to improve the numerical accuracy. The technique presented in [35] has been further analyzed and extended for multi-term and multi-order systems in [36,40–43]. Ford and Simpson [44] presented a numerical scheme for FDEs of order greater than 1. In this formulation, the FDEs of order greater than 1 are reduced to FDEs of order less than 1, and the resulting systems are solved using a numerical scheme. In all these techniques, the unknown functions between the node points are approximated using a linear function. Kumar and Agrawal [45] present numerical scheme for FDEs of order greater than 1. The scheme allows the solution yðtÞ and its derivatives to be continuous at the time node points.
2603
This paper also presents a numerical scheme to find approximate solution of a FDE. The scheme is based on the fact that a class of FDEs can be converted to a Volterra type integral equations. In particular, we use quadratic approximating functions to develop the algorithm, and show that the method can be used to find numerical solution of a class of FDEs. Two examples, linear and nonlinear, are solved to demonstrate the efficiency and accuracy of the solution. It is demonstrated that the numerical scheme is stable. 2. The numerical algorithm Several definitions of a FD have been proposed. These definitions include Riemann–Liouville, Grunwald–Letnikov, Weyl, Caputo, Marchaud, and Riesz FDs [18,20,22,46]. Here, we formulate the problem in terms of the Caputo FD, which is defined as [20], n Z t 1 d Da yðtÞ ¼ ðt tÞna1 yðtÞ dt, Gðn aÞ 0 dt ðn 1oaonÞ,
ð1Þ
where a is the order of the derivative, and n is the smallest integer greater than a. According to Butzer and Westphil [46] and Diethelm et al. [36], Eq. (1) had already been introduced in a 19th century paper of Liouville, and it was used by Rabotnov a year before Caputos paper was published. However, in the literature, the FD defined by Eq. (1) is known as Caputo FD. In the discussion to follow, we consider the following initial value problem containing Caputo FD: find the solution of Da yðtÞ ¼ f ðt; yðtÞÞ
(2)
subject to the following initial conditions: yðkÞ ð0Þ ¼ yðkÞ 0 ;
k ¼ 0; 1; . . . ; n 1,
(3)
where f is an arbitrary function, yðkÞ ðtÞ is the kth derivative of y, and yðkÞ 0 , k ¼ 0; 1; . . . ; n 1 are the specified initial conditions. It is assumed that this function is continuous with respect to both its arguments in the domain of integration, and that it satisfies Lipschitz condition with respect to its second argument. In pure mathematics, Riemann–Liouville (RL) derivative is more commonly used than Caputo derivative. However, Caputo FD is considered here because Riemann–Liouville based FDEs require
ARTICLE IN PRESS P. Kumar, O.P. Agrawal / Signal Processing 86 (2006) 2602–2610
2604
FDs and integrals of y at t ¼ 0. Generally, the physical meanings of these conditions are not known, and for practical applications, they are unavailable. Lorenzo and Hartley [47] have discussed the problem of finding the correct form of the initial conditions in a more general setting (not necessarily assuming that the entire history of the process can be observed [35]). Ref. [48] does the same for FDEs containing Caputo FDs. It is demonstrated in Diethelm and Ford [42] that the problem described by Eqs. (2) and (3) is equivalent to the following Volterra-integral equation: yðtÞ ¼ gðtÞ þ
1 GðaÞ
Z
t
ðt tÞa1 f ðt; yðtÞÞ dt,
(4)
0
where gðtÞ is given as gðtÞ ¼
n1 X k¼0
yðkÞ 0
tk . k!
(5)
In order to explain the quadratic polynomial based numerical method, let us assume that we want to integrate the FDE defined by Eq. (2) from 0 to T. To accomplish this, divide the time T into N equal parts, and let h ¼ T=N be the time interval of each part. The times at the grid points are given as tj ¼ jh, j ¼ 0; 1; . . . ; N. Also assume that the approximate numerical values for yðtÞ have been determined at the grid points tj , j ¼ 0; 1; . . . ; 2m, tj oT. The basic idea of the proposed method is to numerically obtain the values of the function yðtÞ at the next two time points t2mþ1 and t2mþ2 , and repeat the process to advance the integration until the final time T is reached. For simplicity in the discussion to follow, we will use the following notations: yðtj Þ ¼ yðjhÞ yj , gðtj Þ ¼ gðjhÞ gj , and f ðtj ; yðtj ÞÞ ¼ f ðjh; yðjhÞÞ F j . The method developed here will require determination of the integral in Eq. (4) twice at every step. There are two different methods to accomplish this. First, approximate yðtÞ using some approximating functions and determine the integral in Eq. (4) using a numerical technique. This will require initial guesses for y2mþ1 and y2mþ2 without which the integration cannot be performed. Second, approximate yðtÞ and f ðt; yðtÞÞ both using the approximating functions and determine the integral in Eq. (4) explicitly. Note that in this case y2mþ1 and y2mþ2 will appear as coefficients through function f . In this paper, the second approach will be used.
We now provide further details of the algorithm. Let us first determine the value of yðtÞ at t1 and t2 . Using quadratic interpolation functions, yðtÞ and f ðt; yðtÞÞ can be approximated in the interval ½0; t2 as yðtÞ
2 X
f0;k ðtÞyk
(6)
k¼0
and f ðt; yðtÞÞ
2 X
f0;k ðtÞf k ,
(7)
k¼0
where f k f ðtk ; yðtk ÞÞ is the value of the function at the kth time point, fj;k ðtÞ, k ¼ 0; 1; and 2 are the quadratic interpolating functions, and the subscript ðj; kÞ represents the kth approximating function at the step j þ 1, j ¼ 0; . . . ; m. Let us first determine the values of yðtÞ at t1 and t2 . Substituting Eq. (7) into Eq. (4), and integrating, we obtain yðt2 Þ gðt2 Þ þ
2 X
c0;k f k ,
(8)
k¼0
where c0;k
1 ¼ GðaÞ
Z
t2
ðt2 tÞa1 f0;k ðtÞ dt,
(9)
0
which can be exactly computed. Note that Eq. (8) requires the values of f (or indirectly, the values of y) at t1 and t2 . To determine y1 , approximate f ðt; yðtÞÞ on ½t0 ; t1 as X f ðt; yðtÞÞ c0;k ðtÞf k , (10) k¼0;1=2;1
where t1=2 ¼ h=2, f 1=2 ¼ f ðt1=2 ; yðt1=2 ÞÞ, and c0;k ðtÞ, k ¼ 0; 1=2; and 1 are another set of quadratic interpolating functions. Functions f0;k , k ¼ 1, 2, 3, are given as f0;1 ðtÞ ¼
ðt hÞðt 2hÞ , 2h2
f0;2 ðtÞ ¼
ðtÞðt 2hÞ , h2
f0;3 ðtÞ ¼
ðt hÞðtÞ . 2h2
Functions c0;k , k ¼ 1, 2, 3, are defined in a similar way.
ARTICLE IN PRESS P. Kumar, O.P. Agrawal / Signal Processing 86 (2006) 2602–2610
Substituting Eq. (10) into Eq. (4), and integrating, we obtain X yðt1 Þ gðt1 Þ þ d 0;k f k , (11) k¼0;1=2;1
where d 0;k ¼
1 GðaÞ
Z
t1
ðt1 tÞa1 c0;k ðtÞ dt,
(12)
0
which, like c0;k in Eq. (9), can be exactly computed. The value of f 1=2 is determined from Eq. (7), which leads to f 1=2 ¼ 38f 0 þ 34f 1 18f 2 .
(13)
Here, we have explicitly used the properties of the quadratic polynomials. In case of other set of polynomials, the coefficients of f k would be different. Substituting Eq. (13) into Eq. (11), we obtain yðt1 Þ ¼ gðt1 Þ þ d 0 f 0 þ d 1=2 ð38f 0 þ 34f 1 18f 2 Þ þ d 1 f 1 . (14) Note that Eqs. (8) and (14) represent a set of two equations in terms of two unknowns (y1 and y2 ), which can be determined using Newton–Raphson method, fixed point iteration, or another nonlinear solver. In this study, we use Newton–Raphson method to solve these nonlinear equations. This method requires initial guess for y1 and y2 . For a41, the value of the slope at t ¼ 0 could be used to make a better guess for y1 and y2 . However, in this study, for both ao1 and a41, we take the initial guess for these variables as y0 . Note that in each iteration, the time advances by 2h. Let us now assume that the values of y are known at tk , k ¼ 0; 1; . . . ; 2m, and we want to determine y2mþ1 and y2mþ2 . Following the above approach, y2mþ2 and y2mþ1 can be written as Z t2m 1 yðt2mþ2 Þ ¼ gðt2mþ2 Þ þ ðt2mþ2 tÞa1 GðaÞ 0 2mþ2 X f ðt; yðtÞÞ dt þ cm;k f k ð15Þ k¼2m
and Z t2m 1 ðt2mþ1 tÞa1 GðaÞ 0 f ðt; yðtÞÞ dt þ d 2m f 2m þ d 2mþ1=2
yðt2mþ1 Þ ¼ gðt2mþ1 Þ þ
ð38f 2m þ 34f 2mþ1 18f 2mþ2 Þ þ d 2mþ1 f 2mþ1 ,
ð16Þ
2605
where cm;k , k ¼ 2m; 2m þ 1; 2m þ 2 and d m;k , k ¼ 2m; 2m þ 1=2; 2m þ 1 are coefficients which are determined in the same way as c0;k , k ¼ 0; 1; 2 and d 0;k , k ¼ 0; 1=2; 1. Note that the integral in Eqs. (15) and (16) can be determined numerically since the values of yðtÞ at tk , k ¼ 0; 1; . . . ; 2m are known. These equations contain two unknowns, yðt2mþ1 Þ and yðt2mþ2 Þ, which are determined using Newton–Raphson method. In this study, we take y2m as the initial guess for both yðt2mþ1 Þ and yðt2mþ2 Þ. Thus, one can integrate Eq. (1) over the desired interval. As a special case, consider the following nonlinear system: Da yðtÞ ¼ yðtÞ. For this case, f k ¼ yk , and Eqs. (16) and (15) reduce to " #" # " # y2mþ1 a11 a12 b1 ¼ , (17) y a21 a22 b2 2mþ2 where a11 ¼ 1 þ 34d 2mþ1=2 þ d 2mþ1 ,
(18)
a12 ¼ 18d 2mþ1=2 ,
(19)
a21 ¼ c2m;2mþ1 ,
(20)
a22 ¼ 1 þ c2mþ2 ,
(21)
Z t2m 1 ðt2mþ1 tÞa1 yðtÞ dt GðaÞ 0 ðd 2m þ 38d 2mþ1=2 Þy2m ,
b1 ¼ gðt2mþ1 Þ
ð22Þ
and Z t2m 1 b2 ¼ gðt2mþ2 Þ ðt2mþ2 tÞa1 GðaÞ 0 yðtÞ dt cm;2m y2m .
ð23Þ
Eq. (17) is a set of linear simultaneous equation which can be solved using a linear solver. Note the following two additional remarks. First, Eq. (1) contains only one yðtÞ. In case yðtÞ is a vector function, the same algorithm can be applied. However, in that case all yðtÞ and f ðt; yðtÞÞ must be replaced with appropriate vector functions. Second, the scheme requires storage of all past values of y. This is typical of many algorithms for FDEs. This could pose some problems especially when the dimension of y is large as in case of a finite element model of a fractional system. In case, the system has near short term memory, the distance past values of yðtÞ could be neglected to improve the storage space requirements and the computational efficiency.
ARTICLE IN PRESS P. Kumar, O.P. Agrawal / Signal Processing 86 (2006) 2602–2610
2606
3. Numerical results
such that
To demonstrate the effectiveness of this scheme, we consider two examples, one linear and one nonlinear. These examples are considered because closed form solutions are available for them, and they have also been solved using other numerical schemes. This allows one to compare the results obtained using this scheme with the analytical solution and the solutions obtained using other schemes.
yð0Þ ¼ 1;
3.1. Example 1: linear equation As the first example, we consider a linear equation, which is defined as follows: find the numerical solution of the FDE Da yðtÞ ¼ yðtÞ;
0oao2,
(24)
1
y (t)
0.8
0.6
0.4
0.2
0
0
2
4 Time t (sec.)
6
Fig. 1. Comparison of yðtÞ for a ¼ 0:75 (O: Analytical, X: This scheme).
y0 ð0Þ ¼ 0.
(25)
The second initial condition is for a41 only. The solution of Eqs. (24) and (25) is given as yðtÞ ¼ E a ðta Þ,
(26)
where E a ðzÞ ¼
1 X
zk Gðak þ 1Þ k¼0
(27)
is the Mittag–Leffler function of order a. Results are obtained for various values of a and h, some of which are presented here. In each case we take T ¼ 6:4 s. This interval was considered because it is close to the period of the system for a ¼ 2. Fig. 1 shows the comparison of yðtÞ for a ¼ 0:75 obtained analytically and the quadratic numerical scheme developed here. For this case, we take h ¼ 0:1 s. Observe that the two results practically overlap. To highlight the convergence, the numerical results are presented for a ¼ 0:75 and h ¼ 0.4, 0.2, 0.1, 0.05, and 0.025 in Table 1. Note that as the step size is reduced, the error is also reduced, as expected. The ratio of the error R ¼ Eð2hÞ=EðhÞ is for most part very close to 3.37, which indicates that the error is of the order of 1.75 (or EðhÞ ¼ Oðh1:75 Þ). Fig. 2 shows the numerical results for yðtÞ for h ¼ 0:025 and a ¼ 0:25, 0.5, 0.75, 0.95, and 1. The analytical results are not plotted in this figure because they almost coincide with the numerical results. For a ¼ 1, the exact solution is given as yðtÞ ¼ et . Note that as a approaches close 1, the numerical solution converges to the analytical solution yðtÞ ¼ et , i.e. in the limit, the solution of the FDE approaches to that of the integer order DE.
Table 1 Error in function yðtÞ for a ¼ 0:75 and different values of h t#
h ¼ 0:4
h ¼ 0:2
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
0 0.8 1.6 2.4 3.2 4 4.8 5.6 6.4
0 9:832E 4 3:108E 4 2:813E 4 2:100E 4 1:564E 4 1:191E 4 9:287E 5 7:411E 5
0 3:906E 4 2:113E 4 1:243E 4 8:003E 5 5:513E 5 3:998E 5 3:018E 5 2:352E 5
0 1:509E 4 7:031E 5 3:955E 5 2:486E 5 1:687E 5 1:211E 5 9:069E 6 7:029E 6
0 4:730E 5 2:139E 5 1:190E 5 7:433E 6 5:022E 6 3:592E 6 2:685E 6 2:078E 6
0 1:423E 5 6:382E 6 3:537E 6 2:205E 6 1:487E 6 1:063E 6 7:939E 7 6:139E 7
ARTICLE IN PRESS P. Kumar, O.P. Agrawal / Signal Processing 86 (2006) 2602–2610 1
0.8
y (t)
0.6
0.4
0.2
0
0
2
4 Time t (sec.)
6
Fig. 2. Comparison of yðtÞ for different as (O: a ¼ 0:25, X: a ¼ 0:5, þ: a ¼ 0:75, D: a ¼ 0:95, *: a ¼ 1.).
1
y (t)
0.5
0
-0.4 0
2
4 Time t (sec.)
6
Fig. 3. Comparison of yðtÞ for a ¼ 1:5 (O: Analytical, X: This Scheme).
2607
We further present a set of results for a41. The results for ao1 and a41 are presented separately because there is a jump in the slope of yðtÞ at a ¼ 1. Fig. 3 shows the comparison of yðtÞ for a ¼ 1:5 and h ¼ 0:4 obtained analytically and the numerical scheme developed here. Once again the two results practically overlap. To highlight the convergence, Table 2 presents the numerical results for a ¼ 1:5 and h ¼ 0.4, 0.2, 0.1, 0.05, and 0.025. As observed above, for this case also, as the step size is reduced, the error is also reduced. For this case, the ratio of the error R is approximately 5.5 which indicates that the error is of the order of 2.5 (or EðhÞ ¼ Oðh2:5 Þ). Thus, observing the convergence results for ao1 and a41, it can be suggested that the order of error is 1 þ a (or EðhÞ ¼ Oðh1þa Þ), i.e. the order of the error depends not only on h but also on the order of the derivative a as demonstrated in [35]. Fig. 4 shows the numerical results for yðtÞ for h ¼ 0:025 and a ¼ 1:25, 1.5, 1.75, 1.95, and 2. Once again, the analytical results are not plotted in this figure because they almost coincide with the numerical results. For a ¼ 2, the exact solution is given as yðtÞ ¼ cos ðtÞ. Note that once again as a approaches close 2, the numerical solution converges to that of integer order DE. The convergence results shown in Figs. 2 and 4 are important as they demonstrate that the FDEs and their results approach, in the limit, to integer order DEs and their analytical solutions. Table 3 compares the errors in the solution at t ¼ 1:0 reported in [35] and obtained using this scheme for a ¼ 0:5 and 1.25, and for h ¼ 0:1, 0.05 and 0.025. Observe that this scheme gives lower error. This is because the method presented here uses higher order scheme. Similar trend is observed for other values of a and h.
Table 2 Error in function yðtÞ for a ¼ 1:5 different values of h t#
h ¼ 0:4
h ¼ 0:2
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
0 0.8 1.6 2.4 3.2 4 4.8 5.6 6.4
0 1:253E 3 4:999E 4 1:346E 4 3:702E 4 2:915E 4 1:026E 4 4:302E 5 9:386E 5
0 1:664E 4 8:669E 5 6:038E 6 3:438E 5 3:696E 5 2:108E 5 4:275E 6 5:095E 6
0 2:615E 5 1:598E 5 3:683E 6 3:707E 6 5:525E 6 4:008E 6 1:647E 6 6:904E 8
0 4:415E 6 2:920E 6 8:685E 7 4:766E 7 9:064E 7 7:327E 7 3:572E 7 4:760E 8
0 7:671E 7 5:267E 7 1:716E 7 7:061E 8 1:554E 7 1:317E 7 6:824E 8 1:266E 8
ARTICLE IN PRESS P. Kumar, O.P. Agrawal / Signal Processing 86 (2006) 2602–2610
2608
1.8
1
1.6 1.4
0.5
y (t)
y (t)
1.2
0
1 0.8 0.6
-0.5
0.4 Increasing α
0.2
-1
0
2
4 Time t (sec.)
6
Fig. 4. Comparison of yðtÞ for different as (O: a ¼ 1:25, X: a ¼ 1:5, þ: a ¼ 1:75, D: a ¼ 1:95, *: a ¼ 2.).
Table 3 Comparison of error in yðtÞ reported in [35] and obtained using this scheme a ¼ 0:5 h
Ref. [35]
a ¼ 1:5 This scheme
Ref. [35]
This scheme
0.1 1:30E 03 3:98E 04 5:46E 04 2:45E 05 0.05 3:93E 04 1:43E 04 1:28E 04 4:22E 06 0.025 1:26E 04 5:06E 05 3:04E 05 7:40E 07
3.2. Example 2: nonlinear equation As the second example, we consider a nonlinear equation defined as follows: find the numerical solution of Da yðtÞ ¼
40320 8a Gð5 þ a=2Þ 4a=2 t t 3 Gð9 aÞ Gð5 a=2Þ 9 þ Gða þ 1Þ 4 3 3 a=2 4 þ t t ½yðtÞ3=2 2
ð28Þ
such that yð0Þ ¼ 0;
y0 ð0Þ ¼ 0.
(29)
As before, the second initial condition is for a41 only. The exact solution of Eqs. (28) and (29) is given as [35], yðtÞ ¼ t8 3t4þa=2 þ 94ta .
(30)
0 0
0.4 Time t (sec.)
0.8
1
Fig. 5. Comparison of yðtÞ for as ¼ 0:25, 0.75, 1.25, and 1.75 (O: Analytical, X: This Scheme).
Note that for ao1, the slope of the solution at t ¼ 0 goes to infinity. Therefore, one can expect a large numerical error near t ¼ 0. Numerical results are obtained for various values of a and h, some of which are shown below. Fig. 5 shows the analytical and numerical results for h ¼ 0:05 and a ¼ 0:25, 0.75, 1.25, and 1.75. It can be observed that (1) the analytical and numerical results practically overlap. Same observation is made for other values of a, (2) in spite of a very large slope at t ¼ 0, the scheme gives very accurate results, and (3) as expected, at t ¼ 1 the values of yðtÞ for all values of a converges to 0.25 (see Eq. (30)). Table 4 shows the errors in the numerical solutions for a ¼ 0:75 and 1.5, and h ¼ 0:1, 0.05, and 0:025 s. Note that errors decrease as the step size is decreased. The same trend is observed for other values of a and h. For the tried values of a, the error ratios R ¼ Eð2hÞ=EðhÞ suggests no particular order of convergence. However, for a ¼ 0:75 and 1.5, the average values of error of convergence obtained are close to 12 and 15. Table 5 compares the errors in the solution at t ¼ 1:0 reported in [35] and obtained using this scheme for a ¼ 0:25 and 1.25, and h ¼ 0:05. Here we take the values obtained using Richardson extrapolation reported in [35]. Observe that once again this scheme gives lower errors. For a ¼ 0:25, this scheme gives much lower error than Richardson extrapolation scheme does. This may be because for a ¼ 0:25 the slope of yðtÞ near t ¼ 0 changes very rapidly and the linear scheme is unable to capture it accurately. It should be pointed out that this scheme
ARTICLE IN PRESS P. Kumar, O.P. Agrawal / Signal Processing 86 (2006) 2602–2610
2609
Table 4 Error in function yðtÞ for a ¼ 0:75 and 1.5, and for different values of h a ¼ 0:75
a ¼ 1:5
t
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0 1:967E 4 1:374E 4 3:060E 4 2:519E 4 1:117E 4 1:220E 4 7:829E 4 6:911E 4 3:243E 3 2:867E 3
0 6:138E 6 1:246E 5 1:752E 5 1:988E 5 1:692E 5 4:758E 6 2:171E 5 6:861E 5 1:425E 4 2:488E 4
0 5:584E 7 1:037E 6 1:401E 6 1:531E 6 1:213E 6 1:381E 7 2:093E 6 5:953E 6 1:194E 5 2:043E 5
0 1:447E 4 1:058E 4 2:066E 4 1:594E 4 3:803E 5 7:468E 5 6:058E 4 2:928E 4 2:054E 3 1:161E 3
0 4:013E 6 7:124E 6 9:188E 6 9:576E 6 7:308E 6 1:152E 6 1:033E 5 2:879E 5 5:630E 5 9:566E 5
0 2:719E 7 4:651E 7 5:838E 7 5:837E 7 3:968E 7 6:166E 8 8:899E 7 2:201E 6 4:134E 6 6:872E 6
Table 5 Comparison of error in yðtÞ reported in [35] and obtained using this scheme
h
a ¼ 0:25
a ¼ 1:25
Ref. [35]
This scheme Ref. [35]
0.1 2:50E 01 2:51E 03 0.05 1:50E 01 3:47E 04 0.025 4:09E 02 4:26E 05
the step size is reduced. In the limit, as a approaches and integer value, the scheme provides solution for the integer order system. Results also suggest that the scheme is numerically stable.
This scheme
5:53E 03 1:14E 04 2:80E 04 3:14E 05 1:63E 05 3:12E 06
gives more accurate results also for other values of a and h. 4. Conclusions A scheme for approximate numerical solution of a class of fractional differential equations (FDEs) was presented. The FDEs were expressed in terms of Caputo type fractional derivative and the properties of this derivative were used to reduce the FDE into a Volterra type integral equation. The time interval was divided into a set of grid points, and the unknown and known functions yðtÞ and f ðt; yðtÞÞ were approximated over three successive node points using quadratic polynomials. These polynomials were substituted into the Volterra type equations to obtain a set of algebraic equations, and schemes were discussed to solve these equations and to obtain the solution for yðtÞ. Two examples, one linear and the other nonlinear, were solved to demonstrate the performance of the scheme. Results obtained using this scheme agrees with analytical solutions and the numerical results obtained using other schemes. It is shown that results converge as
References [1] R.L. Bagley, P.J. Torvik, A theoretical basis for the application of fractional calculus to viscoelasticity, J. Rheology 27 (3) (1983) 201–210. [2] R.L. Bagley, P.J. Torvik, Fractional calculus—a different approach to the analysis of viscoelastically damped structures, AIAA J. 21 (5) (1983) 741–748. [3] R.L. Bagley, P.J. Torvik, Fractional calculus in the transient analysis of viscoelastically damped structures, AIAA J. 23 (6) (1985) 918–925. [4] M. Ichise, Y. Nagayanagi, T. Kojima, An analog simulation of non-integer order transfer functions for analysis of electrode processes, J. Electronical. Chem. Interfacial Electrochem. 33 (1971) 253–265. [5] H.H. Sun, B. Onaral, Y. Tsao, Application of positive reality principle to metal electrode linear polarization phenomena, IEEE Trans. Biomed. Eng. BME-31 (10) (October 1984) 664–674. [6] H.H. Sun, A.A. Abdelwahab, B. Onaral, Linear approximation of transfer function with a pole of fractional order, IEEE Trans. Automat. Control AC-29 (5) (1984) 441–444. [7] B. Mandelbrot, Some noises with 1/f spectrum, a bridge between direct current and white noise, IEEE Trans. Inform. Theory 13 (2) (April 1967) 289–298. [8] R.L. Bagley, R.A. Calico, Fractional order state equations for the control of viscoelastic structures, J. Guid. Control Dynam. 14 (2) (March–April 1991) 304–311. [9] R.C. Koeller, Application of fractional calculus to he theory of viscoelasticity, J. Appl. Mech. 51 (June 1984) 299–307. [10] R.C. Koeller, Polynomial operators. Stieltjes convolution and fractional calculus in hereditary mechanics, Acta Mech. 58 (1986) 251–264.
ARTICLE IN PRESS 2610
P. Kumar, O.P. Agrawal / Signal Processing 86 (2006) 2602–2610
[11] S.B. Skaar, A.N. Michel, R.K. Miller, Stability of viscoelastic control systems, IEEE Trans. Automat. Control AC33 (4) (April 1988) 348–357. [12] T.T. Hartley, C.F. Lorenzo, H.K. Qammar, Chaoes in a fractional order chua system, IEEE Trans. Circuits Systems I 42 (8) (August 1995) 485–490. [13] F. Mainardi, Fractional calculus: some basic problem in continuum and statistical mechanics, in: A. Carpinteri, F. Mainardi (Eds.), Fractals and Fractional Calculus in Continuum Mechanics, Springer, Wein, New York, 1997, pp. 291–348. [14] Y.A. Rossikhin, M.V. Shitikova, Applications of fractional calculus to dynamic problems of linear and nonlinear hereditary mechanics of solids, Appl. Mech. Rev. 50 (1997) 15–67. [15] R.L. Magin, Fractional calculus in bioengineering, Crit. Rev. Biomed. Eng. 32 (1) (2004) 1–104. [16] R.L. Magin, Fractional calculus in bioengineering—part 2, Crit. Rev. Biomed. Eng. 32 (2) (2004) 105–193. [17] R.L. Magin, Fractional calculus in bioengineering—part 3, Crit. Rev. Biomed. Eng. 32 (3/4) (2004) 194–377. [18] K.B. Oldham, J. Spanier, The Fractional Calculus, Academic Press, New York, 1974. [19] S.G. Samko, A.A. Kilbas, O.I. Marichev, Fractional Integrals and Derivatives—Theory and Applications, Gordon and Breach Science Publishers, Longhorne, PA, 1993. [20] I. Podlubny, Fractional Differential Equations, Academic Press, New York, 1999. [21] R. Gorenflo, F. Mainardi, Fractional calculus: integral and differential equations of fractional order, in: A. Carpinteri, F. Mainardi (Eds.), Fractals and Fractional Calculus in Continuum Mechanics, Springer, Wein, New York, 1997, pp. 223–276. [22] K.S. Miller, B. Ross, An Introduction to the Fractional Calculus and Fractional Differential Equations, Wiley, New York, 1993. [23] L. Gaul, P. Klein, S. Kempfle, Impulse response function of an oscillator with fractional derivative in damping description, Mech. Res. Comm. 16 (5) (1989) 4447–4472. [24] L. Gaul, P. Klein, S. Kempfle, Damping description involving fractional operators, Mech. Systems Signal Process. 5 (2) (1991) 8–88. [25] L.E. Suarez, A. Shokooh, An eigenvector expansion method for the solution of motion containing fractional derivatives, ASME. J. Appl. Mech. 64 (1997) 629–635. [26] M. Enelund, B.L. Josefson, Time domain finite element analysis of viscoelastic structures with fractional derivatives constitutive relations, AIAA J. 35 (10) (1997) 1630–1637. [27] M. Enelund, G.A. Lesieutre, Damping described by fading memory—analysis and application to fractional derivative models, Internat. J. Solid Structures 36 (29) (1999) 939–970. [28] M. Enelund, P. Olsson, Time domain modeling of damping using anelastic displacement fields and fractional calculus, Internat. J. Solids Structures 36 (6) (1999) 939–970. [29] R. Gorenflo, Fractional calculus: some numerical methods, in: A. Carpinteri, F. Mainardi (Eds.), Fractals and Fractional Calculus in Continuum Mechanics, Springer, Wein, New York, 1997, pp. 277–290. [30] C.G. Koh, J.M. Kelly, Application of fractional derivatives to seismic analysis of base-isolated models, Earthquake Eng. Structural Dynamics 19 (1990) 229–241.
[31] J. Padovan, Computational algorithms and finite element formulation involving fractional operators, Comput. Mech. 2 (1987) 271–287. [32] P. Ruge, N. Wagner, Time-domain solutions for vibration systems with feding memory, European Conference of Computational Mechanics, Munchen, Germany, August 31–September 3, 1999. [33] A. Shokooh, L.E. Suarez, A comparison of numerical methods in applied to a fractional model of damping, J. Vibration Controls 5 (1999) 331–354. [34] L. Yuan, O.P. Agrawal, A numerical scheme for dynamic systems containing fractional derivatives, Trans. ASME. J. Vibration and Acoustics 124 (2002) 321–324. [35] K. Diethelm, N.J. Ford, A.D. Freed, A predictor–corrector approach for the numerical solution of fractional differential equations, Nonlinear Dynamics 29 (1–4) (2002) 3–22. [36] K. Diethelm, N.J. Ford, A.D. Freed, Y. Luchko, Algorithms for the fractional calculus: a selection of numerical methods, Comput. Methods Appl. Mech. Eng. 194 (2005) 743–773. [37] P. Kumar, O.P. Agrawal, A cubic scheme for numerical solution of fractional differential equations, in: Proceedings of the Fifth EUROMECH Nonlinear Dynamics Conference, ENOC-2005, pp. 1479–1487. [38] P. Linz, Analytical and Numerical Methods for Volterra Equations, SIAM, Philadelphia, 1985. [39] S. Shaw, J.R. Whiteman, Applications and numerical analysis of partial differential Volterra equations: a brief survey, Comput. Method Appl. Mech. Eng. 150 (1997) 397–409. [40] K. Diethelm, Efficient solution of multi-term fractional differential equation using PðECÞm E methods, Computing 71 (2003) 305–319. [41] K. Diethelm, N.J. Ford, Multi-order fractional differential equations and their numerical solution, Appl. Math. Comput. 154 (3) (2004) 621–640. [42] K. Diethelm, N.J. Ford, Analysis of fractional differential equations, J. Math. Anal. Appl. 265 (2002) 229–248. [43] J.T. Edwards, N.J. Ford, C.A. Simpson, The numerical solution of linear multi-term fractional differential equations: systems of equations, J. Comput. Appl. Math. 148 (2) (2002) 401–418. [44] N.J. Ford, C.A. Simpson, The approximate solution of fractional differential equations of order greater than 1, Internet download, URL hhttp://www.chester.ac.uk/maths/ nevillepub.htmli May 2003. [45] P. Kumar, O.P. Agrawal, Numerical scheme for the solution of fractional differential equations, in: Proceedings of the 2005 ASME Design Engineering Technical Conferences and Computer and Information Engineering Conference, Long Beach, California, September 24–28, 2005. [46] P.L. Butzer, U. Westphil, An introduction to fractional calculus, in: R. Hilfer (Ed.), Applications of Fractional Calculus in Physics, World Scientific, New Jersey, 2000, pp. 1–85. [47] C.F. Lorenzo, T.T. Hartley, Initialized fractional calculus, Internat. J. Appl. Math. 3 (2000) 249–265. [48] B.N. Achar, C.F. Lorenzo, T.T. Hartley, Initialization issue of the Caputo fractional derivative, in: Proceedings of the 2005 ASME Design Engineering Technical Conferences and Computer and Information Engineering Conference, Long Beach, California, September 24–28, 2005.