Determining the Spectrum of the Nonlinear Local Lyapunov Exponents in a Multidimensional Chaotic System

The estimation of the average rates of divergence (or convergence) of initially nearby trajectories in phase space has been studied with (global) Lyapunov exponents, which are used to quantify the average predictability properties of a chaotic system (Eckmann and Ruelle, 1985; Wolf et al., 1985). The sum of all the positive Lyapunov exponents is an estimate of the Kolmogorov entropy K, and the inverse of K is a measure of the total predictability of the system (Kolmogorov, 1941; Fraedrich, 1987, 1988). Considering that the (global) Lyapunov exponents only provide a measure of the total predictability of a system, various local or finite-time Lyapunov exponents (Nese, 1989; Houtekamer, 1991; Yoden and Nomura, 1993; Ziehmann et al., 2000) have been subsequently proposed to measure the local predictability around a point x₀ in phase space. However, the existing global or local Lyapunov exponents have limitations because they all satisfy the assumption that the initial perturbations are sufficiently small and the tangent linear model (TLM) of a nonlinear system could approximately govern their evolution (Lacarra and Talagrand, 1988; Feng and Dong, 2003; Mu and Duan, 2003; Duan and Mu, 2009). If an initial perturbation is large enough to invalidate the TLM, it is no longer possible to apply the existing global or local Lyapunov exponents in predictability studies of chaotic systems (Kalnay and Toth, 1995).

In view of the limitations of the existing global or local Lyapunov exponents, (Ding and Li, 2007) introduced the concept of the nonlinear local Lyapunov exponent (NLLE). The NLLE measures the average growth rate of the initial errors of nonlinear dynamical models without linearizing the governing equations. The experimental results of (Ding and Li, 2007) show that, compared with a linear local or finite-time Lyapunov exponent, the NLLE is more appropriate for the quantitative determination of the predictability limit of a chaotic system. Based on observational or reanalysis data, the NLLE method has been used to investigate the atmospheric predictability at various timescales (Ding et al., 2008, 2010, 2011, 2016; Li and Ding, 2008, 2011, 2013).

However, recall that the NLLE defined by (Ding and Li, 2007) only characterizes the nonlinear growth rate of the initial perturbations along the fastest growing direction, which is insufficient to describe the expanding or contracting nature of the initial perturbations along different directions in phase space. This explains why the Lyapunov exponent spectrum, rather than only the largest Lyapunov exponent, was introduced in the traditional Lyapunov theory. Therefore, for an n-dimensional chaotic system, it is necessary to extend the definition of the NLLE from one- to n-dimensional spectra, as this would allow us to investigate nonlinear evolution behaviors of initial perturbations along different directions in phase space.

In this paper, we first introduce the definition of the NLLE spectrum, and then propose a method to compute the NLLE spectrum. Finally, we demonstrate the validation and usefulness of the NLLE spectrum in characterizing the nonlinear evolutionary behaviors of initial perturbations along different directions, and in measuring the predictability limit of chaotic systems, by applying it to three chaotic systems.

For an n-dimensional nonlinear dynamical system, each initial error vector tends to fall along the fastest growing direction (Wolf et al., 1985). Therefore, for any initial error vector δ'(t₀), after a short time τ, the error vector δ'(t₀+τ) will capture the fastest growing direction. If this short time τ is not taken into consideration, instead, taking δ'(t₀+τ) as the initial error δ(t₀), the first NLLE along the fastest growing direction can be defined as \begin{equation} \label{eq1} \lambda_1({x}(t_0),{\delta}(t_0),\tau)=\frac{1}{\tau}\ln\frac{\|{\delta}(t_0+\tau)\|}{\|{\delta}(t_0)\|} ,\ \ (1) \end{equation} where Λ₁(x(t₀),δ(t₀),τ) depends on the initial state x(t₀) in phase space, the magnitude of the initial error vector δ(t₀)=\|δ(t₀)\|, and the evolution time τ (Ding and Li, 2007). Note that the direction of the fastest growing initial error vector δ(t₀) is determined by the local state x(t₀) in phase space, and Λ₁ is therefore independent of the direction of δ(t₀) and dependent only on its magnitude δ(t₀).

Once the first NLLE, Λ₁(x(t₀),δ(t₀),τ), has been obtained, the mth NLLE along the mth fastest growing direction can be successively determined from the growth rate of the volume V_m of m-dimensional subspace spanned by the orthogonal vector set of initial errors \(\pmb\Omega_m(t_0)=(\delta_1(t_0),\delta_2(t_0),\cdots\cdots,\delta_m(t_0))\) (m=2,3,……,n): \begin{equation} \sum_{i=1}^m\lambda_i=\frac{1}{\tau}\ln\frac{V_m({\pmb{\Omega}}_m(t_0+\tau))}{V_m({\pmb{\Omega}}_m(t_0))} , (m=2,3,\cdots\cdots,n) ,\ \ (2) \end{equation} where Λ_i=Λ_i(x(t₀),δ(t₀),τ) is defined as the ith NLLE of the dynamical system; δ₁(t₀) is equal to δ(t₀) in Eq. (2), which captures the local direction of the most rapid growth; and δ_i(t₀) (i=2,3,……,m) is the ith fastest growing initial error vector. These initial error vectors, δ₁(t₀),δ₂(t₀),……,δ_m(t₀), have the same magnitude δ(t₀), but they are orthogonal to each other. After δ₁(t₀) is obtained, δ_i(t₀) (i=2,3,……,m) can be successively determined by comparing their error growth rates.

It has already been pointed out that each error vector in a chaotic system generally converges towards the local direction of most rapid growth (Wolf et al., 1985). Due to the finite precision of computer calculations, the collapse towards a common direction causes the orientation of all error vectors to become indistinguishable. As a result, the calculation of the volume V_m of m-dimensional subspace becomes extremely difficult. This problem can be overcome by the repeated use of the Gram-Schmidt reorthogonalization (GSR) procedure on the vector frame (Björck, 1967, 1994). In the GSR procedure, the first vector, whose direction is never affected, tends to freely seek out the fastest growing direction (vectors along other directions are either growing less rapidly or shrinking); the second vector is projected onto the subspace orthogonal to the first vector and, iteratively, the nth vector is projected onto the subspace orthogonal to the previous n-1 vector. That is, given a set of vectors {v₁,v₂,……,v_n}, the GSR provides the following orthogonal set {v'₁,v'₂,……,v'_n}: \begin{eqnarray} {v}'_1&=&{v}_1 ;\nonumber\\ {v}'_2&=&{v}_2-\frac{({v}_2,{v}'_1)}{({v}'_1,{v}'_1)}{v}'_1 ;\nonumber\\ &&\vdots\nonumber\\ \label{eq2} {v}'_n&=&{v}_n-\frac{({v}_n,{v}'_{n-1})}{({v}'_{n-1},{v}'_{n-1})}{v}'_{n-1}-\cdots\cdots- \frac{({v}_n,{v}'_1)}{({v}'_1,{v}'_1)}{v}'_1 . \ \ (3)\end{eqnarray}

Note that error vectors in the GSR are not rescaled (i.e., dividing the vectors by their magnitudes) because of the dependence of the nonlinear growth of error vectors on their magnitudes. This is an important difference between the calculations of nonlinear and linear Lyapunov exponent spectra. The linear Lyapunov exponent spectral calculation requires the reorthonormalization of error vectors at each step of the GSR procedure (Wolf et al., 1985). It can be seen that the subspace spanned by the first m vectors is unaffected by GSR, so that the volume defined by these m vectors is not changed before and after the GSR procedure. In practice, as v'₁,v'₂,……,v'_n are orthogonal to each other, representing the vectors along the directions from most rapid growth to most rapid decline, we may determine the ith NLLE, Λ_i (i=1,2,……,n), directly from the growth rate of vector v'_i relative to the initial vector v_i. In addition to the GSR method (also referred to as the Wolf method), the spectrum of traditional Lyapunov exponents can be computed using the singular value decomposition (SVD) of the Jacobian matrix (referred to the Jacobian method) (Barana and Tsuda, 1993). However, because the error growth equations governing the NLLE are not linearized, the NLLE spectrum cannot be computed using the SVD of the Jacobian matrix. This is another difference between the calculations of nonlinear and linear Lyapunov exponent spectra.

Figure 1.  The Lorenz system’s mean NLLEs \(\bar\lambda_i\) [(a) i=1, (b) i=2, (c) i=3], and the logarithm of the corresponding mean RGIEs \(\ln\overline{E}_i\) [(d) i=1, (e) i=2, (f) i=3], as a function of time τ. The magnitudes δ(t₀) of the initial error vectors are all 10^-6, but their directions are orthogonal to each other.

To study the dynamic characteristics of the whole system, the ensemble average of the ith NLLE, \(\bar\lambda_i(\delta(t_0),\tau)=\langle\lambda_i(x(t_0),\delta(t_0),\tau)\rangle\) (i=1,2,……,n, and the symbol \(\langle\ \ \rangle\) denotes the ensemble average over a great number of different initial states) should be introduced, and then the ith mean relative growth of initial error (RGIE), \(\overline{E}_i(\delta(t_0),\tau)=\exp(\bar\lambda_i(\delta(t_0),\tau)\tau)\), can be obtained. In light of dynamical systems theory and probability theory, (Ding and Li, 2007) proved a saturation theorem of \(\overline{E}_1(\delta(t_0),\tau)\); that is, for a chaotic system, \(\overline{E}_1(\delta(t_0),\tau)\) will necessarily reach a saturation value in a finite time interval. According to this theorem, the average predictability limit of an error vector along the fastest growing direction could be quantitatively determined as the time at which \(\overline{E}_1(\delta(t_0),\tau)\) reaches its saturation level. In addition to \(\overline{E}_1(\delta(t_0),\tau)\) (corresponding to the fastest growing directions), it would be possible for several \(\overline{E}_i(\delta(t_0),\tau)\) (i≥ 2, corresponding to other growing directions) to finally reach their respective saturation levels in a high-dimensional chaotic system. Therefore, we can also quantitatively determine the predictability limits of these several \(\overline{E}_i(\delta(t_0),\tau)\) according to their respective saturation levels.

Figure 2.  The hyperchaotic Lorenz system’s mean NLLEs \(\bar\lambda_i\) [(a) i=1, (b) i=2, (c) i=3, (d) i=4], and the logarithm of the corresponding mean RGIEs \(\ln\overline{E}_i\) [(e) i=1, (f) i=2, (g) i=3, (h) i=4], as a function of time τ. The magnitudes δ(t₀) of the initial error vectors are all 10^-6, but their directions are orthogonal to each other. Due to the finite precision of computer calculations, the error along the most rapidly shrinking direction will become indistinguishable, and so the possible evolutions of \(\bar\lambda_4\) and \(\ln\overline{E}_4\) are denoted by the dashed lines in (d) and (h), respectively.

Once the NLLE spectrum and corresponding RGIEs are obtained, we can examine the evolutionary behaviors of the NLLE spectrum and corresponding RGIEs in a multidimensional chaotic system. Our method for computing the NLLE spectrum is tested on three dynamical systems with different complexity: the 3-variable Lorenz system (Lorenz, 1963), the 4-variable hyperchaotic Lorenz system (Li et al., 2005; Wang and Liu, 2006), and the 40-variable Lorenz96 model (Lorenz, 1996). The 3-variable Lorenz system is \begin{equation} \label{eq3} \left\{ \begin{array}{l} \dot{X}=-\sigma X+\sigma Y \\ \dot{Y}=rX-Y-XZ \\ \dot{Z}=XY-bZ \end{array} \right. , \ \ (4)\end{equation} where X, Y, Z are state variables, σ=10, r=28, and b=8/3, for which the well-known butterfly attractor exists. Figure 1 shows the Lorenz system’s \(\bar\lambda_i\) and \(\ln\overline{E}_i\) (i=1,2,3) as a function of time τ. As can be seen, \(\bar\lambda_1\) initially remains constant at around 0.91, then decreases rapidly, and finally approaches zero as τ increases (Fig. 1a). (Ding and Li, 2007) proved that \(\bar\lambda_1\) decreases asymptotically in a similar manner to O(1/τ) as τ tends to infinity. Correspondingly, \(\overline{E}_1\) initially increases exponentially, then its growth rate slows, and finally stops increasing and reaches a saturation value (Fig. 1d). \(\bar\lambda_2\) initially remains constant near zero, and then decreases gradually as τ increases (Fig. 1b). Accordingly, \(\overline{E}_2\) initially remains almost constant, and then decreases gradually towards zero (Fig. 1e). \(\bar\lambda_3\) initially moves rapidly towards a large negative value (-14.5), and subsequently remains almost constant at this value (Fig. 1c). Correspondingly, \(\overline{E}_3\) decreases exponentially from the beginning and rapidly approaches zero (Fig. 1f).

Three Lyapunov exponents of the Lorenz system under given parameters are 0.906, 0, and -14.572 (Wolf et al., 1985). Note that these values of the Lyapunov exponents are close to those of the NLLEs at the initial phase (i.e., the phase of the NLLEs remaining almost constant) in the Lorenz system (Figs. 1a-c), indicating that when the errors are small enough to validate the TLM, the error growth rates measured by the NLLE spectrum are close to the Lyapunov exponents. However, as the errors increase, error evolution enters the nonlinear phase, the TLM is no longer valid, and the error growth rates measured by the NLLE spectrum show different behaviors from those during the linear phase. These results indicate that the NLLE spectrum may realistically reflect the time-varying characteristics of error growth rates along different directions from the linear to nonlinear phases of error growth. This represents an improvement over the traditional Lyapunov exponent spectrum, which only characterizes the error growth rates during the linear phase of error growth.

The hyperchaotic Lorenz system is formulated by introducing an additional state variable into the 3-variable Lorenz system. The equations that describe the hyperchaotic Lorenz system are \begin{equation} \label{eq4} \left\{ \begin{array}{l} \dot{x}_1=a(x_2-x_1)\\ \dot{x}_2=bx_1 +cx_2-x_1 x_3 +x_4\\ \dot{x}_3=-dx_3 +x_1x_2\\ \dot{x}_4=-kx_1 \end{array} \right. ,\ \ (5) \end{equation} where x_i (i=1,2,3,4) are state variables, a, b, c, d, and k are system parameters. Here, taking a=35, b=7, c=12, d=3, and k=5, the hyperchaotic attractor exists (Li et al., 2005; Wang and Liu, 2006). As can be seen from Fig. 2, \(\bar\lambda_1\), \(\bar\lambda_3\), and \(\bar\lambda_4\) of the hyperchaotic Lorenz system have similar time-varying characteristics to \(\bar\lambda_1\), \(\bar\lambda_3\), and \(\bar\lambda_4\), respectively, of the Lorenz system. Therefore, a detailed description of the time evolution of \(\bar\lambda_i\) and \(\overline{E}_i\) (i=1,3,4) of the hyperchaotic Lorenz system is not given here. Next, we place more emphasis on the evolutions of \(\bar\lambda_2\) and \(\overline{E}_2\), and their influences on the predictability estimate of the hyperchaotic Lorenz system.

In the hyperchaotic Lorenz system, \(\bar\lambda_2\) initially remains a positive constant, and then decreases gradually and changes from positive to negative as τ increases (Fig. 2b). Correspondingly, \(\overline{E}_2\) initially increases, and then decreases gradually after reaching its maximum value (Fig. 2f). Because the maximum value of \(\overline{E}_2\) is far below the saturation value of \(\overline{E}_1\), \(\overline{E}_2\) no longer plays an important role in the error growth of the hyperchaotic Lorenz system. As mentioned earlier, the traditional Lyapunov theory uses the inverse of the sum of all the positive Lyapunov exponents (an estimate of the Kolmogorov entropy) as a measure of the total predictability of chaotic systems (Kolmogorov, 1941; Fraedrich, 1987, 1988). In this case, the second Lyapunov exponent will play an important role in limiting the predictability of the hyperchaotic Lorenz system (the traditional Lyapunov exponents of the hyperchaotic Lorenz system under given parameters are 0.41, 0.20, 0.00, and -26.7). These results indicate that the traditional Lyapunov exponents, based on linear error dynamics, may be insufficient to estimate the predictability time because they only consider the linear phase of error growth, without considering the nonlinear phase during which error growth either slows down or decreases gradually after the error initially increases exponentially. From this point of view, the NLLE spectrum is more suitable for characterizing the predictability of chaotic systems.

Our method for computing the NLLE spectrum is also applied to the Lorenz96 model. The model has 40 state variables, X₁,X₂,……,X₄₀, which are governed by the equation dX_i/dt=(X_i+1-X_i-2)X_i-1-X_i+F, where the index 1≤ i≤ 40 is arranged cyclically, and F is a fixed forcing. When F=8.0, the model displays sensitive dependence on the initial conditions (Lorenz, 1996). The Lorenz96 model has been used as a low-order proxy for atmospheric prediction and assimilation studies. Similar to the cases of the 3-variable and 4-variable Lorenz systems presented above, the NLLE spectrum of the Lorenz96 model initially remains close to its Lyapunov exponent spectrum, but depart from each other afterwards as error growth enters the nonlinear phase (Fig. 3a). Note that in Fig. 3a only the first 12 NLLEs [\(\bar\lambda_i(\delta(t_0),\tau)\), i≤ 12] are shown; the remaining NLLEs, \(\bar\lambda_i(\delta(t_0),\tau)\), are not shown due to space limitations. The first 12 RGIEs [i.e., \(\overline{E}_i(\delta(t_0),\tau)\), i≤ 12] finally reach saturation (Fig. 3b), and the rest, \(\overline{E}_i(\delta(t_0),\tau)\) (13≤ i≤ 40), decrease gradually after reaching a maximum value, or decrease directly from the beginning (not shown).

Figure 3.  The Lorenz96 model’s (a) mean NLLEs \(\bar\lambda_i\) (i≤ 12), and (b) the logarithm of the corresponding mean RGIEs \(\ln\overline{E}_i\) (i≤ 12), as a function of time τ. The magnitudes δ(t₀) of the initial error vectors are all 10^-6, but their directions are orthogonal to each other.

Figure 4.  Mean NLLEs \(\bar\lambda_i\) (i=2,4,6,8,10,12) of the Lorenz96 model, as a function of time τ, obtained from noise-free (red line) and noise-contaminated (blue and green lines) data. The magnitudes of initial error vectors are 10^-6, and the magnitudes of Gaussian white noise are 10^-7 (blue line) and 10^-6 (green line), respectively.

Following the work of (Dalcher and Kalnay, 1987), we determine the predictability limit as the time at which the error reaches 98% of its saturation level. The predictability limits determined from the curves of these 12 RGIEs are shown in Fig. 4. the results show that a more rapid initial growth of error vectors generally corresponds to a lower predictability limit, and the predictability limit of \(\overline{E}_i(\delta(t_0),\tau)\) (i≤ 12) increases more and more quickly with increasing i. As a result, the predictability limits show significant variations among error vectors of different directions, ranging from 11.5 to 50.1, suggesting that the predictability limit is highly sensitive to the direction of the initial error vector in the Lorenz96 model. These results further indicate that the NLLE method can effectively separate the slowly and rapidly growing perturbations, which is very important for studies of predictability and error growth dynamics.

The three examples presented above are all noise-free dynamical systems. However, noise is inevitably present in experimental and natural systems. The effect of noise on the estimation of the linear Lyapunov exponent spectrum has been noted in previous studies (e.g., Wolf et al., 1985; Brown et al., 1991). It is of interest and importance to explore the effects of noise on the calculation of the NLLE spectrum. For this purpose, we add measurement noise to the Lorenz96 model; that is, a small Gaussian white noise is added to each element of the time series of 40 variables of the Lorenz96 model after the entire series of all 40 variables have been generated. The NLLE spectrum is then computed based on these noise-contaminated data (Fig. 5). For noise of relatively small amplitude (magnitude of noise smaller than that of the initial error), the NLLE spectrum is not greatly affected by the noise. In contrast, for noise of sufficiently large amplitude (magnitude of noise close to that of the initial error), the NLLE spectrum cannot be accurately determined. The results suggest that the effects of noise should be considered for an accurate estimation of the NLLE spectrum. (Wolf et al., 1985) pointed out that low-pass filtering may be a feasible approach to reduce the effects of noise, which provides some direction for improving the estimation of the NLLE spectrum in the presence of noise.

Figure 5.  Predictability limits (dimensionless) determined from the curves of the first 12 RGIEs [i.e., \(\overline{E}_i(\delta(t_0),\tau)\), i≤12] in the Lorenz96 model. The magnitudes δ(t₀) of initial error vectors are set to 10^-6.

To conclude, for an n-dimensional chaotic system, we extend the definition of the NLLE from one- to n-dimensional spectra. Our experimental results from three dynamical systems with different complexity show that the NLLEs converge to the Lyapunov exponents during the linear phase of error growth, but they deviate from the Lyapunov exponents during the nonlinear phase of error growth. This time dependence of the NLLE spectrum realistically characterizes the growth rates of initial error vectors along different directions from the linear to nonlinear phases of error growth. This represents an improvement over the traditional Lyapunov exponent spectrum, which only characterizes the error growth rates during the linear phase of error growth. In addition, because the NLLE spectrum can effectively separate the slowly and rapidly growing perturbations, it is shown to be more suitable for estimating the predictability of chaotic systems, as compared to the traditional Lyapunov exponent spectrum.

Our results present a preliminary application of the NLLE spectrum in studies of error growth and predictability in several relatively simple systems. For a more complex system, such as weather and climate, there are many types of instabilities (Toth and Kalnay, 1993, 1997; Trevisan and Legnani, 1995; Norwood et al., 2013). These instabilities generally show different error growth rates and thereby have different predictability times. It would be interesting to extend the investigation to more realistic weather and climate systems, which we intend to do in future research.