Eva Kaslik

Stefan Balint
email: balint@balint.uvt.ro

Existence of oscillatory solutions in longitudinal flight dynamics

Keywords: longitudinal flight, oscillation, periodic solution, coincidence degree theory, continuation theorem, 76G99

By means of coincidence degree theory and Mawhin's continuation theorem, a theoretical proof is given for the existence of oscillatory solutions of the simplified dynamical system which governs the motion around the center of mass in a longitudinal flight with constant forward velocity of a rigid aircraft, when the automatic flight control system is decoupled.

Introduction

Interest in oscillation susceptibility of aircrafts has been generated by the crashes of high-performance fighter airplanes such as YF-22A and B-2, due to oscillations that were not predicted during the aircraft development process [START_REF] Mehra | Bifurcation and limit cycle analysis of nonlinear pilot induced oscillations[END_REF].

Flight quality criteria for oscillation prediction are based on linear analysis and quasi-linear extensions [START_REF]Unified pilot induced oscillation theory[END_REF]. However, these criteria cannot, in general, predict the presence or the absence of oscillations, because of the large variety of nonlinear phenomena that have been identified as factors contributing to oscillations and which are neglected in the linear approach. Sources of these factors include pilot behavioral transitions, actuator rate limiting [START_REF] Kish | A limit flight test investigation of pilot induced oscillation due to rate limiting[END_REF][START_REF] Klyde | Pilot induced oscillation analysis and prediction with actuator rate limiting[END_REF][START_REF] Mehra | Global stability and control analysis of aircraft at high angles of attack[END_REF] and changes in aircraft dynamics caused by transitions in operating conditions [START_REF] Ionita | Delay induced oscillations[END_REF],

gain scheduling and mode switching [START_REF] Shama | Guaranteed properties of gain schedule control for linear parameter varying plants[END_REF]. The analysis of non-linear oscillations

A c c e p t e d m a n u s c r i p t

involves the computation of nonlinear phenomena including Hopf bifurcation that lead sometimes to large changes in the stability of the pilot-vehicle-system [START_REF] Mehra | Bifurcation and limit cycle analysis of nonlinear pilot induced oscillations[END_REF].

More recently, theoretical bifurcation studies have been undertaken for longitudinal flight dynamics, using the elevator deflection and mass of the vehicle as bifurcation parameters [START_REF] Liaw | Analysis of longitudinal flight dynamics: A bifurcationtheoretic approach[END_REF][START_REF] Liaw | Two-parameter bifurcation analysis of longitudinal flight dynamics[END_REF][START_REF] Pereira | On nonlinear dynamics and an optimal control design to a longitudinal flight[END_REF]. The occurrence of saddle-node and Hopf bifurcations has been pointed out in the case of the F-8 aircraft, and it has been emphasized that these bifurcations may result in jump behavior and pitch oscillations of flight dynamics. Moreover, system controllability with respect to the variation of the elevator deflection angle has been discussed in [START_REF] Pereira | On nonlinear dynamics and an optimal control design to a longitudinal flight[END_REF][START_REF] Liaw | Control design for longitudinal flight dynamics[END_REF].

However, these bifurcation studies can only explain locally the appearance of oscillatory behavior (associated with supercritical Hopf bifurcations), and they do not represent a tool for understanding the global nature of longitudinal flight dynamics. More precisely, a supercritical Hopf bifurcation that occurs at the critical value δ e of the elevator deflection, can only explain the appearance of asymptotically stable limit cycles for values of δ e close to the critical value δ e , i.e. for δ e in a neighborhood of the form (δ eε, δ e ) or (δ e , δ e + ε).

Nevertheless, Hopf bifurcations are not the only type of bifurcation phenomena leading to oscillatory behavior. In [START_REF] Balint | Oscillation susceptibility analysis along the path of longitudinal flight equilibriums[END_REF][START_REF] Kaslik | Numerical analysis of the oscillation susceptibility along the path of longitudinal flight equilibria of a reentry vehicle[END_REF], it has been shown that in a longitudinal flight with constant forward velocity, equilibria exist for the AD-MIRE aircraft and the ALFLEX reentry vehicle only if the elevator deflection δ e belongs to a closed and bounded interval J. When the elevator deflection is at the boundary of the interval J, a countable infinity of saddle-node bifurcation points is present. When the elevator deflection exceeds these critical values and is outside the interval J, numerical simulations show that the angle of attack and pitch rate oscillate with the same period, while the pitch angle increases or decreases infinitely. Hence, the orbit of the system is spiraling.

In this paper, by means of coincidence degree theory and Mawhin's continuation theorem, a theoretical proof is given for the existence of oscillatory solutions of the simplified dynamical system which governs the motion around the center of gravity in a longitudinal flight with constant forward velocity of a rigid outside the interval which corresponds to the existence of equilibria.

Preliminary results

The simplified system of differential equations which governs the motion around the center of gravity in a longitudinal flight with constant forward velocity of a rigid aircraft (assuming that the angle of attack α is small), when the automatic flight control system is decoupled, is given by:

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ α = a 11 (α -α 0 ) + q + g V (cos θ -cos θ 0 ) + b 1 (δ e -δ e0 ) q = a 21 (α -α 0 ) + a 22 q + b 2 (δ e -δ e0 ) θ = q (1)
In this system, the state parameters are: angle of attack α, pitch rate q and Euler pitch angle θ. The control parameter is the elevator angle δ e . V is the forward velocity of the vehicle, considered constant, and g is the gravitational acceleration. The constants α 0 , θ 0 and δ e0 represent the trim angle of attack, trim pitch angle and trim elevator angle, respectively. The appendix contains the formulas for all the parameters appearing in system [START_REF] Mehra | Bifurcation and limit cycle analysis of nonlinear pilot induced oscillations[END_REF] and their numerical values corresponding to the ALFLEX reentry vehicle [START_REF] Goto | Bifurcation analysis for the control of a reentry vehicle[END_REF][START_REF] Goto | Bifurcation analysis for the inertial coupling problem of a reentry vehicle[END_REF][START_REF] Goto | Problems of analysis methods for studying nonlinear flight dynamics[END_REF].

A detailed description of how to derive system (1) from the general system of differential equations [START_REF] Etkin | Dynamics of Flight: Stability and Control[END_REF][START_REF] Cook | Flight Dynamics Principles[END_REF] which describes the motion around the center of gravity of a rigid aircraft, with respect to an xyz body-axis system, where xz is the plane of symmetry, has been presented in [START_REF] Balint | Oscillation susceptibility analysis along the path of longitudinal flight equilibriums[END_REF].

The following notations will be used throughout the paper:

τ = -a 11 -a 22 δ = a 11 a 22 -a 21 γ = (a 21 b 1 -a 11 b 2 )(δ e -δ e0 ) -a 21 g V cos θ 0 ε = -a 21 g V

A c c e p t e d m a n u s c r i p t

The following assumption is taken into consideration:

(A) τ > 0 , δ > 0 , τ 2 > 4δ
One can easily verify that this assumption is fulfilled in the case of the numerical data presented in the appendix.

The following proposition addresses the existence of equilibrium states for system [START_REF] Mehra | Bifurcation and limit cycle analysis of nonlinear pilot induced oscillations[END_REF].

Proposition 1. Consider Λ = ga21 V (a11b2-a21b1
) and let J be the closed interval

J = ⎧ ⎨ ⎩ [δ e0 -Λ(1 + cos θ 0 ), δ e0 + Λ(1 -cos θ 0 )] , if Λ > 0 [δ e0 + Λ(1 -cos θ 0 ), δ e0 -Λ(1 + cos θ 0 )] , if Λ < 0 a.
If δ e / ∈ J then, for system (1), there are no equilibria corresponding to δ e .

b. If δ e ∈ J then, for system (1), there exists a countable infinity of equilibria corresponding to δ e , namely, for any n ∈ Z:

α ± n = α 0 - b 2 a 21 (δ e -δ e0 ) , q ± n = 0 , θ ± n = 2nπ ± arccos cos θ 0 + 1 Λ (δ e -δ e0 ) T
When δ e ∈ ∂J, then the equilibria corresponding to δ e are saddle-node bifurcation points.

Proof. The equilibrium states of system (1) are the solutions of the following algebraic system:

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ a 11 (α -α 0 ) + q + g V (cos θ -cos θ 0 ) + b 1 (δ e -δ e0 ) = 0 a 21 (α -α 0 ) + a 22 q + b 2 (δ e -δ e0 ) = 0 q = 0
Hence, (α, q, θ) T is an equilibrium state of system (1) if and only if q = 0 and α and θ satisfy:

⎧ ⎨ ⎩ a 11 (α -α 0 ) + g V (cos θ -cos θ 0 ) + b 1 (δ e -δ e0 ) = 0 a 21 (α -α 0 ) + b 2 (δ e -δ e0 ) = 0

A c c e p t e d m a n u s c r i p t

From the second equation we obtain that α = α 0 -b2 a21 (δ eδ e0 ), and by eliminating α from the first equation, we have In the followings, we will investigate the case when the control parameter δ e is at the boundary of the interval J. If (α, q, θ) T is an equilibrium state corresponding to δ e = δ e0 -Λ(1 + cos θ 0 ), then there exists n ∈ Z such that θ = (2n+1)π. On the other hand, if (α, q, θ) T is an equilibrium state corresponding

cos θ = cos θ 0 + 1 Λ (δ e -δ e0 ) As cos θ ∈ [-1, 1],
to δ e = δ e0 + Λ(1cos θ 0 ), then there exists n ∈ Z such that θ = 2nπ. The jacobian matrix of system (1) at all these equilibrium states is

⎛ ⎜ ⎜ ⎜ ⎝ a 11 1 0 a 21 a 22 0 0 1 0 ⎞ ⎟ ⎟ ⎟ ⎠
Hence, it is clear that the jacobian matrix has one zero eigenvalue, and the other two eigenvalues are the roots of the equation λ 2 + τλ + δ = 0. Therefore, it is easy to see that the equilibria corresponding to δ e ∈ ∂J are saddle-node bifurcation points.

Lemma 1. If (α(t), q(t), θ(t)) T is a solution of system (1), then θ(t) is the solution of the following third order differential equation:

...

θ (t) + τ θ(t) + δ θ(t) = γ -ε cos θ(t) ( 2 )
Viceversa, if θ(t) is a solution of (2) then (α(t), q(t), θ(t)) T , where

α(t) = α 0 + 1 a 21 [ θ(t) -a 22 θ(t) -b 2 (δ e -δ e0 )] and q(t) = θ(t),
is a solution of system [START_REF] Mehra | Bifurcation and limit cycle analysis of nonlinear pilot induced oscillations[END_REF].

A c c e p t e d m a n u s c r i p t

Proof. If (α(t), q(t), θ(t)) T is a solution of system (1), taking into account that θ(t) = q(t), from the second equation of (1) we obtain:

θ(t) = a 21 (α(t) -α 0 ) + a 22 θ(t) + b 2 (δ e -δ e0 )
From here, we can express α with respect to θ as follows:

α(t) = α 0 + 1 a 21 θ(t) -a 22 θ -b 2 (δ e -δ e0 )
From here, computing α(t) and replacing α(t) and α(t) in the first equation of ( 1) and making all the necessary simplifications, leads us to equation ( 2).

Viceversa, if θ(t) is a solution of ( 2), it can be verified by straightforward computations that (α(t), q(t), θ(t)) T , where

α(t) = α 0 + 1 a 21 [ θ(t) -a 22 θ(t) -b 2 (δ e -δ e0 )] and q(t) = θ(t),
is a solution of system (1).

Remark 1. Proposition 1 translates into the following necessary and sufficient condition for the existence of equilibrium states for both (1) and ( 2):

|γ| ≤ |ε|

At γ = ±ε, saddle-node bifurcations occur.

Definition 1. We say that θ(t) is a monotonic oscillatory solution of equation (2) if the derivative θ(t) is a strictly positive or strictly negative periodic function.

In this paper, we will show that if γ > |ε| then equation ( 2) has at least one increasing oscillatory solution. In a similar way, we will show that if γ < -|ε| then equation ( 2) has at least one decreasing oscillatory solution. However, before proceeding to the proof of the main results, we will give some useful preliminary lemmas.

Lemma 2. (a) If there exists an increasing oscillatory solution θ(t) of ( 2), such that T > 0 is the period of θ(t), then there exists n ∈ N such that θ(t + T ) =

A c c e p t e d m a n u s c r i p t

θ(t) + 2nπ and there exists a 2nπ-periodic solution x(s) of the equation Proof. (a) Assume that there exists an increasing oscillatory solution θ(t) of ( 2), such that θ(t) is T -periodic. From (2) we deduce that cos(θ(t + T )) = cos(θ(t)), and hence, there exists n ∈ N such that θ(t + T ) = θ(t) + 2nπ.

x (s) = 2x (s) 2 -τx (s)e x(s) + δe 2x(s) -(γ -ε cos s)e 3x(s) (3)
The function θ : R → R is a strictly increasing bijection. The inverse function ψ = θ -1 is also strictly increasing and therefore, we can define x(s) = ln ψ(s).

As θ(t+T ) = θ(t)+2nπ, making t = ψ(s) in this equation, and then applying

ψ, we obtain ψ(s + 2nπ) = ψ(s) + T . Hence, x(s) is 2nπ-periodic. Moreover,
we can easily see that

2nπ 0 e x(s) ds = ψ(2nπ) -ψ(0) = T .
By successive differentiation, taking into account that θ(ψ(s)) = s, we obtain:

θ(ψ(s)) = 1 ψ(s) = e -x(s) θ(ψ(s)) = - 1 ψ(s) x (s)e -x(s) = -x (s)e -2x(s) ... θ (ψ(s)) = - 1 ψ(s) [x (s) -2x (s) 2 ]e -2x(s) = [2x (s) 2 -x (s)]e -3x(s)
Making t = ψ(s) in (2) and using the equations from above, we finally obtain that x(s) is a solution of (3).

(b) Suppose that equation ( 3) has a 2nπ-periodic solution x(s). Consider θ 0 ∈ R. It can be easily seen that the function ψ : R → R defined by

ψ(s) = s θ 0 e x(u) du
is strictly increasing and bijective. Consider the inverse function θ = ψ -1 : R → R which is also strictly increasing. One has

θ(t) = 1 ψ(θ(t)) = e -x(θ(t)) > 0 A c c e p t e d m a n u s c r i p t θ(t) = -θ(t)x (θ(t))e -x(θ(t)) = -x (θ(t))e -2x(θ(t))
...

θ (t) = -θ(t)x (θ(t))e -2x(θ(t)) + 2 θ(t)x (θ(t)) 2 e -2x(θ(t)) = [2x (θ(t)) 2 -x (θ(t))]e -3x(θ(t))
Hence ...

θ (t) + τ θ(t) + δ θ(t) = e -3x(θ(t)) [2x (θ(t)) 2 -x (θ(t)) -τx (θ(t))e x(θ(t)) + δe 2x(θ(t)) ] = = γ -ε cos θ(t)
so θ(t) is a solution of (2). Moreover, as x(s) is 2nπ-periodic, it follows that

ψ(s + 2nπ) = s+2nπ θ 0 e x(u) du = θ 0 +2nπ θ 0 e x(u) du + s+2nπ θ 0 +2nπ e x(u) du = = 2nπ 0 e x(u) du + s θ 0 e x(u) du = ψ(s) + T where T = 2nπ 0 e x(u) du. Making s = θ(t) in the equality ψ(s+2nπ) = ψ(s)+T ,
we get ψ(θ(t) + 2nπ) = t + T and applying θ to this equality, we obtain:

θ(t + T ) = θ(t) + 2nπ
Therefore, we obtain that θ is a strictly positive T -periodic function and hence, θ(t) is an increasing oscillatory solution of (2).

Therefore, in order to prove that equation ( 2) has an increasing oscillatory solution, it is sufficient to prove that there exists a 2nπ-periodic solution of equation (3).

Remark 2. Denoting x 1 = x and x 2 = τe -xx e -2x , equation ( 3) is replaced by the system:

⎧ ⎨ ⎩ x 1 (s) = τe x1(s) -x 2 (s)e 2x1(s) x 2 (s) = (γ -ε cos s)e x1(s) -δ (4)
Hence, in order to show that equation ( 2) has an increasing oscillatory solution, it is sufficient to study the existence of a 2nπ-periodic solution for the system (4).

A c c e p t e d m a n u s c r i p t

Lemma 3. (a) If there exists a decreasing oscillatory solution θ(t) of ( 2), such that T > 0 is the period of θ(t), then there exists n ∈ N such that θ(t + T ) = θ(t) -2nπ and there exists a 2nπ-periodic solution x(s) of the equation Proof. (a) Assume that there exists a decreasing oscillatory solution θ(t) of ( 2), such that θ(t) is T -periodic. From (2) we deduce that cos(θ(t + T )) = cos(θ(t)), and hence, there exists n ∈ N such that θ(t

x (s) = 2x (s) 2 + τx x(s)e x(s) + δe 2x(s) + (γ -ε cos s)e 3x(s) (5
+ T ) = θ(t) -2nπ.
The function θ : R → R is a strictly decreasing bijection. The inverse function ψ = θ -1 is also strictly decreasing and therefore, we can define x(s) = ln(-ψ(s)).

As θ(t+T ) = θ(t)-2nπ, making t = ψ(s) in this equation, and then applying ψ, we obtain ψ(s -2nπ) = ψ(s) + T . Hence, x(s) is 2nπ-periodic. Moreover, we can easily see that

θ(ψ(s)) = 1 ψ(s) = -e -x(s) θ(ψ(s)) = 1 ψ(s) x (s)e -x(s) = -x (s)e -2x(s) ... θ (ψ(s)) = - 1 ψ(s) [x (s) -2x (s) 2 ]e -2x(s) = [x (s) -2x (s) 2 ]e -3x(s)
Making t = ψ(s) in (2) and using the equations from above, we finally obtain that x(s) is a solution of (5). = -e -x(θ(t)) < 0 θ(t) = θ(t)x (θ(t))e -x(θ(t)) = -x (θ(t))e -2x(θ(t))

...

θ (t) = -θ(t)x (θ(t))e -2x(θ(t)) + 2 θ(t)x (θ(t)) 2 e -2x(θ(t)) = [x (θ(t)) -2x (θ(t)) 2 ]e -3x(θ(t))
Hence ...

θ (t) + τ θ(t) + δ θ(t) = e -3x(θ(t)) [x (θ(t)) -2x (θ(t)) 2 -τx (θ(t))e x(θ(t)) -δe 2x(θ(t)) ] = = γ -ε cos θ(t)
so θ(t) is a solution of (2). Moreover, as x(s) is 2nπ-periodic, it follows that

-ψ(s + 2nπ) =
θ(t + T ) = θ(t) -2nπ
Therefore, we obtain that θ is a negative T -periodic function and hence, θ(t) is a decreasing oscillatory solution of (2).

Therefore, in order to prove that equation (2) has a decreasing oscillatory solution, it is sufficient to prove that there exists a 2nπ-periodic solution of equation ( 5).

Remark 3. Denoting x 1 = x and x 2 = τe -x + x e -2x , equation ( 5) is replaced by the system:

⎧ ⎨ ⎩ x 1 (s) = -τe x1(s) + x 2 (s)e 2x1(s) x 2 (s) = (γ -ε cos s)e x1(s) + δ (6)
Hence, in order to show the existence of a decreasing oscillatory solution of equation [START_REF]Unified pilot induced oscillation theory[END_REF], it is sufficient to study the existence of a 2nπ-periodic solution for the system (6)

A c c e p t e d m a n u s c r i p t

We introduce the following notations:

[f ] M = max s∈[0,2nπ] f (s), [f ] L = min s∈[0,2nπ] f (s),
where f is a continuous 2nπ-periodic function. The following result holds:

Lemma 4. If f is a smooth 2nπ-periodic function then [f ] M ≤ [f ] L + 1 2 2nπ 0 | ḟ (s)|ds
Proof. As f is a smooth 2nπ-periodic function, there exists ζ ∈ [0, 2nπ] and

η ∈ [ζ, ζ + 2nπ] such that [f ] L = f (ζ) and [f ] M = f (η).
On one hand, we have:

f (η) = f (ζ) + η ζ ḟ(s)ds ≤ f (ζ) + η ζ | ḟ (s)|ds
On the other hand:

f (η) = f (ζ + 2nπ) + η ζ+2nπ ḟ(s)ds ≤ f (ζ) + ζ+2nπ η | ḟ (s)|ds
Adding up the above inequalities, we obtain:

2f (η) ≤ 2f (ζ) + ζ+2nπ ζ | ḟ (s)|ds = 2f (ζ) + 2nπ 0 | ḟ(s)|ds which leads to [f ] M ≤ [f ] L + 1 2 2nπ 0 | ḟ (s)|ds.
The existence of 2nπ-periodic solutions of systems ( 4) and ( 6) is obtained via Mawhin's continuation theorem [START_REF] Mawhin | Equivalence theorems for nonlinear operator equations and coincidence degree theory[END_REF][START_REF] Gaines | Coincidence Degree and Nonlinear Differential Equations[END_REF], which is presented below.

Mawhin's continuation theorem

Let X, Y be two infinite dimensional Banach spaces. A linear mapping

L : DomL ⊂ X → Y is

Existence of monotonic oscillatory solutions of (2)

4.1. Existence of an increasing oscillatory solution of (2) when γ > |ε| Theorem 6. If assumption (A) holds and γ > |ε|, then system (4) has at least one 2nπ-periodic solution.

Proof. Consider X = Y = {x = (x 1 , x 2 ) T ∈ C(R, R 2 ) : x(s) = x(s + 2nπ)} and the norm x = max s∈[0,2nπ] |x 1 (s)| + max s∈[0,2nπ]
X = Y = KerL ⊕ KerP = ImL ⊕ ImQ Considering K P Q = L -1 P (I -Q) we obtain that K P Q y = s 0 y(u)du - 1 2nπ s 0 2nπ 0 y(v)dvdu
Using the fact that N takes bounded sets into bounded sets, we deduce that N is L-compact on Ω for any open bounded set Ω ⊂ X.

In the followings, we will prove that assumption (a) from Theorem 5 holds. Indeed, assume that x(s) = (x 1 (s), x 2 (s)) T ∈ DomL is a solution of the operator equation Lx = λN x for a certain λ ∈ (0, 1). Therefore, x(s) = (x 1 (s), x 2 (s)) T is a smooth 2nπ-periodic function satisfying

⎧ ⎨ ⎩ x 1 (s) = λ[τe x1(s) -x 2 (s)e 2x1(s) ] x 2 (s) = λ[(γ -ε cos s)e x1(s) -δ] (7)
As x 1 and x 2 are continuous 2nπ-periodic functions, there exist

ζ i , η i ∈ [0, 2nπ] such that [x i ] L = x i (ζ i ) and [x i ] M = x i (η i ) for i ∈ {1, 2}
Integrating the second equation of system [START_REF] Shama | Guaranteed properties of gain schedule control for linear parameter varying plants[END_REF] on the interval [0, 2nπ] we obtain

A c c e p t e d m a n u s c r i p t

Using (8), we obtain

[x 1 ] L ≤ ln δ γ ≤ [x 1 ] M (9)
From ( 8) we also obtain 9) it results that

(γ -|ε|)
[x 2 ] L ≤ τγ δ ≤ [x 2 ] M (11)
Using ( 8) we get As τ > 0, δ > 0 and τ 2 > 4δ it results that the equation z 2τz + δ = 0 has two positive real roots 0 < z 1 < z 2 .

We have shown that y(ζ 1 ) = 0. Assume that there exists ξ > ζ 1 such that

y(ξ) = z 1 and y(s) < z 1 for any s ∈ [ζ 1 , ξ). As ẏ(ξ) = -λ(γ -ε cos ξ)e 2x1(ξ) < 0
it follows that y is strictly decreasing in a neighborhood of ξ, which is absurd, because y(s) < z 1 = y(ξ) for any s ∈ [ζ 1 , ξ).

Hence, we have shown that y(s) < z 1 for any s ∈ R. It results that x 2 (s)e x1(s) > τz 1 = z 2 > 0 and therefore, x 2 (s) > 0 for any s ∈ R.

A c c e p t e d m a n u s c r i p t

Integrating the first equation of system [START_REF] Shama | Guaranteed properties of gain schedule control for linear parameter varying plants[END_REF] on the interval [0, 2nπ] we obtain Using ( 9) and ( 10), it follows that

[x 1 ] M ≤ [x 1 ] L + 1 2 2nπ 0 |x 1 (s)|ds < [x 1 ] L + τ 2nπ 0 e x1(s) ds ≤ ln δ γ + 2nπδτ γ -|ε| (13)
and

[x 1 ] L ≥ [x 1 ] M - 1 2 2nπ 0 |x 1 (s)|ds > [x 1 ] M -τ 2nπ 0 e x1(s) ds ≥ ln δ γ - 2nπδτ γ -|ε| (14)
More, as x 2 (s) > z 2 e -x1(s) for any s ∈ [0, 2nπ] it follows that

[x 2 ] L > z 2 e -[x1]M = z 2 γ δ e -2nπδτ γ-|ε| (15)
We have obtained, based on inequalities ( 12)-( 15) that, if x(s) ∈ DomL is a solution of the operator equation Lx = λN x then, for any s ∈ R we have

ln δ γ - 2nπδτ γ -|ε| < x 1 (s) < ln δ γ + 2nπδτ γ -|ε| and z 2 γ δ e -2nπδτ γ-|ε| < x 2 (s) < τγ δ + 2nπδ
Therefore, considering the intervals

I 1 = ln δ γ - 2nπδτ γ -|ε| , ln δ γ + 2nπδτ γ -|ε| I 2 = z 2 γ δ e -2nπδτ γ-|ε| , τγ δ + 2nπδ
considering the open bounded set In order to check assumption (b) of Theorem 5, we observe that the equation

Ω = {x = (x 1 , x 2 ) T ∈ X : x 1 (s) ∈ I
QN x = 0 is equivalent to ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 2nπ 0 [τe x1(s) -x 2 (s)e 2x1(s) ]ds = 0 2nπ 0 (γ -ε cos s)e x1(s) ds = 2nπδ
The only solution of this system belonging to

KerL = R 2 is x = (x 1 , x 2 ) T = ln δ γ , τγ δ T ∈ Ω.
As for assumption (c) of Theorem 5, taking into account that ImQ = KerL = R 2 , we can chose the isomorphism J = Id. Considering the func-

tion ϕ : R 2 → R 2 defined by ϕ = QN | KerL , we obtain that ϕ(x 1 , x 2 ) = ⎡ ⎣ τe x1 -x 2 e 2x1 γe x1 -δ ⎤ ⎦ As ϕ(x) = 0 if and only if x = (x 1 , x 2 ) T = ln δ γ , τ γ δ T ∈ Ω, we evaluate det[J ϕ (x 1 , x 2 )] = γe 3x 1 = δ 3 γ 2 > 0
Hence the Brouwer degree is

deg B (JQN, Ω ∩ KerL, 0) = sign(det[J ϕ (x 1 , x 2 )]) = 1
In conclusion, there exists at least one 2nπ-periodic solution of system (4).

Therefore, the conclusion of this section can be summarized in the following Proposition:

Proposition 2.
If assumption (A) holds and γ > |ε|, then equation ( 2) has at least one solution θ(t) such that its derivative θ(t) is a positive periodic function (i.e. θ(t) is an increasing oscillatory solution).

A c c e p t e d m a n u s c r i p t

Proof. Consider the same Banach spaces endowed with the same norm as in the proof of Theorem 6.

Define the linear mapping L :

DomL = X ∩ C 1 (R, R 2 ) → Y by Lx = x = (x 1 , x 2 ) T .
Consider the nonlinear mapping N : X → Y defined by

N x = ⎡ ⎣ -τe x1 + x 2 e 2x1 (γ -ε cos s)e x1 + δ ⎤ ⎦
In the same way as in the proof of Theorem 6, it follows that L is a Fredholm operator of index zero and N is L-compact on Ω, for any open bounded set Ω ⊂ X.

In the followings, we will prove that assumption (a) from Theorem 5 holds.

Indeed, assume that x(s) = (x 1 (s), x 2 (s)) T ∈ DomL is a solution of the operator equation Lx = λN x for a certain λ ∈ (0, 1). Therefore,

x(s) = (x 1 (s), x 2 (s)) T is a smooth 2nπ-periodic function satisfying ⎧ ⎨ ⎩ x 1 (s) = λ[-τe x1(s) + x 2 (s)e 2x1(s) ] x 2 (s) = λ[(γ -ε cos s)e x1(s) + δ] (16)
As x 1 and x 2 are continuous 2nπ-periodic functions, there exist

ζ i , η i ∈ [0, 2nπ] such that [x i ] L = x i (ζ i ) and [x i ] M = x i (η i ) for i ∈ {1, 2}
Integrating the second equation of system [START_REF] Goto | Problems of analysis methods for studying nonlinear flight dynamics[END_REF] on the interval [0, 2nπ] we obtain

2nπ 0 (γ -ε cos s)e x1(s) ds = -2nπδ (17)
Taking into account that γ < -|ε|, we can see that γε cos s < γ + |ε| < 0 and

e [x1]M 2nπ 0 (γ-ε cos s)ds ≤ 2nπ 0 (γ-ε cos s)e x1(s) ds ≤ e [x1]L 2nπ 0 (γ-ε cos s)ds
Using [START_REF] Etkin | Dynamics of Flight: Stability and Control[END_REF], we obtain

[x 1 ] L ≤ ln δ |γ| ≤ [x 1 ] M (18)
From [START_REF] Etkin | Dynamics of Flight: Stability and Control[END_REF] we also obtain

(γ -|ε|)
As x 1 (ζ 1 ) = [x 1 ] L it follows that x 1 (ζ 1 ) = 0 and from the first equation of system [START_REF] Goto | Problems of analysis methods for studying nonlinear flight dynamics[END_REF] we obtain that

x 2 (ζ 1 ) = τe -[x1]L and therefore τe -[x1]L ≤ [x 2 ] M .
In a similar way we can show that τe

-[x1]M ≥ [x 2 ] L . By means of (18) it results that [x 2 ] L ≤ τ |γ| δ ≤ [x 2 ] M (20)
Using [START_REF] Etkin | Dynamics of Flight: Stability and Control[END_REF] we get Using ( 18) and ( 19), it follows that

[x 1 ] M ≤ [x 1 ] L + 1 2 2nπ 0 |x 1 (s)|ds < [x 1 ] L + τ 2nπ 0 e x1(s) ds ≤ ln δ |γ| + 2nπδτ |γ| -|ε| (22)
and

[x 1 ] L ≥ [x 1 ] M - 1 2 2nπ 0 |x 1 (s)|ds > [x 1 ] M -τ 2nπ 0 e x1(s) ds ≥ ln δ |γ| - 2nπδτ |γ| -|ε| (23)
More, as x 2 (s) > -z 1 e -x1(s) for any s ∈ [0, 2nπ] it follows that

[x 2 ] L > -z 1 e -[x1]M = |z 1 γ| δ e -2nπδτ |γ|-|ε| (24)
We have obtained, based on inequalities ( 21

I 2 = |z 1 γ| δ e -2nπδτ |γ|-|ε| , τ |γ| δ + 2nπδ
considering the open bounded set Ω = {x = (x 1 , x 2 ) T ∈ X : x 1 (s) ∈ I 1 and x 2 (s) ∈ I 2 for any s ∈ R} it follows that assumption (a) of Theorem 5 is fulfilled.

In order to check assumption (b) of Theorem 5, we observe that the equation

QN x = 0 is equivalent to ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩
x = (x 1 , x 2 ) T = ln δ |γ| , τ |γ| δ T ∈ Ω, we evaluate det[J ϕ (x 1 , x 2 )] = γe 3x 1 = δ 3 γ 2 > 0
Hence the Brouwer degree is

deg B (JQN, Ω ∩ KerL, 0) = sign(det[J ϕ (x 1 , x 2 )]) = 1
In conclusion, there exists at least one 2nπ-periodic solution of system [START_REF] Ionita | Delay induced oscillations[END_REF].

The conclusion of this section can be summarized in the following Proposition:

Proposition 3. If assumption (A) holds and γ < -|ε|, then equation ( 2) has at least one solution θ(t) such that its derivative θ(t) is a negative periodic function (i.e. θ(t) is a decreasing oscillatory solution).

Oscillations in longitudinal flight regime

Taking into account the relationship between system (1) and equation ( 2), the following result is obtained: Proposition 4. If Assumption (A) is fulfilled and |γ| > |ε|, then there exists at least one oscillatory solution for system (1). More precisely:

a. If γ > |ε|, then there exists at least one solution (α(t), q(t), θ(t)) T of system

(1) such that α(t) and q(t) = θ(t) are periodic functions and q(t) is strictly positive.

A c c e p t e d m a n u s c r i p t

b. If γ < -|ε|, then there exists at least one solution (α(t), q(t), θ(t)) T of system (1) such that α(t) and q(t) = θ(t) are periodic functions and q(t) is strictly negative.

Numerical example

Numerical simulations have been carried out for the ALFLEX reentry vehicle (see numerical data in the appendix). In this case, the interval J given by Proposition 1 is J = [2.9944 • , 3.8727 • ]. When the elevator deflection δ e is set to values outside this interval, oscillatory behavior is expected, according to the previous theoretical results.

The case δ e < 2.9944 • corresponds to γ > |ε|. A large number of simulations have been completed, using a wide set of initial conditions. Fig. 1 shows that setting δ e = 2.5 • the solutions of system (1) with the initial condition (8.18 • , 0, -9.16 • ) T , converge towards a presumably asymptotically stable oscillatory orbit (α(t), q(t), θ(t)) T , where α(t) and q(t) are periodic functions and θ(t) is strictly increasing. The same phenomenon can be noticed for δ e = 2.8 • , observing that the period of α(t) and q(t) becomes larger. Numerical simulations show that the period increases infinitely as δ e approaches the critical value δ e = 2.9944 • . This suggests that at the critical value δ e = 2.9944 • , the existing saddle-node bifurcation points (given by Proposition 1) are linked by a chain of heteroclinic orbits. More, this allows us to believe that when δ e = δ e = 2.9944 • , a global bifurcation takes place which involves this whole chain of heteroclinic orbits, resulting in the birth of an oscillatory orbit when δ e < δ e .

The case δ e > 3.8727 • corresponds to γ < -|ε|. Fig. 2 shows that setting δ e = 4 • the solutions of system (1) with the initial condition (8.18 • , 0, -9.16 • ) T , converge towards a presumably asymptotically stable oscillatory motion (α(t), q(t), θ(t)) T , where α(t) and q(t) are periodic functions and θ(t) is strictly decreasing. The same phenomenon can be noticed for δ e = 5 • , observing that the period of α(t) and q(t) becomes smaller. Numerical simulations show that the period increases infinitely as δ e approaches the critical value δ e = 3.8727 • . The same

A c c e p t e d m a n u s c r i p t

phenomenon can be noticed as in the previous paragraph.

Conclusions

A theoretical proof has been given for the existence of oscillatory solutions in the case of a simplified mathematical model of longitudinal flight dynamics, whenever the value of the elevator deflection δ e belongs to (-∞, δ e ) ∪ (δ e , ∞)

(where the interval J = [δ e , δ e ] corresponds to the existence of equilibria in system (1)). The mathematical tools involved in this proof are the coincidence degree theory and Mawhin's continuation theorem.

It has been emphasized that at the critical values δ e and δ e of the elevator deflection, an infinity of saddle-node bifurcations take place. In the case of system (1), Hopf bifurcations are not present, and hence, a Hopf bifurcation analysis such as in [START_REF] Liaw | Analysis of longitudinal flight dynamics: A bifurcationtheoretic approach[END_REF][START_REF] Liaw | Two-parameter bifurcation analysis of longitudinal flight dynamics[END_REF][START_REF] Pereira | On nonlinear dynamics and an optimal control design to a longitudinal flight[END_REF] is not an appropriate tool for depicting oscillatory behavior.

Numerical simulations suggest that the saddle-node bifurcation points are linked by a chain of heteroclinic orbits, allowing us to believe that the onset of oscillatory behavior, when δ e exceeds the critical values δ e and δ e , is caused by a global bifurcation which involves this whole chain of heteroclinic orbits. A more thorough investigation of this phenomenon, as well as extending these results

to more general systems describing longitudinal flight dynamics, may constitute directions for future research.

u

  = first, second, third order derivative of u in Newton's notation u , u = first, second order derivative of u in Lagrange's notation S = closure of the set S ∂S = boundary of the set S N = set of positive non-zero integers Z = set of integers R = set of real numbers DomL = domain of the operator L ImL = codomain (range) of the operator L KerL = kernel of the operator L

A c c e p t e d m a n u s c r i p t

  aircraft, when the automatic flight control system is decoupled and the elevator deflection exceeds the bifurcation values. Sufficient conditions are obtained for the existence of oscillatory solutions for any value of the elevator deflection δ e

satisfying 2nπ 0 e

 0 x(s) ds = T . (b) If for n ∈ N , there exists a 2nπ-periodic solution x(s) of (3), then there exists an increasing oscillatory solution θ(t) of (2) satisfying θ(t + T ) = θ(t) + 2nπ, where T = 2nπ 0 e x(s) ds.

) satisfying 2nπ 0 e

 0 x(s) ds = T .(b) If for n ∈ N , there exists a 2nπ-periodic solution x(s) of (5), then there exists a decreasing oscillatory solution θ(t) of (2) satisfying θ(t+T ) = θ(t)-2nπ, where T = 2nπ 0 e x(s) ds.

  s) ds = -ψ(2nπ) + ψ(0) = T . By successive differentiation, taking into account that θ(ψ(s)) = s, we obtain:

A c c e p t e d m a n u s c r i p t θ = ψ - 1 :

 1 (b) Suppose that equation (5) has a 2nπ-periodic solution x(s). Consider θ 0 ∈ R. It can be easily seen that the function ψ : R → R defined by ψ(s) = -s θ 0 e x(u) du is strictly decreasing and bijective. Consider the inverse function R → R which is also strictly decreasing. One has θ(t) = 1 ψ(θ(t))

  u) du = -ψ(s) + T where T = 2nπ 0 e x(u) du. Making s = θ(t) in the equality ψ(s + 2nπ) = ψ(s) -T , we get ψ(θ(t) + 2nπ) = t -T and applying θ to this equality, we obtain: θ(t -T ) = θ(t) + 2nπ and hence:

A c c e p t e d m a n u s c r i p t by K P Q = L - 1 PTheorem 5 (

 15 called a Fredholm mapping if KerL has finite dimension and ImL is closed and has finite codimension. The index of a Fredholm mapping L is the integer i(L) = dim KerLcodim ImL. In the followings, consider L : DomL ⊂ X → Y a Fredholm mapping of index zero which is not injective. Let P : X → X and Q : Y → Y be continuous projectors such that KerQ = ImL, ImP = KerL, X = KerL ⊕ KerP and Y = ImL ⊕ ImQ. It is clear that the operator L P = L| DomL∩KerP : DomL ∩ KerP → ImL is an isomorphism. Consider the operator K P Q : Y → X defined (I -Q). Let Ω ⊂ X be an open bounded set and N : Ω → Y be a continuous mapping. We say that N is L-compact if K P Q N is compact, QN is continuous and QN ( Ω) is a bounded set in Y . Since ImQ is isomorphic to KerL, there exists an isomorphism J : ImQ → KerL. Continuation theorem [19, 20]). Let Ω ⊂ X be an open bounded set, let L be a Fredholm mapping of index zero and let N be L-compact on Ω. Assume (a) Lx = λN x for any λ ∈ (0, 1) and x ∈ ∂Ω ∩ DomL; (b) QN x = 0 for any x ∈ KerL ∩ ∂Ω; (c) Brouwer degree deg B (JQN, Ω ∩ KerL, 0) = 0. Then Lx = N x has at least one solution in DomL ∩ Ω.

|x 2

 2 (s)|, x∈ X = Y X and Y endowed with the norm • are Banach spaces. Define the linear mappingL : DomL = X ∩ C 1 (R, R 2 ) → Y by Lx = x = (x 1 , x 2 ) T . Consider the nonlinear mapping N : X → Y defined by N x = ⎡ ⎣ τe x1x 2 e 2x1(γε cos s)e x1δ It can be easily seen that KerL = R 2 , ImL = {y ∈ Y : 2nπ 0 y(s)ds = 0} and L is a Fredholm mapping of index zero. Define the continuous projectors P : X → X and Q : Y → Y by P x = x(0) ds = ȳ. Then ImP = KerL, KerQ = ImL, ImQ = KerL = R 2 and

8 )

 8 ε cos s)e x1(s) ds = 2nπδ (Taking into account that γ > |ε|, we can see that γε cos s > γ -|ε| > 0 and e ε cos s)e x1(s) ds ≤ e [x1]M 2nπ 0 (γ-ε cos s)ds

2nπ 0 e

 0 x1(s) ds ≤ 2nπδ ≤ (γ + |ε|) x 1 (ζ 1 ) = [x 1 ] L it follows that x 1 (ζ 1 ) = 0 andfrom the first equation of system (7) we obtain that x 2 (ζ 1 ) = τe -[x1]L and therefore τe -[x1]L ≤ [x 2 ] M . In a similar way we can show that τe -[x1]M ≥ [x 2 ] L . By means of (

  -ε cos s)e x1(s) -δ| < 2nπ 0 (γ-ε cos s)e x1(s) +2nπδ = 4nπδ and therefore [x 2 ] M ≤ [x 2 ] Denoting y(s) = τx 2 (s)e x1(s) , it results that ẏ(s) = λe x1(s) [y(s) 2τy(s) + δ -(γε cos s)e x1(s) ]

  s)x 2 (s)e 2x1(s)

1 A c c e p t e d m a n u s c r i p t it follows that

 1 and x 2 (s) ∈ I 2 for any s ∈ R} assumption (a) of Theorem 5 is fulfilled.

4. 2 .Theorem 7 .

 27 Existence of a decreasing oscillatory solution of (2) when γ < -|ε| If assumption (A) holds and γ < -|ε|, then system (6) has at least one 2nπ-periodic solution.

2nπ 0 e

 0 x1(s) ds ≤ -2nπδ ≤ (γ + |ε|)

2nπ 0 |x 2 0 x 2 x 2

 02022 (s)|ds = λ 2nπ 0 |(γ-ε cos s)e x1(s) +δ| < 2nπ 0 (γε cos s)e x1(s) +2nπδ = 4nπδ and therefore[x 2 ] M ≤ [x 2 ] s) = -τ + x 2 (s)e x1(s) , it results that ẏ(s) = λe x1(s) [y(s) 2 + τy(s) + δ + (γε cos s)e x1(s) ]As τ > 0, δ > 0 and τ 2 > 4δ it results that the equation z 2 + τz + δ = 0 has two negative real roots z 1 < z 2 < 0.We have shown that y(ζ 1 ) = 0. Assume that there exists ξ < ζ 1 such that y(ξ) = z 2 and y(s) > z 2 for any s ∈ (ξ, ζ 1 ]. As ẏ(ξ) = λ(γε cos ξ)e 2x1(ξ) < 0 it follows that y is strictly decreasing in a neighborhood of ξ, which is absurd, because y(s) > z 2 = y(ξ) for any s ∈ (ξ, ζ 1 ].Hence, we have shown that y(s) > z 2 for any s ∈ R. It results that x 2 (s)e x1(s) > τ + z 2 = -z 1 > 0 and therefore, x 2 (s) > 0 for any s ∈ R.Integrating the first equation of system[START_REF] Goto | Problems of analysis methods for studying nonlinear flight dynamics[END_REF] on the interval [0, 2nπ] we obtain 2nπ (s)e 2x1(s) ds = τ (s)e 2x1(s) ds = 2τ

A c c e p t e d m a n u s c r i p t⎦

  x1(s) + x 2 (s)e 2x1(s) ]ds = 0 2nπ 0 (γε cos s)e x1(s) ds = -2nπδ The only solution of this system belonging to KerL = R 2 is x = (x 1 , x 2 ) (c) of Theorem 5, taking into account that ImQ = KerL = R 2 , we can chose the isomorphism J = Id. Considering the function ϕ : R 2 → R 2 defined by ϕ = QN | KerL , we obtain that ϕ(x 1 , x 2 ) = ⎡ ⎣ -τe x1 + x 2 e 2x1 γe x1 + δ ⎤ As ϕ(x) = 0 if and only if

Figure 2 :

 2 Figure 2: Evolution of the state parameters α, q ad θ of system (1) in time, considering the initial condition (8.18 • , 0, -9.16 • ) T , and choosing δe = 4 • (left) and δe = 5 • (right).

A c c e p t e d m a n u s c r i p t A.

 t Numerical data for the ALFLEX reentry vehicle m = 760 g = 9.81 δ e0 = 3 • α 0 = 8.18 • θ 0 = -9.16 • S = 9.45 b = 3.295 c = 3.154 ρ = 1.156 k = ρS/2 I x = 407 I y = 1366 I z = 1634 I xz = 10.4 V = 73.84 C yβ = -0.6849 C L = 0.2387 C lβ = -0.1774 C nβ = -0.0657 C yδr = 0.1907 C Lα = 2.016 C lp = -0.007 C np = 0.0032 C yr = 0 C Lδe = 0.6355 C lr = 0.004 C nr = -0.006 C mα = -0.0134 C D = 0.0745 C lδa = 0.1488 C nδa = -0.0266 C mq = -0.0474 C Dα = 0.2714 C lδr = 0.0788 C nδr = -0.099 C mδe = -0.2152 C Dδe = 0.1019 C zα = -C Lα cos α 0 + C L sin α 0 -C D cos α 0 -C Dα sin α 0 C zδe = -C Lδe cos α 0 -C Dδe sin α 0 a 11 = kV m C zα = -1.10058 a 21 = ckV 2 C mα = -0.92142 a 22 = ckV 2 C mq = -3.25936 b 1 = kV m C zδe = -0.34151 b 2 = ckV 2 C mδe = -14.79777

  )-(24) that, if x(s) ∈ DomL is a solution of the operator equation Lx = λN x then, for any s ∈ R we have

	ln	δ |γ|	-	2nπδτ |γ| -|ε|	< x 1 (s) < ln	δ |γ|	+	2nπδτ |γ| -|ε|
	and
		|z 1 γ| δ	e -2nπδτ |γ|-|ε| < x 2 (s) <	τ |γ| δ	+ 2nπδ
	Therefore, considering the intervals
	I 1 = ln	δ |γ|	-	2nπδτ |γ| -|ε|	, ln	δ |γ|	+	2nπδτ |γ| -|ε|

Acknowledgements

The authors would like to thank the anonymous referees for their valuable comments and suggestions, which led to an improvement of the paper's presentation.

34C15✩ This work has been supported by the Romanian National Authority for Scientific Research under the grant ID-354-PN-II-7/2007.


Attila Pálfalvi

Efficient solution of a vibration equation involving fractional derivatives

Keywords: vibration, fractional derivative, fractional differential equation, Adomian decomposition

Fractional order (or, shortly, fractional) derivatives are used in viscoelasticity since the late 1980's, 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.

Introduction

It is common to date fractional calculus back to Leibniz, and several textbooks are available on the subject (e.g. [START_REF] Podlubny | Fractional differential equations[END_REF]). Several authors have been examining the possibility of using fractional derivatives in material modelling in the last two-three decades [START_REF] Bagley | A theoretical basis for the application of fractional calculus to viscoelasticity[END_REF][START_REF] Bagley | On the fractional calculus model of viscoelastic behavior[END_REF], and the field is of a growing interest

A c c e p t e d m a n u s c r i p t

nowadays [START_REF] Alvelid | Modelling of constrained thin rubber layer with emphasis on damping[END_REF][START_REF] Ramrakhyani | Modeling of elastomeric materials using nonlinear fractional derivative and continuously yielding friction elements[END_REF][START_REF] Freed | Fractional calculus in biomechanics: A 3D viscoelastic model using regularized fractional derivative kernels with application to the human calcaneal fat pad[END_REF][START_REF] Heymans | Fractional calculus description of non-linear viscoelastic behaviour of polymers[END_REF][START_REF] Hernandez-Santiago | Relaxation modulus in PMMA and PTFE fitting by fractional Maxwell model[END_REF][START_REF] Schmidt | Finite element formulation of viscoelastic constitutive equations using fractional time derivatives[END_REF]. Finite element formulations are also in development [START_REF] Schmidt | Finite element formulation of viscoelastic constitutive equations using fractional time derivatives[END_REF][START_REF] Cortés | Homogenised finite element for transient dynamic analysis of unconstrained layer damping beams involving fractional derivative models[END_REF][START_REF] Enelund | Time domain modeling of damping using anelastic displacement fields and fractional calculus[END_REF].

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 time-history integral. There are several solution methods available in the literature.

The analytical solution is theoretically known [START_REF] Podlubny | Fractional differential equations[END_REF]. 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 Grünwald-Letnikov definition of fractional derivatives, used e.g. by Schmidt and Gaul [START_REF] Schmidt | Finite element formulation of viscoelastic constitutive equations using fractional time derivatives[END_REF]. These authors proceeded to reduce the volume of the time history integral [START_REF] Schmidt | On the numerical evaluation of fractional derivatives in multi-degree-of-freedom systems[END_REF], achieving an almost two-magnitude gain in calculation time on an example.

Suarez and Shokooh [START_REF] Suarez | An eigenvector expansion method for the solution of motion containing fractional derivatives[END_REF] solved a vibration equation with fractional damping of order 1/2 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 problem of the vibration system, which is uncommon in the literature.

Yuan and Agrawal [START_REF] Yuan | A numerical scheme for dynamic systems containing fractional derivatives[END_REF] 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 [START_REF] Schmidt | On a critique of a numerical scheme for the calculation of fractionally damped dynamical systems[END_REF] 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 [START_REF] Pálfalvi | Methods for solving a semi-differential vibration equation[END_REF] for the field of interest of this paper.

The Adomian decomposition method (ADM) has been introduced by George Adomian in the late 1980's. Essentially, it approximates the solution of a non-linear differential equation with a series of functions. The method

A c c e p t e d m a n u s c r i p t

is getting into use for the solution of fractional differential equations [START_REF] Wang | Numerical solutions for fractional KdV-Burgers equation by Adomian decomposition method[END_REF][START_REF] Daftardar-Gejji | Adomian decomposition: a tool for solving a system of fractional differential equations[END_REF][START_REF] Ray | Analytical solution of a dynamic system containing fractional derivative of order one-half by Adomian decomposition method[END_REF][START_REF] Daftardar-Gejji | Solving a multi-order fractional differential equation using Adomian decomposition[END_REF][START_REF] Daftardar-Gejji | Solving multi-term linear and nonlinear diffusion-wave equations of fractional order by Adomian decomposition method[END_REF][START_REF] Shawagfeh | Analytical approximate solutions for nonlinear fractional differential equations[END_REF][START_REF] Momani | A numerical scheme for the solution of multi-order fractional differential equations[END_REF].

For a special case, Saha Ray, Poddar and Bera [START_REF] Ray | Analytical solution of a dynamic system containing fractional derivative of order one-half by Adomian decomposition method[END_REF] have proven that the Adomian solution converges to the analytical solution. Recently, Hu, Luo and Lu [START_REF] Hu | Analytical solution of the linear fractional differential equation by Adomian decomposition method[END_REF] 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 [START_REF] Momani | Numerical approach to differential equations of fractional order[END_REF][START_REF] Momani | Numerical comparison of methods for solving linear differential equations of fractional order[END_REF], the method of Atanackovic and Stankovic [START_REF] Atanackovic | An expansion formula for fractional derivatives and its application[END_REF][START_REF] Atanackovic | On a numerical scheme for solving differential equations of fractional order[END_REF] and the method of Singh and Chatterjee [START_REF] Singh | Galerkin projections and finite elements for fractional order derivatives[END_REF]) that will not be treated in this paper. Also, the Adomian decomposition method has been improved by some authors [START_REF] Jafari | Revised Adomian decomposition method for solving a system of nonlinear equations[END_REF][START_REF] Jafari | Revised Adomian decomposition method for solving systems of ordinary and fractional differential equations[END_REF][START_REF] Ray | Application of modified decomposition method for the analytical solution of space fractional diffusion equation[END_REF]. 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 Grünwald-Letnikov definition will also be calculated.

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.

A c c e p t e d m a n u s c r i p t

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.

The equation to be solved

The differential equation to be solved is the vibration equation with fractional damping, with one degree of freedom:

D 2 x(t) + c m D α x(t) + k m x(t) = f (t) , (1)
where D = d dt is the differential operator. Another common form of Equation (1) is

D 2 x(t) + 2ηω 2-α n D α x(t) + ω 2 n x(t) = f (t) , (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

D -q as a D -q t x(t) = 1 Γ(q) t a (t -τ ) q-1 x(τ ) dτ , ( 3 )
where q > 0 and Γ(x) is the Gamma function

Γ(x) = ∞ 0 e -z z x-1 dz. ( 4 )
For a continuous x(t),

D -p D -q x(t) = D -(p+q) x(t) , (5)
as given in [START_REF] Podlubny | Fractional differential equations[END_REF] (if both p and q are non-negative).

A c c e p t e d m a n u s c r i p t

The fractional integral of a polynomial will be needed later so it is given below:

D -q t ν = Γ(ν + 1) Γ(ν + 1 + q) t q+ν , ( 6 )
With the fractional integral operator, fractional derivatives are easily introduced.

Fractional derivatives

For a real α > 0, D α is defined by the Riemann-Liouville definition [START_REF] Podlubny | Fractional differential equations[END_REF],

using the above fractional integral operator:

a D α t x(t) = d dt n a D -(n-α) t = 1 Γ(n -α) d dt n t a (t -τ ) n-1-α x(τ ) dτ . ( 7
)
Another choice is the Caputo definition

C a D α t x(t) = 1 Γ(n -α) t a (t -τ ) n-1-α d dτ n x(τ ) dτ . (8)
In both cases, (n -1) < α < n.

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

a D α t x(t) = C a D α t x(t) + 1 Γ (n -α) n-1 k=0 Γ(n -α) Γ(k -α + 1) (t -a) k-α x (k) (a) . (9)
In the applications, D practically always means 0 D t , and most authors use the Riemann-Liouville, or the mathematically equivalent Grünwald-Letnikov definition (see [START_REF] Podlubny | Fractional differential equations[END_REF] 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 [START_REF] Heymans | Physical interpretation of initial conditions for fractional differential equations with Riemann-Liouville fractional derivatives[END_REF].

The composition (5) can be extended to an integral and a derivative:

D -p D q x(t) = D q D -p x(t) = D q-p x(t) , (10)
with p, q ≥ 0 and a continuous x(t). (However, the composition of two fractional order derivatives is not straightforward.)

A c c e p t e d m a n u s c r i p t 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 [START_REF] Wang | Numerical solutions for fractional KdV-Burgers equation by Adomian decomposition method[END_REF][START_REF] Daftardar-Gejji | Adomian decomposition: a tool for solving a system of fractional differential equations[END_REF][START_REF] Ray | Analytical solution of a dynamic system containing fractional derivative of order one-half by Adomian decomposition method[END_REF][START_REF] Daftardar-Gejji | Solving a multi-order fractional differential equation using Adomian decomposition[END_REF][START_REF] Daftardar-Gejji | Solving multi-term linear and nonlinear diffusion-wave equations of fractional order by Adomian decomposition method[END_REF][START_REF] Shawagfeh | Analytical approximate solutions for nonlinear fractional differential equations[END_REF][START_REF] Momani | A numerical scheme for the solution of multi-order fractional differential equations[END_REF]. As it will be also used in this paper, it is presented in the followings.

The method has been described very clearly in [START_REF] Chen | An algorithm for Adomian decomposition method[END_REF]. The non-linear 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) = x 0 (t) + x 1 (t) + x 2 (t) + . . . , (12)
using

x 0 (t) = L -1 Lx(t) -x(t) + L -1 f (t) (13)
and

x n+1 (t) = -L -1 Rx n (t) -L -1 A n (t) (14)
with

A n (t) = 1 n! d n dλ n N ∞ i=0 λ i x i (t) λ=0 . ( 15
)
It has been proven recently by Hu, Luo and Lu [START_REF] Hu | Analytical solution of the linear fractional differential equation by Adomian decomposition method[END_REF] 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 followings, it will be shown that it can be very efficient, too.

A c c e p t e d m a n u s c r i p t

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.

3.1. Analytical solution for α = 1/2

Suarez and Shokooh [START_REF] Suarez | An eigenvector expansion method for the solution of motion containing fractional derivatives[END_REF] calculate the analytical solution of Equation ( 2)

for the special case α = 1/2, both for a free vibration and a step function excitation. For the former, they obtain

x IC (t) = 4 j=1 Ψ 4j R j 1 √ πt + 4 j=1 Ψ 4j R j λ j g j (t) , (16)
where

g j (t) = e λ 2 j t 1 + erf λ j √ t (17)
and λ j and Ψ j are the solution of the eigenproblem

AΨ j = λ j BΨ j (18)
with

A = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 -ω 2 n ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ and B = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 2ηω 3/2 n ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
, and Ψ 4j is the 4th coordinate of the jth eigenvector. In the above,

erf(x) = 2 √ π x 0 e -t 2 dt
is the error function, and R j comes from the initial conditions, by

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ R 1 R 2 R 3 R 4 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = Ψ 1 Ψ 2 Ψ 3 Ψ 4 -1 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ v 0 0 x 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 19
)

A c c e p t e d m a n u s c r i p t

The response to a step function force (using the Heaviside function H(t))

f (t) = f 0 H(t) ,
where X 0 and V 0 are considered zero, is

x H (t) = f 0 4 j=1 Ψ 2 4j λ j (g j (t) -1) . ( 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.

Grünwald-Letnikov formulation

It is common in the literature, especially in applications, to use the Grünwald-Letnikov definition of a fractional derivative [START_REF] Podlubny | Fractional differential equations[END_REF]:

D α f (t) = lim N →∞ ⎡ ⎣ t N -α N -1 j=0 Γ(j -α) Γ(-α) Γ(j + 1) f t -j t N ⎤ ⎦ . (21)
This, along with any time stepping scheme, allows to find a numerical solution for Equation [START_REF] Podlubny | Fractional differential equations[END_REF]. For simplicity, an explicit scheme has been used for time stepping in this paper:

v i+1 = v i + a i h (22)
and

x i+1 = x i + v i h + a i h 2 2 , ( 23
)
where

a i = f i - k m x i - c m h -α Γ(-α) i-1 j=0 Γ(j -α) Γ(j + 1) x i-j-1 , (24)
h being the time step size and f i being the force per unit mass at t = ih.

A c c e p t e d m a n u s c r i p t 4. New solutions

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.

Proposed method: Taylor-Adomian series

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

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 = D 2 , R = c m D α + k m and N = 0.
In this case, conforming to Equations ( 13) and ( 14),

x 0 (t) = D -2 D 2 x(t) -x(t) + D -2 f (t) (25)
and

x n+1 (t) = - k m D -2 x n (t) - c m D α-2 x n (t) , (26)
the latter leading to

x n (t) = (-1) n m n n j=0 n j c n-j k j D -((2-α)n+αj) x 0 (t) (27)
using the composition property [START_REF] Cortés | Homogenised finite element for transient dynamic analysis of unconstrained layer damping beams involving fractional derivative models[END_REF] of an integral and a derivative operator.

A c c e p t e d m a n u s c r i p t

In Equation [START_REF] Momani | Numerical approach to differential equations of fractional order[END_REF], the first term describes the contribution of initial conditions:

x IC 0 (t) = D -2 D 2 x(t) -x(t) = X 0 + V 0 t, (28)
where X 0 and V 0 are the initial position and velocity, respectively. With Equations ( 27) and ( 6), one can write

x IC n (t) = (-1) n m n t (2-α)n n j=0 n j c n-j k j t jα X 0 Γ(2n + 1 -(n -j) α) + V 0 t Γ(2n + 2 -(n -j) α) . ( 29
)
This result has also been obtained by Baclic and Atanackovic [START_REF] Baclic | Stability and creep of a fractional derivative order viscoelastic rod[END_REF] using Laplace transformation.

The second term of Equation [START_REF] Momani | Numerical approach to differential equations of fractional order[END_REF] gives the response of the system to the excitation:

x f 0 (t) = D -2 f (t) . ( 30
)
Suppose an excitation of the form of the Taylor series

f (t) = ∞ i=0 T i t i , (31)
T i s being the coefficients of the polynomial. This, together with Equations ( 30) and ( 6), leads to

x f 0,i (t) = T i Γ(i + 1) Γ(i + 3) t i+2 (32)
and

x f n,i (t) = T i (-1) n m n Γ(i + 1) t i+(2-α)n+2 n j=0 n j c n-j k j Γ(i + 3 + 2n -(n -j) α) t jα , (33)
x f n,i (t) being the nth term of the Adomian series of the response to the ith term of the Taylor series of the excitation. The total response is

x f (t) = ∞ n=0 ∞ i=0 x f n,i (t) . (34)
Thus, the solution of Equation ( 1) has been given.

A c c e p t e d m a n u s c r i p t

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, Equations ( 29) and [START_REF] Heymans | Physical interpretation of initial conditions for fractional differential equations with Riemann-Liouville fractional derivatives[END_REF] show that terms t jα 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 n inc terms of the series: if the sum of the next n inc terms is larger than a small limit ('error goal'), then the number of terms is increased by n inc , and the error estimation will be based on the following n inc terms. This is applied separately for both series.

Another point is the required computational effort. The first two terms of the sum in Equation ( 29) give

x IC (t) = X 0 +V 0 t- c m X 0 Γ(3 -α) t 2-α - k m X 0 Γ(3) t 2 - c m V 0 Γ(4 -α) t 3-α - k m V 0 Γ(4) t 3 ±. . . (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 α is restricted to be rational, as α = p q with p and q both integers, the right-hand side of Equation ( 35) becomes a polynomial for t 1/q :

x IC (t) = X 0 + V 0 t 1/q q - c m X 0 Γ(3 -α) t 1/q 2q-p - k m X 0 Γ(3) t 1/q 2q - c m V 0 Γ(4 -α) t 1/q 3q-p - k m V 0 Γ(4) t 1/q 3q ± . . . (36)
This can be computed efficiently.

A c c e p t e d m a n u s c r i p t

Examples

Step function excitation. For an excitation with the Heaviside function H(t)

of the form

f (t) = f 0 H(t) ,
solved by Suarez and Shokooh [START_REF] Suarez | An eigenvector expansion method for the solution of motion containing fractional derivatives[END_REF] for α = 1/2, only T 0 is non-zero of the

T i s, resulting in x f n (t) = f 0 (-1) n m n t (2-α)n+2 n j=0 n j c n-j k j Γ(3 + 2n -(n -j) α) t jα . ( 37
)
This has also been obtained by Saha Ray, Poddar and Bera [START_REF] Ray | Analytical solution of a dynamic system containing fractional derivative of order one-half by Adomian decomposition method[END_REF] for α = 1/2 using Adomian decomposition (with a slightly different selection for L, R, and N).

Sine excitation. For a sine excitation of the form

f (t) = f 0 sin(ω e t) ,
the Taylor coefficients are

T i = ⎧ ⎪ ⎨ ⎪ ⎩ 0 i f i is even (-1) (i-1)/2 f 0 ω i e i! if i is odd . Using i = 2h + 1 gives T h = (-1) h ω 2h+1 e (2h + 1)! ,
which, inserted into Equation ( 33), leads to

x f n,h (t) = f 0 (-1) n+h m n ω 2h+1 e t (2-α)n+3+2h n j=0 n j c n-j k j Γ(4 + 2 (n + h) -(n -j) α) t jα , ( 38
)
also using the property

Γ(n + 1) = n!
of the Gamma function.

A c c e p t e d m a n u s c r i p t 4.2. Analytical solution by Adomian decomposition

Equation [START_REF] Atanackovic | An expansion formula for fractional derivatives and its application[END_REF] gives the nth term of the Adomian series of x(t) if the original differential equation is Equation [START_REF] Podlubny | Fractional differential equations[END_REF]. For a sine excitation

f (t) = f 0 sin(ω e t) ,
this yields to

x n (t) = (-1) n m n n j=0 n j c n-j k j D -((2-α)n+αj) D -2 f 0 sin(ω e t) .
Using the definition (3) of the integral operator, one obtains

x 0 (t) = t 0 (t -τ ) sin(ω e τ ) dτ = 1 ω 2 e
(ω e tsin(ω e t)) .

(39)

This leads to

x n (t) = (-1) n m n f 0 ω 2 e n j=0 n j c n-j k j 1 Γ(q) t 0 (t -τ ) q-1 (ω e τ -sin(ω e τ )) dτ , (40)
where

q = (2 -α) n + αj.
This results in

x n (t) = (-1) n m n f 0 ω 2 e n j=0 n j c n-j k j 1 Γ(q + 2) √ t ω 1/2-q e s 3/2+q,1/2 (ω e t) , (41)
for n > 0, where s μ,ν (t) is the Lommel-function [START_REF] Weisstein | From MathWorld -A Wolfram Web Resource[END_REF]

s μ,ν (t) = 1 F 2 1; 1 2 (μ -ν + 3) , 1 2 (μ + ν + 3) ; - 1 4 t 2 t μ+1 (μ + 1) 2 -ν 2 with 1 F 2 (a 1 ; b 1 , b 2 ; t) being the hypergeometric function 1 F 2 (a 1 ; b 1 , b 2 ; t) = ∞ k=0 (a 1 ) k (b 1 ) k (b 2 ) k t k k! , where (a) k is the Pochhammer symbol (a) k = Γ(a + k) Γ(a) .
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.

A c c e p t e d m a n u s c r i p t

Numerical examples

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.

Software

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 [START_REF]GNU multiple precision arithmetic library[END_REF] and one of its derivatives, MPFR [START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with correct rounding[END_REF].

The Grünwald-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 Equations ( 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 10 -25 for ten terms.

Problems and results

In the followings, numerical examples and their results will be presented and compared. For the latter purpose, let us introduce an absolute error indicator ε as

ε = 1 N N i=1 x num i -x ref 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.

A c c e p t e d m a n u s c r i p t

The treated time interval is from 0 to 5 seconds. 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.

5.2.1. Derivative of order 1/2, 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 [START_REF] Suarez | An eigenvector expansion method for the solution of motion containing fractional derivatives[END_REF] to reproduce the same results (Figure 1). Here, and in the followings, only the reference solutions are plotted; other solutions are very near when not practically identical to them.

Figure 2 shows the error indicator ε 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 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 η), which require different numbers of terms in the Adomian series.

As one would expect, the calculation times of the method based on the Grünwald-Letnikov formulation were not affected by the numerical parameters.

5.2.2. Derivative of order 1/2, harmonic excitation (problems 1 and 2)

Here, the reference solution is the solution by Adomian decomposition (Section 4.2). Results are shown in Figure 3.

Error indicators and CPU times are shown in Figure 4. Observations are the same as above: for a calculation time of the 10-second order, the Taylor-Adomian method is much more precise than the scheme with the Results are shown in Figure 5, while Figure 6 shows the values of the error indicator and CPU times. It is immediate to see the same as before:

the Taylor-Adomian method is clearly much more precise while keeping the calculation time low.

Simulating a longer time interval

As seen above, the terms of Equations ( 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 10 -10 . Figure 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 1 N N i=1

x ix steady i 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 seconds. Even problem 2 settles within a CPU time of less than 10 minutes on an average 64-bit system. Contrarily, the

A c c e p t e d m a n u s c r i p t

transient of problem 4 required somewhat more than 6 hours 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.

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

A c c e p t e d m a n u s c r i p t Grünwald-Letnikov derivative. 5 . 2 . 3 .

 523 Derivative of order 1/5, harmonic excitation (problems 3 and 4)To leave derivatives of order 1/2, another value has been chosen for α. The reference solution for these problems is the solution by Adomian decomposition (Section 4.2).

A c c e p t e d m a n u s c r i p t A c c e p t e d m a n u s c r i p t A c c e p t e d m a n u s c r i p tFigure 1 :

 1 Figure 1: Solution curves for problems 1u and 2u.

Figure 2 :

 2 Figure 2: Comparison of methods for problems 1u and 2u (α = 1/2, step function excitation).

Figure 3 :

 3 Figure 3: Solution of problems 1 and 2.

Figure 4 :

 4 Figure 4: Comparison of methods for problems 1 and 2 (α = 1/2, harmonic excitation).

Figure 5 :

 5 Figure 5: Solution of problems 3 and 4.

Figure 6 :

 6 Figure 6: Comparison of methods for problems 3 and 4 (α = 1/5, harmonic excitation).

Figure 7 : 23 AFigure 1 of 7 AFigure 2 of 7 AFigure 3 of 7 AFigure 4 of 7 AFigure 5 of 7 AFigure 6 of 7 AFigure 7 of 7

 72317273747576777 Figure 7: Calculation times for longer time intervals.

Acknowledgement

The author would like to thank Prof. László Szabó for his encouragement and most helpful comments.

X 0 0.25 0.25 0.25 0.25 0.25 0.25