Numerical Integration

From VDrift
Revision as of 10:51, 4 April 2009 by Venzon (talk | contribs)
Jump to: navigation, search

Note: the old version of this page is here: Old 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.

Integrator Order

The order of a numerical integrator is a convenient way to approximate the accuracy of an integrator. All integrators presented here have an order of at least one, which means that they eventually converge to the correct solution if the timestep is small enough. First order methods converge to the exact solution in a linear relation to the timestep. Higher order methods converge to the correct solution in a power relation to the timestep.

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 explicit method for numerical integration for ordinary differential equations.

 real a = system.GetAcceleration(state);
 state.x += state.v * dt;
 state.v += a * dt;

SUVAT

The "suvat" method takes its name from the variables in the equations of motion. This is a first order method that is very similar to the Euler method. The only difference is the addition of an acceleration term to the equation for the change in position.

 real a = system.GetAcceleration(state);
 state.x += state.v * dt + a*dt*dt*0.5;
 state.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 systems of undamped oscillators. Due to this property, it preserves energy better than the standard Euler method and so is often used in simulations of orbital mechanics. It is a first order integrator.

 real a = system.GetAcceleration(state);
 state.v += a * dt;
 state.x += state.v * dt;

Basic Verlet/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. It is a second order integrator. Unfortunately, this method is unsuitable for simulations where the acceleration is dependent on velocities, such as with a damper. The modified verlet variant below uses some tricks to try to get around this limitation, but this method will only provide first-order accuracy.

 if (!have_oldaccel)
     oldaccel = system.GetAcceleration(state);
 
 state.x += state.v*dt + 0.5*oldaccel*dt*dt;
 state.v += 0.5*oldaccel*dt;
 real a = system.GetAcceleration(state);
 state.v += 0.5*a*dt;
 
 oldaccel = a;
 have_oldaccel = true;

Improved Euler/Trapezoidal/Bilinear/Predictor-Corrector/Heun

The so-call "Improved Euler" method, also known as the trapezoidal or bilinear or predictor/corrector or Heun Formula method, is a second order integrator.

 STATE predictor(state);
 predictor.x += state.v * dt;
 predictor.v += system.GetAcceleration(state) * dt;
 
 STATE corrector(state);
 corrector.x += predictor.v * dt;
 corrector.v += system.GetAcceleration(predictor) * dt;
 
 state.x = (predictor.x + corrector.x)*0.5;
 state.v = (predictor.v + corrector.v)*0.5;

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 fourth-order explicit Runge Kutta algorithm. The code snippet shown below is high level, and the actual implementation is a bit more complicated. Note that the evaluate() function calls the system's acceleration function.

 DERIVATIVE a = evaluate(state, 0, DERIVATIVE(), system);
 DERIVATIVE b = evaluate(state, dt*0.5, a, system);
 DERIVATIVE c = evaluate(state, dt*0.5, b, system);
 DERIVATIVE d = evaluate(state, dt, c, system);
 const float dxdt = 1.0/6.0 * (a.dx + 2.0*(b.dx + c.dx) + d.dx);
 const float dvdt = 1.0/6.0 * (a.dv + 2.0*(b.dv + c.dv) + d.dv);
 state.x = state.x + dxdt*dt;
 state.v = state.v + dvdt*dt;

Performance Comparison

In most simulations, calculation of the system's acceleration function is computationally very expensive. Therefore, we can rank the integrators by how many times they evaluate the system's acceleration function:

Method Acceleration Evaluations
Euler 1
SUVAT 1
NSV 1
Modified Verlet 1
Improved Euler 2
RK4 4

Detailed Comparison: oscillating spring-mass

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 analytically 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 analytic 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

File:M250a1k200dt01t20.pdf-cropped.png

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.

File:M250a1k10000dt01t20.pdf.png

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

File:M250a1k10000dt01t200.pdf.png

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.

File:M250a1k100000dt01t20.pdf.png

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.

File:M250a1k100000dt01t200.pdf.png

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.

File:M250a1k1000000dt01t20.pdf.png

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

Detailed Comparison: spring-mass-damper

For realtime dynamics simulations, damping forces are usually applied. The damping force is proportional to the velocity state, while the spring force is proportional to the position state. Acceleration is calculated as:

 a = (-k*x - c*v)/m

where:

 c = 2*sqrt(k*m)

The analytic solution for the position is:

 (A + B*t)*exp(-w*t)

where:

 w = sqrt(k/m);
 B = vo + w*xo;

The constants were set to:

 m = 250.0
 A = 1.0
 xo = 1.0
 vo = 0.0
 dt = 0.1

Because the Velocity Verlet algorithm shown above isn't technically correct since due to the damper the acceleration depends on the velocity, a modified Velocity Verlet algorithm was added to the comparison which is purported to give better results for these sorts of cases:

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

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

oldaccel = a

File:Damped-m250a1k10000dt01t10.png

Cutting right to the chase, k=10,000 is where the unmodified Velocity Verlet algorithm starts to fall apart. The other integrators are similar, although note that the RK4 solution is right on top of the analytic solution.

File:Damped-m250a1k15000dt01t10.png

Increasing k to 15,000 results in the NSV starting to show major inaccuracies. Surprisingly, the Euler algorithm is still stable and doesn't show the same undesirable behavior as the NSV algorithm.

File:Damped-m250a1k18000dt01t10.png

At k = 18,000, the modified Verlet algorithm becomes unstable and we're left with only the Euler and RK4 algorithms. The Euler algorithm starts to have major error at k = 30,000, and the RK4 does the same at k = 80,000.

Summary

The results from the spring oscillator and critically damped spring are shown below as rankings for each integrator, with 1 being the best and 4 the worst. The rankings are based on when the method becomes unstable, with accuracy being used as a tie-breaker where necessary.

Method Oscillating Spring Critically Damped Spring
Euler 4 2
NSV 3 4
RK4 2 1
Velocity Verlet 1 3

It should be noted that for the oscillating spring case, the RK4 integrator resulted in poor accuracy (zero response) while all of the other algorithms except the Euler method were still giving plausible results. Although the RK4 algorithm stayed stable until the same high k values as the Velocity Verlet algorithm, the velocity verlet algorithm's better accuracy at lower k values makes it clearly better for this case.

The performance of the Euler, NSV, and Velocity Verlet integrators is similar, with the Euler and NSV methods resulting in slightly less arithmetic. The RK4 has worse performance and is more complicated. In fact, the RK4 results in 4 evaluations of the acceleration function per timestep. If this function is the performance bottleneck, then the other methods could be run with a smaller timestep (0.25*dt) with similar performance, and with a smaller timestep all methods (except the basic Euler) outperform the RK4 algorithm for both cases.

The RK4 algorithm is very stable, but it comes at a performance and accuracy cost. The RK4 integrator loses energy and requires much more arithmetic. Better stability and accuracy could be achieved at the same performance level by running the Velocity Verlet algorithm with a smaller timestep.

Based on the stability, accuracy, and performance of these integrators, I believe that the best algorithm for general realtime dynamics simulations is the modified Velocity Verlet algorithm.