Efficient solution of a vibration equation involving fractional derivatives

ARTICLE IN PRESS International Journal of Non-Linear Mechanics 45 (2010) 169–175

Contents lists available at ScienceDirect

International Journal of Non-Linear Mechanics journal homepage: www.elsevier.com/locate/nlm

Efficient solution of a vibration equation involving fractional derivatives Attila Pa´lfalvi  + Department of Applied Mechanics, Faculty of Mechanical Engineering, Budapest University of Technology and Economics, 1111 Budapest Muegyetem rkp. 5., Hungary

a r t i c l e in fo

abstract

Article history: Received 26 March 2009 Received in revised form 19 October 2009 Accepted 19 October 2009

Fractional order (or, shortly, fractional) derivatives are used in viscoelasticity since the late 1980s, and they grow more and more popular nowadays. However, their efficient numerical calculation is nontrivial, because, unlike integer-order derivatives, they require evaluation of history integrals in every time step. Several authors tried to overcome this difficulty, either by simplifying these integrals or by avoiding them. In this paper, the Adomian decomposition method is applied on a fractionally damped mechanical oscillator for a sine excitation, and the analytical solution of the problem is found. Also, a series expansion is derived which proves very efficient for calculations of transients of fractional vibration systems. Numerical examples are included. & 2009 Elsevier Ltd. All rights reserved.

Keywords: Vibration Fractional derivative Fractional differential equation Adomian decomposition

1. Introduction It is common to date fractional calculus back to Leibniz, and several textbooks are available on the subject (e.g. [1]). Several authors have been examining the possibility of using fractional derivatives in material modelling in the last two–three decades [2,3], and the field is of a growing interest nowadays [4–9]. Finite element formulations are also in development [9–11]. However, the wide spread of these models is obstructed among others by the difficulty concerning their efficient numerical calculation. This comes from the fact that they need the evaluation of a timehistory integral. There are several solution methods available in the literature. The analytical solution is theoretically known [1]. However, it requires the computation of a convolution of the two-parameter Mittag-Leffler function, which is numerically challenging. Therefore, other calculation methods are sought and used. The most commonly used method is based on the well-known ¨ Grunwald–Letnikov definition of fractional derivatives, used e.g. by Schmidt and Gaul [9]. These authors proceeded to reduce the volume of the time history integral [12], achieving an almost twomagnitude gain in calculation time on an example. Suarez and Shokooh [13] solved a vibration equation with fractional damping of order 12 analytically, taking advantage of the derivative order. Therefore, their solution is restricted to cases where the order of damping is one half. Still, the presented numerical results are widely used as reference solutions. Moreover, Suarez and Shokooh presented a solution for the initial value

 Tel.: + 36 1 463 1369; fax: + 36 1 463 3471.

E-mail address: [email protected] 0020-7462/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijnonlinmec.2009.10.006

problem of the vibration system, which is uncommon in the literature. Yuan and Agrawal [14] have rewritten the definition of a fractional derivative, and turned a fractional differential equation to a system of linear differential equations. However, Schmidt and Gaul [15] have shown that in some cases, this method loses the advantages of fractional calculus over integer-order calculus. Later, this has also been checked by the present author [16] for the field of interest of this paper. The Adomian decomposition method (ADM) has been introduced by George Adomian in the late 1980s. Essentially, it approximates the solution of a non-linear differential equation with a series of functions. The method is getting into use for the solution of fractional differential equations [17–23]. For a special case, Saha Ray et al. [19] have proven that the Adomian solution converges to the analytical solution. Recently, Hu et al. [24] have shown that this is true for any linear fractional differential equation. It will be shown in the present paper that this series expansion is very efficient in many cases. There are some further methods (e.g. the variational iteration method [25,26], the method of Atanackovic and Stankovic [27,28] and the method of Singh and Chatterjee [29]) that will not be treated in this paper. Also, the Adomian decomposition method has been improved by some authors [30–32]. However, the convergence to the analytical solution has only been proven for the original method. The present paper aims to provide a computationally efficient solution method for the fractionally damped vibration equation using the Adomian decomposition method and Taylor series. The obtained solution will be compared to analytical solutions, either known previously or developed here in a series form. To show the ¨ connection with classical methods, solutions using the Grunwald– Letnikov definition will also be calculated.

ARTICLE IN PRESS ´lfalvi / International Journal of Non-Linear Mechanics 45 (2010) 169–175 A. Pa

170

In this paper, Section 2 presents the fractionally damped vibration equation, also giving a brief introduction to fractional derivatives. It also describes the Adomian decomposition method which will be used afterwards. In Section 3, some solutions existing in the literature are shown. Next, in Section 4, new solutions are presented: in Section 4.1, the proposed method is developed, while in Section 4.2, an analytical solution is calculated for reference. Finally, Section 5 shows the numerical examples, and Section 6 concludes the paper.

Another choice is the Caputo definition  n  Z t 1 d n1a C a D xðtÞ ¼ ðt  t Þ xð t Þ dt: a t Gðn  aÞ a dt In both cases, ðn  1Þ o a on:

Actually, the two definitions only differ in the consideration of conditions at the start of the interval: a

a Dt

xðtÞ ¼ Ca Dat xðtÞ þ

2. Preliminaries This section gives the equation to be solved, including a basic introduction to fractional derivatives. The Adomian decomposition method, which is extensively used in this paper, is also presented. 2.1. The equation to be solved The differential equation to be solved is the vibration equation with fractional damping, with one degree of freedom: D2 xðtÞ þ

c a k D xðtÞ þ xðtÞ ¼ f ðtÞ; m m

ð1Þ

where D ¼ d=dt is the differential operator. Another common form of Eq. (1) is a Da xðtÞ þ o2 xðtÞ ¼ f ðtÞ; D2 xðtÞ þ 2Zo2 n n

ð2Þ

with a¼ 2Zo2 n

c m

and

o2n ¼

k : m

2.2. Fractional integrals Fractional integrals and derivatives are deduced from the generalisation of the integer-order operations. It is usual to define the integral operator Dq as Z t 1 q ðt  tÞq1 xðtÞ dt; ð3Þ a Dt xðtÞ ¼ GðqÞ a where q 40 and GðxÞ is the Gamma function Z 1 GðxÞ ¼ ez zx1 dz:

ð4Þ

D

D

xðtÞ ¼ D

ðp þ qÞ

xðtÞ;

ð5Þ

as given in [1] (if both p and q are non-negative). The fractional integral of a polynomial will be needed later so it is given below: Dq t n ¼

Gðn þ 1Þ q þ n t : Gðn þ1 þ qÞ

ð6Þ

With the fractional integral operator, fractional derivatives are easily introduced. 2.3. Fractional derivatives For a real a 4 0, Da is defined by the Riemann–Liouville definition [1], using the above fractional integral operator: a

a Dt

xðtÞ ¼

 n  n Z t d 1 d DtðnaÞ ¼ ðt  tÞn1a xðtÞ dt: dt a Gðn  aÞ dt a

ð9Þ

In the applications, D practically always means 0 Dt , and most authors use the Riemann–Liouville, or the mathematically ¨ equivalent Grunwald–Letnikov definition (see [1] for precise conditions of equivalence). Also, since the Riemann–Liouville definition has a singularity for non-zero initial conditions, the initial conditions are often considered zero. For a physical interpretation of this singularity, see [33]. The composition (5) can be extended to an integral and a derivative: Dp Dq xðtÞ ¼ Dq Dp xðtÞ ¼ Dqp xðtÞ;

ð10Þ

with p; qZ 0 and a continuous xðtÞ. (However, the composition of two fractional order derivatives is not straightforward.) 2.4. Adomian decomposition method The decomposition method has been elaborated by George Adomian in the late 1980s, originally for non-linear differential equations. However, as mentioned in the Introduction, it has also been used by several authors for fractional differential equations [17–23]. As it will be also used in this paper, it is presented in the following. The method has been described very clearly in [34]. The nonlinear differential equation is written as LxðtÞ þ RxðtÞ þ NxðtÞ ¼ f ðtÞ;

ð11Þ

where L is a linear operator which can be inverted easily, R is the remaining linear part and N is a non-linear operator. The function xðtÞ is derived as a series expansion: xðtÞ ¼ x0 ðtÞ þ x1 ðtÞ þ x2 ðtÞ þ    ;

ð12Þ

x0 ðtÞ ¼ ½L1 LxðtÞ  xðtÞ þL1 f ðtÞ

ð13Þ

and

For a continuous xðtÞ, q

n1 X 1 Gðn  aÞ ðt  aÞka xðkÞ ðaÞ: Gðn  aÞ k ¼ 0 Gðk  a þ 1Þ

using

0

p

ð8Þ

ð7Þ

xn þ 1 ðtÞ ¼  L1 Rxn ðtÞ  L1 An ðtÞ;

ð14Þ

with An ðtÞ ¼

" !# 1 X 1 dn i N l x ðtÞ i n! dln i¼0

:

ð15Þ

l¼0

It has been proven recently by Hu et al. [24] that for any linear fractional differential equation, the solution given by the decomposition method converges to the analytical solution. Thus, Adomian decomposition is an adequate tool for solving such equations. In the following, it will be shown that it can be very efficient, too.

3. Known solutions This section presents solutions for the fractionally damped vibration equation (Section 2.1) which are well-established in the literature, and will be used here as reference solutions.

ARTICLE IN PRESS ´lfalvi / International Journal of Non-Linear Mechanics 45 (2010) 169–175 A. Pa

3.1. Analytical solution for a ¼ 12

where

Suarez and Shokooh [13] calculate the analytical solution of Eq. (2) for the special case a ¼ 12, both for a free vibration and a step function excitation. For the former, they obtain

ai ¼ fi 

4 X 1 x ðtÞ ¼ C4j Rj pffiffiffiffiffi þ C4j Rj lj gj ðtÞ; pt j ¼ 1 j¼1 4 X

IC

i1 k c ha X Gðj  aÞ x ; xi  m m GðaÞ j ¼ 0 Gðj þ1Þ ij1

ð24Þ

h being the time step size and fi being the force per unit mass at t ¼ ih.

ð16Þ 4. New solutions

where l2j t

gj ðtÞ ¼ e

pffiffi ð1þ erfðlj t ÞÞ

ð17Þ

and lj and Wj are the solution of the eigenproblem AWj ¼ lj BWj ; with 2

0 6 60 A¼6 61 4 0

0

ð18Þ

1

1

0

0

0

0

0

0

3

7 0 7 7 0 7 5 o2n

2

and

0 6 60 B¼6 60 4 1

0

0

1

0

1

0

1

0

0

0

0

3

3=2 n

2Zo

p

The main goal of this paper is to propose a method for efficient calculation of the numerical solution of Eq. (1) with non-zero initial conditions. In the following, this method will be presented.

0

is the error function, and Rj comes from the initial conditions, by 2 3 2 3 R1 v0 6 7 607 6 R2 7 7 6 7 ¼ ½W1 jW2 jW3 jW4 1 6 ð19Þ 6 7: 6 R3 7 4 x0 5 4 5 R4 0 The response to a step function force (using the Heaviside function HðtÞ) f ðtÞ ¼ f0 HðtÞ;

4 C2 X 4j j¼1

4.1.1. Description In the proposed calculation method, Adomian decomposition is used, separately for initial conditions (starting displacement and velocity) and for the excitation. This can be done, as the differential operator defined in Section 2.3 is linear. The excitation is written in terms of Taylor series, and the resulting method will be referred to as Taylor–Adomian series. For our equation, a possible choice is L ¼ D2 ;

R¼

c a k D þ m m

and

N ¼ 0:

In this case, conforming to Eqs. (13) and (14), x0 ðtÞ ¼ ½D2 D2 xðtÞ  xðtÞ þ D2 f ðtÞ

where X0 and V0 are considered zero, is x ðtÞ ¼ f0

In this section, new solutions will be presented: first, an accurate and efficient series expansion, which is proposed for use; second, an analytical solution using the Adomian decomposition method as reference solution. 4.1. Proposed method: Taylor–Adomian series

7 7 7; 7 5

and C4j is the 4th coordinate of the j th eigenvector. In the above, Z x 2 2 erfðxÞ ¼ pffiffiffiffi et dt

H

171

lj

ðgj ðtÞ  1Þ:

and ð20Þ

For an arbitrary function f ðtÞ and arbitrary initial conditions, the proposal of Suarez and Shokooh is to consider the initial conditions according to the above. As for the force, they suggest discretising the function and handling it as a series of step loads, one in each time step. However, with this method, the time history integral of the fractional derivative is replaced by a sum of displacement responses to step loads, which seems actually slower to calculate.

It is common in the literature, especially in applications, to use ¨ the Grunwald–Letnikov definition of a fractional derivative [1]: 2 3  a NX   1 t Gðj  aÞ t 5 a 4 f tj : ð21Þ D f ðtÞ ¼ lim ! N GðaÞGðj þ1Þ N N1 j¼0 This, along with any time stepping scheme, allows to find a numerical solution for Eq. (1). For simplicity, an explicit scheme has been used for time stepping in this paper:

ai h2 ; 2

k 2 c D xn ðtÞ  Da2 xn ðtÞ; m m

ð26Þ

the latter leading to n ð1Þn X xn ðtÞ ¼ n m j¼0

! n nj j ðð2aÞn þ ajÞ x0 ðtÞ c kD j

ð27Þ

using the composition property (10) of an integral and a derivative operator. In Eq. (25), the first term describes the contribution of initial conditions: ð28Þ

where X0 and V0 are the initial position and velocity, respectively. With Eqs. (27) and (6), one can write ð1Þn ð2aÞn t mn !   n X n X0 V0 t þ : cnj kj t ja  j Gð2n þ1  ðn  jÞaÞ Gð2n þ 2  ðn  jÞaÞ

xIC n ðtÞ ¼

j¼0

ð29Þ

ð22Þ

This result has also been obtained by Baclic and Atanackovic [35] using Laplace transformation. The second term of Eq. (25) gives the response of the system to the excitation:

ð23Þ

xf0 ðtÞ ¼ D2 f ðtÞ:

and xi þ 1 ¼ xi þ vi h þ

xn þ 1 ðtÞ ¼ 

2 2 D xðtÞ  xðtÞ ¼ X0 þ V0 t; xIC 0 ðtÞ ¼ ½D

3.2. Gr¨ unwald–Letnikov formulation

vi þ 1 ¼ vi þ ai h

ð25Þ

ð30Þ

ARTICLE IN PRESS ´lfalvi / International Journal of Non-Linear Mechanics 45 (2010) 169–175 A. Pa

172

Suppose an excitation of the form of the Taylor series 1 X

f ðtÞ ¼

Ti t i ;

ð31Þ

i¼0

f ðtÞ ¼ f0 HðtÞ;

Ti ’s being the coefficients of the polynomial. This, together with Eqs. (30) and (6), leads to xf0;i ðtÞ ¼ Ti

Gðiþ 1Þ i þ 2 t Gðiþ 3Þ

ð32Þ

and xfn;i ðtÞ

n X ð1Þn ¼ Ti Gði þ 1Þti þ ð2aÞn þ 2 mn j¼0

n

!

j

cnj kj t ja ; Gði þ 3 þ 2n  ðn  jÞaÞ

ð33Þ xfn;i ðtÞ

being the n th term of the Adomian series of the response to the i th term of the Taylor series of the excitation. The total response is xf ðtÞ ¼

4.1.3. Examples Step function excitation. For an excitation with the Heaviside function HðtÞ of the form

1 X 1 X

xfn;i ðtÞ:

ð34Þ

n¼0i¼0

solved by Suarez and Shokooh [13] for a ¼ 12, only T0 is non-zero of the Ti ’s, resulting in ! n n ð1Þn ð2aÞn þ 2 X cnj kj xfn ðtÞ ¼ f0 t ja : t ð37Þ n j Gð3 þ2n  ðn  jÞaÞ m j¼0 This has also been obtained by Saha Ray et al. [19] for a ¼ 12 using Adomian decomposition (with a slightly different selection for L, R, and N). Sine excitation. For a sine excitation of the form f ðtÞ ¼ f0 sinðoe tÞ; the Taylor coefficients are 8 if i is even; <0 i Ti ¼ ði1Þ=2 oe : ð1Þ if i is odd: f0 i! Using i ¼ 2hþ 1 gives

Thus, the solution of Eq. (1) has been given.

Th ¼ ð1Þh 4.1.2. Technical considerations The implementation of the method raises several questions. First, this is a double series expansion, which means that both the number of terms in the Taylor series and the number of terms in the Adomian series have to be defined. Moreover, Eqs. (29) and (33) show that terms t ja become large for a relatively small t, which means a need for computing the difference of large numbers. Thus, the evaluation of xðtÞ requires a computational precision better than usual. In the calculation code prepared for this paper, the number of terms is increased dynamically throughout the calculation, based on error estimations from the next ninc terms of the series: if the sum of the next ninc terms is larger than a small limit (‘error goal’), then the number of terms is increased by ninc , and the error estimation will be based on the following ninc terms. This is applied separately for both series. Another point is the required computational effort. The first two terms of the sum in Eq. (29) give

þ1 o2h e

ð2h þ 1Þ!

;

which, inserted into Eq. (33), leads to xfn;h ðtÞ ¼ f0

n ð1Þn þ h 2h þ 1 ð2aÞn þ 3 þ 2h X oe t mn j¼0

n j

!

cnj kj

Gð4 þ 2ðnþ hÞ  ðn  jÞaÞ

t ja ;

ð38Þ also using the property

Gðn þ1Þ ¼ n! of the Gamma function. 4.2. Analytical solution by Adomian decomposition Eq. (27) gives the n th term of the Adomian series of xðtÞ if the original differential equation is Eq. (1). For a sine excitation f ðtÞ ¼ f0 sinðoe tÞ; this yields to

xIC ðtÞ ¼ X0 þ V0 t  

c X0 k X0 2 t 2a  t m Gð3  aÞ m Gð3Þ

c V0 k V0 3 t 3a  t 7  : m Gð4  aÞ m Gð4Þ

xn ðtÞ ¼ ð35Þ

This shows that the exponent of t in the series expansion is either an integer or a real number, the latter being computationally expensive. However, if a is restricted to be rational, as p q

n ð1Þn X n m j¼0

n j

! cnj kj Dðð2aÞn þ ajÞ ðD2 f0 sinðoe tÞÞ:

Using the definition (3) of the integral operator, one obtains Z t 1 ðt  tÞ sinðoe tÞ dt ¼ 2 ðoe t  sinðoe tÞÞ: ð39Þ x0 ðtÞ ¼

oe

0

This leads to

a¼ ; xn ðtÞ ¼

with p and q both integers, the right-hand side of Eq. (35) becomes a polynomial for ðt1=q Þ:

n ð1Þn f0 X n 2 m oe j ¼ 0

n

!

j

cnj kj

1 GðqÞ

Z

t

ðt  tÞq1 ðoe t  sinðoe tÞÞ dt;

0

ð40Þ where

IC

x ðtÞ ¼ X0 þ V0 ðt 

c X0 ðt 1=q Þ2qp Þ  m Gð3  aÞ

1=q q

q ¼ ð2  aÞn þ aj:

k X0 1=q 2q c V0 k V0 1=q 3q ðt Þ  ðt1=q Þ3qp  ðt Þ 7    : m Gð3Þ m Gð4  aÞ m Gð4Þ

ð36Þ This can be computed efficiently.

This results in xn ðtÞ ¼

n ð1Þn f0 X mn o2e j ¼ 0

n j

! cnj kj

1

Gðq þ 2Þ

pffiffi 1=2q t oe s3=2 þ q;1=2 ðoe tÞ; ð41Þ

ARTICLE IN PRESS ´lfalvi / International Journal of Non-Linear Mechanics 45 (2010) 169–175 A. Pa

for n 4 0, where sm;n ðtÞ is the Lommel-function [36]   1 1 1 tm þ 1 sm;n ðtÞ ¼ 1 F2 1; ðm  n þ 3Þ; ðm þ n þ 3Þ;  t 2 ; 2 2 4 ðm þ 1Þ2  n2 with 1 F2 ða1 ; b1 ; b2 ; tÞ being the hypergeometric function 1 F2 ða1 ; b1 ; b2 ; tÞ ¼

1 X

ða1 Þk t k ; ðb 1 Þk ðb2 Þk k! k¼0

where ðaÞk is the Pochhammer symbol ðaÞk ¼

Gða þkÞ : GðaÞ

173

Table 1 Numerical parameters of the examples. Parameters

1

1u

2

2u

3

4

m

1 1=2 10 0.5 1 4p 0.25 0

1 1=2 10 0.5 1 – 0.25 0

1 1=2 10 0.05 1 4p 0.25 0

1 1=2 10 0.05 1 – 0.25 0

1 1=5 10 0.5 1 4p 0.25 0

1 1=5 10 0.05 1 4p 0.25 0

a on Z f0

oe X0 V0

Problems 1u and 2u use a step function excitation.

This can be calculated numerically by any computer algebra system. To the author’s knowledge, this solution has not been calculated in the literature before.

0.25 0.2 0.15

So far, existing (Section 3) and new (Section 4) methods for the solution of the fractionally damped vibration equation have been presented. In this section, they will be tested and compared for some parameter sets.

Displacement

0.1

5. Numerical examples

0.05 0 −0.05 −0.1 −0.15 −0.2

5.1. Software

Problem 1u Problem 2u

−0.25 0

1

2

In the following, numerical examples and their results will be presented and compared. For the latter purpose, let us introduce an absolute error indicator e as N 1X e¼ jxnum  xref i j: Ni¼1 i

Six problems have been analysed, two using a unit step force and four with a harmonic excitation. Table 1 resumes the parameters of the different examples. They lead to steady-state amplitudes between 0.012 and 0.021. The treated time interval is from 0 to 5 s. With the Taylor–Adomian method, 10 000 values of xðtÞ have been calculated for each problem. The floating-point precision was adjusted by hand as the required error level needed it; it was either 128 or 192 bits for all calculations on this time interval.

4

5

Fig. 1. Solution curves for problems 1u and 2u.

0.01 0.0001 1e−06 1e−08 1e−10 1e−12 1e−14 1e−16 Grünwald−Letnikov (1u) Taylor−Adomian (1u) Grünwald−Letnikov (2u) Taylor−Adomian (2u)

1e−18 5.2. Problems and results

3 Time

Error (ε)

First, some details are given on software used for the calculations. The code for the Taylor–Adomian series is written in C. The computational precision is assured by the GNU Multiple precision library [37] and one of its derivatives, MPFR [38]. ¨ The Grunwald–Letnikov formulation is also implemented in C, but it uses the conventional double-precision arithmetics. The analytical solution of Suarez and Shokooh (described in Section 3.1) has been calculated using the Maple computer algebra system. Also, for the analytical solution described in Section 4.2, the Maple computer algebra system has been used with Eqs. (12), (29) and (41). Precision was higher than double precision, usually 50 or 75 decimal digits have been set. The run-time increase of the number of Adomian terms, described for the Taylor–Adomian method, has also been used here, with an error goal of 1025 for 10 terms.

1e−20 1e−22 0.01

0.1

1

10 100 CPU time [s]

1000

10000

Fig. 2. Comparison of methods for problems 1u and 2u (a ¼ 12, step function excitation).

5.2.1. Derivative of order 12, step function excitation (problems 1u and 2u) Here, the reference solution is from Suarez and Shokooh (Section 3.1), and the same numerical parameters have been used as in [13] to reproduce the same results (Fig. 1). Here, and in the following, only the reference solutions are plotted; other solutions are very near when not practically identical to them. Fig. 2 shows the error indicator e and the elapsed CPU time for the two problems. It is immediate to see that, for the same computational effort, the Taylor–Adomian method offers a much

ARTICLE IN PRESS ´lfalvi / International Journal of Non-Linear Mechanics 45 (2010) 169–175 A. Pa

greater precision. (However, an engineering precision is assured by any of the methods.) The difference of the calculation times for the Taylor–Adomian method comes from the difference of numerical conditions (mostly, the dissipation factor Z), which require different numbers of terms in the Adomian series. As one would expect, the ¨ calculation times of the method based on the Grunwald–Letnikov formulation were not affected by the numerical parameters. 5.2.2. Derivative of order 12, harmonic excitation (problems 1 and 2) Here, the reference solution is the solution by Adomian decomposition (Section 4.2). Results are shown in Fig. 3. Error indicators and CPU times are shown in Fig. 4. Observations are the same as above: for a calculation time of the 10-s order, the Taylor–Adomian method is much more precise ¨ than the scheme with the Grunwald–Letnikov derivative.

0.25 0.2 0.15 0.1 Displacement

174

0.05 0 −0.05 −0.1 −0.15 −0.2

Problem 3 Problem 4

−0.25 0

1

2

3

4

5

Time 5.2.3. Derivative of order 15, harmonic excitation (problems 3 and 4) To leave derivatives of order 12, another value has been chosen for a. The reference solution for these problems is the solution by Adomian decomposition (Section 4.2). Results are shown in Fig. 5, while Fig. 6 shows the values of the error indicator and CPU times. It is immediate to see the same as

Fig. 5. Solution of problems 3 and 4.

0.01 0.0001 1e−06 1e−08 Error (ε)

0.25 0.2 0.15

Displacement

0.1

1e−10 1e−12 1e−14 1e−16

0.05

1e−18 0

1e−20 −0.05

1e−22 0.01

−0.1

Grünwald−Letnikov (3) Taylor−Adomian (3) Grünwald−Letnikov (4) Taylor−Adomian (4)

−0.15 −0.2 0

1

2

3

4

5

Time Fig. 3. Solution of problems 1 and 2.

1e−06

Error (ε)

1e−08 1e−10 1e−12 1e−14 1e−16 Grünwald−Letnikov (1) Taylor−Adomian (1) Grünwald−Letnikov (2) Taylor−Adomian (2)

1e−22 0.01

0.1

1

1000

10000

before: the Taylor–Adomian method is clearly much more precise while keeping the calculation time low.

As seen above, the terms of Eqs. (29) and (33) become very large as t increases. It is easy to see that a major limitation of the Taylor–Adomian method will be this property: beyond a certain interval length, a huge number of terms are required, which increases the calculation time horrendously. To check the extents of this phenomenon, calculations have been run with the Taylor–Adomian method up to the end of the transient, but at least to t ¼ 25, with 200 calculated points per time unit. The error goal was 1010 . Fig. 7 shows elapsed CPU time versus simulated time. The points are at the end of the transient period, which is considered to be the end of the last period of excitation for which the error indicator

0.0001

1e−20

10 100 CPU time [s]

5.3. Simulating a longer time interval

0.01

1e−18

1

Fig. 6. Comparison of methods for problems 3 and 4 (a ¼ 15, harmonic excitation).

Problem 1 Problem 2

−0.25

0.1

N 1X jx  xsteady j i Ni¼1 i

10 100 CPU time [s]

1000

10000

Fig. 4. Comparison of methods for problems 1 and 2 (a ¼ 12, harmonic excitation).

is above 1% of the steady-state amplitude. It is immediate to see that problems 1 and 3, i.e. the highly damped equations settle very fast, and their CPU time for a very high precision is in the order of some or some 10 s. Even problem

ARTICLE IN PRESS ´lfalvi / International Journal of Non-Linear Mechanics 45 (2010) 169–175 A. Pa

100000

CPU time [s]

10000

1000

100

Problem 1 Problem 2 Problem 3 Problem 4

10

1 5

10 Simulated time

20

30

40

50

Fig. 7. Calculation times for longer time intervals.

2 settles within a CPU time of less than 10 min on an average 64bit system. Contrarily, the transient of problem 4 required somewhat more than 6 h to calculate on the same computer. This shows that the presented Taylor–Adomian method may be very efficient for the accurate solution of the fractionally damped vibration equation, especially when damping is not very weak.

6. Conclusion In this paper, the main goal was to present a method to solve the fractionally damped, inhomogeneous 1-DOF vibration equation with initial conditions. The Taylor–Adomian method has been proposed and tested on step-function excited and harmonically excited cases. It has proven very efficient in calculating the solution for a reasonably long time interval, practically with an arbitrarily low error. This makes the method immediately usable in providing quick reference results for other methods. The engineering application may be limited, however, due to the necessity of high-precision arithmetics and the fast increase of calculation time as the simulated time interval gets longer. Moreover, the analytical solution of the same equation (with a harmonic excitation) has been calculated by the Adomian decomposition method. This solution served as a reference for the Taylor–Adomian method.

Acknowledgement The author would like to thank Prof. La´szlo´ Szabo´ for his encouragement and most helpful comments. References [1] I. Podlubny, Fractional Differential Equations, in: Mathematics in Science and Engineering, Academic Press, New York, 1999. [2] R.L. Bagley, P.J. Torvik, A theoretical basis for the application of fractional calculus to viscoelasticity, J. Rheol. 27 (1983) 201–210. [3] R.L. Bagley, P.J. Torvik, On the fractional calculus model of viscoelastic behavior, J. Rheol. 30 (1986) 133–155. [4] M. Alvelid, M. Enelund, Modelling of constrained thin rubber layer with emphasis on damping, J. Sound Vib. 300 (2007) 662–675. [5] D.S. Ramrakhyani, G.A. Lesieutre, E.C. Smith, Modeling of elastomeric materials using nonlinear fractional derivative and continuously yielding friction elements, Int. J. Solids Struct. 41 (2004) 3929–3948.

175

[6] A.D. Freed, K. Diethelm, Fractional calculus in biomechanics: a 3D viscoelastic model using regularized fractional derivative kernels with application to the human calcaneal fat pad, Biomech. Model. Mechanobiol. 5 (2006) 203–215. [7] N. Heymans, Fractional calculus description of non-linear viscoelastic behaviour of polymers, Nonlinear Dyn. 38 (2004) 221–231. [8] J. Hernandez-Santiago, A. Macias-Garcia, A. Hernandez-Jimenez, J. SanchezGonzalez, Relaxation modulus in PMMA and PTFE fitting by fractional Maxwell model, Polym. Test. 21 (2002) 325–331. [9] A. Schmidt, L. Gaul, Finite element formulation of viscoelastic constitutive equations using fractional time derivatives, Nonlinear Dyn. 29 (2002) 37–55. [10] F. Corte´s, M.J. Elejabarrieta, Homogenised finite element for transient dynamic analysis of unconstrained layer damping beams involving fractional derivative models, Comput. Mech. 40 (2007) 313–324. [11] M. Enelund, G.A. Lesieutre, Time domain modeling of damping using anelastic displacement fields and fractional calculus, Int. J. Solids Struct. 36 (1999) 4447–4472. [12] A. Schmidt, L. Gaul, On the numerical evaluation of fractional derivatives in multi-degree-of-freedom systems, Signal Process. 86 (2006) 2592–2601. [13] L.E. Suarez, A. Shokooh, An eigenvector expansion method for the solution of motion containing fractional derivatives, J. Appl. Mech. 64 (1997) 629–635. [14] L. Yuan, O.P. Agrawal, A numerical scheme for dynamic systems containing fractional derivatives, J. Vib. Acoust. 124 (2002) 321–324. [15] A. Schmidt, L. Gaul, On a critique of a numerical scheme for the calculation of fractionally damped dynamical systems, Mech. Res. Commun. 33 (2006) 99– 107. [16] A. Pa´lfalvi, Methods for solving a semi-differential vibration equation, Period. Polytech. Ser. Mech. Eng. 51 (2007) 1–5. [17] Q. Wang, Numerical solutions for fractional KdV-Burgers equation by Adomian decomposition method, Appl. Math. Comput. 182 (2006) 1048–1055. [18] V. Daftardar-Gejji, H. Jafari, Adomian decomposition: a tool for solving a system of fractional differential equations, J. Math. Anal. Appl. 301 (2005) 508–518. [19] S. Saha Ray, B.P. Poddar, R.K. Bera, Analytical solution of a dynamic system containing fractional derivative of order one-half by Adomian decomposition method, J. Appl. Mech. 72 (2005) 290–295. [20] V. Daftardar-Gejji, H. Jafari, Solving a multi-order fractional differential equation using Adomian decomposition, Appl. Math. Comput. 189 (2007) 541–548. [21] V. Daftardar-Gejji, S. Bhalekar, Solving multi-term linear and non-linear diffusion-wave equations of fractional order by Adomian decomposition method, Appl. Math. Comput. 202 (2008) 113–120. [22] N.T. Shawagfeh, Analytical approximate solutions for nonlinear fractional differential equations, Appl. Math. Comput. 131 (2002) 517–529. [23] S. Momani, A numerical scheme for the solution of multi-order fractional differential equations, Appl. Math. Comput. 182 (2006) 761–770. [24] Y. Hu, Y. Luo, Z. Lu, Analytical solution of the linear fractional differential equation by Adomian decomposition method, J. Comput. Appl. Math. 215 (2008) 220–229. [25] S. Momani, Z. Odibat, Numerical approach to differential equations of fractional order, J. Comput. Appl. Math. 207 (2007) 96–110. [26] S. Momani, Z. Odibat, Numerical comparison of methods for solving linear differential equations of fractional order, Chaos Solitons Fractals 31 (2007) 1248–1255. [27] T.M. Atanackovic, B. Stankovic, An expansion formula for fractional derivatives and its application, Fractional Calc. Appl. Anal. 7 (2004) 365–378. [28] T.M. Atanackovic, B. Stankovic, On a numerical scheme for solving differential equations of fractional order, Mech. Res. Commun. 35 (2008) 429–438. [29] S.J. Singh, A. Chatterjee, Galerkin projections and finite elements for fractional order derivatives, Nonlinear Dyn. 45 (2006) 183–206. [30] H. Jafari, V. Daftardar-Gejji, Revised Adomian decomposition method for solving a system of nonlinear equations, Appl. Math. Comput. 175 (2006) 1–7. [31] H. Jafari, V. Daftardar-Gejji, Revised Adomian decomposition method for solving systems of ordinary and fractional differential equations, Appl. Math. Comput. 181 (2006) 598–608. [32] S. Saha Ray, K.S. Chaudhuri, R.K. Bera, Application of modified decomposition method for the analytical solution of space fractional diffusion equation, Appl. Math. Comput. 196 (2008) 294–302. [33] N. Heymans, I. Podlubny, Physical interpretation of initial conditions for fractional differential equations with Riemann–Liouville fractional derivatives, Rheol. Acta 45 (2006) 765–771. [34] W. Chen, Z. Lu, An algorithm for Adomian decomposition method, Appl. Math. Comput. 159 (2004) 221–235. [35] B.S. Baclic, T.M. Atanackovic, Stability and creep of a fractional derivative order viscoelastic rod, Bull. Acad. Serb. Sci. Arts CXXI (2000) 115–131. [36] E.W. Weisstein, Lommel function, From MathWorld—A Wolfram Web Resource /http://mathworld.wolfram.com/LommelFunction.htmlS. [37] Free Software Foundation, Inc., GNU Multiple Precision Arithmetic Library, Version 4.2.4, 2008 /http://www.gmplib.org/manual/S. [38] L. Fousse, G. Hanrot, V. Lefe vre, P. Pe´lissier, P. Zimmermann, MPFR: a multiple-precision binary floating-point library with correct rounding, ACM Trans. Math. Softw. 33 (2) (2007) 13:1–13:15.