computational physics - Ljapunov exponent of driven damped pendulum

I have written a computer simulation of the driven damped pendulum, pretty much as the one shown here, only that I did it Python. Next, I have found some parameters for which the pendulum behaves chaotically. Now I want to extract the Lyapunov exponent from the system. To do this, I let the system run two times, one time at the found parameters for chaos, a second time with the same parameters plus a very small deviation in the initial position of the pendulum. Then I plot the difference in the displacement of the two pendulums as a function of time, as shown in the figure below (taken from Giordano, Computational Physics).

I can qualitatively extract the Lyapunov exponent, by doing this several times and fitting a line to the data.

The question is, however, how would I do this quantitatively or even analytically?

Answer

Let's start out using notation similar to the example you linked to:

$$ \ddot{y}+b\dot{y}+\sin(y)=a\cos(ct)+d $$

As in the example, we'll write this in autonomous form:

$$ \mathbf{x}=\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\begin{pmatrix}y\\\dot{y}\\ct\end{pmatrix} $$

the differential equation can now be written:

$$ \mathbf{\dot{x}}=\begin{pmatrix}\dot{x_1}\\\dot{x_2}\\\dot{x_3}\end{pmatrix}=\mathbf{f(x)}=\begin{pmatrix}x_2\\-bx_2-\sin{x_1}+a\cos{x_3}+d\\c\end{pmatrix} $$

Now, suppose we study the evolution of two different starting points in this system that are arbitraily close. Let's denote the trajectory of one of them with $\mathbf{x}(t)$ and the other one with $\mathbf{x}(t)+\mathbf{z}(t)$, where $\mathbf{z}(t)$ is the vector difference between the two points. The time derivative of $\mathbf{z}(t)$ can be written:

$$ \mathbf{\dot{z}}=\mathbf{f(x+z)}-\mathbf{f(x)} $$

If $\mathbf{z}$ is sufficiently small this can be written:

$$ \mathbf{\dot{z}}=\mathbf{Jz} $$

where $\mathbf{J(x)}$ is the Jacobian of $\mathbf{f(x)}$:

$$ \mathbf{J}=\begin{pmatrix} 0 & 1 & 0\\ -\cos{x_1} & -b & -a\sin{x_3}\\ 0 & 0 & 0 \end{pmatrix} $$

Over a short period $\delta t$, the change in $\mathbf{z}$ is given by:

$$ \mathbf{z}(t+\delta t)=\mathbf{z}(t)+\delta t\mathbf{Jz}(t)=(\mathbf{I} +\delta t\mathbf{J})\mathbf{z}(t) $$

Assuming $\mathbf{z}$ and $\delta t$ are small enough for the linearizations to hold we can decribe the change in $\mathbf{z}$ from time $t_1$ to $t_2$ with:

$$ \mathbf{z}(t_2)=\left[\prod_i(\mathbf{I} +\delta t_i\mathbf{J}(\mathbf{x}_i))\right]\mathbf{z}(t_1) $$

where $\delta t_i$ are sufficiently small time intervals from $t_1$ to $t_2$ and $\mathbf{x}_i$ is $\mathbf{x}$ at those times.

What you need to do to calculate the Lyaponov exponents of the this system is to numerically solve the trajectory $\mathbf{x}(t)$ for an arbitrary starting point $\mathbf{x}_0$ over a long time interval $T$. Then you need to calculate the matrix $\mathbf{A}$:

$$ \mathbf{A}=\prod_i(\mathbf{I} +\delta t_i\mathbf{J}(\mathbf{x}_i)) $$

using the values for $\delta t_i$ and $\mathbf{x}_i$ from the solved trajectory. You then need to calculate the eigenvalues of $\mathbf{A}$. Let's call them $\alpha_1$, $\alpha_2$ and $\alpha_3$. The lyaponov exponents are given by:

$$ \lambda_i=\frac{1}{T}\log{\alpha_i} $$

assuming $T$ is long enough. You need to experiment with different $\mathbf{x}_0$ and $T$. Hopefully, the variations in the Lyaponov exponents you calculate for different $\mathbf{x}_0$ will decrease with larger $T$ and converge for sufficiently large $T$.

One problem you will almost certainly run into is that the values of $\mathbf{A}$ will become too big to handle numerically. The way to get around this is to keep track of the magnitude of $\mathbf{A}$ as it is iteratively calculated and normalize it whenever its highest value reaches some threshold. The normalizing factors needs to be saved each time this is done. Assuming this has been done you will in the end have a matrix $\mathbf{B}$ and a series of normalizing factors $N_j$ (in my notation $N_j$ are the values you have divided by when normalizing) relating to $\mathbf{A}$ as:

$$ \mathbf{A}=\left(\prod_j N_j\right)\mathbf{B} $$

If $\beta_i$ are the eigenvalues of $\mathbf{B}$, the Lyaponov exponents can be calculated as:

$$ \lambda_i=\frac{1}{T}\left(\log{\beta_i}+\sum_j\log{N_j}\right) $$

It is usually the largest of the Lyaponov exponents that is of interest as it will dominate the exponential growth in $\mathbf{z}$. At least one of the Lyaponov exponents must be positive for the system to be chaotic.

It is important to understand that $\mathbf{z}$ should be thought of as a vector that remains very small throughout the whole period $T$. This might be a bit counter-intuitive since it grows exponentially. However, for any finite $T$, the starting value of $\mathbf{z}$ can always be chosen small enough to make sure it reamins sufficiently small thoughout the time period.