Predictability of the onset of spiking and bursting in complex chemical reactions

Abstract

For three complex chemical reactions displaying intricate dynamics, we assess the effectiveness of a recently proposed quantitative method to forecast bursting and large spikes, i.e. extreme events. Specifically, we consider predicting extreme events in (i) a copper dissolution model where Bassett and Hudson experimentally observed homoclinic (Shilnikov) chaos, (ii) a model derived from the mass action law of chemical kinetics, and (iii) an autocatalator model. For these systems, we describe how the alignment of Lyapunov vectors can be used to predict the imminence of large-amplitude events and the onset of complex dynamics in chaotic time-series of observables.

---

1. Introduction

The prediction of catastrophic events, outliers, is an enduring problem that continues to attract attention in virtually all scientific disciplines.^1–3 Such extreme events are usually quite disruptive and manifest themselves as large-amplitude fluctuations in time-series of observables, either measured or computed. A plethora of methods to anticipate the onset of catastrophic events exist. Recently, we proposed a quantitative criterion to predict larger-amplitude events in chaotic systems.^4,5 The criterion is obtained by studying the alignment of Lyapunov vectors (LVs), sometimes referred to as “covariant” Lyapunov vectors. The norms of the LVs, called Lyapunov exponents (LEs), measure the temporal expansion (or contraction) of the distance between two initially very close trajectories. Positive LEs are the signature of chaotic trajectories. We observed empirically^4,5 that a tangency between LVs occurs every time that an observable under consideration tends towards its maximum. This tangency is an alignment of Lyapunov vectors. In some cases the value of the maximum is inversely proportional to the degree of the vector alignment.

In the present paper, as a further application of the proposed criterion,^4,5 we investigate the reliability of the alignment between Lyapunov vectors as a predictive tool in the context of complex chemical phenomena. Specifically, we use Lyapunov vectors to predict the onset of bursting and spiking in three familiar reactions. The first reaction was considered in pioneering work in 1988 by Bassett and Hudson⁶ to reproduce aspects of the experimental electrodissolution of a copper rotating disk in a H₂SO₄/NaCl solution. The second reaction is a mass-action model used to study the interplay of three interacting chemical species.^7,8 The third reaction is a popular autocatalator model used to study complex oscillations in isothermal chemical systems.⁹ In all three examples, we find the alignment between Lyapunov vectors to be an effective predictor of extreme events. In addition, we point out some pitfalls of the method and discuss how to apply it effectively.

Lyapunov vectors provide a detailed step-by-step record of what exactly happens with the angles between stable and unstable manifolds during the whole evolution of a given dynamical process. The computation of Lyapunov vectors is no more complicated than the computation of standard Lyapunov exponents and is well-described in the literature.^4,5,10–15 In an Appendix, we summarize the steps involved in their computation. As discussed below, it was observed heuristically^4,5 that the onset of large peaks in a given physical variable corresponds to the alignment of LVs along the flow direction. Therefore, since Lyapunov vector alignment precedes large peaks, it can be used to predict the latter. It has been shown⁴ that predictability using LVs is closely related to alignment conditions, obtained by defining thresholds of how close LVs become tangent to the flow direction. When the alignment conditions are satisfied at a given time, called the prediction time, an extreme spike is expected to occur. The appropriate definition of the alignment conditions depends sensitively on the specific dynamics of the system.

For systems with long quiescent intervals between large spikes, the alignment conditions are very useful to anticipate unexpected large spikes. However, when numerous smaller spikes occur between large spikes, the choice of the appropriate alignment conditions can be tricky and some preparatory trial-error checks may be needed. In general, while small and large spikes can be predicted by the alignment of LVs, large spikes correspond to the prediction times where LVs are strongly parallel to the flow direction. Thus, the amplitude of the predicted peaks is related to how close the LV direction aligns along the flow direction.⁴

Apart from the chemical phenomena described by the models, complex chemical reactions are governed by sets of differential equations which are significantly different from the equations so far considered and no knowledge concerning the interplay of their Lyapunov vectors exists. Thus, it is of interest to verify the predictive power of Lyapunov vectors in foretelling the imminence of large-amplitude events in chemical time-series. Here, we check predictions for three typical systems governed by nonlinear differential equations involving different combinations of polynomial terms. We computed LVs using the methodology explained in detail in ref. 4 and 5. See also ref. 10–14 and the survey in ref. 15.

2. Electrodissolution of copper

In 1988, Bassett and Hudson⁶ performed a series of experiments on the electrodissolution of a copper rotating disk in a H₂SO₄/NaCl solution. They reported evidence of homoclinic (Shilnikov) chaos in their experimental time-series. Among this evidence, they presented a comparison of their time-series with those obtained from a mathematical model considered previously by Argoul et al.,¹⁶ in a discussion concerning homoclinic chaos in the Belousov–Zhabotinski reaction. The model is defined by the following three-dimensional flow

ẋ = y, (1)

ẏ = z, (2)

ż = −z − 1.3y − μx + x² − 1.425y² + 0.2xz − 0.01x²z. (3)

Here, (x, y, z) represent chemical concentrations and (ẋ, ẏ, ż) are their corresponding time derivatives. For μ < 1.3 the origin is a stable steady state and for μ = 1.3 a subcritical Hopf bifurcation occurs where the stable state becomes unstable and Shilnikov chaos appears.

Fig. 1 shows the chaotic attractor for μ = 1.4 and initial conditions (−0.1, −0.1, 1.0). The attractor is obtained by solving eqn (1)–(3) numerically with the standard fourth-order Runge–Kutta algorithm and a fixed time-step h = 0.005. The first 10⁷ steps were discarded as transient, with LVs computed subsequently for approximately 9.6 × 10⁶ steps. The Lyapunov exponents λ_i (i = 1, 2, 3) are λ₁ = 0.234, λ₂ = −0.0013, and λ₃ = −2.739, where the subindex i = 1 refers to the unstable direction, i = 2 to the flow direction and i = 3 to the stable direction. Accordingly, θ₁₂, represented in colors in Fig. 1, is the angle between the unstable manifold and the flow direction. We observe that along the attractor the angle θ₁₂ can change substantially.

Fig. 1 Representative attractor for the model of dissolution of copper, eqn (1)–(3), with colors representing the values of the angle θ₁₂ between the unstable invariant manifold and the flow direction.

In order to get a more quantitative insight concerning the ability of the alignment of LVs to anticipate spikes, Table 1 summarizes some statistical results. For this end, we consider the chaotic trajectory from Fig. 1 which has 2192 spikes over the integrated times, and say that there is a spike whenever x < −9.0 and an extreme spike when x < −15.0. To explain the quantities in Table 1, let us start by discussing the third line, which is the alignment condition displayed in Fig. 2. From the sample of 2192 spikes, N_pp = 1693 of them were predicted, which is around 76.3% of the total number of spikes. However, N_fa = 686 (41.0%) from the 1693 predicted spikes are false alarms, meaning that 1007 spikes are not extreme. False alarms are those alarms which predicted spikes with small amplitudes, i.e. predicted non-extreme spikes. From the total number 2192 − 1693 = 499 of non-predicted spikes, N_ext = 181 of them were extreme spikes. Thus, for this alignment condition the method did not predict 181 extreme spikes. Such a large value of non-predicted extreme spikes has to be avoided. All predicted extreme spikes have an associated prediction time interval Δt_p, which can be large or short. The prediction time interval is defined as the absolute value between the prediction time and the time for which the extreme spike occurs. While large prediction time intervals are most welcome to allow preparation to meet the extreme events, short prediction time intervals can bring difficulties for such preparations. In other words, for very short prediction time intervals there may not be enough time to prepare for the impingent extreme events. To quantify this, the last row displays the mean prediction time interval .

Table 1 Statistics for the chaotic trajectory with 2192 peaks (x < −9.0) for distinct values of thresholds θ^(min)₁₂ and θ^(max)₁₂. The quantities are: N_pp, the number of predicted spikes; N_fa, the number of false alarms, i.e. predicted spikes with smaller amplitudes (x > −15.0); N_ext, the number of non-predicted extreme spikes (x < −15.0) and , the average prediction time interval

θ ^(min)12	θ ^(max)12	N pp	N fa	N ext
0.05	3.09	759 (35%)	315 (41%)	725	19.1
0.10	3.04	1205 (55%)	484 (40%)	447	17.4
0.14	3.0	1693 (76%)	686 (41%)	181	15.2
0.20	2.94	1933 (88%)	795 (41%)	30	14.1
0.30	2.84	2146 (98%)	982 (45%)	4	15.1
0.50	2.64	2186 (99%)	1022 (47%)	4	18.2
0.70	2.44	2187 (99%)	1022 (47%)	3	20.4

---

Fig. 2 A time window displaying the evolution of x (black curve), of θ₁₂ (red), and two straight lines (blue) indicating the values of θ^(min)₁₂ = 0.14 and θ^(max)₁₂ = 3.0 for the alignment condition.

To better understand the relevance of θ₁₂, Fig. 2 shows a time-window displaying the variable x (black curve) from the chaotic attractor from Fig. 1. The angle θ₁₂ is plotted in the red color. Apart from some small amplitude oscillations, θ₁₂ tends to 0 or π before spikes appear in variable x. This means that the peaks in x can be predicted in time if we define thresholds for θ₁₂, namely θ^(min)₁₂ and θ^(max)₁₂. The two horizontal blue lines in Fig. 2 are illustrative thresholds for this trajectory, namely θ^(min)₁₂ = 0.14 and θ^(max)₁₂ = 3.0. Instants of time where θ₁₂ crosses these thresholds (i.e., θ₁₂ < 0.14 or θ₁₂ > 3.0) are called prediction times and a large spike in x is expected to occur. Three such examples are indicated with black arrows in Fig. 2, the arrows chosen subjectively, for illustration purposes. For the first two arrows on the left we have θ₁₂ > 3.0 and for the arrow on the right θ₁₂ < 0.14, and they precede the large spikes with x < −15. The thresholds define the alignment condition for the LV along the flow direction and allows us to predict the spikes. The effectiveness to predict spikes may change for different alignment conditions, as discussed below in more detail.

The overall message from Table 1 is that for stronger alignment conditions the number of non-predicted spikes increases, specially the extreme spikes, N_ext. If the alignment conditions are relaxed, the number of predicted peaks increases but the number of false alarms also increases. Comparing the stronger alignment conditions with the more relaxed condition, N_pp increases by around 65.2% but N_fa increases by only 5.2%. In fact, more relaxed alignment conditions lead to longer prediction time intervals. The selection of a suitable threshold depends very much on what is needed for practical purposes. For example, if we want to avoid all extreme spikes, we should chose an alignment condition consistent with N_ext → 0. In this case the appropriate alignment condition would be θ^(min)₁₂ = 0.7 and θ^(max)₁₂ = 2.44. On the other hand, if the large peaks provoke no major damage, one may accept some such peaks together with the condition to prevent a repetition of false alarms. This means, a compromise must be found between allowing some large peaks while at the same time not allowing the number of false alarms to dominate. To find such a compromise, one may consult Table 1 and select the suitable alignment condition providing an adequate balance. In such case, an adequate alignment condition would be θ^(min)₁₂ = 0.20 and θ^(max)₁₂ = 2.94.

3. Mass-action kinetics

In the next example, we revisit a chemical model introduced by Gaspard and Nicolis^7,8 to study homoclinic phenomena. The dimensionless governing equations are

ẋ = (βx − fy − z + g)x, (4)

ẏ = (x + sz − α)y, (5)

εż = x − az³ + bz² − cz. (6)

As before, of interest is the temporal evolution of the concentrations (x, y, z), which depend on the control parameters a, b, ε, f, g, s, c, α, β. As usual, we kept fixed the constants a = 0.5, b = 3, ε = 0.01, f = 0.5, g = 0.6, s = 0.3 and c = 4.8 and change the pair (α, β). Fig. 3 displays the chaotic attractor for α = 0.7825 and β = 0.55, starting from the initial condition (x, y, z) = (1, 1, 1). Colors represent the values of the angle θ₁₂ for each point on the attractor. The attractor consists of a chaotic trajectory with Lyapunov exponents λ₁ ∼ 0.13, λ₂ ∼ 0.0005 and λ₃ ∼ −544. Accordingly, the angle θ₁₂ is the angle between the unstable manifold and the flow direction. The trajectory from Fig. 3 and all subsequent simulations in this section were obtained by solving eqn (4)–(6) numerically with the standard fourth-order Runge–Kutta algorithm with a fixed time-step h = 0.00005. The first 10⁷ steps were discarded as transient, with LVs computed subsequently for approximately 10⁷ steps.

Fig. 3 Chaotic attractor governed by the mass-action kinetics, eqn (4)–(6), with colors representing the angle θ₁₂ as indicated by the colorbar.

The trajectory forming the attractor in Fig. 3 starts parallel to the xy plane with a chaotic spiraling-out motion, with the angles θ₁₂ oscillating between π/4 (cyan) and π/2 (green), while the reinjection loops are shown in black and red, associated with θ₁₂ approaching either 0 or π, respectively.

To better understand the time evolution of θ₁₂, Fig. 4 presents the evolution of z (black lines) and the associated LV angle θ₁₂ (red lines) for a regular trajectory in Fig. 4(a) and for a chaotic trajectory in Fig. 4(b). By ‘regular’ trajectory we mean a trajectory which has no positive Lyapunov exponent. In both cases, regular and chaotic, the LV angle θ₁₂ oscillates fast in time but tends to approach 0 or π before the spikes in z are observed. This means that the peaks in z can be predicted in time when introducing suitable thresholds for θ₁₂, namely θ^(min)₁₂ and θ^(max)₁₂. As exemplified in Fig. 4, the two horizontal blue lines indicate thresholds for these trajectories, namely θ^(min)₁₂ = 0.1 and θ^(max)₁₂ = 3.04 for the regular case and θ^(min)₁₂ = 0.5 and θ^(max)₁₂ = 2.64 for the chaotic case.

Fig. 4 Time window showing the variable z (black lines) and θ₁₂ (red lines) for part of (a) a regular trajectory with parameters α = 0.7825 and β = 0.39213 and (b) the chaotic trajectory from Fig. 1. Two horizontal blue lines indicate thresholds to predict spikes and the two black arrows show two prediction times where θ₁₂ crosses the thresholds.

In order to choose adequate thresholds for these trajectories, a statistical analysis of the spikes must be performed, as done in the next paragraph. When θ₁₂ crosses these thresholds, spikes in z are to be expected. We note that θ₁₂ tends to approach 0 or π not only before the large spikes, but also before the small spikes. The crucial distinction is that, preceding large spikes, the angle comes closer to 0 or π than before the small spikes. From a physical point of view, our results show that each time the unstable (chaotic case) and the stable (regular case) manifolds align along the direction of motion, a spike in the z concentration is expected.

Fig. 4 also shows that, far from spikes, the angle θ₁₂ oscillates in time. Such oscillations occur when the eigenvalues of the Jacobian matrix are complex and no invariant manifolds exist. In fact, manifolds rotate and the angles oscillate in time. This oscillatory behavior was observed and explained for specific parameter values of the Hénon map.⁵ Analogous angle oscillations were reported recently for interesting critical transition problems,¹⁷ although no explanation concerning the origin of the oscillations was given there.

The results above suggest that for times called prediction times, as exemplified in Fig. 4 by black arrows, it is possible to anticipate the imminent occurrence of a spike. It also becomes clear that prediction time intervals, defined as the absolute values between the prediction times and the times for which the extreme spike occurs, depend on the thresholds chosen.

Similar to before, Table 2 summarizes some statistical results. For this we consider the chaotic trajectory which has 549 spikes during the integration interval. We say that there is a spike when z > 0.9 and an extreme spike when z > 3.0. Let us start to discuss the first line, which has a very strong alignment condition. From the total of 549 spikes, N_pp = 526 of them were predicted. However, N_fa = 33 from the predicted spikes are false alarms, meaning that the 33 spikes are not extreme (z < 3.0). From the 549 − 526 = 23 non-predicted spikes, N_ext = 16 are extreme spikes. The last row displays the mean prediction time intervals .

Table 2 Statistics for the chaotic trajectory with 549 peaks (z > 0.9) for distinct values of thresholds θ^(min)₁₂ and θ^(max)₁₂. The number of extreme spikes (z > 3.0) for this trajectory is 509. The quantities are: N_pp, the number of predicted spikes; N_fa, the number of false alarms, i.e. predicted spikes with smaller amplitudes (z < 3.0); N_ext, the number of non-predicted extreme spikes (z > 3.0) and , the average prediction time interval

θ ^(min)12	θ ^(max)12	N pp	N fa	N ext
0.01	3.13	526 (95.8%)	33 (6.3%)	16	14.7
0.02	3.12	543 (98.9%)	37 (6.8%)	3	16.9
0.03	3.11	549 (100%)	40 (7.3%)	0	17.4
0.04	3.10	549 (100%)	40 (7.3%)	0	17.6
0.50	2.64	549 (100%)	40 (7.3%)	0	18.6

---

The general behavior presented in Table 2 is similar to the one in Table 1. For strong alignment conditions the number of non-predicted spikes increases, especially the number N_ext of extreme spikes. If the alignment conditions are relaxed, there is an increase in the number of predicted peaks as well as in the number of false alarms. However, the percentage of false alarms is much smaller here. Comparing the stronger alignment conditions with the more relaxed condition, one sees that N_pp increases by 4.2% while N_fa increases by only 1.0%. More relaxed alignment conditions require longer prediction time intervals, excluding almost all shorter (<1.0) prediction time intervals.

As discussed in the previous section, the appropriate alignment condition depends on what is sought in practical situations. If we again want to find a compromise between allowing some large peaks while at the same time not allowing the number of false alarms to dominate, consulting Table 2 the suitable alignment condition providing an adequate balance is θ^(min)₁₂ = 0.03 and θ^(max)₁₂ = 3.11.

4. The autocatalator model

Next, we consider the use of LVs to predict outliers in the complex oscillations generated by the isothermal chemical system governed by the following four-parameter autonomous polynomial model⁹where α, β, γ are dimensionless concentrations and μ, κ, σ, δ are control parameters of the reaction.

Fig. 5 displays the attractor for the typical parameters κ = 2.5, δ = 1, μ = 0.2983 and σ = 0.013 using the initial condition (α, β, γ) = (0.5, 2.5599, 0.5). As before, colors represent the values of the angle θ₁₂ for each point on the attractor. The attractor is formed by a chaotic trajectory with Lyapunov exponents λ₁ ∼ 0.755, λ₂ ∼ 0.002 and λ₃ ∼ −11.117. The attractor in Fig. 5 and in subsequent simulations was obtained by solving eqn (7)–(9) numerically with the standard fourth-order Runge–Kutta algorithm and a fixed time-step h = 0.05. The first 10⁷ steps were discarded as transient, with LVs computed during the next 6 × 10⁶ steps. For better visualization, only 12% of these points were plotted in Fig. 5. The angle θ₁₂ changes along the attractor. The rightmost green ridge clearly visible in Fig. 5 is related to angles close to π/2 and corresponds to strong hyperbolic motion.

Fig. 5 Attractor for the autocatalator, eqn (7)–(9), with colors representing the angle θ₂₃.

Fig. 6 presents a time-window of the evolution of β (black lines) and associated LVs for the chaotic trajectory from Fig. 5. Fig. 6(a) presents the angle θ₁₂ between the unstable manifold and flow direction, while Fig. 6(b) displays the angle θ₂₃ between the stable manifold and flow direction. Besides the very fast variations of these angles, they assume values close to 0 and π for almost all integrated times.

Fig. 6 Time evolution of the concentration β (×3) (black lines) for a chaotic trajectory and (a) the angle θ₁₂ (red line) and (b) θ₂₃ (green line).

Next, we consider a regular trajectory obtained for κ = 2.5, δ = 1, μ = 0.29776 and σ = 0.013 and initial conditions (α, β, γ) = (0.5, 0.25599, 0.5). The Lyapunov exponents for this case are λ₁ ∼ 0.003, λ₂ ∼ −0.914 and λ₃ ∼ −16.16. Therefore the angles θ₁₂ and θ₁₃ are the angles between stable manifolds along the flow direction. The time evolution of the concentration β is shown in Fig. 7 (black lines) with colors representing the angle θ₁₂ obtained along the trajectory. Results are analogous to the chaotic case discussed above, also for θ₂₃. Since for the autocatalator the variables θ₁₂ and θ₂₃ show essentially the same behavior, we use θ₂₃ to show that the prediction method also works when using angles between other manifolds.

Fig. 7 Time evolution of the concentration β (black lines) for a regular trajectory and (a) the angle θ₁₂ (red line) and (b) θ₁₃ (green line).

In sharp contrast to the reaction considered in Sections 2 and 3, here the variable of interest, namely the concentration β, never relaxes to random quiescent intervals. It always increases towards a peak (large or not) and decreases very fast to a minimum. After each minimum, β promptly starts to increase towards the next maximum, while θ₁₂ and θ₁₃ tend to approach 0 or π, usually right after the minimum. This is the alignment of both stable and unstable manifolds along the direction of motion and precedes the appearance of a peak. This is in agreement with our recent proposal⁴ that the alignment of LVs along the flow direction precedes peaks in a chaotic dynamical system. Therefore all peaks in the concentration β can be easily predicted by the angles θ₁₂ and θ₁₃. However, since the concentration β always goes, after the minima, towards the peaks, the alignment of LVs to predict peaks is not necessary in this case. In fact, the dynamics is apparently strongly non-hyperbolic for all the duration of the integration.

For all three models considered here, the alignment of LVs is able to predict all peaks during the time evolution, large or not. Therefore, while the comparatively large number of false alarms found in Section 2 (∼46.7%) is compatible with this detection power, false alarms found in Section 3 (∼7.3%) are acceptable. However, this seemingly negative aspect of the method can be essential to predict spikes in systems where the quiescent motion does not give any indication that a spike will appear. For example, in the quiescent intervals shown in Fig. 2 and 4, the plotted quantities x and z, respectively, present some smaller spikes with increasing amplitudes which indicate that an extreme spike might appear. All these smaller spikes are detected by the alignment of LVs and lead to false alarms. Sometimes, the time evolution of the observed quantities x, z and β in the examples above provides clues that a large spike will appear. But the difficult situation is when such clues are not found in the time evolution. In this circumstance, LV alignment is quite helpful to predict large spikes. Comparing results for the two systems in Sections 2 and 3, one sees that when the alignment conditions are relaxed, the number of predicted peaks increases much more than the number of false alarms. In addition, the number of unpredicted extreme spikes tends to zero very fast. All in all, for the reactions considered, the LV alignment is able to foretell the imminence of spikes.

5. Conclusions and outlook

This work assessed the alignment of LVs as a quantitative tool to predict large-amplitude spikes in three chemical models governed by distinct polynomial dynamics. For the models in Section 2 and 3, the time evolution of the concentrations presents random intervals of quiescent motion. In these cases, large events are unpredictable by analyzing just the time evolution of the concentrations. Here, LVs are very helpful to predict spikes. In the model presented in Section 4 the concentrations oscillate continuously between minima and maxima. There is no relaxation and the arrival of spikes can be recognized by the time evolution of the concentrations. The three examples considered show that the alignment of Lyapunov vectors along the flow direction provides a straightforward and reliable quantity to predict small and large events for typical chemical reactions.

The LV alignment as a spike prediction method is expected to be very general since it is able to detect all spikes. It is particularly useful when random intervals of quiescent motion are present in the time evolution of the concentrations. Due to the lack of analytical results, the relationship between the LV alignment and extreme events in dynamical systems must be tested for each particular case. An open challenge is to devise a way to apply LVs (i) to observational data, namely, for data which lack an underlying mathematical framework, and (ii) to realistic systems where the presence of noise is inevitable.

Conflicts of interest

There are no conflicts to declare.

Appendix – numerical procedure

The most unstable direction of a chaotic trajectory can be determined from the direction of the eigenvector related to the largest LE, obtained as an eigenvalue of the Jacobian matrix of the system. The procedure for obtaining it is well known and widely used in the nonlinear and complex systems community and an explicit numerical code for the determination of Lyapunov exponents is given by Wolf et al.¹⁸ It involves a Gramm–Schmidt (GS) orthonormalization procedure which, unfortunately, destroys the information about the direction of all eigenvectors not related to the largest Lyapunov exponent. To compensate for this loss, a method proposed by Ginelli et al.,¹⁴ based on forward and backward time integrations, allows one to recover the direction of all unstable and stable invariant manifolds. This procedure is somewhat time consuming but very precise. In order to obtain stable and unstable directions in phase space we need to perform (i) forward and (ii) backward time evolutions:

(i) The procedure summarized here is well known and uses ordinary differential equation (ODE) integrators and the aforementioned numerical code,¹⁸ based on the GS orthonormalization procedure. For the forward time evolution we integrate the equations of motion for variables (ẋ, ẏ, ż) (eqn (1)–(3), for example), together with their corresponding linearized equations

with the Jacobian matrix

Eqn (10) governs the time evolution in the tangent space. Choosing Δx(0) = a⁽¹⁾₀, Δy(0) = a⁽²⁾₀ and Δz(0) = a⁽³⁾₀, as an initially orthonormal basis G₀ = (a⁽¹⁾₀,a⁽²⁾₀,a⁽³⁾₀), the evolution in tangent space up to a time t = nh (n = 0, 1, 2, 3,…) is then obtained from

G̃_n = J_nG₀, (12)

where J_n = J_n−1,J_n…J₀ is the product of the Jacobian matrix of the map evaluated for every orbital point at times t = nh. The subscript represents the number of iterations with time step h. It is implicit here that at the initial iteration (n = 0) we have already computed the orbit for a sufficiently long time (forward transient time) so that the orthonormal basis of the tangent space converged already to the asymptotic Gram–Schmidt vectors.¹¹ In the numerical code,¹⁸ there are three sets of linearized equations, one for G₀ = (1, 0, 0), one for G₀ = (0, 1, 0) and the last one for G₀ = (0, 0, 1). To avoid divergences, the matrix G̃_n is renormalized for each step. As usual, this can be done by the QR decomposition

G̃_n = G_nR_n. (13)

The matrix R_n is upper-triangular

and contains the information obtained in the GS orthonormalization procedure of G̃_n. Here, G_n = (a⁽¹⁾_n,a⁽²⁾_n,a⁽³⁾_n) are the GS vectors before the orthonormalization at iteration n and (g⁽¹⁾_n,g⁽²⁾_n,g⁽³⁾_n) are the GS vectors after the orthonormalization. Since g⁽¹⁾_n, g⁽²⁾_n and g⁽³⁾_n are orthonormal by construction, they can only provide information about the local rates of expansion (contraction) of the vectors, stored in the diagonal elements of R_n, from which the standard Lyapunov vectors can be obtained. The correct direction of the eigenvectors of the Jacobian matrix were corrupted by the orthonormalization procedure. They can be recovered by the backward integration, as shown below.

(ii) The essential property is that for the backward time evolution the unstable direction is the stable direction from the forward time evolution. Thus, in the backward evolution the Lyapunov vector related to the largest Lyapunov exponent points along the stable direction from the forward evolution. In this way we define the Lyapunov vectors as

v⁽¹⁾_n = c^(1,1)_ng⁽¹⁾_n, (15)

v⁽²⁾_n = c^(1,2)_ng⁽¹⁾_n + c^(2,2)_ng⁽²⁾_n, (16)

v⁽³⁾_n = c^(1,3)_ng⁽¹⁾_n + c^(2,3)_ng⁽²⁾_n + c^(3,3)_ng⁽³⁾_n. (17)

Since g⁽¹⁾_n, g⁽²⁾_n and g⁽³⁾_n were computed in the forward time evolution, the coefficients c^(j,i)_n can be determined from the dynamics in tangent space for the backward direction. In the matrix form we have C_j−1 = R_j⁻¹C_j, where R_j⁻¹ is the inverse of the R_j matrix. The time j starts to count after the backward transient, which is the time in the backward evolution sufficiently long to converge the tangent initial conditions close to the LVs. The LVs have normalized length so that the columns of C_j must be normalized to 1. The initial condition for C_n, before starting the backward evolution, can be a generic nonsingular upper triangular matrix, which is the GS basis. Finally, the angles between LVs are calculated from θ_ij = cos⁻¹(v⁽ⁱ⁾_n·v^(j)_n).

Since in this work we use a fourth-order Runge–Kutta integrator, the error in the variables is around h⁵, which is negligible for the time steps h that we use. As convenient rules of thumb, one has to chose time steps small enough compared to the relevant oscillations of the system. For example, the peaks in Fig. 4 occur for small intervals of time, around Δt ∼ 0.4. Thus, to compute many points along such a short peak, we used h = 5 × 10⁻⁵ so to have around 8000 for each peak (not all plotted here). Once the above compromises are satisfied, the predictions are independent of the time step. On the other hand, the number of steps is more related to how much computational time is available or, in our case, how many peaks are needed to have a reliable statistical analysis. The precision of Lyapunov vectors is more dependent on the transient times than on other parameters.

It is worth mentioning that in general there is no need to realize the backward time evolution. It is enough to determine the unstable direction from the forward evolution and compare it to the flow evolution direction. Below, we present in schematic form the main steps used to determine Lyapunov vectors.

Acknowledgements

The authors were supported by CNPq, Brazil. Open Access funding provided by the Max Planck Society.

References

1. For a survey see, e.g., N. Akhmediev, et al., Roadmap on optical rogue waves and extreme events, J. Opt., 2016, 18, 063001 and references therein.
2. E. Mercier, A. Even, E. Mirisola, D. Wolfersberger and M. Sciamanna, Numerical study of extreme events in a laser diode with optical feedback, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 042914.
3. Extreme Events in Nature and Society, ed. S. Albeverio, V. Jentsch, and H. Kantz, Springer, New York, 2006.
4. M. W. Beims and J. A. C. Gallas, Alignment of Lyapunov Vectors: A Quantitative Criterion to Predict Catastrophes?, Sci. Rep., 2016, 6, 37102.
5. M. W. Beims and J. A. C. Gallas, Manifold angles, the concept of self-similarity, and angle-enhanced bifurcation diagrams, Sci. Rep., 2016, 6, 18859.
6. M. R. Bassett and J. L. Hudson, Shilnikov chaos during copper electrodissolution, J. Phys. Chem., 1988, 92, 6963–6966.
7. P. Gaspard and G. Nicolis, What can we learn from homoclinic orbits in chaotic dynamics?, J. Stat. Phys., 1983, 31, 499.
8. J. A. C. Gallas, The structure of infinite periodic and chaotic hub cascades in phase diagrams of simple autonomous flows, Int. J. Bifurcation Chaos Appl. Sci. Eng., 2010, 20, 197–211.
9. V. Petrov, S. K. Scott and K. Showalter, Mixed-mode oscillations in chemical systems, J. Chem. Phys., 1992, 97, 6191.
10. Y. Pomeau, A. Pumir and P. Pelce, Intrinsic stochasticity with many degrees of freedom, J. Stat. Phys., 1984, 37, 39.
11. C. L. Wolfe and R. M. Samelson, An efficient method for recovering Lyapunov vectors from singular vectors, Tellus A, 2007, 59, 355.
12. D. Pazó, J. M. López and M. A. Rodriguez, On the angle between the first and the second Lyapunov vectors in spatio-temporal chaos, J. Phys. A: Math. Theor., 2013, 46, 254014.
13. A. Norwood, E. Kalnay, K. Ide, S. C. Yang and C. Wolfe, Lyapunov, singular and bred vectors in a multi-scale system an empirical exploration of vectors related to instabilities, J. Phys. A: Math. Theor., 2013, 46, 254021.
14. F. Ginelli, H. Chaté, R. Livi and A. Politi, Covariant Lyapunov vectors, J. Phys. A: Math. Theor., 2013, 46, 254005 , and references therein.
15. P. V. Kuptsov and U. Parlitz, Theory and computation of covariant Lyapunov vectors, J. Nonlinear Sci., 2012, 22, 727–762.
16. F. Argoul, A. Arneodo and P. Richetti, Experimental evidence for homoclinic chaos in the Belousov-Zhabotinsky reaction, Phys. Lett. A, 1987, 120, 269–275.
17. N. Sharafi, M. Timme and S. Hallerberg, Critical transitions and perturbation growth directions, Phys. Rev. E, 2017, 96, 032220.
18. A. Wolf, J. B. Swift, H. L. Swinney and J. A. Vastano, Determining Lyapunov exponents from a time series, Phys. D, 1985, 16, 285.

---