Linearized perturbation dynamics of Henon Attractor

Question

I am reading the following paper:

ergodic theory of chaos and strange attractors, by J.-P. Eckmann (can be easily downloaded)

My question is around equation (2.111), p. 624:

Given $a = 1.4$, $b = 0.3$. One can finds numerically that $$\delta x(t) \approx \delta x(0) e^{\lambda t}, \ \ \ \ \lambda = 0.42$$

I have question on how to obtain $e^{\lambda t}$ and $\lambda = 0.42$.

I think this comes from linearization. From the formula on the same page: $$\delta x(t) = (D_x f^t)\delta x(0).$$ I think $(D_x f^t) = (\partial f_i/\partial x_j)$, please see p.619, is a Jacobian matrix: $$(D_x f^t) = \begin{bmatrix}-2ax_1 & 1 \\ b & 0 \end{bmatrix}.$$ How to obtain $e^{\lambda t}$ and $\lambda = 0.42$?

Answer 1

Your assumption about $D_x f^t$ is wrong. Note that $f^t$ is defined as

$$ x(t) = f^t x(0) $$

that means that $f^t$ tells you the state of the system at any time $t$, and you do not know that. What you know is how to get $x(t + 1)$ from $x(t)$.

That being said $\delta(t)$ measures the separation between orbits. Imagine you start with two orbits very close to each other, but already in the attractor. Call the initial separation $\delta(0)$. If the system exhibits divergence of the initial conditions it is reasonable to assume that after some time $t$, the distance between the orbits will be modeled by something like

$$ \delta(t) \approx \delta(0)e^{\lambda t} $$

for some $\lambda > 0$. $\lambda$ is called a Lyapunov exponent, you can follow this link to learn how to calculate it.

Here's a very inefficient python code to calculate it, with this I calculated

$$ \lambda = 0.42220711662906674 $$

import numpy as np    
# henon map
def henon(x, **kwargs):

    x1 = 1 + x[1] - kwargs['a'] * x[0]**2
    y1 = kwargs['b'] * x[0]

    return np.array([x1, y1])


# iterates for a few reps
def transient(x, **kwargs):

    for k in range(100):
        x = henon(x, **kwargs)

    return x


# reinitializes the orbit
def reinit(x1, x2, d0, d1):

    xp = x1 + d0 * (x2 - x1) / d1

    return xp

def exponent(x0, eps = 1e-10, **kwargs):

    x1 = transient(x0, **kwargs)
    x2 = x1 + eps

    e = []

    for k in range(4000):

        # update
        d1 = np.linalg.norm(x2 - x1)

        x1 = henon(x1, **kwargs)
        x2 = henon(x2, **kwargs)

        d2 = np.linalg.norm(x2 - x1)

        # reinit
        x2 = reinit(x1, x2, d1, d2)

        e.append(np.log(d2 / d1))

    return np.mean(e)

print(exponent(np.array([0, 0]), a = 1.4, b = 0.3, eps = 1e-10))
>> 0.42220711662906674

Here's a plot of $\lambda$ as a function of $a$. You can see there are some region where $\lambda < 0$, in those cases paths do not diverge, so this plot is very useful to determine chaotic regions