ASTR615 Fall 2015 Problem Set #5

Due Wed November 18th, 2015

Write a program to integrate any number of coupled differential equations using the Euler method, fourth-order Runge-Kutta, and Leapfrog (note: Leapfrog only applies to special cases). You will be using this program in a future assignment, so make sure it's well documented. It's recommended that you use double precision throughout.

1. Use your program to solve the following differential equation for $x(t)$:
   
   $$\frac{d^2x}{dt^2} + x = 0,$$
   
   
   with initial conditions $x(0) = 0$, $\dot{x}(0) = 1$. Note the analytical solution is $x = \sin(t)$.
   1. Integrate the equation for $0 \le t \le 15$ using each of the methods, and step sizes of 1, 0.3, 0.1, 0.03, and 0.01.
   2. Plot your integration results against the analytical solution for each case. (Hint: do all the Euler plots on one page, with one plot per timestep; then all the Leapfrog plots on another page, etc.) Comment on the results.
   3. Plot $\log\vert x_\mathrm{numerical}(15) - x_\mathrm{exact}(15)\vert$ as a function of $\log(\mathit{stepsize})$ in each case and comment. (Hint: does the error have the expected dependence on the stepsize? Remember you're integrating over many steps, not just one.)

2. Now try the two-dimensional orbit described by the potential:
   
   $$\Phi = - \frac{1}{\sqrt{1 + 2x^2 + 2y^2}} ,$$
   
   
   where we are assuming unit mass for the particle in this potential. Show analytically that the orbits are given by the coupled differential equations:
   
   $$\frac{d^2x}{dt^2} = - \frac{2x}{(1 + 2x^2 + 2y^2)^{3/2}}$$
   
   
   
   
   $$\frac{d^2y}{dt^2} = - \frac{2y}{(1 + 2x^2 + 2y^2)^{3/2}}$$
   
   
   and then reduce these to 4 coupled first-order equations.
   1. Integrate this system for $0 \le t \le 100$ with the initial conditions $x = 1$, $y = 0$, $\dot{x} = 0$, $\dot{y} = 0.1$. Try this with Leapfrog and Runge-Kutta, and step sizes of 1, 0.5, 0.25, and 0.1. Plot $x$ vs. $y$ for these integrations.
   2. Plot the energy $E = (\dot{x}^2 + \dot{y}^2)/2 + \Phi(x,y)$ as a function of time for your integrations.

3. Plot phase diagrams ($x$ vs. $y$) for the Lotka-Volterra Predator-Prey model:

   $\dot{x} =$	$Ax - Bxy - dx$
   $\dot{y} =$	$-Cy + Dxy - ey$

   where $x$ is the prey density (rabbits), $y$ the predator density (foxes), $A = 1$ (rabbit reproduction rate), $B = 0.1$ (rabbit consumption rate by foxes), $C = 1.5$ (fox death rate by natural causes), $D = 0.03$ (fox population growth rate due to consumption of rabbits), and $d = e = 0$ (hunting rate of foxes and rabbits, respectively). Use only your Runge-Kutta integrator, with $t$ = 0 to 100 and timesteps of 1, 0.5, 0.25, and 0.1 to solve this system, starting with $x = 30$ and $y = 3$. If $d = e = q$, for roughly what value of $q$ do both populations drop below $10^{-9}$ by $t = 100$ for a timestep of 0.1?