Finding and Visualizing of Limit Cycles Yuri K.Dem'yanovich Aleksey A.Fefelov Department of Electrical Departmen of Mathematics Emperor Alexander I and Mechanics St. Petersburg State St. Petersburg State University Transport University Saint Petersburg, Russia Saint Petersburg, Russia fefaleksey@mail.ru Yuri.Demjanovich@gmail.com Abstract of various devices reduced to the definition of steady states for a certain system of differential equations. The article reflects a study aimed at using Inadequate investigation of the sustainability can lead a parallel computing system for automated to catastrophic consequen\-ces for designed devices retrieval and visualization of cycles of a and apparatuses. Hopping from a planned state of quadratic system of two differential equations. resistance to an unplanned condition repeatedly led to The study was conducted in the seven- the destruction of bridges and structures, railway dimensional space of parameters --- system accidents, the total destruction of aircraft, etc. coefficients and initial data of the Cauchy Let a pair of functions x (t), y (t) be a solution to problem. It is very important for sustainable the Cauchy problem for a system of two differential operation of transport systems. To implement equations of the form the calculations, supercomputers of Moscow x'=P(x,y) (1) State University were remotely used. y'=Q(x,y), (2) Visualization of the results carried out x(0)=x0,y(0)=y0. (3) on Hewlett Packard personal computers. The The behavior of the trajectory of the point developed software model is applicable to (x (t), y (t)) is important in the plane (x, y) with weaning and visualization of cycles for increasing parameter t ( t usually represents the time different systems of two differential equations. of the corresponding physical system). The mentioned trajectories are called phase trajectories, and the plane 1. Introduction (x, y) is the phase plane. Steady state defined by stability points and limit cycles (the word "cycle" A.Poincare examined the geometric pattern means a closed curve). Limit cycles (attractors) are solutions of the differential equation. In 1900, characterized by the fact that they are approached D.Hilbert set the task of research limit cycles (wind on them) by phase trajectories (i.e. they are (attractors) in two-dimensional quadratic systems attraction cycles). Repulsion points and cycles may (see [Per01], [Li03], [Leo14], [Ruz15], [Leo15], also exist: phase trajectories are wound from them. [Leo17]). In the fifties A. N. Kolmogorov suggested Note that between two nested limit cycles (they estimating the number of these cycles. correspond steady states) there is always a cycle In his book [Arn05] V. I. Arnold wrote: "To repulsion (it corresponds to an unstable state). estimate the number of limit cycles quadratic vector The aforementioned tragedies show an urgent need fields on the plane, A. N. Kolmogorov distributed to develop reliable methods of finding attraction cycles several hundred such fields (with random selected and repulsion cycles. coefficients of polynomials of the second degree) In the second half of the twentieth century, a large to several hundred students at the mechanical and number of theoretical papers in which the existence of mathematical Faculty of Moscow State University as a limit cycles and the boundaries of the parameters are mathematical workshop. Each student had to find the indicated, where they should be searched (we skip the number of limit cycles of their field.The results of this review of these works). experiment were completely unexpected: there were no Unfortunately, in the vast majority of cases, there limit cycles at all." are no analytical methods (formulas) for determining It is well known that the definition of sustainability these cycles. In view of this, for finding cycles have Copyright c by the paper's authors. Use permitted under Creative become widely used in modern computing methods Commons License Attribution 4.0 International (CC BY 4.0). In: A. and computers. Khomonenko, B. Sokolov, K. Ivanova (eds.): Selected Papers of the Models and Methods of Information Systems Research Workshop, St. Petersburg, Russia, 4-5 Dec. 2019, published at http://ceur-ws.org 108 The exception of unplanned states of transport showed significant advantage of this method (relative devices, and bridge and tunnel structures come down to computational speed) compared to the high- to finding all attraction and repulsion cycles for phase precision Gear's method used by other researchers. In trajectories of the Cauchy problem this work, the cycles of attraction and repulsion cycles x'=x*x+x*y+y, (4) are automatically determined. y'=a*x*x+b*x*y+c*y*y+alpha*x+beta*y (5) Note that between every two nested attractors there x(0)=x0, y(0)=y0. (6) is repulsion cycle (the cycle of unstable equilibrium). The selection of the seven parameters appearing The definition of the location of these cycles is very here a, b, c, alpha,beta,x0,y0 should exclude device important in calculating stability in the case of jump from a planned attractor to an unplanned one. designing mechanisms and structures (unstable Such a jump may cause the device to malfunction and equilibrium cycle unsafe for designed devices). even to its complete destruction. For reliable results in Initial testing of algorithms and programs was points of the selected region of the seven-dimensional conducted in uniprocessor mode, and then with parameter space all cycles (attractive and repulsive) parallelization emulation with an MPI interface on a have to be defined, and they have to be presented with laptop and on a parallel cluster. After that, work was a video monitor. Problem (4) - (6) is a special case of a carried out remotely on supercomputers "Chebyshev" more general problem (1) - (3). and "Lomonosov-1" of Supercomputer Research The problem of finding and visualizing cycles was Computing Center of Moscow State University. The solved, thanks to the use of modern computing tools most interesting calculation results in the latter case and high-speed computing systems. To solve this were automatically saved on the supercomputer, problem, Professor G.A. Leonov formed a group then sent and autonomously visualized on the HP with researchers who conducted a series of numerical 27-p251ur All-in-One and on the HP Pavilon x360 experiments using various methods on computers of Convertable Notebook (Figures 1 -- 3 show the results various types. of some visualizations). Due to difficulties in processing seven-dimensional The results can be used in calculating and parameter spaces to these studies, the authors of this designing various devices, as well as for work were also involved in the organization of parallel checking reliability of created designs. The simple computing on a super\-computer. modification of algorithms and programs allows you This article reflects a study aimed at using to use the program in the case of solving similar a parallel computing system for automated retrieval problems for other autonomous systems of differential and visualization of cycles of a quadratic system of equations. two differential equations in the seven-dimensional space of parameters. The parameters are the 3 Results coefficients systems and initial data of the Cauchy problem for the mentioned system. When implementing calculations, supercompu- 3.1 First Series of Values Parameter beta ters "Chebyshev" and "Lomonosov-1" of Super- computer Research Computing Center of Moscow State University are remotely used. The visualization Search for limit cycles (attractors) for a set of of the obtained cycles is carried out on HP 27-p251 parameters a=-10.0, b=2.7, c=0.4, alpha=-473.5, ur All-in-One and HP Pavilon x360 Convertable beta=0.003-epsilon, epsilon=s*0.0000000001, Notebook PC. The developed software model is s=0,1,…, 1000, in each of these options led to 3 applicable to the finding and visualization of cycles attractors (the case epsilon=0 with a gradual expansion for different systems of two differential equations. of the study area, see Fig. 1 -- 3). 3.2 Second Series of Values Parameter beta 2 Methods and Algorithms To solve problem (3) - (4), the authors use the Runge- When searching for limit cycles (attractors) in another Kutta method of fourth order precision with the set parameters, namely, for a=-10.0, b=2.7, c=0.4, automatic choice of step. A numerical experiment alpha=-173.5,beta=0.004+ 0.0001*s, s=0,1,…,9950. 109 In each of the options listed, exactly one attractor space of parameters. The developed software model appeared. is obviously applicable to weaning and visualization of cycles for different systems of two differential equations. 3.3 Pair (b, c) gets 32 Million Values Here we consider 32 million of the pairs of parameter (b,c) according to the next formulas a=(b-1)*(b-1)/(4*(c-1)+1), b=2.1+0.0001*s, s=0,1,…,8000, c=0.5+0.0001*p, p=0,1,…,4000, alpha=a*(2+b)/(b*c-1) +|a*(2+b)/( b*c-1)|/2, beta=0. Here, in each variant, three attractors appeared, but for some parameters (and for small perturbation of the parameter \beta) a fourth attractor arose. The occurrence of the fourth attractor cannot be considered reliable because rounding errors are occurring when floating point is used. Therefore we will not discuss the fourth attractor. Figure 1: Visualization of wound trajectories 3.4 Wavelet Decomposition To speed up data transfer was considered a first-order wavelet decomposition with the following parameters a=-10.0, b=2.7, c=0.4, alpha=-173.5,beta=0.003. 1000 Cauchy problems were solved with the initial data (x_0, y_0) = (j, 0) ; here j = 1,2, \ldots, 1000. The resulting sets of values were saved and then divided into the main and wavelet streams (the main stream turned out to be about 2 times less source). Next, to another computer via ssh- protocol source and main streams are transferred. Let T0 be the transmission time of the original flow, let T1 be the transmission time of the main flow, and let k be their ratio, i.e. k = T0 / T1. In the described numerical experiment, the coefficient k turned out to be 1.92. Thus this indicates Figure 2: Expansion of the area. More wound savings significant resource when the wavelet trajectories decomposition is used. 4 Conclusion This study showed that the problem of finding and visualization of cycles in multidimensional (in seven- dimensional) space of parameters can be solved using modern computational algorithms and high-speed parallel computing systems. In particular, the application of the Runge-Kutta method proved to be very effective. The usage of a parallel computing system for automated search and visualization of cycles quadratic system of two differential equations gives very precise results in the seven-dimensional Figure 3: Visualization of three cycles 110 5 Acknowledgments This work was supported by Botan Investments and [Kuz15] Kuznetsov, N.V., Leonov, G.A., Shumafov, M.M.: Supercomputer Research Computing Center of A Short Survey on Pyragas Time-delay Feedback Stabilization and Odd Number Limitation. In: Moscow State University. Science\-Direct, 48-11, IFAC-PapersOnLine, 2015. P. 706-709. References [Leo15] Leonov, G. A., Zvyagintseva, K. A.: Pyragas Stabilization of Discrete Systems with Delayed [Per01] Perko, L.: Differential Equations and Dynamical Feedback and Pulse Periodic Gain. In: ISSN 10634541, Systems, Springer, N.Y., 2001. Vestnik St. Petersburg University. Mathe\-matics, [Li03] Li, J.: Hilbert’s 16th problem and bifurcations of 2015. Volume 48. Issue 3. P. 147–156. planar polynomial vector field. In: Internat. J. [Leo17] Leonov, Gennady A., Moskvin Alexander V.: Bifurcation Chaos, 2003. Volume 13. Issue 1. P. 47– Stabilizing Unstable Periodic Orbits of Dynamical 106. Systems Using Delayed Feedback Control with [Arn05] Arnold, V.I.: Experimental Mathematics. Fazis, Periodic Gain. In: International Journal Dynamics, Moscow, 2005. Control DOI 10.1007/s40435-017-0316-8 © Springer- Verlag Berlin Heidelberg, 2017. [Leo14] 8. Leonov, G.A.: Pyragas Stabilizability via Delayed Feedback with Periodic Control Gain. In: System \& Control Letters, 2014. Volume 69. P. 34-37. 111