Bound Orbit (numerical solution)

[brainpushups]

Homework Statement

Consider a particle with mass m and angular momentum l in the field of a central force F=$\frac{-k}{r^{5/2}}$. To simplify your equations, choose units for which m=l=k=1. a) find the value $r_{0}$ of r at which $U_{eff}$ is a minimum and make a plot of $U_{eff}$(r) for 0<r<5$r_{0}$ b) Assuming now that the particle has energy E=-0.1 find an accurate value of $r_{min}$, the particle's distance of closest approach to the force center. c) Assuming that the particle is a r=$r_{min}$ when ψ=0 use a computer program to solve the transformed radial equation and find the orbit in the form r=r(ψ) for 0≤ψ≤7π. Plot the orbit. Does it appear to be closed?

Homework Equations

The transformed radial equation: u''(ψ)=-u(ψ)-$\frac{μ}{l^{2}u(ψ)^{2}}$F comes from the radial equation of the orbit μ$r^{..}$=F(r)+$\frac{l^{2}}{μ r^{3}}$ by making the substitution u=1/r.

The Attempt at a Solution

I have completed everything but the final plot of the orbit and I'm not sure why it is giving me trouble.

To do part a) simply find where the derivative of the effective potential is equal to zero and solve for r. Using the suggested simplifications for the units it is easy to find that

$r_{0}$=1.

I plotted the effective potential in Mathematica along with the energy of the particle by defining $U_{eff}$(r)=-$\frac{2}{3 r^{3/2}}$+$\frac{1}{2r^{2}}$ and then using

Plot[{Ueff[r],-0.1},{r,0,5}].

From the plot it is easy to see that the distance of closest approach is around 0.7 and if you find the roots you will see that it is $r_{0}$=0.6671.

Now, for part c) you need to write F(r) as F(u) using the substitution mentioned above. After simplification u''(ψ)=$\sqrt{u(ψ)}$-2u(ψ).

If this can be solved for u(ψ) then you can easily find r(ψ). I tried using Mathematica again to numerically solve and then plot the reciprocal of my solution. Part of the problem is that I wasn't sure what to use for my second initial condition, u'(0), but I tried several different numbers just to see if I could get a reasonable plot and I did not.

Anyway, here is what I entered

s=NDSolve[{u''[ψ]==$\sqrt{u[ψ]}$-2u[ψ], u[0]=1.499, u'[0]=1}, u, {ψ, 0, 7π}]

Here, u[0] comes from the fact that u=1/r and r0=0.667. Once again, I just picked a number for u'[0] but I have tried various selections.

To plot the solution

PolarPlot[Evaluate[1/u[ψ]/.s] ,{ψ,0,7π}]

The result is not the orbit (not closed) that is pictured in the back of the book, it is more or less a straight line. Not sure where I went wrong here and any help would be appreciated. To check that the technique was correct I did check my work with other examples that I knew the solution to (a free particle, and a particle under the influence of an inverse square force). I was able to get correct plots using the method described above for these.

Homework Statement

Homework Equations

The Attempt at a Solution

 

[brainpushups]

After playing around with Mathematica I found a plot that seemed to match the answer (this is problem 8.25 out of Classical Mechanics by Taylor). If you enter as the transformed radial equation

u''[ψ]=$\sqrt{u[ψ]}$-u[ψ] instead of u''[ψ]=$\sqrt{u[ψ]}$-2u[ψ] and use an initial condition that u'[0]=0 the orbit is indistinguishable from the one in the back of the book (on the given scale, at least).

However, it certainly doesn't appear that this is the correct transformed radial equation if you work through the math.

We have u''=-$\frac{μ}{l^{2}u^{2}}$F(u)-u

u''=-$\frac{μ}{l^{2}u^{2}}$(-u$^{5/2}$+$u^{3}$)-u

With µ=1 l=1 we have

u''=$u^{1/2}$-u-u

u''=$u^{1/2}$-2u

What gives?