HIGHER-ORDER MELNIKOV FUNCTIONS FOR SINGLE-DOF MECHANICAL OSCILLATORS: THEORETICAL TREATMENT AND APPLICATIONS

HIGHER-ORDER MELNIKOV FUNCTIONS FOR SINGLE-DOF MECHANICAL OSCILLATORS: THEORETICAL TREATMENT AND APPLICATIONS

STEFANO LENCI AND GIUSEPPE REGA Received 30 October 2003

A Melnikov analysis of single-degree-of-freedom (DOF) oscillators is performed by tak- ing into account the first (classical) and higher-order Melnikov functions, by considering Poincar´e sections nonorthogonal to the flux, and by explicitly determining both the dis- tance between perturbed and unperturbed manifolds (“one-half ” Melnikov functions) and the distance between perturbed stable and unstable manifolds (“full” Melnikov func- tion). The analysis is developed in an abstract framework, and a recursive formula for computing the Melnikov functions is obtained. These results are then applied to various mechanical systems. Softening versus hardening stiffness and homoclinic versus hetero- clinic bifurcations are considered, and the influence of higher-order terms is investigated in depth. It is shown that the classical (first-order) Melnikov analysis is practically inac- curate at least for small and large excitation frequencies, in correspondence to degenerate homo/heteroclinic bifurcations, and in the case of generic periodic excitations.

1. Introduction

The exact- or reduced-order nonlinear dynamics of many mechanical models or infinite- dimensional systems can be described by single-degree-of-freedom (DOF) oscillators, that is, by second-order ordinary nonlinear differential equations. Examples of the first kind are the mathematical pendulum [23], the inverted pendulum with lateral barriers [10,11], rocking rigid blocks [9], and many others, whereas buckled beams [16], shal- low arches [21], and cables [4] are just some samples of the second type. These systems are usually conservative plus damping, excitations, and possibly other kinds of perturba- tions which can be considered small in the first approximation. They are described by the equation

x¨= f(x) +εg1(x, ˙x,t) +ε²g2(x, ˙x,t) +···, x∈R, (1.1) whereεis a dimensionless smallness parameter measuring the amplitude of perturba- tions. Furthermore, in practical applications one often deals with smooth (or at least

Copyright©2004 Hindawi Publishing Corporation Mathematical Problems in Engineering 2004:2 (2004) 145–168 2000 Mathematics Subject Classification: 37C29, 34C37 URL:http://dx.doi.org/10.1155/S1024123X04310045

piecewise smooth) systems, and is interested in periodic excitations. This yields the fol- lowing hypothesis on (1.1):

(H1) f(x) andg_i(x, ˙x,t),i≥2, are sufficiently smooth and bounded on bounded sets, andgi(x, ˙x,t) areT-periodic int.

In (1.1) the mechanical differences between various systems are taken into account by considering different nonlinear stiffnesses, that is, different restoring forces f(x). These have strong consequences in terms of the dynamical response. For example, for soften- ing versus hardening systems, we have left versus right bending of the nonlinear reso- nance curve, escape versus scattered chaotic attractor for large excitation amplitude, and so forth. In spite of these important distinctions, there are some common dynamical fea- tures which permit a unified approach to the analysis of various oscillators.

Among others, we consider the saddles embedded in the system dynamics, physically representing unstable equilibrium positions, and their stable and unstable invariant man- ifolds. Indeed, it is well known that they play a central role in the nonlinear behaviour, being at the heart of such complex phenomena as multistability, chaotic dynamics [5,26], fractal basin boundaries [14], safe basin erosion and escape from potential wells [13,22], and so on.

Since the aim of this paper is to investigate some properties of the invariant manifolds of the system (1.1), we make the second hypothesis:

(H2) forε=0, the system possesses an orbitx^u(t) backward-asymptotic to a hyperbolic saddle pointp0=[p0, 0]^T, that is, limt→−∞x^u(t)=p0, f(p0)=0, and f,x(p0)= λ²>0.

Here and in the following, f,x(·) means the derivative of f(x) with respect to its argu- mentx.

Remarks. (1) Equation (1.1) is equivalent to the first-order system

˙ x=z,

˙

z=f(x) +εg1(x,z,t) +ε²g2(x,z,t) +···, x

z

∈R², (1.2)

which is sometimes useful from both a theoretical and a numerical point of view.

(2) Hypothesis (H1) guarantees the existence of the solution on a bounded region of the phase planeR², where we restrict our analysis.

(3) Forε=0, the system is called “unperturbed.” It is Hamiltonian withH(x, ˙x)=

˙

x²/2 + f(x)dx−C, where the constantC is chosen in such a way thatH(p0, 0)=0.

Thus, the orbitx^u(t) can be computed by solving the first-order equation

˙ x= ±

2

C−

f(x)dx

(1.3) with boundary condition limt→−∞x^u(t)=p0.

(4) The functiont→q0(t)=[x^u(t), ˙x^u(t)]^Tparametrically describes in the phase space theunstablemanifoldW^u(p0) ofp0.

t x˙

x

Σt0

p₀ p_ε

S

n t q0 q^u_ε

T t0

W^u(p₀) Wu(p_ε)

Figure 1.1. Unperturbed and perturbed manifolds and relevant Poincar´e sections.

(5) Since the unperturbed system is autonomous,x^u(t+ ˆt) is a solution for all ˆt∈R. This property is used in the variational approach to homo/heteroclinic bifurcations [2].

tˆis a parameter which permits to choose freely the pointq0=q0(0)=[x^u(0), ˙x^u(0)]^T∈

t0∩W^u(p0) around which we develop our analysis (see Figure 1.1andSection 2for further details). Note that ˆt=0 is usually assumed, but other choices can be made [9].

(6)x^u(t) solves ¨x^u(t)= f(x^u(t)). By taking the derivative of both sides with respect tot, we get...

x^u(t)= f,x(x^u(t)) ˙x^u(t), namely, ˙x^u(t) is one solution of the homogeneous variational equation ¨y=f_,x(x^u(t))y.

(7) Sincep0is a hyperbolic fixed point, ˙x^u(t) vanishes exponentially fort→ −∞, that is, ˙x^u(t)^∼=e^λt, whereλis defined in hypothesis (H2).

(8) The analysis is initially developed by referring to theunstablemanifoldW^u(p0).

The study of thestablemanifoldW^s(p0) is analogous and is summarized at the end of Section 2. Further results on stable and unstable manifolds of the unperturbed equa- tion (1.1) can be found in [8]. Combining results for stable and unstable manifolds per- mits detecting perturbed homo or heteroclinic distance (see the end ofSection 2) in both smooth [13] and nonsmooth systems [9].

(9) The extension to multi-DOF (and even infinite-dimensional [6]) dynamical sys- tems requires more sophisticated mathematical tools [25], but in some cases this entails only technicalities.

The object of this work is detecting the effects of the perturbationsεg1(x, ˙x,t) +ε²g2(x,

˙

x,t) +··· on the invariant manifolds. This is usually done by the classical Melnikov method [15], which is a perturbative technique permitting us to calculate the first-order distance (inε) between stable and unstable manifolds (see remark (8)). This method is specially conceived for determining homo/heteroclinic bifurcation thresholds, an issue which is very important from both a theoretical and a practical point of view because it permits the enlightenment of the associated dynamical phenomena and the skillful pursuit of their elimination or, possibly, their enhancement. This question has been sys- tematically investigated in the recent past (see [9,10,11,13,18]).

The present analysis mimics the spirit of the classical Melnikov’s method, which is illustrated, for example, in [5, Section 4.5] and [26, Section 4.5], two works which are our basic references. The perturbative approach is similar to that used by Vakakis [24], some results of which are also extended. More recent interesting results on the Melnikov function can be found in [3], while other works related to the present analysis are, for example, [7,19,20].

Three specific points, which seem not to have been investigated in classical analyses [5, 26], are studied in this paper. (i) We consider Poincar´e sections which arenotorthogonal to the unperturbed unstable manifoldW^u(p0) (resp.,W^s(p0)). (ii) We also consider the distance between perturbed and unperturbed unstable (stable) manifolds (related to the so-called “one-half ” Melnikov function). (iii) Moreover, we do not restrict to a first-order perturbative analysis but consider higher-order terms.

The first two extensions are motivated by the necessity of dealing with piecewise smooth systems (nonsmooth in the following for conciseness), where the Poincar´e sec- tion can be chosen on the discontinuity manifold, thus circumventing the difficulties of nonsmoothness. In this respect, it is worth stressing that the smoothness required in hy- pothesis (H1) must hold up to the discontinuity manifold. Among various applications, this permits the detection of heteroclinic bifurcations in the nonlinear dynamics of a generic rigid block [9], an issue which represents one of the motivations for this work.

The last point, on the other hand, is suggested by the fact that Melnikov method pro- vides a reliable approximation of the true distance for small ε, but when ε increases, or when the homo/heteroclinic bifurcation is degenerate (i.e., two distinct, simultane- ous points of tangency occur, seeSection 6), the error becomes appreciable. An example, which actually represents another motivation of this work, is reported at the end of [13, Section 4.3]. This paper is thus also aimed at overcoming this drawback by considering further terms in theε-series of the distance.

Some results of the present work were anticipated in [12].

2. The perturbative analysis

Before proceeding with the detection of the perturbed unstable manifolds, we need few preliminary comments. After having chosen one representation of the backward-asymp- totic orbitx^u(t), that is, after having chosen ˆt, see remark (5), we consider

x^u t−t0

, (2.1)

wheret0is the time where we fix the stroboscopic Poincar´e section_t0(seeFigure 1.1).

Another Poincar´e section is required. It is denoted bySand it is the plane passing throughq0and orthogonal to the vectorn=[cos(α), sin(α)]^T (in the phase space (x, ˙x)), as shown inFigure 1.1. In the classical Melnikov analysis,Sis orthogonal to the vector t=[ ˙x^u(0),f(x^u(0))]^T tangent to the unperturbed flux atq0(see (1.2) and remarks (4) and (5)), but in this work this restriction is removed.Sis only constrained to betransverse to the flux, that is,

n·t=cos(α) ˙x^u(0) + sin(α)f x^u(0)=0. (2.2)

Under the action of the perturbationsεg1(x, ˙x,t) +ε²g2(x, ˙x,t) +···, the hyperbolic fixed pointp0becomes a hyperbolic periodic orbitp_ε(t)ε-close top0, which on_t0cor- responds to the hyperbolic saddle pointpε, as shown inFigure 1.1[5, Lemma 4.5.1]. Fur- thermore, Lemma 4.5.2 of [5] guarantees that the unstable manifold ofp0converts to the unstable manifold ofp_ε, whose intersection with_t0is also shown inFigure 1.1, where it is denoted byW^u(pε). It intersectsSin the pointq_ε^u=W^u(pε)∩S∩

t0(Figure 1.1).

Thanks to [5, Lemma 4.5.2], we have that on the time interval−∞< t≤t0, the solution of (1.1) ending atq^u_ε can be expressed in the perturbative form

q^u_ε(t)=x^u t−t0

+εx^u₁(t) +ε²x^u₂(t) +···. (2.3) To determine the correction termsx^u_i(t), we insert (2.3) in (1.1), expand inε-series, and get (in the followingx^umeansx^u(t−t0) andx_i^ustands forx^u_i(t))

x¨^u−f x^u+εx¨₁^u−f_,x x^ux₁^u−g1 x^u, ˙x^u,t +ε²x¨^u₂−f,x x^ux₂^u−1

2f,xx x^u x₁^u2−g1,x x^u, ˙x^u,tx^u₁

−g1, ˙x x^u, ˙x^u,tx˙1^u−g2 x^u, ˙x^u,t+··· =0.

(2.4)

The first term is trivially zero by the definition ofx^u(t). Equating to zero the terms multi- plying successive powers ofε, we obtain the variational equations permitting the compu- tation ofx_i^u(t). They are the following linear second-order ordinary differential equations with nonconstant coefficients:

¨

x_i^u(t)=f_,x x^ux_i^u(t) +k^u_i(t), k₁^u(t)=g1 x^u, ˙x^u,t, k₂^u(t)=1

2f,xx x^u x₁^u2+g1,x x^u, ˙x^u,tx₁^u+g1, ˙x x^u, ˙x^u,tx˙^u₁+g2 x^u, ˙x^u,t,. . . . (2.5) First, we note that y^u1(t)=x˙^u(t−t0) is one solution of the homogenous version of (2.5), as shown in remark (6). The other can be computed by integrating the Wronskian expression

˙

y₁^u(t)y^u₂(t)−y˙^u₂(t)y₁^u(t)=1. (2.6) This is a first-order linear equation in the unknowny₂^u(t), whose solution, computed by the variation of constants method, can be written, at least formally, in the form

y₂^u(t)= −y^u₁(t) dτ

y₁^u(τ)². (2.7)

Remark 2.1. Since y^u₁(t)^∼₌e^λt,λ >0, for t→ −∞(see remark (7)), we have that y^u₂(t)

∼=e^−λt. Thus,y1^u(t) vanishes andy^u2(t) is unbounded fort→ −∞.

Remark 2.2. Note thaty^u1(t) andy^u2(t) donotdepend on the order of the approximation or on perturbations (the homogeneous equation is always the same) so that they must be computed one time only.

The knowledge of y^u₁(t) and y^u₂(t) permits to compute the general solution of (2.5), which is given by [27, equation (5.6.6), with changing the sign to negative] (note that equation (5.6.6) is correct up to the sign):

x_i^u(t)=y^u₁(t)

N_i^u+ _t

t0

y₂^u(τ)k_i^u(τ)dτ

+y₂^u(t)

M^u_i − _t

t0

y₁^u(τ)k_i^u(τ)dτ

. (2.8)

For the following purposes, it is also useful to compute its derivative, which is

˙

x_i^u(t)=y˙^u₁(t)

N_i^u+ _t

t0

y₂^u(τ)k_i^u(τ)dτ

+ ˙y₂^u(t)

M^u_i − _t

t0

y₁^u(τ)k_i^u(τ)dτ

. (2.9)

The constantsN_i^uandM_i^ucan be determined by boundary conditions. The first one is the boundedness ofx^u_i(t) fort→ −∞, which, becausey^u2(t) is unbounded fort→ −∞

(seeRemark 2.1), requires that the term multiplyingy₂^u(t) in (2.8) vanishes fort→ −∞, namely,

M_i^u=M_i^u t0

= _−∞

t0

y^u₁(τ)k^u_i(τ)dτ

= − _t0

−∞x˙^u τ−t0

k_i^u(τ)dτ

= − ₀

−∞x˙^u(τ)k_i^u τ+t0

dτ.

(2.10)

OnceM_i^uhave been computed, it remains to determineN_i^u. They follow from the require- ment that

q^u_ε =q^u_ε t=t0

=

x^u(0), ˙x^u(0)^T+εx₁^u t0

, ˙x₁^u t0

T

+ε²x₂^u t0

, ˙x₂^u t0

T

+··· (2.11)

belongs toS(Figure 1.1), namely, 0=

q^u_ε−q0

·n=εx₁^u t0

cos(α) + ˙x^u₁ t0 sin(α) +ε²x₂^u t0

cos(α) + ˙x^u₂ t0

sin(α)+···. (2.12)

Thus, the other boundary condition for (2.4) is x_i^u t0

cos(α) + ˙x^u_i t0

sin(α)=0 (2.13)

which, by taking into account (2.8) and (2.9) (fort=t0), yields y^u₁ t0

N_i^u+y₂^u t0

M_i^ucos(α) +y˙^u₁ t0

N_i^u+ ˙y₂^u t0

M_i^usin(α)=0, (2.14) N_i^u=N_i^u t0,α= −M_i^u t0

y^u₂ t0

cos(α) + ˙y₂^u t0

sin(α) y^u₁ t0

cos(α) + ˙y₁^u t0

sin(α). (2.15)

Note that the denominator of (2.15) is justn·t, which is different from zero due to the transversality condition (2.2). Three particular cases of (2.15) are worthy of attention: (i) nparallel tot, a case which is the classical assumption of the Melnikov theory [5,26] and leads toN_i^u= −M^u_i[y2^u(t0)y1^u(t0) + ˙y^u2(t0) ˙y1^u(t0)]/{[y1^u(t0)]²+ [ ˙y^u1(t0)]²}; (ii)nparallel to the ˙x-axis, that is,α=π/2, which yieldsN_i^u= −M_i^uy˙₂^u(t0)/y˙^u₁(t0); (iii)nparallel to the x-axis, that is,α=0, which givesN_i^u= −M^u_iy^u2(t0)/ y1^u(t0). Vakakis [24] notes thatN_i^u=0 can be viewed as a time shift, becausex^u(t−t0) +εN1^uy1^u(t) +··· =x^u(t−t0) +εN1^ux˙^u(t− t0) +··· ∼=x^u(t−t0+εN₁^u), and assumes thatN_i^u=0. The condition (iii) is used in [9].

Thanks to (2.12), the vectorq^u_ε−q0 has component only along the unit vectorm= [−sin(α), cos(α)]^Torthogonal ton. This component, which actually represents thesigned distance on_t0∩S between perturbed and unperturbed unstable manifolds, is finally given by

d^u t0,α= q^u_ε−q0

·m

=ε−x^u₁ t0

sin(α) + ˙x^u₁ t0

cos(α) +ε²−x^u₂ t0

sin(α) + ˙x^u₂ t0

cos(α)+···

=ε− y^u₁ t0

N₁^u+y₂^u t0

M₁^usin(α) +y˙^u₁ t0

N₁^u+ ˙y₂^u t0

M^u₁cos(α) +ε²−

y₁^u t0

N₂^u+y₂^u t0

M^u₂sin(α) +{y˙^u₁ t0

N₂^u+ ˙y₂^u t0

M^u₂cos(α)+···

= y^u₁ t0

˙ y^u₂ t0

−y˙^u₁ t0

y^u₂ t0

εM₁^u+ε²M₂^u+···

y^u₁ t0

cos(α) + ˙y₁^u t0

sin(α)

= −

εM1^u+ε²M^u2+···

x˙^u(0) cos(α) +fx^u(0)sin(α),

(2.16)

where use is made of (2.8) and (2.9) fort=t0and of (2.6). We explicitly report the ex- pressions of the first two coefficientsM_i^uwhich determine the distanced^u:

M₁^u t0

= − ₀

−∞x˙^u(τ)g1

x^u(τ), ˙x^u(τ),τ+t0

dτ, M₂^u t0

= −1 2

0

−∞x˙^u(τ)f_,xxx^u(τ)x₁^u τ+t0

2

dτ

− 0

−∞x˙^u(τ)g_1,xx^u(τ), ˙x^u(τ),τ+t0

x^u₁ τ+t0

dτ

− 0

−∞x˙^u(τ)g1, ˙x

x^u(τ), ˙x^u(τ),τ+t0

x˙^u₁ τ+t0

dτ

− 0

−∞x˙^u(τ)g2

x^u(τ), ˙x^u(τ),τ+t0 dτ.

(2.17)

The functionsM^u_i are the “one-half ” Melnikov functions of orderi.

Remark 2.3. The Melnikov functions have been obtained by setting equal to zero the coefficients of the unbounded part of (2.8). Generalizations of this idea in a functional analysis abstract framework are reported in [17].

The analysis can be repeated analogously for thestablemanifolds. In this case hy- pothesis (H2) requires the existence of an orbitx^s(t) forward-asymptotic to p0, that is, lim_t→+∞x^s(t)=p0, and it is justt→q0(t)=[x^s(t), ˙x^s(t)]^T that parametrically describes thestablemanifoldW^s(p0) ofp0. On the time intervalt0≤t <∞, the solution of equa- tion (1.1) starting fromq^s_ε=W^s(pε)∩S∩

t0can be expressed in the perturbative form (2.3), with a simple substitution of the apex “s” instead of “u”. The remaining part of the analysis is basically identical, up to this change of label. In particular, the correction termsx_i^s(t) are determined by solving “the same” variational problems (2.5). A slight dif- ference is that the vanishing ofy₁^s(t) and the unboundedness ofy^s₂(t) fort→+∞, instead oft→ −∞, are used in remark (7) andRemark 2.1. Thus, instead of (2.10), we have

M^s_i=M_i^s t0

= _∞

t0

y^s₁(τ)k^s_i(τ)dτ

= _∞

t0

˙ x^s τ−t0

k^s_i(τ)dτ= _∞

0 x˙^s(τ)k^s_i τ+t0

dτ,

(2.18)

while thesigneddistance on _t0∩Sbetween perturbed and unperturbed stable mani- folds is

d^s t0,α= q^s_ε−q0

·m= −

εM^s1+ε²M2^s+···

x˙^s(0) cos(α) +fx^s(0)sin(α). (2.19) The functionsM^s_i, which are conceptually and numerically distinct fromM_i^u, are the other “one-half ” Melnikov functions of orderi. The explicit expressions of the first two are

M₁^s t0

= _∞

0 x˙^s(τ)g1

x^s(τ), ˙x^s(τ),τ+t0

dτ, M₂^s t0

=1 2

_∞

0 x˙^s(τ)f_,xxx^s(τ)x₁^s τ+t0

2

dτ +

_∞

0 x˙^s(τ)g_1,xx^s(τ), ˙x^s(τ),τ+t0

x₁^s τ+t0

dτ +

_∞

0 x˙^s(τ)g1, ˙x

x^s(τ), ˙x^s(τ),τ+t0

x˙₁^s τ+t0

dτ +

_∞

0 x˙^s(τ)g2

x^s(τ), ˙x^s(τ),τ+t0

dτ.

(2.20)

With (2.16) and (2.20), it is then possible to compute the distance, onS, between stable and unstable perturbed manifolds:

d t0

=

q^s_ε−q^u_ε·m= q^s_ε−q0

− q^u_ε−q0

·m=d^s t0

−d^u t0

= −ε

M₁^s

˙

x^s(0) cos(α) +fx^s(0)sin(α)⁻

M₁^u

˙

x^u(0) cos(α) +fx^u(0)sin(α)

−ε²

M2^s

˙

x^s(0) cos(α) +fx^s(0)sin(α)⁻

M2^u

˙

x^u(0) cos(α) +fx^u(0)sin(α)

− ···. (2.21) Remark 2.4. The functions x^u(t) and x^s(t) are the restrictions of the same homo/

heteroclinic orbitx_h(t) to ]− ∞,t0] and [t0, +∞[, respectively. However, for nonsmooth systems [9],x^s(0)=x^u(0) and ˙x^s(0)=x˙^u(0) in general, and expression (2.21) cannot be further simplified. For smooth systems, on the other hand, x^s(0)=x^u(0)=xh(0) and

˙

x^s(0)=x˙^u(0)=x˙_h(0), and we have d t0

= −

εM1+ε²M2+···

˙

xh(0) cos(α) +fxh(0)sin(α), (2.22) whereMi(t0)=M_i^s(t0)−M_i^u(t0) are the “full” Melnikov functions of orderi. In particular, we have

M1 t0

= _∞

0 x˙^s(τ)g1

x^s(τ), ˙x^s(τ),τ+t0

dτ+ ₀

−∞x˙^u(τ)g1

x^u(τ), ˙x^u(τ),τ+t0

dτ

= _∞

−∞x˙h(τ)g1

xh(τ), ˙xh(τ),τ+t0 dτ,

(2.23)

which is the classical Melnikov function (see [5,26]).

Remark 2.5. Expression (2.21) applies to both homoclinic (the same unperturbed sad- dle for stable and unstable manifolds) and heteroclinic (different saddles for stable and unstable manifolds) orbits. In the case of heteroclinic orbits, however, the simple inter- section of perturbed stable and unstable manifolds is not sufficient to give “chaos” in the Smale’s horseshoe sense [26], and we must have a couple of intersecting manifolds (the so-called heteroclinic loop).

3. A different expression for correction terms

The differential equations (2.5), which are at the base of the previous perturbative anal- ysis, have the disadvantage that the known termsk^u_i(t) do not vanish fort→ −∞, and this may be unpleasant in numerical calculations of the integrals giving the closed-form expression of the solution. To overcome this point, the following trick is suggested.

First, we note that fort→ −∞, the perturbed solutionq^u_ε(t) must approach the saddle pε(t), which can also be determined by a perturbative analysis,

pε(t)=p0+ε˜x1(t) +ε²x˜2(t) +···, (3.1)

where, similarly to what has been done forq_ε(t), the terms ˜x_i(t) are the periodic solutions of (remember thatλ²=f_,x(p0), see remark (7))

x¨˜i(t)=λ²x˜i(t) + ˜ki(t), k˜1(t)=g1 p0, 0,t, k˜2(t)=1

2f_,xx p0

x˜²₁(t) +g_1,x p0, 0,tx˜1(t) +g1, ˙x p0, 0,tx˙˜1(t) +g2 p0, 0,t,. . . . (3.2)

On the basis of the previous considerations, we can write the corrections termsx^u_i(t) ofq^u_ε(t) in the form

x^u_i(t)=x˜i(t) + ˆx^u_i(t), (3.3) where ˜xi(t) is the asymptotically oscillating part ofx^u_i(t), and ˆx_i^u(t) is the decaying part.

By comparing (2.5) with (3.2), it is easy to see that ˆx^u_i(t) are solutions of x¨ˆ_i^u(t)=f,x x^uxˆ^u_i(t) + ˆk_i^u(t),

kˆ^u_i(t)=

f,x x^u−λ²x˜i(t) +k^u_i(t)−k˜i(t). (3.4) Equations (3.4) are similar to (2.5) so that the same expressions for the solutions de- veloped in the previous section apply, but now the known terms ˆk^u_i(t), contrary to the k_i^u(t), vanish fort→ −∞, and are therefore more easily handled. The same reasoning applies to the stable manifold.

4. Preliminary examples

To illustrate the features of the proposed perturbative approach for computing perturbed manifolds, in this section we consider two preliminary examples where computations can be done analytically. Two more interesting cases, however requiring numerical computa- tions of the involved integrals, will be discussed in the following sections.

We initially investigate the most simple case, that is, the linear equation

x¨+εδx˙−λ²x=εγsin(ωt), (4.1) which corresponds to f(x)=λ²x,g1(x, ˙x,t)= −δx˙+γsin(ωt), andgi(x, ˙x,t)=0,i≥2.

In this casex^u(t)=e^λt,y1(t)=λe^λ(t−t0),y2(t)=e^{−λ(t−t0)}/(2λ²). After some computations we get

x^u₁(t)= −δ 2 t−t0

e^λ(t−t0)−γsin(ωt)−sin ωt0

e^λ(t−t0)

λ²+ω² , (4.2)

which, in the special caset0=0 andα=0, provides q_ε(t)=x_h(t) +εx^u₁(t) +··· =e^λt+ε

−δ

2te^λt−γsin(ωt) λ²+ω²

+···. (4.3)

As expected, (4.3) contains the first two terms of theε-power series of the exact per- turbed solution

q_ε(t)=e^zt−εγ λ²+ω²sin(ωt) +εδωcos(ωt)−e^zt

(εδω)²+ λ²+ω²² , z=−εδ+(εδ)²+ 4λ²

2 .

(4.4) In this case the pointq0isq0=[x^u(0), ˙x^u(0)]^T=[1,λ]^T. By choosingα=0, that is, the sectionSperpendicular to thex-axis, the distance onSinq0is given by

d^u t0

= −εδ

2 ⁻εγωcos ωt0

−λsin ωt0 λ²+ω² +(εδ)²

8λ + (εδ)(εγ)8λωcos ωt0

+ 4 ω²−λ²sin ωt0

λ²+ω²² +···.

(4.5)

Note that the coupling between damping and excitation occurs only at second order.

The second preliminary example is the escape (Helmholtz) oscillator subjected to a damping which, though small, is larger than the excitation. To take this fact into account, we scale to different powers ofεthe two perturbation terms and consider

x¨+εδx˙−x+x²=ε²γsin(ωt). (4.6) In this case we have f(x)=x−x²,g1(x, ˙x,t)= −δx,˙ g2(x, ˙x,t)=γsin(ωt), andgi(x, ˙x,t)= 0,i≥3. Easy computations give the expression of thehomoclinicorbit and of the auxiliary functionsy1andy2:

xh=3 2

1

cosh²(t/2), y1= −3 2

sinh(t/2) cosh³(t/2), y2=5

6+cosh²(t/2)

3 ⁻

5

2 cosh²(t/2)+ 5tsinh(t/2) 4 cosh³(t/2).

(4.7)

If we look for the perturbations around the point described by ˆt=0, that is,q0=[3/2, 0]^T, and measure the distance on the sectionSperpendicular to the y-axis, that is,α=π/2, we obtain

x₁^u(t)= δ

10e^t e^3t+ 15e^2t−65e^t−15+ 30t e^t−1

e^t+ 1³ ,

x₁^s(t)= δ

10e^{−t −}e^−3t−15e^−2t+ 65e^−t+ 15+ 30t e^−t−1

e^−t+ 1³ .

(4.8)

The functionsx_h(t) (unperturbed stable and unstable manifolds),x_h(t) +εx^u₁(t) (per- turbed unstable manifold), andx_h(t) +εx₁^s(t) (perturbed stable manifold) are depicted in

4 2

−2 0

−4

0.2 0.4 0.6 0.8 1 1.2 1.4

x

t (a)

0 0.2 0.4 0.6

−0.2

−0.4

−0.6

0.2 0.4 0.6 0.8 1 1.2 1.4 x

x˙

(b)

Figure 4.1. The functionsxh(t) (thin line),xh(t) +εx1^u(t) (thick continuous line), andxh(t) +εx1^s(t) (thick dashed line) forεδ=0.05: (a) time histories and (b) phase space.

Figure 4.1forεδ=0.05. The distances, on the other hand, are given by d^u t0

=4 5εδ+

49 75⁻

8 5ln 2

(εδ)² +ε²γ

cos ωt0

4πω² sinh(πω) +2 sin ωt0

−1+ω²

ψiω 2 +ψ⁻iω

2 ⁻ψ 1

2+iω 2

−ψ 1

2⁻ iω

2

+···, d^s t0

= −4 5εδ+

49 75⁻

8 5ln 2

(εδ)² +ε²γ

−cos ωt0

4πω² sinh(πω) +2 sin ωt0

−1+ω²

ψiω 2 +ψ⁻iω

2 ⁻ψ 1

2+iω 2

−ψ 1

2⁻ iω

2

+···, d t0

= −8

5εδ−ε²γcos ωt0 8πω²

sinh(πω)+···,

(4.9) whereψ(·) is the digamma function [1]. Owing to the considered powers ofε, there is no coupling at second order between damping and excitation.

5. The mathematical pendulum

We consider the damped mathematical pendulum subjected to horizontal harmonic ex- citation:

¨

x+εδx˙+ sin(x) +εγcos(ωt) cos(x)=0. (5.1)