Can anyone help me with calculating Lyapunov exponent of 2D map, for example Henon map?

If there is simple code, that would be great!

=================

=================

2 Answers

2

=================

As noted in Nasser’s answer to 17593, code to solve this problem already exists and can be downloaded from the web site of Dr Marco Sandri, who also wrote an article explaining his algorithms. Once downloaded, lce.m should be stored in some accessible place on the disk, for instance, the subdirectory MyPackages of Mathematica. The calculation proceeds as follows.

AppendTo[$Path, FileNameJoin[{$UserDocumentsDirectory, “Mathematica/MyPackages”}]];

Quiet@Needs[“lce`”]

Then, define the Henon recurrence function.

Î± = 1.4; Î² = .3;

h[{x_, y_}] := {1 – Î± x^2 + y, Î² x}

If desired, the Henon Map can be plotted.

lst = NestList[h, {.1, 0}, 40000];

ListPlot[lst, PlotStyle -> Red, AxesLabel -> {“x”, “y”}]

As noted in a comment above by thils, further discussion of the Henon Map is at Wolfram MathWorld. The Lyapunov exponent now can be found with

LCEsD[h, {0.1, 0}, 1000, 100, LCEsPlot -> False]

(* {{0.445647, -1.64962}, 1.27015} *)

the first two numbers being numerical estimates of the exponents, and the third the approximate dimension of the attractor. The convergence rate of the approximation can be obtained, if desired, from

LCEsD[h, {0.1, 0}, 1000, 100, LCEsPlot -> True]

Addendum

Lyapunov exponents for a range of parameters can be computed quickly, for instance,

ly = Quiet@Table[h[{x_, y_}] := {1 – Î± x^2 + y, Î² x};

{Î±, LCEsD[h, {0.1, 0}, 1000, 100, LCEsPlot -> False][[1, 1]]}, {Î±, 1.0, 1.4, 0.005}]

ListPlot[ly, AxesLabel -> {“Î±”, None}]

Negative exponents correspond to stable windows within an otherwise chaotic band. (LCEsD is unable to compute the dimension of the attractor, which is zero in these cases, and instead generates a number of error messages. Quiet is used to suppress them.) As an example, the Henon Map for Î± = 1.23 is finite number of discrete points, representing a normal, although complicated, attractor. (PointSize[.01] is used in this plot to improve visibility of the points.)

This paper provides a reasonable background to evaluation of the Lyaponov exponent. Description of the Henon map is also provided. For the case of the Henon map, we include two parameters (Î±, Î²), and we can start by using

n = 100;

max = 20;

r = 0.005;

Manipulate[

ListPlot[Flatten[

Table[({Î±, #1} & ) /@

Take[Quiet[

NestList[{#1[[2]] +

1 – Î±*#1[[1]]^2, Î²*#1[[1]]} & , {-0.1, 0.},

n][[All, 1]]], -max], {Î±, 0.9, 1.5, r}], 1]],

{Î², 0.1, 0.25}]

n is the no. of iterations, r is resolution of the parameter Î± and max is the maximum period. All these variables determine the computational time.

This doesn’t seem to address the OP’s question other than by providing an link to another site. The OP didn’t ask for code that draws the bifurcation diagram.

– m_goldberg

Oct 4 ’15 at 13:06

@m_Goldberg I ran out time, hopefully will attach the 2nd part which shows explicit calc of lyapunov exponent which can compare with the Ice.m package (@bbgodfrey response).

– thils

Oct 4 ’15 at 19:44