The last post detailed how specific types of information could be gleaned from the equation of motion for a $$1/r$$ central potential. Three conserved quantities were discovered: the energy, the angular momentum, and the eccentricity vector. In addition, the radial equation, which gives the trajectory’s distance from the scattering center as a function of polar angle along the trajectory, presented itself after a clever manipulation.
This column presents specific techniques and results for modeling individual scattering trajectories for both repulsive ($$\alpha < 0 $$) and attractive ($$\alpha > 0$$) forces.
The first step in calculating specific, individual trajectories is the transformation of the equation of motion from the usual Newton form
\[ \ddot {\vec r} + \frac{\alpha \vec r}{r^3} = 0 \]
into a state space form. The most straightforward method is to use Cartesian coordinates, which yields the coupled first-order equations
\[ \frac{d}{dt} \left[ \begin{array}{c} x \\ y \\ z \\ v_x \\ v_y \\ v_z \end{array} \right] = \left[ \begin{array}{c} v_x \\ v_y \\ v_z \\ -\frac{\alpha x}{r^3} \\ -\frac{\alpha y}{r^3} \\ -\frac{\alpha z}{r^3} \end{array} \right] \; .\]
Numerical integration of the above system results in an ephemeris representation of the trajectory – a listing of discrete times and corresponding position components $$x, y, z$$ and velocity components $$v_x, v_y, v_z$$. Nonlinear combinations of the state will be used to determine the conserved quantities both as a check on the numerical method but also to emphasize the power affording by finding and exploiting conserved quantities.
All of the numerical modeling was done using Python 2.7, numpy, and scipy under the Jupyter notebook format. Matplotlib is used for all of the plots.
The core function used to integrate the equations of motion is
def TwoBody(state,t,params):
#unpack the state
x = state[0]
y = state[1]
z = state[2]
vx = state[3]
vy = state[4]
vz = state[5]
#unpack params
alpha = params['alpha']
#define helpers
r = sp.sqrt(x**2 + y**2 + z**2)
q = alpha/r**3
#state derive
ds = [vx, vy, vz, -q*x, -q*y, -q*z]
return ds
Note its structural similarity to the coupled equation above – specifically, the function TwoBody returns the right-hand side of the system of differential equations. The scipy.integrate function odeint, which provides the numerical integration technique, takes this function, arrays specifying the initial conditions and the time span, and any additional arguments that need to be passed. In this case, the additional argument is a Python dictionary that contains the value of $$\alpha$$.
The odeint function returns an array of state values at the times contained in the array that specifies the time span. Together the two arrays comprise the ephemeris of the system.
One additional, useful ingredient is the function conserved_quantities
def conserved_quantities(soln,params):
pos = soln[:,0:3]
vel = soln[:,3:6]
r = np.sqrt((pos*pos).sum(axis=1)).reshape(-1,1) #makes the broadcast shape correct
E = 0.5*(vel*vel).sum(axis=1).reshape(-1,1) - params['alpha']/r #other reshape needed to avoid NxN shape for E
L = np.cross(pos,vel)
Ecc = np.cross(vel,L)/params['alpha'] - pos/r
return E, L, Ecc
that takes the state array and returns arrays corresponding to the energy, the angular momentum vector, and the eccentricity vector. The conserved quantities serve two purposes. The first is that they help make sense out of the jumble of state variables. This was the point of last week’s column. The second is they provide indicators of the goodness of the numerical techniques. If they are not well conserved, then the numerical techniques are not doing what they should.
The first test case run was for the initial conditions
\[ \bar S_0 = \left[ \begin{array}{c} -30 \\ 0.1 \\ 0 \\ 2.25 \\ 0 \\ 0 \end{array} \right] \]
and a value of $$\alpha = -0.4$$. Since $$\alpha < 0 $$, the scattering is repulsive. The resulting trajectory is
The black dot indicates the force or scattering center, the blue dot the point of closest approach, and the red dot location of the eccentricity vector.
While the trajectory looks like what we would expect, the is no obvious way to tell, from this figure, whether the ephemeris is correct. Based on the initial conditions, the conserved quantities should be
\[ E = 2.54458 \; ,\]
\[ \vec L = \left[ \begin{array}{c} 0 & 0 & -0.225 \end{array} \right]^T \; ,\]
and
\[ \vec E = \left[ \begin{array}{c} 1 & -1.268958 & 0 \end{array} \right]^T \; .\]
The figure below shows that error in each quantity, where the error is defined as the difference between the instantaneous value and the initial value (e.g. $$\epsilon_E = E(t) – E_0$$).
It is clear that the integration does an excellent job of honoring the conservation of the energy, angular momentum, and the eccentricity vector. So now we can believe that the ephemeris is accurate and can test the radial equation. The next figure shows the distance of the trajectory as a function of time in the neighborhood of the closest approach.
From the figure, the closest approach occurs at a distance of $$r = 0.2055838$$. The radial equation, which reads
\[ r(\nu) = \frac{L^2/|\alpha|}{e \cos\nu – 1} \; ,\]
predicts that the closest approach, which occurs at $$\nu = 0$$, of $$r_{CA} = \frac{L^2/|\alpha|}{e-1} = 0.2055838$$.
In the same fashion, we can generate a trajectory for $$\alpha = 0.4$$ with the same initial conditions as above. The resulting trajectory is
Based on the initial conditions, the conserved quantities should be
\[ E = 2.51792 \; ,\]
\[ \vec L = \left[ \begin{array}{c} 0 & 0 & -0.225 \end{array} \right]^T \; ,\]
and
\[ \vec E = \left[ \begin{array}{c} 1 & 1.262292 & 0 \end{array} \right]^T \; .\]
The corresponding figure for the time evolution of the conserved quantities is
The radial equation predicts $$r_{CA} = 0.04848406$$ compared to a value of $$0.04848647$$ from the plot
Next week, I’ll be looking at the behavior of an ensemble of trajectories and will be deriving the classical scattering cross section for Rutherford scattering and examining how well the numerical experiments agree.