Mathematica Activity 1:
Visualising dynamical systems

by Phillip Kent, METRIC Project, 
for M3A16/M4A16, 1999-2000.
Version 1.1, 26 October 1999.

These activities are based on a "toolkit" of Mathematica functions, designed to help you to explore various (mostly graphical) aspects of dynamical systems. In this first activity we focus on two standard examples, the simple pedulum and the Kepler problem.

The Mathematica notebook containing this document, plus related materials can be accessed at

http://metric.ma.ic.ac.uk/dynamics/

Note that you can save yourself a LOT of typing by downloading this notebook!

This activity forms part of your assessed coursework. You will be required to submit a "portfolio" of graphics and explanatory text, preferably in the form of a printed Mathematica notebook. The different pieces of "portfolio work" that we would like you to produce are clearly indicated by the text cells with grey background. Please show all the Mathematica code that you use and explain fully the steps in your calculations.

This coursework should be submitted by THURSDAY 18th NOVEMBER.

0. Loading the dynamics toolkit

On the Mathematics Department Windows NT network, you can load the toolkit into Mathematica with the following commands:

AppendTo[$Path,  "//OPHELIA/OPHELIA_D"]

<< dynapack.m

To install the toolkit file locally in NT, or on your own computer, see the instructions on the web page.

Please send any questions, comments or bug reports to p.kent@ic.ac.uk . Note that some of the toolkit commands can easily cause Mathematica to crash, so save your notebooks regularly.

1. The simple pendulum - Trajectories in the phase plane

The basic Mathematica object which this toolkit deals with is a “Trajectory”. This is not part of regular Mathematica, and the toolkit contains various commands to create, plot and manipulate this new kind of object. There is a reference summary for the major commands in the toolkit in Section 3.

Making a Trajectory

The dynamical system for the simple pendulum (with mass set to 1 and gravitational constant set to 1) is described by the equations:

dθ/dt = \(TraditionalForm\`θ\& . \),  d\(TraditionalForm\`θ\& . \)/dt = - sin(θ).

A Trajectory is created with the MakeTrajectory command, for example:

traj1 = 
      MakeTrajectory[{θ'[t] == θdot[t], θdot'[
              t] == -Sin[θ[t]]}, {θ[0] == 
            19  π/20, θdot[0] == 
            0.0}, {θ, θdot}, {t, 12}]

traj1 now represents a numerical solution of the simple pendulum for the given initial conditions (θ = 19π/20, "θdot" = 0) and time interval (12 units).

Note that MakeTrajectory makes use of Mathematica's built-in numerical solver, NDSolve.

Plotting Trajectories

Several commands exist for making plots of Trajectories. For example, to see the variation of θ against t (sometimes called the "time series" for θ):

TimeSeries[traj1, θ]

And to see the "phase plane":

TrajPlot[traj1, {θ, θdot}, AspectRatio -> Automatic]

There is an option you can set so that the direction of time evolution is indicated by a "blob" at the start, and an arrow at the end of the Trajectory:

TrajPlot[traj1, {θ, θdot}, AspectRatio -> Automatic, 
      DirectionQ -> True]

You can plot any number of Trajectories on the same axes by specifying a list of them as the first argument of TrajPlot. For example:

traj2 = 
      MakeTrajectory[{θ'[t] == θdot[t], θdot'[
              t] == -Sin[θ[t]]}, {θ[0] == 
            2  π/3, θdot[0] == 0.0}, {θ, θdot}, {t, 
          10}]

TrajPlot[{traj1, traj2}, {θ, θdot}, 
      AspectRatio -> Automatic, DirectionQ -> True]

You can combine Trajectory plots with regular Mathematica plots. Here, for example, you might want to compare the numerically-calculated Trajectory "traj2" with the corresponding level curve for the pendulum system:

phaseplot = 
      TrajPlot[traj2, {θ, θdot}, AspectRatio -> Automatic, 
        DirectionQ -> True, PlotStyle -> RGBColor[1, 0, 0]]

levelcurve = 
      Plot[{(Sqrt[2 ((1/2) + Cos[θ])]), 
          -(Sqrt[2 ((1/2) + Cos[θ])])}, {θ, 
          -2 π/3, 2  π/3}, AspectRatio -> Automatic]

Show[levelcurve, phaseplot, 
      AxesLabel -> {"θ", "θdot"}]

The agreement is very good. But you will see a significant "drift" if you integrate the system over longer times, say 200 units or more (this requires a larger value for the option MaxSteps in MakeTrajectory, say MaxSteps->2000). (Note: this effect can be reduced by making NDSolve perform more accurately, but this only works up to the "machine accuracy", about 16 digits, of the computer.)

Energy conservation

Conserved quantities are very important in Hamiltonian dynamics, and you can check for conservation along Trajectories using the CheckQuantity command. For example (note: the quantity here should not involve use of the t variable, so θ and not θ[t], etc):

energyvals1 = 
      CheckQuantity[traj1, θ((dot^2)/2) - Cos[θ]]

ListPlot[energyvals1]

So, you should see that energy is conserved to about 5 decimal places, which is generally OK for "qualitative" work such as we're doing here.

Portfolio Work 1.1

Produce a phase plane plot for the simple pendulum containing a selection of  bounded (oscillating) and unbounded ("loop-the-loop") orbits, and, in a contrasting plot style, the separatrix* curves.

* Separatrix = an exceptional orbit that tends to the unstable equilibrium point at θ =π for both  t ->-∞ and t -> +∞; in the phase plane, it separates regions of bounded and unbounded orbits.

Period as a function of amplitude

An important difference between the simple pendulum and the simple harmonic oscillator is that for the pendulum the period of an oscillating orbit depends on the amplitude (as you have probably already noted from the previous task).

Let's create a set of initial points with increasing amplitude:

startpoints = Table[{a π, 0}, {a, 0.1, 0.9, 0.1}]

Then, we'll "shoot" Trajectories, for the same time interval, from each of these initial points. First, for the SHO:

shootingsSHO = 
      Map[MakeTrajectory[{x'[t] == xdot[t], xdot'[t] == -x[t]}, {
              x[0] == #[[1]], xdot[0] == #[[2]]}, {x, xdot}, {t, 
              5}]&, startpoints]

And for the pendulum:

shootingsPend = 
      Map[MakeTrajectory[{θ'[t] == θdot[t], 
              θdot'[t] == -Sin[θ[t]]}, {
              θ[0] == #[[1]], θdot[0] == #[[2]]}, {
              θ, θdot}, {t, 5}]&, startpoints]

Compare:

TrajPlot[shootingsSHO, {x, xdot}, DirectionQ -> True, 
      AspectRatio -> Automatic, PlotStyle -> {Hue[0], Hue[0.15], Hue[0.4]}]

TrajPlot[shootingsPend, {θ, θdot}, DirectionQ -> True, 
      AspectRatio -> Automatic, PlotStyle -> {Hue[0], Hue[0.15], Hue[0.4]}]

You'll see more about the relation between period and amplitude later in the course (under "action-angle variables").

Plotting in cylindrical coordinates [OPTIONAL]

An alternative visualisation of the pendulum Trajectories is given by plotting them in cylindrical coordinates. For this, it is useful to construct a "wire frame" cylinder:

Needs["Graphics`Shapes`"]

cylinderWire = 
      Show[WireFrame[
          First[ParametricPlot3D[{Cos[θ], Sin[θ], z, 
                GrayLevel[0.5]}, {θ, 0, 2  π}, {z, -3, 3}]]], 
        Axes -> True]

(The effect of GrayLevel[0.5] is to grey out the wires so that the Trajectories will show more clearly.) It's also useful to construct a θ = 0 line and a θdot = 0 circle:

thetazero = 
      ParametricPlot3D[{1, 0, 
          z, {Thickness[0.012], RGBColor[0, 0, 1]}}, {z, -3, 
          3}]

thetadotzero = 
      ParametricPlot3D[{Cos[t], Sin[t], 
          0, {RGBColor[0, 0, 1], Thickness[0.012]}}, {t, 0, 
          2 π}]

To draw a Trajectory in the cylindrical system, the θ and θdot components are first extracted using the TrajComponent command:

traj3 = 
      MakeTrajectory[{θ'[t] == θdot[t], θdot'[
              t] == -Sin[θ[t]]}, {θ[0] == 
            19  π/20, θdot[0] == 
            0.0}, {θ, θdot}, {t, 16}]

cyltraj1 = 
      ParametricPlot3D[{Cos[TrajComponent[traj3, θ][t]], 
          Sin[TrajComponent[traj3, θ][t]], TrajComponent[
              traj3, θdot][t], Thickness[0.012]}, {t, 0, 
          16}]

And here's the full cylindrical phase plot:

cylphaseplot = 
      Show[cylinderWire, thetazero, thetadotzero, cyltraj1]

You should add further bounded and unbounded Trajectories to this plot.

To view this plot from different directions, you have a couple of choices. One is to use the ViewPoint option to specify the viewing point by {x,y,z} coordinates. For example

Show[cylphaseplot,  ViewPoint -> {1, 1, 1}]

The other approach is to use a new, experimental feature in Mathematica version 4, for "real time" interactive 3-D graphics. We have to say that its performance is limited, and a big problem with it is that it does not handle plot styles such as colours and line thicknesses (which is rather crucial for these cylindrical plots). Anyway, to try it out enter:

<< Experimental`

<< RealTime3D`

The back quotes ` here are essential. Now evaluate:

Show[thetazero, thetadotzero, cyltraj1]

The output you will get is not a regular graphic but an interactive one: click and drag the mouse to rotate it. You can zoom in or out by pressing (in Windows) the "Shift" and "Ctrl" keys together and then dragging the mouse up or down. All 3-D graphics will come out in "real time"  versions until you switch back to regular 3D graphical output by entering:

<< Default3D`

Again, note the back quote `.

2. The Kepler problem

From your course notes -- or from Problem Sheet 2, Question 5 -- you should have that the Hamiltonian for the Kepler problem is:

\(\(keplerHam = 
        pr[t]\^2\/\(2  m\) + 
          pθ[t]\^2\/\(2  m\ r[t]\^2\) - \(μ\ m\)\/r[t];\)\)

We'll use "pr" and "pθ" instead of the conventional subscripted forms, (Subscript[p, r]) and (Subscript[p, θ]), since Mathematica cannot handle these very well in differential equations.

Then, the equations of motion can be calculated:

keplerEqns = {r'[t] == D[keplerHam, pr[t]], θ'[t] == 
          D[keplerHam, pθ[t]], pr'[
            t] == -D[keplerHam, r[t]], pθ'[t] == 0}

Here is a sample elliptic orbit. Note the use of "/." to set values for m and μ:

keptraj1 = 
      MakeTrajectory[
        keplerEqns /. {m -> 1, μ -> 1}, {r[0] == 5, θ[0] == 0, 
          pr[0] == 0.1, pθ[0] == 1}, {r, θ, pr, pθ}, {t, 
          25}]

pθ should be a constant of the motion:

TimeSeries[keptraj1, pθ]

Here is the projection of the (4-dimensional) orbit in the (r, pr) plane:

TrajPlot[keptraj1, {r, pr}, AspectRatio -> Automatic, 
      DirectionQ -> True]

To see the orbit in Cartesian coordinates, the Transform option can be used. First, we define a transform function from (r, θ) to (x, y):

keplertrans[{r_, θ_}] := {r Cos[θ], 
        r Sin[θ]}

And then apply that transformation to an (r, θ) plot:

TrajPlot[keptraj1, {r, θ}, Transform -> keplertrans, 
      AspectRatio -> Automatic, DirectionQ -> True]

A consequence of the conservation of angular momentum is Kepler's Second Law, which states that the areal velocity of a particle moving in an elliptical orbit under gravity is a constant -- or, that equal areas are swept out in equal times. As an illustration of this, the following code generates an animation in which two lines join the origin (one of the foci of the ellipse) to the position of the particle at time t, and at time t - 1.

keptraj2 = 
      MakeTrajectory[
        keplerEqns /. {m -> 1, μ -> 1}, {r[0] == 
            1.0, θ[0] == 1.0, pr[0] == 0.9, 
          pθ[0] == 1}, {r, θ, pr, pθ}, {t, 
          76}]

orbit = 
      TrajPlot[keptraj2, {r, θ}, Transform -> keplertrans, 
        AspectRatio -> Automatic]

orbitPoint[traj_Trajectory, 
        t_] := 
TrajComponent[traj, r][
          t]*{Cos[TrajComponent[traj, θ][t]], 
          Sin[TrajComponent[traj, θ][t]]}

Table[
      Show[orbit, 
Graphics[{RGBColor[1, 0, 0], 
            Line[{{0, 0}, orbitPoint[keptraj2, t]}], 
            Line[{{0, 0}, orbitPoint[keptraj2, t - 1]}]}], 
        AspectRatio -> Automatic, 
        PlotRange -> {{-8.75, 0.75}, {-1, 6}}], {t, 1, 
        76}]

To animate these frames, the first thing you have to do is double-click on the cell bracket that contains them all (on the right-hand side of the Mathematica window). They should all disappear apart from the first frame in the sequence. Then, select Animate Selected Graphics from the Cell menu.

You can change the speed or direction of the movie using the buttons at the bottom left of the screen. To stop, click the mouse button anywhere in the notebook. If you want to resize the frames in the animation, resize the first frame and then reevaluate the code for the animation.

Portfolio Work 1.2

(cf. Problem Sheet 2, Question 5)

Use ContourPlot to plot the level curves H = constant in the (r, pr) plane. For simplicity, take m = 1, μ = 1 and pθ = 1. (To get smooth looking curves, you will need to increase the value of PlotPoints, say PlotPoints->50.)

Portfolio Work 1.3

Make a set of Trajectories for some elliptic, parabolic and hyperbolic motions. Again, for simplicity take m = 1, μ = 1 and pθ = 1. 

Make plots of these Trajectories in the (r, pr) projection and in (x, y) space. 

For the unbounded orbits, can you think of an easy way to generate both the t -> +∞ and t ->-∞ sections of the orbit?

3. Toolkit reference section

To find out about the toolkit commands whilst using Mathematica, recall the use of ? to bring up a "usage message":

?TrajPlot3D

MakeTrajectory

Syntax

Examples

lorenz[σ_, r_, b_] := {x'[t] == σ (y[t] - x[t]), 
        y'[t] == (r - z[t]) x[t] - y[t], 
        z'[t] == x[t] y[t] - b z[t]}

lorenz1 = 
      MakeTrajectory[
        lorenz[10, 28, 8/3], {x[0] == 10.0, y[0] == 10.0, z[0] == 10.0}, {x, 
          y, z}, {t, 10}]

lorenz2 = 
      MakeTrajectory[
        lorenz[10, 28, 8/3], {x[0] == 12.0, y[0] == 10.0, z[0] == 10.0}, {x, 
          y, z}, {t, 10}, MaxSteps -> 1500]

Options

MakeTrajectory has no options of its own, but you may need to pass options to NDSolve, which is called internally. Some Trajectories will exceed the default value of MaxSteps.

Note: Trajectory objects have three pieces: a list of the InterpolatingFunction objects returned by NDSolve; a rule, indVar->{...}, which contains the list of independent variable values at which the interpolating functions are determined; and a rule, varLabels->{...}, which is the list of variable names (the first being the independent variable).

TimeSeries

Plots the variation of one component of a Trajectory against the independent variable ("time").

TimeSeries[lorenz1, x]

TimeSeries[{lorenz1, lorenz2}, x, 
      PlotStyle -> {Automatic, RGBColor[1, 0, 0]}]

You can specify all the usual options for Plot with TimeSeries - for example, PlotStyle, Axes, AspectRatio .

TimeSeries3D

Plots the variation of two components of a Trajectory against the independent variable ("time").

TimeSeries3D[lorenz1, {x, z}, BoxRatios -> {3, 1, 1}]

You can specify all the usual options for ParametricPlot3D with TimeSeries - for example, PlotStyle, Axes, BoxRatios .

TrajPlot (or TrajPlot2D)

Plot a graph of two components of a Trajectory against each other - i.e. the phase plane for a 2-D system, or a phase space projection for 3-D and higher systems.

TrajPlot[lorenz1, {x, z}, DirectionQ -> True]

TrajPlot[{lorenz1, lorenz2}, {x, z}, DirectionQ -> True, 
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 0, 1]}]

You can specify all the usual options for Plot with TrajPlot - for example, PlotStyle, AxesLabel, AspectRatio .

With DirectionQ->True (default False), the direction of "time" evolution is indicated by a "blob" at the initial point and an arrow at the end of the Trajectory. A plot containing a lot of Trajectories may look nicer if you remove the initial blobs: the option to do this is StartObject->({}&).

Transform->transformfunc allows the (variable1, variable2) coordinate system to be transformed. transformfunc must take a list of two numbers as input, and return a list of two numbers. See the example in Section 2.

TrajPlot3D

Plot a graph of three components of a Trajectory against each other.

TrajPlot3D[lorenz1, {x, y, z}]

TrajPlot3D[{lorenz1, lorenz2}, {x, z, y}, DirectionQ -> True, 
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 0, 1]}]

You can specify all the usual options for ParametricPlot3D with TrajPlot3D - for example, PlotStyle, AxesLabel, BoxRatios .

With DirectionQ->True (default False), the direction of "time" evolution is indicated by a "blob" at the initial point and an arrow at the end of the Trajectory. 

Transform->transformfunc allows the (variable1, variable2,variable3) coordinate system to be transformed. transformfunc must take a list of three numbers as input, and return a list of three numbers.

TrajComponent

Selects the component of a Trajectory given by the specified variable.

TrajComponent[lorenz1, x]

The output of TrajComponent is an "InterpolatingFunction" object, which is handled like a function in one variable, e.g.

TrajComponent[lorenz1, x][0]

Plot[TrajComponent[lorenz1, x][t], {t, 3, 10}]

xcomp[t_] = TrajComponent[lorenz1, x][t]

xcomp[5.34]

Plot[xcomp[t], {t, 0, 10}]

CheckQuantity

Generates a list of numerical values of the given expression along a Trajectory. The expression must be given in terms of the dependent variables belonging to the Trajectory.

Useful for checking on whether quantities are conserved or not - see the energy example in Section 1.

CheckQuantity[lorenz1, (Sqrt[(x^2) + (y^2) + (z^2)])]

ListPlot[CheckQuantity[lorenz1, (Sqrt[(x^2) + (y^2) + (z^2)])]]

CutTrajectory

Outputs a new Trajectory containing the section of the input Trajectory lying in the range {t1, t2} of the independent variable.

lorenz1cut = CutTrajectory[lorenz1,  {2,  5}]

TrajPlot[CutTrajectory[lorenz1,  {5, 8}], {x, z}, DirectionQ -> True]

You can use CutTrajectory to make animations of the time evolution of a Trajectory:

Table[
      TrajPlot3D[CutTrajectory[lorenz1, {0, T}], {x, y, z}, Axes -> True, 
        DirectionQ -> True, AxesLabel -> {"x", "y", "z"}, 
        PlotRange -> {{-20, 20}, {-25, 25}, {5, 50}}], {T, 0.5, 10, 
        0.5}]

TrajRange is useful for determining the overall plot ranges here.

TrajRange

Calculates the minimum and maximum values for the given variable along the input Trajectory.

TrajRange[lorenz1, x]

TrajRange[lorenz1, y]

TrajRange[lorenz1, z]