Equilibria, periodic orbits and computing them

TL;DR

This paper explains how to find equilibria and periodic orbits using the flow map and Jacobian-free Newton-Krylov method, providing practical code examples for complex systems like pipe flow and the Lorenz system.

Contribution

It introduces a Jacobian-free Newton-Krylov approach for computing equilibria and periodic orbits that is easy to implement with existing time stepping codes.

Findings

Method successfully applied to Lorenz system

Code is problem-independent and scalable to large systems

Enables efficient computation of complex dynamical structures

Abstract

In this short exposition, we describe equilibria and periodic orbits in terms of the flow map, {\Phi}, and discuss the essentials of the Jacobian-free Newton-Krylov (JFNK) method that can be used to find them. This method requires little more than calls to an existing time stepping code, which {\Phi} can be considered to represent. Fortran90 / MATLAB code is available to try it out for yourself, where, in the template/example the method is applied to the Lorenz system. This code is problem-independent and can be applied to large systems, having initially been developed to find periodic orbits in simulations of pipe flow.

Equations

$$\mbox{\boldmath$x$}_{t}=\mbox{\boldmath$\Phi$}^{t}(\mbox{\boldmath$x$}_{0})\,.$$

$$\mbox{\boldmath$x$}_{t^{\prime}+t}=\mbox{\boldmath$\Phi$}^{t}(\mbox{\boldmath$x$}_{t^{\prime}})\,.\qquad\mbox{{\color[rgb]{0,0,1}[flow-map]}}$$

$$\dot{X}$$

$$\dot{Y}$$

$$\dot{Z}$$

$$T(x,z,t)=\theta(x,z,t)-(z/H)\,\Delta T,\qquad u_{x}=-\frac{\partial\psi}{\partial z},\quad u_{z}=\frac{\partial\psi}{\partial x}\,,~{}~{}~{}\mbox{(stream function $\psi$)}$$

$$\psi$$

$$\theta$$

$$\mbox{\boldmath$x$}_{0}=\mbox{\boldmath$\Phi$}^{t}(\mbox{\boldmath$x$}_{0})\,,\qquad\mbox{{\color[rgb]{0,0,1}[equilibrium = fixed point]}}$$

$$\mbox{\boldmath$x$}_{p}=\mbox{\boldmath$\Phi$}^{T}(\mbox{\boldmath$x$}_{p})\,,\qquad\mbox{{\color[rgb]{0,0,1}[periodic orbit]}}$$

$$\mbox{\boldmath$x$}_{0}=g(-ct)\,\mbox{\boldmath$\Phi$}^{t}(\mbox{\boldmath$x$}_{0})\,,\qquad\mbox{{\color[rgb]{0,0,1}[relative equilibrium = travelling wave]}}$$

$$\mbox{\boldmath$x$}_{p}=g(-\bar{c}T)\,\mbox{\boldmath$\Phi$}^{T}(\mbox{\boldmath$x$}_{p})\,,\qquad\mbox{{\color[rgb]{0,0,1}[relative periodic orbit]}}$$

$$\mbox{\boldmath$x$}_{p}=\sigma\,\mbox{\boldmath$\Phi$}^{T/2}(\mbox{\boldmath$x$}_{p})\,.\qquad\mbox{{\color[rgb]{0,0,1}[pre-periodic orbit]}}$$

$$\mbox{\boldmath$x$}_{p}=\mbox{\boldmath$\Phi$}(\mbox{\boldmath$x$}_{p})\,.$$

$$\mbox{\boldmath$x$}_{0}=g(-ct)\,\mbox{\boldmath$\Phi$}^{t}(\mbox{\boldmath$x$}_{0})~{}~{}\to~{}~{}\hat{\mbox{\boldmath$x$}}_{0}=\hat{\mbox{\boldmath$\Phi$}}^{t}(\hat{\mbox{\boldmath$x$}}_{0})\,.$$

$$\mbox{\boldmath$x$}_{p}=g(-\bar{c}T)\,\mbox{\boldmath$\Phi$}^{T}(\mbox{\boldmath$x$}_{p})~{}~{}\to~{}~{}\hat{\mbox{\boldmath$x$}}_{p}=\hat{\mbox{\boldmath$\Phi$}}^{T}(\hat{\mbox{\boldmath$x$}}_{p})\,,$$

$$(5-\sqrt{5})/8\to (5+\sqrt{5})/8\to (5-\sqrt{5})/8\to ...$$

$$\frac{||\mbox{\boldmath$x$}_{t}-\mbox{\boldmath$x$}_{t-\Delta t}||}{||\mbox{\boldmath$x$}_{t-\Delta t}||}\,.$$

$$\min(||\mbox{\boldmath$x$}_{t}-\mbox{\boldmath$x$}_{t-\Delta t}||,~{}||\mbox{\boldmath$x$}_{t}-\sigma\,\mbox{\boldmath$x$}_{t-\Delta t}||)~{}~{}~{}\mbox{or}~{}~{}~{}\min_{l}(||\mbox{\boldmath$x$}_{t}-g(-l)\mbox{\boldmath$x$}_{t-\Delta t}||)~{}.$$

$$x_{i+1}=x_i-f(x_i)/f^{\prime }(x_i).$$

$$x_{i+1}=x_i+\delta x_i\text{mbox}where{f}^{\prime }(x_i)\delta x_i=-f(x_i).$$

$$\mbox{\boldmath$F$}(\mbox{\boldmath$x$}_{p})=\mbox{\boldmath$0$}~{}~{}~{}\mbox{where}~{}~{}~{}\mbox{\boldmath$F$}(\mbox{\boldmath$x$})=\mbox{\boldmath$\Phi$}(\mbox{\boldmath$x$})-\mbox{\boldmath$x$}.$$

$$(a)~{}~{}\mbox{\boldmath$x$}_{i+1}=\mbox{\boldmath$x$}_{i}+\mbox{\boldmath$\delta x$}_{i}\,~{}~{}~{}\mbox{where}~{}~{}(b)~{}~{}\left.\frac{\mbox{\boldmath$\partial F$}}{\mbox{\boldmath$\partial x$}}\right|_{\mbox{\boldmath$x$}_{i}}\mbox{\boldmath$\delta x$}_{i}=-\mbox{\boldmath$F$}(\mbox{\boldmath$x$}_{i})\,.$$

$$\left.\frac{\mbox{\boldmath$\partial F$}}{\mbox{\boldmath$\partial x$}}\right|_{\mbox{\boldmath$x$}_{i}}~{}~{}=~{}~{}\left.\frac{\mbox{\boldmath$\partial\Phi$}}{\mbox{\boldmath$\partial x$}}\right|_{\mbox{\boldmath$x$}_{i}}-I~{}~{}=~{}~{}J-I$$

$$\mbox{\boldmath$x$}=(x_{1},\,x_{2},\,x_{3})\,,\quad\mbox{\boldmath$\Phi$}(\mbox{\boldmath$x$})=\left[\begin{array}[]{c}\Phi_{1}\\ \Phi_{2}\\ \Phi_{3}\end{array}\right],\quad J=\left[\begin{array}[]{ccc}\frac{\partial\Phi_{1}}{\partial x_{1}}&\frac{\partial\Phi_{1}}{\partial x_{2}}&\frac{\partial\Phi_{1}}{\partial x_{3}}\\[2.0pt] \frac{\partial\Phi_{2}}{\partial x_{1}}&\frac{\partial\Phi_{2}}{\partial x_{2}}&\frac{\partial\Phi_{2}}{\partial x_{3}}\\[2.0pt] \frac{\partial\Phi_{3}}{\partial x_{1}}&\frac{\partial\Phi_{3}}{\partial x_{2}}&\frac{\partial\Phi_{3}}{\partial x_{3}}\end{array}\right]\,.$$

$$A\,\mbox{\boldmath$\delta x$}=\mbox{\boldmath$b$}\,,$$

$$\left.\frac{\mbox{\boldmath$\partial F$}}{\mbox{\boldmath$\partial x$}}\right|_{\mbox{\boldmath$x$}_{i}}\,\mbox{\boldmath$\delta x$}~{}~{}\approx~{}~{}\frac{1}{\epsilon}\,(\mbox{\boldmath$F$}(\mbox{\boldmath$x$}_{i}+\epsilon\,\mbox{\boldmath$\delta x$})-\mbox{\boldmath$F$}(\mbox{\boldmath$x$}_{i}))\,.$$

$$\min_{\mbox{\boldmath$\delta x$}_{i}:\ ||\mbox{\boldmath$\delta x$}_{i}||<\delta}\ \left|\left|\left.\frac{\mbox{\boldmath$\partial F$}}{\mbox{\boldmath$\partial x$}}\right|_{\mbox{\boldmath$x$}_{i}}\mbox{\boldmath$\delta x$}_{i}+\mbox{\boldmath$F$}(\mbox{\boldmath$x$}_{i})\right|\right|\ .$$

$$\mbox{\boldmath$F$}(\mbox{\boldmath$x$},T)=\mbox{\boldmath$\Phi$}^{T}(\mbox{\boldmath$x$})-\mbox{\boldmath$x$}=\mbox{\boldmath$0$}\,,$$

$$\tilde{\mbox{\boldmath$x$}}_{i}=(\mbox{\boldmath$x$}_{i},\,T_{i})\,~{}~{}~{}\mbox{and}~{}~{}~{}\tilde{\mbox{\boldmath$b$}}=(-\mbox{\boldmath$F$}(\tilde{\mbox{\boldmath$x$}}_{i}),\,0)\,.$$

$$A\,\tilde{\mbox{\boldmath$\delta x$}}_{i}=\tilde{\mbox{\boldmath$b$}}\,,$$

$$A\,\tilde{\mbox{\boldmath$\delta x$}}=\left(\left.\frac{\mbox{\boldmath$\partial F$}}{\mbox{\boldmath$\partial\tilde{x}$}}\right|_{\tilde{\mbox{\boldmath$x$}}_{i}}\,\tilde{\mbox{\boldmath$\delta x$}},\,\langle\dot{\mbox{\boldmath$x$}_{i}}|\mbox{\boldmath$\delta x$}\rangle\right)\,.$$

Full text

**Equilibria, periodic orbits and computing them.

**

EPSRC Summer School on Modal decompositions in fluid mechanics.

DAMTP, Cambridge, 5-8 August 2019.

Ashley P. Willis,

School of Mathematics and Statistics, **University of Sheffield, U.K.

**a.p.willis@sheffield.ac.uk, openpipeflow.org.

In this short exposition, we describe equilibria and periodic orbits in terms of the flow map, $\Phi$, and discuss the essentials of the Jacobian-free Newton–Krylov (JFNK) method that can be used to find them. This method requires little more than calls to an existing time stepping code, which $\Phi$ can be considered to represent. Fortran90 / MATLAB code is available to try it out for yourself, where, in the template/example the method is applied to the Lorenz system. This code is problem-independent and can be applied to large systems, having initially been developed to find periodic orbits in simulations of pipe flow.

1 Using the flow-map

Let the point \mbox{\boldmathx}_{0} be an $n$-vector representing the state of a system. For a dynamical system, as time $t$ progresses, the point \mbox{\boldmathx}_{t} traces out a trajectory, a one-dimensional curve, in an $n$-dimensional phase space $\mathcal{M}$.

We can describe the trajectory for \mbox{\boldmathx}_{t} using the flow-map denoted \mbox{\boldmath\Phi}^{t}, which takes a point \mbox{\boldmathx}_{0} and evolves it by a time period $t$:    \mbox{\boldmath\Phi}^{t}:\mbox{\boldmathx}_{0}\to\mbox{\boldmathx}_{t}   i.e.

[TABLE]

More generally, the flow-map simply advances a point along its trajectory:

[TABLE]

1.1 Example: Lorenz’s model for convection

Lorenz (1963) derived the following system for three-time dependent amplitudes, $X(t)$, $Y(t)$ and $Z(t)$:

[TABLE]

At each instant in time, the current state \mbox{\boldmathx}=(X,Y,Z) is a point in the phase space $\mathcal{M}=\mathbb{R}^{3}$. As time progresses, \mbox{\boldmathx}_{t}=(X(t),Y(t),Z(t)) traces out a trajectory, i.e. a curve, in $\mathbb{R}^{3}$. Lorenz focussed on parameter values $r=28$, $b=8/3$, $\sigma=10$, which result in chaotic trajectories. The flow-map takes us along this trajectory, see figure 1.

We should not forget that each point \mbox{\boldmathx}_{t} in phase-space corresponds to a whole convection flow pattern! Here it corresponds to a two-dimensional flow between two flat plates a distance $H$ apart, with a temperature difference $\Delta T$ between the top and bottom plates (figure 2).

The amplitudes $X(t)$, $Y(t)$, $Z(t)$ correspond to modulated variations in the temperature and velocity fields:

[TABLE]

[TABLE]

where the $c_{i}$ are scalar constants.

1.2 Invariant solutions

Equilibria and periodic orbits are topological features of the phase space $\mathcal{M}$. Irrespective of the coordinates and measure of distance used to visualise the phase space, an equilibrium remains a point, and a periodic orbit remains a closed loop. Properties, such as their eigenvalues, are also invariant to the measure of distance used.

An equilibrium \mbox{\boldmathx}_{0} is a fixed point of the flow-map that satisfies

[TABLE]

for any time $t$. A point \mbox{\boldmathx}_{p} on a periodic orbit satisfies

[TABLE]

where $T$ is the period of the orbit. In terms of solutions of the flow-map, we can consider an equilibrium to be a special case of a periodic orbit where $T$ may be arbitrarily chosen.

If a system has a homogeneous dimension, $x$, then it can have travelling wave solutions. (See figure 3a.) In a frame moving at some phase speed $c$, the solution looks steady. Equivalently, we can keep shifting the solution so that it looks steady. Let $g(l)$ be an operator that shifts a state by a distance $l$ in the $x$ direction. Then, a travelling wave satisfies

[TABLE]

for any time $t$. The state \mbox{\boldmathx}_{0} is an equilibrium solution of the slightly modified map for any $t$, hence we also call a travelling wave a relative equilibrium. Similarly, a periodic solution that recurs up to a spatial shift,

[TABLE]

for some $\bar{c}$, we call a relative periodic orbit.

Suppose the dimension $x$ has a mirror symmetry about $x=0$. Let $\sigma$ be the flip operator: \sigma\,\mbox{\boldmathx}(x)=\mbox{\boldmathx}(-x). The second half of an orbit satisfying (1.5) might just be a reflection of the first half:

[TABLE]

(See figure 3b.) Such orbits are called pre-periodic orbits. The shortest/simplest periodic orbits of a system with discrete symmetries are typically pre-periodic. Figure 4 shows the shortest periodic orbit of the Lorenz system, where the second half of the orbit is related to the first half by the 180 degree rotation symmetry $(X,Y,Z)\to(-X,-Y,Z)$.

1.3 Poincaré sections

Let \mbox{\boldmathx}^{\prime} be a point and \mbox{\boldmatht}^{\prime} be a normal vector that together define a hypersurface $\mathcal{P}$. Crossings of $\mathcal{P}$ can be defined by times $t$ when \langle\mbox{\boldmathx}_{t}-\mbox{\boldmathx}^{\prime}|\mbox{\boldmatht}^{\prime}\rangle=0. (See figure 5a.) We might restrict to when crossings occur in a particular direction, say when the inner product goes from negative to positive.

We can now let $\Phi$ be the map that takes one point on $\mathcal{P}$ to the next crossing point on $\mathcal{P}$. If a periodic orbit has a point \mbox{\boldmathx}_{p} on $\mathcal{P}$, then it satisfies

[TABLE]

The advantage is that we no longer need to worry about the period $T$ for periodic orbits. The disadvantage is that we know nothing about what happens to the orbit off $\mathcal{P}$, and in general, not all periodic orbits cross a single $\mathcal{P}$.

1.4 Slicing

For a homogeneous spatial dimension $x$, the freedom of a pattern to appear at any location is awkward when we want to compare states. Slicing is an automatic shifting procedure that removes this degree of freedom.

Here we will discuss the simplest form of slicing – ‘Fourier’ slicing. (See Budanur et al., 2015; Willis et al., 2016) When a system has a homogeneous dimension, it is commonplace to work with a periodic domain of length $L=2\pi/\alpha$. In this case, construct \mbox{\boldmathx}^{\prime}=\mbox{\boldmathx}_{c}\cos\alpha x+\mbox{\boldmathx}_{s}\sin\alpha x, where \mbox{\boldmathx}_{c} and \mbox{\boldmathx}_{s} are arbitrary states independent of $x$; they must not both be zero.

Any state $x$ may be projected onto a plane $(a_{1},a_{2})$ via a_{1}=\langle\mbox{\boldmathx}|\mbox{\boldmathx}^{\prime}\rangle and a_{2}=\langle\mbox{\boldmathx}\,|\,g(L/4)\,\mbox{\boldmathx}^{\prime}\rangle. (See figure 6a.) In this projection, the set of shifted states \{g(l)\mbox{\boldmathx}\mbox{ for all }l\} lie on a circle centred on the origin. By shifting the states, we can rotate all points on the circle to the unique point on the circle where it crosses the positive $a_{1}$ axis. All possible shifted versions of $x$ are then mapped to the unique version \hat{\mbox{\boldmathx}}=g(-l)\,\mbox{\boldmathx}, where $l=(\theta/2\pi)L$ and $\theta$ is the polar angle to $(a_{1},a_{2})$.

This operation is a symmetry-reduction, and we say that the symmetry-reduced state \hat{\mbox{\boldmathx}}, indicated by the hat, lies on a slice. Arbitrary shifts have been eliminated, so the slice has dimension one less than that of the original system. The slice $\hat{\mathcal{M}}$ is a hypersurface within the original space of states $\mathcal{M}$.

The slice is different from a Poincaré section because the symmetry reduction can be applied to \mbox{\boldmathx}_{t} for all times $t$. We can compute sliced dynamics with trajectories \hat{\mbox{\boldmathx}}_{t} that lie within the slice. (See figure 5b.) Meanwhile, trajectories only pierce a Poincaré section.

Relative equilibria (travelling waves) are reduced to equilibria automatically:

[TABLE]

All possible shifts of a state are reduced by shifting to one particular version on the slice $\hat{\mathcal{M}}$, i.e. the travelling wave is ‘pinned’ by the shifting. (See figure 6b.) Here, \hat{\mbox{\boldmath\Phi}} is the flow-map of the symmetry reduced dynamics. Note that all travelling waves of the system are reduced to equilibria, even though they typically will have a range of different phase speeds $c$.

Similarly, a relative periodic orbit becomes a closed periodic orbit, because the start and end point are shifted to a single point on $\hat{\mathcal{M}}$. The relative periodic orbit now satisfies the simpler form:

[TABLE]

for any symmetry-reduced point \hat{\mbox{\boldmathx}}_{p} on the orbit.

2 Periodic orbits

All periodic orbits (POs) in a chaotic attractor must be unstable, otherwise the behaviour would eventually be attracted to the orbit and become periodic, not chaotic. So why is it useful to find POs when they’re all unstable!?

2.1 Why periodic orbits?

Firstly, the dynamics is always very close to a PO! A chaotic attractor is dense in POs: For any chaotic point \mbox{\boldmathx}_{0} and $0<\epsilon\ll 1$, there exists a periodic point \mbox{\boldmathx}_{p} within $\epsilon$ of \mbox{\boldmathx}_{0}. See figure 7. Also…

• •

Unlike equilibria, POs exhibit dynamics! — they capture the time-dependent dynamic processes of the system and organise the chaotic set.

• •

The chaotic dynamics tends to follow the least unstable POs.

• •

A PO is a closed loop in state space. It will appear as a closed loop irrespective of the coordinates used to visualise the state space.

• •

The shortest, most fundamental, POs provide an alphabet for symbolic dynamics.

• •

Statistical properties of a chaotic attractor can be calculated in terms of sums over the POs, using their relative stability.

2.2 Examples of periodic orbits

• •

The logistic map $x_{t+1}=r\,x_{t}\,(1-x_{t})$ is chaotic for the case $r=4$ but has the (unstable) period-$2$ orbit

[TABLE]

• •

In a remarkable methodical and computational feat, Viswanath (2003) calculated 111011 periodic orbits of the Lorenz attractor, that is all periodic orbits with itineraries of up to length 20 (number of windings around the left/right ‘wings’), with an accuracy of 14 decimal digits. A few of them are shown in figure 8. The periodic orbit $A^{14}B$ was found to be the least unstable among all orbits calculated. (It is the periodic orbit with the smallest Lyapunov/Floquet exponent.)

• •

For an $n=$ 154755-dimensional model of turbulent flow in a pipe, Willis et al. (2016) calculated periodic orbits and travelling wave solutions. Travelling waves are stationary in a moving frame, so they correspond to fixed points in the sliced phase space. Two ‘clouds’ of solutions were observed, shown in figure 9, and the relationship between them found to be the reflection symmetry. The phase-space visualisation reveals that trajectories don’t like to switch orientation with respect to the symmetry, due to the presence of a strongly repelling (unstable) fixed point, marked ‘A’, that lives between the two clouds. For beautiful examples from Couette flow, see Cvitanović and Gibson (2010).

2.3 Searching for recurrences

At present, the standard approach to finding near-recurrences is pretty crude. It involves recurrence plots, with distance measures tailored for the case at hand, and, if not automated, visual inspection of the plot for close recurrences. Nevertheless, this is the usual means to find points close to periodic orbits that can be refined to exact recurrences using the Newton method (described in the next section).

For a recurrence plot we compute something like

[TABLE]

Figure 10 is an example for pipe flow.

We then look for local minima in the plot that provide candidate recurrent points, \mbox{\boldmathx}_{p}\approx\mbox{\boldmathx}_{t-\Delta t} and $T\approx\Delta t$. The normalisation factor might be chosen to depend on both ||\mbox{\boldmathx}_{t}|| and ||\mbox{\boldmathx}_{t-\Delta t}||, or might not be necessary at all. A complication is that we might need to minimise over discrete symmetries, such as the flip operator $\sigma$, or over shifts (section 1.2) e.g.

[TABLE]

Minimisation over shifts can be avoided if slicing is applied (section 1.4).

The norm itself might need tinkering with. For example, in the sheared flow of fluid, perturbations in the streamwise dimension are typically an order of magnitude larger than the crossflow components. A ‘compensatory’ norm helps in this case, where the components are scaled to be more similar in magnitude.

3 The Newton–Krylov method

The Jacobian-free Newton–Krylov (JFNK) method is a variant of the Newton–Raphson method. In its raw form, the Newton–Raphson method for an $n$-dimensional system involves an $n\times n$ Jacobian matrix, which can be tricky to evaluate. It is possible to avoid this evaluation using a Krylov-subspace method (Knoll and Keyes, 2004).

3.1 The Newton–Raphson method

To find roots $x$ such that $f(x)=0$ in one dimension, given an initial guess $x_{0}$, the Newton-Raphson method generates improvements using the iteration

[TABLE]

Re-arranging, we may re-express the iteration as

[TABLE]

Our task is to find fixed points of the map such that \mbox{\boldmathx_{p}}=\mbox{\boldmath\Phi}(\mbox{\boldmathx}_{p}), i.e.

[TABLE]

(The fixed points could correspond to equilibria, periodic orbits, or their relative equivalents. Augmentations, if necessary, to find a period $T$ or spatial shift $l$ are delayed to section 3.4.) The extension of Newton’s method (3.2) to an $n$-dimensional system is then

[TABLE]

In order to apply the update (3.4a), the linear system (3.4b) needs to be solved for the unknown \mbox{\boldmath\delta x}_{i}.

In (3.4b), the matrix part is given by

[TABLE]

where $J$ is the Jacobian matrix for \mbox{\boldmath\Phi}(\mbox{\boldmathx}) and $I$ is the identity matrix. For the case $n=3$,

[TABLE]

3.2 Jacobian-Free method

The $n\times n$ Jacobian matrix $J$ is usually difficult to evaluate. We might not even have sufficient computer memory to store it for a high dimensional system. The problem (3.4b), however, is in the form

[TABLE]

where $A$ is an $n\times n$ matrix and $\delta x$ and $b$ are $n$-vectors. This can be solved for $\delta x$ using the Krylov-subspace method GMRES(m). The GMRES algorithm does not need to know the matrix $A$ itself, only the result of multiplying a given vector by $A$. The method seeks a solution for $\delta x$ in \mathrm{span}\{\mbox{\boldmathK}_{1},\,\mbox{\boldmathK}_{2},\dots,\,\mbox{\boldmathK_{m}}\}, i.e. \mbox{\boldmath\delta x}=c_{1}\,\mbox{\boldmathK}_{1}+c_{2}\,\mbox{\boldmathK}_{2}+...+c_{m}\,\mbox{\boldmathK}_{m}. It is common to start with \mbox{\boldmathK}_{1}=\mbox{\boldmathb}/||\mbox{\boldmathb}||. The next vector is generated by evaluating \tilde{\mbox{\boldmathK}}_{i+1}=A\,\mbox{\boldmathK}_{i}, then \mbox{\boldmathK}_{i+1} is obtained by orthonormalising \tilde{\mbox{\boldmathK}}_{i+1} against the previous \mbox{\boldmathK}_{j} ($j\leq i$) using the Gram-Schmidt method. Next, \mathrm{error}=||A~{}\mbox{\boldmath\delta x}-\mbox{\boldmathb}\,|| is minimised over the coefficients $c_{j}$ ($j\leq i+1$) and the process repeated if $\mathrm{error}$ is too large.

Iterations of the GMRES algorithm for the problem (3.4b) involve calculating matrix-vector products with given $\delta x$ that may be approximated:

[TABLE]

$\epsilon$ is a small scalar value; a typical value is $\epsilon$ such that (\epsilon||\mbox{\boldmath\delta x}||)\,/\,||\mbox{\boldmathx}_{i}||=10^{-6}. The important point is that we do not need to know the Jacobian — only a routine for evaluating \mbox{\boldmathF}(\mbox{\boldmathx}) is required.

Note that provided that each step of the Newton method, $\delta x$, takes $x$ in approximately the correct direction, the method is expected to converge. Therefore the tolerance specified in the accuracy of the solution for $\delta x$ in each Newton step (calculated via the GMRES method) typically need not be so stringent as the tolerance placed on the Newton method itself for the solution $x$. For example, we might seek a relative error for the Newton solution ||\mbox{\boldmathF}(\mbox{\boldmathx})||/||\mbox{\boldmathx}||=O(10^{-8}), but a relative error for the GMRES solution ||A\,\mbox{\boldmath\delta x}-\mbox{\boldmathb}||/||\mbox{\boldmath\delta x}||=O(10^{-3}) is likely to be sufficient for calculation of the steps $\delta x$.

3.3 Hookstep approach

To improve the domain of convergence of the Newton method, it is commonplace to limit the size of the step taken. One approach is simply to take a ‘damped’ step in the direction of the solution to 3.4(b), i.e. step by \alpha\ \mbox{\boldmath\delta x}_{i}, where $\alpha\in(0,1]$. In the ”’hookstep approach”’, we minimise subject to the condition that the magnitude of the Newton step is limited, ||\mbox{\boldmath\delta x}_{i}||<\delta, where $\delta$ is the size of the ”’trust region”’:

[TABLE]

Given the minimisation, the hookstep \mbox{\boldmath\delta x}_{i} is expected to produce a better result than a simple damped step of the same size. It is also expected to perform much better in ’valleys’, where it produces a bent/hooked step to a point along the valley, rather than jumping from one side of the valley to the other; see figure 11.

The hookstep can be calculated with little extra work to the GMRES method, provided that the size of Krylov-subspace, m, is chosen sufficiently large to solve to the desired accuracy within m GMRES iterations; for details see Viswanath (2007) [particularly v1 on arxiv.org].

For a given \mbox{\boldmath\delta x}_{i}, the reduction in error predicted by the linearisation (3.9) can be compared with the actual reduction in ||\mbox{\boldmathF}(\mbox{\boldmathx}_{i}+\mbox{\boldmath\delta x}_{i})||. According to the accuracy of the prediction, the size of the trust region $\delta$ can be adjusted automatically; see Dennis and Schnabel (1996).

3.4 Adding constraints

3.4.1 Time constraint

When looking for a periodic orbit, the period $T$ is an extra unknown. One way to eliminate needing to find $T$ is to work within a Poincaré section, as described in section 1.3. We can then attempt to solve the function \mbox{\boldmathF}(\mbox{\boldmathx})=\mbox{\boldmath\Phi}(\mbox{\boldmathx})-\mbox{\boldmathx} as it stands to find a point \mbox{\boldmathx}_{p} on the Poincaré section that corresponds to a periodic orbit.

We might not want to restrict ourselves to Poincaré sections. We must then solve

[TABLE]

for (\mbox{\boldmathx},T). We augment the whole system. Let

[TABLE]

We now want to solve a system of the form

[TABLE]

for \tilde{\mbox{\boldmath\delta x}}_{i}=(\mbox{\boldmath\delta x}_{i},\,\delta T_{i}), but need an extra constraint because we have an extra unknown. We choose that the update \mbox{\boldmath\delta x}_{i} has no component that points along the trajectory, i.e. \langle\dot{\mbox{\boldmathx}_{i}}|\mbox{\boldmath\delta x}_{i}\rangle=0. Following the ethos of matrix-free methods, that we do not need to know the matrix $A$ itself, we only need to state the result of multiplication by $A$:

[TABLE]

We use the approximation (3.8) to evaluate the first part of the result. The augmented system (3.12) can be solved for \tilde{\mbox{\boldmath\delta x}}_{i} using the GMRES algorithm by applying multiplications (3.13). The update for both the state and the period is then \tilde{\mbox{\boldmathx}}_{i+1}=\tilde{\mbox{\boldmathx}}_{i}+\tilde{\mbox{\boldmath\delta x}}_{i}.

3.4.2 Shift constraints

For relative equilibria (travelling waves) and relative periodic periodic orbits we need to solve

[TABLE]

where $l$ is an unknown spatial shift in the homogeneous $x$-dimension.

One way to avoid the extra unknown is to work with the sliced dynamics (section 1.4) so that arbitrary shifts are automatically eliminated. Alternatively, we can augment the system again. Let

[TABLE]

We now want to solve a system of the form

[TABLE]

for \tilde{\mbox{\boldmath\delta x}}_{i}, but need another constraint to match the extra unknown. This time we choose that the update \mbox{\boldmath\delta x}_{i} has no component that just corresponds to a spatial shift, i.e. \langle\partial_{x}\mbox{\boldmathx}_{i}|\mbox{\boldmath\delta x}_{i}\rangle=0. We assert that multiplication by $A$ is:

[TABLE]

3.5 Preconditioning

The good news is that you can ignore preconditioning and skip this section if you are combining Newton–Krylov with timestepping to evaluate the flow-map \mbox{\boldmath\Phi}^{T}(\mbox{\boldmathx}).

3.5.1 Exponentiation and timestepping

The GMRES algorithm is closely related to another Krylov-subspace method, the Arnoldi method, which is used to calculated eigenvalues of a matrix $A$. It tends to find the eigenvalues most separated in the complex plane first, but those might be of little interest. For example, the Laplacian $\nabla^{2}$ has a spectrum of very negative eigenvalues corresponding to high frequency oscillations that rapidly decay. Basically, we do not wish to build a Kyrlov-subspace involving such modes.

It may be better to work with $\tilde{A}=\mathrm{e}^{A}=1+A+\frac{1}{2!}A^{2}+...$, corresponding to the eigenproblem \mathrm{e}^{\sigma}\mbox{\boldmathx}=\mathrm{e}^{A}\mbox{\boldmathx}. This problem shares the same eigenvectors as the problem \sigma\,\mbox{\boldmathx}=A\,\mbox{\boldmathx}, but has more suitable eigenvalues, $\tilde{\sigma}=\mathrm{e}^{\sigma}$. The negative eigenvalues $\sigma$ then correspond to eigenvalues $\tilde{\sigma}$ bunched close to the origin. The Arnoldi method then favours the $\tilde{\sigma}$ most distant from the origin, corresponding to the $\sigma$ with largest real parts.

Note that for the system \partial_{t}\mbox{\boldmathx}=A\,\mbox{\boldmathx}, time integration corresponds to exponentiation: Taking eigenvector $x$ with growth rate $\sigma$ as an initial condition, the result of time integration from [math] to $T$ is \mathrm{e}^{\sigma T}\mbox{\boldmathx}. We therefore have that \mathrm{e}^{\sigma T}\mbox{\boldmathx}=\mbox{\boldmathx}+\int_{0}^{T}A\,\mbox{\boldmathx}\,dt=\mathrm{e}^{AT}\mbox{\boldmathx}, which can be written \tilde{\sigma}\,\mbox{\boldmathx}=B\,\mbox{\boldmathx} where $\tilde{\sigma}=\mathrm{e}^{\sigma T}$ is the eigenvalue of the time integration operator $B=\mathrm{e}^{AT}$.

3.5.2 Explicit preconditioning

GMRES is likely to find it easier to solve M^{-1}A\,\mbox{\boldmathx}=M^{-1}\mbox{\boldmathb} than the original system, if $M^{-1}$ is an approximate inverse for $A$. For example, if $A$ is dominated by its diagonal elements, we might take $M$ to be the banded matrix consisting of the diagonal and the first sub- and super-diagonals of $A$. Each GMRES iteration applied to the modified system now requires a muliplication by $A$ then by $M^{-1}$. This is fine, as, for a banded matrix, it is quick and easy to solve M\mbox{\boldmathx}^{\prime}=\mbox{\boldmathx} for \mbox{\boldmathx}^{\prime}. Like $A$, we don’t need to know the matrix $M^{-1}$ itself, only the result of multiplication by the matrix.

4 Try it yourself!

Application of the Newton–Krylov method to the Lorenz system

Given that the Newton–Krylov method is designed to cope with high-dimensional systems (the same code has been used to find travelling waves in pipe flow), this is somewhat overkill, but it helps illustrates how we can use the solver as a black box…

**Please cite openpipeflow.org (Willis, 2017) if you use this code in your research. Thanks!

**

• •

Download the Template/Example (Fortran90 / MATLAB / Octave)

http://www.openpipeflow.org/index.php?title=Newton-Krylov_method

• •

For MATLAB, the unpacked tgz/zip file has separate .m files for each function.

Take a look at

• –

Lorenz_f.m: Lorenz evolution rule \dot{\mbox{\boldmathx}}=\mbox{\boldmathf}(\mbox{\boldmathx}).

• –

steporbit.m: Evaluate \mbox{\boldmath\Phi}^{T}(\mbox{\boldmathx}), i.e. step $f$ by ndts_ timesteps, where the input x(1)$=T$, and x(2:4)$=(X,Y,Z)$. The timestep size is dt$=T/$ndts_.

• –

saveorbit.m: Output at end of each Newton iteration. relative_err\,=||\mbox{\boldmathF}(\mbox{\boldmathx})||\,/\,||\mbox{\boldmathx}||.

• –

MAIN.m: Set up initial guess \mbox{\boldmathx}_{0} and call the black box NewtonHook.m.

• •

Other functions are called by NewtonHook.m, and are unlikely to need changing for a problem of this type, where shifts and other spatial symmetries are ignored:

• –

getrhs.m: Evaluate right-hand side \tilde{\mbox{\boldmathb}} (3.11) i.e. \mbox{\boldmathF}(\tilde{\mbox{\boldmathx}})=\mbox{\boldmath\Phi}^{T}(\mbox{\boldmathx})-\mbox{\boldmathx}.

• –

multJ.m : Evaluate multiplication (3.13), i.e. multiplication by the Jacobian.

• –

multJp.m: Preconditioner for multiplication (here an empty function).

• –

dotprd.m: Evaluate inner product \langle\mbox{\boldmatha}|\mbox{\boldmathb}\rangle.

• –

GMRESm.m: Method of section 3.2.

• –

GMREShook.m: Calculate hookstep, section 3.3.

• •

The following data are points on the periodic orbits of figure 8, taken from Viswanath (2003). $Z=27$ in call cases.

[TABLE]

• •

In MATLAB, call MAIN. It will plot the result of timestepping the initial guess for the AB orbit (green), call the NewtonHook subroutine, then plot the converged solution (blue). Scroll back through the output, and compare relative_err for the initial guess at iteration 0 with the final relative error.

• •

Comment/uncomment other initial guesses new_x\,=\mbox{\boldmathx}_{0}, or experiment with your own. How do they affect the number of Newton iterations taken?

[Typically convergence takes $O(10)$ iterations, otherwise it will never converge.]

• •

Uncomment the initial guess for an equilibrium. Here we assume a short fixed $T$, too short for a PO; $T$ is not permitted to change, otherwise ||\mbox{\boldmath\Phi}^{T}(\mbox{\boldmathx})-\mbox{\boldmathx}|| could be reduced by simply taking $T\to 0$. Check that MAIN can find the analytic equilibrium solution $(\pm\alpha,\pm\alpha,r-1)$, where $\alpha=\sqrt{(r-1)\,b}$.

4.1 Adapting the code for your own use

• •

For a very large system, for which you might consider parallelization (see final comment), you should probably use the Fortran90 version.

• •

Experiment with the Template/Example first, to get used to how the code is set up. The initial guess is put in new_x.

• •

Note that at present, new_x(1)$\,=T$ (the period), and new_x(2:end)\,=\mbox{\boldmathx} (the state).

• •

The place to start is then steporbit. If you already have an existing timestepping code, it could do something as simple as call it externally via system calls:

function y = steporbit(ndts_,x) persistent dt

if ndts_ ~= 1 % Set timestep size dt=T/ndts_ dt = x(1) / ndts_ ; % If only doing one step to calc \dot{x}, end % then use previously set dt.

a = x(2:end) ;

WRITE DATA TO FILES: dt timestep size ndts_ number of steps to take a initial condition

LOAD STATE, TIMESTEP, SAVE STATE: system(’run_my_code.exe’)

LOAD TIMESTEPPED STATE: --> a

y = zeros(size(x)) ; y(2:end) = a ; end

• •

saveorbit is called at the end of each Newton iteration. Add code here to save the current state new_x.

• •

If your inner product corresponds to \langle\mbox{\boldmatha}|\mbox{\boldmathb}\rangle\,=\,\mbox{\boldmatha}^{T}W\mbox{\boldmathb} where $W$ is a diagonal matrix of positive weights, and here $T$ is the transpose, then pass \mbox{\boldmathx}^{\prime}=W^{\frac{1}{2}}\mbox{\boldmathx} to the code. The existing functions that take inner products then need no modification.

• •

For parallel use with MPI+Fortran, the NewtonHook and GMRES codes do not need changing: Split vectors over threads and let each thread pass its section to NewtonHook. The only place where an MPI call is required is an MPI_Allreduce in the dotprod function. To avoid all threads outputting information, set info=1 on rank 0, and info=0 on all other ranks.

• •

Further information at openpipeflow.org.

Acknowledgements

AW would like to thank Rich Kerswell, Predrag Cvitanović (chaosbook.org), John Gibson (channelflow.org), Marc Avila and many others for their generous support in many forms. Developed under EPSRC grants EP/K03636X/1, EP/P000959/1.