Phys 262 Math session--Jan 26

Matrix Manipulations

Suppose we want to solve the simultaneous equations

eq1 = 3 x - 5 y == 5
eq2 = x - y == 2

Your task: Solve[] for x and y.

Next, phrase the problem in matrix form.  I'll type some random numbers here---you fix it up to match the above:

matrix = {{1,2},{3,4}}

MatrixForm[matrix]

Here is how to get every output in matrix form:

$Post = MatrixForm

b = {1,0}

Then here is the system we want to solve, written on neatly on the screen: (this is only for formatting purposes; the “equation” doesn’t make sense to Mathematica in this form)

MatrixForm[matrix].MatrixForm[{x,y}] == MatrixForm[b]

Once you have fixed up the matrix and the vector b, the solution can be obtained using the Inverse function:

Inverse[matrix]

Inverse[matrix].b

For good measure, perhaps you want to check your homework by repeating for the system

eq1 = x1 + 2x2 == 9
eq2 = 3 x1 + 4 x2 == 23

And here is how you turn off the MatrixForm output:

$Post = .

Celestial orbits

Here are the incantations which have appeared repeatedly in class already.
Start by writing down the effective potential, which gives the diffeq for r.
We choose specific numerical parameters so we can use numerical integration, by default all are 1 in some units:

l = m = α=  1

Ueff[r_] =  -α/r + l^2/ ( 2 m  r^2)

reqn = m r''[t] == - D[Ueff[r[t]],r[t]]

Next find the mininum of this potential

sUmin = FindMinimum[Ueff[r],{r,.8}]

and save the information in the symbols r0 and e0

r0 = r /. sUmin[[2]]

e0 = sUmin[[1]]

Choose an energy for our motion somewhat bigger than the minimum:

e = e0+ .05

And finally make our diagram:

Plot[{Ueff[r],e},{r,.1 r0,3 r0}]

We might as well locate the turning points:

r1 = (r/.FindRoot[Ueff[r]==e,{r,.5r0}])

r2 = (r/.FindRoot[Ueff[r]==e,{r,2r0}])

Now to find the actual motion we need to set up initial conditions. 
Let us agree to launch at r[0]=r0 with whatever speed we need to make
the energy right:

eeqn = (1/2) m(rdot)^2 + Ueff[r0] == e

rdot

rdot0 = (rdot/.(Solve[eeqn,rdot])[[2]])

rdot0=0

Now we are ready to solve:

nsoln = NDSolve[{reqn,r'[0]==rdot0,r[0]==r0},r[t],{t,0,100}]

And save a copy of the solution as "myr[t]"

myr[t_] = (r[t]/.nsoln[[1]])

Plot[myr[t],{t,0,100}]

Now that we have r[t], find θ[t].

θeqn = θ'[t] == l/(m (myr[t])^2);

θsoln = NDSolve[{θeqn,θ[0]==0},θ[t],{t,0,100}]

myθ[t_] = θ[t]/.θsoln[[1]]

Plot[myθ[t],{t,0,10}]

And now we can draw the curve.  Note that you might want to play with the aspect ratio.

ParametricPlot[{myr[t]*Cos[myθ[t]],myr[t]*Sin[myθ[t]]},{t,0,100},AspectRatio->Automatic]

And to make the movie:

orbit[tt_] := Show[{ParametricPlot[{myr[t]*Cos[myθ[t]],myr[t]*Sin[myθ[t]]},{t,0,tt},DisplayFunction->Identity,PlotRange->{{-3,3},{-3,3}},AspectRatio->1],
ListPlot[{{myr[tt]*Cos[myθ[tt]],myr[tt]*Sin[myθ[tt]]}},PlotStyle->{PointSize[.03],RGBColor[1,0,0]},DisplayFunction->Identity]},DisplayFunction->$DisplayFunction]

Animate[orbit[t],{t,0,20}]

Your optional tasks:

Repeat for a variety of energies.  Choose a small one and satisfy yourself that r[t] looks sinusoidal; choose a large one and see that it is not.

Repeat for a slightly perturbed potential, and notice that the orbits no longer close

Ueff[r_] = - c/r + l^2/ ( 2 m  r^2) + (1/100)/r^3

Repeat for an α r^k potential for a variety of k's.  Goal: reproduce the plots on p.303 of Morin.

Homework Task: Morin 7.12

We are asked: suppose the earth's tangential velocity suddenly dropped to zero.
How long would it take to fall in to the Sun?
The straighforward way of formulating the problem is to say that at time t=0 we start at r(0)=R,
where R is the usual 1.5x10^11 meters to the sun, and we obey the F=ma equation which says:

fma = m r''[t] == -α/r[t]^2

With great hope, one might try to DSolve, but Mathematica gets us only halfway there (i.e. it can find t(r), but not r(t)).
Instead choosing the numerical parameters as above (α=1,m=1) and using the appriate r0 (which is then 1),
try NDSolve for the motion.  Find the time at which r[t] hits zero, i.e. we crash in to the Sun.
(This is treating the Sun as a point, of course.) 
Tricky point #1: you will get complaints from NDSolve that something funny happens at a particular time.  That time is of course when r=0 the force blows up.  Cure #1: ignore the complaint; Cure #2: adjust your integration time to just short of the singularity.
Tricky point #2: you've found a numerical value for the time, but what are the units?  Answer:  units such that the original orbit at r0=1 takes one year.  And what was the numerical value of the period above?  (The period being the time it takes for θ to advance by 2π.)
Tricky point #3: the book problem hints that you can do this by thinking of falling in as the limiting case of a very eccentric ellipse.  That is, we want (half of) the period of an ellipse with major axis r0.  Since any orbit with major axis 2r0 (such as Earth) takes one year, Kepler's third tells us that (half of) the period of the skinny ellipse is ... (what?!)

Your tasks: Construct and print cells showing the NDSolve and the resulting time, as well as a quick "N[some sqrt to some power] involving 2's and π's"  to demonstrate that the result of the NDSolve agrees with point #3.

(BONUS segment): Halley's Comet

Let' s make a movie of Halley' s comet.  Start by carefully drawing the Earth' s orbit, including nonzeor eccentricity.
  Choose parameters by convention :

α = 1
m =1

Take as inputs the period T (measured in years) and the eccentricity ϵ, e.g. for Earth,

T = 1; ϵ=.017;

To convert to the parameters L and E we use Kepler's law:

s  = Solve[T^2 == π^2 m α^2/(-2 e^3),e]

We want the real solution:

e0 = e/. s[[3]]

To get L

s=Solve[ϵ == Sqrt[1 + 2 e0 L^2/(m α^2)],L]

L0 = L/.s[[2]]

With L in hand we can define the effective potential:

Ueff[r_] = - α/r +L0^2/ ( 2 m  r^2)

We can draw a line on the graph with:

Plot[{Ueff[r],e0},{r,.1,2}]

Find the minimum:

s=FindMinimum[Ueff[r],{r,.5}]

Choose as our starting r0 that radius:

r0 = r/.s[[2]]

And figure out the appropriate v0, i.e. the one consistent with the energy (e0) we have.

s=Solve[(1/2) m v^2 + Ueff[r0]== e0,v]

v0 = v/.s[[1]]

Now that we have the boundary conditions, let's solve F=ma:

reqn = m r''[t] == - D[Ueff[r[t]],r[t]]

s= NDSolve[{reqn,r[0]==r0,r'[0]==v0},r[t],{t,0,80}]

rE = r[t]/.s[[1]]

Plot the r[t] we've got and see if the period really is one year:

Plot[rE,{t,0,2}]

And with r[t] in hand, we can NDSolve for θ[t] using the fact that

θeqn = θ'[t] == L0/(m (rE)^2)

s=NDSolve[{θeqn,θ[0]==0},θ[t],{t,0,80}]

θE = θ[t]/.s[[1]]

Then make a position vector function:

pE[t_]= {rE*Cos[θE], rE*Sin[θE]}

And plot

ParametricPlot[pE[t],{t,0,1},AspectRatio->1]

Now repeat all these steps for Halley's comet:

T = 76; ϵ=.967;

s  = Solve[T^2 == π^2 m α^2/(-2 e^3),e]

e1 = e/.s[[2]]

s=Solve[ϵ == Sqrt[1 + 2 e1 L^2/(m α^2)],L]

L1 = L/.s[[2]]

Ueff[r_] = - α/r +L1^2/ ( 2 m  r^2)

Plot[{Ueff[r],e1},{r,.1,2}]

s=FindMinimum[Ueff[r],{r,.4}]

r1 = r/.s[[2]]

s=Solve[(1/2) m v^2 + Ueff[r1]== e1,v]

v1 = v/.s[[2]]

reqn = m r''[t] == - D[Ueff[r[t]],r[t]]

s= NDSolve[{reqn,r[0]==r1,r'[0]==v1},r[t],{t,0,80}]

rH = r[t]/.s[[1]]

θeqn = θ'[t] == L0/(m (rH)^2)

s=NDSolve[{θeqn,θ[0]==0},θ[t],{t,0,80}]

θH= θ[t]/.s[[1]]

pH[t_]= {rH*Cos[θH], rH*Sin[θH]}

ParametricPlot[pH[t],{t,0,76}]

And now define a function to draw dots at the locations of Earth and Halley's comet:

drawstate2[t_]  := Show[
ListPlot[{pE[t],pH[t]},AspectRatio->1,PlotStyle->{Hue[1],PointSize[.02]},PlotRange->{{-6,6},{-1,11}}],
ParametricPlot[{pE[tprime],pH[tprime]},{tprime,0,t}]]

And animate

?*Animate*

Morin 7.17

This problem asks about the angular momentum (Lz), given that a satellite is in parabolic orbit with a given distance of closest approach (L).  Here is a plot, choosing L=1, with a red dot at the focus, and a generic blue point.

p1 = Plot[x^2/4,{x,0,4}]; p2 = ListPlot[{{0,1}},PlotStyle->{Red,PointSize[.03]}];
p3 =  ListPlot[{{x,x^2/4}/.x->3.8},PlotStyle->{Blue,PointSize[.02]}];
Show[p1,p2,p3]

Part (a) of the problem asks you to compute Lz as the satellites passes through the origin, noting that you know the speed from the fact that the total energy is exactly zero.  Part (b) asks you to express Lz using formulae for r0 (k) and ϵ.  Part (c) asks you to compute Lz at a generic point (x0,y0) on the orbit, simplifying various expressions under the assumption that x0 is much larger than L.  Here let’s use Mathematica to compute Lz exactly.

First express y0 in terms of x0

y0 = x0^2/(4 L)

The displacement vector from the red origin to the blue point is:

r = {x0,y0-L}

which has length

dist = Sqrt[r.r]

In terms of this distance we know the kinetic energy (again since total energy is exactly 0), and the speed:

speed = Sqrt[(2/m) α/dist]

And this speed is the magnitude of a velocity vector which is directed at an angle with slope:

tanθ = D[y0,x0]

From tan(θ) we can figure out the sin and cos as:

sinθ = tanθ/Sqrt[1+tanθ^2]

cosθ = 1/Sqrt[1+tanθ^2]

So the velocity vector is:

v=speed {cosθ,sinθ}

Neatened up:

FullSimplify[v,Assumptions->{α>0,x0>0,m>0,L>0}]

and the angular momentum is this cross product:  (note the fancy way we compute this!)

Det[{ r,m v}]

which simplifies to .... (drum roll)

FullSimplify[%,Assumptions->{α>0,x0>0,m>0,L>0}]

Your job on the homework is to do the same steps where you assume x0 is large, and keep only the leading (dominant) contribution in each expression.  For example the exact distance is

FullSimplify[dist,Assumptions->{α>0,x0>0,m>0,L>0}]

which means the potential energy is some factors times 1/x0^2 plus terms of order L^2/x0^2 and higher, which we drop.  Here is how to get Mathematica to express this:

Series[dist,{x0,Infinity,2}]

Similarly, for the speed ...

Series[speed,{x0,Infinity,2}]//(FullSimplify[#,Assumptions->{α>0,x0>0,m>0,L>0}]&)

and the whole velocity vector:

Series[v,{x0,Infinity,2}]//(FullSimplify[#,Assumptions->{α>0,x0>0,m>0,L>0}]&)

Since Lz is exactly the same for any x0, it should come as no surprise that we get the same result for Lz order by order in powers of (1/x0).