im Journal of JET v°lume i°<2oi7)p-p- 57-69 Issue 3, October 2017 Energy Type of a^de 1XI:L Technology www.fe.um.si/en/jet.html ON THECENTRE-FOCUSPROBLEM IN LRCELE L^NARD SUTLEU L R PPROBLEIUCENTISUNFOKUSA V ELKATLPIH LIEEARDOVIH SISLEMIH MayejMenciscer1'2* Keywords:Polynomlal systemnf ODEs, CentercCocusproblem,Complnx benavinr Abstract A large family of planar systems of ODEs arising from Lienard equations is considered. A Lienard equation x + f(x)x + g(x) = 0 is commonly used in practical problems, in particular in (electro)mechanics. It is well-known that a Lienard equation can be transformed into an autonomous planar system of ODEs of the form x' = y — F(x), y' = ~gix), where F'(x) = f{x). In this paper fix) = 2a2x + 3a3x2 + 4a4x3 and g{x) = c1x + c3x3 + csxs + c7x7. In the parameter space (a2,a3,a4,c1,c3,cs,c7)E R7 we consider the center-focus problem and find necessary conditions for the corresponding system having a center at the origin. In the parameter space (a2,a3,a4,c1,c3,cs,c7)E R7 an example with a possible limit cycle and some examples with other complex dynamic behavior are presented. R Corresponding author: Matej Mencinger, Tel.: +386 2 2294 321, Mailing address: University of Maribor, Faculty of Civil Engineering, Transportation Engineering and Architecture, Smetanova ulica 17, 2000 Maribor, Slovenia. Email address: matej.mencinger(S)um.si 1 University of Maribor, Faculty of Civil Engineering, Transportation Engineering and Architecture, Smetanova ulica 17, 2000 Maribor, Slovenia. 2 Institute of Mathematics, Physics and Mechanics, Jadranska 19, 1000, Ljubljana, Slovenia. JET 22 Matej Mencinger JET Vol. 10 (2017) Issue 3 Povzetek Obravnavana je velika družina ravninskih sistemov navadnih diferencialnih enačb (NDE), ki izhaja iz Lienardove enačbe. Lienardova enačba x + f(x)x + g{x~) = 0 se uporablja v praktičnih problemih, še zlasti v (elektro)mehaniki. Znano je, da se Lienardova enačba lahko preoblikuje v avtonomni ravninski sistem NDE oblike x' = y — F{x), y' = —g{x), kjer je F'{x) = f{x). V tem članku sta f{x) = 2a2x + 3a3x2 + 4a4x3 in g(x) =c1x + c3x3 + csxs + c7x7. V prostoru parametrov (a2,a3,a4,c1,c3,cs,c7)ER7 obravnavamo problem centra in fokusa in poiščemo potrebne pogoje, da ima ustrezen sistem center v izhodišču. Prav tako v prostoru parametrov (a2,a3,a4,c1,c3,cs,c7)ER7 predstavimo primer z možnim limitnim ciklom in nekatere primere z drugimi kompleksnimi dinamičnimi pojavi. 1 INTRODUCTION A great number of mathematical models of physical systems give rise to differential equations of the type x + f{x)x + g{x) = 0 (1.1) which is called a Lienard equation. From the mechanical point of view, equation (1.1) can be interpreted as the generalization of the mass-spring-damper system, where f(x)x is the damping term and g(x) represents the (nonlinear) spring term. Applications of equation (1.1) can be found in many important examples including chemical reactions, predator-prey models, vibration analysis, etc. Two famous examples of a Lienard equation are the Van der Pol equation and Duffing's equation. The Van der Pol equation x + e(x2 — l)x + x = 0, £ > 0 describes the circuit of a vacuum tube, whilst Duffing's equation x + Sx + ax + @x3 = ycos(_wt), aims to model certain nonlinearly damped/driven oscillators (i.e. a spring pendulum whose spring's stiffness does not exactly obey Hooke's law). Here x = x{t) represents the displacement of the (pendulum) bob at time t, x represents the first derivative of x with respect to time t, and x is the second time-derivative of x. Parameters a,p,y,S and w >0 are given (real) constants (case p = S = 0 corresponds to simple harmonic motion). There are two conventional transitions from homogeneous ODE (1.1) to a planar dynamical system of ODEs. Namely, setting y = x we obtain x' =y, y' = —f(x)y — g(x). (1.2) In (1.1) another approach is possible via the so-called Lienard coordinates. Substituting y = x + F(x), where f(x) =F'(x), one obtains 58 JET Onthe centre-focus problemin some Lienord systems x' = y — F(x), y' = —g(x). (1.3) Both systems (1.2) and (1.3) are special cases of (continuous) dynamic system called autonomous systems of ODEs, which generally takes the form x' = P(x,y), y' = Q(x,y). (1.4) Even in planar dynamics, there are several open problems: among them the problem of finding the position and the number of limit cycles bifurcating from a non-hyperbolic singular point, which is a part of famous Hilbert's 16th problem. Most real-life problems are related to the centre-focus problem: this is of distinguishing between a centre, where all orbits in the neighbourhood of the singular point are periodic, and a focus, where all orbits are spiralling away or towards to the singular point/origin (for more details see e.g. [1]). Autonomous systems (1.4) cannot contain chaotic dynamics (according to the Poincare-Bendixon theorem). However, note that Duffing's equation (which is a nonhomogeneous ODE of order two) corresponding to a non-autonomous system x' = P(x,y,t), y' = Q(x,y, t)) is an example of a dynamical system that exhibits chaotic behaviour. In contrast Duffing's equation is also a classic example of a dynamical system with a limit cycle, [2,3]. For example, for /(x) = £(x2 — 1) and g(x) =x we obtain which readily contains periodic solutions for x(t) and y(t). In Figs.1-3, there is an example for £ = 1.2 with initial conditions x(0) = —0.4, y(0) = 0.3. In Fig. 4, the corresponding stream-line plot with a clearly visible limit cycle is presented. JET 59 Matec Mencinger JET Vol. 10 (2017) Issue 3 Figure 2: Initial value solution y = y(t) to Duffing's equation: £ = 1.2; x(0) = —0.4, y(0) = 0.3 Figure 3: A solution tending to limit cycle of Duffing's equation: £ = 1.2 and x(0) = —0.4, y(0) = 0.3 \ \\ wv Vw5 V t K \ y X^vlT1-- k. v v ---- 77) > ^J \ \ \ a \ ^ i \ \ \1 \ k * \ ^A \ V V ' < V -V \ \ V * \ 1 1 j v v y \ Aa V K * v V \ v\ v x\x x O V Figure 4: Streamline-plot (Mathematica) of Duffing's equation for £ = 1.2 containing a limit cycle 66 JET Onthe centre-focus problemin some Lienord systems Usually for systems (1.4) and, consequently, also for systems (1.2) and (1.3) singular points xo are hyperbolic; this means that no eigenvalue of the Jacobian matrix /&>) = 3P dPi 3y 3y. , wherei = (x,y) evaluated at xo has the real part equal to zero. In this case the linearized system x' = /(io) may (locally - near singular point xo = (xo,yo)) be used to approximate the original system (1.4). More precisely, the approximation is in terms of a continuous invertible map, which locally (near singularity) takes parametrized solutions of the linearized system x' = /(io) to the parametrized solutions of the original system (1.4). This is the statement of the Hartman—Grobman theorem [4]. For any singular point xo of system (1.4) one of the following five »generic cases«, according to the Jordan canonical form of/(io), appears (i) 0 A.*O(IO[—>„ a«-«!» a^e SI,"Owb 0]. None, except the first one, are within the scope of the Hartman—Grobman theorem (for any w ^ 0) and must be considered case by case/separately. Systems corresponding to the Jacobian matrix of the form (ii) are the most studied planar systems. The addition of nonlinear terms may result either in centre or in focus. In this paper, we will consider a special case of (1.2) of the form x' = y, y' = —(2a2x + 3a3x2 + 4a4x3)y — (c1x + c3x3 + c5x5 + c7x7), (1.5) where a2,a3,a4,c1,c3,c5,c7 El, and analyse the stability for the whole family (1.5). 2 THE ANALYSIS OF SYSTEM (1.5) 2.1 Case >0 If c1 >0 introducing new coordinates X = c1x, Y = c1y one can obtain system X' = Y, Y' = —(2a2X + 3a3X2 + 4a4X3)Y — (X + c3X3 + c5X5 + c7X7). The corresponding Jacobian at singular point (0,0) is 7«M>) = [—1 0] which yields a centre or focus at the origin. In Fig. 5 there is a stream-line plot of a centre. JET 59 Matec Mencinger JET Vol. 10 (2017) Issue 3 O \\ \ " * \ \\ 1 Figure 5: Centre for system (5), with a2 = a3 = 0, a4 = ^,C1 = l,C3 = CS = 0, C7 = 2 Note that the complete analysis of case c1 >0 is done in a separate section in the sequel; it will be seen that a3 = 0 is the sufficient condition to obtain a centre in system (1.5). However, distinguishing between a centre and a focus just from a streamline plot (e.g. in Mathematica) is impossible and much too inaccurate (since the stream-plot of a focus is too similar to the stream-plot of a centre). In the case of a focus, it is much more convenient to consider a single solution, like in Figs. 6-7 in which graphs of x = x(t) and y = y(t) are shown. In Fig. 8, the parametric solution (x(t),y(t)) is shown in the phase plane (x,y). Figure 6: Particular solution x = x(t) of (1.5) with a2 = a4 = 0, a3 = C1 = l, C3 = 2,CS = C7 = 0 66 JET Onthe centre-focus problemin some Lienord systems Figure 7: PaetitrSae ncSreich y — y(t) cf (1.5) with a2 = a4 = 0, a3 = C1 — 1, C3 — 2,C5 =C7 = 0 Figure 8: Trajectory (x(t),y(t)) of (1.5) with a2 = a4 — 0, a3 = 1,C1 — 1, C3 = 2, Cs = C7 = 0 2.2 Case ct <0 If c1 <0 introducing new coordinates X = C1x,Y = C1y one yields X' = Y, Y' = —(2a2X + 3a3X2 + 4a4X3)Y — (—X + c3X3 + csX5 + c7X7). The corresponding Jacobian matrix at a singular point (0,0) is '«m»-G 1 The eigenvalues of J(0,0) are Xi2 = +1, yielding a hyperbolic singular point: a saddle with an unstable singularity. The dynamics near the origin (according to the Hartman—Grobman theorem) are topologically conjugate to the linear system X' = Y, Y' —X (see Fig. 9), and no further analysis is needed in this case. JET 59 Matec Mencinger JET Vol. 10 (2017) Issue 3 Figure 9: Phase portrait of a saddle for (1.5); a2 = a4 = 0, a3 = 1,C1 = —1,C3 = CS = 0,C7 = 1 2.3 Case c1 = 0 If c1 = 0, the system takes the form x' = y, y' = —(2a2x + 3a3x2 + 4a4x3)y —(c3x3 + csxs + c7x7), and (0,0) is a non-elementary (nilpotent) singular point if type (iv), since the corresponding Jacobian at singular point (0,0) is '(00»=i0 a- The blow-up desingularization gives rise to the type of stability in this case (see [5,6] for details). dP Note that the case (v) can never appear for Lienard systems since — =1^0. 3 THE CENTRE-FOCUS ANALYSIS In order to consider the centre-focus problem, we briefly recall an approach [1] for studying the problem for planar polynomial systems (1.4), where P(x,y) = y + h.o.t. and Q(x,y) = —x + h. o. t. Here h. o. t. stands for Y,1i+j=2 PijxiyJ, and "Z1i+j=2 Qi,jxlyJ, respectively. Note that via the (reverse) time transformation t = —t this is equivalent to P(x,y) = y + h.o.t. and Q(x,y) = —x + h.o.t., therefore we make no difference between these two cases. Here P(x,y) and Q(x,y) are polynomials of degree at most n without constant and linear terms. It is convenient to introduce the polar coordinates u = r cos (p ,v = r sin q> and consider the so-called Poincare return map R(r). Introducing the polar coordinates into (1.4) with P(x,y) = y + h.o.t. and Q(x,y) = —x + h.o.t. yields the equation of the trajectories dr r2F(r,cosp ,sinp) dp 1 + rG(r,,cosp ,sinp) ' 66 JET Onthe centre-focus problemin some Lienord systems dr Obviously, r = 0 is a solution to — = ft(r, since ft(0, = 0. The function ft(r, is periodic (with the least period in variable and analytic for (small enough) |r| < r* (and all Thus, ft(r, can be expanded in a convergent power series in r to obtain The continuous dependence on the initial conditions and the fact that r = 0 is a solution for all ^ e [0,2tc] yield that every solution to the above equation in a sufficiently small neighbourhood of the origin intersects every ray ^o e [0,2^]. This means that without loss of generality one can choose the line segment I = {(u, v); v = 0,0 < u < r* }, where r* is chosen to be small enough. Next, we consider the solution of (3.1) with the initial condition r(^ = 0) =r0 and expand it into a power series in r0 to obtain r(^,ro) = w1(^)ro + w2(^>o2 + w3(^)ro3 + ••• (3.2) which is also convergent for all ^ e [0,2tc] and all |r0| < r*. The r(^) from (3.2) is a solution of (3.1) and inserting r(^,ro) into (3.2) yields recurrence differential equations for functions (see [1] for more details). We consider one revolution of r = r(^,ro) beginning on ro e I where ^ is assumed to be 0 and study the return to I (which occurs at ^ = that is, after one revolution). Thus, the Poincaré return map ft(r) is defined by ^(ro) = r(2^,ro) = r +M/2(2H>o2 + W3(2re)ro5 + ••• The coefficients = w;(2h) for 7 > 1, defined in the above equation are called Lyapunnv hUdUcen. From the definition of polar coordinates, we readily conclude that the zeros of the difference function ^(ro) — ro correspond to closed orbits. In particular, isolated zeros correspond to limit cycles, and if ^(ro) — ro = 0, the system has a centre at the origin, which means that for all y> 1 the Lyapunov numbers must vanish. However, computing Lyapunov numbers requires the integration of trigonometric functions, which can be very difficult problems for some cases. Poincare and Lyapunov proved that system (1.4) with P(x,y) = y + h. o. t. and Q(x,y) = — x + h. o. t. has centre at the origin if it admits the first integral of the form $>(x,y) =x2 + y2 + ^ v xly;, ¿+;>3 (3.3) which is an analytic function in a neighbourhood of the origin (0,0). By definition, O is a first integral of system (1.4) if O is a solution to the following PDE dO(x,y) N dO(x,y) „ (3.4) ox oy In agreement with formula (3.4), equation ^ = 0 can be solved only on some special variety (set of zeros) in the (affine) space of parameters defined by coefficients of P(x,y) and Q(x,y). Generally, from (3.4) using step-by-step process of equating the proper coefficients of ^ to zero we obtain JET 59 Matec Mencinger JET Vol. 10 (2017) Issue 3 ¥ = gi{alj,blj){x2 + y2)2 +g6(aiJ,biJ)(x2 + y2)3 +gB(aiJ,biJ)(x2 + y2)4 + ••• yielding ¥ = 0 if and only if 9i{atj,btj) = gb{atj,btj) = • = g2k{atj,btj) = ••• = 0 (3.5) for all k >2. The numbers (polynomials) g2k{aij,bij) are called focus quantities [1]. Finally, the centre-focus problem reduces to find the conditions for vanishing all focus quantities g2k{atj,btj). For case (1.5), we computed the first several polynomials g2k{aij,bij) using and obtained 3a? g* = ~— g6=^(76a22+21c3), gB = -3a2(33Q34a% + 19736a^c3 + 9(133a2. + 229c% - 216cs) - 15728a2a4), g10 = -a^(15467244a% - 11734288a^a4 + 10684399a^c3 - 2776152a2a4c3 - 184320 3a2(484752al - 789569cj + 559144cs) - 3a\(484752a2 - 789569c2 + 559144cs)+9(52064a24 + 5(4929a2c3 + 2901c3 - 6164c3cs + 2568c7))), etc. The expressions for the other focus quantities are too long to be presented here, but can easily be computed using Mathematica or any other appropriate computer algebra system. Obviously, the condition a3 =0 is necessary for (3.5), since a3 is a cofactor of g2k{alj,blj^ for any k. This means that system (1.5) for a3 = 0 and c1 >0 has centre at the origin for arbitrary a2,a4,c3,cs,c7 El and the phase portrait in the neighbourhood of the origin is topologically equivalent to the phase portrait shown in Fig. 5. Figure 10: A eoe-tnivial limit cycle fon (1.2) of the form x' = y,y' = -(x2 - 3)y - 5x3(x2 - 5) 66 JET Onthe centre-focus problemin some Lienord systems ^ 2_ Figure 13: CcmpScx dynamite cf x' =y,y' = (--X2)y----- — X + 2 sin 3t: (*(t),y(t)) JET 59 Matej Mencinger JET Vol. 10 (2017) Issue 3 ^ *3(*2_ Figure 14: Complex dynamics of x' = y,y' = (--x )y----- — X — 3 sin It: (x(t),y(t)) Figure 15: A single trajectory of x' = y,y' = (--X2)y----- — x — 3sin2t2: (x(t),y(t)) 4 CONCLUSIONS The centre-focus problem for system (1.5) is solved for necessary conditions. It is proven that the crucial parameter that distinguishes between the centre and focus in (1.5) is a3. The condition a3 ^ 0 yields focus, while a3 = 0 defines a centre for any choice of other parameters. Note also that there are many nontrivial examples of system (1.2) that probably admit limit cycles, for example system x' = y,y' = —(x2 — 3)y — 5x3(x2 — 5)) shown in Fig.10. According to the well-known Lienard theorem, [3], system (1.2) admits a unique stable limit cycle if f and g are continuously differentiable, g is odd and g(x) >0 for x >0, f is an even function and f* f(u) du < 0 for 0 0 and increasing for x >a; for some a El+. Finally, note that a non-autonomous modification of system (1.2): x' = y, y' = —f(x)y — g(x) + h(t) may exhibit other complex dynamics, as shown in Figs. 11-15. 68 JET Onthe centre-focus problemin some Lienord systems References [1] V.G. Romanovski, D.S. Shafer: Thc Cchtce and Cytlitiee PecUlcme: A Computational AlgcUea Appecat, Birkhauser, Boston, 2009 [2] A. Constantin: On the Oscillation of Solutions of the Liénaed Equation, J. Math. Ann. Appl., Vol. 205, p.p. 207-215, 1997 [3] A. Palit, D.P. Datta: On a Finitc NrmUce cf Limit Cytlce in a Liénaed Systcm, Int. J. Pure and Applied Math., Vol. 59, p.p. 469-488, 2010 [4] N. Li, M. Han, V.G. Romanovski: Cyclicity of some Liénard Systems, Commun. Pure Appl. Anal. Vol. 14, p.p. 2127-2150, 2015 [5] F. Dumortier, J. Llibre, J. C. Artes: Qralieativc Thccey cf Planae Diffcrcntial Systcme, Springer-Verlag, Berlin Heidelberg, 2006 [6] M. Jesús Álvares, A. Ferragut, X. Jarque: A survey on the blow up technique, Internat. J. Bifur. Chaos Appl. Sci. Engrg. Vol. 21, p.p. 3103, 2011 JET 59