Note: If you made rapid progress on the previous exercises, then you're encouraged to give this one a try. But if not (or if this one takes a long time), remember it's ok not to finish.
Consider integrating the equations of motion for a test particle orbiting a star fixed at the origin. In 2D, there are two 2nd order ODEs:
$$\frac{d^2 r_x}{dt^2} = a_x = -GM r_x/r^3$$
$$\frac{d^2 r_y}{dt^2} = a_y = -GM r_y/r^3,$$
where $r = \sqrt(r_x^2+r_y^2)$ and $a_x$ and $a_y$ are the accelerations in the x and y directions.
Numerically, it is generally easier to integrate 1st order ODEs, so we transform these into:
\begin{eqnarray}
d v_x/dt & = a_x & = -GM r_x/r^3 \\
d v_y/dt & = a_y & = -GM r_y/r^3 \\
d r_x/dt & = v_x & \\
d r_y/dt & = v_y & \\
\end{eqnarray}
Euler’s method is a simplistic algorithm for integrating a set of 1st order ODEs with state stored in a vector x:
$$
r(t_n+1) = r(t_n) + dt * \frac{dr}{dt}|_{t_n} \\
v(t_n+1) = v(t_n) + dt * \frac{dv}{dt}|_{t_n}
$$
a) Write a function integrate_euler!(state,dt, duration) that integrates a system (described by state) using Euler’s method and steps of size dt for an interval of time given by duration. It would probably be wise to write at least two additional functions. Have your function return a two dimensional array containing a list of the states of the system during your integration. NOTE: ADD EXAMPLE OF 2D ARRAY SYNTAX.
b) Use the above code to integrate a system with an initial state of (r_x, r_y, v_x, v_y) = (1, 0, 0, 1) for roughly 3 orbits (i.e., 3*2pi time units if you set GM=1) using approximately 200 time steps per orbit. Inspect the final state and comment on the accuracy of the integration.
c) Plot the trajectory of the test particle (e.g., x vs y). Is the calculated behavior accurate?
d) Write a new function integrate_leapfrog!(state,dt, duration) that is analogous to integrate_euler, but uses the “leapfrog” or modified Euler’s integration method: $$ r(t_n+1/2) = r(t_n) + dt/2 * \frac{dr}{dt}|_{t_n} \\ v(t_n+1) = v(t_n) + dt * \frac{dv}{dt}|_{t_n+1/2} \\ r(t_n+1) = r(t_n+1/2) + dt/2 * \frac{dr}{dt}|_{t_n+1} $$
e) Repeat the integration from part 4b, but using the integrate_leapfrog code. Inspect the final end state and make a similar plot. Describe the main difference in the results and explain on what basis you would choose one over the other.
f) Time how long your code for a similar integration of 100 orbits. How long would it take to integrate for 4.5*10^9 orbits (e.g., age of Earth)?
g) Repeat the integration from part 4e, but using a larger dt. Inspect the final end state. How much less accurate is the integration that used a larger dt?
h) Make a log-log plot of the accuracy of the integration versus the integration duration for durations spanning at least a few orders of magnitude. Estimate the slope (by eye is fine). What are the implications for the accuracy of a long-term integration?
i) List a few ideas for how you could accelerate the integrations. Comment on how likely each of your ideas is to make a significant difference. (If you try it out, report the new time for 100 orbits, so we can compare.) Discuss any drawbacks of each idea to accelerate the integrations.