Lyapunov exponent of the logistic map (Mathematica Code)

In a previous post I’d shown a way to get the Lyapunov exponent from the time series data of any map. In this quick tutorial, I’ll show you a cleaner way to get the Lyapunov exponent for the specific case of the logistic map, and then using a really short script in Mathematica, plot it against r.

First the mathematical analysis that simplifies the expression for the Lyapunov exponent for a map, and particularly the logistic map. This discussion follows this article.

Suppose the initial infinitesimal perturbation is \delta x_{0}. Then we have, for n\rightarrow\infty:

\left|\delta x_{n}\right|=\left|\delta x_{0}\right|e^{\lambda n}\Rightarrow e^{\lambda n}=\underset{\delta x_{0}\rightarrow0}{\lim}\left|\frac{\delta x_{n}}{\delta x_{0}}\right|=\left|\frac{dx_{n}}{dx_{0}}\right|.

But we can write:

\frac{dx_{n}}{dx_{0}}=\frac{dx_{n}}{dx_{n-1}}\frac{dx_{n-1}}{dx_{n-2}}\ldots\frac{dx_{1}}{dx_{0}},

where each x_{i}=f(x_{i-1}). So we have:

\frac{dx_{n}}{dx_{0}}=f'(x_{n-1})f'(x_{n-2})\ldots f'(x_{0}).

Therefore,

\;e^{\lambda n}=\left|f'(x_{n-1})f'(x_{n-2})\ldots f'(x_{0})\right|

\Rightarrow\lambda(r;x_{0})=\underset{n\rightarrow\infty}{\lim}\frac{1}{n}\ln\left|f'(x_{n-1})f'(x_{n-2})\ldots f'(x_{0})\right|

=\underset{n\rightarrow\infty}{\lim}\frac{1}{n}\sum\limits_{k=0}^{n-1}\ln\left|f'(x_{k})\right|.

For the logistic map, f'(x)=r(1-2x).
So we have:

\lambda(r;x_{0})=\underset{n\rightarrow\infty}{\lim}\frac{1}{n}\sum\limits_{k=0}^{n-1}\ln\left|r(1-2x_{k})\right|.

We can put the above formula in a short Mathematica script to obtain \lambda as a function of r and plot it. The following is the code:

\[Lambda][r_] := Module[{f, l},
   f[x_] := r x (1 - x);
   l[x_] := Log[Abs[r (1 - 2 x)]];
   Mean[l[NestList[f, 0.1, 1*^2]]]];
Plot[\[Lambda][r], {r, 0, 4}, PlotStyle -> Thickness[.0001], 
 AxesLabel -> {"r", "\[Lambda](r)"}]

And the following is the output:
Lyapunov exponent for logistic map (Mathematica code)In the line that uses Nestlist, we specify the starting point of the trajectories. However, I noticed that the output does not depend on the starting point.

Advertisements

Calculating the Lyapunov Exponent of a Time Series (with python code)

(In a later post I discuss a cleaner way to calculate the Lyapunov exponent for maps and particularly the logistic map, along with Mathematica code.)
I found this method during my Masters while recreating the results of an interesting paper on how some standard tests for chaos fail to distinguish chaos from stochasticity (Stochastic neural network model for spontaneous bursting in hippocampal slices, Biswal and Dasgupta, 2002).
Procedure
The Lyapunov exponent is a measure of sensitive dependence on initial conditions, i.e. how quickly two nearby states diverge.
Now consider two points in the time-series, ti and tj, whose values are very close. That means the system reached near the same state at the ith and jth iterations. Now consider the two sequences  ti,  ti+1,  ti+2 … and  tj,  tj+1,  tj+2 … We wish to know how these two sequences diverge from each other. For this, consider the distance between the two sequences after k steps: d(k) = | ti+ktj+k |. (This is for a 1D time series. For higher dimensions, you can define this to be the Euclidean distance and modify the code accordingly.) If the system is chaotic, d(k) will initially rise exponentially with k. For this, one can plot ln d(k) vs k and apply a linear fit. The slope will be an estimate for the Lyapunov exponent.
(Since the system is bounded, the two nearby states will not diverge indefinitely though. Their exponential divergence will stop after some length. We must fit the straight line only within this region.)
Now, this was for a single pair of initial states. The Lyapunov exponent is an average of this divergence exponent over all nearby initial pairs. So for this, define d(k)>, where is averaging over all starting pairs  ti,  tj, such that the initial distance d(0) = | t– tj | is less than some fixed small value. The program finds all such initial pairs, calculates d(k)>, plots it against k, and the slope of the initial linear part gives us the Lyapunov exponent.
Python Code
The following code takes a text file with the time series, ‘timeseries.txt’, as the argument. The text file must contain only the time series values in a single column, no serial numbers or any other text before or after. It asks for the starting diameter within which to limit the initial pairs. It displays how many such pairs it is finding in the time series, so you can vary the diameter based on this.
It outputs a text file, ‘lyapunov.txt’ with two columns, k and d(k)>, which you can then plot and fit in the correct region by visual inspection.
from math import log

def d(series,i,j):
    return abs(series[i]-series[j])

f=open('timeseries.txt', 'r')
series=[float(i) for i in f.read().split()]
f.close()
N=len(series)
eps=input('Initial diameter bound: ')
dlist=[[] for i in range(N)]
n=0 #number of nearby pairs found
for i in range(N):
    for j in range(i+1,N):
        if d(series,i,j) < eps:
            n+=1
            print n
            for k in range(min(N-i,N-j)):
                dlist[k].append(log(d(series,i+k,j+k)))
f=open('lyapunov.txt','w')
for i in range(len(dlist)):
    if len(dlist[i]):
        print>>f, i, sum(dlist[i])/len(dlist[i])
f.close()

The following is the plot and fit of the resulting data from a logistic map series with an appropriately chosen initial diameter.

Lyapunov Exponent of Logistic Map

I deliberately did not automate the plotting and fitting part, because a. it’s tedious and hard to write the code in a way that runs on most installations, and b. human eyes will do a much more reliable job of identifying where the linear portion ends.