# Numerical integration methods

This topic is 2166 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## Recommended Posts

I think I made these codes up from the example that Gaffer had given. The Gaffer example is a great resource, but it was a bit too object-oriented for me, and I had to pare it down to understand it better. If anyone has any other integration methods that they'd like to share, I'd be very grateful, and I can put them in the same format. Maybe someone could make a sticky thread about integration, I dunno. Maybe there already is.

I'm also kind of putting this here so that I can find it easier later. Is there a way to make notes for myself on gamedev.net, that can be formatted nicely like in these forum posts?

 // Euler integration inline void proceed_Euler(body &b, const double &dt) { b.velocity += acceleration(b.position, b.velocity)*dt; b.position += b.velocity*dt; } // Runge-Kutta order 2 inline void proceed_RK2(body &b, const double &dt) { vector_3 k1_velocity = b.velocity; vector_3 k1_acceleration = acceleration(b.position, k1_velocity); vector_3 k2_velocity = b.velocity + k1_acceleration*dt*0.5; vector_3 k2_acceleration = acceleration(b.position + k1_velocity*dt*0.5, k2_velocity); b.velocity += k2_acceleration*dt; b.position += k2_velocity*dt; } // Runge-Kutta order 4 inline void proceed_RK4(body &b, const double &dt) { static const double one_sixth = 1.0/6.0; vector_3 k1_velocity = b.velocity; vector_3 k1_acceleration = acceleration(b.position, k1_velocity); vector_3 k2_velocity = b.velocity + k1_acceleration*dt*0.5; vector_3 k2_acceleration = acceleration(b.position + k1_velocity*dt*0.5, k2_velocity); vector_3 k3_velocity = b.velocity + k2_acceleration*dt*0.5; vector_3 k3_acceleration = acceleration(b.position + k2_velocity*dt*0.5, k3_velocity); vector_3 k4_velocity = b.velocity + k3_acceleration*dt; vector_3 k4_acceleration = acceleration(b.position + k3_velocity*dt, k4_velocity); b.velocity += (k1_acceleration + (k2_acceleration + k3_acceleration)*2.0 + k4_acceleration)*one_sixth*dt; b.position += (k1_velocity + (k2_velocity + k3_velocity)*2.0 + k4_velocity)*one_sixth*dt; } 

Here's some code testing it against the code that Gaffer gave out:
 // Simple RK4 integration framework // Copyright (c) 2004, Glenn Fiedler // http://www.gaffer.org/articles #include <iostream> using std::cout; using std::endl; #include <cmath> // My adaptation class vector_3 { public: double x, y, z; vector_3(const double &src_x = 0, const double &src_y = 0, const double &src_z = 0) { x = src_x; y = src_y; z = src_z; } vector_3 operator+(const vector_3 &rhs) { return vector_3(x + rhs.x, y + rhs.y, z + rhs.z); } vector_3 operator*(const double &rhs) { return vector_3(x*rhs, y*rhs, z*rhs); } vector_3& operator+=(const vector_3 &rhs) { x += rhs.x; y += rhs.y; z += rhs.z; return *this; } }; class body { public: vector_3 position, velocity; }; vector_3 acceleration(vector_3 const &pos, vector_3 const &vel) { const double k = 10; const double b = 1; return -k*pos.x - b*vel.x; } // Euler integration inline void proceed_Euler(body &b, const double &dt) { b.velocity += acceleration(b.position, b.velocity)*dt; b.position += b.velocity*dt; } // Runge-Kutta order 2 inline void proceed_RK2(body &b, const double &dt) { vector_3 k1_velocity = b.velocity; vector_3 k1_acceleration = acceleration(b.position, k1_velocity); vector_3 k2_velocity = b.velocity + k1_acceleration*dt*0.5; vector_3 k2_acceleration = acceleration(b.position + k1_velocity*dt*0.5, k2_velocity); b.velocity += k2_acceleration*dt; b.position += k2_velocity*dt; } // Runge-Kutta order 4 inline void proceed_RK4(body &b, const double &dt) { static const double one_sixth = 1.0/6.0; vector_3 k1_velocity = b.velocity; vector_3 k1_acceleration = acceleration(b.position, k1_velocity); vector_3 k2_velocity = b.velocity + k1_acceleration*dt*0.5; vector_3 k2_acceleration = acceleration(b.position + k1_velocity*dt*0.5, k2_velocity); vector_3 k3_velocity = b.velocity + k2_acceleration*dt*0.5; vector_3 k3_acceleration = acceleration(b.position + k2_velocity*dt*0.5, k3_velocity); vector_3 k4_velocity = b.velocity + k3_acceleration*dt; vector_3 k4_acceleration = acceleration(b.position + k3_velocity*dt, k4_velocity); b.velocity += (k1_acceleration + (k2_acceleration + k3_acceleration)*2.0 + k4_acceleration)*one_sixth*dt; b.position += (k1_velocity + (k2_velocity + k3_velocity)*2.0 + k4_velocity)*one_sixth*dt; } // My adaptation // Slightly modified Gaffer code struct State { double x; double v; }; struct Derivative { double dx; double dv; }; double acceleration(const State &state) { const double k = 10; const double b = 1; return - k*state.x - b*state.v; } Derivative evaluate(const State &initial) { Derivative output; output.dx = initial.v; output.dv = acceleration(initial); return output; } Derivative evaluate(const State &initial, double dt, const Derivative &d) { State state; state.x = initial.x + d.dx*dt; state.v = initial.v + d.dv*dt; Derivative output; output.dx = state.v; output.dv = acceleration(state); return output; } void integrate(State &state, double dt) { Derivative a = evaluate(state); Derivative b = evaluate(state, dt*0.5, a); Derivative c = evaluate(state, dt*0.5, b); Derivative d = evaluate(state, dt, c); const double dxdt = 1.0/6.0 * (a.dx + 2.0*(b.dx + c.dx) + d.dx); const double 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; } // Slightly modified Gaffer code struct Method { char const* name; void (*integrator)(body &b, const double &dt); body b; }; int main() { // Gaffer State state; state.x = 100; state.v = 0; // My adaptation body b; b.position.x = state.x; b.velocity.x = 0; Method methods[] = { { "RK4 ", proceed_RK4, b }, { "RK2 ", proceed_RK2, b }, { "Euler", proceed_Euler, b} }; double t = 0; double dt = 0.1; cout << "Name | ( position, velocity )"<<endl; while (fabs(state.x)>0.001 || fabs(state.v)>0.001) { integrate(state, dt); //proceed_Euler(b, dt); //proceed_RK2(b, dt); //proceed_RK4(b, dt); for(int i = 0; i < sizeof(methods)/sizeof(Method); ++i) { methods.integrator(methods.b, dt); cout << methods.name << " (" << fabs(state.x - methods.b.position.x) << ", " << fabs(state.v - methods.b.velocity.x) << ")" << endl; } t += dt; } return 0; }  Edited by Washu

##### Share on other sites
Apparently this is verlet integration:
[source lang="cpp"]void Particle::integrate(float timestep){
Vec2 Temp = P.Position;
P.Position += P.Position - P.OldPosition + P.Acceleration*Timestep*Timestep;
P.OldPosition = Temp;
}
[/source]

##### Share on other sites
Modified your code a bit to allow for comparison between all of the different integrators fairly trivially. Also fixed a few bugs while I was at it.

##### Share on other sites
Personally, I disagree with Gaffer -- Runge-Kutta is a terrible way to do integration for game simulation. Semi-implicit euler is a far better system for most games.

##### Share on other sites
So, for semi-implicit Euler update velocity first, and then use that to update position?

 inline void proceed_SemiImplicit_Euler(body &b, const double &dt) { b.velocity += acceleration(b.position, b.velocity)*dt; b.position += b.velocity*dt; } 

It seems like Verlet and velocity Verlet are designed to reverse engineer the derivatives, such as velocity and acceleration. Reduced variable count, reduced possibility of error accumulation, reduced possibility of unstable behaviour.

For the tennis, it looks like semi-implicit Euler handles things almost as well as RK4 (and obviously much faster), and a bit better than Euler. The differences are small to begin with. With drag and spin acceleration, it might be interesting to use the path as a target curved "laser" beam, where the many force-related parameters can be varied slowly and smoothly over time to create animation. Screams render-to-texture to me. Lotsa beams. Twisted beams, with a high spin RPM. Time-dependent forces would be interesting. Edited by taby

##### Share on other sites

Modified your code a bit to allow for comparison between all of the different integrators fairly trivially. Also fixed a few bugs while I was at it.

Envy.

##### Share on other sites
This is not related to the original question, but I wanted to say that passing numbers by const reference is a weird thing to do. Passing by value should be the default, and then turning that into constant references is an optimization that makes sense for classes of some size, to avoid copying. Copying native numeric types is trivial, and when you use the reference you'll end up paying the cost of a dereference, which can actually make your code slower. It also forces the calling code to put the number in memory, which again could make your code slower.

##### Share on other sites
Copying native numeric types is trivial, and when you use the reference you'll end up paying the cost of a dereference, which can actually make your code slower. It also forces the calling code to put the number in memory, which again could make your code slower.
Maybe I seriously misunderstand what you're saying or how references work on some compilers, but... what dereference overhead are you talking about?

References are aliases, not pointers. Certainly they can be implemented via pointers, and sometimes this happens, but when that is the case then it's the most efficient (and safest) way, too. References to "numbers" (and if the compiler can help it, complex objects too!) are treated by any compiler I've ever used just as if you were using that same "number" variable declared in another place, under another name.
To me, advising against references feels a bit like telling someone not to use typedefs (because the extra type indirection makes the code run slower). It's really just about the same thing, a typedef aliases a type, and a reference aliases an object.

##### Share on other sites
David Whysong's higher order symplectic integration algorithms are truly excellent. I've completely left Runge-Kutta and all the similar, large integrators for these. They come with fourth, sixth, and eighth order accuracy!

Source code as text:

There's a paper too, but I can't find a link.

When I don't need the high accuracy, I use plane old symplectic Euler:

v(n+1) = v(n) + a(n) * dt
x(n+1) = x(n) + v(n+1) * dt

I've stopped using the velocity verlet, because it handles differences in velocity very poorly.

Cheers,
Mike Edited by h4tt3n

##### Share on other sites

So, for semi-implicit Euler update velocity first, and then use that to update position?

 inline void proceed_SemiImplicit_Euler(body &b, const double &dt) { b.velocity += acceleration(b.position, b.velocity)*dt; b.position += b.velocity*dt; } 
Yeah, that should be right.

For the tennis, it looks like semi-implicit Euler handles things almost as well as RK4 (and obviously much faster), and a bit better than Euler. The differences are small to begin with. [/quote]The key is that semi-implicit Euler sacrifices some accuracy for dramatically improved stability and better energy conservation (well, bounded energy fluctuation actually). RK4 systems will tend to explode or slow down very easily because of the energy instability. This isn't noticeable for simple motions, but becomes very evident in springs for example.

This is also called sympletectic Euler at times, and there are many other classes of symplectic integrators as h4tt3n pointed out. All of them are much preferable to the Verlet or Runge-Kutta type approaches. But the semi-implicit Euler is dead simple and very good for how easy and cheap it is. Edited by Promit

• 10
• 19
• 14
• 19
• 15