Advertisement

A Comparison of Numeric Integration Schemes

Runge Kutta vs. Euler - On the Influence of time step sizes on the accuracy of numerical simulations

Calculation examples

In the previous section we learned about global error, total error any roundoff error. Now we will compare different integration schemes such as Runge-Kutta or Adams Bashforth on two sample problems. We will compare the results and draw conclusions about the accuracy of the different integration schemes.

Damped harmonic oscillation

The first thing you should do with your integration engine is checking it against a known mathematical model for verification. In order to do this you need a problem with a known analytic solution. I choose the damped harmonic oscillator as the verification model. The behaviour of such a system can be described by the following homogeneous 2nd order differential equation:

damped oscillation damped spring
Image 2: Graph of a damped oscillation. (Courtesy of Wikimedia commons [1] and [2])
\begin{equation} \nonumber\ddot {x} + {\mu \over m} \dot{x} + {k \over m} x = 0 \label{eq:eqn_of_motion2} \end{equation}

mit:

We assume the following initial conditions:

\begin{equation} \nonumber x(0)=x_0 \quad \dot{x}(0)=0 \end{equation}

The result of this differential equation is a function describing the position of the oscilator as a function of time:

\begin{equation} \nonumber x(t) = x_0 e^{-\delta t} (cos(\omega t) + \frac{\delta}{\omega} sin (\omega t)) \end{equation}

with:

\begin{equation} \nonumber \delta = {\mu \over {2m} } \quad,\quad \omega_0 = \sqrt {k \over m} \quad,\quad \omega = \sqrt{\omega_0^2 - \delta^2} \label{eq:omega} \end{equation}

The difference between analytic solution (see Image 3) and simulated solution is the total error. The relation between the time step size and the total error for different integration methods is shown in image 4.

Image 3: Comparison between the analytic and simulated solution for a time step size of h=0.005. All methods except Euler's method match the solution pretty well.
Image 4: Integration step size plotted against the square sum of total error for different integration schemes.

Obviously the theoretical relationship between the total error and the time step size is confirmed. You can clearly see there is an optimal time step size for which the total error is minimal. You can also see that the total error is indeed smaller for the high order schemes with RK5 clearly showing the lowest error of all integration methods. Something you will also notice is that the higher the order of the integration scheme the larger is the step size that produces the best result. If you look at the Euler scheme something else is noteable: Despite being able to use step sizes smaller than in any other method the total error is still significantly larger in comparison to any other scheme, even if they use time step sizes of a order of magnitude larger. And if you look at Image 3 you will see that Euler just doesn't get it right even with a moderate step size. So what do we learn from that? How about that:

Lets have a closer look at the timing of the simulation cause after all what one really wants is getting an accurate result without spending too much time. So lets look at image 5 which illustrates the dependency between time step size and calculation time. The basic (and obvious) rule for all integration schemes is: The larger the timestep the faster the calculation. If you look closely at image 5 you will see a perfect set of straight lines, except for well they arn't really straight because i had to apply logarithmic scaling to both axis in order to get a meaningful diagram. So it's no linear relation between time step size and calculation time its a power law. If you cut the time step size in half you double the number of points in the calculation. The bottom line is that you really pay for small timesteps regardless which integration method you use! So you'd better choose wise...

So let's choose wisely and analyse yet another diagram (image 6). This time we see the calculation time vs. the global error squared. This is the point were it's getting interesting because now the relation between accuracy and the calculation time becomes obvious. So how long does any of the schemes need to produce a result with a given accuracy? The answer is right here in image 6:

Image 5: Integration step size vs. calculation time. Unsuprisingly the calculation time decreases over the the stepsize.
Image 6: Calculation time vs. sum of error squares (global error).

So what do we make of it? A table! At least thats what i make of it. Lets check how long it takes to solve our differential equation and lets say we do not want the square sum of the error to exceed 1e-5. Which scheme does it quickest?

Rank Integration Scheme Calculation time
in ms
Time step width
1 RK 5 0.01 0.02
2 RK 4 0.015 0.01
3 ADB 5
(* see below)
0.03 0.005
4 ADB 4 0.039 0.0035
5 ADB 3 0.09 0.0015
6 Euler mod.
Heun
0.5 0.0003
7 RK3
ADB2
0.7 0.0003
9 Euler 36812 4e-9
Table 2: Integration scheme ranking for an accuracy of 1e-5.

The Winner is the Runge Kutte 5th order scheme! In this test it performs better than any other scheme because it allows using larger timesteps without loosing too much accuracy. A strong second goes to the RK4 scheme which is still better than the 5th order Adams-Bashforth method. ADB 5 is a bit of a disappointment here cause as a multistep method it had all to win this race. The Calculations needed for each timesteps are minimal and yet it came in at the third place. The reason for this result is very simple: The calculations needed for calculating the forces for a damped harmonic oszillator are too simple. Multistep schemes can't play out their advantage over the single Step schemes if there is almost nothing to calculate. This results would change when dealing with more complex calculations such as 3-body problems or an n-body problem. So don't write multistep schemes off based on these results! Also disappointing is the performance of the RK3 scheme which fell behind the modified Euler and Heun methods and came in tied with ADB2. The reason is simple: It needs more calculations but fails to allow much larger timesteps.

On the next page we will compare the performance of different integration schemes on the magnetic pendulum simulation.

back      next

You might also like: