The calculation process of the limit cycle for Lorenz system

: This work focuses on the bifurcation behavior before chaos phenomenon happens. Traditional numerical method is unable to solve the unstable limit cycle of nonlinear system. One algorithm is introduced to solve the unstable one, which is based one of the continuation method is called DEPAR approach. Combined with analytic and numerical method, the two stable and symmetrical equilibrium solutions exist through Fork bifurcation and the unstable and symmetrical limit cycles exist through Hopf bifurcation of Lorenz system. With the continuation algorithm, the bifurcation behavior and its phase diagram is solved. The results demonstrate the unstable periodical solution is around the equilibrium solution, besides the trajectory into the unstable area cannot escape but only converge to the equilibrium solution.


Introduction
Lorenz equation, nowadays its related phenomenon is famous known as 'Butterfly effect' [1][2][3].Originally in 1963, Lorenz's work focused on the numerical simulation of the dynamical model for the atmospheric prediction.His work demonstrated the results were rather sensitive to the initial conditions.Specifically, a slight variation in the initial condition of the equation might lead to a significant difference of the results, the phenomenon of which is called chaos.
Lorenz 's work was a milestone for later researchers.Many works focused on the attractors.Shil'Nikov A L et al. applied the norm forms theory to investigate the Lorenz attractors [4].Doering C R et al. studied the shape and dimension of the Lorenz attractor by the compution of the Lyapunov dimension with using numerical method [5].Stewart I discuss the whether the Lorenz attractors exist [6].Luzzatto S et al, focused on a class of geometric Lorenz flows and found the classical Lorenz attractor is mixing [7].
Moreover, the chaos phenomenon and application had also been widely studied from many different domains.For example, in mechanics, chaos helps prevent undesirable resonances [8].In thermodynamics, it offers one way to seek the relationship between the heat transports [9][10][11].Besides in chemical engineering, using chaotic theory may facilitate the implementation of novel operation/control strategies for the process industry [12].In living organisms, chaotic dynamics in vital functions can make the difference between health and disease [13].In Meteorology, it is common used to predict the weather [14].In biology, it has been suggested that the disappearance of chaos may be the signal of pathological behavior [15].In many such works, they demonstrate the Lorenz equation is the basis of kinds of nonlinear system in reality.
In this work, we focus on the computation of the limit cycle for the Lorenz system on the case before the attractor appearance, especially for the unstable limit cycle which is the limitation between two stable solutions (periodical or trivial ) but cannot be solved by regular algorithm.The unstable limit cycle is not able to be observed in reality for small perturbations, however, it does exists ideally and theoretically.A few works focused on this issue.For example, Gao Xuejun used one of continuation method to solve the unstable cycle nonlinear rail vehicle system [16].Krasnoschechenko V I discussed the stabilization of the unstable limit cycle of one relay chaotic system [17].Therefore, in this work we introduce the DEPAR approach to compute the unstable limit cycle for Lorenz system, to get a better understanding of the bifurcation behavior before chaos happens.

b. The Hopf bifurcation point
The Jacobian matrix of equation (2.1) is According to Hopf bifurcation theorem [18], firstly, Let λ Solve it and we can obtain And another eigenvalue η = −σ − β − 1 < 0 And then calculate the real part with term (2.5), it yields The bifurcation is subcritical.
Figure 1 the limit cycle, L1 denotes the stable one while L2 denotes the unstable one, point O is the stable trivial solution With using the time integration and Poincare map [19] method, the stable solution (or the stable limit cycle if it exists) and the chaos solution can be solved as figure 2 (the unstable cycle cannot be solved by this method but introduced another algorithm in chapter 4).From the figure, the results may be rough and far from the reality.For example, according to the Analytic solution, chaos should happen in ρ >  H = 24.7 (The Hopf bifurcation point).However, when ρ = 21, chaos happens.This difference illustrates the time integration method has low precision with many factors such as the initial value, step size, accepted error, et al.
On the other hand, from the phase portrait (figure 3) and time domain (figure 4), when ρ = 13.9 or 24<ρ H , chaos dose not happen, the trajectory are converge to the stable solution C ± .When ρ = 25>ρ H , points C ± are not stable as shown in figure 5 and 6.The two become strange attractors, which attract the trajectory nearby, but not converge to themselves.Although nowadays it is known the chaos phenomenon may happen in different cases.According to chapter Analytic solution, we know this system has subcritical Hopf bifurcation, which means the unstable limit cycle exists in this system.Besides, the stable limit cycle can be computed by the time integration method, however, the method fails to calculate the unstable one.Therefore, other algorithm may be used to solve this problem.For continuation method is widely applied in rail vehicle for limit cycle computation [16].

Continuation method for the unstable limit cycle
From figure 1, it is known the unstable cycle cannot be solved by the time integration method because it only can converge to the stable solution.But it is able to be calculated by other algorithm such as shooting or continuation ones [20].Considering the shooting method has lower computational efficiency, here we use one of continuation method, which is DERPAR algorithm.
In this paper ρ is the bifurcation parameter.We start the trivial solution which is known to be asymptotically stable at C 0 .Then we increase ρ in small steps and follow the solution and the eigenvalues of its Jacobian by the numerical solution of the system and numerical computation of eigenvalues for each value of ρ.When a bifurcation point is reached, the continuation algorithm is adopted to investigate the stationary and periodic solutions.The majority of continuation algorithms [21] are based on a predictor/corrector scheme.The equations of motion of system can be given as the ODE-formulation y = ẋ= f(x, ρ) (4.1)where ρ is the varied bifurcation parameter.The algorithm starts from a known point (yj-1, ρj-1) on a solution of Equation (4.1), as shown in Figure 7.The task in every step j (j=1,2,…) along the path is to find a neighboring solution point (yj, ρj) for a continued value ρj of the bifurcation parameter in two successive steps: For each ρ，the solution x(ρ) depends on ρ and satisfies f(x(ρ), ρ) = 0 (4.3) x(ρ ) may contain different branches and each branch begins at the point itself.
To solve (4.2), derivate (4.3) with parameter ρ.The the ODEs on x(ρ) can be obtained.If the Jocobian matrix J(x(ρ), ρ) is regular on [ρ 0 , ρ 1 ], the integral result x(ρ) by (4.6) satisfies But at the branch point, the matrix J is singular which make the continuation fail.
To make the algorithm can continue, it needs to revise the method.Furthermore, the revised method can totally depend on x(ρ) at any point and make the solution be a continuous and derivable curve in n+1 dimension (x  To determine z, the initial conditions of (4.9) are written as z = 0, x = x 0 , ρ = ρ 0 (4.11) For each z, (4.9) is one algebraic equations on n+1 unknown elements If Jk is regular for one z and k (1 ≤ k ≤ n + 1), for clarity let Then (4.9) about unknown elements Then the coefficient β i can be solved.For n(n+1) dimension matrix, combined The sign of the derivative dx k dz can be determined according to the direction of curve.Other derivatives can be computed according to (4.14).The whole process can be summarized as follow process.
(1) Input parameters and the initial (start) value ρ (2) In the correction step of solving (4.8), let the correction number NC=1 (3) Calculate f, J. Solve dx i dz and k, β i .(Corrector in figure 7)   8, it shows when ρ < 1, the system (2.1) has only one stable trivial solution (0, 0, 0); when ρ = 1, the fork bifurcation happens till ρ < 24.7 and at this case system (2.1) has one unstable trivial solution and two stable equilibrium solutions; when ρ = ρ H = 24.7, the subcritical Hopf birfurcation happens and each equilibrium solution bifurcates a unstable periodical solution (limit cycle) till the amplitude of the cycle is zero which also means the cycle and the trivial solution interests at ρ = 13.89 .Therefore, the above result at the range ρ ∈ [13.89, 24.7] demonstrates the bifurcation characteristic of system (2.1) before chaos happen (ρ > ρ H = 24.7).From figure 9, it is obvious that in phase plane, for one stable equilibrium solution C -, there is one unstable limit cycle nearby.It means the trajectory can only converge to Cbut not unstable limit cycle, which corresponds to the time domain diagram as shown in figure 4. The case for C + is similar.The phase diagram in 3D space for ρ = 20 Figure 10 shows in phase space when ρ = 20, system (2.1) has one unstable trivial solution (0, 0, 0), two stable equilibrium solutions C ± and two unstable limit cycles, wherein C ± and limit cycles are symmetrical on axis z.  11 shows dynamical phase diagram with parameter ρ increasing from 1 to ρ H .For clarity, it gives sets of curves of limit cycle.Actually from ρ ∈ [1, ρ H ], the limit cycles form torus, which make the trajectories cannot escape from the unstable limit cycle.

Conclusions
Recently, many study on the chaotic phenomenon and its application.This work mainly focuses on the Lorenz equation before the chaos phenomenon happens.Firstly the analytic way is used to obtain two bifurcation behaviours before chaos happen which are Fork and Hopf bifurcation respectively.Then the numerical approach is applied to solve the bifurcation diagram.Moreover, the two stable and symmetrical equilibrium solutions exist through Fork bifurcation and the unstable and symmetrical limit cycles exist through Hopf bifurcation.Finally, one special algorithm to solve the unstable limit cycle is introduced.The global bifurcation behaviour and its phase diagram before chaos of Lorenz system is computed by this algorithm and numerical method.From the result, the unstable periodical solution is around the equilibrium solution, besides the trajectory into the unstable area cannot escape but only converge to the equilibrium solution.It helps get a better understand mechanism of the chaos and the bifurcation behaviour before it occurs.

8 3
and ρ is the bifurcation parameter.Accoeding to (2.5), we can obtain ρ H = 24.7.The principle diagram for the limit cycle is shown in figure1.

Figure 2
Figure 2 The bifurcation and chaos phenomenon of equation (2.1)

Figure 5
Figure 5 The phase portrait x-y plane for ρ = 25 Figure 6 The time domain for ρ = 25

Figure 7
Figure 7 Continuation with a predictor/corrector scheme

5 . 3 ,
step (Predictor in figure 7).Do step (2) Results and discussions With using continuation method illustrated in chapter 4 to solve the bifurcation curves for Lorenz equation (2.1) with parameters σ = 10, β = 8 we can obtain the Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted: 18 February 2021 doi:10.20944/preprints202102.0419.v1following results wherein the bifurcation parameter is ρ. Figure 8 is the bifurcation curve, figure 9 are the phase diagram in 2D space for ρ = 20, and figure 10 is the phase diagram in 3D space for ρ = 20 and figure 11 is the phase diagram with dynamic parameter ρ.

Figure 8
Figure 8  The bifurcation curve From figure8, it shows when ρ < 1, the system (2.1) has only one stable trivial solution (0, 0, 0); when ρ = 1, the fork bifurcation happens till ρ < 24.7 and at this case system (2.1) has one unstable trivial solution and two stable equilibrium solutions; when ρ = ρ H = 24.7, the subcritical Hopf birfurcation happens and each equilibrium solution bifurcates a unstable periodical solution (limit cycle) till the amplitude of the cycle is zero which also means the cycle and the trivial solution interests at ρ = 13.89 .Therefore, the above result at the range ρ ∈ [13.89, 24.7] demonstrates the bifurcation characteristic of system (2.1) before chaos happen (ρ > ρ H = 24.7).

Figure 9
Figure 9 The phase diagram in 2D space for ρ = 20

Figure 10
Figure 10  The phase diagram in 3D space for ρ = 20 Figure10shows in phase space when ρ = 20, system (2.1) has one unstable trivial solution (0, 0, 0), two stable equilibrium solutions C ± and two unstable limit cycles, wherein C ± and limit cycles are symmetrical on axis z.

Figure 11
Figure 11  The phase diagram with dynamic parameter ρ Figure11shows dynamical phase diagram with parameter ρ increasing from 1 to ρ H .For clarity, it gives sets of curves of limit cycle.Actually from ρ ∈ [1, ρ H ], the limit cycles form torus, which make the trajectories cannot escape from the unstable limit cycle.