This is a preprint of: “Rocard’s 1941 chaotic relaxation econometric oscillator”, Jean-Marc Ginoux, Franck Jovanovic, Riccardo Meucci, Jaume Llibre, Internat. J. Bifur. Chaos Appl. Sci. Engrg., vol. 32(3), 2250043 (12 pages), 2022. DOI: [10.1142/S0218127422500432] begindocument/before ROCARD’S 1941 CHAOTIC RELAXATION ECONOMETRIC OSCILLATOR Jean-Marc GINOUXa∗ , Franck JOVANOVICb , Riccardo MEUCCIc & Jaume LLIBREd a Aix Marseille Univ, Université de Toulon, CNRS, CPT, Marseille, France, ginoux@univ-tln.fr, b School of Business Administration, Teluq University, Canada LEO, Université d’Orléans, France, c Istituto Nazionale di Ottica, Consiglio Nazionale delle Ricerche, Largo E. Fermi 6, 50125 Firenze, Italy, d Departament de Matemàtiques, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Spain, In the beginning of the Second World War, the French physicist, Yves Rocard, published a book entitled Théorie des Oscillateurs (Theory of Oscillators). In Chapter V, he designed a mathemat- ical model consisting of a set of three nonlinear differential equations and allowing to account for economic crises. Numerical integration of his model has highlighted a chaotic attractor. Its analysis with classical tools such as bifurcation diagram and Lyapunov Characteristic Exponents has con- firmed the chaotic features of its solution. It follows that Rocard’s 1941 chaotic econometric model has thus most likely preceded Lorenz’ butterfly of twenty-two years. Moreover, apart this historical discovery which upsets historiography, it is also established that this “new old” three-dimensional autonomous dynamical system is a new jerk system whose solution exhibits a chaotic attractor the topology of which varies, from a double scroll attractor to a Möbius-strip and then to a toroidal attractor, according to the values of a control parameter. I. INTRODUCTION According to the historiography, it is generally considered that the very first chaotic attractor has been designed in 1963 by the late Edward Norton Lorenz [15]. But the winding road taken by the theory of nonlinear oscillations sometimes leads to surprises. Long and deep investigations performed in this domain [4] have led us to the “discovery” of a book entitled Théorie des Oscillateurs (Theory of Oscillators) and published by the French physicist Yves Rocard (1903-1992) in 1941. In chapter V: “Les oscillateurs des théories économiques” (Oscillators of economics theories), Rocard designed two mathematical models allowing to account for economic crises. In the first one, presented in subsection B, Rocard proves that economic crises can be modeled by using a Van der Pol’s relaxation oscillator [22] and obtained, to our knowledge, probably the very first relaxation econometric oscillator. Then, by considering that the frequency of oscillations may strongly depends on the amplitude, Rocard modified his first model and obtained a second one which is a chaotic relaxation econometric oscillator presented in subsection C. This latter model, which is the subject of this present work, will be analyzed in Sec. II and III. Then, by using classical analysis tools such as bifurcation diagram and Lyapunov Characteristic Exponents, it will be thus established that this model is a new jerk system the solution of which exhibits not only one but several chaotic attractors according to the value of the bifurcation parameter. A. Rocard’s econometric oscillators Starting from the analogy with a nonholonom oscillator, Rocard [17, p. 126] imagines an econometric model that he describes as follows: “Suppose that y is the price of a commodity, that y1 is the number of consumers of that commodity, or its total consumption, and assume that y2 is the degree of tooling or mechanization, or rationalization, involved in the manufacture of this commodity and tending to lower its price. We will reason less about the quantities themselves than about their differences from a position of equilibrium that will not be quantified.” 2 Then, he obtained the following three-dimensional dynamical model consisting in linear ordinary differential equa- tions:   dy1   = −ay1 + by,   dt  dy2 = K (y + y1 ) , (1)   dt    dy  m = −y2 . dt where a is a positive parameter while b is negative and m and K are two unknown coefficients to be determined. In this model: −a (a > 0) represents the decreasing rate of the number of consumers or decreasing rate of the total consumption, b (b < 0) is the decreasing rate of a commodity price, K corresponds to the commodity price growth rate and to the total consumption growth rate, m is the mass of investments. By taking the second time derivative of the last equation of (1) and by making a linear combination of the two others, he obtains the following third-order linear ordinary differential equation: ... m y + amÿ + K ẏ + K (a + b) y = 0. (2) B. Relaxation econometric oscillator Then, Rocard explains that as long as b is negative, the dynamical system (1) or its corresponding linear differential equation (2) cannot exhibit any self-oscillations, i.e., self sustained oscillations. As a consequence, in order to obtain such kind of oscillations, he suggests to replace b by b[1 − y 2 /y02 ] where y0 is a constant in the above equations (1-2). Let’s notice that Rocard introduces this nonlinear oscillations characteristic “by hand” in his model, i.e., without any economical justification. He obtains the following three-dimensional dynamical consisting in nonlinear ordinary differential equations:      dy1 y2   = −ay 1 + b 1 − ,   dt y02  dy2 (3)  = K (y + y1 ) ,   dt   m dy = −y2 ,  dt that he transforms into:    ... y2 m y + amÿ + K ẏ + K a + b 1 − 2 y = 0. (4) y0 According to Sprott [20], the third-order nonlinear ordinary differential equation (4) is a jerk equation and the dynamical system (3) is jerk system. Then, Rocard [17, p. 128] explains that: “The equations of the system (3) are no more linear, and their mathematical analysis becomes more difficult. However, we have the study of relaxation oscillations to guide us, and we will quickly see that we can conclude to the existence of self-sustaining oscillations of finite amplitude.” Thus, Rocard [17, p. 130] performs a classical analysis of his third-order nonlinear ordinary differential equation, i.e. his jerk equation (4) and plot its solution (see Fig. 1). 3 FIG. 1: Solution of Rocard’s jerk equation (4). Then, Rocard [17, p. 130] explains that the curves represented on Fig. 1 are “very similar to those of relaxation oscillations”. From Fig. 1, he deduces that “for the variation of prices y over time we obtain a fairly characteristic law of slow rise when prices are low, accelerated when they are high, then slowly fall, becoming a little faster when they are. lower, etc . . . ”. Finally, Rocard [17, p. 131] states mathematically that the frequency decreases as the amplitude increases and he considers that it would be interesting to analyze the case of an “oscillator whose frequency depends much on the amplitude”. By posing: K ω2 = ; a = εω ; b = ηω ; y = y0 z, m He obtains the following dimensionless third-order nonlinear ordinary differential equation: ...   z + εωz̈ + ω 2 ż + ω 3 ε + η 1 − z 2 z = 0. (5) Although Le Corbeiller [13] and Hamburger [8, 9] had suggested to apply Van der Pol’s relaxation oscillations to Econometry, this is only in 1951 that Goodwin [6] proposed a prototype nonlinear differential equation exhibiting maintained or self-sustained oscillations including relaxation oscillation. So, it appears that Rocard’s jerk equation (5) which has preceded that of Goodwin of ten years can be considered as the paradigm of relaxation oscillations in Econometry and also upsets the historiography. C. Chaotic relaxation econometric oscillator In the second section of chapter V, Rocard [17, p. 133] designed a second model which is an “oscillator whose frequency depends much on the amplitude”. To this aim, he modified the nonlinear oscillations characteristic, i.e., the last term of Eq. (5), by replacing (1 − z 2 ) with (1 − z 2 − ż 2 /ω 2 ). Thus, he obtains the following dimensionless third-order nonlinear ordinary differential equation:    ... ż 2 z + εωz̈ + ω 2 ż + ω 3 ε + η 1 − z 2 − 2 z = 0. (6) ω Then, Rocard [17, p. 133] explains that: “It would be interesting to provide a case study for which the frequency variation according to the ampli- tude can even be totally abnormal.” This last model, is, to our knowledge the first chaotic relaxation econometric oscillator. It will be analyzed in the next section. We will show that what he considered as “abnormal” is in fact the expression of the chaotic behavior of the solution of his jerk equation (6). 4 II. ROCARD’S 1941 CHAOTIC RELAXATION OSCILLATOR First, let’s notice that, according to D’Alembert [1], the third-order nonlinear ordinary differential equation (6) can be cast in the form of a system of coupled first-order nonlinear differential equations as follows:   dx   = −ω [εx + ωy + ωz] ,   dt   dy    x2 = ω ε + η 1 − z2 − 2 z, (7)   dt ω     dz = x.  dt Although in the second section of his chapter V, Rocard [17, p. 133] did not assign any value to the parameter set of his model, we have performed many preliminary tests to determine the parameters range within which chaotic attractors may appear. This led us to use the following parameter set: ε = 0.5, ω = 2 and η ∈ [−1.34, 0.94] which will be used in the next sections. These parameter values correspond in the non-dimensionless form to a = εω = 1, b ∈ [−2.68, 1, 88] since b = ηω and K = ω 2 m = 4m. Let’s notice that they do have a real economic significance, according to REF since... III. STABILITY ANALYSIS A. Equilibrium points By using the classical nullclines method, we found that Rocard’s system (7) has the following three equilibrium points:  r r   r r  ε+η ε+η ε+η ε+η O (0, 0, 0) ; I1 0, ,− ; I2 0, − , (8) η η η η B. Jacobian matrix The Jacobian matrix of Rocard’s dynamical system (7) reads:   −εω −ω 2 −ω 2    2η ηx2  J = − ω xz 0 − + ω ε + η − 3ηz 2   (9)  ω  1 0 0 By replacing the coordinate of the equilibrium point O (8) in the Jacobian matrix (9) one obtains the following Cayley-Hamilton third degree eigenpolynomial: λ3 + εωλ2 + ω 2 λ − ω 3 (ε + η) = 0. (10) By using the Routh-Hurwitz criterion [10, 18] to state the stability of O, we obtain the following three determinants:  ∆1 = εω,  ∆2 = −ω 3 η,   ∆3 = −ω 6 (η + ε) η. Since with our parameter set, ε = 0.5 and ω = 2, it follows that all these three determinants are strictly positive provided that: −ε < η < 0. (11) 5 Thus, all the real parts of the eigenvalues of the eigenpolynomial (10) are negative and so, O is a stable equilibrium point provided that condition (11) is verified. Now, by replacing the coordinate of the equilibrium points I1 or I2 (8) in the Jacobian matrix (9) one obtains the following Cayley-Hamilton third degree eigenpolynomial: λ3 + εωλ2 + ω 2 λ + ω 3 (ε + η) = 0. (12) Still using the Routh-Hurwitz criterion [10, 18] to state the stability of I1,2 , we obtain the three determinants:  ∆1 = εω,  ∆2 = ω 3 (3ε + 2η) ,   ∆3 = −2ω 6 (3ε + 2η) (ε + η) . Since with our parameter set, ε = 0.5 and ω = 2, it follows that all these three determinants are strictly positive provided that: 3ε − < η < −ε. (13) 2 Thus, all the real parts of the eigenvalues of the eigenpolynomial (12) are negative and so, I1,2 are stable equilibrium points provided that condition (13) is verified. Then, by taking into account both conditions (11) and (13), we find that if: 3ε − <η<0 2 at least one of the three equilbrium points is stable, but this does not imply that other attractors can coexist with this local stable equilibrium, see for instance [21]. The characteristic polynomial (10) becomes the polynomial (λ − κ)(λ − σi)(λ + σi) with κσ ̸= 0, and consequently the equilibrium point at the origin of coordinates O has the possibility of exhibiting a Hopf bifurcation because then the eigenvalues of its linear part are κ and ±σi for ω > 0 if and only if ε = −κ/σ, η = 0 and ω = σ. We can expect to see a small-amplitude limit cycle biburcating from the equilibrium point O. In order to confirm that a Hopf bifurcation appears at O we must compute the first Lyapunov coefficient ℓ1 (O) of the differential system at O. When ℓ1 (O) ̸= 0 the equilibrium point O is a weak focus of the differential system (7) restricted to the central manifold of O and the limit cycle that emerges from O is stable if ℓ1 (O) < 0, and unstable if ℓ1 (O) > 0. In the first case we say that the Hopf bifurcation is supercritical, and in the second case we say that the Hopf bifurcation is subcritical. Here we use the following result presented in the pages 175–180 of the book [12] for computing ℓ1 (O). Lemma 1. Let ẋ = F (x) be a differential system with x ∈ R3 having O as an equilibrium point. Consider the third 1 1 order Taylor approximation of F around O given by F (x) = Ax + B(x, x) + C(x, x, x) + O(|x|4 ), where A is a 2! 3! matrix, B is a bilinear function and C is a trilinear one. Assume that the matrix A has a pair of purely imaginary eigenvalues ±σi. Let q be the eigenvector of A corresponding to the eigenvalue σi, normalized so that q · q = 1, where q is the conjugate vector of q. Let p be the adjoint eigenvector such that AT p = −σip and p · q = 1. If I denotes the 3 × 3 identity matrix, then 1 ℓ1 (O) = Re(p · C(q, q, q) − 2p · B(q, A−1 B(q, q)) + p · B(q, (2σiI − A)−1 B(q, q))). 2σ 6 Some easy but tedious computations show that   κ −σ 2 −σ 2 A=0 0 −κ  , 1 0 0   iσ iκ 1 q = q , q ,q , κ2 2 κ2 2 κ2 2 σ2 + σ + 1 σ σ2 + σ + 1 σ2 + σ + 1   1 iσ 1 p = − q , q ,q , 2 2 κ2 (κ + iσ) σκ2 + σ 2 + 1 (κ + iσ) σκ2 + σ 2 + 1 σ2 + σ 2+1 B((x1 , y1 , z1 ), (x2 , y2 , z2 )) = (0, 0, 0),   6η(x1 x2 z3 + x1 z2 x3 + z1 x2 x3 ) C((x1 , y1 , z1 ), (x2 , y2 , z2 ), (x3 , y3 , z3 )) = 0, − − 6ηωz1 z2 z3 , 0 , ω  3ησ 3 σ 2 + ω 2 ℓ1 (O) = − . 2ω (κ2 + σ 2 ) (κ2 + σ 4 + σ 2 ) In summary, by Lemma 1 it follows that when ℓ1 (O) ̸= 0 the differential system (7) for ε = −κ/σ, η = 0 and ω = σ > 0 exhibits a Hopf bifurcation at the equilibrium point O. C. Bifurcation diagram According to Rocard [17, p. 133], the amplitude of his jerk equation (6) and so, of the jerk system (7) much depends on the parameter η. The same is true for existence of chaotic attractors for the jerk system (7). Thus, in order to highlight how the changes of this control parameter impact the corresponding topology of the attractor, we have built a bifurcation diagram for η ∈ [−1.34, −0.75] (see Fig. 2) and for η ∈ [0, 0.94] (see Fig. 3) since in the interval η ∈ [−0.75, 0] one at least of the three equilibrium points is stable. On Fig. 2, we observe a reverse period doubling cascade which confirms the existence of chaotic attractors for −1.34 ⩽ η ⩽ −0.75. As parameter η increases from −1.34 to −0.75, the chaotic attractor becomes a limit cycle. Let’s notice on Fig. 2 the presence of several “windows” within which the chaotic attractor becomes a limit cycle whose period is determined by the number of branches. As an example, for η ∈ [−1.319, −1.308], the attractor becomes a limit cycle of period 9. On Fig. 3, bifurcation diagram highlights a period doubling cascade for 0 ⩽ η ⩽ 0.94. Starting from η = 0 to η ≈ 0.72 the attractor is a limit cycle and then becomes chaotic. There are also several windows within which the attractor topology changes and we observe limit cycles whose period is given by the number of branches. As an example, for η ∈ [0.87, 0.92], the attractor becomes a limit cycle of period 5. The attractor topology changes according to the control parameter values of η may be represented as follows: −1.34 −−−−−−−−−−→ −0.75 −−−−−−−−−−→ 0 −−−−−−−−−−→ 0.65 −−−−−−−−−−→ 0.94 Chaos / LCn HB SFP HB LC1 Chaos / LCn where LCn means limit cycle of period n, SFP, Stable equilibrium Points and HB, Hopf Bifurcation. For −1.34 ⩽ η ⩽ −0.75, a reverse period doubling cascade occurs (see Fig. 3) till the control parameter η reaches the vale of the Hopf bifurcation ηHopf = −0.75. Thus, we observe for −1.34 ⩽ η ⩽ −1.15 a chaotic double-scroll (see Fig. 4a & 4c). Within this interval, several windows appear on the bifurcation diagram (see Fig. 2) and correspond to limit cycles of period n. As another example, for η = −1.25 the attractor becomes a limit cycle of period 5 (see Fig. 4b). Then, starting for −1.05 ⩽ η ⩽ −0.995, the topology of the attractor changes and it becomes a Möbius-strip (see Fig. 4d-4f). For η ⩾ −0.75, the attractor is a limit cycle of period 1. For −0.75 ⩽ η ⩽ 0, one of the three equilibrium points is stable. As highlighted on the bifurcation diagram (see Fig. 3), from η = 0 to η ≈ 0.72 the attractor is a limit cycle. Then, a period doubling cascade occurs. From η ≈ 0.77, we observe a toroidal chaotic attractor (see Fig. 5a, 5b, 5d & 5e). Again, there are several window within which some limit cycles of period n appear. As an example, Fig. 5c highlights a limit cycles of period 7 for η = 0.825 while a limit cycles of period 5 for η = 0.9 (see Fig. 5e). 7 FIG. 2: Bifurcation diagram xmax as function of η for η < 0. FIG. 3: Bifurcation diagram xmax as function of η for η > 0. In order to confirm the topology of these attractors, Lyapunov exponents have been computed in each case. D. Numerical computation of the Lyapunov exponents ® The algorithm developed by Sandri [19] for Mathematica has been used to perform the numerical calculation of the Lyapunov characteristics exponents (LCE) of dynamical system (7) in each case. LCEs values have been computed within each considered interval (η ∈ [−1.34, −0.75] and [0, 0.94]). Then, following the works of Klein and Baier [11], a classification of (autonomous) continuous-time attractors of dynamical system (12) on the basis of their Lyapunov ® spectrum, together with their Hausdorff dimension is presented in Tab. 1. LCEs values have been also computed with the Lyapunov Exponents Toolbox (LET) developed by Siu for MatLab and involving the two algorithms proposed by Wolf et al. [24] and Eckmann and Ruelle [3] (see https://fr.mathworks.com/matlabcentral/fileexchange/233-let). Results obtained by both algorithms are consistent. We observe on figures 4 that the topology of the attractor of Rocard’s chaotic relaxation econometric oscillator (7) varies. For −1.34 ⩽ η ⩽ −1.15 the attractor is a double scroll (see Fig. 4a & 4c) which may become for particular values of η, a limit cycle (see Fig. 4b), the period of which is given by the number of “branches” observed in the 8 (a) η = −1.34 (b) η = −1.25 (c) η = −1.15 (d) η = −1.05 (e) η = −1 (f) η = −0.995 FIG. 4: Rocard’s chaotic relaxation econometric oscillator (7) in the phase space for various values of η < 0. 9 (a) η = 0.77 (b) η = 0.8 (c) η = 0.825 (d) η = 0.865 (e) η = 0.9 (f) η = 0.925 FIG. 5: Rocard’s chaotic relaxation econometric oscillator (7) in the phase space for various values of η > 0. 10 TABLE I: Lyapunov characteristics exponents of Rocard’s dynamical system (7) for various values of η. η LCE spectrum Dynamics of the attractor Hausdorff dimension −1.34 ⩽ η ⩽ −0.995 (+, 0, −) Chaos 2.04 ⩽ D ⩽ 2.22 −0.95 ⩽ η ⩽ −0.75 (0, −, −) Limit cycle of period 1 D=1 0 ⩽ η ⩽ 0.72 (0, −, −) Limit cycle of period 1 D=1 0.72 ⩽ η ⩽ 0.94 (+, 0, −) Chaos 2.04 ⩽ D ⩽ 2.17 bifurcation diagram (see Fig. 2). Then, for η ≈ −1.09, one the two scrolls of the attractor disappears giving rise to the Möbius-strip (see Figs. 4d, 4e & 4f). This latter disappears on its turn to become a limit cycle according to the reverse period doubling cascade scenario presented on the bifurcation diagram (see Fig. 2). For 0 ⩽ η ⩽ 0.94 the attractor slips from a limit cycle to a Möbius-strip (see Fig. 5a & 5b) via a period doubling cascade route to chaos as highlighted in the bifurcation diagram (see Fig. 3). Moreover, it may become for particular values of η, a limit cycle (see Fig. 5c & 5e), the period of which is given by the number of “branches” observed in the bifurcation diagram (see Fig. 3). Then, for η ≈ 0.825, the attractor becomes toroidal (see Fig. 5d & 5f). IV. CONCLUSIONS Deep investigations of applications of‘nonlinear oscillations theory in the domain of Econometric induced by one of us (F.J.) have led us to the “discovery” of a book entitled Théorie des Oscillateurs (Theory of Oscillators) and published by the French physicist Yves Rocard (1903-1992) in 1941. In chapter V: “Les oscillateurs des théories économiques” (Oscillators of economics theories), Rocard designed two mathematical models allowing to account for economic crises. Each of these two models consists of a set of three nonlinear differential equations which allows to account for economic crises and which can be transformed into a third-order nonlinear ordinary differential equation, that is to say into a jerk equation according to Sprott [20]. With the former model (3), Rocard proved that economic crises can be modeled by using a Van der Pol’s relaxation oscillator [22] and obtained, to our knowledge, probably the very first relaxation econometric oscillator (5). Till recently, the historiography [23] considered that this is only in 1951 that the American mathematician and economist, Richard Goodwin [6], proposed a prototype nonlinear differential equation exhibiting maintained or self-sustained oscillations including relaxation oscillation. Then, Rocard modified his first equation (5) in order to have an “oscillator whose frequency depends much on the amplitude”. Thus, he obtained the third-order nonlinear ordinary differential equation or jerk equation (6) the investigations of which led him to the conclusion that the “frequency variation according to the amplitude can even be totally abnormal”. Nevertheless, he neither further analyzed this jerk equation (6) nor assign any value to the parameter set. So, by considering a realistic range of parameter set, from the econometric point of view, we performed many preliminary tests and determined that for ε = 0.5, ω = 2 and η ∈ [−1.34, 0.94] several chaotic attractors may appear. Then, we have transformed the third-order nonlinear ordinary differential equation or jerk equation (6) into a dynamical system that we have analyzed by using classical tools such as, equilibrium points stability, occurrence of Hopf bifurcations, bifurcation diagram and Lyapunov Characteristic Exponents. Such mathematical and numerical analysis has enabled to confirm that the solution of this “new old” three-dimensional autonomous dynamical system or new jerk system (7) exhibits a chaotic attractor the topology of which varies, from a double scroll attractor to a Möbius-strip and then to a toroidal attractor, according to the values of a control parameter η via a reverse period doubling cascade and period doubling cascade. Thus, it appears that Rocard [17] has stated in 1941, twenty-two years before Edward Norton Lorenz [15], the very first chaotic attractor. This result upsets the historiography [2, 5, 7, 14] who considered till now that Lorenz [15] had been the first to propose a nonlinear dynamical system the solution of which was exhibiting the famous butterfly. So, in this work, we have shown that contrary to what one thought, the very first chaotic attractor has not been designed for modeling atmospheric convection in the domain of Meteorology but for modeling great amplitude variations of relaxation oscillations in Econometry. Acknowledgments ??????????????? 11 The fourth author is partially supported by the Agencia Estatal de Investigación grant PID2019-104658GB-I00, and the H2020 European Research Council grant MSCA-RISE-2017-777911. [1] D’Alembert, J., “Suite des recherches sur le calcul intégral, quatrième partie : Méthodes pour intégrer quelques équations différentielles,” Hist. Acad. Berlin, tome IV (1748) pp. 275–291. [2] Bergé, P., Pomeau, Y. & Vidal, C., Order within Chaos: Towards a Deterministic Approach to Turbulence, New York: John Wiley & Sons, 1984. [3] Eckmann, J.P. & Ruelle, D., “Ergodic Theory of Chaos and Strange Attractors,” Rev. Mod. Phys., 57 (1985) pp. 617-656. [4] Ginoux, J. M., History of Nonlinear Oscillations Theory, Archimede, New Studies in the History and Philosophy of Science and Technology, vol. 49, Springer, New York, 2017. [5] Gleick, J., Chaos: Making a New Science, New York, Viking, 1987. [6] Goodwin, R. M., “The Non-linear Accelerator and the Persistence of Business Cycle,” Econometrica 19 (1951) pp. 1–17. [7] Grassberger, P., & Procaccia, I., “Measuring the strangeness of strange attractors,” Physica D: Nonlinear Phenomena, 9(1-2), (1983) pp. 189-208. [8] Hamburger, L., “Analogie des fluctuations économiques et des oscillations de relaxation,” Indices du Mouvement des Affaires 9 - supplément, no Janvier (1933) pp. 1–35. [9] Hamburger, L., “Note on economic cycles and relaxation-oscillations,” Econometrica 2(1) (1933) pp. 112. [10] Hurwitz, A., “”Uber die Bedingungen, unter welchen eine Gleichung nur Wurzeln mit negativen reellen Theilen besitzt,” Math. Ann. 41 (1893) pp. 403–442. [11] Klein, M. & Baier, G., Hierachies of dynamical systems, In A Chaotic Hierarchy, edited by G. Baier and M. Klein. Singapore: World Scientific, 1991. [12] Kuznetsov, Y., Elements of applied bifurcation theory, Applied Mathematical Sciences, Vol. 112, Springer-Verlag, New York, (2004). [13] Le Corbeiller, Ph., “Les systèmes autoentretenus et les oscillations de relaxation,” Econometrica 1 (1933) pp. 328–332. [14] Li, T. Y. & Yorke, J. A., “Period Three Implies Chaos,” The American Mathematical Monthly, 82, (10) (Dec., 1975) pp. 985-992. [15] Lorenz, E. N., “Deterministic Nonperiodic Flow,” Journal of Atmospheric Sciences, 20(7), (1963) pp. 130–141. [16] Lyapunov, A.M., The General Problem of the Stability of Motion (In Russian), Doctoral dissertation, Univ. Kharkov 1892, English translations: The General Problem of the Stability of Motion, (A. T. Fuller trans.) Taylor & Francis, London 1992, pp. 250. [17] Rocard, Y., Théorie des Oscillateurs (Theory of Oscillators), Édition de la Revue Scientifique, Paris, 1941. [18] Routh, E. J., A Treatise on the Stability of a Given State of Motion: Particularly Steady Motion, Macmillan and co, 1877. [19] Sandri M., “Numerical Calculation of Lyapunov Exponents,” The Mathematica Journal, 6(3) (1996) pp. 78-84. [20] Sprott, J. C., Chaos and Time-Series Analysis, Oxford University Press, Oxford, 2003. [21] Sprott, J. C, Wang, X. and Chen, G., “Coexistence of point, periodic and strange attractors,” Internat. J. Bifur. Chaos Appl. Sci. Engrg. 23 (2013), no. 5, 1350093, 5 pp. [22] Van der Pol, B., On “relaxation-oscillations”, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, (VII), 2 (1926) pp. 978-992. [23] Velupillai, K. V., Nonlinear and Mutisectoral Macrodynamics: Essays in Honour of Richard Goodwin, (ed. K.V. Velupillai), Macmillan, London, 1989. [24] Wolf, A., Swift, J.B., Swinney, H.L. & Vastano, J.A., “Determining Lyapunov Exponents from a Time Series,” Physica D, 16 (1985) pp. 285-317.