Scientific Computing, Spring 1999

Assignment 6.

Given March 23, due April 13.

Objective: To explore numerical methods of ODE’s.

We want to compute the trajectory of a comet. In nondimensionalized variables, the equations of motion are given by the inverse square law:

\frac{d^{2}}{dt^{2}} ( \begin{matrix}x \\ y \end{matrix} ) = \frac{−1}{(x^{2} + y^{2})^{3/2}} ( \begin{matrix}x \\ y \end{matrix} ) .

We will always suppose that the comet starts at t = 0 with x = 10, y = 0, ẋ = 0, and ẏ = v_{0}. If v_{0} is less than the escape velocity (1/\sqrt{5} here), the comet will move in an elliptical orbit, returning to x = 10, y = 0 after a time T (v_{0}). This is the time for a single orbit. For larger values of v_{0}, the trajectory will not be a closed orbit, but will move farther and farther from the sun as time increases.

(1) Compute the position and velocity at time t = 30 with v_{0} = .2 using a fixed Δt and the three methods: Forward Euler, second order Adams Bashforth, and four stage fourth order Runge Kutta. Do a convergence study to verify the expected order of accuracy of each of these methods. Compare the number of time steps needed to achieve 1% accuracy in each case.

(2) Use a fixed time step to compute T (v_{0}) as a function of v_{0} for v_{0} in the range [.01, .25]. Plot the result. For this, you will need some way to pick an appropriate time step, possibly by doing a convergence study. You will also have to safeguard your program against the possibility that T (v_{0}) is infinite^{1}.

(3) Write an adaptive time step program that uses the fourth order Runge Kutta together with a local truncation error strategy such as the one discussed in class. The adaptive program should much more efficient than the fixed step method when v_{0} is small.

1“This is not a drill.”