<jupyter_start><jupyter_text>Two-body ProblemsThe gravitational potential energy and forces involving two masses $a$ and $b$ are $$\begin{eqnarray}U_{ab}&=&-\frac{Gm_am_b}{|\boldsymbol{r}_a-\boldsymbol{r}_b|},\\\nonumberF_{ba}&=&-\frac{Gm_am_b}{|\boldsymbol{r}_a-\boldsymbol{r}_b|^2}\hat{r}_{ab},\\\nonumber\hat{r}_{ab}&=&\frac{\boldsymbol{r}_b-\boldsymbol{r}_a}{|\boldsymbol{r}_a-\boldsymbol{r}_b|}.\end{eqnarray}$$ Here $G=6.67\times 10^{-11}$ Nm$^2$/kg$^2$, and $F_{ba}$ is the forceon $b$ due to $a$. By inspection, one can see that the force on $b$due to $a$ and the force on $a$ due to $b$ are equal and opposite. Thenet potential energy for a large number of masses would be $$\begin{equation}U=\sum_{a<b}U_{ab}=\frac{1}{2}\sum_{a\ne b}U_{ab}.\label{_auto1} \tag{1}\end{equation}$$  Relative and Center of Mass MotionThus far, we have considered the trajectory as if the force iscentered around a fixed point. For two bodies interacting only withone another, both masses circulate around the center of mass. Onemight think that solutions would become more complex when bothparticles move, but we will see here that the problem can be reducedto one with a single body moving according to a fixed force byexpressing the trajectories for $\boldsymbol{r}_1$ and $\boldsymbol{r}_2$ into thecenter-of-mass coordinate $\boldsymbol{R}_{\rm cm}$ and the relativecoordinate $\boldsymbol{r}$, $$\begin{eqnarray}\boldsymbol{R}_{\rm cm}&\equiv&\frac{m_1\boldsymbol{r}_1+m_2\boldsymbol{r}_2}{m_1+m_2},\\\nonumber\boldsymbol{r}&\equiv&\boldsymbol{r}_1-\boldsymbol{r_2}.\end{eqnarray}$$ Here, we assume the two particles interact only with one another, so$\boldsymbol{F}_{12}=-\boldsymbol{F}_{21}$ (where $\boldsymbol{F}_{ij}$ is the force on $i$due to $j$. The equations of motion then become $$\begin{eqnarray}\ddot{\boldsymbol{R}}_{\rm cm}&=&\frac{1}{m_1+m_2}\left\{m_1\ddot{\boldsymbol{r}}_1+m_2\ddot{\boldsymbol{r}}_2\right\}\\\nonumber&=&\frac{1}{m_1+m_2}\left\{\boldsymbol{F}_{12}+\boldsymbol{F}_{21}\right\}=0.\\\ddot{\boldsymbol{r}}&=&\ddot{\boldsymbol{r}}_1-\ddot{\boldsymbol{r}}_2=\left(\frac{\boldsymbol{F}_{12}}{m_1}-\frac{\boldsymbol{F}_{21}}{m_2}\right)\\\nonumber&=&\left(\frac{1}{m_1}+\frac{1}{m_2}\right)\boldsymbol{F}_{12}.\end{eqnarray}$$ The first expression simply states that the center of mass coordinate$\boldsymbol{R}_{\rm cm}$ moves at a fixed velocity. The second expressioncan be rewritten in terms of the reduced mass $\mu$. $$\begin{eqnarray}\mu \ddot{\boldsymbol{r}}&=&\boldsymbol{F}_{12},\\\frac{1}{\mu}&=&\frac{1}{m_1}+\frac{1}{m_2},~~~~\mu=\frac{m_1m_2}{m_1+m_2}.\end{eqnarray}$$ Thus, one can treat the trajectory as a one-body problem where thereduced mass is $\mu$, and a second trivial problem for the center ofmass. The reduced mass is especially convenient when one isconsidering gravitational problems because then $$\begin{eqnarray}\mu \ddot{r}&=&-\frac{Gm_1m_2}{r^2}\hat{r}\\\nonumber&=&-\frac{GM\mu}{r^2}\hat{r},~~~M\equiv m_1+m_2.\end{eqnarray}$$ For the gravitational problem, the reduced mass then falls out and thetrajectory depends only on the total mass $M$.The kinetic energy and momenta also have analogues in center-of-masscoordinates. The total and relative momenta are $$\begin{eqnarray}\boldsymbol{P}&\equiv&\boldsymbol{p}_1+\boldsymbol{p}_2=M\dot{\boldsymbol{R}}_{\rm cm},\\\nonumber\boldsymbol{q}&\equiv&\mu\dot{\boldsymbol{r}}.\end{eqnarray}$$ With these definitions, a little algebra shows that the kinetic energy becomes $$\begin{eqnarray}T&=&\frac{1}{2}m_1|\boldsymbol{v}_1|^2+\frac{1}{2}m_2|\boldsymbol{v}_2|^2\\\nonumber&=&\frac{1}{2}M|\dot{\boldsymbol{R}}_{\rm cm}|^2+\frac{1}{2}\mu|\dot{\boldsymbol{r}}|^2\\\nonumber&=&\frac{P^2}{2M}+\frac{q^2}{2\mu}.\end{eqnarray}$$ The standard strategy is to transform into the center of mass frame,then treat the problem as one of a single particle of mass $\mu$undergoing a force $\boldsymbol{F}_{12}$. Scattering angles can also beexpressed in this frame, then transformed into the lab frame. Inpractice, one sees examples in the literature where $d\sigma/d\Omega$expressed in both the "center-of-mass" and in the "laboratory"frame. Deriving Elliptical OrbitsKepler's laws state that a gravitational orbit should be an ellipsewith the source of the gravitational field at one focus. Deriving thisis surprisingly messy. To do this, we first use angular momentumconservation to transform the equations of motion so that it is interms of $r$ and $\theta$ instead of $r$ and $t$. The overall strategyis to1. Find equations of motion for $r$ and $t$ with no angle ($\theta$) mentioned, i.e. $d^2r/dt^2=\cdots$. Angular momentum conservation will be used, and the equation will involve the angular momentum $L$.2. Use angular momentum conservation to find an expression for $\dot{\theta}$ in terms of $r$.3. Use the chain rule to convert the equations of motions for $r$, an expression involving $r,\dot{r}$ and $\ddot{r}$, to one involving $r,dr/d\theta$ and $d^2r/d\theta^2$. This is quitecomplicated because the expressions will also involve a substitution $u=1/r$ so that one finds an expression in terms of $u$ and $\theta$.4. Once $u(\theta)$ is found, you need to show that this can be converted to the familiar form for an ellipse.The equations of motion give $$\begin{eqnarray}\label{eq:radialeqofmotion} \tag{2}\frac{d}{dt}r^2&=&\frac{d}{dt}(x^2+y^2)=2x\dot{x}+2y\dot{y}=2r\dot{r},\\\nonumber\dot{r}&=&\frac{x}{r}\dot{x}+\frac{y}{r}\dot{y},\\\nonumber\ddot{r}&=&\frac{x}{r}\ddot{x}+\frac{y}{r}\ddot{y}+\frac{\dot{x}^2+\dot{y}^2}{r}-\frac{\dot{r}^2}{r}.\end{eqnarray}$$ Recognizing that the numerator of the third term is the velocity squared, and that it can be written in polar coordinates, $$\begin{equation}v^2=\dot{x}^2+\dot{y}^2=\dot{r}^2+r^2\dot{\theta}^2,\label{_auto2} \tag{3}\end{equation}$$ one can write $\ddot{r}$ as $$\begin{eqnarray}\label{eq:radialeqofmotion2} \tag{4}\ddot{r}&=&\frac{F_x\cos\theta+F_y\sin\theta}{m}+\frac{\dot{r}^2+r^2\dot{\theta}^2}{r}-\frac{\dot{r}^2}{r}\\\nonumber&=&\frac{F}{m}+\frac{r^2\dot{\theta}^2}{r}\\\nonumberm\ddot{r}&=&F+\frac{L^2}{mr^3}.\end{eqnarray}$$ This derivation used the fact that the force was radial,$F=F_r=F_x\cos\theta+F_y\sin\theta$, and that angular momentum is$L=mrv_{\theta}=mr^2\dot{\theta}$. The term $L^2/mr^3=mv^2/r$ behaveslike an additional force. Sometimes this is referred to as acentrifugal force, but it is not a force. Instead, it is theconsequence of considering the motion in a rotating (and thereforeaccelerating) frame.Now, we switch to the particular case of an attractive inverse squareforce, $F=-\alpha/r^2$, and show that the trajectory, $r(\theta)$, isan ellipse. To do this we transform derivatives w.r.t. time toderivatives w.r.t. $\theta$ using the chain rule combined with angularmomentum conservation, $\dot{\theta}=L/mr^2$. $$\begin{eqnarray}\label{eq:rtotheta} \tag{5}\dot{r}&=&\frac{dr}{d\theta}\dot{\theta}=\frac{dr}{d\theta}\frac{L}{mr^2},\\\nonumber\ddot{r}&=&\frac{d^2r}{d\theta^2}\dot{\theta}^2+\frac{dr}{d\theta}\left(\frac{d}{dr}\frac{L}{mr^2}\right)\dot{r}\\\nonumber&=&\frac{d^2r}{d\theta^2}\left(\frac{L}{mr^2}\right)^2-2\frac{dr}{d\theta}\frac{L}{mr^3}\dot{r}\\\nonumber&=&\frac{d^2r}{d\theta^2}\left(\frac{L}{mr^2}\right)^2-\frac{2}{r}\left(\frac{dr}{d\theta}\right)^2\left(\frac{L}{mr^2}\right)^2\end{eqnarray}$$ Equating the two expressions for $\ddot{r}$ in Eq.s ([4](eq:radialeqofmotion2)) and ([5](eq:rtotheta)) eliminates all the derivatives w.r.t. time, and provides a differential equation with only derivatives w.r.t. $\theta$, $$\begin{equation}\label{eq:rdotdot} \tag{6}\frac{d^2r}{d\theta^2}\left(\frac{L}{mr^2}\right)^2-\frac{2}{r}\left(\frac{dr}{d\theta}\right)^2\left(\frac{L}{mr^2}\right)^2=\frac{F}{m}+\frac{L^2}{m^2r^3},\end{equation}$$ that when solved yields the trajectory, i.e. $r(\theta)$. Up to thispoint the expressions work for any radial force, not just forces thatfall as $1/r^2$.The trick to simplifying this differential equation for the inversesquare problems is to make a substitution, $u\equiv 1/r$, and rewritethe differential equation for $u(\theta)$. $$\begin{eqnarray}r&=&1/u,\\\nonumber\frac{dr}{d\theta}&=&-\frac{1}{u^2}\frac{du}{d\theta},\\\nonumber\frac{d^2r}{d\theta^2}&=&\frac{2}{u^3}\left(\frac{du}{d\theta}\right)^2-\frac{1}{u^2}\frac{d^2u}{d\theta^2}.\end{eqnarray}$$ Plugging these expressions into Eq. ([6](eq:rdotdot)) gives anexpression in terms of $u$, $du/d\theta$, and $d^2u/d\theta^2$. Aftersome tedious algebra, $$\begin{equation}\frac{d^2u}{d\theta^2}=-u-\frac{F m}{L^2u^2}.\label{_auto3} \tag{7}\end{equation}$$ For the attractive inverse square law force, $F=-\alpha u^2$, $$\begin{equation}\frac{d^2u}{d\theta^2}=-u+\frac{m\alpha}{L^2}.\label{_auto4} \tag{8}\end{equation}$$ The solution has two arbitrary constants, $A$ and $\theta_0$, $$\begin{eqnarray}\label{eq:Ctrajectory} \tag{9}u&=&\frac{m\alpha}{L^2}+A\cos(\theta-\theta_0),\\\nonumberr&=&\frac{1}{(m\alpha/L^2)+A\cos(\theta-\theta_0)}.\end{eqnarray}$$ The radius will be at a minimum when $\theta=\theta_0$ and at amaximum when $\theta=\theta_0+\pi$. The constant $A$ is related to theeccentricity of the orbit. When $A=0$ the radius is a constant$r=L^2/(m\alpha)$, and the motion is circular. If one solved theexpression $mv^2/r=-\alpha/r^2$ for a circular orbit, using thesubstitution $v=L/(mr)$, one would reproduce the expression$r=L^2/(m\alpha)$.The form describing the elliptical trajectory inEq. ([9](eq:Ctrajectory)) can be identified as an ellipse with onefocus being the center of the ellipse by considering the definition ofan ellipse as being the points such that the sum of the two distancesbetween the two foci are a constant. Making that distance $2D$, thedistance between the two foci as $2a$, and putting one focus at theorigin, $$\begin{eqnarray}2D&=&r+\sqrt{(r\cos\theta-2a)^2+r^2\sin^2\theta},\\\nonumber4D^2+r^2-4Dr&=&r^2+4a^2-4ar\cos\theta,\\\nonumberr&=&\frac{D^2-a^2}{D+a\cos\theta}=\frac{1}{D/(D^2-a^2)-a\cos\theta/(D^2-a^2)}.\end{eqnarray}$$ By inspection, this is the same form as Eq. ([9](eq:Ctrajectory)) with $D/(D^2-a^2)=m\alpha/L^2$ and $a/(D^2-a^2)=A$.Let us remind ourselves about what an ellipse is before we proceed.<jupyter_code>%matplotlib inline

import numpy as np
from matplotlib import pyplot as plt
from math import pi

u=1.     #x-position of the center
v=0.5    #y-position of the center
a=2.     #radius on the x-axis
b=1.5    #radius on the y-axis

t = np.linspace(0, 2*pi, 100)
plt.plot( u+a*np.cos(t) , v+b*np.sin(t) )
plt.grid(color='lightgray',linestyle='--')
plt.show()<jupyter_output><empty_output><jupyter_text>Effective or Centrifugal PotentialThe total energy of a particle is $$\begin{eqnarray}E&=&U(r)+\frac{1}{2}mv_\theta^2+\frac{1}{2}m\dot{r}^2\\\nonumber&=&U(r)+\frac{1}{2}mr^2\dot{\theta}^2+\frac{1}{2}m\dot{r}^2\\\nonumber&=&U(r)+\frac{L^2}{2mr^2}+\frac{1}{2}m\dot{r}^2.\end{eqnarray}$$ The second term then contributes to the energy like an additionalrepulsive potential. The term is sometimes referred to as the"centrifugal" potential, even though it is actually the kinetic energyof the angular motion. Combined with $U(r)$, it is sometimes referredto as the "effective" potential, $$\begin{eqnarray}U_{\rm eff}(r)&=&U(r)+\frac{L^2}{2mr^2}.\end{eqnarray}$$ Note that if one treats the effective potential like a real potential, one would expect to be able to generate an effective force, $$\begin{eqnarray}F_{\rm eff}&=&-\frac{d}{dr}U(r) -\frac{d}{dr}\frac{L^2}{2mr^2}\\\nonumber&=&F(r)+\frac{L^2}{mr^3}=F(r)+m\frac{v_\perp^2}{r},\end{eqnarray}$$ which is indeed matches the form for $m\ddot{r}$ in Eq. ([4](eq:radialeqofmotion2)), which included the **centrifugal** force.The following code plots this effective potential for a simple choice of parameters, with a standard gravitational potential $-\alpha/r$. Here we have chosen $L=m=\alpha=1$.<jupyter_code># Common imports
import numpy as np
from math import *
import matplotlib.pyplot as plt

Deltax = 0.01
#set up arrays
xinitial = 0.3
xfinal = 5.0
alpha = 1.0   # spring constant
m = 1.0   # mass, you can change these
AngMom = 1.0  #  The angular momentum
n = ceil((xfinal-xinitial)/Deltax)
x = np.zeros(n)
for i in range(n):
    x[i] = xinitial+i*Deltax
V = np.zeros(n)
V = -alpha/x+0.5*AngMom*AngMom/(m*x*x)
# Plot potential
fig, ax = plt.subplots()
ax.set_xlabel('r[m]')
ax.set_ylabel('V[J]')
ax.plot(x, V)
fig.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Gravitational force exampleUsing the above parameters, we can now study the evolution of the system using for example the velocity Verlet method.This is done in the code here for an initial radius equal to the minimum of the potential well.  We seen then that the radius is always the same and corresponds to a circle (the radius is always constant).<jupyter_code># Common imports
import numpy as np
import pandas as pd
from math import *
import matplotlib.pyplot as plt
import os

# Where to save the figures and data files
PROJECT_ROOT_DIR = "Results"
FIGURE_ID = "Results/FigureFiles"
DATA_ID = "DataFiles/"

if not os.path.exists(PROJECT_ROOT_DIR):
    os.mkdir(PROJECT_ROOT_DIR)

if not os.path.exists(FIGURE_ID):
    os.makedirs(FIGURE_ID)

if not os.path.exists(DATA_ID):
    os.makedirs(DATA_ID)

def image_path(fig_id):
    return os.path.join(FIGURE_ID, fig_id)

def data_path(dat_id):
    return os.path.join(DATA_ID, dat_id)

def save_fig(fig_id):
    plt.savefig(image_path(fig_id) + ".png", format='png')


# Simple Gravitational Force   -alpha/r
    
DeltaT = 0.01
#set up arrays 
tfinal = 100.0
n = ceil(tfinal/DeltaT)
# set up arrays for t, v and r
t = np.zeros(n)
v = np.zeros(n)
r = np.zeros(n)
# Constants of the model, setting all variables to one for simplicity
alpha = 1.0
AngMom = 1.0  #  The angular momentum
m = 1.0  # scale mass to one
c1 = AngMom*AngMom/(m*m)
c2 = AngMom*AngMom/m
rmin = (AngMom*AngMom/m/alpha)
# Initial conditions
r0 = rmin
v0 = 0.0
r[0] = r0
v[0] = v0
# Start integrating using the Velocity-Verlet  method
for i in range(n-1):
    # Set up acceleration
    a = -alpha/(r[i]**2)+c1/(r[i]**3)
    # update velocity, time and position using the Velocity-Verlet method
    r[i+1] = r[i] + DeltaT*v[i]+0.5*(DeltaT**2)*a
    anew = -alpha/(r[i+1]**2)+c1/(r[i+1]**3)
    v[i+1] = v[i] + 0.5*DeltaT*(a+anew)
    t[i+1] = t[i] + DeltaT
    # Plot position as function of time
fig, ax = plt.subplots(2,1)
ax[0].set_xlabel('time')
ax[0].set_ylabel('radius')
ax[0].plot(t,r)
ax[1].set_xlabel('time')
ax[1].set_ylabel('Velocity')
ax[1].plot(t,v)
save_fig("RadialGVV")
plt.show()<jupyter_output><empty_output><jupyter_text>Changing the value of the initial position to a value where the energy is positive, leads to an increasing radius with time, a so-called unbound orbit. Choosing on the other hand an initial radius that corresponds to a negative energy and different from the minimum value leads to a radius that oscillates back and forth between two values.  Harmonic Oscillator in two dimensionsConsider a particle of mass $m$ in a 2-dimensional harmonic oscillator with potential $$U=\frac{1}{2}kr^2=\frac{1}{2}k(x^2+y^2).$$ If the orbit has angular momentum $L$, we can find the radius and angular velocity of the circular orbit as well as the b) the angular frequency of small radial perturbations.We consider the effective potential. The radius of a circular orbit is at the minimum of the potential (where the effective force is zero).The potential is plotted here with the parameters $k=m=0.1$ and $L=1.0$.<jupyter_code># Common imports
import numpy as np
from math import *
import matplotlib.pyplot as plt

Deltax = 0.01
#set up arrays
xinitial = 1.0
xfinal = 5.0
k = 0.1   # spring constant
m = 0.1   # mass, you can change these
AngMom = 1.0  #  The angular momentum
n = ceil((xfinal-xinitial)/Deltax)
x = np.zeros(n)
for i in range(n):
    x[i] = xinitial+i*Deltax
V = np.zeros(n)
V = 0.5*k*x*x+0.5*AngMom*AngMom/(m*x*x)
# Plot potential
fig, ax = plt.subplots()
ax.set_xlabel('r[m]')
ax.set_ylabel('V[J]')
ax.plot(x, V)
fig.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>$$\begin{eqnarray*}U_{\rm eff}&=&\frac{1}{2}kr^2+\frac{L^2}{2mr^2}\end{eqnarray*}$$ The effective potential looks like that of a harmonic oscillator forlarge $r$, but for small $r$, the centrifugal potential repels theparticle from the origin. The combination of the two potentials has aminimum for at some radius $r_{\rm min}$. $$\begin{eqnarray*}0&=&kr_{\rm min}-\frac{L^2}{mr_{\rm min}^3},\\r_{\rm min}&=&\left(\frac{L^2}{mk}\right)^{1/4},\\\dot{\theta}&=&\frac{L}{mr_{\rm min}^2}=\sqrt{k/m}.\end{eqnarray*}$$ For particles at $r_{\rm min}$ with $\dot{r}=0$, the particle does notaccelerate and $r$ stays constant, i.e. a circular orbit. The radiusof the circular orbit can be adjusted by changing the angular momentum$L$.For the above parameters this minimum is at $r_{\rm min}=1$. Now consider small vibrations about $r_{\rm min}$. The effective spring constant is the curvature of the effective potential. $$\begin{eqnarray*}k_{\rm eff}&=&\left.\frac{d^2}{dr^2}U_{\rm eff}(r)\right|_{r=r_{\rm min}}=k+\frac{3L^2}{mr_{\rm min}^4}\\&=&4k,\\\omega&=&\sqrt{k_{\rm eff}/m}=2\sqrt{k/m}=2\dot{\theta}.\end{eqnarray*}$$ Here, the second step used the result of the last step from part(a). Because the radius oscillates with twice the angular frequency,the orbit has two places where $r$ reaches a minimum in onecycle. This differs from the inverse-square force where there is oneminimum in an orbit. One can show that the orbit for the harmonicoscillator is also elliptical, but in this case the center of thepotential is at the center of the ellipse, not at one of the foci.The solution is also simple to write down exactly in Cartesian coordinates. The $x$ and $y$ equations of motion separate, $$\begin{eqnarray*}\ddot{x}&=&-kx,\\\ddot{y}&=&-ky.\end{eqnarray*}$$ So the general solution can be expressed as $$\begin{eqnarray*}x&=&A\cos\omega_0 t+B\sin\omega_0 t,\\y&=&C\cos\omega_0 t+D\sin\omega_0 t.\end{eqnarray*}$$ The code here finds the solution for $x$ and $y$ using the code we developed in homework 4.<jupyter_code>DeltaT = 0.01
#set up arrays 
tfinal = 10.0
n = ceil(tfinal/DeltaT)
# set up arrays
t = np.zeros(n)
v = np.zeros((n,2))
r = np.zeros((n,2))
radius = np.zeros(n)
# Constants of the model
k = 0.1   # spring constant
m = 0.1   # mass, you can change these
omega02 = sqrt(k/m)  # Frequency
AngMom = 1.0  #  The angular momentum
rmin = (AngMom*AngMom/k/m)**0.25
# Initial conditions as compact 2-dimensional arrays
#x0 =rmin*0.5; y0 = sqrt(rmin*rmin-x0*x0)
x0 = 1.0; y0= 1.0
r0 = np.array([x0,y0]) 
v0 = np.array([0.0,0.0])
r[0] = r0
v[0] = v0
# Start integrating using the Velocity-Verlet  method
for i in range(n-1):
    # Set up the acceleration
    a =  -r[i]*omega02  
    # update velocity, time and position using the Velocity-Verlet method
    r[i+1] = r[i] + DeltaT*v[i]+0.5*(DeltaT**2)*a
    anew = -r[i+1]*omega02  
    v[i+1] = v[i] + 0.5*DeltaT*(a+anew)
    t[i+1] = t[i] + DeltaT
# Plot position as function of time
radius = np.sqrt(r[:,0]**2+r[:,1]**2)
fig, ax = plt.subplots(3,1)
ax[0].set_xlabel('time')
ax[0].set_ylabel('radius squared')
ax[0].plot(t,r[:,0]**2+r[:,1]**2)
ax[1].set_xlabel('time')
ax[1].set_ylabel('x position')
ax[1].plot(t,r[:,0])
ax[2].set_xlabel('time')
ax[2].set_ylabel('y position')
ax[2].plot(t,r[:,1])

fig.tight_layout()
save_fig("2DimHOVV")
plt.show()<jupyter_output><empty_output><jupyter_text>With some work using double angle formulas, one can calculate $$\begin{eqnarray*}r^2&=&x^2+y^2\\\nonumber&=&(A^2+C^2)\cos^2(\omega_0t)+(B^2+D^2)\sin^2\omega_0t+(AB+CD)\cos(\omega_0t)\sin(\omega_0t)\\\nonumber&=&\alpha+\beta\cos 2\omega_0 t+\gamma\sin 2\omega_0 t,\\\alpha&=&\frac{A^2+B^2+C^2+D^2}{2},~~\beta=\frac{A^2-B^2+C^2-D^2}{2},~~\gamma=AB+CD,\\r^2&=&\alpha+(\beta^2+\gamma^2)^{1/2}\cos(2\omega_0 t-\delta),~~~\delta=\arctan(\gamma/\beta),\end{eqnarray*}$$ and see that radius oscillates with frequency $2\omega_0$. Thefactor of two comes because the oscillation $x=A\cos\omega_0t$ has twomaxima for $x^2$, one at $t=0$ and one a half period later.The following code shows first how we can solve this problem using the radial degrees of freedom only.<jupyter_code>DeltaT = 0.01
#set up arrays 
tfinal = 10.0
n = ceil(tfinal/DeltaT)
# set up arrays for t, v and r
t = np.zeros(n)
v = np.zeros(n)
r = np.zeros(n)
E = np.zeros(n)
# Constants of the model
AngMom = 1.0  #  The angular momentum
m = 0.1
k = 0.1
omega02 = k/m
c1 = AngMom*AngMom/(m*m)
c2 = AngMom*AngMom/m
rmin = (AngMom*AngMom/k/m)**0.25
# Initial conditions
r0 = rmin
v0 = 0.0
r[0] = r0
v[0] = v0
E[0] = 0.5*m*v0*v0+0.5*k*r0*r0+0.5*c2/(r0*r0)
# Start integrating using the Velocity-Verlet  method
for i in range(n-1):
    # Set up acceleration
    a = -r[i]*omega02+c1/(r[i]**3)    
    # update velocity, time and position using the Velocity-Verlet method
    r[i+1] = r[i] + DeltaT*v[i]+0.5*(DeltaT**2)*a
    anew = -r[i+1]*omega02+c1/(r[i+1]**3)
    v[i+1] = v[i] + 0.5*DeltaT*(a+anew)
    t[i+1] = t[i] + DeltaT
    E[i+1] = 0.5*m*v[i+1]*v[i+1]+0.5*k*r[i+1]*r[i+1]+0.5*c2/(r[i+1]*r[i+1])
    # Plot position as function of time
fig, ax = plt.subplots(2,1)
ax[0].set_xlabel('time')
ax[0].set_ylabel('radius')
ax[0].plot(t,r)
ax[1].set_xlabel('time')
ax[1].set_ylabel('Energy')
ax[1].plot(t,E)
save_fig("RadialHOVV")
plt.show()<jupyter_output><empty_output><jupyter_text>Stability of OrbitsThe effective force can be extracted from the effective potential, $U_{\rm eff}$. Beginning from the equations of motion, Eq. ([2](eq:radialeqofmotion)), for $r$, $$\begin{eqnarray}m\ddot{r}&=&F+\frac{L^2}{mr^3}\\\nonumber&=&F_{\rm eff}\\\nonumber&=&-\partial_rU_{\rm eff},\\\nonumberF_{\rm eff}&=&-\partial_r\left[U(r)+(L^2/2mr^2)\right].\end{eqnarray}$$ For a circular orbit, the radius must be fixed as a function of time,so one must be at a maximum or a minimum of the effectivepotential. However, if one is at a maximum of the effective potentialthe radius will be unstable. For the attractive Coulomb force theeffective potential will be dominated by the $-\alpha/r$ term forlarge $r$ because the centrifugal part falls off more quickly, $\sim1/r^2$. At low $r$ the centrifugal piece wins and the effectivepotential is repulsive. Thus, the potential must have a minimumsomewhere with negative potential. The circular orbits are then stableto perturbation.The effective potential is sketched for two cases, a $1/r$ attractivepotential and a $1/r^3$ attractive potential. The $1/r$ case has astable minimum, whereas the circular orbit in the $1/r^3$ case isunstable.If one considers a potential that falls as $1/r^3$, the situation isreversed and the point where $\partial_rU$ disappears will be a localmaximum rather than a local minimum. **Fig to come here with code**The repulsive centrifugal piece dominates at large $r$ and the attractiveCoulomb piece wins out at small $r$. The circular orbit is then at amaximum of the effective potential and the orbits are unstable. It isthe clear that for potentials that fall as $r^n$, that one must have$n>-2$ for the orbits to be stable.Consider a potential $U(r)=\beta r$. For a particle of mass $m$ withangular momentum $L$, find the angular frequency of a circularorbit. Then find the angular frequency for small radial perturbations.For the circular orbit you search for the position $r_{\rm min}$ where the effective potential is minimized, $$\begin{eqnarray*}\partial_r\left\{\beta r+\frac{L^2}{2mr^2}\right\}&=&0,\\\beta&=&\frac{L^2}{mr_{\rm min}^3},\\r_{\rm min}&=&\left(\frac{L^2}{\beta m}\right)^{1/3},\\\dot{\theta}&=&\frac{L}{mr_{\rm min}^2}=\frac{\beta^{2/3}}{(mL)^{1/3}}\end{eqnarray*}$$ Now, we can find the angular frequency of small perturbations about the circular orbit. To do this we find the effective spring constant for the effective potential, $$\begin{eqnarray*}k_{\rm eff}&=&\partial_r^2 \left.U_{\rm eff}\right|_{r_{\rm min}}\\&=&\frac{3L^2}{mr_{\rm min}^4},\\\omega&=&\sqrt{\frac{k_{\rm eff}}{m}}\\&=&\frac{\beta^{2/3}}{(mL)^{1/3}}\sqrt{3}.\end{eqnarray*}$$ If the two frequencies, $\dot{\theta}$ and $\omega$, differ by aninteger factor, the orbit's trajectory will repeat itself each timearound. This is the case for the inverse-square force,$\omega=\dot{\theta}$, and for the harmonic oscillator,$\omega=2\dot{\theta}$. In this case, $\omega=\sqrt{3}\dot{\theta}$,and the angles at which the maxima and minima occur change with eachorbit. Code example with gravitional forceThe code example here is meant to illustrate how we can make a plot of the final orbit. We solve the equations in polar coordinates (the example here uses the minimum of the potential as initial value) and then we transform back to cartesian coordinates and plot $x$ versus $y$. We see that we get a perfect circle when we place ourselves at the minimum of the potential energy, as expected.<jupyter_code># Simple Gravitational Force   -alpha/r
    
DeltaT = 0.01
#set up arrays 
tfinal = 8.0
n = ceil(tfinal/DeltaT)
# set up arrays for t, v and r
t = np.zeros(n)
v = np.zeros(n)
r = np.zeros(n)
phi = np.zeros(n)
x = np.zeros(n)
y = np.zeros(n)
# Constants of the model, setting all variables to one for simplicity
alpha = 1.0
AngMom = 1.0  #  The angular momentum
m = 1.0  # scale mass to one
c1 = AngMom*AngMom/(m*m)
c2 = AngMom*AngMom/m
rmin = (AngMom*AngMom/m/alpha)
# Initial conditions, place yourself at the potential min
r0 = rmin
v0 = 0.0  # starts at rest
r[0] = r0
v[0] = v0
phi[0] = 0.0
# Start integrating using the Velocity-Verlet  method
for i in range(n-1):
    # Set up acceleration
    a = -alpha/(r[i]**2)+c1/(r[i]**3)
    # update velocity, time and position using the Velocity-Verlet method
    r[i+1] = r[i] + DeltaT*v[i]+0.5*(DeltaT**2)*a
    anew = -alpha/(r[i+1]**2)+c1/(r[i+1]**3)
    v[i+1] = v[i] + 0.5*DeltaT*(a+anew)
    t[i+1] = t[i] + DeltaT
    phi[i+1] = t[i+1]*c2/(r0**2)
# Find cartesian coordinates for easy plot    
x = r*np.cos(phi)
y = r*np.sin(phi)
fig, ax = plt.subplots(3,1)
ax[0].set_xlabel('time')
ax[0].set_ylabel('radius')
ax[0].plot(t,r)
ax[1].set_xlabel('time')
ax[1].set_ylabel('Angle $\cos{\phi}$')
ax[1].plot(t,np.cos(phi))
ax[2].set_ylabel('y')
ax[2].set_xlabel('x')
ax[2].plot(x,y)

save_fig("Phasespace")
plt.show()<jupyter_output><empty_output><jupyter_text>Try to change the initial value for $r$ and see what kind of orbits you get.In order to test different energies, it can be useful to look at the plot of the effective potential discussed above.However, for orbits different from a circle the above code would need modifications in order to allow us to display say an ellipse. For the latter, it is much easier to run our code in cartesian coordinates, as done here. In this code we test also energy conservation and see that it is conserved to numerical precision. The code here is a simple extension of the code we developed for homework 4.<jupyter_code># Common imports
import numpy as np
import pandas as pd
from math import *
import matplotlib.pyplot as plt

DeltaT = 0.01
#set up arrays 
tfinal = 10.0
n = ceil(tfinal/DeltaT)
# set up arrays
t = np.zeros(n)
v = np.zeros((n,2))
r = np.zeros((n,2))
E = np.zeros(n)
# Constants of the model
m = 1.0   # mass, you can change these
alpha = 1.0
# Initial conditions as compact 2-dimensional arrays
x0 = 0.5; y0= 0.
r0 = np.array([x0,y0]) 
v0 = np.array([0.0,1.0])
r[0] = r0
v[0] = v0
rabs = sqrt(sum(r[0]*r[0]))
E[0] = 0.5*m*(v[0,0]**2+v[0,1]**2)-alpha/rabs
# Start integrating using the Velocity-Verlet  method
for i in range(n-1):
    # Set up the acceleration
    rabs = sqrt(sum(r[i]*r[i]))
    a =  -alpha*r[i]/(rabs**3)
    # update velocity, time and position using the Velocity-Verlet method
    r[i+1] = r[i] + DeltaT*v[i]+0.5*(DeltaT**2)*a
    rabs = sqrt(sum(r[i+1]*r[i+1]))
    anew = -alpha*r[i+1]/(rabs**3)
    v[i+1] = v[i] + 0.5*DeltaT*(a+anew)
    E[i+1] = 0.5*m*(v[i+1,0]**2+v[i+1,1]**2)-alpha/rabs
    t[i+1] = t[i] + DeltaT
# Plot position as function of time
fig, ax = plt.subplots(3,1)
ax[0].set_ylabel('y')
ax[0].set_xlabel('x')
ax[0].plot(r[:,0],r[:,1])
ax[1].set_xlabel('time')
ax[1].set_ylabel('y position')
ax[1].plot(t,r[:,0])
ax[2].set_xlabel('time')
ax[2].set_ylabel('y position')
ax[2].plot(t,r[:,1])

fig.tight_layout()
save_fig("2DimGravity")
plt.show()
print(E)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17216389154028353, 'fraction_numerical': 0.027998727330575884, 'mean_word_length': 4.986878306878307, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 34, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '127560', 'license': 'CC0-1.0', 'n_tokens_mistral': 10470, 'n_tokens_neox': 9996, 'n_words': 2960, 'path': 'doc/src/LectureNotes/_build/jupyter_execute/testbook/chapter6.ipynb', 'repo_name': 'anacost/MachineLearning', 'size': '28287'}