Objects.nb

This notebook accompanies the article 
Using a Computer Algebra System to Create a Language for Mathematical "Objects", published in the TALUM Newsletter number 9 (September 1998)

by Phillip Kent (p.kent@ic.ac.uk)

© Copyright Imperial College of Science, Technology and Medicine 1999

http://www.bham.ac.uk/ctimath/talum/newsletter/talum9.htm

http://metric.ma.ic.ac.uk/articles/talum/

This notebook contains graphs and animations to illustrate the article. If you are viewing this with the MathReader application, you may only view the contents and play animations. If you have a full version of Mathematica (version 3 or 4), consult the "Code setup" section for installing and loading the necessary Mathematica programs.

Code setup (for full Mathematica only)

Downloading the extra Mathematica packages

Before you can use the code in this notebook, you will need two extra code packages. They may be placed, for example, either in the main Mathematica folder, typically inside "AddOns/Applications", or, put them in the same folder as this notebook itself.

ExtendGraphics

These packages accompany the book Mathematica Graphics by Tom Wickham-Jones. Note that I only use a few of the packages from the full set, i.e. "Geometry" and "Geometry3D". Download the zip file either from:

http://metric.ma.ic.ac.uk/articles/talum/ExtendGraphics30.zip

or from:

http://www.mathsource.com/Content/Enhancements/Graphics/3D/0208-976

You must unzip this file into a new folder called "ExtendGraphics" - most modern unzip programs will do this automatically.

The ExtendGraphics programs are copyrighted by Springer-Verlag - see the README.TXT file.

dynapack.m

This contains functions written by me for working with "Trajectory" objects, as demonstrated in Sections 3 and 4 of this notebook.

http://metric.ma.ic.ac.uk/articles/talum/dynapack.m

EVALUATE ALL THE CODE IN THIS SECTION!

(Suppress annoying spell warnings)

Off[General::spell1]

My package for dynamics:

<< dynapack.m

The Wickham-Jones packages:

<< ExtendGraphics`Geometry`

And an extra definition for convenience:

IntersectionPoint[l1_Line, l2_Line] := 
      IntersectionPoint[l1, ToImplicit[l2]]

<< ExtendGraphics`Geometry3D`

The convenience Draw function (MyFilterOptions is defined in dynapack.m):

Draw[object_, opts___] := 
      Show[Graphics[object], 
        MyFilterOptions[Graphics, Join[{opts}, Options[Draw]]]]

\(\(\( (*not\ sure\ if\ this\ is\ required, \ 
      but\ does\ no\ harm ... *) \)\(
\)\(Draw[objects_List, opts___] := 
      Show[Graphics[objects], 
        MyFilterOptions[Graphics, 
          Join[{opts}, Options[Draw]]]]\)\)\)

Options[Draw] = {AspectRatio -> Automatic};

The functions for Poincaré sections:

TestForIntersection[{p1_List, p2_List}, pl_Plane]:=
 If[
  Not[SideOfPlaneQ[p1, pl]] && 
    (SideOfPlaneQ[p2, pl] || PointInPlaneQ[p2, pl]),
   IntersectionPoint[Line[{p1, p2}], pl],
   {}
 ]

PoincareSection[traj_Line, pl_Plane]:=
 Map[Point,
 DeleteCases[
   Map[TestForIntersection[#, pl]&, 
    Transpose[{Drop[traj[[1]], -1], Rest[traj[[1]]]}]],
   {}]
 ]

\(PoincarePoints[traj_Trajectory, pl_Plane] := 
      Module[{trajline}, \((
			trajline = \ 
            Line[Transpose[
                TrajToPoints[traj]]]; \ 
\ \ \ \ \ \ \ PoincareSection[
            trajline, \ pl]
		)\)]\)

(Turn spell warnings back on)

On[General::spell1]

1. Circle example

figure1 = Draw[Circle[{1, 3}, 1.5], Axes -> True]

(This is the usual "raw" Mathematica version:

Show[Graphics[Circle[{1, 3}, 1.5]], AspectRatio -> Automatic, 
      Axes -> True]

)

2. Wickham-Jones example

P = {0, 0}; Q = {6, 0}; R = {4, 2};

triangle1 = Polygon[{P, Q, R}];

OutCircle[triangle1]

Draw[Boundary[triangle1]]

figure2 = Draw[{Boundary[triangle1], OutCircle[triangle1]}]

An algebraic calculation with a geometrical "object":

triangle2 = Polygon[{{0, 0}, {a, 0}, {a, b}}];

OutCircle[triangle2]

3. Lorenz Trajectory

Define the Lorenz ODEs for convenience:

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]};

Make a Trajectory "object":

traj1 = 
      MakeTrajectory[{x'[t] == 
            10 (y[t] - x[t]), y'[t] == (28 - z[t]) x[t] - 
              y[t], z'[t] == x[t] y[t] - (8/3) z[t]}, {x[0] == 
            10.0, y[0] == 10.0, z[0] == 10.0}, {x, y, z}, {t, 20}, 
        MaxSteps -> 2500]

The "time series" plot:

figure3 = 
      TimeSeries[traj1, x, AxesLabel -> {"t", "x"}]

(The figure in the article shows only a portion of the above plot. One way to achieve this is with the PlotRange option:

TimeSeries[traj1, x, AxesLabel -> {"t", "x"}, 
      PlotRange -> {{0, 11.5}, Automatic}]

)

The "phase space" plot:

figure4 = 
      TrajPlot3D[traj1, {x, y, z}, AxesLabel -> {"x", "y", "z"}, 
        DirectionQ -> True]

The figure in the article restricts the time range to (0,11.5). Another function defined in dynapack.m allows Trajectories to be "cut":

\(TrajPlot3D[
	CutTrajectory[traj1, {0, 11.5}], 
	{x, y, z}, 
      AxesLabel -> {"x", "y", "z"}, 
      DirectionQ -> True]\)

4.  Poincaré sections for the Hénon-Heiles system

NOTE: This is code under development and is given here for demonstration purposes only. This is not intended for use by students, yet. The main fudge is that I do not construct the Poincaré section in the full 4-D phase space. Instead, I manually delete the p1-coordinate from the Trajectory data, and use the intersection plane "q1 = 0" in the 3-D q1-q2-p2 space. I haven't explained how the code works at all --- if you're interested, please get in touch.

(The Hénon-Heiles system was first investigated by Michel Hénon and Carl Heiles in the early 1960s [Astronomical Journal, vol. 69, p. 73, 1964]. As with the Lorenz system, it became a classic of the nonlinear dynamics literature and is discussed in most books on "chaos" [e.g. at a popular level, Chaos by James Gleick, 1987, pages 144-150].)

This is the Hamiltonian (H) for the Hénon-Heiles system:

HHHam = (1/2)*(p1[t]^2 + p2[t]^2) + (1/2)*(q1[t]^2 + q2[t]^2)  +  q1[t]^2*q2[t] - q2[t]^3/3;

From that, I can construct the equations of motion (see any undergraduate textbook on classical mechanics for how and why this works):

hheqns = {q1'[t] == D[HHHam, p1[t]], 
q2'[t] == D[HHHam, p2[t]], p1'[t] == -D[HHHam, q1[t]],
p2'[t] == -D[HHHam, q2[t]]}

The following is a useful function for setting up initial conditions that satisfy a given total energy (i.e. specify E and three of the dependent variables, and solve the H equation for the fourth variable):

\(hhinitial[q1i_, q2i_, p2i_, E_] := 
	{q1[0] == q1i, q2[0] == q2i, 
        p2[0] == p2i, 
		p1[0] == 
          First[Select[

              Map[\(First[#]\)[\([2]\)] &, \ 
                NSolve[\(HHHam\  == \ E\  /. \ 
                      t -> 0\)\  /. \ 
\ \ {q1[0]\  -> \ q1i, \ 
                      q2[0]\  -> \ q2i, 
\ \ p2[0]\  -> \ p2i}, \ 
                  p1[0]]], 
\ \((#\  > \ 0)\)\  &, \ 1]]}\)

This is the plane "q1 = 0" in q1-q2-p2 space:

pplane = Plane[{0,0,0}, {1,0,0}];

This is the Trajectory that I want to take a section of (you'll see why I call it "bananasection" shortly):

bananasection = 
      MakeTrajectory[hheqns, 
        hhinitial[0.0, -0.121, 0.0, 1/12], {q1,  q2,  p1,  
          p2},  {t,  200}, MaxSteps -> 1500]

PoincareSection calculates the points of intersection (such that p1>0 only) between the Trajectory and the Plane q1=0:

PoincareSection[
        Line[Transpose[Drop[TrajToPoints[bananasection], {3}]]], pplane] // 
      Chop

(An abbreviated listing of the above was given  in the article.)

This is a phase space plot of the Trajectory in q1-q2-p2 space:

figure5a = 
      TrajPlot3D[bananasection, {q1, q2, p2}, 
        PlotRange -> {{-0.5, 0.5}, {-0.5, 0.5}, {-0.35, 
              0.35}}]

This combined figure shows the Trajectory, the Plane and the Points of the Poincaré section:

\(figure5b = 
      Show[Graphics3D[{ToPolygon[pplane], 
            Line[Transpose[
                Drop[TrajToPoints[bananasection], {3}]]], {PointSize[0.02], 
              PoincareSection[
                Line[Transpose[Drop[TrajToPoints[bananasection], {3}]]], 
                pplane]}
		}], Boxed -> False, Axes -> True, 
        Lighting -> True, AmbientLight -> RGBColor[1, 1, 0.3], 
        PlotRange -> {{\(-0.5\), 0.5}, {\(-0.5\), 0.5}, {\(-0.35\), 
              0.35}}]\)

I can use the output from PoincareSection to get a 2-D plot of the points in the q2-p2 plane (and hence the name "bananasection"):

Draw[{PointSize[0.02], 
        Chop[PoincareSection[
              Line[Transpose[Drop[TrajToPoints[bananasection], {3}]]], 
              pplane]] /. Point[{x_, y_, z_}] -> Point[{y, z}]}, 
      Axes -> True, AxesLabel -> {"q2", "p2"}]
\(		\)

These are the other two plots shown in the article, using different choices for the section plane:

\(figure5c = 
      Show[Graphics3D[{ToPolygon[Plane[{0, 0, 0}, {1, 1, 0}]], 
            Line[Transpose[
                Drop[TrajToPoints[bananasection], {3}]]], {PointSize[0.02], 
              PoincareSection[
                Line[Transpose[Drop[TrajToPoints[bananasection], {3}]]], 
                Plane[{0, 0, 0}, {1, 1, 0}]]}
		}], Boxed -> False, 
        Axes -> True, Lighting -> True, AmbientLight -> RGBColor[1, 1, 0.3], 
        PlotRange -> {{\(-0.5\), 0.5}, {\(-0.5\), 0.5}, {\(-0.35\), 
              0.35}}]\)

\(figure5d = 
      Show[Graphics3D[{ToPolygon[Plane[{0, 0, 0}, {0, 1, 0}]], 
            Line[Transpose[
                Drop[TrajToPoints[bananasection], {3}]]], {PointSize[0.02], 
              PoincareSection[
                Line[Transpose[Drop[TrajToPoints[bananasection], {3}]]], 
                Plane[{0, 0, 0}, {0, 1, 0}]]}
		}], Boxed -> False, 
        Axes -> True, Lighting -> True, AmbientLight -> RGBColor[1, 1, 0.3], 
        PlotRange -> {{\(-0.5\), 0.5}, {\(-0.5\), 0.5}, {\(-0.35\), 
              0.35}}]\)

(The following is how to put together the four parts of the figure given in the article:

figure5 = 
      Show[GraphicsArray[{{figure5a, figure5b}, {figure5c, 
              figure5d}}]]

)

5. A Poincaré section movie

(To run a Mathematica animation, click on the second blue bracket to the right-hand-side of the graphic--the one with the hook on the bottom--and then either select Animate Selected Graphics from the Cell menu or hold down the "Ctrl" key and press "Y". You can change the speed or direction of the movie using the buttons at the bottom left of the Mathematica window. To stop, click the mouse button.)

The following code creates a movie of the section Plane sweeping through 360° (if you are running this, it may take a while):

\(Table[
	
      p = Plane[{0, 0, 0}, {Cos[\(-α\)], Sin[\(-α\)], 0}]; 
      Show[Graphics3D[{ToPolygon[p], 
            Line[Transpose[
                Drop[TrajToPoints[bananasection], {3}]]], {PointSize[0.02], 
              PoincareSection[
                Line[Transpose[Drop[TrajToPoints[bananasection], {3}]]], 
                p]}
		}], Boxed -> False, Axes -> True, Lighting -> True, 
        AmbientLight -> RGBColor[1, 1, 0.3], 
        PlotRange -> {{\(-0.5\), 0.5}, {\(-0.5\), 0.5}, {\(-0.35\), 
              0.35}}], {α, 0, 2  π, π/8}]\)