Computing Lyapunov Exponents using Weighted Birkhoff Averages

In this paper we will test several weight functions $g$ (8) to see how they influence the convergence of the time average (7) for the Lyapunov spectrum.

The standard choice for $g$ is the $C^{\infty}$ bump function of [20] ¹¹1Other possibilities with varying “widths” were explored in [33]. For time averages of functions on phase space, (9) was found to have the best convergence properties.

We also consider the skewed bump function

This is still $C^{\infty}$, but is no longer symmetric about $\tau=\tfrac{1}{2}$—it has a maximum at $\tau=\tfrac{1}{\sqrt{2}}$. As a third example, we will use the function

which is $C^{\infty}$ and monotone increasing on $[0,1)$ but has a discontinuity at $t=1$. This will test whether the computations of $\mu_{T}$ are more sensitive to initial transient behavior in the stretching of the vector $v_{t}$. These functions are shown in Fig. 1.

In [35], a modification of the average (3) was introduced to give a chaos indicator that they entitle the “mean exponential growth rate for nearby orbits” (MEGNO). As reviewed in [10], a generalized MEGNO can be formulated, labeled by a pair integers $(m,n)$. It is obtained by first computing a weighted average of the stretching factor (4) for the first $t$ iterates: ²²2To be consistent with the definition of the stretch (4), and the concept of a bump function, we shift the indices by one step from [10, Eqs. 4.38-39].

The $(m,n)$ MEGNO is obtained as an additional time average of this quantity:

Note that we can reorder this double sum to obtain an expression that closely resembles the weighted average (7):

where the MEGNO “weight” is effectively

Though $\bar{Y}_{(m,n)}$ now has a form similar to (7), the weight function is not normalized: MEGNO does not attempt to compute an accurate value for $\mu_{T}$. We propose that a normalized weight function would be more appropriate, so we rescale (12) to define

where $N_{T}$ is the normalization constant as in (8). Then the MEGNO-weighted average for the Lyapunov exponent is defined using this weight in (7).

According to [10], the most useful cases of MEGNO correspond to $(m,n)=(1,-1)$ and $(2,0)$. We will compare our results only with the second case since it was shown to converge more rapidly. When $n=0$ the sum in (13) is trivial. Normalizing then gives

Note that this function vanishes at $t=0$ and $T$: it is a bump function like those in §II.1. To compare more directly with these, we rescale time and set the interval to $[0,1]$ to obtain

and then $w_{T}^{meg}$ is given by the normalization (8). This weight is $C^{1}$ at the endpoint $\tau=0$, but only $C^{0}$ at $\tau=1$. This function is the orange curve in Fig. 1.

Recall that the weighted Birkhoff average $\mathit{WB}_{T}(h)$ for a function $h\in C^{\infty}(M,{\mathbb{R}})$ on phase space converges slowly when an orbit is chaotic but for “regular” orbits (those smoothly conjugate to incommensurate rotations) it converges at a rate rate determined by the smoothness of the weight function—for a $C^{\infty}$ weight, such as (9), this can be super-polynomial [20, 21, 22, 34]. By contrast, the standard unweighted average nominally convergences at best as $T^{-1}$ [32].

Here we similarly observe that the Lyapunov spectrum also converges slowly whenever the orbit is chaotic, regardless of the method used. But for a regular orbit, the weighted averages do typically enhance the convergence of the Lyapunov spectrum.

As a first example, consider the three-dimensional “discrete Lorenz map” [36]

Here we fix two of the parameters ($\nu_{1}=-0.85$ and $\nu_{3}=0.95$, in the notation of [36]) and allow only the parameter $\nu_{2}$ to vary. For (16), the determinant of the Jacobian is independent of the point: $\det(Df)=\nu_{1}\nu_{3}=-.8075$. This implies that the sum of the exponents should be

However, in the computations below, we do not use this since we want to test convergence of the individual exponents (5).

Figure 2 shows the three Lyapunov exponents for this map as a function of $\nu_{2}$, computed using the standard WBA weight (9) with the iterative Gram-Schmidt method (5). The orbit is arbitrarily chosen to start at $(0,-0.01,0.0001)$, discarding the first $4000$ iterates to remove transients. The exponents are computed using (7) for the next $T=2(10)^{4}$ iterates. In all cases, $|\mu^{(1)}+\mu^{(2)}+\mu^{(3)}-d^{*}|<10^{-14}$, consistent with the constant Jacobian determinant. This excellent convergence follows from the fact that, up to floating point error, $\ln(r_{t}^{(1)})+\ln(r_{t}^{(2)})+\ln(r_{t}^{(3)})=d^{*}$ for each $t$, so the sum of the averages is the average of the sum, with or without a weight function.

As shown in [36], the fixed point $(x,y,z)=(0,0,0)$ of the map (16) is stable up to $\nu_{2}=1.85$ where it undergoes a pitchfork bifurcation—at this point the largest exponent (blue in the figure) hits zero. The newly created pair of fixed points lose stability at $\nu_{2}\approx 1.8645$ in Neimark-Sacker bifurcations. The resulting pair of attracting circles or periodic orbits have basins of attraction limited by the stable and unstable manifolds of the origin, and there can be additional attractors. At $\nu_{2}=1.87$ there is a chaotic, Lorenz-like attractor, as shown in Fig. 3(a). This attractor has a single positive Lyapunov exponent and a tangential exponent of zero. Using the exponential weight (9) with $T^{*}=2(10)^{6}$ gives the exponents

As $\nu_{2}$ varies, the Lorenz-like attractor can collapse onto an attracting circle; for example this occurs at $\nu_{2}=1.8785$, see Fig. 3(c). For this attracting circle the maximal Lyapunov exponent is zero, and again using the exponential weight (9), we find (to higher accuracy)

for the same $T^{*}$. These cases are marked with vertical dashed lines in Fig. 2.

The convergence of the largest exponent using the five weight functions of §II is shown in Fig. 3, panels (b) and (d), for $\nu_{2}=1.87$ and $1.8785$, respectively. These plots show the errors for $100$ logarithmically spaced values of $T\in[900,1.5(10)^{6}]$. In each case the first $4000$ iterates of the initial condition are discarded to remove transients. The error is estimated by comparing to $\mu^{(1)}_{T^{*}}$ with $T^{*}=2(10)^{6}$.

As seen in Fig. 3(b), all of the methods perform poorly for the chaotic attractor at $\nu_{2}=1.87$: the convergence is at best like $T^{-1}$. Consistent with this, one might believe the results to 5 or 6 digits at $T=10^{6}$. The convergence of the second and third Lyapunov exponents (not shown) is essentially indistinguishable from that shown in panel (b).

For the invariant circle at $\nu_{2}=1.8785$, the WBA and skew weights, which are $C^{\infty}$, far outperform the other methods. The best convergence is for the skew weight; it essentially reaches machine precision by $T=3(10)^{5}$. The MEGNO weight (15) also gives increased convergence at a rate nearing $T^{-2}$. The left and constant (regular) weights still converge as $T^{-1}$: there is no improvement since these weight functions are not continuous.

The behavior seen in Fig. 3 appears to be typical: we have seen similar performance for convergence of the Lyapunov spectrum in many other simulations for chaotic and nonchaotic orbits for various maps including the discrete Lorenz map for other parameter values, the classic Hénon map, the Derived from Anosov (DA) map [37], and the 2D torus map studied in [20, Section 3.7.3]. Similar behavior is also seen for Poincaré maps of flows, such for a periodically forced, double-well Duffing oscillator [38].

We now describe several cases where the convergence does not follow the pattern seen in §III.1.

Dynamics with Shear: Integrable symplectic maps have families of invariant tori whose rotation vectors generically vary across tori: they have shear. This results in the linear growth of the length of vectors transverse to the tori, and this gives the slow convergence

to zero for (3). It has long been known [39] that shear causes a similar slow convergence even when the map is not integrable, whenever the orbit lies on an invariant torus [28]. It appears to be better for the weighted cases, as we describe below and see in Fig. 4.

For example, consider the 3D volume-preserving map for $(x,y,z)\in{\mathbb{T}}^{2}\times{\mathbb{R}}$ [22]:

The initial point $(x,y,z)=(0,0,-0.05)$ appears to lie on an invariant two-torus that is a graph over $(x,y)$ on which the dynamics has the incommensurate rotation vector $\omega\simeq(0.544519,0.411571)$ [22]. All three of the Lyapunov exponents for this orbit are zero: the two tangential exponents vanish because the invariant set is a two-torus, and the transverse exponent is then zero because of volume-preservation. Nevertheless, since the rotation vector varies with $z$, the map has shear and the length of a vector transverse to the torus will grow linearly. This should result in slow convergence of the exponents.

The convergence of a Birkhoff average of the function $h=\cos(2\pi x)$ and of the largest Lyapunov exponent are shown in Fig. 4 for the five weight functions of §II. The figure shows the averages for $100$ values of $T\in[800,1.5(10)^{7}]$. Since the map is volume preserving, no transient removal is needed. For the $C^{\infty}$ weight functions the convergence rate of the Birkhoff average is excellent, as we would expect for a quasiperiodic orbit [22]. The errors for the Lyapunov exponent in panel (b) were computed by comparing to $\mu^{(1)}=0$. Note that this convergence is very slow—even though the orbit is nonchaotic. The other two Lyapunov exponents (not shown) also have the same convergence rates.

We observe similar behavior for other parameters and orbits of the map (17) as well as for nonchaotic orbits of the 2D Chirikov standard map.

Weak Chaos: A dynamical system has weak chaos when the Lyapunov exponents are not positive, but it still has sensitive dependence on initial conditions. Such dynamics can lead to strange nonchaotic attractors (SNA), where the orbit lies on a geometrically strange (fractal) set, [40, 41]. In this case we have previously observed that averages converge slowly for both Lyapunov exponents and functions on phase space, and adding a weight function does not improve convergence [33, 42].

Noninvertible maps: Noninvertibility makes the computation of Lyapunov exponents more delicate. Denote the set of points where the Jacobian of the map is singular by

If an invariant set intersects $J_{0}$, it may be nonsmooth and even self-intersecting, see [43, 44, 45] and references therein. If an ergodic component intersects $J_{0}$, even if it is smooth, some of the exponents will be undefined, since a zero determinant implies that $r_{t}^{(j)}=0$ for some $j$, so that $\ln r_{t}^{(j)}$ is infinite. Moreover, the average (5) will also be undefined for initial conditions on the dense, countably infinite set of preimages of $J_{0}$. Of course it is still possible for the Lyapunov exponents to exist almost everywhere. However, from the numerical standpoint, if an orbit nears the dense set of singular points, then at least one $r_{t}^{(j)}\approx 0$, and this will lead to significant floating point errors.

This phenomenon is shown in Fig. 5 for a so-called Tinkerbell map [7]

In this case $J_{0}$ is a circle centered at $(-0.2,-0.65)$ with radius $\sqrt{0.125}$ that intersects the attractor, see Fig. 5(a). The Lyapunov exponents for this attractor are nonpositive, $\mu\simeq(0,-0.14292)$, and as seen in Fig. 5(b), the convergence of $\mu^{(1)}$ for the $C^{\infty}$ bump functions is excellent: since $J_{0}$ intersects the orbit transversally, $r_{t}^{(1)}$ is never near zero. However, $r_{t}^{(2)}$ does get arbitrarily close to zero infinitely often on the orbit. As seen in Fig. 5(c), this results in slow convergence of the second Lyapunov exponent for all methods. For this case, the attractor is smooth, and we observe that the convergence of $\mathit{WB}_{T}(h)$ using (9) for functions on phase space (such as $h=\cos(2\pi x)$ and the rotation number) is excellent (not shown). The implication is that for this case, the numerical difficulities are restricted to the smaller Lyapunov exponent.

Constant Jacobian: When the map $f$ has a constant, hyperbolic Jacobian, the convergence rates of $\mathit{WB}_{T}(S)$, (7) are enhanced for smooth weight functions. Indeed, when $Df=A$ is constant, the method given in (5) is equivalent to a computing the singular values $\sigma_{j}$ of $A$ via normalized simultaneous power iteration. As long as $0\leq\sigma_{j+1}<\sigma_{j}<\dots<\sigma_{1}$ (the strict inequality is the generic case), we have $r_{t}^{(j)}\to\sigma_{j}$ with error ${\cal O}(|\sigma_{j+1}/\sigma_{j}|^{t})$ [46]. As a result, the errors in the weighted average (5), occur only for the initial transients. The implication is that the left (11) and skew (10) weights perform better than the others, since they suppress the initial portion of the average. In contrast, when the orbit is chaotic, $\mathit{WB}_{T}(h)$ for functions $h$ converges slowly for any weight function.

We have verified this improved convergence of $\mu^{(j)}_{T}$ and slow convergence of other averages for examples including Arnold’s cat map and the skinny baker map [7], both of which are uniformly hyperbolic.