Numerical Integration

The term "Numerical Integration" seems to encompass quite a few different but related topics. But for the purposes of my project, we are refering to a family of algorithms for approximately solving ordinary differential equations. Many different methods exist, such as Euler, Verlet and a whole family of methods called Runge-Kutta.

Euler Integration is the most basic method. It involves multiplying a differential by a change (or delta) in the independent variable to work out the change in the dependant variable over a "step" or interval. This value is then added to the old value from the last step.

E.g. At time 0.1 seconds, a particle was moving at 10m/s and the acceleration of the particle was 5m/s^2. How fast is the particle moving at 0.2 seconds?

The step measurement here is in time so let's say the distance between each step is 0.1 seconds.

And differential here is acceleration (the change in velocity with respect to time):

A(t) = d V(t) / d t

We want to know how much the velocity has changed this step. So it's intuitive that we multiply this differential by the size of the time step (delta t):

Change in V(t) = A(t) * dt

The new V(t), then, is this change plus the old value from the last step:

V(t) = V(old_time) + A(t) * dt

Punching in the numbers from above:-

V(t) = 10 + (5 * 0.1)
V(t) = 10.5

...we get the velocity to be 10.5m/s after 0.2 seconds.

This is basic Euler Integration. It is a fast and simple method. It's major drawback is the accuracy of the results it produces. All Numerical Integration techniques are designed to give an approximate solution but the Euler method is the most inaccurate of all. It's huge downfall is that it assumes the differential to be constant over the step interval. What happens in the example above when the acceleration wildly fluctuates between 0.1 seconds and 0.2 seconds? The result at 0.2 seconds could be hugely inaccurate. What's more is the error accumulates since each result is based on the value from the last step.

-----------------------------------------------

Verlet Integration is a scheme often used in computer simulations. It evaluates the differential, not as a result of integrating, but as the difference in the dependant variable since the last step.

For example:

Si+1 = Si + (Si - Si-1) + a* dt*dt

Where Si+1 is the new position, Si is the last position, a is acceleration, dt is the size of the time step and (Si - Si-1) is velocity.

In this example, not evaluating velocity by integrating acceleration means greater accuracy and more predictable results. It breeds avantages such as being able to instantaneoulsy place objects, with their velocity implicitly defined as a result of their movement. However, as with Euler, the Verlet scheme also assumes the differential (acceleration in the example) to be constant over the time step.

-------------------------------------------------------

The family of Runge-Kutta methods evaluate the derivative at multiple points along a time step interval to converge to a more accurate solution.

One known as the Midpoint Method (or second-order Runge-Kutta) uses the initial derivative to find the derivative half-way between this step and the next. Then this "half-way derivative" is used to integrate over the whole step.

These methods are superior to Euler and Verlet in that they accept changes in the derivative over the time step.

A Higher-order method, such as the Runge-Kutta 4 method, calculates four derivatives between the current step and the next. Each successive calculation of the derivative uses the old derivative as input. An example of how the four derivatives are calculated follows.

Suppose that velocity (derivative of position wrt time) is given by a known function of time and position:

dS/dt = f (t, s)

Derivatives per time step are calculated as follows:

• a = f (t0, position from last interval)
• b = f (t0 + h/2, position half-way to next interval (calculated from derivative a))
• c = f (t0 + h/2, position half-way to next interval (calculated from derivative b))
• d = f (t0 + h, position at next interval (calculated from derivative c))

Where t0 is time at the beginning of the last interval and h is step-size.

Then the Runge-Kutta 4 method evaulates the final result as follows:

[new value] = [value from last interval] + h/6(a + 2b + 2c + d)

---------------------------------------------------------------------------------

Within the field of numerical integration is a concept of Adaptive Time-Stepping. Adaptive Time Stepping is designed to improve the accuracy of the solution whilst maintaining as much efficiency as possible. It works on the principle that with smaller time-steps, results become more accurate. This fact makes all these algorithms adaptive in terms of accuracy, simply by modifying the time-step. This is an advantage but the problem is, decreasing the time step means increasing the computation cost. So adaptive time stepping works on the idea that the ideal would be to take small time-steps when the results coming back are inaccurate and larger time steps when they are of an acceptable accuracy.

Adaptive time stepping adapts the integration scheme through a feedback loop, which per step calculates the error in the last result, makes a decision whether this is "acceptable" and then modifies the time-step accordingly.

Sources:

Numerical Recipes In C : The Art of Scientific Computing, by W.Press, B.Flannery, S.Teukolsky and W.Vetterling.
http://gafferongames.wordpress.com/game-physics/integration-basics/
http://www.gamedev.net/reference/programming/features/verlet/default.asp
http://www.gamedev.net/reference/programming/features/verlet/default.asp
http://mathworld.wolfram.com/Runge-KuttaMethod.html
http://www.myphysicslab.com/runge_kutta.html
http://numericalmethods.eng.usf.edu/mcquizzes/08ode/runge4th.pdf