A Numerical-Analytical Method for Constructing Periodic Solutions of the Lorenz System
aa r X i v : . [ m a t h . NA ] F e b dxdt ✻ ✛✲ ❄ DIFFERENTIAL EQUATIONSANDCONTROL PROCESSESN. 4, 2020Electronic Journal,reg. N Φ C77-39410 at 15.04.2010ISSN 1817-2172http://diffjournal.spbu.ru/e-mail: jodiff@mail.ru
Numerical methods
A Numerical-Analytical Method for Constructing PeriodicSolutions of the Lorenz System
Alexander N. PchelintsevTambov State Technical University,ul. Sovetskaya 106, Tambov, 392000, Russiae-mail: [email protected]
Abstract.
This article describes a method for constructing approximations toperiodic solutions of dynamic Lorenz system with classical values of the systemparameters. The author obtained a system of nonlinear algebraic equations ingeneral form concerning of the cyclic frequency, constant terms and amplitudesof harmonics that make up harmonic approximations to the desired solutions.The initial approximation for the Newton method is selected, which convergesto a solution describing a periodic solution different from the equilibrium posi-tion. The results of a computational experiment are presented. The results areverified using high-precision calculations.
Keywords:
Attractor, Lorenz Attractor, Trigonometric Polynomial, Newton’sMethod. ifferential Equations and Control Processes,
N. 4, 2020
Let us consider the nonlinear system of differential equations introduced by E.Lorenz in [1] ˙ x = σ ( x − x ) , ˙ x = rx − x − x x , ˙ x = x x − bx , (1)where σ = 10, r = 28 and b = 8 / X ( t ) = [ x ( t ) x ( t ) x ( t )] T . It is proved in the article[1] that there exists a number C > X ( t ) of thesystem (1), starting at time moment, | X ( t ) | < C , and the divergence of thevector velocity field of the system (1) is negative everywhere in R for classicalvalues of the system parameters. Then [1] there exists a limit set, called theLorenz attractor, to which all trajectories of the dynamical system are attractedwhen time tends to infinity. Thus the attractor determines the behavior of thesolutions of a dynamical system over large segments of time.W. Tucker in his work [2] proved that the attractor is hyperbolic in thesystem (1), that is, the attractor consists of cycles everywhere dense on it alongwhich the near trajectories diverge exponentially. This creates their chaoticbehavior.As know [3, 4], the symbolic dynamics is used to track cycles in the Lorenzsystem. The region in the phase space containing the attractor is divided intoa finite number of subdomains. Denoting each partition element by a symbol,the trajectories on the attractor passing through the corresponding regions areencoded by sequences of such symbols. If the sequence has regularity (repeata-bility of groups of characters), then we have a cycle. However, the return oftrajectories in a neighborhood of its part does not mean its closure. A critiqueof the results of such computational experiments can be found, for example, in[5]. In 2004, D. Viswanath published the paper [6], in which he presented theinitial conditions and periods for three cycles in the Lorenz attractor with ahigh accuracy. The calculation algorithm is based on the Lindstedt-Poincar´e(LP) method, which (unlike numerical integration methods) is not affected bythe stability of the cycle to which approximations are constructed. Electronic Journal. http://diffjournal.spbu.ru/ ifferential Equations and Control Processes, N. 4, 2020
An analysis of the Viswanath’s articles [6, 7] showed that the author givesa general description of the algorithm without reference to the computer im-plementation (in MATLAB as indicated in his works). Moreover, it is not clearhow the obtained inhomogeneous linear system of differential equations withperiodic coefficients is symbolically solved by the LP-method. For example,this can be done for the Van der Pol equation without any special problems.In the article [6] Viswanath showed data that can be verified by solving theCauchy problem with high-precision numerical methods (for example, [8]), butthe details of the algorithm are not disclosed.Therefore, it is important here to obtain the values of the initial conditionsand the period with a given accuracy, having described in detail the implemen-tation of the cycles search algorithm in the system (1).The goal of this article is to develop a numerical-analytical method forconstructing approximations to periodic solutions of the Lorenz system, whichis simpler to implement than the LP-method. In this case, a system of nonlinearalgebraic equations concerning of the cyclic frequency, constant terms, andamplitudes of harmonics making up the desired solution will be obtained ingeneral form.
Attempts to construct approximate periodic solutions in the system (1) withwere made before Viswanath (for example, [9]) by the method of harmonicbalance, but with low accuracy in representing real numbers, while in the article[9] initial conditions and periods of found cycles are not indicated (only drawingswith cycles are given). Now this method is actively developing in the works of[10, 11, 12] A. Luo to find periodic solutions of nonlinear systems of differentialequations.Next, we describe a numerical-analytical method for constructing approxi-mations to periodic solutions of the system (1). We make for this an approxi-mation of the phase coordinates on the period T by trigonometric polynomialsin general form with an unknown cyclic frequency ω (since we do not know the Electronic Journal. http://diffjournal.spbu.ru/ ifferential Equations and Control Processes, N. 4, 2020 value of T ; in the general case, it can be an irrational number): x ( t ) ≈ ˜ x ( t ) = x , + h X i =1 ( c ,i cos( iωt ) + s ,i sin( iωt )) ,x ( t ) ≈ ˜ x ( t ) = x , + h X i =1 ( c ,i cos( iωt ) + s ,i sin( iωt )) ,x ( t ) ≈ ˜ x ( t ) = x , + h X i =1 ( c ,i cos( iωt ) + s ,i sin( iωt )) , where h is given number of harmonics. If i > h , then we assume c ,i = s ,i = c ,i = s ,i = c ,i = s ,i = 0 . (2)By the right-hand side of the system (1), we compose the residuals δ ( t ) = ˜ x ′ ( t ) − σ [˜ x ( t ) − ˜ x ( t )] ,δ ( t ) = ˜ x ′ ( t ) − [ r ˜ x ( t ) − ˜ x ( t ) − ˜ x ( t )˜ x ( t )] ,δ ( t ) = ˜ x ′ ( t ) − [˜ x ( t )˜ x ( t ) − b ˜ x ( t )] , where the prime denotes the time derivative of the function. If we make calcu-lations in an analytical form, then for each residual you need the following:1. Differentiate by time the corresponding trigonometric polynomial;2. Where there are products of phase coordinates, multiply the correspond-ing trigonometric polynomials, converting the products of trigonometricfunctions into sums;3. Give similar terms for each function cos() and sin() with the correspondingargument;4. By virtue of the equalities (2), to cut off the higher-order harmonics fromthe resulting residual;5. Set the resulting residual to zero, i.e., coefficients at its harmonics.If we put together the found algebraic equations for each residual, we obtaina still unclosed system of nonlinear equations concerning of unknown amplitudes c ,i , s ,i , c ,i , s ,i , c ,i and s ,i ( i = 1 , h ), constant terms x , , x , and x , andthe cyclic frequency ω . The number of unknown variables in the system is3(1 + 2 h ) + 1 = 6 h + 4, but the equations are one less. Electronic Journal. http://diffjournal.spbu.ru/ ifferential Equations and Control Processes, N. 4, 2020
An additional equation can be taken from the following considerations. Itis known (see [4, 6]) that the desired cycles intersect the plane passing throughthe equilibrium positions of the system (1) O (cid:16) − p b ( r − , − p b ( r − , r − (cid:17) , O (cid:16)p b ( r − , p b ( r − , r − (cid:17) (3)and parallel to the plane x Ox (a Poincare section). Then the third coordinatein the initial condition for the desired cycles is equal to r −
1, whence ˜ x (0) = r − x , + h X i =1 c ,i −
27 = 0 . The author did not find in literature of other additional information on theperiodic solutions in the Lorenz system. Note that for the three cycles foundby Viswanath, in the initial condition for the third coordinate, the number 27was taken.
Electronic Journal. http://diffjournal.spbu.ru/ ifferential Equations and Control Processes, N. 4, 2020
Next, we give an example of a system of equations for h = 2: ωs , − c , + 10 c , = 0 , − s , + 10 s , − c , ω = 0 , ωs , − c , + 10 c , = 0 , − s , + 10 s , − c , ω = 0 , x , − x , = 0 ,c , x , + c , x , + s , s , s , s , ωs , + c , c , c , c , c , − c , = 0 ,s , x , + s , x , + c , s , − c , s , s , + c , s , − c , s , − s , − c , ω = 0 ,c , x , + c , x , − s , s , ωs , + c , c , c , − c , = 0 ,s , x , + s , x , + c , s , s , − s , + c , s , − c , ω = 0 ,x , x , + x , − x , + s , s , s , s , c , c , c , c , , − c , x , − c , x , + ωs , − s , s , − s , s , c , − c , c , − c , c , , − s , x , − s , x , + 8 s , − c , s , c , s , − c , s , c , s , − c , ω = 0 , − c , x , − c , x , + 2 ωs , + s , s , c , − c , c , , − s , x , − s , x , + 8 s , − c , s , − c , s , − c , ω = 0 , x , − x , x , − s , s , − s , s , − c , c , − c , c , ,x , + c , + c , −
27 = 0 . Note that for any h a similar system has solutions x , = x , = ± p b ( r − , x , = r − , c k,i = 0 , s k,i = 0 ,ω is any number , k = 1 , , i = 1 , h, corresponding to the equilibrium positions (3).Therefore the resulting nonlinear system of algebraic equations has a non-unique solution. To find its approximate solutions, we will use the Newtonnumerical method, whose a convergence to the desired solution (i.e., describinga periodic solution of the system (1) different from its the equilibrium positions)depends on the choice of the initial approximation. Electronic Journal. http://diffjournal.spbu.ru/ ifferential Equations and Control Processes, N. 4, 2020
Thus, to obtain an approximation to the periodic solution, we must obtaina nonlinear system concerning of unknown decomposition coefficients and fre-quencies. As shown in the previous section, even for two harmonics, the systemhas a bulky form. Therefore, we consider the algorithm for performing symboliccalculations to obtain it.When developing software [13], the Maxima math package (a computeralgebra system) was chosen. The program for obtaining the amplitudes andconstant terms of the residuals for h = 2 is presented below. /* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*//* [wxMaxima: input start ] */display2d:false$x1:x10+c1c1*cos(1*omega*t)+s1c1*sin(1*omega*t)+c1c2*cos(2*omega*t)+s1c2*sin(2*omega*t)$x2:x20+c2c1*cos(1*omega*t)+s2c1*sin(1*omega*t)+c2c2*cos(2*omega*t)+s2c2*sin(2*omega*t)$x3:x30+c3c1*cos(1*omega*t)+s3c1*sin(1*omega*t)+c3c2*cos(2*omega*t)+s3c2*sin(2*omega*t)$assume(omega > 0)$delta1:trigreduce(diff(x1,t)-(10*(x2-x1)),t)$delta2:trigreduce(diff(x2,t)-(28*x1-x2-x1*x3),t)$delta3:trigreduce(diff(x3,t)-(x1*x2-8/3*x3),t)$expand(diff(delta1,cos(1*omega*t)));expand(diff(delta1,sin(1*omega*t)));expand(diff(delta1,cos(2*omega*t)));expand(diff(delta1,sin(2*omega*t)));expand(integrate(delta1,t,0,2*%pi/omega)*omega/(2*%pi));expand(diff(delta2,cos(1*omega*t)));expand(diff(delta2,sin(1*omega*t)));expand(diff(delta2,cos(2*omega*t)));expand(diff(delta2,sin(2*omega*t)));expand(integrate(delta2,t,0,2*%pi/omega)*omega/(2*%pi));expand(diff(delta3,cos(1*omega*t)));expand(diff(delta3,sin(1*omega*t)));expand(diff(delta3,cos(2*omega*t))); Electronic Journal. http://diffjournal.spbu.ru/ ifferential Equations and Control Processes, N. 4, 2020 expand(diff(delta3,sin(2*omega*t)));expand(integrate(delta3,t,0,2*%pi/omega)*omega/(2*%pi));/* [wxMaxima: input end ] */
The expression display2d:false$ turns off multi-line drawing of fractions,degrees, etc. The sign $ allows to calculate the result of an expression, butnot display it (instead of ; ). The function trigreduce(expression,t) col-lapses all products of trigonometric functions concerning of the variable t ina combination of sums. Differentiation of residuals according to harmonicfunctions is necessary to obtain the corresponding amplitudes. The function expand(expression) expands brackets (performs multiplication, exponentia-tion, leads similar terms).To find the constant terms of the residuals, their integration over the periodis applied, i.e. the constant term of the k -residual is ω Z πω δ k ( t ) dt π . So that during symbolic integration the package does not ask a questionabout the sign of the frequency, a command is given assume(omega > 0)$ .A file with package commands is generated similarly for any number of h harmonics by a computer program written in C++ [13]. After executing thisprogram, the package will output symbolic expressions to the console for theleft side of the system of algebraic equations, which will be solved in it by theNewton method.Note that the most time-consuming operation here is symbolic integration.For example, for 120 harmonics, the system formation time is more than 2 days.We can here parallelize the computational process on three computers, but thiswill not have a significant effect. Therefore, a system of algebraic equationsmust be formed immediately. Next, we get a general form of this system. Notethat when solving the system of nonlinear equations by the Newton method,the Jacobi matrix for the left side of the system does not invert. The Maximapackage uses LU decomposition to solve a system of linear equations at eachiteration of the method. Electronic Journal. http://diffjournal.spbu.ru/ ifferential Equations and Control Processes, N. 4, 2020
Since the right-hand side of the (1) system contains nonlinearities in the formof products of phase coordinates, let us obtain relations expressing the coeffi-cients of trigonometric polynomials obtained by multiplying the approximations˜ x ( t )˜ x ( t ) and ˜ x ( t )˜ x ( t ).We consider two functions f ( t ) and F ( t ) represented by Fourier series f ( t ) = a + ∞ X i =1 ( a i cos( iωt ) + b i sin( iωt )) ,F ( t ) = A + ∞ X i =1 ( A i cos( iωt ) + B i sin( iωt )) . Let f ( t ) F ( t ) = α + ∞ X i =1 ( α i cos( iωt ) + β i sin( iωt )) . Following the book [14, pp. 123-125], we have the relations: α = a A + 12 ∞ X m =1 ( a m A m + b m B m ) ,α i = a A i + 12 ∞ X m =1 ( a m ( A m + i + A m − i ) + b m ( B m + i + B m − i )) , (4) β i = a B i + 12 ∞ X m =1 ( a m ( B m + i − B m − i ) − b m ( A m + i − A m − i )) . (5)We assume that for i > ha i = b i = A i = B i = 0 . Since for our problem we find for an approximation up to and including the h -harmonic, we zero all the amplitudes in the product for i > h , i.e. α i = β i = 0 . Thus, we pass from the product of series to the product of trigonometricpolynomials. Also in the relations (4) and (5) we will assume [14, p. 124] that A m − i = A i − m , B m − i = − B i − m , B = 0 . Electronic Journal. http://diffjournal.spbu.ru/ ifferential Equations and Control Processes, N. 4, 2020
Then we get α = a A + 12 h X m =1 ( a m A m + b m B m ) ,α i = a A i + 12 ∞ X m =1 a m A m + i + 12 ∞ X m =1 a m A m − i + 12 ∞ X m =1 b m B m + i + 12 ∞ X m =1 b m B m − i == a A i + 12 h − i X m =1 a m A m + i + 12 a i A + 12 i − X m =1 a m A i − m + 12 h X m = i +1 a m A m − i ++ 12 h X m =1 b m B m + i + 12 b i B − i − X m =1 b m B i − m + 12 h X m = i +1 b m B m − i == a A i + a i A + 12 h − i X m =1 ( a m A m + i + b m B m + i ) + 12 i − X m =1 ( a m A i − m − b m B i − m ) ++ 12 h X m = i +1 ( a m A m − i + b m B m − i ) ,β i = a B i + 12 ∞ X m =1 a m B m + i − ∞ X m =1 a m B m − i − ∞ X m =1 b m A m + i + 12 ∞ X m =1 b m A m − i == a B i + 12 h − i X m =1 a m B m + i + 12 i − X m =1 a m B i − m − h X m = i +1 a m B m − i −− h − i X m =1 b m A m + i + b i A + 12 i − X m =1 b m A i − m + 12 h X m = i +1 b m A m − i == a B i + b i A + 12 h − i X m =1 ( a m B m + i − b m A m + i ) + 12 i − X m =1 ( a m B i − m + b m A i − m ) ++ 12 h X m = i +1 ( − a m B m − i + b m A m − i ) . Applying the obtained formulas to calculate the products of trigonometricpolynomials to the residuals, we can write the equations for the i -th harmonics( i = 1 , h is the number of harmonics, k = 1 , Electronic Journal. http://diffjournal.spbu.ru/ ifferential Equations and Control Processes, N. 4, 2020 k = 1: iωs ,i − c ,i + 10 c ,i = 0 , − iωc ,i − s ,i + 10 s ,i = 0 , the equation corresponding to the constant term for the first residual is x , − x , = 0 ,k = 2: iωs ,i − c ,i + c ,i + x , c ,i + c ,i x , + 12 h − i X m =1 ( c ,m c ,m + i + s ,m s ,m + i ) ++ 12 i − X m =1 ( c ,m c ,i − m − s ,m s ,i − m ) ++ 12 h X m = i +1 ( c ,m c ,m − i + s ,m s ,m − i ) = 0 , − iωc ,i − s ,i + s ,i + x , s ,i + s ,i x , + 12 h − i X m =1 ( c ,m s ,m + i − s ,m c ,m + i ) ++ 12 i − X m =1 ( c ,m s ,i − m + s ,m c ,i − m ) ++ 12 h X m = i +1 ( − c ,m s ,m − i + s ,m c ,m − i ) = 0 , the equation corresponding to the constant term for the second residual is − x , + x , + x , x , + 12 h X m =1 ( c ,m c ,m + s ,m s ,m ) = 0 ,k = 3: iωs ,i − x , c ,i − c ,i x , − h − i X m =1 ( c ,m c ,m + i + s ,m s ,m + i ) −− i − X m =1 ( c ,m c ,i − m − s ,m s ,i − m ) −− h X m = i +1 ( c ,m c ,m − i + s ,m s ,m − i ) + 83 c ,i = 0 , Electronic Journal. http://diffjournal.spbu.ru/ ifferential Equations and Control Processes, N. 4, 2020 − iωc ,i − x , s ,i − s ,i x , − h − i X m =1 ( c ,m s ,m + i − s ,m c ,m + i ) −− i − X m =1 ( c ,m s ,i − m + s ,m c ,i − m ) −− h X m = i +1 ( − c ,m s ,m − i + s ,m c ,m − i ) + 83 s ,i = 0 , the equation corresponding to the constant term for the third residual is − x , x , − h X m =1 ( c ,m c ,m + s ,m s ,m ) + 83 x , = 0 , the additional system equation is x , + h X i =1 c ,i −
27 = 0 . As a result of numerous computational experiments, the initial approximationwas chosen for the cyclic frequency, constant terms, and amplitudes at h = h = 5: ω = 4 , x , = x , = x , = 0 , c ,i = − , i = 1 , ,s ,j = 0 , j = 1 , , , , s , = 1 . This result is remarkable in that the Newton method converges to a solutiondifferent from the equilibrium positions. Therefore, to improve the accuracy ofthe approximate periodic solution, we consider a system of algebraic equationsfor the value of h equal to some h > h .The obtained numerical solution of the system for h = h is taken as theinitial approximation for amplitudes with indices i ≤ h for a system with h = h , and the values of the initial approximation for amplitudes with indices i > h are assumed to be zero.Tables 1–3 show the result of solving the system for h = 35; the accu-racy of the Newton method is 10 − . The period value is obtained equal to T = 1 . Electronic Journal. http://diffjournal.spbu.ru/ ifferential Equations and Control Processes, N. 4, 2020
Table 1: The amplitudes of harmonics for ˜ x ( t ), x , = 0. i c ,i s ,i − . . . . . − . − . − . − . . . . . − . − . · − − . · −
16 0 017 − . · − . · −
18 0 019 4 . · − . · −
20 0 021 1 . · − − . · −
22 0 023 − . · − − . · −
24 0 025 − . · − . · −
26 0 027 1 . · − . · −
28 0 029 5 . · − − . · −
30 0 031 − . · − − . · −
32 0 033 − . · − . · −
34 0 035 5 . · − . · − Electronic Journal. http://diffjournal.spbu.ru/ ifferential Equations and Control Processes, N. 4, 2020
Table 2: The amplitudes of harmonics for ˜ x ( t ), x , = 0. i c ,i s ,i − . . . − . − . − . − . . . . . − . − . − . − . . . · − . . · − − . · −
20 0 021 − . · − − . · −
22 0 023 − . · − . · −
24 0 025 5 . · − . · −
26 0 027 2 . · − − . · −
28 0 029 − . · − − . · −
30 0 031 − . · − . · −
32 0 033 2 . · − . · −
34 0 035 1 . · − − . · − Electronic Journal. http://diffjournal.spbu.ru/ ifferential Equations and Control Processes, N. 4, 2020
Table 3: The amplitudes of harmonics for ˜ x ( t ), x , = 23 . i c ,i s ,i . − . − . − . − . . . . . − . − . − . − . . . . · −
17 0 018 2 . · − − . · −
19 0 020 − . · − − . · −
21 0 022 − . · − . · −
23 0 024 1 . · − . · −
25 0 026 1 . · − − . · −
27 0 028 − . · − − . · −
29 0 030 − . · − . · −
31 0 032 1 . · − . · −
33 0 034 6 . · − − . · −
35 0 0
Electronic Journal. http://diffjournal.spbu.ru/ ifferential Equations and Control Processes, N. 4, 2020 -30 -20 -10 0 10 20 30 -30-20-10 0 10 20 30 0 5 10 15 20 25 30 35 40x3 x1 x2x3
Figure 1: The cycle obtained by described method.
Electronic Journal. http://diffjournal.spbu.ru/ ifferential Equations and Control Processes, N. 4, 2020 solution is ˜ x (0) = − . , ˜ x (0) = 2 . , ˜ x (0) = 27 . (6)The initial values (6) were checked on the period in a computer program thatimplements the numerical integration of the system (1) by the modified powerseries method [8] with an accuracy of estimating the common term of the series10 − , 100 bits for mantissa real number and machine epsilon 1 . · − .With such parameters of the method, the approximate values of the phasecoordinates obtained by numerical integration were also verified by the samenumerical method, but in reverse time. The values in the reverse time coincidewith (6) up to the 9th character inclusive after the point. The resulting valuesof x ( T ), x ( T ) and x ( T ) coincide with (6) up to the 8th character inclusive.The cycle corresponding to (6) is shown in Fig. 1. Note that the cyclefound coincides with the first Viswanath cycle in [6], all signs after the pointfor T also coincide with the data from [6]. The reported study was funded by RFBR according to the research project20-01-00347.
References [1] Lorenz, E. N. Deterministic Nonperiodic Flow, Journal of the AtmosphericSciences, vol. 20, no. 2 (1963), pp. 130-141.[2] Tucker, W. A Rigorous ODE Solver and Smale’s 14th Problem, Founda-tions of Computational Mathematics, vol. 2, no. 1 (2002), pp. 53-117.[3] Rabinovich, M. I. Stochastic Self-Oscillations and Turbulence, SovietPhysics Uspekhi, vol. 21, no. 5 (1978), pp. 443-469.[4] Galias, Z., Tucker, W. Validated Study of the Existence of Short Cyclesfor Chaotic Systems Using Symbolic Dynamics and Interval Tools, Interna-tional Journal of Bifurcation and Chaos, vol. 21, no. 2 (2011), pp. 551-563.[5] Lozi, R. Can We Trust in Numerical Computations of Chaotic Solutionsof Dynamical Systems?, Topology and Dynamics of Chaos. In Celebration
Electronic Journal. http://diffjournal.spbu.ru/ ifferential Equations and Control Processes, N. 4, 2020 of Robert Gilmore’s 70th Birthday. - World Scientific Series in NonlinearScience Series A, vol. 84 (2013), pp. 63-98.[6] Viswanath, D. The Fractal Property of the Lorenz Attractor, Physica D:Nonlinear Phenomena, vol. 190, no. 1-2 (2004), pp. 115-128.[7] Viswanath, D. The Lindstedt-Poincare Technique as an Algorithm forComputing Periodic Orbits, SIAM Review, vol. 43, no. 3 (2001), pp. 478-495.[8] Pchelintsev, A. N. Numerical and Physical Modeling of the Dynamics ofthe Lorenz System, Numerical Analysis and Applications, vol. 7, no. 2(2014), pp. 159-167.[9] Neymeyr, K., Seelig, F. Determination of Unstable Limit Cycles in ChaoticSystems by Method of Unrestricted Harmonic Balance, Zeitschrift f¨urNaturforschung A, vol. 46, no. 6 (1991), pp. 499-502.[10] Luo, A. C. J., Huang, J. Approximate Solutions of Periodic Motions inNonlinear Systems via a Generalized Harmonic Balance, Journal of Vibra-tion and Control, vol. 18, no. 11 (2011), pp. 1661-1674.[11] Luo, A. C. J. Toward Analytical Chaos in Nonlinear Systems, John Wiley& Sons, Chichester, ISBN: 978-1-118-65861-1, 2014, 258 pp.[12] Luo, A. C. J., Guo, S. Analytical Solutions of Period-1 to Period-2 Mo-tions in a Periodically Diffused Brusselator, Journal of Computational andNonlinear Dynamics, vol. 13, no. 9, 090912 (2018), 8 pp.[13] Pchelintsev, A. N. The Programs for Finding of Pe-riodic Solutions in the Lorenz Attractor, GitHub, https://github.com/alpchelintsev/periodic_sols [14] Tolstov, G. P. Fourier Series, Dover Publications, New York (1962), 336pp.