Henon-Heiles Model

These equations model a star moving in the field of its galaxy. The model possesses chaotic orbits. The dynamical system is continuous, and is governed by Hamilton's equations from physics for the motion of a particle.

Clear[x, y, px, py]
Hamiltonian[x_, y_, px_, 
        py_]  =  (px^2 + py^2 + x^2 + y^2)/
          2 + (x^2  y - y^3/3)

Clear[x, y, px, py]
eq1  =  y'[t]  == 
        D[Hamiltonian[x[t], y[t], px[t], py[t]], 
          py[t]]
eq2  =  x'[t]  == 
        D[Hamiltonian[x[t], y[t], px[t], py[t]], 
          px[t]]
eq3  =  px'[t]  == - 
          D[Hamiltonian[x[t], y[t], px[t], py[t]], 
            x[t]]
eq4  =  py'[t]  == - 
          D[Hamiltonian[x[t], y[t], px[t], py[t]], y[t]]

Mathematica is unable to solve this system analytically.

DSolve[{eq1, eq2, eq3, eq4}, {px[t], py[t], y[t], x[t]}, t]

However, we can approximate it using Euler's method, and then proceed to analyze the system.

f1[x_, y_, px_, py_]  =  
      D[Hamiltonian[x, y, px, py], px]
f2[x_, y_, px_, py_]  =  
      D[Hamiltonian[x, y, px, py], py]
f3[x_, y_, px_, py_]  = - 
        D[Hamiltonian[x, y, px, py], x]
f4[x_, y_, px_, py_]  = - 
        D[Hamiltonian[x, y, px, py], y]

Clear[x, y, px, py]
x[0] :=  0.1
y[0] :=  0.05
px[0] :=  0.0
py[0] :=  0.11

    
num = 50000;
to  =  0;
tf  =  100;
h  =  
      N[(tf - to)/num]

x[n_] :=  x[n]  =  
        h f1[x[n - 1], y[n - 1], px[n - 1], py[n - 1]] + 
          x[n - 1]
y[n_] :=  y[n]  =  
        h f2[x[n - 1], y[n - 1], px[n - 1], py[n - 1]] + 
          y[n - 1]
px[n_] :=  px[n]  =  
        h f3[x[n - 1], y[n - 1], px[n - 1], py[n - 1]] + 
          px[n - 1]
\(py[n_] := \ \(py[n]\  = \ 
        h\ f4[x[n - 1], y[n - 1], px[n - 1], py[n - 1]] + 
          py[n - 1]\)
\)
data  =  
        Table[{h*n, x[n], y[n], px[n], py[n]}, {n, 0, 
            num}];

xyphase  =  
        Table[{Part[data, n + 1, 2], Part[data, n + 1, 3]}, {n, 0, 
            num}];
phaseplot = 
      ListPlot[xyphase, 
        PlotRange -> {{-0.15, 0.15}, {-0.15, 0.15}}, 
        PlotJoined -> True, Frame -> True, 
        FrameLabel -> {"x", "y"}, 
        PlotStyle -> {RGBColor[1, 0, 
              0]}]

xtime  =  
        Table[{Part[data, n + 1, 1], Part[data, n + 1, 2]}, {n, 0, 
            num}];
xplot = 
      ListPlot[xtime, PlotRange -> All, PlotJoined -> True, 
        Frame -> True, FrameLabel -> {"Time", "x"}, 
        PlotStyle -> {RGBColor[1, 0, 0]}]

Here is another trajectory in the Henon Heiles model:

Clear[x, y, px, py]
x[0] :=  0.2
y[0] :=  0.1
px[0] :=  0.04
\(\(\(py[0] := \ \(-0.05\)\)\(
\)
    \)\)
num = 50000;
to  =  0;
tf  =  140;
h  =  N[(tf - to)/num]

    
x[n_] :=  x[n]  =  
        h f1[x[n - 1], y[n - 1], px[n - 1], py[n - 1]] + 
          x[n - 1]
y[n_] :=  y[n]  =  
        h f2[x[n - 1], y[n - 1], px[n - 1], py[n - 1]] + 
          y[n - 1]
px[n_] :=  px[n]  =  
        h f3[x[n - 1], y[n - 1], px[n - 1], py[n - 1]] + 
          px[n - 1]
\(\(\(py[n_] := \ \(py[n]\  = \ 
        h\ f4[x[n - 1], y[n - 1], px[n - 1], py[n - 1]] + py[n - 1]\)\)\(
\)
    \)\)
data  =  
        Table[{h*n, x[n], y[n], px[n], py[n]}, {n, 0, 
            num}];

    
xyphase  =  
        Table[{Part[data, n + 1, 2], Part[data, n + 1, 3]}, {n, 0, 
            num}];
phaseplot1 = 
      ListPlot[xyphase, PlotRange -> {{-0.4, 0.4}, {-0.4, 0.4}}, 
        PlotJoined -> True, Frame -> True, 
        FrameLabel -> {"x", "y"}, 
        PlotStyle -> {RGBColor[1, 0, 0]}]

    
xtime  =  
        Table[{Part[data, n + 1, 1], Part[data, n + 1, 2]}, {n, 0, 
            num}];
xplot1 = 
      ListPlot[xtime, PlotRange -> All, PlotJoined -> True, 
        Frame -> True, FrameLabel -> {"Time", "x"}, 
        PlotStyle -> {RGBColor[1, 0, 0]}]

The Rossler Attractor

This attractor I will solve numerically using NDSolve.

Clear[x, y, z, f1, f2, f3, t, c, sol]

    
tf  =  210;
sol  = 
      NDSolve[{
x'[t] == -(y[t] + z[t]), 
          x[0] == -0.04, 
y'[t] == 
            x[t] + 0.2 y[t], 
          y[0] == -0.3, z'[t] == 
            0.2 + x [t] z[t] - 5.7 z[t], z[0] == 0.2}, {x, y, z}, {t, 
          0, tf}, Method -> RungeKutta, 
        MaxSteps -> 2000]

    
plot1 = 
      ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. sol], {t, 0, tf}, 
        PlotPoints -> 6000, 
        PlotRange -> All]

Clear[x, y, z, f1, f2, f3, t, c, 
      sol]

tf  =  210;
sol  = 
      NDSolve[{
x'[t] == -(y[t] + z[t]), 
          x[0] == -0.04, 
y'[t] == 
            x[t] + 0.2 y[t], 
          y[0] == -0.3, z'[t] == 
            0.2 + x [t] z[t] - 5.7 z[t], z[0] == 0.2}, {x, y, z}, {t, 
          0, tf}, Method -> RungeKutta, 
        MaxSteps -> 2000]

plot1 = 
      ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. sol], {t, 0, tf}, 
        PlotPoints -> 6000, 
        PlotRange -> All]

An nimation of the Rossler attractor:

Do[{Clear[x, y, z, f1, f2, f3, t, c, 
          sol], 
c  =  5.7, 
sol  = 
          NDSolve[{
x'[
                  t] == -(y[t] + z[t]), 
              x[0] == -0.04, 
y'[t] == 
                x[t] + 0.2 y[t], 
              y[0] == -0.3, z'[t] == 
                0.2 + x [t] z[t] - c z[t], z[0] == 0.2}, {x, y, 
              z}, {t, 0, tf}, Method -> RungeKutta, 
            MaxSteps -> 2000], 
plot1 = 
          ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. sol], {t, 0, tf}, 
            PlotPoints -> 4000, 
            PlotRange -> {{-10, 15}, {-12, 10}, {0, 25}}]}, {tf, 
        10, 110, 1}]

Clear[x, y, z, f1, f2, f3, t, c, sol]

    
tf  =  200;
sol  = 
        NDSolve[{
x'[
                t] == -(y[t] + z[t]), 
            x[0] == -0.04, 
y'[t] == 
              x[t] + 0.2 y[t], 
            y[0] == -0.3, z'[t] == 
              0.2 + x [t] z[t] - 5.7 z[t], z[0] == 0.2}, {x, y, 
            z}, {t, 0, tf}, Method -> RungeKutta, 
          MaxSteps -> 2000];

    
xone[t_]  =  {x[t]} /. sol;

    
plot1 = 
      Plot[Evaluate[xone[t]], {t, 0, tf}, PlotPoints -> 6000, 
        PlotRange -> All, 
        PlotStyle -> {RGBColor[0, 0, 1]}]

How we know we have chaos in the Rossler attractor:

Clear[x, y, z, f1, f2, f3, t, c, sol]

    
tf  =  200;
sol  = 
        NDSolve[{
x'[
                t] == -(y[t] + z[t]), 
            x[0] == -0.045, 
y'[t] == 
              x[t] + 0.2 y[t], 
            y[0] == -0.3, z'[t] == 
              0.2 + x [t] z[t] - 5.7 z[t], z[0] == 0.2}, {x, y, 
            z}, {t, 0, tf}, Method -> RungeKutta, 
          MaxSteps -> 2000];

    
plot2 = 
      Plot[Evaluate[{x[t]} /. sol], {t, 0, tf}, PlotPoints -> 6000, 
        PlotRange -> All, 
        PlotStyle -> {RGBColor[0, 0, 1]}]

    
plot3 = 
      Plot[Evaluate[xone[t]  - {x[t]} /. sol], {t, 0, tf}, 
        PlotPoints -> 6000, PlotRange -> All, 
        PlotStyle -> {RGBColor[1, 0, 0]}]

The Henon Map

Clear[x, y]

x[0]  :=  -1
y[0]  :=  0

x[n_]  :=  x[n]  =  
        y[n - 1]  + 1 - 1.4  x[n - 1]^2
y[n_]  :=  y[n] = 
        0.3  x[n - 1]

points  =  
        Table[{x[n], y[n]}, {n, 0, 
            12000}];

plot1  =  
      ListPlot[points, Frame -> True, Axes -> False, 
        FrameLabel -> {"x", "y"}]

Clear[x, y]

x[0]  :=  1
y[0]  :=  -0.3

x[n_]  :=  x[n]  =  
        y[n - 1]  + 1 - 1.4  x[n - 1]^2
y[n_]  :=  y[n] = 
        0.3  x[n - 1]

points  =  
        Table[{x[n], y[n]}, {n, 0, 
            12000}];

plot2  =  
      ListPlot[points, Frame -> True, Axes -> False, 
        FrameLabel -> {"x", "y"}]

Show[plot1, plot2]

Fractals

Here is the Julia set, on two different scales, to show self similarity (fractal):

Clear[z]
c  =  I
z[0]  :=  0.5
z[n_] := z[n] = 
        If[Random[Integer] == 1, 
          Sqrt[z[n - 1] - c], -Sqrt[
              z[n - 1] - 
                c]]

JuliaSet  =  
        Table[{Re[z[n]], Im[z[n]]}, {n, 0, 
            500000}];
ListPlot[Take[JuliaSet, 10000], AspectRatio -> 1,  
      Frame -> True,  Axes -> False, 
      PlotRange -> {{-1.5, 1.5}, {-1.5, 
            1.5}}]
ListPlot[JuliaSet, AspectRatio -> 1,  Frame -> True,  
      Axes -> False, 
      PlotRange -> {{-0.15, 0.15}, {-0.15, 0.15}}]

Clear[z]
c  = I
z[0]  :=  -1.0 - 1.0  I
z[n_] := z[n] = 
        If[Random[Integer] == 1, 
          Sqrt[z[n - 1] - c], -Sqrt[
              z[n - 1] - 
                c]]

JuliaSet  =  
        Table[{Re[z[n]], Im[z[n]]}, {n, 0, 
            500000}];
ListPlot[Take[JuliaSet, 10000], AspectRatio -> 1,  
      Frame -> True,  Axes -> False, 
      PlotRange -> {{-1.5, 1.5}, {-1.5, 
            1.5}}]
ListPlot[JuliaSet, AspectRatio -> 1,  Frame -> True,  
      Axes -> False, 
      PlotRange -> {{-0.15, 0.15}, {-0.15, 0.15}}]

References

Fractals:

http://mathforum.org/advanced/robertd/index.html#frac

Attractors, with animation:

http://www.cs.gsu.edu/~cscrwh/dynam/dynam.html

Fractal music:

http://www.thinks.com/sounds/fractal.htm