Thomas J Waters

Stability of a 2-dimensional Mathieu-type system with quasiperiodic coefficients

Keywords: 2-d Mathieu, quasiperiodic, Floquet, harmonic balance, multiple scales NODY9599_source.tex, 30/10/2009, 10:15 p. 2

In the following we consider a 2-dimensional system of ODE's containing quasiperiodic terms. The system is proposed as an extension of Mathieu-type equations to higher dimensions, with emphasis on how resonance between the internal frequencies leads to a loss of stability. The 2-d system has two 'natural' frequencies when the time dependent terms are switched off, and it is internally driven by quasiperiodic terms in the same frequencies. Stability charts in the parameter space are generated first using numerical simulations and Floquet theory. While some instability regions are easy to anticipate, there are some surprises: within instability zones small islands of stability develop, and unusual 'arcs' of instability arise also. The transition curves are analyzed using the method of harmonic balance, and we find we can use this method to easily predict the 'resonance curves' from which bands of instability emanate. In addition, the method of multiple scales is used to examine the islands of stability near the 1:1 resonance.

Introduction

The study of systems of ode's with periodic coefficients arises naturally in many branches of applied mathematics. Aside from being interesting in their own right, in applications such systems arise when considering the orbital stability of a periodic solution to some dynamical system (see, for example, Betounes [START_REF] Betounes | Differential Equations: theory and applications[END_REF]). The archetypal Thomas J. Waters Dept. Applied Mathematics NUI Galway, Ireland E-mail: thomas.waters@nuigalway.ie equation for this class is the Mathieu equation, ẍ + (δ + ε cos t)x = 0 [START_REF] Betounes | Differential Equations: theory and applications[END_REF] where a dot denotes differentiation w.r.t. time. This equation, though linear, is sufficiently complicated by its explicit time dependence that closed form solutions do not exist. Instead, the common method of attack is to consider the stability of the origin using Floquet theory (or other methods). The (δ, ε) parameter space is clearly divided into regions of stability/instability (see Jordan and Smith [START_REF] Jordan | Nonlinear Ordinary Differential Equations[END_REF] for a detailed discussion). The appearance of these instability zones is intuitively understood in the following way: consider the time-dependent term in (1) as a perturbation. When ε = 0, x has periodic solutions with frequency √ δ (the 'natural' frequency of the system). As we switch on the perturbation, there will be a resonance between the 'driving' frequency, 1, and the natural frequency √ δ, which leads to instability.

There are many extensions to this equation, of which we describe only a few. Firstly the equation can be extended to many dimensions, by replacing x with a vector and δ and ε with matrices. This problem was analyzed in Hansen [START_REF] Hansen | Stability diagrams for coupled Mathieu-equations[END_REF] using Floquet theory, however (in contrast with the present work) the time dependent terms only depended on one frequency. Secondly the equation may be extended to contain nonlinear terms; for example El-Dib [START_REF] El-Dib | Nonlinear Mathieu equation and coupled resonnace mechanism[END_REF] consider an extension up to cubic order with each term containing subharmonics in the periodic term, and Younesian et al [START_REF] Younesian | Existence of periodic solutions for the generalized form of Mathieu equation[END_REF] consider also cubic order with two small parameters and terms in sin and cos.

Thirdly, and most importantly for the present work, we may develop the time dependent term to contain two periodic terms in different frequencies, and as the two frequencies need not be commensurable this system is referred to as the 'quasiperiodic Mathieu equation'. This type of system was developed and studied in great detail by Rand and coworkers: in Rand and Zounes [START_REF] Zounes | Transition curves for the quasiperiodic Mathieu equation[END_REF] and Rand et al [START_REF] Rand | A quasiperiodic Mathieu equation[END_REF] the following quasiperiodic Mathieu equation is examined ẍ + [δ + ε(cos(t) + cos(ωt))]x = 0 [START_REF] El-Dib | Nonlinear Mathieu equation and coupled resonnace mechanism[END_REF] and transition curves in the (δ, ω) parameter space are developed. Specific resonances are examined in Rand and Morrison [START_REF] Rand | 2:1:1 resonance in the quasiperiodic Mathieu equation[END_REF] and Rand et al [START_REF] Rand | 2:1:1 resonance in the quasiperiodic Mathieu equation[END_REF]. In Sah et al [START_REF] Sah | Autoparametric quasiperiodic excitation[END_REF] a form of the quasiperiodic Mathieu equation is derived and examined in the context of the stability of motion of a particle constrained by a set of linear springs, and finally in Rand and Zounes [START_REF] Zounes | Global behaviour of a nonlinear quasiperiodic Mathieu equation[END_REF] a nonlinear version of the quasiperiodic Mathieu equation is studied. In these papers a combination of numerical and analytical techniques are used to examine how resonances between the internal driving frequencies lead to a loss of stability in solutions, and it is very much in the spirit of these papers that we present the following work. We seek to extend consideration of quasiperiodic Mathieu-type systems to higher dimensions. The motivation for this is clear: as mentioned previously, linear systems with (quasi)periodic coefficients often arise when a nonlinear dynamical system is linearized around a (quasi)periodic solution, the so-called 'variational equations' [START_REF] Betounes | Differential Equations: theory and applications[END_REF]. As most dynamical systems of interest will have more than one degree of freedom (for example the motion of a particle in space under the influence of a collection of forces), the associated variational equations will also have more than one degree of freedom (unless only certain modes are being considered as in [START_REF] Sah | Autoparametric quasiperiodic excitation[END_REF]). We seek a 2-dimensional system containing more than one frequency and a small parameter ε to enable a perturbation analysis. However, it is beneficial to keep the number of free parameters low in order to make the behavior clear. With this in mind, we shall consider the following system:

ẍ + [α 2 + ε cos(βt)]x + ε cos(αt)y = 0 ÿ + ε cos(βt)x + [β 2 + ε cos(αt)]y = 0 (3) which in general is of the form d 2 X dt 2 + [P 0 + εP 1 (t)]X = 0 (4)
with X = (x, y) and

P 0 = α 2 0 0 β 2 , P 1 (t) = cos(βt) cos(αt) cos(βt) cos(αt) . ( 5
)
Here a dot denotes differentiation w.r.t. time, ε is small and α, β ∈ R.

We can predict to some degree how solutions of this system will behave: when ε = 0, the solutions are periodic in the natural frequencies α and β. As the perturbation, εP 1 (t)X, is switched on we would expect resonance between the natural and driving frequencies, as in the Mathieu equation [START_REF] Betounes | Differential Equations: theory and applications[END_REF]. As the natural and driving frequencies are the same, this means we expect instability to arise where α and β are commensurable. And this is to a degree what we do see, with two unexpected features. Firstly, curved 'arcs' of instability can be seen which clearly do not arise from integer ratios of the parameters. Secondly, within the bands of instability small islands of stability appear. These two features will be addressed in more detail in what follows.

While this system may seem a little artificial, we will show how it may arise as an approximate variational equation. Consider the following system of nonlinear ODE's:

ü = -α 2 u -uv, v = -β 2 v -uv. ( 6
)
This is an example of a system which is not dissipative but nonetheless is not Hamiltonian; as such it falls into the more general class of Louiville systems (see below). The author's direct experience with systems of this type is in the problem of the solar sail in the circular restricted 3-body problem (see for example the author and McInnes [START_REF] Waters | Periodic orbits above the ecliptic in the solar sail restricted 3-body problem[END_REF]), although other applications exist.

The origin of ( 6) is a fixed point, and letting u = ũ + x, v = ṽ + y we have the linear system

d 2 X dt 2 = -α 2 -v -u -v -β 2 -u u=ũ,v=ṽ X (7)
where X = (x, y). Letting (ũ, ṽ) = (0, 0) we find two linear oscillators which we write as ε cos(αt), ε cos(βt) (letting the amplitudes be equal and setting the phases to zero), and then letting ũ = ε cos(αt), ṽ = ε cos(βt) we get system (3) above. Note that the autonomous system (6) is not Hamiltonian and therefore [START_REF] Rand | 2:1:1 resonance in the quasiperiodic Mathieu equation[END_REF] is not symmetric; neither, for the same reason, is P 1 (t) in (3). Whether or not the system under consideration is Hamiltonian only becomes an issue in dimension 2 and above; in one dimension we can always find a (time dependent) Hamiltonian, a fact taken advantage of in Yeşiltaş and S ¸imşek [START_REF] Yeşiltaş | Multiple scale and Hamilton-Jacobi analysis of extended Mathieu equation[END_REF]. Note also that (u, v) = (ε cos(αt), ε cos(βt)) is not a solution to [START_REF] Rand | 2:1:1 resonance in the quasiperiodic Mathieu equation[END_REF], rather a solution to the linearized system [START_REF] Rand | 2:1:1 resonance in the quasiperiodic Mathieu equation[END_REF]. As such (3) is not a true variational equation, rather an approximation to one. We make this approximation, rather than seek more accurate solutions to [START_REF] Rand | 2:1:1 resonance in the quasiperiodic Mathieu equation[END_REF], so as to enable an analytical rather than purely numerical treatment. The Mathieu equation itself can be seen in this light, as an approximate variational equation to the system ü = -uu 3 .

The fact that (3) is only approximately a variational equation will be of relevance in what follows.

The rest of the paper is laid out as follows: In Section 2 we numerically integrate system (3) to examine the stability properties of various values of α and β. As computers can only handle finite precision we must choose rational values of α and β and therefore (3) is periodic, rather than quasiperiodic (a more appropriate term might be multiply-periodic). As the system is periodic, albeit of possibly very long period, we may use Floquet theory to enhance the numerical integration. We are able to generate stability charts in the (α, β) parameter space for various values of ε and we discuss the main features of same. The fact that we have lost the true quasiperiodic nature of the system in this section is made up for by the following: every irrational number has rational numbers arbitrarily close to it which thus give a very good approximation; the numerically generated stability charts allow us to isolate the interesting features of the system which informs later analysis; and the remaining analysis in the paper holds for arbitrary α and β and thus recaptures the quasiperiodic nature of the system.

In Section 3 we use the method of harmonic balance to find the transition curves in parameter space. Writing solutions in terms of truncated Fourier series we approximate the transition curves by the vanishing of the Hill's determinant. As the truncation size becomes large the matrices of vanishing determinant also become large making this method cumbersome. Also, there are convergence issues with the infinite Hill's determinant appropriate to (3) which we will elaborate on. However, we show that we can use the method of harmonic balance to quickly find approximations to transition curves by understanding the following mechanism of instability: when ε = 0 there are certain 'resonance curves' in parameter space (for the Mathieu equation (1) these are simply the points δ = n 2 /4 for n ∈ Z). As ε grows, these curves widen into bands of instability. The method of harmonic balance allows us to find these resonance curves quickly, and we demonstrate the usefulness of this technique on a second system closely related to [START_REF] Goldstein | Classical Mechanics[END_REF].

In Section 4 we use the method of multiple scales to examine in detail the fine structure near the α ≈ β resonance. From Section 1 we can see small islands of stability arising in the middle of a large band of instability, and we examine this feature more closely. We separate out the 'slow flow' using three time variables, and find we can analytically predict the appearance of these pockets of stability. Finally we give some discussion and suggestions for future work.

Numerical integration

We wish to fix the values of the parameters (α, β, ε) and numerically integrate the system (3) to examine the stability of the origin. In deciding what range of values to choose for the parameters, we note that (3) has the following scale invariance: consider the system ẍ + [(mα

) 2 + m 2 ε cos(mβt)]x + m 2 ε cos(mαt)y = 0 ÿ + m 2 ε cos(mβt)x + [(mβ) 2 + m 2 ε cos(mαt)]y = 0 (8)
By defining a new time coordinate t = mt we recover the original system (3). Thus the following parameter values are equivalent as regards stability:

(α, β, ε), (mα, mβ, m 2 ε).
This means we need only consider α and β on the unit square (for example), and the stability of other values of α and β are found by varying the value of ε. For example, the stability of system (3) when α = 3, β = 2, ε = 0.1 is the same as when α = 3/4, β = 1/2, ε = 0.0125.

As mentioned in the introduction, to integrate numerically we must use finite precision values of α and β and thus the system is (multiply-)periodic. As such, we may take advantage of Floquet theory. Writing (3) as a first order system, the fundamental matrix Φ(t) solves the following initial value problem (see, for example, Yakubovich and Starzhinskii [START_REF] Yakubovich | Linear differential equations with periodic coefficients[END_REF] for a detailed treatment of Floquet theory)

dΦ(t) dt = A(t)Φ(t), Φ(0) = I (9)
where

A(t) = ⎛ ⎜ ⎜ ⎝ 0 0 1 0 0 0 0 1 -[α 2 + ε cos(βt)] -ε cos(αt) 0 0 -ε cos(βt) -[β 2 + ε cos(αt)] 0 0 ⎞ ⎟ ⎟ ⎠ (10)
and A(t + T ) = A(t) for some fundamental period T .

The stability of the solutions to this system are determined by the monodromy matrix1 , given by Φ(T ), since it can be shown that

x(T ) y(T ) = Φ(T ) x(0) y(0) .
More precisely, stability is determined by the eigenvalues of Φ(T ), the so-called (characteristic) multipliers λ i . If |λ i | ≤ 1 the solutions remain bounded, implying stability, and if any |λ i | > 1 the solution is unbounded.

From the structure of ( 10) we can make some comments about the multipliers:

(a) ( b) ( c) ( d)
Fig. 1 Possible positions of multipliers in the complex plane with respect to the unit circle for a t-invariant system of order 4.

1. Louiville's theorem [START_REF] Jordan | Nonlinear Ordinary Differential Equations[END_REF] gives

det[Φ(T )] = λ 1 λ 2 λ 3 λ 4 = exp T 0 T r{A(t)}dt = 1.
This however does not mean that the eigenvalues must appear in reciprocal pairs, merely that all four multiplied together must be 1.

That they must appear in reciprocal pairs comes

from the t-invariant nature of the system [START_REF] Yakubovich | Linear differential equations with periodic coefficients[END_REF]. These are systems with the property

A(-t)G + GA(t) = 0 (11)
for some non-singular constant matrix G. Since A(t) is even, i.e. A(-t) = A(t), this holds with

G = I 0 0 -I , ( 12
)
the entries understood to be 2 × 2 matrices. In this case, we can show

Φ(T ) -1 = G.Φ(T ).G -1 , ( 13
)
in other words, the monodromy matrix and its inverse are similar and thus have the same eigenvalues; thus if λ i is a multiplier then so is 1/λ i .

As Φ(T ) is real, complex eigenvalues must appear in conjugate pairs, and since they must also be reciprocals this means complex eigenvalues must be on the unit circle (with one rare exception, see point 4 below). Therefore, complex multipliers represent stability, and real multipliers instability. 3. If the system (3) were a variational equation, that is a (nonlinear autonomous) system linearized around a periodic solution, we would expect a unit eigenvalue (see chap. 7 of Betounes [START_REF] Betounes | Differential Equations: theory and applications[END_REF]) and therefore (from the previous point) 2 unit eigenvalues. This would make life easier, as stability would be ensured if

T r{Φ(T )} ≤ 4 ( 14
)
and thus we would not need to calculate the multipliers themselves. Unfortunately we cannot say this is the case (see the Introduction), and so there is no guarantee of unit eigenvalues. As the system can be viewed as an approximation to a variational equation then we would expect (for small ε) to have two multipliers close to, but not equal to, one. 4. There are a number of possibilities for the positions of the multipliers in the complex plane, shown in Figure 1. The most common is for the set

1 + e, 1 1 + e , a + bi, a -bi (15)
with e ∈ R + and a 2 + b 2 = 1 (since λ = 1/λ). We may describe this as 'saddle-centre' (Figure 1

(b)).

Stability is given by 'centre-centre' configurations (Figure 1(c)), however we must also accept the possibility that the multipliers may not be on the unit circle (Figure 1(d)). This is possible when the eigenvalues of the centre-centre configuration come together and then bifurcate off the unit circle (a Krein collision). Note, the intermediate stage here means the eigenvalues are equal and this is a degenerate case; this will not be pursued further.

Based on these considerations, we will calculate the multipliers of the monodromy matrix and find their norms; if the norm is less than some cut-off value we say the system is stable, otherwise unstable.

We focus on (α, β) = (0, 1)×(0, 1) (as we may do due to the scale invariance of the system, mentioned previously). Splitting this interval into n equal segments, we let α = i/n and β = j/n for i, j = 1, . . . , n. The fundamental period of system ( 10) is given by

T = 2πn GCD[i, j]
where GCD is the greatest common divisor. We then integrate the system (9) for T units of time, evaluate Φ(T ), and calculate its eigenvalues.

The most important issue with regards accuracy is the long integration time. If i and j do not have a common divisor, than T = 2πn, which is large for n large. We may halve this integration time by taking advantage where the norm of the maximum multiplier is greater than 1.025 (and therefore the system is unstable), and black where the norm is less than 1.025 (the system being stable). The dashed circle shows the region of parameter space where the solutions in Figure 3 3 Long-time solutions for very close together points in α, β parameter space. On the left is the solution for α = 435/800, β = 425/800; this point is in the stability island indicated in Figure 2. On the right is the solution for α = 436/800, β = 425/800; this solution is outside the stability island. The integration times here are 1600π, and in both columns we present the x-and y-solution versus time and the x-y solution. The initial conditions for both are x(0) = 1, y(0) = 1, ẋ(0) = 0, ẏ(0) = 0. of the t-invariant nature of the system (see above). If A(t) is such that (11) holds, then we may show that

Φ(T ) = G.Φ(T /2) -1 .G.Φ(T /2)
since G = G -1 . Hence we need only integrate over one half period to calculate Φ(T ). This improves the accuracy greatly, with one caveat: we must invert Φ(T /2). However, since Louiville's theorem holds for all t, this means the determinant of Φ(T /2) should be one, and thus inverting may be done accurately.

In fact, we may use the condition det[Φ(t)] = 1 to monitor the accuracy of the integration, however this proves very costly with regards integration time. Indeed, for some values of α, β and ε we would expect the solution to be grossly unstable. Monitoring the accuracy of such solutions is inefficient as our concern is the stability-instability divide, rather than how large these very large multipliers may be. As such, we are satisfied that if det[Φ(t)] grows very large this is due to the huge growth in very unstable combinations of α and β; we will not waste effort with the accuracy of these solutions.

In Figure 2 we give the stability diagram for system (3) with ε = 0.1. Here, we set n = 800 and thus divide the α, β unit square parameter space into 800 2 data points. For each we calculate the maximum multiplier norm, with white regions being where the max norm is larger than 1.025, black regions where the max norm is less than 1.025. We may make some observations:

-The most striking features are the large, linear bands of instability. These are aligned along lines of constant slope in the parameter space, in other words where α/β =constant. These are clearly the dominant, low order resonances between α and β, and the lower the order of resonance the wider the band of instability. There are also higher order resonances, and these are of narrower width. This behavior is very much what we would expect, based on the classical Mathieu equation and its variations. -There are two interesting and unexpected features.

The first is the large 'arc' of instability seen in the corners. This curved band of instability is not predicted by simple resonance between α and β. -Second and perhaps more interesting are the small islands of stability which appear in the large α ≈ β resonance band. We would fully expect the bands of instability to grow with ε, and this we do observe. However, for small pockets of stability to develop when surrounding them is very large regions of instability is most unexpected. These islands are precarious: we show in Figure 3 the solution for two very nearby points in parameter space, one inside a small island of stability and one outside.

-The corresponding plot for other values of ε can be inferred from Figure 2 in the following way. Larger values of ε generate stability charts which 'zoom in' towards the origin of Figure 2, smaller values 'zoom out'; this is due to the scale invariance of the system as described at the beginning of this section. For example, the stability chart for α, β on the unit square and ε = 0.4 is the same as the lower left quadrant of the same chart for ε = 0.1, since the parameter values (α, β, ε) = (1/2, 1/2, 0.1) have the same stability properties as (1, 1, 0.4). This is why there is a cluster of unstable points near the origin: moving towards the origin is equivalent to increasing ε; the resonance bands widen and overlap. As we decrease ε, regions beyond the unit square in Figure 2 are drawn in; since the resonance bands narrow as we move away from the origin this is equivalent to saying they narrow as ε is decreased, as expected, and the chart is increasingly made up primarily of stable points. The exception is the α ≈ β resonance band; it does not appear to taper as we move along it away from the origin. Some insight into this is gained in Section 4.

While the plots generated using numerical integration are informative, they give us no feel for why certain bands of instability develop or which are dominant. Nor do they give us an analytic expression for the values at which the system bifurcates from stable to unstable. In the next two sections, we use analytical techniques which do just this, and what's more recapture the true quasiperiodic nature of the original system.

Transition curves using the method of harmonic balance

In Floquet theory, t-invariant systems must have multipliers on the unit circle or real line (with the rare exception of Figure 1(d)). As the system passes from stable to unstable, the multipliers must pass through the point (1, 0) or (-1, 0) in the complex plane. These multipliers represent solutions which have the same period, or twice the period, of the original system. This means solutions transiting from stable to unstable have a known frequency and thus can be written as a Fourier series. For this series to be a solution to the system puts a constraint on the parameters involved, and this in turn tells us the transition curves in parameter space.

In quasiperiodic systems, such as the one under consideration in this paper, which we give again as

ẍ + [α 2 + ε cos(βt)]x + ε cos(αt)y = 0 ÿ + ε cos(βt)x + [β 2 + ε cos(αt)]y = 0, (16)
Floquet theory does not hold. However, based on the discussion above we can use the following ansatz: solutions on transition curves from regions of stability to instability in parameter space have the same quasiperiodic nature as the system itself. This was justified and used by Rand et al in [START_REF] Zounes | Transition curves for the quasiperiodic Mathieu equation[END_REF][START_REF] Rand | A quasiperiodic Mathieu equation[END_REF], and while lacking in rigorous theoretical proof we find this assertion works well in predicting the transition curves in parameter space. Since any quasiperiodic function can be written as an infinite Fourier series in its frequencies (see, for example, Goldstein et al [START_REF] Goldstein | Classical Mechanics[END_REF]), we write the solutions to system (16) in the following form,

x y = ∞ n=-∞ ∞ m=-∞ A nm B nm exp it 2 (nα + mβ) , ( 17
)
where the 2 appears in the exponent to account for solutions with twice the period of the system (corresponding to the -1 multiplier). We substitute (17) into (16), use the identity cos(θt) = 1 2 (e iθt +e -iθt ), and shift indices to cancel out the exponential terms. What remains are the following two recurrence relations:

α 2 - (nα + mβ) 2 4 A nm + ε 2 (A n,m-2 + A n,m+2 + B n-2,m + B n+2,m ) = 0, ( 18
)
β 2 - (nα + mβ) 2 4 B nm + ε 2 (A n,m-2 + A n,m+2 + B n-2,m + B n+2,m ) = 0, (19)
with n, m = -∞, . . . , ∞. This is a set of linear, homogeneous equations in A nm , B nm of the form Cx = 0, and for non-trivial solutions to exist we require det(C) = 0; these are the well know infinite Hill's determinants.

In practice, we truncate at some order N so as n, m = -N, . . . , N. Truncating at larger values of N will produce more accurate determinants, provided the determinant converges. However, close inspection of the recurrence relations (18,19) reveals that this is an issue for the system in question.

Writing the coefficient of A nm in (18) as γ nm and dividing across gives the following:

A nm + ε 2γ nm (A n,m-2 + . . .) = 0 (20)
with a similar expression for (19). In most systems, as N → ∞ the 'off-diagonal' terms here grow progressively smaller, ensuring convergence. However, we see that if we let, for example, n = 2, m = 0 then

γ 20 = 0.
Thus the infinite determinant does not converge. This is a problem particular to systems with a driving frequency that is equal to a natural frequency. As a way around this problem we notice that a solution with period 2T also has period 4T , and indeed 6T, 8T, . . . and so on. Thus we may write

x y = ∞ n=-∞ ∞ m=-∞ A nm B nm exp it 2M (nα + mβ) , (21)
where M = 1, 2, 3, . . .. In this case we see

γ nm = α 2 - (nα + mβ) 2 (2M ) 2
and the singular term does not appear in the Hill determinant until N = 2M . Thus we may work up to order N = 2M -1 and calculate the determinant. It will be a function of the parameters (α, β, ε) and by setting it equal to zero we may therefore find analytic approximations to the transition curves in parameter space.

For example, with M = 3 we generate the matrix appropriate to 5th order, find its determinant, and set it equal to zero. Factorizing the determinant we see it contains terms like the following:

(4α 2 -β 2 + ε)(4α 2 -β 2 -ε) = 0, (9α 4 -82α 2 β 2 + 9β 4 -64e 2 ) × (9α 4 -82α 2 β 2 + 9β 4 -16e 2 ) = 0.
The first two clearly indicate a band of instability with width O(ε) around the 2:1 resonance curve, which grows to either side of the resonance as ε grows. The second two terms indicate a band of width O(ε 2 ) around the 3:1 resonance curve, which grows away from the resonance as ε grows. This implies that stability for resonant values of α and β is not as straightforward as we would initially expect: for α/β = a/b with a, b integers the band of instability may move away from the resonance curve as ε grows and thus commensurable values of α, β will in fact be stable.

We give the zeroes of some of the main factors of the determinant in Figure 4. We can see in the figure some of the main features we expect from our numerical simulations, however some features are notable absent; foremost among these is the 1:1 resonance. Going to higher orders will capture more of the transition curves, however the dimension of the problem quickly becomes very large. If we truncate the series at N = 2M -1 we must find the determinant of a 2(2N + 1) 2 square matrix, which for the example described above (M = 3) is 242 × 242. Calculating the determinant can be facilitated in the following way: these large matrices are very sparse as is usual for recurrence relations. We may separate out 2M submatrices, each describing the T, 2T, . . . , 2M T solutions, as in

C = C 1 C 2 . . . C 2M ,
and set the determinants of each to zero. However, the improvements in accuracy from increasing M do not justify the increasingly cumbersome and lengthy matrices. Thus for this system it is worthwhile to study the resonances of greatest import individually; this we do in the next section. ******* Aside from these complications, we note that the method of harmonic balance can be used very quickly and directly to predict, in broad terms, the locations of the bands of instability. As mentioned previously, we may understand the mechanism of instability so: resonance curves in the (α, β) parameter space grow with ε into bands of instability due to resonance between the natural and driving frequencies. Predicting the location of said resonance curves will give approximately the regions of parameter space which will become unstable. As these resonance curves are the bands of instability in the limit of ε → 0, we may simply let ε = 0 in the recurrence relations given in (18,19). This leads to the following definition of the resonance curves: or

α 2 - (nα + mβ) 2 4 = β 2 - (nα + mβ) 2 4 = 0, (22)
α = ±mβ 2 ± n , β = ±nα 2 ± m (23)
for n, m = -∞, . . . , ∞, unless n(m) = ±2 in which case β(α) = 0. These are precisely the straight lines in the α, β plane where the frequencies are in ratio.

As an illustration of how useful this technique can be we can examine the following problem, closely related to (16),

ẍ + [α + ε cos(βt)]x + ε cos(αt)y = 0 ÿ + ε cos(βt)x + [β + ε cos(αt)]y = 0. ( 24
)
The time dependent terms here are of the same form as in ( 16), and thus we write solutions on the transition curves as in (17). Now the recurrence relations, in the limit of ε → 0, yield the following resonance curves:

α - (nα + mβ) 2 4 = β - (nα + mβ) 2 4 = 0, ( 25
)
for n, m = -∞, . . . , ∞ (we see that in this case there is no issue with convergence of the infinite Hill's determinants, as the driving frequencies do not equal the natural frequencies). The resonance curves for the first few values of n, m are given in Figure 5, which can be compared with numerically generated stability charts (using methods already discussed in Section 2) for nonzero ε as shown in Figure 6. We see that the resonance curves as predicted by (25) give a good approximation to the location of the bands of instability in parameter space, with very little effort. 24) with ε = 0.1 and the α, β unit square divided into an 800 × 800 grid. White represents points where the norm of the maximum multiplier is greater than 1.025, and black where the norm is less than 1.025. We see the strong correspondence between bands of instability in this plot and the resonance curves in Figure 5.

In fact, we may generalize this technique to arrive at the following proposition:

Proposition: For the system

d 2 X dt 2 + (P 0 + εP 1 (t)) X = 0 (26) with X ∈ R n , P 0 = diag(λ 2 1 , λ 2 2 , . . . , λ 2 n ) and P 1 (t) of the form P 1 (t) = ⎛ ⎜ ⎝ cos(ω 1 t) . . . cos(ω n t) . . . . . . . . . cos(ω 1 t) . . . cos(ω n t) ⎞ ⎟ ⎠ ,
bands of instability will grow with ε from resonance curves given by

λ 2 j - (n.ω) 2 4 = 0, j = 1, . . . , n,
where n is a vector of length n whose elements are all integers in the range -∞, . . . , ∞, that is n ∈ Z n .

Proof: On transition curves, the solution has the same quasiperiodic structure as P 1 (t) and so we can write it as a Fourier series of the form

X(t) = n A n exp it 2 (n.ω)
where A n is an n-dimensional vector of Fourier coefficients. Subbing this into (26) we find

- (n.ω) 2 4 n A n exp it 2 (n.ω) + (P 0 + εP 1 (t)) n A n exp it 2 (n.ω) = 0.
The term in ε in this equation is complicated and would involve shifting indices were we looking for the infinite Hill's determinant to determine the transition curves; however we seek only the resonance curves which are found by letting ε = 0. Canceling the exponentials

λ 2 j - (n.ω) 2 4 = 0, n ∈ Z n .
Now we may simply plot these curves for some elements of n as in Figure 5.

We note this method predicts the resonance points for the Mathieu equation ( 1) to be δ-n 2 /4 = 0, and the resonance curves for the quasiperiodic Mathieu equation (2) to be δ -(n + mω) 2 /4 = 0 (as in [START_REF] Zounes | Transition curves for the quasiperiodic Mathieu equation[END_REF]), with n, m ∈ Z.

The multiple scales method near the 1:1 resonance

Based on the previous section we see that the method of harmonic balance is unable to describe behavior near the 1 : 1 resonance, that is close to α = β. We use instead the method of multiple scales, which is a perturbative technique that allows us to derive analytic representations near resonant values of α and β.

The typical way to proceed is to expand α, β, x and y in powers of ε, then define different time scales and separate out the behaviour on slower time scales. We could let

α = α 0 + εα 1 + ε 2 α 2 , β = α 0 + εβ 1 + ε 2 β 2
so α and β only differ at first order. Then, the differential equations would contain terms involving α 0 t + α 1 εt + α 2 ε 2 t which suggests we define three time scales, t, τ = εt, T = ε 2 t (T here is not to be confused with the fundamental period T from Section 2). Note the ε 2 time scale is necessary as there are no secular terms at first order. Now we may expand x and y in powers of ε where each term is a function of these three times, i.e.

x = x 0 (t, τ, T ) + εx 1 (t, τ, T ) + ε 2 x 2 (t, τ, T ) + . . . (27)
and similarly for y. Subbing into the differential equations and collecting powers of ε we find a set of equations with inhomogeneous terms involving partial derivatives of lower order solutions. Removing secular terms gives differential equations in the slow-time dependence of low order terms. However, we find that the crucial system of equations only describes flow on the 'very' slow time T , and is independent of α 1 and β 1 (see the Appendix for details).

Therefore, it is sufficient to define

α = α, β = α + ε 2 β 2 . ( 28
)
Now there is only one slow time, T = ε 2 t, however we are still expanding x and y in powers of ε (rather than ε 2 ). To second order we find the following system of equations:

ε 0 : ∂ 2 x 0 ∂t 2 + α 2 x 0 = 0, ∂ 2 y 0 ∂t 2 + α 2 y 0 = 0, ε 1 : ∂ 2 x 1 ∂t 2 + α 2 x 1 = -x 0 cos(αt + β 2 T ) -y 0 cos(αt), ∂ 2 y 1 ∂t 2 + α 2 y 1 = -y 0 cos(αt + β 2 T ) -y 0 cos(αt), ε 2 : ∂ 2 x 2 ∂t 2 + α 2 x 2 = -2 ∂ 2 x 0 ∂t∂T -x 1 cos(αt + β 2 T ) -y 1 cos(αt), ∂ 2 y 2 ∂t 2 + α 2 y 2 = -2 ∂ 2 y 0 ∂t∂T -y 1 cos(αt + β 2 T ) -y 1 cos(αt) -2αβ 2 y 0 ,
where x 0 = x 0 (t, T ) and so on. Solving one at a time, we find at zero order:

x 0 = A(T ) cos(αt) + B(T ) sin(αt), (29a)
y 0 = C(T ) cos(αt) + D(T ) sin(αt). (29b)
There are no secular terms on the right hand side at first order so we find

x 1 = a(T ) cos(αt)+b(T ) sin(αt)+ A 6α 2 cos(2αt+β 2 T )+. . .
and a similar expression for y 1 . Now when we look at second order, we see there are secular terms in the right hand side (that is, terms involving cos(αt) and sin(αt)). Setting these terms to zero results in four linear differential equations in the slow time dependence of x 0 and y 0 , that is the four functions A(T ), B(T ), C(T ) and D(T ) defined in (29). Letting Y = (A, B, C, D) we may write these four equations as dY /dT = Q(T )Y where

Q = ⎛ ⎜ ⎜ ⎝ -sin(β 2 T ) + 3 sin(2β 2 T ) -2 + cos(β 2 T ) + 3 cos(2β 2 T ) 7s i n ( β 2 T ) 1 + c o s ( β 2 T ) 2 + 5 cos(β 2 T ) + 3 cos(2β 2 T ) -7 sin(β 2 T ) -3 sin(2β 2 T ) 5+5c o s ( β 2 T ) s i n ( β 2 T ) -sin(β 2 T ) + 3 sin(2β 2 T ) -2 + cos(β 2 T ) + 3 cos(2β 2 T ) 7s i n ( β 2 T ) 1+2 4 α 3 β 2 + cos(β 2 T ) 2 + 5 cos(β 2 T ) + 3 cos(2β 2 T ) -7 sin(β 2 T ) -3 sin(2β 2 T ) 5 -24α 3 β 2 + 5 cos(β 2 T ) s i n ( β 2 T ) ⎞ ⎟ ⎟ ⎠ 0.05 0.1 0.15 0.2 Μ 1 5 10
Fig. 7 The largest multiplier of system (30) for μ small.

The stability of this system will determine the stability of the original system.

This system is still Louiville (T r(Q) = 0) and tinvariant (as defined in [START_REF] Yakubovich | Linear differential equations with periodic coefficients[END_REF] with G = diag(1, -1, 1, -1)). A nice simplification happens if we define a new time coordinate t * = β 2 T , and we may write the system in the form

dX dt * = [A + μ(A 0 + A 1 cos(t * ) + A 2 cos(2t * ) +B 1 sin(t * ) + B 2 sin(2t * ))] X ( 30
)
where μ = 1/(24α 3 β 2 ) and the constant coefficient matrices are given by

A = ⎛ ⎜ ⎜ ⎝ 0 0 0 0 0 0 0 0 0 0 0 1 0 0 -1 0 ⎞ ⎟ ⎟ ⎠ , A 0 = ⎛ ⎜ ⎜ ⎝ 0 -2 0 1 2 0 5 0 0 -2 0 1 2 0 5 0 ⎞ ⎟ ⎟ ⎠ ,
with A 1 , A 2 , B 1 , B 2 in the Appendix. This is a one parameter 2π-periodic system, with the μ parameter multiplying terms very much like a Fourier series. It is straightforward to analyze the stability of this system using the methods already outlined, and we show the maximum multiplier in Figure 7.

The appearance of a small region of stability can be understood in the following way. In system (30), when μ = 0 the system has natural frequency 1. As we switch on μ, the driving frequency is also 1 and resonance causes the system to be unstable. However, as μ grows, the natural frequencies are now the eigenvalues of A + μA 0 , which become sufficiently distinct from the driving frequency for stability to resume. This is transient however and as the 'perturbation' μ grows the system becomes unstable again. This window of stability occurs for (μ -=) 0.18106 μ 0.189 (= μ + ), which tells us the values of β 2 where stability occurs and this in turn gives us the following band of stability in (α, β) parameter space

α ± ε 2 1 24α 3 μ + < β < α ± ε 2 1 24α 3 μ -
This result is in excellent agreement with the numerical results generated earlier, as shown in Figure 8, at least for values of α, β over 0.5. As discussed previously, this is because small values of α, β are equivalent to large values of ε and the multiple scales method would not be expected to work well in this regime. Note that the symmetric nature of the system is recovered by the expansion

β = β, α = β + ε 2 α 2 , ( 31
)
analogous to (28), which results in the same stability region as in Figure 8 reflected about the line α = β.

Conclusions

In this paper we have begun to analyze the rich dynamical structure of two dimensional systems with quasiperiodic coefficients. Through a combination of numerical and perturbative techniques we are able to predict and determine the transition curves in parameter space which demarcate the regions leading to stable and unstable solutions. In particular, we find that using the method of harmonic balance we can quickly approximate the resonance curves in parameter space from which bands of instability emanate, with the caveat that the band of instability may spread out on either side of the resonance curve or move to its side. The multiple scales method has proven to be particularly informative, if a little cumbersome, for examining the vicinity of specific resonances. There are three obvious issues which need to be understood further with regards the particular system under consideration in this work. The first is the 'arcs' of instability which are commented on in §2 and can be seen in Figure 2. These arcs do not arise in the harmonic balance analysis, nor do they correspond to a linear relationship between the frequencies and thus the method of multiple scales would be inappropriate. The exact mechanism which leads to the development of these arcs is not understood, and is worthy of further investigation.

The second issue is the possible resonances between first order terms in the expansion of α, β as described in the Appendix. Perhaps these could explain why the multiple scales method predicts a continuous band of stability near the 1:1 resonance, whereas numerically we observe a broken series of islands of stability (see Figure 2 and8). Finally, all bands of instability visible in Figure 2 taper as we increase α and β, as predicted by the harmonic balance method, with the exception of the large region of instability surrounding the 1:1 resonance. While the multiple scales analysis could be interpreted as telling us that the true band of instability around α = β is in fact very narrow and exists between the predicted bands of stability, as shown in Figure 8, it still leaves an open question as to the origin of this large region of instability.

There are numerous extensions to the system in question, however it would be wise to keep the number of free parameters small to enable an informative analysis. For example, a general form could consist of two natural frequencies and four driving frequencies. While this would undoubtedly lead to exceedingly rich stability properties, six free parameters would be difficult to present in a clear fashion. A 'detuning' parameter in the driving frequencies would instead lead to three parameters and would provide an interesting extension. There are of course other extensions involving nonlinear terms, more dimensions and so on. Finally, there are numerous applications involving quasiperiodic oscillations about fixed points whose analysis would benefit from the present work.

Appendix

As is typical in multiple scales analysis, the equations at high orders are quite lengthy and involve very many terms; as such we will only give an overview to justify ignoring first order terms in α, β to arrive at (28).

Let us write

α = α 0 + εα 1 + ε 2 α 2 , β = α 0 + εβ 1 + ε 2 β 2
and expand x, y in the three time scales t, τ and T , as in (27). The second time derivative becomes (up to second order in ε)

d 2 dt 2 = ∂ 2 ∂t 2 + ε 2 ∂ 2 ∂t∂τ + ε 2 2 ∂ 2 ∂t∂T + ∂ 2 ∂τ 2 ,
leading to the following system of equations in x:

ε 0 : ∂ 2 x 0 ∂t 2 + α 2 0 x 0 = 0 ε 1 : ∂ 2 x 1 ∂t 2 + 2 ∂ 2 x 0 ∂t∂τ + α 2 0 x 1 + 2α 0 α 1 x 0 + x 0 cos(α 0 t + β 1 τ + β 2 T ) + y 0 cos(α 0 t + α 1 τ + α 2 T ) = 0 ε 2 : ∂ 2 x 2 ∂t 2 + 2 ∂ 2 x 1 ∂t∂τ + 2 ∂ 2 x 0 ∂t∂T + ∂ 2 x 0 ∂t 2 + α 0 x 2 + 2α 0 α 1 x 1 + (α 2 1 + 2α 0 α 2 )x 0 + x 1 cos(α 0 t + β 1 τ + β 2 T ) + y 1 cos(α 0 t + α 1 τ + α 2 T ) = 0
For the corresponding system in y, simply swap x with y, and α 1 , α 2 with β 1 , β 2 . We examine the equations one order at a time. At zero order, the solutions are x 0 = a(τ, T ) cos(α 0 t) + b(τ, T ) sin(α 0 t), y 0 = c(τ, T ) cos(α 0 t) + d(τ, T ) sin(α 0 t).

At first order, which we write as At second order, which again we write as ∂ 2 x 2 ∂t 2 + α 0 x 2 = inhomogeneous terms, we use the zero and first order solutions to calculate the right hand side and look for secular terms, that is terms multiplying cos(α 0 t), sin(α 0 t). The full expression for the right hand side is very lengthy and complicated, and involves multiplying many cos and sin terms of different arguments together, however we are only interested in those with α 0 in the t-dependence and so can ignore those containing 2α 0 or 3α 0 . This gives us four equations: the coefficients of cos(α 0 t), sin(α 0 t) in the x-equation, and the same for the y-equation, which must be set equal to zero to prevent resonance in x 2 , y 2 . These four equations are of the form

-2α 0 ∂Q ∂τ -2α 0 α 1 P + f 1 (t, τ, T ) = 0, (33a)
2α 0 ∂P ∂τ -2α 0 α 1 Q + f 2 (t, τ, T ) = 0, (33b)
and again in R, S and f 3 , f 4 with α 1 replaced by β 1 . These define the τ -dependence of P, Q, R, S (from the first order solutions x 1 , y 1 in (32)), and we see these will be oscillatory if we remove the cos(α 1 τ ),sin(α 1 τ ) terms in f 1 , f 2 , and the cos(β 1 τ ),sin(β 1 τ ) terms in f 3 , f 4 . This leads to eight equations involving the eight functions describing the T -dependence of the zero order solutions, that is

A, A , B, B , C, C , D, D

and their derivatives with respect to T . We may write this as a first order system dZ/dT = Λ(T )Z, where Λ has block form

Λ(T ) = Λ 1 Λ 2 (T ) Λ 3 (T ) Λ 4 ,
with Λ being Louiville (vanishing trace) and t-invariant, as defined in [START_REF] Yakubovich | Linear differential equations with periodic coefficients[END_REF] with

G = diag(1, -1, -1, 1, 1, -1, -1, 1).
The stability of our original system now hangs on the stability of this one, and crucially the parameters α 1 , β 1 do not appear in the matrix Λ. Thus the difference between α and β at first order is irrelevant and we may set α 1 = β 1 = 0.

We note that there are a number of possible resonances in the first order terms in α, β which we have ignored. The f i functions described above contain terms like, for example, cos(2α

1 τ + β 1 τ + α 2 T -β 2 T ). ( 34
)
In general, this will not be secular in equations (33). If however it was the case that β 1 = -α 1 then this would be secular. There are six combinations which lead to additional secular terms, they are

α 1 = ±β 1 , α 1 = ±3β 1 , 3α 1 = ±β 1 .
Each would need to be considered in turn, and would lead to a different coefficient matrix Λ given above. We will not pursue this issue further in the present work. Finally, we give the coefficient matrices from (30),

1 ΒFig. 2

 12 Fig.2Stability diagram for system (3) with ε = 0.1 and the α, β unit square divided into an 800 × 800 grid. White represents points where the norm of the maximum multiplier is greater than 1.025 (and therefore the system is unstable), and black where the norm is less than 1.025 (the system being stable). The dashed circle shows the region of parameter space where the solutions in Figure3are taken from.

  Fig.2Stability diagram for system (3) with ε = 0.1 and the α, β unit square divided into an 800 × 800 grid. White represents points where the norm of the maximum multiplier is greater than 1.025 (and therefore the system is unstable), and black where the norm is less than 1.025 (the system being stable). The dashed circle shows the region of parameter space where the solutions in Figure3are taken from.

Fig. 4

 4 Fig.[START_REF] Hansen | Stability diagrams for coupled Mathieu-equations[END_REF] Some transition curves on the unit square for ε = 0.1 using the method of harmonic balance. The curves are where the determinant of the coefficient matrix C vanishes.

1 ΒFig. 5

 15 Fig. 5 Resonance curves in the α, β plane for the system (24), with |n|, |m| ≤ 2.

1 ΒFig. 6

 16 Fig.[START_REF] Rand | 2:1:1 resonance in the quasiperiodic Mathieu equation[END_REF] Stability diagram for system (24) with ε = 0.1 and the α, β unit square divided into an 800 × 800 grid. White represents points where the norm of the maximum multiplier is greater than 1.025, and black where the norm is less than 1.025. We see the strong correspondence between bands of instability in this plot and the resonance curves in Figure5.

  P 0 A n = 0 and since P 0 = diag(λ 2 i ) for i = 1, . . . , n this leaves us with the n equations

1 ΒFig. 8

 18 Fig.8The multiple scales prediction for a band of stability overlaid onto the stability chart generated numerically in Section 2.

∂ 2 x 1

 1 ∂t 2 + α 0 x 1 = inhomogeneous terms, the right hand side includes2α 0 sin(α 0 t) ∂a ∂τ α 1 b + 2α 0 cos(α 0 t) -∂b ∂τ α 1 ain the x-equation, and similarly for c, d in the y-equation.Setting these secular terms equal to zero gives us the τ -dependence at zero order, i.e.a(τ, T ) = A(T ) cos(α 1 τ ) + A (T ) sin(α 1 τ ), b(τ, T ) = B(T ) cos(α 1 τ ) + B (T ) sin(α 1 τ ),and similarly for c, d. The solution at first order will then involve the homogeneous solutionx 1 = P (τ, T ) cos(α 0 t) + Q(τ, T ) sin(α 0 t), (32a)y 1 = R(τ, T ) cos(α 0 t) + S(τ, T ) sin(α 0 t),(32b)and the inhomogeneous solution, made up of sin and cos terms with arguments (α 0 ± α 0 )t + (α 1 /β 1 )τ + (α 2 /β 2 )T which we will not list in full (by α i /β i we mean linear combinations of the form k 1 α i + k 2 β i with k ∈ Z and |k 1 | + |k 2 | ≤ 3).

In some of the literature, A(t) is known as the monodromy matrix and Φ(T ) as the matrizant, propagator or characteristic matrix.

NODY9599_source.tex; 30/10/2009; 10:15 p. 9


A Halanay
email: ha-lanay@mathem.pub.ro

A Ionit ¸ȃ

C-A Safta

HOPF BIFURCATIONS THROUGH DELAY IN PILOT REACTION IN A LONGITUDINAL FLIGHT

Keywords: delay differential equations, Hopf bifurcation, Lyapunov coefficient, pilot induced oscillations AMS 2000 subject classification: 34K18, 34K20, 34G15

The paper is devoted to the study of pilot induced oscillations in the landing transition between the approach task and flare to touch-down. These oscillations are proved to appear in a longitudinal flight model when the delay in pilot's reactions exceeds a certain threshold for which the stability of equilibria is lost and a Hopf bifurcation appears. The formulae needed to compute the Lyapunov coefficient and an approximation of the solution are developed for the delay differential equations that model the pilot-vehicle interaction in landing task. These are applied for a concrete model.

§1. Introduction

The so-called Pilot Induced Oscillations (PIO) phenomena are now recognized to belong to a certain subclass of Aircraft Pilot Coupling (APC) (see [START_REF]Flight Control Design -Best Practices[END_REF]). The oscillatory APC event is defined as "inadvertent, unwanted aircraft motion that is associated with anomalous interactions between the aircraft and the pilot" ( [START_REF]Aviation Safety and Pilot Control -Understanding and Preventing Unfavorable Pilot-Vehicle Interactions[END_REF]). In what follows a generic model (see [START_REF] Forssel | The Aero-Data Model In a Research Environment (ADMIRE) for Flight Control Robustness Evolution[END_REF]) will be considered that outlines the interaction between effective aircraft dynamics and pilot's characteristics in the longitudinal motion since PIO are dominant in terminal flight conditions. Famous PIO in landing phase were reported for YF22, Saab Grippen, C-17A, Airbus A321, Boeing 777 aircraft, some of which are using advanced highly augmentation flight control systems. Yet, the Hopf bifurcations of equilibria and the study of the stability of the resulting limit cycle taking as parameter the time delay in pilot's action were not considered until much later (see [START_REF] Mehra | Bifurcation and limit cycle analysis of nonlinear pilot[END_REF], [START_REF] Ph | Stability Analysis Through Bifurcation Theory[END_REF], [START_REF] Garteur Fm | AG12 -Analysis of nonlinear pilot-vehicle systems using modern control theory[END_REF], [START_REF] Ionit | Delay Induced Oscillations[END_REF]).

In the study of PIO based on the qualitative theory of delay differential equations, two points of view are met: one is pessimistic with respect to the use of such a model due to its inherent complexity: functional initial conditions, an infinite number of roots for the characteristic equation, etc. A second point of view, more optimistic, hopes to use such models in order to choose closed-loop system's parameters such that the bifurcating value be high with respect to certain requirements (see [START_REF] Ionit | Delay Induced Oscillations[END_REF]).

The time-delay model that will be studied below is a conventional one and can be seen as an equivalent system obtained from the high order augmentation dynamics of an airplane.

The paper is organized as follows: in Section 2 the model of the pilot-vehicle interaction in landing tasks is introduced. In Section 3 the system of equations for equilibria is analyzed in connection with Hopf bifurcation theorem in [START_REF] Hale | Theory of Functional Differential Equations[END_REF]. Section 4 is devoted to the study of the stability of the limit cycle provided by the Hopf bifurcation theorem. The first Lyapunov coefficient, as well as others constants involved in the description of the limit cycle, depend, of course, on the equilibrium point where bifurcation appears. In Section 5 a case study is presented. Some concluding remarks are given in Section 6. §2. The Pilot-Vehicle Model in the Landing Task

The ADMIRE generic model (see [START_REF] Forssel | The Aero-Data Model In a Research Environment (ADMIRE) for Flight Control Robustness Evolution[END_REF]) is a single engine aircraft with a deltacanard configuration. The coefficients in the nonlinear aerodynamic model are calculated on the basis of the Generic Aerodynamic Model (GAM). The aircraft dynamics is completed by a flight control system in order to provide stability and affordable handling qualities within the operational envelope. The longitudinal controller provides the pitch rate control and angle of attack control. It consists of two parts

u(t) = u 1 (t) + u 2 (t -τ )
where u 1 is the contribution of the flight control system with stability augmentation and u 2 is the contribution of pilot's action. So

u 1 (t) = k α α(t) + k q q(t), u 2 (t) = kθ(t -τ ).
The model for pilot's action is in the frame of classical crossover model: pure pilot gain k and pilot delay τ . The mathematical model for a longitudinal landing flight becomes

(2.1) α = m 11 α + m 12 q + c cos θ + b 1 k θ(t -τ ) q = m 21 α + m 22 q + c m 0 cos θ -c c 1 sin θ + b 2 k θ(t -τ ) θ = q.
The state vector consists of the incidence angle α, the pitch attitude θ and the pitch rate q. The coefficients are supposed constant but depend on flight conditions and aerodynamic characteristics such as mass or moment of inertia of the generic aircraft.

S 3. Hopf bifurcation of equilibria

Equilibria for (2.1) are obtained from the system (3.1)

q = 0 m 11 α + b 1 k θ + c cos θ = 0 m 21 α + b 2 k θ + c m 0 cos θ -c c 1 sin θ = 0.
Let (α 0 , 0, θ 0 ) be an equilibrium point for (2.1). Perform a translation into zero through

(3.2) α 1 = α -α 0 , q 1 = q, θ 1 = θ -θ 0 System (2.1) becomes (3.3) α1 = m 11 α 1 + m 12 q 1 + c cos θ 0 cos θ 1 - -c sin θ 0 sin θ 1 + b 1 k θ 1 (t -τ ) + c cos θ 0 q1 = m 21 α 1 + m 22 q 1 + (c m 0 cos θ 0 -c c 1 sin θ 0 ) cos θ 1 - -(c m 0 sin θ 0 + c c 1 cos θ 0 ) sin θ 1 + b 2 k θ 1 (t -τ )- -(c m 0 cos θ 0 -c c 1 sin θ 0 ) θ1 = q
The characteristic equation for the trivial equilibrium of (3. Recall from [START_REF] Hale | Theory of Functional Differential Equations[END_REF] the following theorem Theorem 3.1 ([8], Ch11, §11.1) Suppose that there exists τ c > 0 such that the characteristic equation (3.4), Δ(λ, τ c ) = 0, has a pair of simple purely imaginary roots λ 0 = iν 0 , λ0 = -iν 0 and all other roots have negative real parts. Suppose also that (3.6) Re λ 0 (τ c ) = 0.

Then a Hopf bifurcation occurs for τ = τ c , that is, for τ close to τ c every system (3.3) (thus system (2.3) too) has a periodic solution.

In order to find ν 0 and τ c we follow [START_REF] Cooke | On zeros of some transcendental equations[END_REF] and begin with the equation (P and Q are given in (3.5)). We suppose that conditions of Theorem 1 in [START_REF] Cooke | On zeros of some transcendental equations[END_REF] are fulfilled so, if ν 0 > 0 is a simple root of (3.7) it is a simple root of

|P (iy)| = |Q(iy)|.
Introduce the real and imaginary parts of P (iy) and Q(iy) through

P = P 1 + iP 2 , Q = Q 1 + iQ 2 so, P 1 (y) = Re P (iy), P 2 (y) = Im P (iy), Q 1 (y) = Re Q(iy), Q 2 (y) = Im Q(iy).
Remark that the solution of (3.1) is present in the coefficients of P .

Being a solution of (3.4) is equivalent for λ = iν 0 to

(3.8) cos τ ν 0 = - P 1 (ν 0 ) Q 1 (ν 0 ) + P 2 (ν 0 )Q 2 (ν 0 ) |Q(iν 0 )| 2 sin τ ν 0 = - P 1 (ν 0 ) Q 2 (ν 0 ) -P 2 (ν 0 ) Q 1 (ν 0 ) |Q(iν 0 )| 2
(by the first hypothesis in Theorem 1 in [START_REF] Cooke | On zeros of some transcendental equations[END_REF] and by (3.7), Q(iν 0 ) = 0). If τ c is a solution of (3.8), λ(τ c ) crosses the imaginary axis at λ = iν 0 in the direction given by (3.9) s = sgn Re dλ dτ

⏐ ⏐ ⏐ λ=iν0 = sgn F (ν 0 ).
As remarked in [START_REF] Cooke | On zeros of some transcendental equations[END_REF] the condition that ν 0 is a simple root of F (y) = 0 assures that, for τ = τ c , iν 0 is a simple root of (3.4).

A particular case of the formulae given in [START_REF] Cooke | On zeros of some transcendental equations[END_REF] implies that if (3.7) has only one positive root, this one crosses the imaginary axis, as τ increases to τ c , from left to right. So, if for τ = 0 the roots of the characteristic equation (3.4) are in the left half plane they will remain there until τ = τ c , thus equilibria given by (3.1), (α 0 , 0, θ 0 ), are stable for τ < τ c . The stability properties of the limit cycle are determined using the first Lyapunov coefficient. §4. Stability

of limit cycles

The basic reference for the formulae to compute the first Lyapunov coefficient is [START_REF] Hassard | Theory and Applications of Hopf Bifurcations[END_REF]. Similar computations were carried on in [START_REF] Adimy | Stability of limit cycles in a pluripotent stem cell dynamics model[END_REF].

With

(4.1) μ = τ -τ c equation (3.
3) can be rewritten as

(4.2) ẋ = G μ (x t ), t ≥ 0
where x t (s) = x(t + s) and G μ is a functional that acts on the space

C := C [-μ -τ c , 0] , C 3 of continuous functions defined on [-μ -τ c , 0] with values in C 3 . G μ is defined through (4.3) G (1) μ ⎛ ⎝ ϕ 1 ϕ 2 ϕ 3 ⎞ ⎠ = m 11 ϕ 1 (0) + m 12 ϕ 2 (0) + c cos θ 0 cos ϕ 3 (0)- -c sin θ 0 sin ϕ 3 (0) + b 1 k ϕ 3 (-μ -τ c ) + a 0 G (2) μ ⎛ ⎝ ϕ 1 ϕ 2 ϕ 3 ⎞ ⎠ = m 21 ϕ 1 (0) + m 22 ϕ 2 (0) + d 1 cos ϕ 3 (0)- -d 2 sin ϕ 3 (0) + b 2 k ϕ 3 (-μ -τ c ) + b 0 G (3) μ ⎛ ⎝ ϕ 1 ϕ 2 ϕ 3 ⎞ ⎠ = ϕ 2 (0)
where (4.4)

d 1 = c m 0 cos θ 0 -c c 1 sin θ 0 , d 2 = c m 0 sin θ 0 + c c 1 cos θ 0 .
The linearized equation is given by the Fréchet derivative of G μ in zero.

L μ = G μ (0) is explicitly given by (4.5) L μ ⎛ ⎝ ϕ 1 ϕ 2 ϕ 3 ⎞ ⎠ = = ⎛ ⎝ m 11 ϕ 1 (0) + m 12 ϕ 2 (0) -(c sin θ 0 ) ϕ 3 (0) + b 1 k ϕ 3 (-μ -τ c ) m 21 ϕ 1 (0) + m 22 ϕ 2 (0) -d 2 ϕ 3 (0) + b 2 k ϕ 3 (-μ -τ c ) ϕ 2 (0) ⎞ ⎠ . Introducing (4.6) F μ := G μ -L μ equation (4.2) becomes (4.7) ẋ = L μ x t + F μ (x t ), t ≥ 0. Obviously F μ (0) = 0, F μ (0) = 0. Define as in [10], [1], X 0 (θ) = 0, -μ -τ c ≤ θ < 0 1, θ = 0
, and then define (4.8)

A μ ϕ = ϕ + X 0 (L μ ϕ -ϕ (0)) , ϕ ∈ C.
Equation (4.7) can be written as

(4.9) ẋt = A μ x t + X 0 F μ (x t ). The initial data is x(s) = ϕ(s), s ∈ [-μ -τ c , 0]
. More details on operator A μ and initial problem (4.9) can be found in [START_REF] Hale | Theory of Functional Differential Equations[END_REF], [START_REF] Hassard | Theory and Applications of Hopf Bifurcations[END_REF], [START_REF] Faria | Normal forms for retarded functional differential equations and applications to Bogadanov-Takens singularity[END_REF].

An eigenvector of A 0 with respect to λ = iν 0 is given by γ

(θ) = e i ν0 θ ⎛ ⎝ γ 1 γ 2 γ 3 ⎞ ⎠ where (4.10) γ 1 = c sin θ 0 + b 1 k e -i ν0 τc -i ν 0 m 12 m 11 -i ν 0 , γ 2 = i ν 0 , γ 3 = 1.
For ϕ ∈ C and ψ ∈ C [0, τ c + μ] , C 3 define, according to [START_REF] Hale | Theory of Functional Differential Equations[END_REF], [START_REF] Hassard | Theory and Applications of Hopf Bifurcations[END_REF], [START_REF] Faria | Normal forms for retarded functional differential equations and applications to Bogadanov-Takens singularity[END_REF], the bilinear form

(4.11) ψ, ϕ = 3 j=1 ψ j (0) ϕ j (0) - 3 j=1 -μ-τc 0 ψ j (ξ + τ c + μ) ϕ j (ξ)dξ.
With respect to (4.11) the adjoint of A μ given in (4.8) is

A * μ ⎛ ⎝ ψ 1 ψ 2 ψ 3 ⎞ ⎠ = - ⎛ ⎝ ψ 1 ψ 2 ψ 3 ⎞ ⎠ + +X * 0 ⎡ ⎣ ⎛ ⎝ m 11 ψ 1 (0) + m 21 ψ 2 (0) m 12 ψ 1 (0) + m 22 ψ 2 (0) + ψ 3 (0) -(c sin θ 0 )ψ 1 (0) + b 1 k ψ 1 (μ + τ c ) -d 1 ψ 2 (0) -b 2 k ψ 2 (μ + τ c ) ⎞ ⎠ - - ⎛ ⎝ ψ 1 (0) ψ 2 (0) ψ 3 (0) ⎞ ⎠ ⎤ ⎦ . It follows that γ * (θ) = e i ν0 θ d ⎛ ⎝ γ1 γ2 γ3 ⎞ ⎠ is an eigenvector of A * 0 for (4.12) γ1 = 1, γ2 = i ν 0 -m 11 m 12 , γ3 = -m 12 + (i ν 0 -m 22 ) (i ν 0 -m 11 ) m 21
and for every d ∈ C. We choose d such that the normalizing condition γ * , γ = 1 be satisfied, thus d • (

3 j=1 γ j γj ) 1 + τ c e -i ν0 τc = 1 so (4.13) d = ⎡ ⎣ ( 3 j=1 γ j γj ) 1 + τ c e -i ν0 τc ⎤ ⎦ -1
.

Remark that, from Δ (i ν 0 , τ c ) = 0 and from (4.11) it follows that

(4.14) γ * , γ = 0.
The vectors γ and γ * are used to define the restriction of equation (4.9) to the section C 0 of the center manifold that corresponds to ±i ν 0 and to μ = 0 (see [START_REF] Hassard | Theory and Applications of Hopf Bifurcations[END_REF], [START_REF] Adimy | Stability of limit cycles in a pluripotent stem cell dynamics model[END_REF], [START_REF] Carr | Applications of Center Manifold Theory[END_REF], [START_REF] Chafee | A bifurcation problem for a functional differential equation of finitely retarded type[END_REF]). If x t is a solution of (4.9) define, for t ≥ 0, (4.15)

z(t) = γ * , x t
z and z will be used as local coordinates on C 0 in the directions γ * and γ * respectively. (4.11) implies that

z(t) = d ⎛ ⎝ 3 j=1 γj x j (t) - 3 j=1 γj -τc 0 e -i ν0 (ξ+τc) • x j (t + ξ) dξ ⎞ ⎠ .
Define, for t ≥ 0 and s ∈ [-τ c , 0],

(4.16) w(t, s) = x t (s) -2Re [z(t) γ(s)] = ⎛ ⎝ w (1) (t, s) w (2) (t, s) w (3) (t, s) ⎞ ⎠ .
On C 0 we have w(t, s) = W [z(t), z(t), s] where

(4.17) W (z, z, s) = ⎛ ⎜ ⎝ w (1) 20 (s) w (2) 20 (s) w (3) 20 (s)
⎞ ⎟ ⎠ z 2 2 + ⎛ ⎜ ⎝ w (1) 11 (s) w (2) 11 (s) w (3) 11 (s)
⎞ ⎟ ⎠ z z + ⎛ ⎜ ⎝ w (1) 02 (s) w (2) 02 (s) w (3) 02 (s)
⎞ ⎟ ⎠ z2 2 + . . .
For real x t , w is real so w 02 = w20 . From (4.14) and (4.15) it follows that γ * , w = γ * , x tγ * , γ * , x t γ + γ * , γ * , x t γ = 0. On the section C 0 , (4.18)

ż(t) = γ * , A 0 x t + X 0 [F 0 (x t )] = i ν 0 z(t) + γ * , F 0 [x(t)] (4.16) = = i ν 0 z(t) + γ * , F 0 [w(t, 0) + 2Rez(t) γ(0)] = = i ν 0 z(t) + γ * (0) T • f 0 (z(t), z(t)) := := i ν 0 z(t) + g [z(t), z(t)]
where dot means the scalar product in C 3 , T means transposed and f 0 (z, z) is given by the Taylor expansion of (3.3) around zero. Namely

f 1 (α 1 , q 1 , θ 1 ) = m 11 α 1 + m 12 q 1 + (c cos θ 0 ) • 1 - θ 2 1 2! + θ 4 1 4! -. . . - -(c sin θ 0 ) • θ 1 - θ 3 1 3! + . . . + a 0 + k b 1 θ 1 (t -τ c ) f 2 (α 1 , q 1 , θ 1 ) = m 21 α 1 + m 22 q 1 + d 1 • 1 - θ 2 1 2! + θ 4 1 4! -. . . - -d 2 • θ 1 - θ 3 1 3! + . . . + b 0 + k b 2 θ 1 (t -τ c ) (4.18
) is called the normal form obtained by restricting the flow to the center manifold. The Liapunov coefficients are computed using the coefficients of the Taylor expansion of g around (0, 0). Define (4.21)

L 1 (0) = i 2ν 0 g 20 g 11 -2|g 11 | 2 - 1 3 |g 02 | 2 + 1 2 g 21
The first Lyapunov coefficient is

(4.22) l 1 (0) = ReL 1 (0)
So, for l 1 (0), we need g 20 , g 11 and g 21 . It follows straightly from (4.20) that

d dt w(t, •) = d dt x t (•) -d dt [z(t)γ(•) + z(t)γ(•)] = = A 0 w(t, •) + X 0 F 0 [w(t, •) + 2Re (z(t)γ(•))] - -2Re [g(z(t), z(t))γ(•)] = = A 0 w(t, •) + H [z(t), z(t), •]
where, for s ∈ [-τ c , 0] we define

(4.26) H(z, z, s) = -2Re [g(z, z)γ(s)] + X 0 (s)F 0 [W (z, z, s) + 2Re [zγ(s)]] Thus, for s ∈ [-τ c , 0) H(z, z, s) = -2Re [g(z, z)γ(s)] = H 20 (s) z 2 2 + H 11 (s)z z + H 02 (s) z2 2 + . . . with (4.27) H 20 (s) = -g 20 γ(s) -ḡ02 γ(s) = -2Re [g 20 γ(s)] H 11 (s) = -2Re [g 11 γ(s)] H 02 (s) = H 20 (s)
. By (4.25) and (4.17)

(4.28) A 0 w(t, s) + H [z(t), z(t), s] = d dt w(t, s) = = d dt w 20 (s) z 2 (t) 2 + w 11 (s)z(t)z(t) + w 02 (s) z2 (t) 2 + . . . = = w 20 (s)z(t) ż(t) + w 11 (s) [ ż(t)z(t) + z(t) ż(t)] + w 02 (s)z(t) ż(t) + . . . = = w 20 (s)z [i ν 0 z + g(z, z)] + w 11 (s) [z (i ν 0 z + g(z, z)) + + z -i ν 0 z + g(z, z) + w 02 (s)z -i ν 0 z + g(z, z) (A 0 is acting on variable s).
Identifying the terms corresponding to z 2 , z z, z2 in (4.28) and recalling (4.17) we get (4.29)

(A 0 -2i ν 0 ) w 20 (s) = -H 20 (s) A 0 w 11 (s) = -H 11 (s) (A 0 + 2i ν 0 ) w 02 (s) = -H 02 (s) . So, ẇ20 = 2i ν 0 w 20 + g 20 e i ν0 s ⎛ ⎝ γ 1 i ν 0 1 ⎞ ⎠ + ḡ02 e -i ν0 s ⎛ ⎝ γ1 -i ν 0 1 ⎞ ⎠ therefore (4.30) ẇ(3) 20 (s) = 2i ν 0 w (3)
20 (s) + 2Re g 20 e i ν0 s .

Making s = 0 in (4.25) and identifying again the coefficient of z 2 one gets (4.31)

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ i ν 0 w (1) 20 (0) = m 11 w (1) 20 (0) + m 12 w (2) 20 (0) -(c 1 sin θ 0 )w (3) 20 (0)+ +kb 1 w (3) 20 (-τ c ) + H (1) 20 (0) i ν 0 w (2) 20 (0) = m 21 w (1) 20 (0) + m 22 w (2) 20 (0) -d 1 w (3) 20 (0)+ +kb 2 w (3) 20 (-τ c ) + H (2) 20 (0) i ν 0 w (3) 20 (0) = w
g 20 iν 0 + ḡ20 3iν 0 e 2 i ν0 s - g 20 iν 0 e i ν0 s - ḡ20 3iν 0 e -i ν0 s Then (4.32) w (3) 20 (-τ c ) = w (3) 20 (0) + g 20 iν 0 + ḡ20 3iν 0 e -2 i ν0 τc - g 20 iν 0 e -i ν0 τc - ḡ20 3iν 0 e i ν0 τc
and this value is to be introduced in (4.31).

NODY9605_source.tex; 30/10/2009; 8:33 p. 11

Identification in (4.28) of coefficients of z z when s = 0 leads to (4.33)

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ m 11 w (1)
11 (0) + m 12 w

(2)

11 (0) -(c sin θ 0 )w (3) 11 (0)+ +kb 1 w (3) 11 (-τ c ) + H (1) 11 (0) = 0 m 21 w (1) 11 (0) + m 22 w (2) 11 (0) -d 1 w (3) 11 (0)+ +kb 2 w (3) 11 (-τ c ) + H (2) 11 (0) = 0 w (2) 11 (0) + H (3) 11 (0) = 0 .
We infer from the second equation in (4.29) and from (4.27) that

ẇ11 (s) = g 11 γ(s) + ḡ11 γ(s) = g 11 e i ν0 s ⎛ ⎝ γ 1 γ 2 γ 3 ⎞ ⎠ + ḡ11 e -i ν0 s ⎛ ⎝ γ1 γ2 γ3 ⎞ ⎠ so w 11 (s) = g 11 e i ν0 s -1 i ν 0 ⎛ ⎝ γ 1 γ 2 γ 3 ⎞ ⎠ -ḡ11 e -i ν0 s -1 i ν 0 ⎛ ⎝ γ1 γ2 γ3 ⎞ ⎠ + w 11 (0).
Then, using (

(-τ c ) = g 11 e -i ν0 τc -1 i ν 0 -ḡ11 e i ν0 τc -1 i ν 0 + w 4.10), it follows that (4.34) w (3) 11
and this value will be introduced in (4.33).

In order to solve system (4.31) to find w

(3) 20 (0) and system (4.33) to find w From definition (4.26), using (4.19), one gets

H(z, z, 0) = -2Re [g(z, z) γ(0)] + F 0 [W (z, z, 0) + 2Re [z γ(0)]] = = -2Re ⎡ ⎣ ⎛ ⎝ γ 1 γ 2 γ 3 ⎞ ⎠ 1 2 g 20 z 2 + g 11 z z + 1 2 g 02 z2 + 1 2 g 21 z 2 z + . . . ⎤ ⎦ + + ⎛ ⎝ -c cos θ0 2 w (3) (0) + z + z 2 + c sin θ0 6 w (3) (0) + z + z 3 + . . . -d1 2 w (3) (0) + z + z 2 -d2 6 w (3) (0) + z + z 3 + . . . 0 ⎞ ⎠ With w (3) (0) = w (3) 20 (0) z 2 2 + w (3)
11 (0)z z + w

μ 2 = - l 1 (0) Re [λ (τ c )] , T 2 = - Im [L 1 (0)] + μ 2 Im [λ (τ c )] ν 0
we have Theorem 4.1 ([10]) If the Lyapunov coefficient l 1 (0) defined in (4.22) is negative, periodic solutions (limit cycles) exist, for equation (3.3), if τ > τ c , τ close to τ c , and are orbitally stable. They exist for τ < τ c and are unstable if l 1 (0) > 0. Their period increases if T 2 > 0 and decreases for T 2 < 0 (T 2 given in (4.37)).

By [START_REF] Hassard | Theory and Applications of Hopf Bifurcations[END_REF] the periodic solutions are approximated by In Figure 1, the dependence of τ c on k is depicted. The higher k is the smaller is the critical time for pilot reaction. The limit cycles are depicted in Figure 2.

(4.38) ϕ(t, μ) = 2 μ μ 2 1 2 Re ⎡ ⎣ e i ν0 t ⎛ ⎝ γ 1 γ 2 γ 3 ⎞ ⎠ ⎤ ⎦ + ⎛ ⎝ α 0 0 θ 0 ⎞ ⎠ (μ is defined in (4. 1
For k = 4, equilibria in (2.1) are α 0 = 0, 0954029; θ 0 = -0.00789153. The solution of (3.8) is τ c = 0.552144 and ν 0 = 2.04975. It follows that Re λ 0 (τ c ) = 1.8379 so Theorem 3.1 applies.

Computations on the lines in §4 yield l 1 (0) = -1.63799 so the limit cycle provided by Theorem 3.1 is a stable one (see Figures 3a and3b). Also μ 2 = 0.891232 and T 2 = 2.07495. The qualitative analysis of the mathematical model for a longitudinal motion with flight control system revealed that high pilot gains reduce pilot's allowed time delay, that is, in highly demanding tasks the pilot's reaction must be quick enough. It follows from the analysis of the stability of equilibria that, after a certain threshold τ c , the asymptotic stability is lost and a limit cycle appears through a Hopf bifurcation. Its characteristics depend on pilot's gain k. Numerical calculations of the Lyapunov coefficient show that in some situations the limit cycle is stable (for example when k ∈ {3, 4, 5}) while in other cases is unstable (k ∈ {1, 2, 6, 7}).

Numerical simulations performed for k = 4 and τ close to τ c , τ > τ c reveal a more complex situation: for constant functions as initial data not very close to the limit cycle (see (4.38)) the solution approaches another limit cycle, encircling the one given by Hopf bifurcation theorem (see Fig. 4). This is much alike the toric bifurcation (see [START_REF] Hassard | Theory and Applications of Hopf Bifurcations[END_REF]) and deserves further study both theoretical and numerical.

12 c

 12 sin θ 0b 1 k e -λτ -m 21 λm 22 c m 0 sin θ 0 + c c 1 cos θ 0b 2 k e -λτ 0 , τ ) = P (λ) + Q(λ)e -λτ where P, Q ∈ R[X], deg P = 3 and deg Q = 1.

(3. 7 )

 7 F (y) := |P (iy)| 2 -|Q(iy)| 2 = 0

2 g

 2 c cos θ 0 ) -γ2 d d1 02 = g 20 g 11 = d γ1 (-c cos θ 0 ) -γ2 d d 1 and (4.24) g 21 = d γ1 (-c cos θ 0 ) w (3) 20 (0) + 2w

  remark that it follows from (4.9), (4.16) and (4.18) that(4.25)

e

  2i ν0 (s-θ) Re g 20 e i ν0 θ dθ = = w(3) 20 (0)e 2 i ν0 s + e 2 i ν0 s

  we need H 20 (0) and H 11 (0).

2 ,

 2 identification of coefficients for z 2 and z z gives

  31) and (4.33) can be solved now and (4.24) gives g 21 and calculation of the Lyapunov coefficient l 1 (0) is completed. Define (4.37)

  )). §5. Case study (M = 0,25; H = 500 m)The formulae developed in paragraphs 3 and 4 are applied for a specific model based on ADMIRE (AeroData Model In a Research Environment, see[START_REF] Forssel | The Aero-Data Model In a Research Environment (ADMIRE) for Flight Control Robustness Evolution[END_REF]). Specifically m 11 = -1.3892; m 12 = 1.6688; c = 0.1161; b 1 = -0.5209; m 21 = 4.2859; m 22 = -14.4124; c 1 = 0.7153; b 2 = -6.3859; m 0 = -5.2641

Figure 1 :Figure 3 :

 13 Figure 1: Dependence of τ c on k

Figure 4 :

 4 Figure 4: Complex situation: two limit cycles

NODY9605_source.tex; 30/10/2009; 8:33 p. 3

NODY9605_source.tex; 30/10/2009; 8:33 p. 5

NODY9605_source.tex; 30/10/2009; 8:33 p. 10

NODY9605_source.tex; 30/10/2009; 8:33 p.[START_REF] Faria | Normal forms for retarded functional differential equations and applications to Bogadanov-Takens singularity[END_REF]

NODY9605_source.tex; 30/10/2009; 8:33 p. 15 the-Loop Oscillations -analysis and test techniques for their prevention" phase II 2004-2007.

Acknowledgement

This research was partially supported by Grant CNCSIS ID-919/2007 and by the National Authority for Scientific Research of Romania under contract C143 and was performed in the framework of GARTEUR FM/AG-15 "Pilot-in-

NODY9605_source.tex; 30/10/2009; 8:33 p. [START_REF] Hale | Theory of Functional Differential Equations[END_REF] We infer from (4.6) that (4.19)

w (1) (t, 0) + (z + z) γ 1 w (2) (t, 0) + (z + z)