Morin 1.7 (Two masses, one swinging)

As you can see in the book, Morin 1.7 concerns the standard issue pulley problem, with two equal masses
connected by a rope draped over a pulley.  While one mass moves up and down only, the other one swings
back and forth, obeying these equations of motion:

eqn17r= 2 r''[t] == r[t]( θ'[t])^2 - g (1-Cos[θ[t]])
eqn17θ = θ''[t] == - 2 r'[t] θ'[t]/r[t] - g Sin[θ[t]]/r[t]

To solve using Euler's integration method, assign starting values to r, r1 (meaning r-dot), q and q1 (meaning θ and θ-dot).

r=1
r1=0
q = 3.14/18
q1 = 0

Choose a step size, define g, and initialize the counter.

ϵ = .001
g = 9.8
i = 0

Learn about the While function:

?While

The "body" can consist of a whole string of Mathematica expressions, seperated by semicolons.
For extra clarity I added some parentheses around the body, though they are not necessary:

While[r<2.0,
	(i = i+1;
	r2 = (r*q1^2 - g*(1-Cos[q]))/2.;
	r = r+ ϵ*r1;
	r1 = r1 + ϵ*r2;
	q2 = -2*r1*q1/r - g*Sin[q]/r;
	q = q+ϵ*q1;
	q1 = q1 + ϵ*q2;)
]
i*ϵ

The loop runs until r exceeds 2, and we extract the result by counting the steps.

Your job: repeat for step sizes .0001 and .00001, comparing to the book if available.
Note that Morin chooses to write π as 3.14, and g as 9.8, so we are looking for an accuracy of  maybe a percent or better.
In other words, the stepsize is small enough when we get 3 digits stable.

Also: use NDSolve[] to solve the above eqn17r and eqn17θ with the given boundary conditions.
Again find the time at which r[t]==2, perhaps defining a myr[t] and using FindRoot[].

Morin 1.17 (Pendulum)-- assigned homework

For the simple pendulum (mass m at the end of a string L) the equation of motion for the angle is:

eqn117 = θ''[t]== (-g/L) Sin[θ[t]]

Your job: write a program like the one above to answer the question: how long does it take to swing
from θ=π/2 down to vertical (θ=0), choosing L=1, g=9.8.  Also use NDSolve to get a more accurate answer.
Print the few cells you need to show your work, and hand in with the homework.  (Don't print this whole notebook!)

Morin 7.4

Let's find some orbits with a central force of the form β r^k.
Start with the effective potential.  If you defined L above, we clear it (and others) now with:

L=.
r=.
k=.

Ueff[r_] =  L^2/ ( 2 m  r^2) + β r^k

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

Can we DSolve?!?!  You may want to learn about the "Abort Evaluation" menu item under "Evaluation" above.

DSolve[reqn,r[t],t]

So let's choose specific numbers

p = {m->1,L->1,k->7,β->1}

And Plot the potential:

Plot[Ueff[r]/.p,{r,.4,1.3}]

Next find the mininum of this potential

sUmin = FindMinimum[Ueff[r]/.p,{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+ .2

And add this chosen energy to the plot:

Plot[{Ueff[r],e}/.p,{r,.4,1.2}]

We might as well locate the turning points:

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

r2 = (r/.FindRoot[(Ueff[r]/.p)==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)/.p

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

Now we are ready to solve:

nsoln = NDSolve[{reqn/.p,r'[0]==rdot0,r[0]==r0},r[t],{t,0,10}]

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

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

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

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

θeqn = θ'[t] == (L/(m (myr[t])^2))/.p

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

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.

p1=ParametricPlot[{myr[t]*Cos[myθ[t]],myr[t]*Sin[myθ[t]]},{t,0,10},AspectRatio->1,PlotRange->{{-1.5,1.5},{-1.5,1.5}}]

And of course we need a movie:

orbit[tt_] := Show[{ParametricPlot[{myr[t]*Cos[myθ[t]],myr[t]*Sin[myθ[t]]},{t,0,tt},DisplayFunction->Identity,PlotRange->{{-1.5,1.5},{-1,2}},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,10}]

Your job: repeat for different values of the constant k and the energy e, looking for orbits which close upon themselves.
The book claims (and you are to show in your own words in the homework) that for small amplitude oscillations 
(i.e. e close to e0) there are certain special values of k (e.g. k=7), for which the frequency of radial oscillation is a rational multiple of the orbital freqency.
And that means the orbit eventually returns to it starting position.  Which value of k gives the "spring ellipse," even for large values of energy?