ASTR415 Spring 2007 Problem Set #5

Solution

The code for my integration package is given at the end. It consists of a source file integ.c with routines for Euler, LeapFrog, and 4th-order Runge-Kutta (the latter adapated from NRiC's rk4()), a header file integ.h with prototypes and a macro for specifying float or double precision, and a Makefile for building the implementation specific to this problem set. The routines in integ.c roughly follow the NRiC rk4() calling style. Each function expects a derivatives function to be supplied (the prototype is the same for each). I have left provision for arbitrary data to be passed to this function via a generic (void *) pointer (this will be useful for PS6).

1. To reduce this 2nd-order ODE to 1st order, define a new variable $v = dx/dt$ so the system of equations becomes:
   
   $$\frac{dx}{dt} = v, \hspace{0.2in} \frac{dv}{dt} = -x.$$
   
   
   The code for integrating this system is given at the end (integ1.c). The code uses various #defines at the start to specify the range of $t$ and the step sizes $h$ to use. Since $t$ will be incremented by $h$ many times, both are stored in double precision. Even so, a small offset is used when comparing $t$ to the end time ($t$ = 15) to ensure the final step is performed (this is important as the value of $x$ at this step is used to determine the error). Outputs are handled by combining a prefix identifying the integration scheme with the step size to form a unique filename. Internally the variables $x$ and $v$ are stored in a single array x[] (as defined in main()), in that order, and their derivatives $\dot{x}$ and $\dot{v}$ are stored in dxdt[].
   1. Plots for the three methods and five step sizes are given at the end. The solid red line in each plot is the analytical solution $x = \sin(t)$ while the crosses show the numerical solution at each step. The Euler integrator never seems stable, in the sense that all step sizes appear (by eye) to give diverging solutions. The other methods, LeapFrog and 4th-order Runge-Kutta (RK4), appear stable for all step sizes tested. However, plots of $x - \sin(t)$ (not shown) indicate that all methods diverge from the true solution to some extent. But we know LeapFrog is symplectic, so what must be happening in the LeapFrog case is phase drift, not amplitude drift (Euler and RK4 show both phase and amplitude errors, but over the time interval tested only Euler goes amplitude unstable, in the sense of deviating from the allowed range of the solution, i.e., -1 to 1, by more than $\sim$ 1%). Another way to see this is to sort the numerical values $x$ and analytical values $\sin(t)$ and plot their difference. This separates the phase error from the amplitude error. Euler always diverges but LeapFrog and RK4 remain bounded.

   2. Plots of the position error are shown with the previous results. The error in the Euler method drops roughly as $h$, LeapFrog as $h^2$, and RK4 as $h^4$ (note that more than 6 digits of output are required to see the proper trend at the smallest step sizes for RK4). At first it might seem these slopes are too shallow by one power of $h$ in each case, since a method of order $m$ has a truncation error of order $m+1$, where $m$ = 1 for Euler, 2 for LeapFrog, and 4 for RK4. But the formal truncation error is only for a single step. Consider instead the approximate solution arrived at by taking $n$ steps of size $h$. Its leading error term will be of order $n(h^{m+1})$ (where we have invoked the fact that the corresponding derivative in this term can be treated as constant over each step to order $h^{m+1}$). But for a fixed integration interval $\tau$ (in this case, $\tau = 15$), the number of steps $n = \tau/h$. Hence the error term becomes of order $\tau h^m$, and since $\tau$ is fixed, the overall error after $n$ steps should then go as $h^m$, not $h^{m+1}$.

2. The equations of motion for this system are derived by equating the acceleration with the negative gradient of the (specific) potential:
   
   $$\ddot{\mathbf{r}} = -\mbox{\boldmath$\nabla$}\Phi.$$
   
   
   For the $x$ component:
   
   $$\frac{d^2x}{dt^2} = - \frac{\partial \Phi}{\partial x} = -\frac{2x}{(1 + 2x^2 + 2y^2)^{3/2}} .$$
   
   
   Similarly for the $y$ component:
   
   $$\frac{d^2y}{dt^2} = - \frac{\partial \Phi}{\partial y} = -\frac{2y}{(1 + 2x^2 + 2y^2)^{3/2}} .$$
   
   

   To reduce these equations to first order, introduce $p = dx/dt$ and $q = dy/dt$, so
   

   $$\frac{dx}{dt} = p,   \frac{dp}{dt} = - \frac{2x}{(1 + 2x^2 ... ... q,   \frac{dq}{dt} = - \frac{2y}{(1 + 2x^2 + 2y^2)^{3/2}} .$$
   
   

   The code for solving this system is given at the end (integ2.c). It is analogous to integ1.c except now x[] and dxdt[] have four components each (e.g., x[] contains $x$, $p$, $y$, and $q$, in that order).

   1. The plots for these cases are given at the end. The result is a bound precessing orbit.

   2. The energy plots are given separately at the end. Note how $E$ remains bounded for LeapFrog (at least for the 3 smallest step sizes) but steadily decays for RK4. Also note that over a given interval the variation of $E$ is roughly inversely proportional to $h^2$ for LeapFrog and $h^4$ for RK4, similar to the first problem, though here we are comparing a function of the solution (namely energy), not the solution itself, so it is harder to derive the expected behavior.

3. In this case the equations are already first-order, so no transformations are required. The code (integ3.c) is given at the end, followed by the phase diagrams for each of the 4 values of $h$. Notice the spiral instability at large $h$ (the color coding is as follows: the yellow X indicates the starting point, red indicates the first third of the simulation, green the middle third, blue the final third). With the smaller $h$ values it is evident that the two populations are stable and cyclical, so long as $q < 1$ (the plots are for $q = 0$). However, for $1 \le q \lesssim 1.25$, the foxes get ``wiped out'' (population drops below $10^{-9}$) but the rabbits survive. For $q \gtrsim 1.25$, the rabbits get wiped out too (the solution goes from the X to the origin). The critical $q$ values were arrived at by trial and error, not root finding, using a stepsize of 0.1. Note that the code checks to make sure no populations go negative (which could happen due to numerical overshooting).