Computational Physics Lab 2

Exploring ODE Algorithms

Whenever you model a dynamic system numerically, two very important questions invariably come up:

Which ODE algorithm to use?
How to pick the time step?

The purpose of this lab is to develop some intuition and gather experience needed to answer these questions.

Part 1. Implementing an ODE Solver

In this part of the lab you will implement and test an ODE solving algorithm. The 2^nd order Runge-Kutta represents a good case to study.

Download the following python modules:

v3.py -- A 3d vector class
cpode.py -- Includes a base class for ODE solvers and a few algorithms (Euler, RK4, etc.)
cpforces.py -- A vector model for physical forces

Note the new implementation of the "EulerSolver" class in the cpode module. The features specific to the Euler method were extracted into a small class which now inherits from a much larger "OdeSolver" class (compare with the old code). Classes which implement other ODE solving algorithms (e.g., Euler-Cromer) also inherit from "OdeSolver" and modify just the part of the system specific to those algorithms. This kind of design conforms to the Don't Repeat Yourself (DRY) principle and allows for very efficient implementation of new ODE solvers.

The following examples illustrate the usage of the above modules. You can download them, make them executable, and run.

golf_ball.py -- Code in this file simulates the motion of the golf ball using the air drug and Magnus force models discussed by Giordano and Nakanishi in section 2.5. The ball trajectories generated by this script should be similar to the "normal drive" curve in Figure 2.13 of the textbook.

oscillations.py -- This example simulates the behavior of a simple harmonic oscillator using several different ODE solvers.

Write a python module which implements the 2^nd order Runge-Kutta algorithm. For a dynamic ODE in the standard form

the 2^nd order Runge-Kutta stepping scheme looks like this:

where

In your module you should implement a class which inherits from "cpode.OdeSolver", just like the "EulerSolver" and other classes inside "cpode.py". Modify the "oscillations.py" example so that you can see the curve produced by your algorithm when you know exactly what to expect. Modify the "golf_ball.py" example and use your algorithm instead of "EulerSolver" in order to make sure that your code works correctly with both scalar and vector forces.

Part 2. Exploring ODE Algorithm Performance

The purpose of this part of the lab is to study the precision of several ODE solving algorithms. We are going to model frictionless oscillatory motion. For this type of motion, conservation of energy provides an extremely usefull crosscheck of the algorithm performance. Proceed as follows:

1) Write a class which implements a 1-dimensional force model which corresponds to the following potential:

Your class should inherit from "cpforces.BasicForce", just like various classes in the cpforces module. In your calculations, don't forget that the derivative of |x| is sgn(x).

2) Replace the "RestoringForce" in the "oscillations.py" example with your model. Check that you get the same behavior as the original in case p = 2. See what happens for the case p = 8. (Note: the Euler method might behave very poorly or even diverge to infinity with the default time step values. Do you know why? Reduce Δt until it works reasonably well.) Does the oscillation period depend on the amplitude?

3) Simulate the motion for 10 periods for the cases p = 2 and p = 8 (you can determine the period for the p = 8 case by eyeballing the plots). For the last oscillation period, calculate the average relative energy error. For each time step i, the relative energy error is

and you need to average it over all steps in the period. The solver keeps the whole history of x and v values, you just need to use the last 1/10^th of it. The energy is conserved, so that E_exact can be calculated from the initial conditions.

You might want to plot e_i as a function of time. This will help you visualize the error trend and get a detailed impression about the algorithm performance.

4) Make a function which would perform average relative energy error calculation automatically for an arbitrary ODE solving algorithm and an arbitrary time step. Plot the relative energy precision as the function of the time step for several algorithms: Euler, Euler-Cromer, 2^nd order Runge-Kutta, 4^th order Runge-Kutta. Use logarithmic scales for both the time step and the precision (you can use the loglog function to do this) -- then the plot will tell you directly the number of significant digits in your results. Useful Δt range is from 1 to 10^-5 or so, for smaller time steps Python solver becomes too slow. What is the empirical order of convergence of these algorithms for p = 2 and p = 8? At what Δt value the precision of the 4^th order Runge-Kutta algorithm is dominated by the round-off errors?

Please submit your lab report before 2 pm 01/29/08. In your report, include the script which generates the plots from Part 2 item 4) and place a small text file named "lab2_report.txt" with your conclusions about convergence rates in the "doc" directory.