We apply a new method for the determination of periodic orbits of general dynamical systems to the Lorenz equations. The accuracy of the expectation values obtained using this approach is shown to be much larger and have better convergence properties than the more traditional approach of time averaging over a generic orbit. Finally, we discuss the relevance of the present work to the computation of unstable periodic orbits of the driven Navier–Stokes equations, which can be simulated using the lattice Boltzmann method.
Periodic orbits are interesting features of a wide variety of dynamical systems. Driven dissipative dynamical systems, exemplified by the Lorenz equations , are well known to possess attractors that are replete with unstable periodic orbits (UPOs). Indeed, these attractors may be thought of as the closure of the set of all their UPOs .
It has been pointed out  that UPOs provide a countable sequence of orbits whose limiting behaviour may be used to understand chaotic dynamics on an attractor. Cvitanović refers to them as the ‘skeleton of classical and quantum chaos’ . For certain attractors, such as those arising from the Lorenz and Duffing equations, the UPOs are reasonably well understood via symbolic dynamics and template theory .
In view of the exceptional nature of UPOs, it is remarkable that knowledge of these orbits and their properties can be used to efficiently calculate expectation values on generic orbits. This may be accomplished by means of the dynamical zeta function (DZF) formalism , whose basic ideas we will outline here.
Given an observable A that is additive along periodic orbits and has value Ap on the prime (i.e. non-composite) periodic orbit p, one definition of the DZF for discrete-time dynamical systems may be written as 1.1 where z is a complex variable, β is a real variable and the product index p ranges over all prime periodic orbits. Λp represents the stability eigenvalue of the p UPO and is the largest eigenvalue of the stability matrix J over the UPO. This matrix is computed by solving numerically the differential equation 1.2 where D is the Jacobian, and ρ the trajectory of the UPO, with an initial condition .
If we suppose that the function defined by equation (1.1) has a pole at position z=sA(β) in the complex z-plane, then the dynamical average of A over generic orbits is given by 1.3 This formalism thus reduces the problem of computing dynamical averages of a given observable to the computation of the prime periodic orbits of the corresponding dynamical system. Fortunately for us, the DZFs have very good convergence properties and including even just a few of the smaller prime orbits in our description can already furnish very accurate results .
For the driven Navier–Stokes equations (NSEs), an infinite-dimensional dynamical system, it has been proved that the attracting set is finite dimensional, with a dimension that scales as a power law in Reynolds number . This theorem thus suggests that a cleverly chosen finite-dimensional representation ought to be sufficient to capture the essential features of driven Navier–Stokes turbulence. In the present work, however, we merely discuss our method for locating UPOs in the Lorenz equations, which can be seen as a sort of ‘toy model’ for fluid turbulence.
The methodology we use to locate the UPOs is based on a space–time variational principle, which builds on previous work by Lan & Cvitanović . Our method is briefly described in Fazendeiro et al.  and, with much greater detail, in Boghosian et al. , which should be seen as a companion piece to the present work. This methodology has also been applied to the identification of UPOs in the NSEs [9,11,12], simulated using the lattice Boltzmann method , something for which petascale resources are required.
The work is organized as follows. We begin with a very brief description of the Lorenz equations and discuss the application of the variational principle to the identification of UPOs in this system. Several of the numerical hurdles involved in this work are discussed. We then show how knowledge of the system’s UPOs leads to the computation of dynamical averages via the DZF formalism. Finally, we make some assessment of the computational capabilities that are needed to numerically solve for UPOs of the Navier–Stokes equations using this technique.
2. Locating unstable periodic orbits
In this work, we apply a variational principle in order to locate UPOs in the Lorenz equations , which can be written as 2.1 In this work, we take the constants to be σ=10, R=28 and b=8/3, for which the equations are well known to exhibit chaotic orbits.
In order to simulate the gradient descent equations obtained by discretizing equation (2.1) using a variational principle, it is necessary to start with an initial condition reasonably close to a UPO. In this work, we adapt a method previously suggested by Auerbach et al.  that consists of plotting the quantity Δ(t,T), defined by 2.2 versus T, for various values of t, and an appropriate choice of norm where r(t) represents a state of the system. Note that Δ(t,0)=0. For the Lorenz system, using the l2 norm and initial condition (0.695295, 1.27032, 14.132), which lies very close to the attractor, such a contourplot is shown in figure 1 for 50≤t≤100 and 0≤T≤10. The time integration was accomplished via fourth-order Runge–Kutta with a time step of 0.001 in these units. Low values of Δ are pictured in dark blue, transitioning through light blue, green and red for the highest values. The region near the horizontal axis is coloured dark blue because Δ(t,0)=0. It is evident that there are many regions of dark blue far from the horizontal axis, corresponding to UPOs.
As is well known, the Lorenz attractor is distinguished by two prominent ‘lobes’, and each UPO traverses those lobes in a certain sequence. We denote the lobe with x,y>0 by 1, and the lobe with x,y<0 by 0, so that a traversal may be represented by a binary sequence. If we consider as equivalent two binary sequences that may be brought into correspondence by a cyclic shift, then the correspondence between UPOs and binary sequences is one-to-one. Using this terminology, an initial condition for the space–time relaxation can be generated from a symbolic description of a UPO.
To numerically compute UPOs for the Lorenz equations, we minimize the functional 2.3 where T stands for the period of an orbit (see [9,10] for more details). In this definition, the following change of variables is assumed: 2.4
To compute the value of numerically, an orbit ρ is discretized by a sequence of points, and the curve is approximated using polynomials of order seven. The conjugate gradient method is used to minimize . Within this method, a Brent solver, which does not require any knowledge of the gradient of , is used to find the minimum along a line search. In this work, iterations were ceased when a value of was reached.
This approach was used in the numerical study presented in figure 2. In this case, an analysis of the Δ(t,T) values yielded an initial condition for the UPO of the Lorenz equations, which consists of two loops in the left lobe, and one loop in the right lobe. The evolution of the functional is presented as a function of the number of iteration steps using three different algorithms. The value of seems to decrease as a power law for the Ginzburg–Landau-type minimization (see  for an example of this method) and with gradient descent, whereas it decays exponentially with the conjugate gradient method. Generally speaking, the conjugate gradient method is, therefore, numerically superior and should be preferred over the other methods.
A major obstacle to fast convergence are low-frequency disturbances, which take a long time to travel around the orbit. To avoid these, a multi-grid approach is chosen. The orbit is first discretized with a few points to which the minimization procedure is applied. The number of points is then doubled, using interpolations based on order seven polynomials, and minimization applied once more. This procedure is iterated until grid convergence is obtained, i.e. until further grid refinement does not modify the value of .
3. Dynamical zeta functions
We now describe how the DZF formalism can be used to compute the average 〈A〉 of an observable A in the Lorenz equations. This quantity can be defined as an average over a generic orbit (other than a UPO or a fixed point) with asymptotically infinite length, 3.1 where ρ(τ) is a solution to the dynamical equation. The value of 〈A〉 can also be derived from the DZF, 3.2 where the product runs over all UPOs of period smaller than or equal to Tmax, Tp is the period length of the UPO and Ap is the integral value of A over one traversal of the UPO, 3.3 where ρ′(τ) is now the solution of the dynamical equation on the UPO. As in equation (1.1), Λp represents the stability eigenvalue of the UPO. We note that equation (3.2) may not converge in all situations. However, we can use it confidently in the present case since for the classical values of the parameters σ, R, b of the Lorenz equations considered in this work, it has been rigorously proved  that the system does possess an attractor.
The dynamical average 〈A〉 is found from the relation 3.4 where the function s(β) is defined implicitly through 3.5 where s(β) is assumed to be the largest existing root.
To compute the integral of an observable along a UPO, we use a Gauss–Legendre quadrature formula to evaluate the integral between two consecutive points. This formula is exact for polynomials of order seven, and is thus consistent with the numerical model of the UPO, which also uses order seven polynomials between points.
A first attempt at enumerating all UPOs up to a certain period length consists of listing all symbolic sequences up to a fixed length. All words with two 1 or 0 characters are listed, then all words with three characters, and so on. One obvious deficiency of this approach is that the same orbits are listed several times (the orbit associated with 10 is, for example, identical to the one associated with 01). This is overcome by solving the ‘necklace problem’, i.e. by listing only sequences corresponding to distinct necklaces.
A second difficulty stems from the observation that there is no simple correlation between the number of characters in the symbolic description and the length of the period in a UPO. It is observed that for each number of characters, there are two orbits (only one is visible, as they are superimposed), which have a distinctly shorter period than the others. These extremal orbits consist of the repetition of a single character, followed by a single occurrence of the other one. It is clear that the period length of an extremal orbit of n characters is a lower bound to all orbits with n characters or more. One may, therefore, list all UPOs up to a character length k, and then use the period length of an extremal orbit of character length k+1 as a cut-off value.
Figure 3 shows the efficiency of the extremal orbits approach as a function of the number of characters. The efficiency is defined as the number of orbits left after cut-off, divided by the total number of computed orbits. Although more numerical data would be desirable, there seems to be a clear trend for the convergence of the efficiency towards the optimal value of 1/2.
Figure 4a shows how the value of 〈z〉 converges as a function of the cut-off value up to which UPOs are taken for evaluation of the DZF, using a reordering of equation (3.2), which is basically a sum over the orbits, with each one appearing only once at most. This is in clear contrast with figure 4b, where 〈z〉 has been computed, using equation (3.1), along an arbitrary orbit: the displayed results are noisy and do not exhibit any clear convergence. We should note that the averages 〈x〉 and 〈y〉 vanish for symmetry reasons because the Lorenz equations are symmetric under the transformation . The formalism of the DZF respects this symmetry, and these two averages are found to vanish exactly when they are computed from the DZF, regardless of the value of Tmax, as opposed to what happens in the more traditional time-averaging approach, exemplified by figure 4b.
We have illustrated our new variational principle by applying it to the Lorenz attractor. Several UPOs were thus found, and very accurate dynamical averages computed from them, using the DZF formalism. This can be seen as a proof of principle for the much harder problem of performing a similar calculation for the Navier–Stokes equations in a weakly turbulent regime. This is indeed a very computationally demanding task, as it requires the storage of the computational lattice in both space and time. Nevertheless, modern petascale facilities such as JUGENE, Intrepid and Ranger already allow such calculations to be performed.
We believe that, by facilitating the discovery and classification of UPOs for dynamical systems exhibiting turbulent behaviour, this methodology can allow us to make very accurate predictions of important dynamical quantities, such as topological entropy, from first principles and in a systematic fashion. This ambition can best be exemplified by comparing figure 4, in which we contrast the systematic DZF approach with the more traditional time averaging, for the Lorenz equations.
We note that the patterns exhibited in figure 1, for the Lorenz equations, appear to be anything but random. The parallelograms of dark blue regions suggest that orbits that fall near certain UPOs, with a range of periods T, can undergo transitions to different subsets of UPOs with different periods. Such a transition is evident in the vicinity of T=88. This approach suggests the possibility of a Markovian model in which orbits transition from the vicinity of certain subsets of UPOs to others. The complete description of the UPOs in the (t,T)-plane, even for the Lorenz attractor, is not yet known, and we leave that description to future work.
This work was supported by EPSRC grants EP/D051754/1 and EP/E045111/1, together with US TeraGrid NSF MRAC grants TG-DMR070013N and TG-DMR070014N. B.M.B. and J.L. were partially supported by AFOSR grant FA9550-04-1-0176 and ARO grant W911NF0410334; graphical visualization was carried out on equipment purchased by NSF MRI grant 0619447. We acknowledge a DOE INCITE allocation grant funded by the US Department of Energy, under contract DE-AC02-06CH11357, which allowed us to use Intrepid at ALCF and the PRACE Early User PRA020 grant, which provided access to JUGENE at the Jülich Supercomputing Centre in Germany.
One contribution of 26 to a Theme Issue ‘Discrete simulation of fluid dynamics: methods’.
- This journal is © 2011 The Royal Society