Exercise 3

Forward integration

In class01, we wrote a julia program orbit.jl that simulates a two-dimensional gravitating system. It describes a test mass \(\mathsf{m}\) in orbit around a heavy mass \(\mathsf{M}\), following the usual \(\mathsf{1/r^2}\) law and also subject to a velocity-dependent drag force:

\[\mathsf{m\frac{d^2\vec{r}}{dt^2}} = \mathsf{m\ddot{\vec{r}}} = \mathsf{-\frac{GMm}{r^2}\hat{r} - C_0 m\vec{v}} = \mathsf{-\frac{GMm}{r^3}\vec{r} - C_0 m\vec{v}.}\]

Here, \(\mathsf{\vec{r}(t)} = \mathsf{x\hat{x} + y\hat{y}}\) is the position vector, and \(\mathsf{\vec{v}(t)} = \mathsf{\dot{\vec{r}}} = \mathsf{\dot{x}\hat{x} + \dot{y}\hat{y}} = \mathsf{v_x\hat{x} + v_y\hat{y}}\) is the vector velocity. In units where GM = 1 and C_0 = 0.01, the initial conditions are \(\mathsf{\vec{r}(0)} = \mathsf{\hat{x}}\) and \(\mathsf{\vec{v}(0)} = \mathsf{\hat{y}}\) (corresponding to the planet starting at the “eastern” edge of the orbit and moving “north”). These are chosen to produce a near-circular orbit that slowly spirals in toward the central mass \(\mathsf{M}\).

In the code, we don’t use vectors. Rather, we define individual floating-point variables defined for the \(\mathsf{x}\) and \(\mathsf{y}\) components of the position (x and y), velocity (vx and \mathsf{vy}), and acceleration (ax and ay). Breaking the vector differential equation into components gives

\[\begin{aligned} \mathsf{\dot{x}} &= \mathsf{v_x},\\ \mathsf{\dot{y}} &= \mathsf{v_y},\\ \mathsf{\ddot{x}} &= \mathsf{a_x} = \mathsf{-GM \frac{x}{(x^2+y^2)^{3/2}} - C_0 v_x}, \\ \mathsf{\ddot{y}} &= \mathsf{a_y} = \mathsf{-GM \frac{y}{(x^2+y^2)^{3/2}} - C_0 v_y.} \end{aligned}\]

Integrating over small finite steps of duration \(\mathsf{\Delta t}\) produces

\[\begin{aligned} \mathsf{\int_{t}^{t+\Delta t}\!dt'\,\dot{\vec{r}}(t')} &= \mathsf{\vec{r}(t+\Delta t) - \vec{r}(t)},\\ \mathsf{\int_{t}^{t+\Delta t}\!dt'\,\ddot{\vec{r}}(t')} &= \mathsf{\dot{\vec{r}}(t+\Delta t) - \dot{\vec{r}}(t)} \\ &= \mathsf{-\int_{t}^{t+\Delta t}\!dt'\,\biggl(\frac{GM\vec{r}(t')}{r(t')^3} + C_0 \vec{v}(t') \biggr)}\\ &= \mathsf{-\biggl(\frac{GM\vec{r}(t^*)}{r(t^*)^3} + C_0 \vec{v}(t^*) \biggr)\Delta t}, \end{aligned}\]

where the final equality is guaranteed to hold (by a theorem of calculus, very easy to justify geometrically) for some instant in time satisfying \(\mathsf{t < t^* < t + \Delta t}\).

Verlet algorithm

A quick and dirty (but surprisingly effective) approximation is to put \(\mathsf{t^* = t + \Delta t/2}\) at the interval midpoint and to update the position in two half-steps; this leads to the famous Verlet algorithm, viz.,

\[\begin{aligned} \mathsf{x(t^*)} &:= \mathsf{x(t) + v_x(t)(\Delta t/2)}\\ \mathsf{y(t^*)} &:= \mathsf{y(t) + v_y(t)(\Delta t/2)}\\ \mathsf{v_x(t+\Delta t)} &:= \mathsf{v_x(t) - \biggl(\frac{GMx(t^*)}{[x(t^*)^2+y(t^*)^2]^3} + C_0 v_x(t) \biggr)\Delta t},\\ \mathsf{v_y(t+\Delta t)} &= \mathsf{v_y(t) - \biggl(\frac{GMy(t^*)}{[x(t^*)^2+y(t^*)^2]^3} + C_0 v_y(t) \biggr)\Delta t},\\ \mathsf{x(t+\Delta t)} &= \mathsf{x(t^*) + v_x(t+\Delta) (\Delta t/2)}\\ \mathsf{y(t+\Delta t)} &= \mathsf{y(t^*) + v_y(t+\Delta) (\Delta t/2)}\\ \end{aligned}\]

Some things to try

Try generalizing the program from two to three dimensions, and give the test mass a small out-of-plane initial velocity, \(\mathsf{z(0) = 0.1}\). Let’s call that code orbit3D.jl.

Try converting your julia program to python. Confirm that orbit3D.jl and the new orbit3D.py both produce the same trajectory.

$ julia orbit3D.jl > traj1.dat
$ python3 orbit3D.py > traj2.dat
$ gnuplot
> plot for [n=1:2] "traj".n.".dat" using 1:2 with lines

Notice that the stopping condition (r < R0) terminates the program when the test mass hits the surface of the “planet” (we model the large mass \(\mathsf{M}\) as a sphere or radius \(\mathsf{R_0}\)). What if that stopping condition is removed? In physical terms, what if \(\mathsf{M}\) is a point mass with no spatial extent? See what the outcome is and try to relate it to your understanding of the floating-point representation.