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

A Melnikov analysis of single-degree-of-freedom (DOF) oscillators is performed by taking into account the first (classical) and higher-order Melnikov functions, by considering Poincaré sections nonorthogonal to the flux, and by explicitly determining both the distance between perturbed and unperturbed manifolds (“one-half” Melnikov functions) and the distance between perturbed stable and unstable manifolds (“full” Melnikov function). 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 heteroclinic 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 inaccurate at least for small and large excitation frequencies, in correspondence to degenerate homo/heteroclinic bifurcations, and in the case of generic periodic excitations.


Introduction
The exact-or reduced-order nonlinear dynamics of many mechanical models or infinitedimensional 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], shallow 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 perturbations which can be considered small in the first approximation.They are described by the equation where ε is a dimensionless smallness parameter measuring the amplitude of perturbations.Furthermore, in practical applications one often deals with smooth (or at least piecewise smooth) systems, and is interested in periodic excitations.This yields the following hypothesis on (1.1): (H1) f (x) and g i (x, ẋ,t), i ≥ 2, are sufficiently smooth and bounded on bounded sets, and g i (x, ẋ,t) are T-periodic in t.
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 softening versus hardening systems, we have left versus right bending of the nonlinear resonance curve, escape versus scattered chaotic attractor for large excitation amplitude, and so forth.In spite of these important distinctions, there are some common dynamical features 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 manifolds.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 orbit x u (t) backward-asymptotic to a hyperbolic saddle point p 0 = [p 0 ,0] T , that is, lim t→−∞ x u (t) = p 0 , f (p 0 ) = 0, and f ,x (p 0 ) = λ 2 > 0.
Here and in the following, f ,x (•) means the derivative of f (x) with respect to its argument x.
Remarks.(1) Equation (1.1) is equivalent to the first-order system ẋ = z, 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 plane R 2 , where we restrict our analysis.
(4) The function t → q 0 (t) = [x u (t), ẋu (t)] T parametrically describes in the phase space the unstable manifold W u (p 0 ) of p 0 .(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 point q 0 = q 0 (0) = [x u (0), ẋu (0)] T ∈ t0 ∩W u (p 0 ) around which we develop our analysis (see Figure 1.1 and Section 2 for further details).Note that t = 0 is usually assumed, but other choices can be made [9].(6) x u (t) solves ẍu (t) = f (x u (t)).By taking the derivative of both sides with respect to t, we get ... x u (t) = f ,x (x u (t)) ẋu (t), namely, ẋu (t) is one solution of the homogeneous variational equation ÿ = f ,x (x u (t))y.
(7) Since p 0 is a hyperbolic fixed point, ẋu (t) vanishes exponentially for t → −∞, that is, ẋu (t) ∼ = e λt , where λ is defined in hypothesis (H2).(8) The analysis is initially developed by referring to the unstable manifold W u (p 0 ).The study of the stable manifold W s (p 0 ) is analogous and is summarized at the end of Section 2. Further results on stable and unstable manifolds of the unperturbed equation (1.1) can be found in [8].Combining results for stable and unstable manifolds permits detecting perturbed homo or heteroclinic distance (see the end of Section 2) in both smooth [13] and nonsmooth systems [9].(9) The extension to multi-DOF (and even infinite-dimensional [6]) dynamical systems 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 εg 1 (x, ẋ,t) + ε 2 g 2 (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 systematically 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é sections which are not orthogonal to the unperturbed unstable manifold W u (p 0 ) (resp., W s (p 0 )).(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é section 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 hypothesis (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 provides a reliable approximation of the true distance for small ε, but when ε increases, or when the homo/heteroclinic bifurcation is degenerate (i.e., two distinct, simultaneous points of tangency occur, see Section 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].

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-asymptotic orbit x u (t), that is, after having chosen t, see remark (5), we consider where t 0 is the time where we fix the stroboscopic Poincaré section t0 (see Figure 1.1).Another Poincaré section is required.It is denoted by S and it is the plane passing through q 0 and orthogonal to the vector n = [cos(α),sin(α)] T (in the phase space (x, ẋ)), as shown in Figure 1.1.In the classical Melnikov analysis, S is orthogonal to the vector t = [ ẋu (0), f (x u (0))] T tangent to the unperturbed flux at q 0 (see (1.2) and remarks (4) and ( 5)), but in this work this restriction is removed.S is only constrained to be transverse to the flux, that is, S. Lenci and G. Rega 149 Under the action of the perturbations εg 1 (x, ẋ,t) + ε 2 g 2 (x, ẋ,t) + •••, the hyperbolic fixed point p 0 becomes a hyperbolic periodic orbit p ε (t) ε-close to p 0 , which on t0 corresponds to the hyperbolic saddle point p ε , as shown in Figure 1.1 [5,Lemma 4.5.1].Furthermore, Lemma 4.5.2 of [5] guarantees that the unstable manifold of p 0 converts to the unstable manifold of p ε , whose intersection with t0 is also shown in Figure 1.1, where it is denoted by W u (p ε ).It intersects S in the point q u ε = W u (p ε ) ∩ S ∩ t0 (Figure 1.1).Thanks to [5, Lemma 4.5.2],we have that on the time interval −∞ < t ≤ t 0 , the solution of (1.1) ending at q u ε can be expressed in the perturbative form To determine the correction terms x u i (t), we insert (2.3) in (1.1), expand in ε-series, and get (in the following x u means x u (t − t 0 ) and x u i stands for x u i (t)) The first term is trivially zero by the definition of x u (t).Equating to zero the terms multiplying successive powers of ε, we obtain the variational equations permitting the computation of x u i (t).They are the following linear second-order ordinary differential equations with nonconstant coefficients: x u , ẋu ,t x u 1 + g 1, ẋ x u , ẋu ,t ẋu 1 + g 2 x u , ẋu ,t ,.... (2.5) First, we note that y u 1 (t) = ẋu (t − t 0 ) 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 This is a first-order linear equation in the unknown y u 2 (t), whose solution, computed by the variation of constants method, can be written, at least formally, in the form (2.7) Remark 2.1.Since y u 1 (t) ∼ = e λt , λ > 0, for t → −∞ (see remark ( 7)), we have that y u 2 (t) ∼ = e −λt .Thus, y u 1 (t) vanishes and y u 2 (t) is unbounded for t → −∞.Remark 2.2.Note that y u 1 (t) and y u 2 (t) do not depend 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.150 Higher-order Melnikov functions for SDOF oscillators The knowledge of y u 1 (t) and y u 2 (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): For the following purposes, it is also useful to compute its derivative, which is (2.9) The constants N u i and M u i can be determined by boundary conditions.The first one is the boundedness of x u i (t) for t → −∞, which, because y u 2 (t) is unbounded for t → −∞ (see Remark 2.1), requires that the term multiplying y u 2 (t) in (2.8) vanishes for t → −∞, namely, (2.10) Once M u i have been computed, it remains to determine N u i .They follow from the requirement that belongs to S (Figure 1.1), namely, (2.12) Thus, the other boundary condition for (2.4) is which, by taking into account (2.8) and (2.9) (for t = t 0 ), yields ) S. Lenci and G. Rega 151 Note that the denominator of (2.15) is just n • t, which is different from zero due to the transversality condition (2.2).Three particular cases of (2.15) are worthy of attention: (i) n parallel to t, a case which is the classical assumption of the Melnikov theory [5,26] and leads to 2 (t 0 )/y u 1 (t 0 ).Vakakis [24] notes that N u i = 0 can be viewed as a time shift, because x u (t , and assumes that N u i = 0.The condition (iii) is used in [9].Thanks to (2.12), the vector q u ε − q 0 has component only along the unit vector m = [−sin(α),cos(α)] T orthogonal to n.This component, which actually represents the signed distance on t0 ∩S between perturbed and unperturbed unstable manifolds, is finally given by where use is made of (2.8) and (2.9) for t = t 0 and of (2.6).We explicitly report the expressions of the first two coefficients M u i which determine the distance d u : (2.17) The functions M u i are the "one-half " Melnikov functions of order i.
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 the stable manifolds.In this case hypothesis (H2) requires the existence of an orbit x s (t) forward-asymptotic to p 0 , that is, lim t→+∞ x s (t) = p 0 , and it is just t → q 0 (t) = [x s (t), ẋs (t)] T that parametrically describes the stable manifold W s (p 0 ) of p 0 .On the time interval t 0 ≤ t < ∞, the solution of equation (1.1) starting from q s ε = W s (p ε ) ∩ S ∩ t0 can 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 terms x s i (t) are determined by solving "the same" variational problems (2.5).A slight difference is that the vanishing of y s 1 (t) and the unboundedness of y s 2 (t) for t → +∞, instead of t → −∞, are used in remark ( 7) and Remark 2.1.Thus, instead of (2.10), we have while the signed distance on t0 ∩S between perturbed and unperturbed stable manifolds is The functions M s i , which are conceptually and numerically distinct from M u i , are the other "one-half " Melnikov functions of order i.The explicit expressions of the first two are With (2.16) and (2.20), it is then possible to compute the distance, on S, between stable and unstable perturbed manifolds: S. Lenci and G. Rega 153 (2.21) Remark 2.4.The functions x u (t) and x s (t) are the restrictions of the same homo/ heteroclinic orbit x h (t) to ] − ∞, t 0 ] and [t 0 ,+∞[, respectively.However, for nonsmooth systems [9], x s (0) = x u (0) and ẋs (0) = ẋ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) = x h (0) and ẋs (0) = ẋu (0) = ẋh (0), and we have where M i (t 0 ) = M s i (t 0 ) − M u i (t 0 ) are the "full" Melnikov functions of order i.In particular, we have which is the classical Melnikov function (see [5,26]).
Remark 2.5.Expression (2.21) applies to both homoclinic (the same unperturbed saddle for stable and unstable manifolds) and heteroclinic (different saddles for stable and unstable manifolds) orbits.In the case of heteroclinic orbits, however, the simple intersection 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).

A different expression for correction terms
The differential equations (2.5), which are at the base of the previous perturbative analysis, have the disadvantage that the known terms k u i (t) do not vanish for t → −∞, 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.
On the basis of the previous considerations, we can write the corrections terms x u i (t) of q u ε (t) in the form where xi (t) is the asymptotically oscillating part of x u i (t), and xu i (t) is the decaying part.By comparing (2.5) with (3.2), it is easy to see that xu i (t) are solutions of (3.4) Equations (3.4) are similar to (2.5) so that the same expressions for the solutions developed in the previous section apply, but now the known terms ku i (t), contrary to the k u i (t), vanish for t → −∞, and are therefore more easily handled.The same reasoning applies to the stable manifold.

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 computations of the involved integrals, will be discussed in the following sections.
We initially investigate the most simple case, that is, the linear equation which corresponds to f (x) = λ 2 x, g 1 (x, ẋ,t) = −δ ẋ + γ sin(ωt), and In this case x u (t) = e λt , y 1 (t) = λe λ(t−t0) , y 2 (t) = e −λ(t−t0) /(2λ 2 ).After some computations we get which, in the special case t 0 = 0 and α = 0, provides S. Lenci and G. Rega 155 As expected, (4.3) contains the first two terms of the ε-power series of the exact perturbed solution (4.4) In this case the point q 0 is q 0 = [x u (0), ẋu (0)] T = [1,λ] T .By choosing α = 0, that is, the section S perpendicular to the x-axis, the distance on S in q 0 is given by 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 In this case we have f (x) = x − x 2 , g 1 (x, ẋ,t) = −δ ẋ, g 2 (x, ẋ,t) = γ sin(ωt), and g i (x, ẋ,t) = 0, i ≥ 3. Easy computations give the expression of the homoclinic orbit and of the auxiliary functions y 1 and y 2 : sinh(t/2) cosh 3 (t/2) , If we look for the perturbations around the point described by t = 0, that is, q 0 = [3/2,0] T , and measure the distance on the section S perpendicular to the y-axis, that is, α = π/2, we obtain   + ε 2 γ cos ωt 0 4πω 2 sinh(πω) where ψ(•) is the digamma function [1].Owing to the considered powers of ε, there is no coupling at second order between damping and excitation.

The mathematical pendulum
We consider the damped mathematical pendulum subjected to horizontal harmonic excitation: ẍ + εδ ẋ + sin(x) + εγ cos(ωt)cos(x) = 0.In this case we have f (x) = −sin(x), g 1 (x, ẋ,t) = −δ ẋ − γ cos(ωt)cos(x), and g i (x, ẋ,t) = 0, i ≥ 2. Simple computations give the expression of the heteroclinic orbit and of the auxiliary functions y 1 and y 2 : 2) The saddle is p 0 = −π, and we focus attention on the point q 0 = [0,2] T which is obtained for t = 0.The section S is perpendicular to the x-axis, that is, α = 0, and is also orthogonal to the tangent vector t.After some computations, we obtain where the expressions of functions a(s), b(s,ω), and c(s,ω) are reported in the appendix.
On the basis of these results we can compute the first two "one-half " Melnikov functions, which are given by respectively.The functions M 1a (ω), M 1b (ω), M 2a (ω), M 2b (ω), M 2c (ω), M 2d (ω), and M 2e (ω) are depicted in Figure 5.1, and their expressions are reported in the appendix.Figure 5.1 shows that for ω → ∞, all the functions M i j (ω) tend to zero so that M u 1 → 4δ and M u 2 → 2δ 2 .This means that the unstable perturbed and unperturbed manifolds keep well disjoint for large frequencies.Physically, this can be justified by noting that when the excitation oscillations are very fast, basically they do not affect the system, and only the manifolds separation due to damping survives.Note that the phenomenon is captured by the first-order term and then confirmed by the second-order term.A different situation is observed in the opposite case of quasistatic excitation, ω → 0. Here, in fact, we still have M u 1 → 4δ, but now M u 2 → 2δ 2 − 0.5γ 2 + 2δγ.Thus, while the excitation amplitude does not contribute to the first-order term, it appears in the second-order term, which therefore adds not only a quantitative but also a qualitative difference, and it is not negligible in this range of frequencies.In particular, if we take into account only the first-order term, we miss the property that the distance may vanish for sufficiently large excitation amplitudes.
The second-order term has an appreciable effect not only for small ω.To illustrate this by an example, we consider ω = 1.The associated distance between perturbed and unperturbed unstable manifolds is given by 4(εδ) + 1.252(εγ)cos t 0 + 0.297(εγ)sin t 0 + 2(εδ) 2 + 0.219(εγ) 2 cos 2 t 0 − 0.262(εγ) 2 sin 2 t 0 + 0.139(εγ) 2 sin t 0 cos t 0 + 0.522(εδ)(εγ)cos t 0 + 1.166(εδ)(εγ)sin t 0 . (5.5) The distance (5.5) is reported in Figure 5.2 for various values of εδ and εγ.This picture clearly shows the modifications due to the second-order terms.Note that the differences are not negligible even for the relatively small values of damping and excitation amplitude employed in the picture.The other "one-half " Melnikov functions M s 1 (t 0 ) and M s 2 (t 0 ) can be computed analogously.They share the same properties of M u 1 (t 0 ) and M u 2 (t 0 ), and are not reported.

The escape oscillator
Other aspects of the importance of higher-order Melnikov terms can be highlighted by considering the damped asymmetric escape (Helmholtz) oscillator subjected to a generic periodic excitation: We remark that, contrary to the Helmholtz oscillator analyzed in Section 4, here the damping and the excitation have the same order of magnitude ε.Another difference is that the excitation is no longer harmonic.The functions f (x), x h , y 1 , and y 2 are reported in Section 4, g 1 (x, ẋ,t) = −δ ẋ + γ(t), and g i (x, ẋ,t)=0, i≥2.The saddle is p 0 = 0, we focus attention on the point q 0 = [3/2,0] T , and we consider α = π/2, that is, n orthogonal to the x-axis and S perpendicular to the tangent vector t.When compared with that of Section 5, this case considers homoclinic instead of heteroclinic orbits, and we focus directly on the "full" Melnikov functions instead of dealing with "one half " as in Section 5.
After some algebra, the first two Melnikov functions can be expressed in the form β jk (ω)γ j γ k sin jωt 0 + ψ j sin kωt 0 + ψ k . (6. 2) The expression of M 1 (t 0 ) is a known result, see [13], while the functions χ(s), β 11 (ω), β 12 (ω), and β 21 (ω), whose expressions are reported in the appendix, are depicted in Figure 6.1.Note that β kk (ω) = β 11 (kω), and this is why β 22 (ω) is not reported in Figure 6.1(b).In (6.2) we do not have the term δ 2 according to the fact that d(t 0 ) must be an odd function in δ.In fact, in the absence of the excitation, changing the sign of δ means changing the direction of time (if x(t) is a solution for a given δ, then x(−t) is the solution for −δ).This implies that stable and unstable manifolds exchange, that is, d u (t 0 ) becomes d s (t 0 ) and vice versa.Accordingly, d(t 0 ) = d s (t 0 ) − d u (t 0 ) changes sign, namely, it is odd in δ.This of course does not require d s (t 0 ) and d u (t 0 ) to be odd in δ (indeed, see (5.4)), but simply that the coefficients of δ 2i in d u (t 0 ) and d s (t 0 ) are the same so that they disappear in d(t 0 ) = d s (t 0 ) − d u (t 0 ).Expressions (6.2) share the same asymptotic properties of (5.4).In fact, M 1 → −1.2δ and M 2 → 0 for ω → ∞, thus showing that stable and unstable perturbed manifolds keep well disjoint for large frequencies.Contrary to (5.4), however, the second-order terms do not affect the analysis.For ω → 0, on the other hand, we have that γ j sin ψ j .( Thus, apart from the special case ψ j = 0 for which M 2 → 0, we have that second-order terms entail qualitative corrections to the first order as the superharmonics amplitudes appear in the expression of the distance.The Melnikov functions are usually of interest because the vanishing of the distance d(t 0 ) for some t 0 implies transverse homoclinic intersection and, roughly, chaos (see [5,26]).The smallest excitation value for which this can occur is the homoclinic bifurcation threshold (actually, it corresponds to homoclinic tangency), whose detection has recently been investigated in the case of classical Melnikov function [13].
The curves εγ 1,cr (ω) obtained considering only first-order term and first-plus secondorder terms are reported in Figure 6.2.This picture clearly shows how the classical Melnikov analysis (thin lines) is accurate only for medium values of ω around the chaotic resonance ω = 0.6096, where the second-order correction basically does not modify the classical results.For low and especially for high excitation frequencies, the error becomes dramatically large, and the classical results largely overestimate the true εγ 1,cr , even for small values of the damping.As a matter of fact, it is believed that for very large ω, the analysis based on second-order terms is still inaccurate, and further higher-order Melnikov functions should be taken into account.
Contrary to what is suggested by Figure 6.2, the second-order Melnikov functions may have nonnegligible effects also for medium excitation frequencies.This is seen by considering generic periodic excitations.In this case we can illustrate by an example the differences between theoretical predictions based on first-order Melnikov function and numerical simulations, a discrepancy which can be eliminated by considering M 2 .We refer to the example of [13,Figure 15], which considers εδ = 0.1, ψ 1 = 0, ω = 0.7, and N = 2, and is therefore the simplest nonharmonic case.The explanation of the inconsistency of [13,Figure 15] actually constitutes one of the initial motivations of this work.
We initially consider εγ 2 = 0.8 × εγ 1 and ψ 2 = 0.This choice arises in the framework of optimal chaos control [13], but this is irrelevant to the purposes of this work.If we let εγ 1 increase from zero, computer simulations of (6.1) show the first homoclinic tangency to occur at εγ num 1,cr = 0.0794, as shown in Figure 6.3(a), where the first touching point A is emphasized.However, the first-order Melnikov function is always less than zero, as shown by the thin line in Figure 6.4(a), and thus we theoretically predict no tangency, contrary to the numerical evidence.If, however, we add the second-order Melnikov function (thick line of Figure 6.4(a)), we theoretically recover the homoclinic tangency, that is, εM 1 (t 0 ) + ε 2 M 2 (t 0 ) = 0 for a certain t 0 .
The theoretical analysis based on εM 1 (t 0 ) predicts homoclinic tangency for εγ 1 ord 1,cr = 0.0818, as shown by the thin line in Figure 6.4(b).In this case, however, we numerically observe homoclinic intersection (Figure 6.3(b)), with the point A beyond the unstable manifold and the point B in front of W u (p ε ), a circumstance which is theoretically retrieved if we consider also the second-order terms (thick line in Figure 6.4(b)).An important aspect highlighted by Figure 6.4(b) is that a degenerate homoclinic bifurcation is expected by considering only first-order analysis, as the tangency is predicted 162 Higher-order Melnikov functions for SDOF oscillators to occur simultaneously in two distinct primary intersection points (PIP, see [26]).This codimension-two event is highly sensitive to small parameters variations, and this shows the inaccuracy of the first-order analysis and the necessity of considering at least secondorder terms.
If we want to recover the degenerate homoclinic bifurcations in real manifolds (this has interest in the field of optimal chaos control [13]), we must correct the excitation εγ(t) = 0.0818 × [sin(ωt) + 0.8 × sin(2ωt)] of Figures 6.3(b) and 6.4(b) by a term of the form −µ cos(2ωt).When µ increases, the point A is seen to move back while the point B approaches the unstable manifold: at last, for µ = 0.0818 × 0.048, there is a simultaneous homoclinic tangency.This is seen numerically in Figure 6. with first-order theoretical analysis are confirmed.The singular character of the degenerate excitation of Figures 6.3(c) and 6.4(c) is demonstrated by the fact that, by further increasing µ, the point A detaches from the unstable manifold while B intersects it, as shown numerically in Figure 6.3(d) and theoretically in Figure 6.4(d).
We note that since the predictions based on εM 1 (t 0 ) are corrected up to the ε order, the difference between the theoretical and numerical "degenerate" excitations must be of ε 2 order.This agrees very well with the obtained results, since the amplitude of the second-order "correction" terms is at least one order of magnitude smaller than the amplitudes of the first-order terms.However, we have shown that in the neighborhood of a codimension-two homoclinic bifurcation, the ensuing corrections in the response amplitudes, although small, are no longer negligible, even in the range of medium frequencies where the first-order analysis is accurate with harmonic excitation.On the basis of these results, we conjecture that higher-order analysis is always required in the case of degenerate homo/heteroclinic bifurcations.
To investigate in more detail the effects on the homoclinic bifurcation threshold of the first superharmonic εγ 2 cos(2ωt + ψ 2 ) added to the principal harmonic εγ 1 cos(ωt + ψ 1 ), we have reported in Figure 6.5(a) the critical threshold in the plane (εγ 1 ,εγ 2 ) in the case εδ = 0.1, ω = 0.7 (the same values of Figures 6.3 and 6.4), and ψ 1 = ψ 2 = 0.The points inside the curves correspond to no homoclinic intersection, which instead occurs for the excitations corresponding to outer points.The line εγ 2 = 0.8 × εγ 1 along which we developed the analysis of the previous point is also represented by a dashed line.
We see that first-and second-order analyses (represented by thin and thick lines, resp.)basically coincide along the lines εγ 2 = 0 and εγ 1 = 0, that is, for harmonic excitations with frequencies ω = 0.7 and ω = 1.4,respectively.This is consistent with the findings of Figure 6.2.The noncoincidence, on the other hand, does not occur only in a neighborhood of the degenerate excitations, which correspond to the dashed line, but it extends to a very large range of combinations of εγ 1 and εγ 2 .This shows that higher-order analyses become always necessary in the case of generic excitations, where the simultaneous presence of different harmonics needs a more accurate analysis in order to detect their competitions and interactions.
The differences between the two analyses of Figure 6.5(a) are quantitatively modest (of the order of 5%), although they have strong consequences in terms of homoclinic bifurcations, as shown in Figures 6.3 and 6.4, and although they are important in the field of optimal chaos control [13].However, if we increase the excitation frequency, still remaining in the central range of medium frequencies, the differences become large, even reducing the damping, as shown by the example reported in Figure 6.5(b).They rapidly become huge by further increasing ω.
As a final comment, we note that first-order analysis systematically overestimates the "true" homoclinic threshold, a fact that can be justified by noticing that ε 2 M 2 (t 0 ) contains S. Lenci and G. Rega 165 only oscillating terms (see (6.2)), and thus its addition is likely to entail increasing the maximum of the whole oscillating part of the distance, and thus enhances the possibility to have homoclinic bifurcation.This does not occur only if the phases of the harmonics are properly synchronized, a phenomenon which is not robust and which does not occur in Figure 6.5.

Conclusions
A higher-order Melnikov analysis for single-DOF oscillators possessing homo-or heteroclinic loops has been developed in an abstract context.Classical Melnikov analysis has been extended in various directions: (i) by considering higher-order Melnikov functions, (ii) by measuring the distance on Poincaré sections nonperpendicular to the flux, and (iii) by computing the distance between perturbed and unperturbed manifolds ("onehalf " Melnikov functions).
The last two points are required for dealing with nonsmooth systems, and their relevance is shown by the application to the rigid block dynamics, which is in [9] and which shows how the present results permit to overcome the difficulties ensuing from nonsmoothness.
Concerning the first point, on the other hand, a recursive formula furnishing closedform expression of higher-order Melnikov functions is obtained without referring to any specific system.Then, the contributions of the second-order terms, which appear to be nonnegligible even for relatively small values of the perturbation parameters ε, have been investigated in given examples of smooth systems aimed at covering the main features of the problem (homo-versus heteroclinic loops, "one-half " versus "full" Melnikov functions, and so on).It has thus been proved that higher-order terms are sometimes necessary to extend to a range of practical interest the mathematical results, like the homoclinic bifurcation theorem [5, Theorem 4.5.3],holding for "ε sufficiently small."Detecting some of these cases constitutes one of the main achievements of the present work.
In fact, the analysis of the Helmholtz oscillator permits to draw some conclusions, which could be conjectured to hold in general.In particular, it has been shown that the classical (first-order) Melnikov analysis is inaccurate for detecting homoclinic bifurcation under harmonic excitation for low and high excitation frequencies, where the relevant threshold is overestimated.Furthermore, in the case of degenerate bifurcation ensuing from generic periodic excitations, it is inaccurate also for medium frequencies.More generally, it has been shown that higher-order terms are required for accurately measuring the interactions between various harmonics simultaneously exciting the system.
We do hope that the present investigation will permit to significantly enlarge the range of application of the classical Melnikov method, and to improve all investigations and achievements which are based on it.

Figure 6 . 2 .
Figure 6.2.(a)The theoretical homoclinic bifurcation thresholds εγ 1,cr .The thin lines correspond to εγ 1,cr calculated with first-order terms only (classical Melnikov analysis); the thick lines take into account also the second-order terms.(b) The ratio between the second-order and first-order theoretical predictions of εγ 1,cr .