Execises (4)

Assigned: February 14, 2013
Due: (Thursday) February 21, 2013

Setup your Solution Notebook (in-class)

Just like last week's exercises, you must set-up your notebook in order to solve these execises. Copy the definitions for the equations of motion for Hyperion orbiting Saturn.

Equation of Motion for Center of Mass

Since the separation between the two dumbell masses is small, the equations of motion for Hyperion's center-of-mass (COM) orbit about Saturn are identical to the motion of planets about the Sun. Below, we normalize time to the mean period of orbit (21 days) and distance to the mean separation between the center of Saturn and the COM of Hyperion (0.01 AU).

Off[General::spell1];

pCOM = {x[t], y[t], vx[t], vy[t]};

dpCOM = {vx[t], vy[t],
-4 (π^2) x[t]/r^3, -4 (π^2) y[t]/r^3}/.  r-> Sqrt[x[t]^2 + y[t]^2]

p0COM = {x0, y0,vx0, vy0};

Equation of Motion for Rotation

Hyperion's rate of change of angular rotation is determined by the torque applied by Saturn.

Some simple vector algebra

a = {ax,ay,az};
b = {bx, by, bz};

a \[Cross] b

a . b

The torque

Let "d" be the distance from the COM to one of the two masses.

r1 = {d Cos[θ[t]], d Sin[θ[t]],0};
r2 =  {d Cos[θ[t]+Pi], d Sin[θ[t]+Pi],0};
rCOM = {x[t],y[t],0};

f1 = -4(π^2) m(r1 + rCOM)/((r1 + rCOM).(r1 + rCOM) )^(3/2)

f2= -4(π^2) m(r2 + rCOM)/((r2 + rCOM).(r2 + rCOM) )^(3/2)

torque = r1 \[Cross] f1 + r2 \[Cross] f2 // Simplify

Let's consider only the case where "d" is very small...

Series[torque[[3]],{d,0,2}] // Simplify

We see the leading order term is proportional to (d^2).

The Leading-Order Expression for Rotational Dynamics

moment = 2 m (d^2); (* Moment of Inertial *)
pROT = {θ[t],ω[t]};
dpROT = {ω[t],
		-(24 m (d^2)(π^2) (Sin[θ[t]] x[t]-Cos[θ[t]] y[t]) (Cos[θ[t]] x[t]+Sin[θ[t]] y[t]))/(((((x[t])^2)+((y[t])^2)))^(5/2))/moment};
p0ROT = {θ0, ω0};

Combined Equations of Motion

p = pCOM ~ Join ~ pROT;
p0 = p0COM ~ Join ~ p0ROT;
dp = dpCOM ~ Join ~ dpROT;

(Thread[...] is an easy way of writing a single equation for a list of variables.)

eqs = Thread[D[p,t] == dp]~ Join ~Thread[{x[0], y[0], vx[0], vy[0], θ[0],ω[0]} == p0]

h1 = NDSolve[eqs /.x0 -> 1 /. y0 -> 0 /. vy0->2 π/. vx0 -> 0/.θ0 -> 0 /. ω0 -> 1, p, {t, 0,10}]//Flatten

h2 = NDSolve[eqs /.x0 -> 1 /. y0 -> 0 /. vy0->2 π /1.25/. vx0 -> 0/.θ0 -> 0 /. ω0 -> 1, p, {t, 0,10}]//Flatten

Plot[Mod[Evaluate[θ[t]-ArcTan[x[t],y[t]]/.h2]+π,2 π],{t,0,10},PlotLabel->"Radial Orientation of Hyperion",AxesLabel->{"time","Δθ"}]

Eccentricity (in-class)

A command is defined below that can be used to find the eccentricity of an orbit. To use this command, you must make a table of the square of the radius of the orbit for at least one orbital period. (Remember, you can find this using "(x[t]^2 + y[t]^2 /. solution)"Use this command to compute the eccentricity of the circular and elliptical orbits demonstrated in class (and called "h1" and "h2").

eccentricty[r2_List] := Module[{mn,mx},
		mn = Sqrt[ Min[r2]];
		mx = Sqrt[Max[r2]];
		(mx-mn)/(mx + mn) (* e = c/a *)]

r1List = Table[Evaluate[x[t]^2 + y[t]^2 /. h1], {t, 0, 2, 0.05}];

r2List = Table[Evaluate[x[t]^2 + y[t]^2 /. h2], {t, 0, 2, 0.05}];

r2List // Shallow

eccentricty[r1List]

eccentricty[r2List]

Lyapunov (in-class)

The command below takes a list of the differences (of the logrythm of the magnitude) between two orbits and fits these list of points to a straight line. Demonstrate the use of this function by recreating a table of the logrythm of the magnitude of two solutions with slightly different initial conditions. (Use the solution "h2a" illustrated in class.)

lyapunov[dt_List] := Module[{ft},
		ft = Fit[dt,{1,t}, t];
		Coefficient[ft,t]]

h2a = NDSolve[eqs /.x0 -> 1 /. y0 -> 0 /. vy0->2 π /1.25/. vx0 -> 0/.θ0 -> 0.001 /. ω0 -> 1, p, {t, 0,10}]//Flatten

dl[t_] = Log[Abs[(θ[t] /. h2a) - (θ[t] /. h2)]]

dl[1.2]

dList = Table[{t, dl[t]}, {t, 0, 10, 0.1}];

dList // Shallow

ListPlot[dList]

ListPlot[Select[dList,Last[#1]<0&]]

We only want to send the "linear" part of the data list to the lyapunov[...] command. This is the part of the data list that can be fit to a straight line.

lyapunov[Select[dList, (Last[#]<0)&]]

Measure Chaos versus Eccentricity (in-class)

When Hyperion is in a circular orbit, then the orientation of the long axis is stable. When the orbit is sufficiently elliptical, the orientation becomes unstable and "flips" randomly.

In this assignment, you are to use the computer to estimate how elliptical the orbit must be to create chaotic tumbling. This is when the difference between orientations of two nearly identical initial conditions grows exponentially fast.

Execute the commands described by this notebook for different eccentricities of Hyperion's orbit. As the eccentricity increases, calculate the Lyapunov number using the method previously demonstrated. Can you estimate for what values of eccentricity the moon tumbles chaotically? If you can, create a plot of the Lyapunov exponent as a function of eccentricity.

To help your study, the "getLyapunov" command is defined below. This is command first solves two orbits beginning with the initial velocity, "v0". The second orbit begin "nearly" the same as the first. Second, the command computes the eccentricity of the orbit. Finally, it compares the two closely-spaced orbits, and computes a Lyapunov exponent. The command returns a list of two numbers: the eccentricity and an estimate of the Lyapunov exponent.

getLyapunov[v0_Real] :=Module[{h1,h2,r2,dt},h1 = (* solution 1 *)NDSolve[eqs /.x0 -> 1 /. y0 -> 0 /. vy0->v0/. vx0 -> 0/.θ0 -> 0 /. ω0 -> 1, p, {t, 0,5}] // Flatten;	h2 = (* solution 2 *)NDSolve[eqs /.x0 -> 1 /. y0 -> 0 /. vy0->v0/. vx0 -> 0/.θ0 -> 0.01 /. ω0 -> 1, p, {t, 0,5}] // Flatten;
r2 = (* list of radii for eccentricity *)Table[x[t]^2 + y[t]^2 /. h1, {t, 0, 5, 0.1}];dt = (* how far apart? *)
Table[{t,Log[Abs[(θ[t] /. h1) - (θ[t] /. h2)]]},{t,0,5,0.05}];
{eccentricty[r2],lyapunov[Select[dt, (Last[#] < 0)&]]}]

When you are finished, you should find Hyperion becomes chaotic (i.e. unstable) when the eccentricity exceeds a critical value (about 0.35 for our model problem).

getLyapunov[N[2 π]]

getLyapunov[N[2 π]/1.1]

el = Table[getLyapunov[N[2 π]/a], {a, 1, 1.5, 0.01}];

ListPlot[el,PlotStyle->PointSize[0.02`],Frame->True,FrameLabel->{"Eccentricity","Lyapunov"}]

Chirikov Map

The Russian mathematician and physicist, Chirikov, studied the dynamics of coupled oscillators in 1979. Instead of solving a differential equation, he defined a "map" or a "stepper" that advances the system according to a model. 
His map is often called the "circle map", and it is relatively easy to model on the computer. The map has a single parameter, "k", that controls the interaction between the drive and the natural oscillations. When k ~ 0, then the motion is purely periodic. As k exceeds unity, the motion is chaotic. 
In Mathematica, we use the Chirikov map with NestList[...] instead of NDSolve[…], and this is computationally efficient.

One simple model for the Chirikov's map is:

chirikov[k_][{th_,p_}] := Module[{thp, pp},
		pp = Mod[p - k Sin[th], 2 Pi];
		thp = Mod[th + pp, 2 Pi];
		{thp, pp}];

Note that this map has a parameter, "k". To use the map, you must assign a value to k.

chirikov[0.5][{π/2, 3 π/4}]

NestList[chirikov[0.5], {π/2, 3 π/4}, 5]

ListPlot[NestList[chirikov[0.5],{(π/2),((3 π)/4)},1000],PlotRange->{{0,2 π},{0,2 π}},Frame->True,FrameLabel->{"θ","p"}]

ListPlot[NestList[chirikov[1.5],{(π/2),((3 π)/4)},1000],PlotRange->{{0,2 π},{0,2 π}},Frame->True,FrameLabel->{"θ","p"}]

In this exercise, you are to write a module that takes a parameter k and an initial condition and computes the Lyapunov exponent for the Chirikov map. When the Lyapunov exponent becomes large, the system is chaotic. Determine, approximately, the value of k when this system becomes chaotic. (With some very  sophisticated mathematics, it can be shown that the critical value of k leading to global chaos is k = 0.9716.)

Perhaps the best way to illustrate the onset of chaos is to compute the Lyapunov exponent as k slowly increases. If you choose this method, you should start two trajectories very close together (differing by perhaps a few hundreds of a percent) and compute how quickly (or slowly) these points move away from each other. Note, however, the important point: the trajectories can never be more than about π apart! Thus, look for exponential divergence only when the separation of the two trajectories are less than unity.

Solution

Phase - Space Trajectories Depend on Initial Condition

Show[Table[ListPlot[NestList[chirikov[0.85],i{(π/5),(π/5)},1000],PlotRange->{{0,2 π},{0,2 π}},Frame->True,FrameLabel->{"θ","p"}, PlotStyle -> {PointSize[Tiny],Hue[i/10.0]}], {i, 1, 9}], PlotLabel -> "Chirikov Parameter = 0.85"]

Show[Table[ListPlot[NestList[chirikov[1.5],i{(π/5),(π/5)},5000],PlotRange->{{0,2 π},{0,2 π}},Frame->True,FrameLabel->{"θ","p"}, PlotStyle -> {PointSize[Tiny],Hue[i/10.0]}], {i, 1, 9}], PlotLabel -> "Chirikov Parameter = 1.5"]

Example with k = 0.5...

l1 = NestList[chirikov[0.5], {π/2.0,3π/4.0}, 500];
l2 = NestList[chirikov[0.5], {π/2.0  1.00001,3π/4.0}, 500];
d = l1 - l2;
d = Sqrt[#1^2 + #2^2]& @@@ d;

ListPlot[Log[Abs[d]]]

ListPlot[Select[Log[Abs[d]],(# < -1)&]]

Fit[Select[Log[Abs[d]],(# < -1)&],{1,t}, t]

Example with k = 1.5...

l1 = NestList[chirikov[1.5], {π/2.0,3π/4.0}, 500];
l2 = NestList[chirikov[1.5], {π/2.0  1.00001,3π/4.0}, 500];
d = l1 - l2;
d = Sqrt[#1^2 + #2^2]& @@@ d;

ListPlot[Log[Abs[d]]]

ListPlot[{Select[Log[Abs[d]],(# < -1)&],Select[Log[Abs[d]],(# < -3)&]}, PlotStyle -> {{},{Red}}]

Fit[Select[Log[Abs[d]],(# < -3)&],{1,t}, t]

Table of Results

lyapunovC[k_][{th0_, p0_}] := Module[{l1, l2, d, ft},
l1 = NestList[chirikov[k], {th0,p0}, 500];
l2 = NestList[chirikov[k], {th0 1.00001,p0}, 500];
d = l1 - l2;
d = Sqrt[#1^2 + #2^2]& @@@ d;
ft = Fit[Select[Log[Abs[d]],(# < -3)&],{1,t}, t];
Coefficient[ft,t]]

lyapunovC[0.5][{π/2, 3π/4}]

lTable = Table[{k,lyapunovC[k][{π/2, 3π/4}]}, {k, 0.5, 1.5, 0.005}];

ListPlot[lTable,PlotStyle-> PointSize[0.01], Frame -> True,PlotLabel-> "Lyapunov Exponent versus k", AxesOrigin-> {0.5, 0.0}]

The Lyapunov number increases rapidly near the onset of gloabl chaos at the predicted value of k = 0.9716…