# Difference between revisions of "Numerical Integration"

Numerical integration constitutes a broad family of algorithms for calculating the numerical value of a definite integral. This the backbone of physics simulations because it allows calculation of velocity and position from forces (and therefore acceleration) applied to a rigid body.

## Criteria

In realtime simulations, the most important integrator criteria are performance, stability, and accuracy. Performance refers to how long it takes to perform the integration for a given timestep. Stability refers to how well the integrator copes with stiff constraints such as high spring constants before errors becomes unacceptably large. Accuracy refers to how well the integrator matches the expected result, and includes whether or not the integrator damps out (loses energy) over time.

## Euler Integration

The Euler method is a first order numerical procedure for solving ordinary differential equations (ODEs) with a given initial value. It is the most basic kind of explicit method for numerical integration for ordinary differential equations. The performance is excellent, but both stability and accuracy are poor. Because it is so simple, it unfortunately gets used very often. Vamos uses Euler integration.

``` a = acceleration(state, t+dt)
x += v*dt
v += a*dt
```

## Newton-Stormer-Verlet (NSV) / Symplectic Euler / Euler–Cromer algorithm

the Euler–Cromer algorithm or symplectic Euler method or Newton-Stormer-Verlet (NSV) method is a modification of the Euler method for solving Hamilton's equations, a system of ordinary differential equations that arises in classical mechanics. It is a symplectic integrator, which is a class of geometric integrators that is especially good at simulations of dynamics and hence it yields much better results than the standard Euler method. The performance is excellent, and stability is fair, and accuracy is excellent.

``` a = acceleration(state, t+dt)
v += a*dt
x += v*dt
```

## Velocity Verlet

Verlet integration is a numerical integration method originally designed for calculating the trajectories of particles in molecular dynamics simulations. The velocity verlet variant directly calculates velocity. The performance is fair, and both stability and accuracy are excellent. Unfortunately, calculating the velocity depends on knowing the acceleration for the current iteration, which poses a problem when the acceleration depends on the velocity. Using the velocity from the last iteration to calculate the acceleration may have implications on the accuracy of the method.

``` if (not oldaccel)
oldaccel = acceleration(state, t+dt)

x += v*dt + 0.5*oldaccel*dt*dt
a = acceleration(state, t+dt)
v += 0.5*(a + oldaccel)*dt

oldaccel = a;
```

## Runge Kutta 4

The Runge–Kutta methods are an important family of implicit and explicit iterative methods for the approximation of solutions of ordinary differential equations. The Runge Kutta 4 (or RK4) is a well-known 4th-order explicit Runge Kutta algorithm. The code snippet shown below is high level, and the actual implementation is a bit more complicated. Performance is poor, since the acceleration must be evaluated for all objects 4 times per iteration, stability is excellent, and accuracy is fair.

``` Derivative a = evaluate(state, t)
Derivative b = evaluate(state, t, dt*0.5f, a)
Derivative c = evaluate(state, t, dt*0.5f, b)
Derivative d = evaluate(state, t, dt, c)

const float dxdt = 1.0f/6.0f * (a.dx + 2.0f*(b.dx + c.dx) + d.dx)
const float dvdt = 1.0f/6.0f * (a.dv + 2.0f*(b.dv + c.dv) + d.dv)

state.x = state.x + dxdt*dt
state.v = state.v + dvdt*dt
```

## Detailed Comparison

For comparing these algorithms I used a simple spring-mass oscillator, because it can be difficult to integrate when the spring is very stiff, but it can be exactly solved easily so I have something to compare the integrators to. In addition, its force calculation depends only on position, which allows the Velocity Verlet algorithm to work as it is commonly used. For this simulation the instantaneous acceleration input into all integrators is calculated as:

``` a = -k*x/m;
```

where k is the spring constant, x is the position, and m is the mass.

The exact solution is calculated as:

``` A * cos (sqrt(k/m)*t)
```

where A is the amplitude (and the initial position) and t is the time in seconds.

The constants were set to:

``` A = 0.5
m = 250.0
dt = 0.1
```

This is the first simulation, with k set to 200 and run for 20 seconds. The Euler integrator is already unstable, with quickly increasing error as time goes on. All of the other methods are similar for this simulation.

The k constant has been increased to 10,000. At this value all of the non-Euler methods are initially similar, but....

This is the same k constant of 10,000 after ~200 seconds. The RK4 integrator is losing energy, while the NSV and Velocity Verlet methods have amplitudes similar to the exact answer.

The k constant has been increased to 100,000. The NSV integrator is unstable at this level. The RK4 integrator is almost uniformly zero. The only integrator that is still close to the exact value is the Velocity Verlet integrator.

This is the same k constant of 100,000 after ~200 seconds. The Velocity Verlet integrator is doing pretty well here, mostly preserving energy. The RK4 integrator is zero.

The k constant is now 1,000,000. At this level both the RK4 and Velocity Verlet integrators quickly become unstable.

## Summary

Each method is ranked by the three criteria as Excellent, Fair, or Poor:

Method Stability Accuracy Performance
Euler Poor Poor Excellent
NSV Fair Excellent Excellent
RK4 Excellent Fair Poor
Velocity Verlet Excellent Excellent Fair

Both the NSV and Velocity Verlet algorithms are good choices for realtime dynamics simulations. For the Velocity Verlet algorithm to perform as well as it did in this example, the acceleration must be derived from the position only and not the velocity. If force calculations require velocity (such as a damper), then one must make sure one understands the implications of using the last velocity to calculate forces, as it may degrade the algorithm's accuracy.