A reduced variational approach for searching cycles in high-dimensional systems

In practice, not all directions are equally important for orbit adjustment when approaching a periodic orbit from a guess loop. On a hyperplane perpendicular to the periodic orbit, the vicinity of an orbit point is stretched in a few directions and compressed in others when moving along the orbit. The correction in the compressed directions is automatically made during the orbit evolution and the deviation in the stretching or neutral directions has to be corrected by the variational scheme. Therefore, it is essential to find these important directions and rewrite the original variational equation in reduced coordinate frames.

Assuming that along the important directions, we already built a family of local coordinate frames $\{\bm{e}_{k}(s_{i})\}$ $(k=1,2,...,n;~{}i=1,2,...,N$ and $n\leq d)$ for the projection ( Fig.1), the reduced variational equation may be derived as follows.

Firstly, we write an approximation of Eq.(II) in the subspace $\{\bm{e}_{k}\}$,

and a dot-product of both sides of Eq.(3) with $\bm{e}_{k}(s_{i})$ results in

By feeding Eq. (9) into Eq. (10), we obtain the projection of the first term in Eq.(10),

where $\Delta\tilde{x}^{j}_{i}$ and $\bm{e}^{j}_{i}$ labels $\Delta\tilde{x}_{j}(s_{i})$ and $\bm{e}_{j}(s_{i})$ respectively. According to Eq. (11), every constant block in the differential operator $\hat{\bm{D}}$ in Eq. (7) is replaced by an $n\times n$ matrix which is defined by $\bm{e}^{k}_{i}\cdot\bm{e}^{j}_{i^{\prime}}$, where $(k,j)$ marks the block position and $(i,i^{\prime})$ the position within the $(k,j)-$block. The second term $\frac{\partial\bm{v}}{\partial\bm{x}}\frac{\partial\tilde{\bm{x}}}{\partial\tau}$ in Eq.(10) can be written as

where $k\,,j$ marks the new coordinates and $\{\hat{\bm{e}}_{p}\}$ make up the original Cartesian coordinate system which are constant at every point in the phase space. $v_{k}=\bm{\bm{v}}\cdot\bm{e}_{k}(s_{i})$ ($\{\bm{e}_{k}\}$ is an orthogonal basis), $e_{kp}(s_{i})=\bm{e}_{k}(s_{i})\cdot\hat{\bm{e}}_{p}$, and $\Delta\tilde{\bm{x}}(s_{i})=\sum\limits_{j}\Delta x_{ij}\bm{e}_{j}(s_{i})=\Delta\tilde{\bm{x}}_{i}$. Therefore, the size of the velocity gradient matrix $\frac{\partial\bm{v}_{k}}{\partial\bm{x}_{j}}$, which is obtained by numerical differentiation, is much smaller than that of the original velocity gradient $A=\partial v/\partial x$ of Eq.(3) if $n\ll d$. About other terms in Eq.(10), the vectors can be just projected into the local subspace $\{\bm{e}_{k}\}$ and the matrix $\tiny{\begin{pmatrix}\vec{\bm{A}}&~{}-\vec{\bm{v}}\\ \vec{\bm{a}}&~{}~{}~{}0\end{pmatrix}}$ in Eq. (8) is then reduced to $(Nn+1)\times(Nn+1)$. Because $n\ll d$ actually holds in many nonlinear systems, the reduced matrix order $Nn+1$ will be much smaller than the original one $Nd+1$. Therefore the storage and computing time required to deal with the reduced form of Eq. (8) are greatly cut down. The larger the dimension $d$, the more prominent the benefits of this reduction.

In this section, we will discuss how to build the local coordinate system $\{\bm{e}_{k}\}$ mentioned above at each lattice point $s_{i}$. However, before that, Let’s first explain the concept of pseudo-evolution. Two adjacent points on the loop are approximately related by the evolution $\bm{x}(s_{i})\longrightarrow\bm{x}(s_{i+1})$ with Eq. (1) if the guess is reasonably good. Similarly, the corresponding nearby points are assumed evolving as $\bm{x}(s_{i})+\bm{h}(s_{i})\longrightarrow\bm{x}(s_{i+1})+\bm{h}(s_{i+1})$ governed by Eq.(1) either, where $\bm{h}(s_{i})$ and $\bm{h}(s_{i+1})$ represent respectively tiny offset vectors at $\bm{x}(s_{i})$ and $\bm{x}(s_{i+1})$. Then we have approximately

where $\Delta t$ is the time interval between points $s_{i}$ and $s_{i+1}$. Eq. (13) is the pseudo-evolution equation of points near the guess loop, which obviously is rough when the guess loop is far away from the periodic orbit we are looking for. As the loop gradually approaches the periodic orbit, it becomes more and more accurate.

In order to decently depict the evolution along the orbit, we employ a second-order numerical scheme based on Eq. (13),

where $\bm{h}_{j}(s_{i})$ is a tiny vector in the $j$-$th$ direction $\bm{e}_{j}(s_{i})$, $\bm{h}^{\prime}_{j}(s_{i+1})$ is the tiny vector at $s_{i+1}$ after one step of evolution, while $\bm{h}_{j}(s_{i+1})$ is the vector at point $s_{i+1}$ before the evolution. In order to ensure that the new coordinate system is orthogonal, after each evolution it is necessary to carry out the Schmidt orthogonalization and normalization of $\{\bm{h}_{j}^{\prime}(s_{i})\}$ to get the new set $\{\bm{e}_{j}(s_{i})\}$ at each $s_{i}$,

where $c_{k}=\bm{e}_{k}(s_{i})\cdot\bm{h}^{\prime}_{j}(s_{i})$. As shown in Fig. 2, the evolution of the local coordinate axis in the $j$-$th$ direction is plotted. For each evolution, a very small vector $\bm{h}_{j}$ (black arrow) in some $j$-th direction at each lattice point $s_{i}$ evolves to $\bm{h}_{j}^{\prime}(s_{i+1})$ (red arrow) at $s_{i+1}$. After orthomalizing the vector $\bm{h}_{j}^{\prime}(s_{i+1})$, we get a new direction vector $\bm{e}_{j}(s_{i+1})$ (not displayed) at $s_{i+1}$. Repeating the procedure in all involved directions we will finally build up a new coordinate frame at $s_{i+1}$.

What is a proper dimension of the local coordinate system? Obviously, according to previous arguments it should not be less than the number of positive Lyapunov exponents of the system. The Lyapunov exponent equals $0$ in the direction of the velocity field, but the sampling points $x(s_{i})$ will approach the desired periodic orbit only along stable directions. Therefore, we include the velocity vector in the reduced coordinate subspace, as well as $2\sim 3$ other not-so-unstable or neutral directions to accelerate the convergence. To start with, at each lattice point we may select the same first few coordinate directions in the full phase space to build the initial local frame. Obviously, such a coordinate system is not in the stretching direction of the orbit, and we should have evolved it for some time to achieve that. Nevertheless, it does not hinder us from computing corrections in this frame and modifying it along the way. As the guess loop gradually approaches the periodic orbit by repeatedly solving Eq. (8)) , modifications of the coordinate system are made with the scheme illustrated in Fig. 2. In each modification, each local coordinate frame moves along the loop for $n_{1}$ steps (roughly $5\sim 15$ steps in the current work). If the total number of modifications during the whole search is $n_{2}$, in order to largely cover the expanding directions, $n_{1}n_{2}\Delta t\geq 1/\lambda_{L}$ should be satisfied, where $\lambda_{L}$ is some kind of average Lyapunov exponent. In the following, we find that $n_{1}\times n_{2}>N/2$ could do the job quite well, where $N$ is the number of lattice points depicting the periodic orbit.

We start the cycle search by numerical integration of the dynamical system under investigation and locate possible orbit segments that nearly repeat themselves - the recurring segments. Numerical experiments reveal regions where a trajectory spends most of its life, giving us the first hunch where to initialize a loop. After an extensive search for the recurring orbits, to start the variational scheme, we first take a spatial FFT of one such segment and keep only the lowest wavenumber components, of which an inverse Fourier transform back to the phase space yields a smooth loop , which will be used as our initial guess. In the following, we will apply the reduced scheme of the variational method discussed above to several examples.

The incompressible Navier-stokes equation on the torus $T^{2}=[0,2\pi]\times[0,2\pi]$ is written as[19, 20]

where $\bm{u}$ and $p$ is the velocity field and the pressure respectively, and $f$ is the external periodic driving. Under certain conditions, Eq. (16) could be expanded with $9$ Fourier modes[21], which reads

where $R_{e}=58$ is the Reynolds number.

For the system described by Eq. (IV.1), we calculated numerically its Lyapunov exponents which are displayed in Table 1. It can be seen that the first two values are greater than $0$. According to the previous discussion in Sec.III.1, the dimension of the local coordinate system is set to $n=4$ (the value of $\lambda_{3}$ and $\lambda_{4}$ are closest to $0$, which correspond to the orbit direction and the translation symmetry of Eq. (IV.1)). The results by the reduced method are displayed in Fig. 3 in which (a) and (b) depict the projections of the initial guess loop constructed by the FFT, while (c) and (d) plot the projection of the final cycle. The initial guess loop does not seem to have any symmetry, but the cycle is obviously symmetric. As can be seen from diagrams (a) and (c), the initial loop roughly possesses the topology of the periodic orbit it is about to approach. Therefore, even if the initial conditions are not very good, our method still drives the guess loop to a periodic orbit.

Here we further investigate the influence of the number of lattice points on the convergence. The same periodic orbit (as shown in Fig. 3) is used as the target with different numbers of lattice points $N=1000,~{}900,~{}800,~{}700,~{}600$ and the accuracy of the initial loop being kept the same. With the guess loop iterated for $100$ times, the convergence of our reduced method is shown in Fig. 4. It is clear that in the first $40$ iterations, the convergence are similar among different representations, but becomes different in the ensuing iterations. The larger the number of orbit points, the better the convergence. The fast convergence at the first stage stems from the initial exponential decay of the error in the projected subspace and thus is similar for different discretization. However, the relatively slow convergence at the second stage is enabled by the modifications of the coordinate frames which proceed distinctively for different $N$’s.

In the flowing we test our reduced method in case of the coupled Lorenz systems[22] $\sum_{i}:\{\dot{x}_{i}(t),\dot{y}_{i}(t),\dot{z}_{i}(t)\}$, described by the following equations

where the parameters $\sigma=10$, $\rho=28$, $\mu=8/3$ with a periodic boundary condition: $y_{\bar{n}+1}=y_{1}$. Here, $\bar{n}$ is the number of the coupled lorenz systems and $k_{i}$ is the coupling strength. In the uncoupled case $k_{i}=0$, all the subsystem perform chaotic motions on the well-known Lorenz attractor. Linear coupling is imposed via the $y_{i+1}$ components of the adjacent subsystem, controlled by $k_{i}$. Here, we take the case of $\bar{n}=33$ as an example, where the total dimension of the whole system is $d=99$. In order to make the coupling of the system non-uniform, we set the coupling strength as $k_{i}=[0.5+0.5\frac{(i-1)}{(\bar{n}-1)})]k_{0}$ (with $k_{0}=0\sim 1$).

For such a system, we may employ an alternative way of constructing an interesting initial loop. Since the system consists of $33$ coupled lorentz oscillators, we may use periodic orbits of the original lorentz system to construct the initial loop. For example, we take the same cycle for each lorenz system (as shown in Fig.5 (a)$\sim$(c), the periodic orbit is symmetrical with respect to $x=0$). When the coupling strength is $k_{i}=0$, the initial guess is already a true periodic orbit of the whole system. Recursively, we search a nearby new cycle upon gradually increasing the coupling constant, starting with the already determined cycle at the previous coupling. This provides an alternative way to start the search and a family of cycles are obtained when $k_{0}$ goes from $0$ to $1$. We have numerically calculated the Lyapunov exponents displayed in the Table 2, and found that the values of the first $33$ Lyapunov exponents are positive and away from zero corresponding to the main stretching directions. Such a result is also easy to understand, because there is an expanding direction in each Lorentz system, therefore the $33$ coupled Lorentz systems should have $33$ unstable directions as long as the coupling is small. According to the previous discussion in Sec.III.1, the dimension of the reduced space is set as $n=34$.

The final result is shown in Figs.5 (d) $\sim$ (f). By comparing the graphs in (d) and (e), we see that the periodic orbit of each individual Lorentz system does not seem to change much in the $x-y$ plane. However, the relation between corresponding variables in different subsystems exhibit nontrivial features (see Figs. 5 (f)). In this example, from the construction of the initial loop to the search of periodic orbit, the whole process is completed in local subspaces without evolution in the full phase space, which certainly cuts down the computation complexity and provides a good example to deal with high-dimensional systems with similar spatial structures. As shown in Fig.6, the algorithm converges faster initially and slows down when the guess loop is about to approach the true periodic orbit for the same reason given in the previous example. But on the whole, the convergence is good and remains exponential.

In previous example, the unstable directions are quite few so that the dimension of the reduced local coordinate subspace is also low, the convergence may be expected. Nevertheless, in the current example, there are many unstable directions but the convergence is still good with a proper construction of the local subspaces. As a matter of fact, the local dimension can be further reduced even down to $n=9$ and the method is still effective though with a slow convergence. How is it possible? During the evolution, all the local coordinate frames are essentially different. Each modification rotates these frames to new dimensions and the variational method will bring down errors in these directions. Finally, when all the directions are covered, a cycle may be identified. Nevertheless, the number of the included dimensions could not be too small to make the algorithm unstable. We still have no idea of what the necessary number of dimension is.

Here, we directly apply the current scheme to a nonlinear partial differential equation - the Kuramoto-Sivashinsky equation (KSe) to check the influcence of high dimensions and possible complications arising from possible stiffness due to the existence of inertial manifold. The KSe was independently proposed by Kuramoto and Tsuzuki in the study of phase turbulence in reaction-diffusion systems [23] and Sivashinsky in the study of flame combustion propagation models [24]. The KSe is an intrinsic phase equation describing the slow change of vibrations in a spatially extended system and also is an ideal model for studying high-dimensional systems. The KSe in one-dimensional space reads

where the subscript $x$ and $t$ represent the partial derivative of $u$ with respect to $x$ or $t$, with $t\geq 0$ and $x\in[0,2\pi]$. The $\nu$ is a viscous damping parameter which is the dissipation rate of the system, and the first term on the right hand side of Eq. (21) is a nonlinear convection term, which induces interaction between different spatial modes and energy transfer from low wave-number modes to high ones. As $\nu$ decrease, the system undergoes a series of bifurcations, leading to increasingly turbulent dynamics.

If we choose to study only the odd solutions $u(t,-x)=-u(t,x)$ with the periodic boundary condition $u(t,x+2\pi)=u(t,x)$, the state variable $u(x,t)$ can be expanded in spatial Fourier series [25],

where $a_{-k}=-a_{k}\in\mathbb{R}$. In terms of the Fourier components, Eq. (21) is rewritten as an infinite ladder of ODEs:

In numerical computation, a Galerkin truncation of the Fourier series has to be applied. As the amplitude of $a_{k}$ decreases with the increase of $k$ near the strange attractor, the high wavenumber modes in the asymptotic regime are negligible. Thus we can reduce the equation to a finite but large number of ODEs by the Galerkin truncation. In the current exploration, we work with two different truncations $d=16$ and $64$ in the following. However, in Eq. (23), the fourth power of $k$ arising from the dissipation could be quite large for large $k$, which makes the system stiff and brings trouble to numerical computation.

The Lyapunov exponents of the system are computed and filled in Table. 3. We show the values of the first eight Lyapunov exponents at $d=16$ ($\nu=0.015$) and $64$ ($\nu=0.0049$), where the numbers of positive values are $2$ and $4$. The exponents $\lambda_{3}=-0.0092$ and $\lambda_{5}=-0.0013$ correspond to the orbit direction of the system with $d=16$ and $d=64$, respectively. According to the previous discussion in Sec.III.1, we select $n=5$ and $n=8$ directions to form the local coordinate frames along the loop, including the velocity direction.

As shown in Figs. 7 and 8, the projections of the determined periodic orbits and the corresponding initial guess loops of with $d=16$ and $64$ are shown. It is found that the amplitudes of high-wavenumber modes are much smaller than the low ones ($\bm{a}_{11}/\bm{a}_{63}\sim 10^{4}$).

Nevertheless, the high-wavenumber modes have intricate features, which makes them sensitive to the initial guess, especially when it is rough. This is not surprising considering the stiffness problem mentioned after Eq. (23). Therefore, extra care should be exercised during the modification of high-wavenumber modes the corrections of which could be multiplied by someweight ranging from $1$ to $0.1$ to slow down the variation. In addition, it is not difficult to find that the topological structures of the high-wavenumber modes change a great deal during loop evolution when the guess loop is not good, which prevents a smooth convergence and possibly induces numerical instability. With further increase of the dimensions of the system, this problem becomes prominent (such as in the case of $d=64$).

In order to promote the robustness of the reduced method the average position of the neighboring points on the loop is used as the new value at $s_{i}$ (namely, $\bm{X}(s_{i})\longrightarrow[\bm{X}(s_{i})+\bm{X}(s_{i+1})]/2$) after each correction. Despite all the details dealing with a bad initial guess, in the Kuramoto-Sivashinsky system with $d=64$ Fourier mode truncation, for a subspace with dimension $n=8$, our reduced method has good convergence. As shown in Fig. 9, when the initial guess loop is close to a periodic orbit (entering the linearization neighborhood of the orbit), the error function Err. will decreases exponentially.

The space-time evolutions of $u(x,t)$, the unstable spatiotemporally solutions corresponding to the orbits in Fig. 7 and Fig. 8, are plotted in Fig. 10. As $u(x,t)$ is antisymmetric in $[-\pi,\pi]$, it suffices to display the solutions in the $[0,\pi]$ interval. There is little difference between the spatiotemporal profiles in Fig. 10 (a) and (b), which indicates that the guess loop is close to the periodic orbit. While big difference between (c) and (d) indicates that the initial guess is bad. A comparison of (b) and (d) reveals that the structure displayed in (d) has more features than that in (b) which indicates that the higher the dimension of the attractor, the finer the space-time evolution of the $u(x,t)$.

Often than not, the original variational method (VM) consumes a lot of storage resource (SR) and computer time which will become more and more striking with the increase of system dimension, but the reduced variational method (RM) introduced in the current manuscript partly solves the problem. As shown in Fig. 11, we compared the required storage and computing time for the VM and the RM. In Fig. 11 (a), the black-quare curve depicts the storage resource required by the VM, which increases exponentially with the increase of system dimension $d$, while the black-dot curve for the RM grows only linearly. The difference between them becomes more and more pronounced with the increase of system dimension $d$ (the ratio of storage resources required is $SR_{VM}/SR_{RM}\sim 6$ for $d=64$). In Fig. 11 (b), the black-square curve is the time required by the VM, and the black-dot curve by the RM. Because for different system dimensions, the number of lattice points is different, it is difficult to compare directly. A quantity $t_{VM}/t_{RM}$ (that is, the ratio of the time $t_{VM}$ required by VM to the time $t_{RM}$ required by RM)is defined and plotted with the red-circle curve in Fig. 11 (b). It is clear that with the increase of system dimension, the ratio $t_{VM}/t_{RM}$ becomes larger and larger, and reaches $t_{VM}/t_{RM}\sim 28$ at $d=64$. Certainly, when the dimension of the system is not high, the advantage of the RM is not so prominent. Henceforth, the caveat that lies in the variational method is essentially filled in the reduced scheme and we expect its fruitful application to high dimensional systems.