The study of chaotic systems began with the first chaotic system, the Lorenz equations, in 1963.^[1] As the first chaotic model, the Lorenz system is an important milestone in the history of chaos. Studying the Lorenz system has greatly promoted the development of chaos. To date, researchers have strictly proven the existence of the Lorenz attractor via Conley index theory and strict numerical calculation methods.^[2–4]

Many computer-assisted numerical results for Lorenz flow have been presented in subsequent years. A high-precision method of calculating periodic orbits on the Lorenz attractor was discussed in Ref. [5], in which the establishment of symbolic dynamics based on the return map was exploited. In Ref. [6], the bound of the generalized Lorenz system was obtained, based on the Lyapunov function and the Lagrange multiplier method. The complex dynamical behavior of a fractional-order Lorenz-like system and its control were investigated in Ref. [7]. Global bifurcations of the Lorenz manifold were explored in Ref. [8], in which 512 generic heteroclinic orbits were found that are given as the intersection curves of two-dimensional (2D) manifolds. In Ref. [9], time-delayed feedback control was used as a method of stabilizing unstable periodic orbits in the Lorenz equations. Recently, a novel hidden chaotic attractor in the classical Lorenz system was revealed in Ref. [10]. Time series of orbits, Lyapunov exponents, and bifurcations of the hidden chaotic attractor have also also reported. In Ref. [11] topological phase transitions when varying the parameters were studied. For special parameter values, the Lorenz flow will be topologically equivalent to an Anosov flow on a knot complement, and different knots appear for different parameter values. Dynamical analysis of the hyperchaos Lorenz system was investigated in Ref. [12], where the ultimate bound on the trajectories based on Lyapunov stability theory was investigated. Sparrow’s monograph provided a detailed discussion of the dynamic behavior of the Lorenz system.^[13]

In chaotic systems, owing to the randomness of motion, individual events cannot be predicted with precision in the long run. However, due to the ergodicity of dynamical motion, the average behavior under many initial conditions is still predictable. Cycle expansion is a powerful tool for calculating averages of physical observables in nonlinear dynamical systems by the periodic orbit theory.^[14,15] It was first applied to a spatially extended nonlinear system in 1997,^[16] in which a multi-point shooting method was used for finding unstable periodic orbits up to certain topological length and yielded several global averages characterizing the chaotic dynamics. For the Lorenz system, periodic orbit theory will continue to play an important role. The key in the application of periodic orbit theory is to establish appropriate symbolic dynamics; otherwise, even when some cycles are found, we still do not know the number of them that should exist and how they are related to each other. Ignoring a short periodic orbit will also impact the accuracy of calculations. In this paper, we put forth a general scheme to locate systematically the cycles of the Lorenz flow, and to build the one-dimensional (1D) symbolic dynamics effectively based on the orbital topology in the phase space for the classification of all periodic orbits.

The rest of this paper is organized as follows. In Section 2, we review basic dynamical properties of the Lorenz system. The extraction of periodic orbits in dissipative chaotic flows requires forceful numerical algorithms, and in Section 3 we present the variational method and apply it at our location for unstable periodic orbits of the Lorenz flow. In Section 4, we discuss the establishment of one-dimensional symbolic dynamics, which plays a key role in the classification of all short cycles up to a certain topological length. The evolution of periodic orbits with the change of parameters is also discussed. The results in the paper are summarized in the last section.

The Lorenz system consists of three ordinary differential equations, which are used to describe the thermal convective instability take the forms

Many studies have been done on the bifurcation behavior of the Lorenz system. The Hopf bifurcation will appear when the parameter c = 24.74, and the two fixed points S₋ and S₊ will become a strange attractor. At this moment the two equilibria comprise an unstable center, and when the phase space orbit approaches the equilibrium, it will be excluded, extending from the inside to the outside with the spiral. Moreover, due to the dissipative property of the system, the phase space must shrink as a whole. When the phase-space orbit expands to a certain extent, it will accordingly return to one of the fixed points and expand from the inside to outside around it. When it extends to a certain extent, it will suddenly randomly enter one of the fixed points, which causes the trajectories to hover between the two centers and never stop. Figure 1(a) shows the complex trajectory of the Lorenz system in two-dimensional (2D) phase space when we take the parameters for a = 10, b = 8/3, and c = 28.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) (a) 2D projection of a trajectory in the Lorenz flow at time t = 100. (b) Return map for the Poincaré section z = 27 after reduction to x > 0 by virtue of inversion symmetry.

The stability matrix of system (1) under the parameters taken is

Together with the eigenvalues of the stability matrix, the equilibria yield quite detailed information about the flow. The eigenvalues are obtained from the stability matrix (3):

In the vicinity of the S₀ fixed point, the extremely strong contraction λ⁽³⁾ = −22.83 along the e⁽³⁾ direction restricts the hyperbolic dynamics near S₀ to the plane spanned by the unstable eigenvector e⁽¹⁾ with λ⁽¹⁾ = 11.83, and the slowest contraction rate eigenvector e⁽²⁾ with λ⁽²⁾ = −2.67. In the plane, the strong expansion along the e⁽¹⁾ direction prevails over the slow λ⁽²⁾ = −2.67 contraction down the z axis, making it extremely unlikely for a random orbit to approach S₀. Therefore, trajectories close to S₀ are rare.

The S₋ eigenvalues indicate that the rotation period is T_S− = 2π/ω⁽¹⁾ = 0.62, and the turnover timescale is of order 1. The associated multipliers per spiral-out turn are Λ⁽¹⁾ = exp(μ⁽¹⁾T_S−) = 1.06 and Λ⁽³⁾ = exp (λ⁽³⁾T_S−) = 1.96 × 10⁻⁴. In the S₋ neighborhood, the unstable manifold trajectories slowly spiral out with a very small radial per-turn expansion multiplier Λ⁽¹⁾ = 1.06 and a very strong contraction multiplier Λ⁽³⁾ = 10⁻⁴ onto the unstable manifold.^[18] It is obvious that periodic orbits densely covering the strange attractor exist, and the Poincaré section technique can be used in cycle searching, which will probably reduce the 3D flow to a unimodal return map.

The Poincaré section is required since it cannot be tangential to the trajectories; thus, we choose the Poincaré section for a plane z = 27. Once the particular section has been selected, we can show the return map. Figure 1(b) displays the return map of the Lorenz system, and from the figure we can see that the return map is unimodal, which indicates the possibility for the establishment of 1D symbolic dynamics by finite Markov partitions,^[19] which will help us to encode all the cycles. In order to calculate the unstable periodic orbits systematically, in the next section we introduce a powerful method for extracting them in a chaotic system.

The variational method was proposed by Lan in 2004 to explore the periodic orbit in high-dimensional chaotic systems.^[20] This method not only preserves the robustness characteristics of most other methods, but it also has the characteristics of fast convergence when the search process is close to the real periodic orbit in practice. The physical idea of the method is that, first, a rough loop guess must be made based on the entire topology for the cycle to be searched, and then the initial guess loop must be driven to evolve toward the real periodic orbit by the variational algorithm. To achieve stability of the numerical method, the Newton descent method is used instead of the Newton–Raphson iteration method.

The variational method uses partial differential equation (5) to describe the evolution process of the initial guess loop toward the real cycle:

The important feature of Eq. (5) is that a minimizing cost functional exists monotonically as the loop guess evolves toward the periodic orbit:

For each iteration, the difference between the direction of loop velocity and flow velocity in the dynamical system decreases. The two become consistent when τ → ∞, so the real periodic orbit of the system is captured. The period of the orbit can be calculated by

Numerical implementation of the variational method requires use of finite-difference methods to discretize the loop guess. In order to ensure the numerical stability, we accurately discretize the loop derivatives:

A five-point approximation is used in the numerical work:

The discretized version of Eq. (5) with a fictitious time Euler step δτ is

The iteration of Eq. (11) can yield quickly and robustly the unstable periodic orbits for standard models of low-dimensional dissipative flows, which have an important feature, namely the existence of strange attractors. The method can also be applied to detect cycles in Hamiltonian flows and for spatio-temporal periodic solutions of partial differential equations. The method is effective for locating periodic orbits in high-dimensional phase space, especially when the low-dimensionality of the strange invariant set is incorporated.

In the previous work, we used the variational method effectively in calculating the periodic orbits in low- and high-dimensional systems, such as the Kuramoto–Sivashinsky equation and its steady-state solutions,^[22,23] the Rössler flow,^[24] and the Rydberg atom in crossed electromagnetic fields.^[25] It is obvious that this method is applicable to the Lorenz flow.

Here, we first discuss how to classify all the short periodic orbits of the Lorenz system when a = 10, b = 8/3, and c = 28 through the variational method, and then the evolution rule of the cycles with parameter changes is investigated. Two orbital segments can be used as basic building blocks, according to which we established 1D symbolic dynamics successfully to find all the periodic orbits up to certain topological length. The methodology is quite effective for this system and apparently applicable to other analogous dynamical systems, such as the Chen and Lü systems.^[26,27] In the numerical examples that follow, the convergence condition is F² ˂ 10⁻⁶.

4.2. Topological classification of periodic orbits of the Lorenz system

Using the initialization approach as previously mentioned, we found several short periodic orbits with simple topological structures in the system after a few trials. Figure 2(a) exhibits the simplest periodic orbit found by the variational method. We initialized the guess loop by numerical integration with the dynamical system, keeping the nearly close orbital segments, and then manually connected it into a loop. Although the guessed orbit is not smooth enough, the variational method will make the loop evolution the true periodic orbit of the system, as shown in the figure. By observing the topological structure of the orbit obtained, we established the one-dimensional symbolic dynamics. Figure 2(b) shows the periodic orbit of Fig. 2(a) projected onto the x–z plane. We can mark the orbital fragment winding near the left equilibrium S₋ with the symbol 0 and that winding near the right equilibrium S₊ with the symbol 1; therefore, the cycle in Fig. 2 is a 01 orbit. Not all conceivable symbol sequences are actually admissible as a dynamical trajectory according to the kneading theory.^[18] We also devised the initial loop guess for a 0 orbit; however, it does not converge to a genuine periodic orbit of the dynamical flow by the variational algorithm. For a 1 orbit, the same situation arises. Therefore, we did not find the presence of the 0 and 1 orbits, which means that they are pruned.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) (a) Shortest periodic orbit of Lorenz flow found by the variational method; blue line is the loop guess and red line is the periodic orbit. (b) 2D projection of the periodic orbit 01; the sequence of symbols corresponding to the orbit fragments is labeled. The two equilibria S− and S+ are marked with “+”.

The two orbital fragments can be used as basic building blocks for locating other complicated cycles. Searching for complex cycles with multiple circuits around the two fixed points requires more delicate initial conditions when using the variational method; otherwise, it will probably lead to non-convergence. We can initialize the loop guess for longer cycles constructed by cutting and gluing the short, known ones. For most systems, such a method yields a quite good systematic initial guess for longer cycles. Even if we deform the orbit manually into a closed loop, the variational method still shows its strength for good convergence.

In this way, we can initialize the loop guess corresponding to a longer sequence of symbols utilizing 1D symbolic dynamics. Figure 3(a) shows a 011 orbit with a topological length of 3 that is made of two 1 orbital fragments and one 0 orbital fragment. Figure 3(b) shows a 0011 orbit with a topological length of 4 that is made up of two circles around the two respective fixed points. Figures 3(c) and 3(d) show the cycles with a topological length of 5. Aided by the 1D symbolic dynamics, we can systematically locate all periodic orbits up to a certain topological length by first constructing the initial loop based on the symbolic sequence, and then evolving the loop to a true cycle using the variational method to check whether the corresponding symbolic dynamics is pruned.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) Four periodic orbits of the Lorenz system: (a) 011 orbit, (b) 0011 orbit, (c) 00001 orbit, and (d) 00101 orbit.

We constructed a 1D symbolic dynamics hierarchy for the Lorenz flow in an exhaustive manner. Altogether, we found all 12 admissible orbits up to a topological length of 5 and list them in Table 1. The Lyapunov exponent of a periodic orbit is defined as ln |Λ|/T_p, where Λ is the magnitude of its leading characteristic multiplier and T_p is its period. From the table we can see that the cycle 001 has the symmetry partner 011, which has the same period, and that the cycle 00101 has the symmetry partner 01011, while the cycle 0011 is conjugate to itself. This is determined by the z-axis symmetry of the system.

Table 1.

Table 1.

Table 1. Cycles up to topological length 5 for the Lorenz flow. Listed are the topological lengths, the itineraries p, periods Tp, the largest Lyapunov exponents λ and the coordinates x, y, z of one point on the periodic orbit. The cycles pruned are marked “−”. . Length p Tp λ x y z 1 0 – – – – – 1 – – – – – 2 01 1.558615 0.995 −8.606282 0.256056 35.754806 3 001 2.305912 0.961 −2.131913 −3.131322 15.653809 011 2.305912 0.961 2.131913 3.131322 15.653809 4 0001 3.023466 0.919 −0.107177 0.719410 19.296022 0011 3.084255 0.969 4.676212 −0.312933 29.448208 0111 3.023466 0.919 0.107177 −0.719410 19.296022 5 00001 3.725386 0.884 −0.663911 −4.837620 26.272703 00011 3.820246 0.946 0.869511 2.619909 21.704587 00101 3.869531 0.978 −5.803396 −9.394599 15.788917 00111 3.820246 0.946 −0.869511 −2.619909 21.704587 01011 3.869531 0.978 5.803396 9.394599 15.788917 01111 3.725386 0.884 0.663911 4.837620 26.272703

Table 1. Cycles up to topological length 5 for the Lorenz flow. Listed are the topological lengths, the itineraries p, periods Tp, the largest Lyapunov exponents λ and the coordinates x, y, z of one point on the periodic orbit. The cycles pruned are marked “−”. .

4.3. Evolution of periodic orbits with a change of parameters

Here, we investigate how the periodic orbits evolve upon changing the values of parameters a, b, and c. We first explore the 01 periodic orbit, which deforms with increasing a value. In our calculations, the 01 cycle calculated previously can be used as a initialization loop guess for the next a value; thus, we are able to continuously deform the 01 cycle. Figure 4(a) shows the deformation of the 01 orbit with four different a values and Table 2 lists the evolution. We find that the period of the 01 orbit decreases gradually with increasing a value.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) (a) Deformation of cycle 01 with four different a values. (b) Four different b values. (c) Four different c values.

Table 2.

Table 2.

Table 2. The periods Tp of the 01 orbit at different a, b, and c values. . a Tp b Tp c Tp 5 1.910341 2.1 1.685825 22 1.877405 7 1.667089 2.3 1.632464 38 1.275800 12 1.533349 2.8 1.538086 58 0.996013 15 1.524488 4.0 1.474188 88 0.793060

Table 2. The periods Tp of the 01 orbit at different a, b, and c values. .

We also studied the 01 orbit which deforms with increasing b and c values separately. Figure 4(b) shows the deformation with four different b values, while figure 4(c) shows that of the 01 orbit with four different c values. We find that the period of the 01 orbit also decreases gradually with increasing b and c values as listed in Table 2.

In this paper, we proposed a topological classification of the periodic orbits in the Lorenz system. We established 1D symbolic dynamics successfully and employed the variational method for cycle searching due to its numerical stability. Two basic building blocks are used to initialize the loop guess, and we calculated the long cycles by cutting and gluing the short, known blocks according to the symbolic sequence, which can locate all the cycles up to a certain topological length in the system. We also discussed the deformation of the 01 orbit when parameters undergo changes, and the evolution rule of orbital period with parameter changes was obtained.

The Chen and Lorenz systems have similar but different topological structure. They are both 3D continuous dynamical systems with two nonlinear terms, and they cannot transform one to the other through topological transformation. The two systems are dual to each other, and the Chen system has more complex topological structure and more dynamical behavior than the Lorenz system. Thus, how to classify the periodic orbits in the Chen system is an open problem. As demonstrated in the present paper, the classification method combining topology and the variational method used here could be a promising candidate for such classification.

The Lü system is a bridge between the Lorenz and Chen systems that realizes the transition from one system to the other. The above three systems can be expressed as a unified chaotic system. There are still many unknown problems in the Lü system, such as the organization of its periodic orbits and the existence of connecting orbits. As with the Lorenz flow, we can track the deformation of the orbit continuously in the Chen and Lü system by utilizing the homotopy evolution, which can help in the analysis of the existing range of the orbits under a change of parameters. Thus, a variety of bifurcation properties merit further study.