Homework #1: ODE's & Particle Dynamics

This assignment uses solvers for ordinary differential equations (ODE's) to simulate the motion of particles. Relatively simple equations give rise to complex and realistic motion. At the most basic level, the ODEs are solved by using linear extrapolation of the state variable using the derivates of the state at discrete time steps. The accuracy and performance of the solver is tied directly to the size of the time step.

The class ODESolver provides a generic interface for the different solvers used in these programs. There are three subclasses of ODESolver, EulerSolver, MidpointSolver, and RK4Solver corresponding to the euler, midpoint, and Runge-Kutta four point methods. Of these the Runge-Kutta is the most accurate, but it is also the most expensive. In general, Runge-Kutta can take a much larger time-step with the same accuracy as the other methods.

ODESystem is the generic interface used by the solver. It includes methods for packing and unpacking the state of a system to and from simple vectors upon which the ODESolver operates. ODESystem also provides a method for calculating derivatives. For this project, there are two ODESystem, CannonBall and Spring. The state maintained by both of these classes is position and velocity. The only difference between them is the derivative calculations.


For the cannon program, a particle is shot from the origin with a muzzle velocity determined by the amount of gunpowder used. The particle is subject to the forces of gravity and viscous drag. Without any drag, the path of the particle is a symmetric arc. When drag is increased the end of the arc becomes foreshortened. If the drag force is very large with respect to the weight of the particle, the path resembles a right triangle. The particle quickly loses all velocity and comes to a stop in mid-air. It then falls very slowly straight down. The following graph shows the trajectory with drag = 50 kg/s and various masses.

This graph shows how the error varies with the time step. With a small time step all of the solvers give the nearly the same solution. As the time step becomes larger however, the Euler and midpoint methods begin to diverge rapidly. With very large time steps Runge-Kutta also diverges.

To compare the efficiency of mid-point and Runge-Kutta to Euler the method I plotted the relative improvement in relative error per the amount of work. Thus if mid-point had an error of .001 and Euler had an error .01 then there would be a relative improvement of 10, but mid-point does twice the work so it is only 5 times more efficient. This graph is weird. For some reason the mid-point efficiency is greater than Runge-Kutta for the first two time steps. Thereafter, Runge-Kutta is always more efficient (though it is hard to see on the compressed scale of this graph.

Stability is more of an issue with the spring-mass system. However this system can demonstrate instability when the drag constant is very high. If the time-step is too large then the particle will change direction within a few steps out of the barrel of the cannon. With the drag constant = 1000 kg/s I saw this happen at .02s+ for Euler and Midpoint and .0279s+ for RK4. The difference is due to the weighting of the sub-steps used in RK4.


To test the stability of the different integration methods on the spring system I set the spring constant very high to 1000 kg/s^2 with a mass of 5 kg and just a little bit of damping at .01 kg/s. I then let the simulation run for 10 seconds. If the time step is sufficiently small the mass should come nearly to rest. If it is too large the mass oscillates with increasingly large amplitude. For no observable increase in amplitude within 10 seconds I obtained the following maximum time steps.
Method Time step
Euler .0001
Midpoint .01
RK4 .2

RK4 is MUCH more stable than the other two methods.

Here are three plots of the simulation using three different solvers at three different time steps. The time steps are .01, .1, and 1. The physical parameters for the system are m = 5 kg, ks = 15 kg/s^2, kd = 1 kg/s, rest length = 5 m. By the time we get to the large time step of 1s, only Runge-Kutta has not exploded.


Here is a screen shot of the cannon program.

Keys for cannon program:
R Reset simulation
V Reset view
Z Auto-zoom (keeps the projectile in the view)
Space Start/Stop

This is a screeshot of the spring program.

This is actually a 2D spring. Use the left button to place the mass. Use the right button to place the anchor. Great fun can be had by dragging around the anchor while the simulation is running.

Source code and executables (145K)
Qt and GLUT DLLs (1.5M)