. Dynamical Systems and Chaos

Paul Abbott

Pre-requisites

To improve the appearance of graphics, TraditionalForm will be used as the default format type for text that appears in graphics,

$FormatType = TraditionalForm;

and the default style for text in graphics is set as follows.

$TextStyle = {FontSize -> 10};

The following options for ParametricPlot, DensityPlot, and ListPlot are desirable.

SetOptions[ParametricPlot, 
        AspectRatio -> Automatic];

SetOptions[DensityPlot, PlotPoints -> 100, 
        Mesh -> False, Compiled -> False, 
        ColorFunction -> Hue];

SetOptions[ListPlot, 
        PlotStyle -> AbsolutePointSize(0.002)];

Consider a particle of mass m that is constrained by a rigid and massless rod to move in a vertical circle of radius l.  The equation of motion is

(Overscript[θ, ¨])(
          t) + (ω^2) sin(θ(t)) = 0,  ω = 
      (g/l),

where θ is the angle measured from the vertical direction with clockwise positive.

Usually, only the harmonic approximation (i.e., sin(θ) is replaced by θ) of this equation is considered:

(Overscript[θ, ¨])(
          t) + (ω^2) θ(t) = 0.

Here, the more general equation (which includes damping when β  ≠  0 and a sinusoidal driving term when α  ≠  0)

\(TraditionalForm\`\(θ\&¨\)(
          t) + β\ \(\(θ\& . \)(
            t)\) + ω\^2\ \(sin(θ(
              t))\) = α\ \(sin(ω\_d\ t)\)\)

will be examined. The units of mass, length, and time are chosen so that they are all equal to unity, i.e., we are working with dimensionless variables.

.. Undamped pendulum

The dimensionless equation for the undamped pendulum in Mathematica notation reads

𝒫=θ''(t)+sin(θ(t))==0;

.. Damped pendulum

The equation for the damped pendulum (where the damping is proportional to the velocity) can be written:

𝒫=θ''(t)+β θ'(t)+sin(θ(t))==0;

.. Damped driven pendulum

The equation for the damped driven pendulum can be written

𝒫=θ''(t)+β θ'(t)+sin(θ(t))==α cos(ω t);

You should see that the resulting phase-space diagram is periodic (in the language of chaos it shows a period-3 limit cycle).

The Poincaré section provides a simplification of the phase diagram by sampling at time intervals equal to the period of the driving frequency.

To avoid transient effects, the data for the first few cycles should be discarded.

.. Series solution

The easiest way to get Mathematica to solve differential equations such as y''(x)+(k^2) y(x)=0 using series methods is direct expansion of the equation, i.e.,

eqn=y''(x)+(k^2) y(x)+((O(x))^7)==0

This is equivalent to expanding y(x) (and y''(x)) into a power series.

yseries = y(x) + ((O(x))^7)

You can see that the expression multiplying each power of x in the expansion of y''(x)+(k^2) y(x) must vanish identically, i.e., y(0) (k^2)+y''(0)=0, y'(0) (k^2)+Derivative[3][y](0)=0, …   To see the equations generated you can use LogicalExpand (this step is not required).

LogicalExpand[eqn]

Note that the initial conditions, y(0) and y'(0) are not determined but that you can solve for the undetermined coefficients (y''(0), Derivative[3][y](0), …) in terms of the initial conditions.

Solve[eqn,Table[Derivative[n][y](0),{n, 2, 8}]]

To obtain the series solution to the original equation you substitute these coefficients into the series expansion of y(x).

yseries /. First[%]

If you now supply the initial conditions — say y(0) = 1 and y'(0)=0 — you get

%/.{y(0) -> 1,y'(0)->0}

.. Series solution for the undamped pendulum

We would like to analytically solve the equation

θ''(t)+(ω^2) sin(θ(t))=0.

The linearized equation results by replacing sin(θ) by θ.

θ''(t)+(ω^2) θ(t)=0.

. Iterative Operations

This section presents four examples of visualizing iterative operations.

.. Julia set

The Julia set associated to a fixed complex number c is found by iterating the map z ↦ (z^2) + c. The set of points, z, that do not escape to infinity (i.e., their absolute value remains bounded, say |z| < 2, after m iterations have been carried out) is the Julia set.

The following Mathematica procedural code (i.e., similar in style to Fortran, C, Pascal, Basic,…) iterates z = (z^2) + c for a particular value of c.  Either the value of z escapes to ∞ which can be tested using |z| < 
        2.0 or it hops around a finite number of times, bounded above by m.

JuliaModule(w_, c_, m_) := 
        Module({z = w, n = 0}, 
          While(\[LeftBracketingBar]z\[RightBracketingBar] < 2 ∧ n < m, 
            z = (z^2) + c; ++n); n);

Alternatively, we can implement this as a pure function (using Function) with {w, c, m} as dummy arguments.

JuliaFunction = 
        Function({w, c, m}, 
          Module({z = w, n = 0}, 
            While(\[LeftBracketingBar]z\[RightBracketingBar] < 2 ∧ 
                n < m, z = (z^2) + c; ++n); n));

A table of values using either implementation is somewhat slow:

Table[
            JuliaModule(x + I y, 1 + 0.3 I, 
              35), {x, -1.5, 1.5, 0.1}, {y, -1.5, 1.5, 0.1}]; // 
        Timing // First

Table[
            JuliaFunction(x + I y, 1 + 0.3 I, 
              35), {x, -1.5, 1.5, 0.1}, {y, -1.5, 1.5, 0.1}]; // 
        Timing // First

Numerical function compilation, using Compile, is a useful way of speeding up a computation. The following code expects three arguments: the point z, the point c, and the maximum number of iterations.

JuliaCompile=Compile[(({{w, _Complex}, {c, _Complex}, {m, _Integer}})),Module[{z = w, n = 0}, 
                While[\[LeftBracketingBar]z\[RightBracketingBar] < 2 ∧ 
                    n < m, z = (z^2) + c; n++]; n]];

(({{w, _Complex}, {c, _Complex}, {m, _Integer}}))

or as a list of lists:

{{w, _Complex}, {c, _Complex}, {m, _Integer}}

Algorithms for computing fractal sets often involve iterating a formula until a specified condition is met. By default, FixedPoint and FixedPointList stop iterating when two successive calculations give the same result. To specify our own stopping condition, we can use the SameTest option. This option takes a pure function of two arguments that returns True if the two arguments can be considered “the same”. For example, here we iterate the square root function, applied to 10, and stop when two successive values are within 0.1 of each other.

FixedPoint(Sqrt, 10.0, 
      SameTest -> (\[LeftBracketingBar]#1 - #2\[RightBracketingBar] < 
              0.1 &))

Here is an elegant alternative implementation of JuliaCompile using a fixed number of iterations (here 35).

JuliaCompile=Compile((({{z, _Complex}, {c, _Complex}})),\(Length(
                FixedPointList(#\^2 + c &, z, 35, 
                  SameTest -> \((\[LeftBracketingBar]#2\[RightBracketingBar] > 2.0 &)\)))\));

Using \(SameTest -> \((|#2| < 2.0 &)\)\) tests whether the last value (i.e., #2) generated by FixedPointList remains bounded (i.e., |z|<2). The length of the result of FixedPointList is the number of iterations.

For fun, run the following code to produce an array of plots with z  ∈  [-1.1 - 1.1  i, 
          1.1 + 1.1  i] and \(TraditionalForm\`c = 
        a + b\ i\ \( ∈ \_\ \)\([\(-2\) - i, 
            1 + i]\)\).

Show(
        GraphicsArray(
          Table(DensityPlot(-JuliaCompile(x + I y, 
                  a + b I), {x, -1.1, 
                1.1}, {y, -1.1, 1.1}, Frame -> False, 
              Axes -> False, PlotPoints -> 20, 
              DisplayFunction -> Identity), {b, -1, 1, 
              (1/3)}, {a, -2, 1, (1/3)})));

.. The Mandelbrot set

For each c in the complex plane, those points that remain bounded upon iterating the map c ↦ (c^2) + c generate the Mandelbrot set. It turns out that points in the Mandelbrot set have connected Julia sets.

.. Newton's Method

Consider using Newton's method to compute the (complex) argument of the n-th roots of unity, i.e., the roots of f(z) = (z^n) - 1. The Newton iteration, z->z-((f(z))/(f'(z))), becomes

z-((f(z))/(f'(z)))/.f -> Function(z, (z^n) - 1)//Simplify

The compiled function

NewtonIteration=Compile((({{n, _Integer}, {c, _Complex}})),Arg(
                FixedPoint(Function(z, ((z ((z^(-n)) + n - 1))/n)), c, 
                  35)));

can be used to show the convergence of points in the complex plane for specific n.

For fun, run the following code to produce an array of graphics (generated using GraphicsArray) that show the convergence of points in the complex plane, for n ranging from 5 to 13.

Show(
        GraphicsArray(
          Partition(
            Table(DensityPlot(
                NewtonIteration(n, x + I y), {x, 
                  1.1, -1.1}, {y, -1.1, 1.1}, PlotPoints -> 30, 
                Compiled -> False, Frame -> False, 
                DisplayFunction -> Identity), {n, 5, 13}), 3), 
          GraphicsSpacing -> 0.05));

.. Feigenbaum Bifurcation Diagram

The Feigenbaum bifurcation diagram shows the behavior of iterates of x ↦ μ x (1 - x) and displays the period-doubling route to chaos. The following code iterates this mapping 300 times using NestList and then deletes the first 100 iterations (during which time the iteration "settles down").

Feigenbaum=Compile((({{μ, _Real}})),({μ, #} &) /@ Drop(
                  NestList(μ # (1 - #) &, 0.2, 300), 100));