Technology Journal of Energy JET Volume 8 (2015) p.p. 43-55 Issue 4, December 2015 Typology of article 1.01 www.fe.um.si/en/jet.html ON THE ISOCHRONICITY OF PERIODIC SOLUTIONS AT A CENTRE MANIFOLD O IZOHRONOSTI PERIODIČNIH REŠITEV NA CENTRALNI MNOGOTEROSTI The problem of isochronicity is discussed from the historical and dynamical systems point of view. The model of Huygen's cycloidal chops is mathematically explained. We consider two dynamical systems arising from a three-dimensional system with a centre manifold. Based on the period function approach we find necessary and sufficient criteria on the coefficients of the system to distinguish between the cases of isochronous and non-isochronous oscillations. Problem izohronosti je obravnavan s stališča dinamičnih sistemov ter iz zgodovinskega stališča. Matematično je razložen Huygensov model ure s cikloidnim nihalom. Na osnovi analize funkcije periode obravnavamo dva dinamična podsistema tridimenzionalnega sistema s centralno mnogoterostjo in poiščemo potrebne in zadostne pogoje za izohronost centra na raznoterosti določeni s koeficienti sistema. R Corresponding author: Matej Mencinger, Tel.: +386 2 229 4321, Mailing address: University of Maribor, Faculty of Civil Engineering, Transportation Engineering and Architecture, Smetanova ulica 17, 2000 Maribor, Slovenia. E-mail address: matej.mencinger@um.si 1 University of Maribor, Faculty of Energy Technology, Hočevarjev trg 1, 8270 Krško, Slovenia 2 Center for applied mathematics and theoretical physics, University of Maribor, Mladinska 3, SI-2000 Maribor, Slovenia 3 University of Maribor, Faculty of Civil Engineering, Transportation Engineering and Architecture, Smetanova ulica 17, 2000 Maribor, Slovenia 4 Institute of Mathematics, Physics and mechanics, Jadranska 19, 1000, Ljubljana, Slovenia Brigita Ferčec12R, Matej Mencinger34 Keywords: centre manifold, centre problem, isochronicity Abstract Povzetek JET 43 arFerta Fercec, Matej Igencinger JET Vol. 8 (2015) Issue 4 1 INTRODUCTION A simple pendulum is a bob of mass m on a massless rigid wire of length l (see Fig. 1). Assuming no damping (i.e. the only force acting on the system is the weight of the bob), the differential equation governing a simple pendulum is d2^ g (1.1) m Figure 1: A simple pendulum of length I driven by the weight F = —m^ of the bob. There are several approximations for a period of a simple pendulum (see Fig.1) of length L One of them (see [1]) is Vsin^o/ ^g (1.2) where ^>o is the initial angle displacement (amplitude) of the bob. Note that substituting x = , y and t = •t into (1.1) one can easily derive the following system dx y3 dy — = —siny« —y + 37—'•■, dT = X. (1.3) The qualitative behaviour of this system is determined by how x(t) and y(x) behave with the change of t. It is convenient to indicate how a solution behaves in the phase plane, that is in the x,y —plane. The qualitative behavior is represented by a family of curves, directed with increasing t. These phase curves are called trajectories or orbits of system (1.1.3). The geometrical representation of the qualitative picture of orbits of system (1.1.3) is called its phase portrait. System (1.1.3) admits a centre at (0,0), that is, all the orbits close to the origin are closed or, in other words, the origin is enclosed with simple closed curves (solutions are periodic). Obviously, 44 JET On the isochTonicity ofperioclic solutions at a centre manifold because of (1.1.2) the centre of (1.1.3) is not isochronous, i.e. not all periodic solutions have the same period. Isochronicity, as will be stated in the next section (see Definition 2), was not known at the time of the great Dutch astronomer, physicist and mathematician Christian Huygens. However, he was aware that if the amplitude of the pendulum's swing changed, the time of swing would also change, which means that the pendulum was not isochronous, [2]. In the limit of small obtained from the equation of a mathematical pendulum (1.1) (again using the change of time However, for larger amplitudes, the period T is a function of the amplitude. The approximation formula (1.1.2) is just one of many existing [1]. Huygens' ingenious idea, which he put into practice, was to vary the effective length of the pendulum by allowing its cord to wrap partially around an obstruction as it swings. What should the shape of this obstruction be in order to ensure that the period is strictly independent of the amplitude? In the middle of the 17th century, the accuracy of clocks was not measured in minutes; even the best of them gained or lost several minutes per day because they were driven by falling weights or springs. Galileo Galilei was probably the first who wanted to design a pendulum driven clock in order to achieve better accuracy. He believed that he had mathematically proven that a (simple) pendulum is isochronous. As already stated, for the case of the pendulum this simply means that the time it takes to complete one full swing is independent of the size of the swing. Nowadays, it is well-known that this is not the case, and that Galileo was mistaken (see approximation formula (1.1.2)). The first scientist who produced a pendulum-driven clock to keep time errors to within one minute per day was Christian Huygens in 1656; within two years the clocks with accuracy of about 10 seconds per day were produced (by Huygens and others). amplitudes, the period (1.1.2) is approximated dx dT (1.4) y (x0 s(t) WO) x Figure 2: A baud ea u wiss ef shops (x(t),y(t)). JET 45 arFerta Fercec, Matej Igencinger JET Vol. 8 (2015) Issue 4 In order to solve the problem of an "isochronous pendulum", a generalization of the simple pendulum in the sense depicted in Fig. 2 was considered. Thus, if we allow the shape connecting (x0,y0) with (0,0) to be different from a semicircle (i.e. a general shape (x(t),y(t))), then Huygens's question (of isochronous swing) becomes for what shape will the bead descend to the origin in the same time regardless of its starting point (x0,y0)? Denoting by s(t) the distance along the string that the bead has traveled at time t and by v(t) = ^ = -J(^) + (y) the velocity of the bead (at time t) which is (in contrast) equal to v(t) = —^2g(y0 — y(t)), where g is the acceleration of a freely falling body due to gravity. Furthermore, assume that the functorial dependence of s(t) in terms of y(t) is s = f(y). Then the time taken by a bead to descend from (x0,y0) to (x(t),y(t)) becomes: , f°dtds f° —1 T(y;y°)=l . f'(y)dy. (15) ■To dy W2g(y° — y) (15) Substituting a dimensionless variable z = — into (1.1.5) and rearranging we obtain yo ^m^r—(1.6) V1 z which is independent of y° (i.e. a constant function of y°) for 0 0. (1.8) y Recall that s = f(y) yielding g=f'(y). Dividing (£) =(§ + by (g yields (2y) (f'(y))2 = (dy) +1 and finally from (1.1.8) we obtain * = /' yo c —y - dy. (1.9) 46 JET On the isochTonicity ofperioclic solutions at a centre manifold The solution to (1.1.9) containing (x,y) = (0,0) can be parametrized (using y = C sin2 J, for 0< 9 < "n:, and the "half angle identities") by CC x(9) =2(9 + sin^), y(9) =2(1_c°s9), (110) which passes through (x0,y0) for m = m0 i tan l2 (satisfying — = 1 cosJ° ) and C= 2y° . I- o v ~r -t-° x° j°+sinj°' 1-cosj The shape of the curve (x(t),y(t)) from (1.10) is called an (inverted) cycloid (see Fig. 3). Definition 1. The cycloid (studied and named by Galileo in 1599) is the locus of a point on the rim of a circle of radius a rolling along a straight line. If the point considered on a circumference of the rolling circle is initially in the origin (0,0), the cycloid has a cusp at the origin and its humps are oriented upward. Its parametric equation is x(t) = a(t- sint), y(t) = a(1 - cost) ............... 2.« y >......... __ Y X Zn-a 0 2na Figure 3: Cycloid: The path followed by a point on the circumference of a circle as that circle rolls along a straight line. An involute is a curve obtained from another given curve by attaching an imaginary taut string to the given curve and tracing its free end as it is (un)wound onto that given curve. The evolute of a curve is the locus of all its centres of curvature. That is to say that when the centre of curvature of each point on a curve is drawn, the resultant shape will be the evolute of that curve. For instance: the evolute of a circle is a single point at its centre, the evolute of a quadratic parabola y = x2 is a semicubical parabola y = i + -V4X2 . The only curve for which the evolute (and involute) is "of the same size" (neglecting translation) as the original curve is the cycloid (see Fig. 3 and Fig. 4). The evolute of an involute is the original curve (fewer portions of zero or undefined curvature). Alternatively, another way to construct the involute of a curve is to replace the taut string by a line segment that is tangent to the curve on one end, while the other end traces out the involute (Fig. 3). JET 47 Brigita Ferčec, Matej Mencinger JET Vol. 8 (2015) Issue 4 Figure 4: Inverted cycloid (black) with its involute (red) which is a (inverted) cycloid with the same parameter a, just translated for a vector (na, —2a). The above results led Huygens to invent a pendulum with cycloidal chops (Fig. 5). Figure 5: The cycloidal chops pendulum invented by C. Huygens. However, it is not possible to determine the centre of oscillation for pendula suspended between cycloids, as seen in Fig. 4, since the motion of a pendulum with cycloidal chops is not planar (which was known by Huygens; such a pendulum is a multiple system, i.e. a system having two degrees of freedom). Therefore, the next step of improvements of the pendulum clocks was the so-called escapement mechanism (which was driven primarily by the gravity force and later by the force of a spring). 2 ISOCHRONICITY OF PLANAR SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS We discuss differential systems with analytical right sides TO TO dx Y"1 ■ ■ dy Y"1 ■ ■ dt = —y+L pijxly],dt=x + Z qijXlyJ, (2.1) i+j = 2 i+j=2 where pij and qij are real constants. There has been a longstanding problem, [3], called the Poincare centre-focus problem; for the system (2.1) to find explicit conditions under which (2.1) has a centre at the origin. The problem is equivalent to an analogue for a corresponding periodic equation 48 JET On the isochTonicity ofperioclic solutions at a centre manifold dr Zr=iAj(9)r d9 l+Zr=iBi(q (2.2) to have periodic solutions. To see this let us note that the phase curves of (2.1) near the origin in polar coordinates x = r cos 9, y=r sin 9 are determined by (2.2), where A;(9) and Bj(9) are polynomials in cos 9 and sin 9 and Rj(9) are 2n —periodic functions of 9 and the series on the right-hand side of (2.2) is convergent for all 9 and all sufficiently small r, [4]. Since the closed orbits of (2.1) correspond to periodic solutions of (2.2), the planar analytic vector field (2.1) has a centre at (0,0) if and only if (2.2) has a centre at r = 0 ; that is, all the solutions nearby are periodic. All methods for determining the nature of the critical point of (2.1) theoretically require determining infinitely many coefficients of some function or as we will see later the corresponding ideal of the infinitely many polynomials, [4,5,6,7]. The most widely used method to resolve the Poincare centre-focus problem are the so-called Lyapunov method (combined with the Bautin method) and Mironenko's method, [8]. Solving the Poincare centre-focus problem for a given system (2.1) may be a tedious job. However, once we know the origin of (2.1) is a centre, the problem of isochronicity becomes sensible, and the definition of the period function makes sense. Definition 2 ([4]). Suppose r* >0 is so small that it lies in the so-called period annulus (the maximal neighbourhood of the origin that is enclosed with simple closed curves). Let us consider the line segment X = {(x,y): 0 1 such that ideal I with infinitely many polynomials Ti(pij, qij) is the same as Im. 3 CENTRE MANIFOLDS The theory of centre manifolds, which originates from the works of V. A. Pliss, [10,11], and was then further developed by many others (see, e.g. [12,13] and references therein) is an extremely effective tool for studying the behaviour of trajectories of high-dimensional systems of ordinary differential equations. Note that one of the important applications of centre manifolds (also called the Pliss reduction principle, [11]) is that it allows the reduction of the system and consequently the study of the stability in lower dimensional phase space (thus, instead of analysing the original high-dimensional phase space we can consider the lower-dimensional subsystem). Consider an m + n-dimensional system of ordinary differential equations of the form x = Ax + u(x,y) y = By + v(x,y), (3.1) where x e!m,y e!n, Re(a(A)) = 0, Re(a(B))^0, a(A) and a(A) are spectrums of A and B, respectively, and u, v are Ck- functions, k>1 which vanish together with their first derivatives 50 JET On the isochTonicity ofperioclic solutions at a centre manifold at the origin. By definitio, a Ck- manifold Wc = WAc(0, U) in a neighbourhood U of 0 is said to be a centre manifold of (3.1) if Wc is invariant under the flow as long as the solution remains in U nd Wc is the graph of a Ck- function y = h(x) which is tangent at (0,0) £ Rmxln to the x-space. The following fundamental result shows that for the system (3.1), there is always a centre manifold in a neighbourhood of the origin. Theorem 1 ([12]). There exists a neighbourhood U of (0,0) £ lm X ln such that there exists a local centre manifold Wc of (3.1) which is the graph of u Ck- function y = h(x). The simplest examples of the centre manifolds can be seen in low dimensions. For instance, consider two-dimensional system x = —x3, y =—y. (3.2) For any constants c1,c2 the curve h(x; c1,c2), where i—^e"^, x <0 0, x = 0 c2e~2x2, x > 0, consists of orbits of (3.2) and is therefore invariant. Furthermore, this curve is tangent to the x-axis at the origin, which implies that it is a centre manifold. From this example, we can also see that the centre manifold is not necessaily unique. Another simple two-dimensional system with a centre manifold is x = —x3, y = —y + x2. Similar as in the previous case, we see that although the system is analytic, it is not difficult to verify that the centre manifold is not analytic. In three dimension, a family of systems with a centre manifold is u = —v +P(u,v,w) v = u + Q(u,v,w) (3.3) w = —A.w + R(u,v, w), where A is a positive real number and P, Q, R are polynomials without constant and linear terms. By Theorem 1, this system has a centre manifold w = f(u,v). There are many systems arising from physics (for instance, the Rikitake system [14], for Earth's magnetic field or the Hide-Acheson Dynamo, [15]) that possess a fixed point at which the linear part has one negative and two purely imaginary eigenvalues and are therefore of the form (3.3). Since system (3.3) has a centre manifold Wc and X>0, the trajectories in a small neighbourhood of the origin tend to the trajectories to the centre manifold as time increases. In systems (3.3), the phase portrait in a neighbourhood of the origin on can be, depending on the added nonlinear terms P and Q, either a centre in which case every trajectory (other than the origin itself) is an oval surrounding the origin, or a focus, in which case every trajectory spirals towards the origin or every trajectory spirals away from the origin as the time increases. The problem of determining the dynamical JET 51 arFerta Fercec, Matej Igencinger JET Vol. 8 (2015) Issue 4 behaviour on Wc, that is, distinguishing between a centre and a focus on the centre manifold for a quadratic polynomial system of the form (3.3) was studied in [16]. 4 THE STUDY OF TWO EXAMPLES In [16], the authors studied the dynamics of trajectories at the centre manifolds for the system u = —v + au2 + av2 + cuw + dvw v = u + bu2 + bv2 + euw + fvw (4.1) w = —w + Su2 + Sv2 + Tuw + Uvw, where the coefficients a, b, c, d, e, f, S, T, U are real. They found five conditions for the existence of a centre on the centre manifold: 1. S = 0; 2. a = b = c + f=8c + T2 — U2=4(e — d) —T2U2 = 2(e + d) + TU and S = 1; 3. a = b = c = f= d + e = 0 and S=1; 4. d + e = c = f= T — 2a = U — 2b = 0 and S = 1; 5. c = d = e = f =0 and S=1. The next question that naturally arises is whether the centre at the centre manifold is isochronous, that is, whether all oscillations have the same period. In this paper we study two subfamilies of system (4.1) satisfying conditions 1 and 2 above. For each system, we compute the period function as described in section 2 and use it to find conditions for the centre to be isochronous. System (4.1) satisfying condition 1 above is written as ù = —v + au2 + av2 + cuw + dvw v =u + bu2 + bv2 + euw + fvw (4.2) w = —w + Tuw + Uvw. It is proven (see [16]) that w = 0 isthe centre manifold for system (4.2). The corresponding 2D system u = —v + au2 + av2 v =u + bu2 +bv2 has a centre at the origin for all a,b£l (see [16]). We now study the isochronicity problem for the above centre. Using the computer algebra software MATHEMATIC, we first turn to the computation of the period function T of the form (2.6). After introducing polar coordinates system, (4.3) becomes 52 JET On the isochTonicity ofperioclic solutions at a centre manifold dr — = r2(acos9 + bsin9) d9 — =1 + r(b cos 9 — asin 9). Following the procedure described in the second section, we find that the first non-zero coefficient of the period function is T2 = 2u(a2 + b2). Thus we see that the necessary condition for the isochronicity of system (4.3) is a = b = 0, which, obviously is also the sufficient condition. In the case of conditions 2, under the same renaming of parameters c,d, e,f,T, U and using a,p as in [16], the system (4.1) takes the form 1 1 2 ù = —v--aBuw--B2vw _ 2 H r (4.5) v = u + 2 a2 uw + ^ apvw w = —w + u2 + v2 + (a + p)uw + (p — a)vw. A search for invariant algebraic surfaces led to the explicit equation of centre manifold Wc given u2+v2 by w = 1 au pv. Inserting the expression for w into system (4.5), we obtain the system ^(au + ^v)(u2 + v2) u = —v — ■ V = u + 2(1 — au — yffv) a(au + ^v)(u2 + v2) 2(1 — au — ^v) Using Taylor series expansion up to the order five, we obtain the system u = —v + ip(au + pv)(u2 + v2) +^p(au + pv)2(u2 + v2) + m1 2 2 (4.6) 1 , , , - „1 where v = u — 2a(au + pv)(u2 + v2) — 2a(au + pv)2(u2 + v2) + m2, m1 = 120(—60a3Pu5 — 180a2p2u4v+10u3v2(—6a3p — 18ap3) + 10u2v3(—18a2p2 — 6p4) — 180ap3uv4 — 60p4vs) and m2 = 120 (60a4 us + 180a3pu4v+10u3v2(6a4 + 18a2p2) +10u2v3(18a3p + 6ap3) + 180a2p2uv4 + 60apV). Further computation following the computational pattern described in section 2 yields T2 = J (a2 + p2), yielding the following result. System (4.6) has an isochronous centre if and only if a = p = 0. JET 53 Brigita Ferčec, Matej Mencinger JET Vol. 8 (2015) Issue 4 5 CONCLUSIONS After more than 300 years, the isochronicity problem remains one of the central problems in the theory of dynamical systems. In addition to the period function approach, there are some other be too complex for some systems. Therefore, in some cases, we need to find other ways to solve the isochronicity problem. It turns out that the isochronicity problem for system (2.1) is equivalent to the linearizability problem in which we look for an analytic change of coordinates that reduces (2.1) to the canonical linear centre X = —y,y = x (see e.g. [4] for more details). The theory of the linearizability of planar systems of ordinary differential equations was applied [17] to system (4.1) corresponding to condition 4 above. In fact, for each of the five conditions above we can find the conditions under which the corresponding system can be reduced to canonical linear centre, and we see that they are the same as conditions obtained using the period function approach. References [1] G.E. Hite: Approximations for the Period of a Simple Pendulum, The Physics Teacher, Vol. 43, p.p. 290, 2005 [2] C. Huygens The pendulum clock or geometrical demonstrations concerning the motion of pendula as applied to clocks. Translated from the Latin and with a preface and notes by Richard J. Blackwell. With an introduction by H. J. M. Bos. Iowa State University Press, Ames, IA, 1986 [3] H. Poincaré: Mémoire sur les courbes définies par une équation différentielle, J. Math. Pures et. Appl. (Sér. 3), Vol. 7, p.p. 375-422, 1881; (Sér. 3), Vol. 8, p.p. 251-296, 1882; (Sér. 4), Voil. 1, p.p. 167-244, 1885; (Sér. 4), Vol. 2, p.p. 151-217, 1886 [4] V.G. Romanovski, D.S. Shafer: The Center and cyclicity Problems: A computational Algebra Approach, Boston: Birkhäuser, 2009 [5] B. Ferčec, X. Chen, V. Romanovski: Integrability conditions for complex systems with homogeneous quintic nonlinearities, Journal of applied analysis and computation, Vol. 1, Iss. 1, p.p. 9 - 20, 2011 [6] B. Ferčec, J. Giné, Y. Liu, V. Romanovski: Integrability conditions for Lotka-Volterra planar complex quartic systems having homogeneous nonlinearities, Acta applicandae mathematicae, Vol. 124, Iss. 1, p.p. 107 - 122, 2013 [7] B. Ferčec, J. Giné, M. Mencinger, R. Oliveira: The center problem for a 1 : -4 resonant quadratic system, Journal of mathematical analysis and applications, Vol. 420, Iss. 2, p.p. 1568 - 1591, 2014 54 JET On the isochTonicity ofperioclic solutions at a centre manifold [8] Z. Zhou: A Nsw Msihed fes Rstsuseh ea ihs Csaiss-Feect Pseblsm ef Differential Systems, Abstract and Applied Analysis, Art. ID 926538, 5 p.p., 2014 [9] D. Cox, J. Little, D. O'Shea: imsult, Vusisiist, uad Algesiihmt: Aa laisemceiiea if Computational Algsbsuie Gsemsiso uad Cemmciuiivs Algsbsu, New York: Springer, 3rd edition, 2007 [10] S.Yu. Pilyugian, G.R. Sell: A Biegsuphieul Sksieh ef Vieies A. Plitt, Journal of Dynamics and Differential Equations, Vol.15, 2003 [11] V.A. Pliss: A ssmceiiea psiaeipls ia ihs ihsesy ef tiubiliio ef meiiea, Izv. Akad. Nauk SSSR Ser. Mat., Vol. 28, p.p. 1297-1324, 1964 [12] S. Chow, J. K. Hale: Msihedt ef BiCcshuiiea Thsesy, Springer-Verlag, New York, 1982 [13] J. Sijbrand: Psepsnist ef esaissmuaifelmt, Trans.~Amer.~Soc., Vol. 289, 431-469, 1985 [14] R. Rikitake: Oteilluiieat ef u tytism ef mite myaumet, Proc. Cambridge Philos. Soc., Vol. 54, p.p. 89-105, 1958 [15] R. Hide, A. C. Skeldon, D.J. Acheson: A ticdy ef iwe aevsl tslC-sxhiiiag tiagls-ditk hemepelus myaumet: ihsesy, Proc. R. Soc. Lond. A , Vol. 452, p.p. 1369-1395, 1996 [16] V.F. Edneral, A. Mahdi, V.G. Romanovski, D.S. Shafer: Ths esaiss pseblsm ea u esaiss muaifelm ia R3, Nonlinear Anal., Vol. 75, p.p. 2614-2622, 2012 [17] V. Romanovski, M. Mencinger, B. Fercec: iavstiiguiiea ef esaiss muaifelmt ef ihsss-mimsatieaul tytismt ctiag eempcissulgsbsu, Programming and computer software, Vol. 39, iss. 2, p.p. 67-73, 2013 JET 55