## Abstract

We present new methods for the determination of periodic orbits of general dynamical systems. Iterative algorithms for finding solutions by these methods, for both the exact continuum case, and for approximate discrete representations suitable for numerical implementation, are discussed. Finally, we describe our approach to the computation of unstable periodic orbits of the driven Navier–Stokes equations, simulated using the lattice Boltzmann equation.

## 1. Introduction

The study of turbulence can be described as one of the most important fundamental problems still faced by current research, sometimes hailed as the last great unsolved problem from classical mechanics [1]. The computational complexity of this problem scales with a power of approximately 9/4 of the number of degrees of freedom, and typical flows that can be seen in everyday life are still out of the reach of even present-day petascale computers. One fruitful way to look at this problem thus consists of aiming to greatly reduce the number of degrees of freedom involved, by identifying some reproducible flows that can represent the actual turbulent state.

Driven dissipative dynamical systems 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 [2]. In many respects, these UPOs are exceptional orbits. The natural measure associated with the dynamics on the attractor may be dramatically different on UPOs than it is on generic orbits, even though the latter are accumulation points of sequences of the former. Indeed, Grebogi *et al*. [3] point out that UPOs can be labelled as both ‘hot spots’ and ‘cold spots’, depending on whether the regions they occupy on the attractor are visited more or less frequently by generic orbits.

The knowledge of UPOs and their properties can be used to efficiently calculate expectation values on generic orbits, by means of the dynamical zeta function (DZF) formalism, an excellent introduction to the subject of which may be found in the book by Cvitanović *et al*. [4]. Even for the driven Navier–Stokes equations, 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 [5]. Though this theorem suggests that a cleverly chosen finite-dimensional representation ought to be sufficient to capture the essential features of driven Navier–Stokes turbulence, this observation has not yet been fully exploited in computer simulation. However, numerical evidence for the presence of UPOs on the attractor for the driven Navier–Stokes equations has been presented by Kawahara & Kida [6] and, more recently, by Gibson *et al*. [7].

The problem of identifying periodic orbits is thus applicable to a wide range of ordinary and partial differential equations, including many of interest in mathematical physics. The values of the dependent variables of the system shall collectively be called the *state* of the system, and the set of all possible states of the system shall be called the *state space*, *Ω*. For example, for the Lorenz equations, for the gravitational three-body problem (three coordinates plus three momenta for each of three particles), and *Ω* is the set of all divergence-free vector fields for the incompressible Navier–Stokes equations of fluid dynamics. Periodic orbits correspond to closed trajectories in *Ω*.

Lan & Cvitanović [8] made important progress on the problem of finding UPOs by introducing a variational method for locating periodic orbits for general flows, , where the dot denotes differentiation with respect to time *t*. The method begins with an initial guess for the entire closed trajectory, with *s*∈[0,2*π*), where *s* is a parameter, and computes the local tangent vector at all points along it. It then seeks to alter so that the vector points in the same direction as the vector ** f**(

**) to the greatest extent possible. This is done by minimizing the quantity 1.1 where**

*r**λ*(

*s*) is a function of parameter

*s*, expressing the proportionality between the two vector fields at the point on the trajectory with parameter

*s*. Once

*F*

^{2}is made to vanish, the function

*λ*(

*s*) is equal to d

*t*/d

*s*, where

*t*is time. The method is remarkably effective; Lan & Cvitanović [8] demonstrate its application to the Hénon–Heiles system and the Kurimoto–Sivashinsky equation for flame-front evolution.

However, unless one begins close to a state on a closed orbit, the above procedure will not be likely to converge owing to the inherent instability of the UPOs. For this reason, attention has been given to methods for locating nearly periodic orbits. One remarkably effective method, introduced by Auerbach *et al*. [9], is to plot the quantity Δ(*t*,*T*) defined by
1.2
versus *T*, for various values of *t*, and some appropriate choice of norm. Note that Δ(*t*,0)=0. If a time *T*>0 is found at which Δ(*t*,*T*) is both at a local minimum and very small in magnitude, then the orbit is nearly periodic. In that case, an orbit with initial condition {** r**(

*s*)|

*t*−

*T*≤

*s*<

*t*} might be an appropriate place to begin the iterative procedure described above. So

*et al*. [10] pointed out that nearly periodic orbits are even more evident in a contour plot of Δ(

*t*,

*T*), and we adopt that procedure in this study.

In this work, we discuss a new variational principle for locating periodic orbits. Like that of Lan & Cvitanović [8], it may be applied to any differential equation that is first order in time.^{1} Rather than insist on a parameter that varies between [0,2*π*) and introduce a separate function that needs to be varied to get the parametrization right, we formulate the variational principle for ** r**(

*t*) directly, where

*t*∈[0,

*T*). Since we have no

*a priori*way to know what the orbit period

*T*is, we leave it as another quantity that must be varied; that is, we define a quantity that is both a functional of

**(**

*r**t*) and a function of

*T*, and we demand that its variation with respect to both of these dependencies vanish.

We begin with a description of the new variational principle, and discuss its application to the Lattice Boltzmann (LB) [11] discretization of the Navier–Stokes equations of viscous fluid dynamics. Because the method is iterative and requires an initial guess for the periodic orbit, we adapt the above-described technique for locating nearly periodic orbits. 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. Variational principle

We seek a periodic solution, ** r**(

*t*), of 2.1 where

**∈**

*r**Ω*may contain many degrees of freedom,

**is a vector field on the set of all possible states of the system,**

*f**Ω*, the dot denotes differentiation with respect to

*t*and we denote the unknown period by

*T*. We may think of equation (2.1) as one differential equation and one algebraic equation for one unknown function,

**(**

*r**t*), and one unknown scalar,

*T*, both of which we seek to determine.

It is convenient to make the following change of variables:
2.2
where ** ρ** is periodic with period one. This in turn yields the equations
2.3
where the prime denotes differentiation with respect to

*χ*.

Our approach to finding solutions to the above equations aims to minimize the functional
2.4
It is manifest that . If we restrict attention to differentiable functions ** r**(

*t*), it is also clear that only for solutions of equation (2.1). In fact, we shall henceforth assume that

**is twice differentiable. Though the above functional seems similar to equation (1.1), we plan to vary this functional with respect to both**

*r***(**

*r**t*) and the unknown period

*T*. By the fundamental theorem of calculus, we have 2.5 so the derivative of with respect to

*T*must vanish as long as the equation of motion continues to hold at

*T*.

To use the above functional, we first use equation (2.2) to change variables in equation (2.4), rewriting the functional in terms of ** ρ** and

*T*, 2.6 where the prime denotes differentiation with respect to

*χ*. The problem of finding periodic solutions to equation (2.1) is thus reduced to that of minimizing the functional , where

**may range over the space of twice differentiable functions with unit period, and**

*ρ**T*is a positive real number. At any such minimum, the Fréchet derivative of with respect to

**(**

*ρ**χ*), and the partial derivative of with respect to

*T*must both vanish. Following this methodology, after some straightforward calculations, we obtain equations for both the period 2.7 and the second-order, nonlinear, integrodifferential boundary-value problem 2.8 with periodic boundary conditions,

**(0)=**

*ρ***(1) and**

*ρ**T*given by equation (2.7). Remarkably, this boundary-value problem must be satisfied along any periodic orbit of equation (2.1).

Based on this approach, an algorithm for solving equation (2.3) will then consist of finding minima of at which . This could be done, e.g. by gradient descent, using the Fréchet derivatives calculated above, or by line searches along the directions thereby defined.

To use the gradient-descent method, one could try to achieve a finite-dimensional representation by approximating the derivatives with respect to *χ* by finite differences. As a rule, however, when differential equations are derivable from variational principles, it is decidedly preferable to instead discretize the functional, and then vary that with respect to each of the finite set of variables used. Adopting this philosophy, we represent ** ρ**(

*χ*,

*s*) by a set of

*N*discrete points,

*ρ*_{0}(

*s*),

*ρ*_{1}(

*s*),…,

*ρ*_{N−1}(

*s*). The functional in equation (2.6) then becomes a

*function*of the

*ρ*_{j}(

*s*), 2.9 Here, we have evaluated

**at the midpoint between**

*f*

*ρ*_{j}and

*ρ*_{j+1}. Note that this is not something that could easily be done in a forward solver for an initial-value problem because it would be a nonlinear implicit problem in that context. Here, however, it does not present any additional difficulty.

The analogues of the Fréchet derivatives for this case are then the partial derivatives
2.10
and
2.11
where we have used, e.g. *ρ*_{k+1/2} to denote (*ρ*_{k}+*ρ*_{k+1})/2.

The corresponding discrete gradient descent equations are
2.12
and
2.13
In practice, these equations may be evolved in *s* using nothing more sophisticated than an Euler stepping scheme. Using this variational principle, we have shown [12], for the case of the Lorenz equations, that UPOs can be found systematically and very accurate averages computed using the DZF formalism.

## 3. The lattice Boltzmann equation

The space–time variational approach presented in §2 has been applied to the LB equation [11,13,14] in the simulation of nearly incompressible three-dimensional weakly turbulent fluid flow, driven by an Arnold–Beltrami–Childress (ABC) force [15–17]. For this work, we used a Lattice–Bhatnagar–Gross–Krook (LBGK) model whose main equation can be written as
3.1
In equation (3.1), *f*_{i}(** r**,

*t*) is the distribution function of the particle located at site

**at time**

*r**t*with velocity index

*i*=1,…,

*Q*,

*c*_{i}represents the (discrete) velocities of each pseudo-particle, is the equilibrium distribution function and

*τ*is the single relaxation time. The equations resulting from the application of the space–time variational approach to the LBGK model can be found in the work of Fazendeiro

*et al*. [15]. Since for this case, the period of the UPOs is an integer (owing to LB being a discrete-time model) we chose to vary instead the inverse of the relaxation time,

*ω*≡1/

*τ*.

In figure 1, we show a Δ(*t*,*T*) plot, using the *f*_{i}(** r**,

*t*) distribution functions of the LBGK D3Q19 model as the dynamical variables, in an ABC flow with Reynolds number equal to 371, after the transients have died down. Interestingly, this shows some similarity, in terms of minima distribution, with the patterns seen on similar plots for the Lorenz equations [12]. The lowest (non-trivial) minima we have found so far for the weakly turbulent flow, are in the region of

*T*>10

^{4}, as seen in figure 1. This imposes very substantial computational requirements, since we require not only all of the time slices of the system to be loaded into the memory of the computer but also extra arrays, with the same dimensions, for the minimization algorithm.

Several of the minima identified in Δ(*t*,*T*) plots similar to figure 1 have been ported to petascale resources. A gradient descent algorithm is then applied for the minimization of a suitably defined functional [15]. Extensive work performed using the conjugate gradient algorithm failed to improve the convergence rate owing to the high instability of these orbits and the extremely high number of degrees of freedom involved.

We are currently developing a slightly different methodology that uses orbit segments, as opposed to invoking all of the LB time slices of the system, and we expect this to greatly improve the convergence rate for these very large space–time orbits. Related to this, we are also implementing a Newton–Krylov hook-step solver, as suggested by the work of Viswanath [18] on plane Couette turbulence, which is expected to converge much faster than gradient descent, while being numerically more stable than the conjugate gradient minimization method.

## 4. Conclusions

We have built upon recent work of Lan & Cvitanović [8] to present a new variational principle for finding UPOs of deterministic dynamical systems. The variational principle presented in this work is distinguished by not requiring a scaling function (*λ*(*s*) in equation (1.1)) to capture the invariance of the periodic orbit under redefinitions of the independent variable, but rather treating the unknown period *T* as one of the quantities to be varied.

The discovery and classification of UPOs for the Navier–Stokes equations using this methodology are very computationally demanding tasks, as they require 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.

Another important aspect raised by this research is the digital curation of the UPOs for turbulent fluid flow. Owing to the sheer size of these UPOs, as well as their inherent instability, new computational methods will be required to effectively label and store them in a way that makes them accessible to the scientific community. The main advantage of storing UPOs to represent a turbulent flow is that it needs to be done only once. After that, the turbulent average of any given quantity can be computed directly from the UPO library with high accuracy and without the need to solve an initial-value problem again. This methodology thus has great potential to become a new paradigm in the study of large driven dissipative dynamical systems and brings us closer to solving one of the main challenges still remaining in classical physics, that of fluid turbulence.

## Acknowledgements

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. was 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. Discussions with Jennifer Tam are also acknowledged.

## Footnotes

↵1 This restriction is not problematic, since one can always introduce additional dependent variables to recast the equations as first order in time. For example, Newton’s equations of motion are second order in time, but they can be recast as the first-order Hamilton’s equations by letting the state of the system include momenta as well as spatial positions.

One contribution of 26 to a Theme Issue ‘Discrete simulation of fluid dynamics: methods’.

- This journal is © 2011 The Royal Society