Can anyone help me with calculating Lyapunov exponent of 2D map, for example Henon map? If there is simple code, that would be great!
-
1Can you explain what it is or provide a link...? – Marius Ladegård Meyer Oct 03 '15 at 22:02
-
http://mathworld.wolfram.com/HenonMap.html may be of help – thils Oct 03 '15 at 22:16
-
Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory [tour] now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. – bbgodfrey Oct 03 '15 at 23:55
-
Search Google for "Lyapunov exponent Mathematica" to find many references. This one.is a possibility. – bbgodfrey Oct 04 '15 at 00:01
-
1Consider also this. – bbgodfrey Oct 04 '15 at 00:10
2 Answers
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.
- 107,779
- 16
- 103
- 257
- 3,228
- 2
- 18
- 35
-
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 04 '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 04 '15 at 19:44




