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.
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.
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.