A high-order, purely frequency based harmonic balance formulation for continuation of periodic solutions: The case of non-polynomial nonlinearities
aa r X i v : . [ phy s i c s . c l a ss - ph ] N ov A high-order, purely frequency based harmonic balanceformulation for continuation of periodic solutions : thecase of non polynomial nonlinearities
Sami Karkar a,b , Bruno Cochelin a,c , Christophe Vergez a a Laboratoire de M´ecanique et d’AcoustiqueCNRS – UPR 705131 chemin J. Aiguier13402 Marseille cedex 20FRANCE b Aix-Marseille University, FRANCE c ´Ecole Centrale MarseillePˆole de l’ ´Etoile, Technopˆole de Chˆateau-Gombert13451 Marseille Cedex 20FRANCE Abstract
In this paper, we extend the method proposed by Cochelin and Vergez in [A highorder purely frequency based harmonic balance formulation for continuation ofperiodic solutions,
Journal of Sound and Vibration , 324:243–262, 2009] to thecase of non polynomial non-linearities. This extension allows for the computa-tion of branches of periodic solutions of a broader class of nonlinear dynamicalsystems.The principle remains to transform the original ODE system into an ex-tended polynomial quadratic system for an easy application of the harmonicbalance method (HBM). The transformation of non polynomial terms is basedon the differentiation of state variables with respect to the time variable, shiftingthe nonlinear non polynomial non-linearity to a time-independent initial condi-tion equation, not concerned with the HBM. The continuation of the resultingalgebraic system is here performed by the Asymptotic Numerical Method (highorder Taylor series representation of the solution branch) using a further dif-ferentiation of the non polynomial algebraic equation with respect to the pathparameter.A one d.o.f. vibro-impact system is used to illustrate how an exponentialnon-linearity is handled, showing that the method works at very high order,1000 in that case. Various kinds of nonlinear functions are also treated, andfinally the nonlinear free pendulum is addressed, showing that very accurateperiodic solutions can be computed with the proposed method.
Keywords: nonlinear dynamical systems, asymptotic numerical method,harmonic balance, nonlinear dynamics, continuation
Email addresses: [email protected] (Sami Karkar), [email protected] (Bruno Cochelin), [email protected] (Christophe
Preprint submitted to Journal of Sound and Vibration August 25, 2018 . Introduction
In a previous issue, Cochelin et Vergez presented
A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions [1],targeting nonlinear dynamical systems. The technique addresses autonomousor forced nonlinear dynamical systems, using an arbitrary high order HarmonicBalance Method (HBM, see [2, 3]) together with the asymptotic numericalmethod (ANM, see [4, 5] for details) for the continuation.The proposed method automatically derives the HBM set of algebraic equa-tions to the desired order, provided that the governing equations are given in aquadratic formalism, ie, dynamical systems must be expressed as a set of first or-der nonlinear ordinary differential equations and nonlinear algebraic equations,all with purely quadratic polynomial nonlinearities at most (see eq. (3) in [1]).The method has been used in various applications, such as the analysis ofmicrobubbles in liquids by Pauzin et al. [6], or the nonlinear oscillations of nano-and micro-electromechanical devices by Kacem et al. [7]. A purely frequency-based stability analysis was also proposed by Lazarus and Thomas [8], as acompanion to the method presented in [1].Whereas classical harmonic balance technique is often limited to very fewharmonics, the proposed method has the advantage of allowing an arbitraryhigh order of the Fourier series. The difficulty is the need to recast the orig-inal system of equations (which is generally not quadratic), into a quadraticframework without changing the nature of the problem, ie, the recast should beonly a change of formlution not a change of the original problem. The recastof polynomial, square root, and rational functions has been treated in the citedpaper, defining intermediate variables and using quadratic polynomial algebraicequations. The authors also set up the basis for the recast of a sine nonlinear-ity, as encountered in the simple, free pendulum equation θ ′′ ( t ) + sin( θ ( t )) = 0(where the prime sign denotes time-derivation). They obtained a set composedof 4 nonlinear ordinary differential equations and 2 time-independent, nonlinearalgebraic equations. But, the calculation was not carried out.In this paper, we propose a generalisation of the method proposed in [1]for the treatment of ODE with non polynomial nonlinearities. This generalisa-tion greatly enlarges the field of application by overcoming the need for strictlyquadratic reformulation. Therefore, the method can be used to compute peri-odic solutions families of most nonlinear dynamical systems. Notice that oncea periodic solution is known, stability and bifurcation can be analyzed by usingclassical methods such as for instance, monitoring the eigenvalue of the mon-odromy matrix. This topic is not addressed hereafter.The present paper is organized as follows : in section 2, the case of theexponential function is discussed as an introductory example and the periodicsolutions of a one d.o.f vibro-impact system are calculated using up to 1000harmonics (with MANLAB software) ; then, in section 3, the method is gen-eralized to a wider class of nonlinearities and in section 4, we show how thisgeneralisation applies to the natural logarithm, trigonometric functions, as well Vergez)
2s non-integer power functions ; and last, in section 5, we use the proposedmethod to compute the periodic solutions of the simple free pendulum and wediscuss the accuracy of the solution.
2. An example using the exponential function: a regularised vibro-impact oscillator
We consider a one-degree-of-freedom, mass-spring oscillator which is limitedto the half-plane x < α to tune the wall stiffness.Let x ( t ) denote the mass position, we look for periodic solutions of : x ′′ ( t ) = − x ( t ) − lx ′ ( t ) − e α ( x ( t ) − , (1)where l is a free parameter (see A for details on this equation).Before embarking in the treatment of (1), we recall that two numerical meth-ods are used in [1] : a high-order HBM technique and a high order Taylor seriesexpansion for the continuation. The automation of these two techniques rely onthe quadratic form of the nonlinearities. To comply with this framework, themethod proposed by Cochelin et Vergez in [1] uses the following formalism: m ( Z ′ ( t )) = c + lc + l ( Z ( t )) + ll ( Z ( t )) + q ( Z ( t ) , Z ( t )) (2)where Z ( t ) is the unknown vector composed of time-dependent state variables, l is the continuation parameter, ( c , c ) are constant vectors, ( l , l , m ) are linearvector-valued operators, and q is a bilinear vector-valued operator.We will now present the successive recasts of system (1) needed to obtain therequired general form (2). It should be noticed in passing that the quadraticform (2) has nothing to do with a second order Taylor expansion. The pas-sage form (1) to (2) should not introduce any approximation, it is just anotherformulation of the same original problem. We introduce an additional variable y ( t ) = x ′ ( t ) and rewrite equation (1)as: x ′ ( t ) = y ( t ) (3a) y ′ ( t ) = − x ( t ) − ly ( t ) − e α ( x ( t ) − (3b) Introducing the additional variable e and its definition equation e ( t ) = exp h α (cid:0) x ( t ) − (cid:1)i (4)we differentiate it with respect to the time variable to obtain e ′ ( t ) = αx ′ ( t ) exp h α (cid:0) x ( t ) − (cid:1)i or equivalently e ′ ( t ) = αe ( t ) y ( t ) (5)3ote that for equation (5) to be exactly equivalent to (4), the following initialcondition must be added: e (0) = exp h α (cid:0) x (0) − (cid:1)i (6)We can now recast (3a)-(3b) into the following set of differential and algebraicequations: x ′ ( t ) = y ( t ) (7a) y ′ ( t ) = − x ( t ) − ly ( t ) − e ( t ) (7b) e ′ ( t ) = αe ( t ) y ( t ) (7c) e (0) =e α ( x (0) − (7d)At this point, the authors wish to underline the fact that the non polynomialnonlinearity has “moved” from an ODE equation (here: 3b) to an algebraicequation (here: 7d) which is not involved in the balance of harmonics: this isthe key point of the method. We will now apply the harmonic balance method (HBM) to the equations(7a)-(7c).Denote Z ( t ) = [ x ( t ) , y ( t ) , e ( t )] T , and the following operators: m ( Z ′ ( t )) =[ x ′ ( t ) , y ′ ( t ) , e ′ ( t )] T c = c =[0 , , T l ( Z ( t )) =[ y ( t ) , − x ( t ) − e ( t ) , T l ( Z ( t )) =[0 , − y ( t ) , T q ( Z ( t ) , Z ( t )) =[0 , , αe ( t ) y ( t )] T The ODEs (7a-7c) now take the required form: m ( Z ′ ) = c + lc + l ( Z ) + ll ( Z ) + q ( Z , Z ) . Eq. (7a-7b) are treated classically, using the method proposed by Cochelinand Vergez in [1], by balancing for the H + 1 harmonics 0 ≤ h ≤ H of Z ( t ),where H is the chosen truncation order of the Fourier series of the vector of(time-dependent) variables Z ( t ): Z ( t ) = Z + H X h =1 (cid:0) Z h − cos( kωt ) + Z h sin( kωt ) (cid:1) (8)For equation (7c), only harmonics 1 ≤ h ≤ H need to be balanced. Thereason is that the balance of the harmonic zero, which consists in equating themean value of each side over a period, is not suitable here: both sides of theequation result from time-differentiation of a periodic quantity (namely e ( t )),and has mean value which is therefore null. Thus, we only apply the HBM forharmonics 1 ≤ h ≤ H , where H is the truncation order.4he unknown (time-independent) variables are gathered in the state vector U : U = (cid:2) { Z i } i =0 .. H , λ, ω (cid:3) , whose size is N u =3(2 H + 1) + 2.The algebraic equations resulting from HBM application to (7a-7b-7c) , arethen automatically put into the standard ANM quadratic formalism: L + L ( U ) + Q ( U , U ) = 0 (9)with a constant vector L , a linear vectir-valued operator L , and a quadraticvector-valued operator Q .Now counting the number of algebraic equations: • H + 1) for the full HBM applied to (7a)-(7b) • H for the partial HBM applied to (7c) • • y (0) = 0)we reach a total number of N e =3(2 H + 1) + 1 equations.Thus, we get N e = N u −
1, as is needed for a 1D family of solutions.The reader shall notice that this algebraic system of N e equations for N u unknowns may now be solved by any continuation algorithm, such as for instancethe very classical first order predictor with Newton corrector and pseudo arc-length parametrization.In the following, we aim at using an ANM continuation method becauseof its robustness as compared to predictor-corrector algorithm. However, thealgebraic system is not fully quadratic, ie, equation (9) is but not (7d), hence,an additional step is necessary. We will now address the last recast that concerns non polynomial, nonlinear, algebraic equation (7d), following the usual procedure presented in [5] and [9].In the continuation process, the N u unknowns are function of a path parameter a . By differentiation of the non polynomial equation with respect to that pathparameter a , quadratic equation can be obtained as shown below.For sake of clarity, let us define two new variables, e t and x t , that represent e (0) and x (0) respectively: e t = e + H X h =1 e h − x t = x + H X h =1 x h − . where e i (respectively x i ) denotes the i -th coefficient of the Fourier series of e ( t )(resp. x ( t )) as defined for Z in the equation (8).By differentiation, equation (7d) is strictly equivalent to: ∀ a ∈ R , d e t ( a ) = α e t ( a ) d x t ( a ) (10a) e t ( a =0) = e α ( x t ( a =0) − (10b)5here d e = d ed a denote differentiation with rexpect to a .Equation (10a) is now quadratic in ( U , dU ) and is suitable for an easy andefficient computation of the Taylor series used for the continuation. The readeris referred to appendix C for details on the formalism used to enter (10a–10b)in the MANLAB software. −4 −3 −2 −1 0 1 2−4−3−2−101234 x y = x ’ (a) −4 −3 −2 −1 0 1 2−4−3−2−101234 x y = x ’ (b) Figure 1:
Phase portrait of the regularised vibro-impact: family of periodic solutionsin the phase plane ( u, v ). (a): non stiff regularisation ( α = 20) with H = 50 harmonics. (b):stiff regularisation ( α = 200) with H = 1000 harmonics. Figure 1 shows samples from the family of periodic solutions computed fortwo cases: a non stiff regularisation, with α = 20, using H = 50 harmonics (a); and a stiff regularisation, with α = 200, using H = 1000 harmonics (b). Inboth cases, the continuation was run with a threshold of 10 − on the residue.The phase portrait cycles are to be compared with those of a free, conserva-tive non-regular vibro-impact system, i.e. a wall modelled with an impact lawusing a restitution coefficient equal to unity : the family of periodic solutions iscomposed of origin-centered circles, when the amplitude is less than unity, andorigin centered arc of circles closed by a vertical segment along the line x = 1,when the amplitude is higher than unity.
3. General treatment of nonlinear functions
Here, we discuss the general method for the recast of most nonlinearitiesinto quadratic formulation. Note that the quadratic recast of rational functionshas been given by Cochelin and Vergez in [1].Let us consider a set of differential and algebraic equations: F ( u, v, g ( u )) = 0 (11)6here v = u ′ , F is at most quadratic in its arguments, and g is any nonlinearfunction of u . We add a new variable w defined by w = g ( u ). By time-derivation, we obtain: w ′ = ∂g∂u ( u ) v (12a) w (0) = g ( u (0)) (12b)If ∂g/∂u can be written as a rational function of ( u, v, w ), then the equation x = ∂g/∂u can be recast into quadratic equations (possibly using additionalvariables). Consequently, the time-dependent equation (12a) becomes w ′ = xv ,which is quadratic in x and v .We apply the HBM to this equation, but only for harmonics 1 ≤ h ≤ H ,while the mean value (harmonic zero) will be constrained by the initial condition(12b). The equation count is then equal to 2 H + 1, as in a standard HBMapplied to a unique equation, which matches the number of variables : the2 H + 1 coefficients of the Fourier series of w ( t ) up to harmonic H . In the case ofautonomous systems, the period (or equivalently the angular frequency) is alsounknown and one need to add a phase equation, as explained by Doedel in [10]as well as Cochelin and Vergez in [1].If additional variables were to be used for the quadratic recast of x , allcorresponding algebraic equations will be treated with full HBM, including thebalance of the harmonic zero (the mean values). If x = ∂g/∂u cannot be expressed as a quadratic polynomial of the currentvariables, we differentiate it with respect to t : x ′ = ∂ g∂u ( u ) v (13a) x (0) = ∂g∂u ( u (0)) (13b)If ∂ g/∂u can be expressed as a rational function of ( u, v, w, x ), then wedefine y = ∂ g/∂u and the previous results apply: the time-dependent equationin (13a) is quadratic in y and v .We then have two differential equations, namely w ′ = xv and x ′ = yv , withtwo associated initial conditions. The differential equations are treated usingHBM for harmonics 1 and higher, while the initial conditions will constrain themean values of w and x (the nonpolynomial, nonlinear, algebraic equation isaddressed by differentiation, as explained in section 2.4).If additional variables were used for the recast of the rational function y , allcorresponding algebraic equations will be treated with full HBM, including thebalance of the harmonic 0 (the mean values).
4. Recast of a few common non-polynomial nonlinearities
For the quadratic recast of the exponential function, the reader is referredto section 2. 7 .1. Natural logarithm
Given w defined as w ( t ) = ln( u ( t )) (assuming u > w ′ = u ′ /u ,or equivalently, using x = w ′ : (cid:26) u ′ = xuw ( t = 0) = ln( u ( t = 0))which is quadratic in u and x . Given w ( t ) = u ( t ) α where α ∈ R is a constant, then w ′ = αu α − u ′ . Using v = u ′ and x = w ′ , one gets: (cid:26) ux = αwvw ( t = 0) = u ( t = 0) α which is quadratic in ( u, v, w, x ). Given s ( t ) = sin( u ( t )) and c ( t ) = cos( u ( t )), we introduce v ( t ) = u ′ ( t ) andtime-derivation of the definition equations of s and c gives: s ′ = cvc ′ = − svs ( t = 0) = sin( u ( t = 0)) c ( t = 0) = cos( u (( t = 0))which are obviously quadratic in ( c, v ) and ( s, v ) respectively.As for w ( t ) = tan( u ( t )), time-differentiation leads to w ′ = (1 + w ) v . Usingan additional variable x = 1 + w , one gets the quadratic equation: (cid:26) w ′ = xvw ( t = 0) = tan( u ( t = 0)) .
5. Periodic solutions of the simple, free, nonlinear pendulum
Denoting θ the angle between the current position and the lower, rest posi-tion, the equation of motion of the free pendulum writes: θ ′′ ( t ) + sin (cid:2) θ ( t ) (cid:3) = 0 (14)Adding a damping parameter l , that will vanish along the family of periodicsolutions, as explained in section A, the equation of motion becomes: θ ′′ ( t ) + lθ ′ ( t ) + sin (cid:2) θ ( t ) (cid:3) = 0 (15)8sing three additional variables, v ( t )= θ ′ ( t ), s ( t )=sin (cid:2) θ ( t ) (cid:3) and c ( t )=cos (cid:2) θ ( t ) (cid:3) ,we recast the equation (15) into the following quadratic, differential-algebraicsystem: θ ′ ( t ) = v ( t ) (16a) v ′ ( t ) = − s ( t ) − lv ( t ) (16b) s ′ ( t ) = c ( t ) v ( t ) (16c) c ′ ( t ) = − s ( t ) v ( t ) (16d) s (0) = sin (cid:0) θ (0) (cid:1) (16e) c (0) = cos (cid:0) θ (0) (cid:1) (16f)Using the proposed method, we apply: • the full HBM to equations (16a-16b) • the balance of the harmonics 1 and higher to equations (16c-16d)and since the frequency is also unknown, we add a phase equation: v (0) = 0 . (17)Finally, the initial conditions (16e-16f) are differentiated with respect tothe path parameter. The equations are then in the right (quadratic) form forapplying the ANM continuation, ie, computing a high orde Taylor series of theunknowns with respect to the path parameter: (cid:26) L + L ( U ) + Q ( U , U ) = ( U . ) = Qh ( U , U . )The total number of equation is then: N e = 2(2 H + 1) | {z } full HBM + 2(2 H ) | {z } HBM h> + 1 | {z } phase eq. | {z } in L0/L/Q + 2 | {z } in Lh/Qh = 4(2 H + 1) + 1while the number of variables is N u = 4(2 H + 1) | {z } Fourier components of θ,v,s,c + 1 |{z} l + 1 |{z} ω = 4(2 H + 1) + 2Figure 2 shows the phase portrait and the amplitude-frequency diagram ofthe periodic solutions family of the nonlinear, free pendulum system, computedwith H = 100 harmonics. The theoretical amplitude-frequency diagram of thissystem (dashes) corresponds to the following analytic formula: T ( θ max ) = 2 π K (cid:16) sin ( θ max / (cid:17) , where K is the complete elliptic integral of the first kind (see [11] for instance).The computed curve (plaine line) is superimposed on the theoretical one.9 θ v = θ ’ (a) θ max ω (b) analyticcomputed Figure 2:
Periodic solutions of the nonlinear free pendulum. (a): phase portrait ( θ, v )showing a few samples of the periodic solution family. (b): amplitude-frequency diagram,the computed curve is superimposed on that of the analytic formula. HBM using H = 100harmonics, ANM-threshold: 10 − , the residue || R || was kept under 10 − . The branch was computed with 29 step of continuation. The computationwas stopped when the relative error between the computed angular frequencyand its theoretical value (given the amplitude) reached 0 . (cid:18) θ max ω (cid:19) = (cid:18) . π rad0 . (cid:19) where the norm of the residual vector is || R || =2 .
76 10 − . Increasing the numberof harmonics only leads to a closer approach to the limit point ( θ max , ω )=( π, .Notice that two additional variables must be introduced for this simple oned.o.f examle. For a multiple d.o.f. system with many kinds of non-polynomialnonlinearities, the quadratic recast has to be applied to each individual non-linear term. This may require a great number of additional variables and thetransformed system may be much larger than the original one. This is howeverthe price to pay with this method, the more complex the original system, thelarger the transformed quadratic system.
6. Conclusion
In this paper, we extended the works of Cochelin and Vergez presented in[1] to the case of general nonlinearities.A vibro-impact with exponential regularization was presented, and its peri-odic solutions were computed with the proposed generalization. The case of avery stiff regularization demonstrated the capabilities of this numerical tool todeal with a very high number of harmonics in the harmonic balance method.Finally, we showed how to apply the method for the continuation of theperiodic solutions of the simple, free, nonlinear pendulum. Our results confirmthat the harmonic balance method with a high number of harmonics is bothaffordable and well suited to the ANM continuation framework. It has been This corresponds to a heteroclinic connection between the two saddle-nodes( θ SN , v SN ) ∈ { ( − π, , ( π, } ppendix A Vibro-impact with exponential wall reaction A.1 Model
We consider a one-degree-of-freedom, mass-spring oscillator which is limitedto the half-plane x < x ( t ) denotes the position of themass.The rigid wall reaction is modelled by an exponential function, with a coef-ficient α to tune the wall stiffness : F w ( x ) = − e α ( x − , Thus, the regularised vibro-impact system is governed by the following equa-tion: x ′′ ( t ) = − x ( t ) − e α ( x ( t ) − , (A.1)where the prime sign denotes time-differentiation. The force F w ( x ) that reflectsthe wall effect derives from a potential energy so that problem (A.1) keeps theproperty of being energy-conservative. A.2 Dissipative recast for continuation
Mu˜noz-Almaraz et al. [13] showed that, in conservative Hamiltonian sys-tems, periodic orbits generally belong to a one-dimensional family of periodicsolutions, parametrised by the value of the first integral (here, the total energy),which is not an explicit parameter of the system. To compute this family ofperiodic solutions in the standard continuation framework R ( x ( t ) , l ) = 0, weperturb the initial equation with a damping term added to the right-hand sideof (A.1). The system is then embedded into a general, dissipative system: x ′′ ( t ) = − x ( t ) − lx ′ ( t ) − e α ( x ( t ) − , (A.2)where l is a free parameter of the continuation.The perturbed system (A.2) possesses periodic solutions that are exactlythose of the unperturbed, conservative system (A.1), if and only if l = 0. Thisway, the additional parameter l allows us to compute the periodic solutionsof the conservative system (A.1) using the classical framework for dissipativesystems possessing an explicit control parameter. Appendix B Classical ANM: quadratic framework
The reader is referred to Cochelin and Vergez [1] (sections 2.4 and 2.5,pp.248-250), for the details concerning the principle and the implementationof the ANM in the classical, quadratic framework.
Appendix C Extended ANM framework
C.1 Principle of the series computation
Given a nonlinear system f ( U ) = 0 with N d equations and N u unknowns,whose differentiated form reads: Lh ( U . ) = Qh ( U , U . ) (C.1)12here Lh : R N u → R N d is a linear, vector-valued operator ; Qh : R N u × R N u → R N d is a bilinear, vector-valued operator.Assuming a known regular solution U of this system, we write the branchof solutions passing through this point as a Taylor series expansion: U ( a ) = U + a U + a U + a U + · · · + a n U n . (C.2)where the branch is parametrised using the pseudo-arclength parameter a de-fined as: a = ( U − U ) t U . (C.3)Differentiating U reads: U . = { U + 2 a U + 3 a U + ... + na n − U n } da. (C.4)Then, substituting both (C.2) and (C.4) in system (C.1) and equating eachpower of a (up to order n ) to zero gives: • power 0: Lh ( U ) − Qh ( U , U ) = 0, which can also be written Jh U . U = 0where Jh U ∈ R N d × N d +1 is the jacobian matrix of f evaluated at U . Thislinear equation in U thus gives the term of order • power 1 ≤ p ≤ n − Jh U . U p +1 = Σ pi =1 p +1 − ip +1 Qh ( U i , U p +1 − i ). Thislinear equation gives the term of order p + 1 of (C.2).The original nonlinear problem is thus replaced by a cascade of n linear systems,which all share the same matrix: Jh U .However, at each order, the linear systems have N d + 1 unknowns and only N d equations. For each linear system, the additional equation is obtained bysubstituting the series (C.2) into the definition (C.3) of the path parameter a : • order 1 : U t U = 1 • order 2 ≤ p ≤ n : U t U p = 0 C.2 Implementation in MANLAB: the example of the vibro-impact
As for the classical framework, the only user input to the MANLAB softwareconsists in M-functions for the operators Lh and Qh , as well as for f (for theresidue computation only), and a starting point U .In the case of the vibro-impact system presented in section 2, some equationsare quadratic (those resulting from the HBM and those for the definition of e t and x t ) while the last one is not. We thus separate the equations in two parts:those resulting from the HBM, that will appear in the L L and Q operators,and the last one, that will appear in the Lh , Qh operators.For the present problem, the state vector is: U = ( x , y , e | {z } Z , x , y , e | {z } Z , x , y , e | {z } Z , · · · ,x H − , y H − , e H − | {z } Z H − , x H , y H , e H | {z } Z H , l, ω, e t , x t )and its size is N u = 3(2 H + 1) + 4. 13or the quadratic part of the problem, the subsystem contains the followingnumber of equations: N e = 2(2 H + 1) | {z } full HBM + 2 H | {z } HBM h =0 + 2 | {z } def. of e t ,x t + 1 | {z } phase eq. = 3(2 H + 1) + 2 . The content of
L0.m , L.m and
Q.m are not listed here, as it is the direct result ofthe harmonic balance as presented in [1]. These three functions return a vectorof size N e − N d = 1 equation. Thecontent of Lh.m , Qh.m and f.m is listed below, with α = 200: function [Lh] = Lh(dU)Lh=zeros(1,1);Lh = dU(end-1);function [Qh] = Qh(U,dU)Qh = zeros(1,1);Qh = 200*U(end-1)*dU(end);function [f] = f(U)f = zeros(1,1);f = exp(200*(U(end)-1)); eferences [1] Bruno Cochelin and Christophe Vergez. A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions. Journal of Sound and Vibration , 324:243–262, 2009.[2] Minoru Urabe. Periodic solutions of differential systems, galerkin’s pro-cedure and the method of averaging.
Journal of Differential Equations ,2(3):265 – 280, 1966.[3] M. Nakhla and J. Vlach. A piecewise harmonic balance technique for deter-mination of periodic response of nonlinear systems.
Circuits and Systems,IEEE Transactions on , 23(2):85 – 91, feb 1976.[4] Bruno Cochelin. A path-following technique via an asymptotic-numericalmethod.
Computers and Structures , 53:1181–1192, 1994.[5] Bruno Cochelin, Noureddine Damil, and Michel Potier-Ferry.
M´ethodeAsymptotique Num´erique . Lavoisier, Paris, 2007. (Asymptotic NumericalMethod).[6] Marie-Christine Pauzin, Serge Mensah, Bruno Cochelin, and Jean-PierreLefebvre. High order harmonic balance formulation of free and encapsulatedmicrobubbles.
Journal of Sound and Vibration , 330(5):987 – 1004, 2011.[7] N. Kacem, S. Baguet, S. Hentz, and R. Dufour. Computational and quasi-analytical models for non-linear vibrations of resonant mems and nemssensors.
International Journal of Non-Linear Mechanics , 46(3):532 – 542,2011.[8] A. Lazarus and O. Thomas. A harmonic-based method for computing thestability of periodic solutions of dynamical systems.
Comptes Rendus deM´ecaniques , 338:510–517, 2010.[9] R´emy Arquier, Bruno Cochelin, Sami Karkar, Arnaud Lazarus, OlivierThomas, and Christophe Vergez.
MANLAB 2.0, an interactive contin-uation software , 2010. http://manlab.lma.cnrs-mrs.fr (last visited28/06/2011).[10] E. J. Doedel. Lecture notes on numerical analysis of nonlinear equa-tions, 2010. http://indy.cs.concordia.ca/auto/notes.pdf (last vis-ited 28/06/2011) 391 slides.[11] A. Belendez, C. Pascual, D. I. Mendez, T. Belendez, and C. Neipp. Ex-act solution for the nonlinear pendulum.
Revista Brasilieira de Ensine deFisica , 29(4):645–648, Between October and December 2007.[12] Isabelle Charpentier and Michel Potier-Ferry. Diff´erentiation automa-tique de la m´ethode asymptotique num´erique typ´ee : l’approche diamant.
Comptes Rendus M´ecanique , 336(3):336 – 340, 2008. (Automatic differen-tiation of asymptotic numerical method : the DIAMANT approach).[13] F. J. Mu˜noz-Almaraz, E. Freire, J. Gal´an, E. Doedel, and A. Vander-bauwhede. Continuation of periodic orbits in conservative and hamiltoniansystems.
Physica D: Nonlinear Phenomena , 181(1-2):1 – 38, 2003.15 ist of Figures Phase portrait of the regularised vibro-impact: family ofperiodic solutions in the phase plane ( u, v ). (a): non stiff regular-isation ( α = 20) with H = 50 harmonics. (b): stiff regularisation( α = 200) with H = 1000 harmonics. . . . . . . . . . . . . . . . . 62 Periodic solutions of the nonlinear free pendulum. (a):phase portrait ( θ, v ) showing a few samples of the periodic so-lution family. (b): amplitude-frequency diagram, the computedcurve is superimposed on that of the analytic formula. HBM us-ing H = 100 harmonics, ANM-threshold: 10 − , the residue || R || was kept under 10 − . . . . . . . . . . . . . . . . . . . . . . . . . 1016 θ t (rad) ω θ maxmax