**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[i].integrator(methods[i].b, dt); cout << methods[i].name << " (" << fabs(state.x - methods[i].b.position.x) << ", " << fabs(state.v - methods[i].b.velocity.x) << ")" << endl; } t += dt; } return 0; }

**Edited by Washu, 14 June 2012 - 09:50 PM.**