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;
}