• Create Account

### #ActualWashu

Posted 14 June 2012 - 09:50 PM

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>

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

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

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


### #7taby

Posted 14 June 2012 - 06:51 PM

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>

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 &pos, vector_3 &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;
}

// 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

int main()
{
// Gaffer
State state;
state.x = 100;
state.v = 0;

body b;
b.position.x = state.x;
b.velocity.x = 0;

double t = 0;
double dt = 0.1;

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

cout << "Distance: " << fabs(state.x - b.position.x) << ' ' << fabs(state.v - b.velocity.x) << endl;

t += dt;
}

return 0;
}


### #6taby

Posted 14 June 2012 - 06:51 PM

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?

The object mercury has a position 3-vector and a velocity 3-vector. The acceleration() function returns an acceleration 3-vector, and takes in a position and velocity (useful for handling gravitation, drag, etc in one place).

// 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>

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 &pos, vector_3 &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;
}

// 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

int main()
{
// Gaffer
State state;
state.x = 100;
state.v = 0;

body b;
b.position.x = state.x;
b.velocity.x = 0;

double t = 0;
double dt = 0.1;

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

cout << "Distance: " << fabs(state.x - b.position.x) << ' ' << fabs(state.v - b.velocity.x) << endl;

t += dt;
}

return 0;
}


### #5taby

Posted 14 June 2012 - 06:50 PM

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?

The object mercury has a position 3-vector and a velocity 3-vector. The acceleration() function returns an acceleration 3-vector, and takes in a position and velocity (useful for handling gravitation, drag, etc in one place).

// 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>

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 &pos, vector_3 &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;
}

// 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

int main()
{
// Gaffer
State state;
state.x = 100;
state.v = 0;

body b;
b.position.x = state.x;
b.velocity.x = 0;

double t = 0;
double dt = 0.1;

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

cout << "Distance: " << fabs(state.x - b.position.x) << ' ' << fabs(state.v - b.velocity.x) << endl;

t += dt;
}

return 0;
}


### #4taby

Posted 14 June 2012 - 06:30 PM

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?

The object mercury has a position 3-vector and a velocity 3-vector. The acceleration() function returns an acceleration 3-vector, and takes in a position and velocity (useful for handling gravitation, drag, etc in one place).

// 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>

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 &pos, vector_3 &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;
}

// 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.5f, a);
Derivative c = evaluate(state, dt*0.5f, b);
Derivative d = evaluate(state, dt, c);

const double dxdt = 1.0f/6.0f * (a.dx + 2.0f*(b.dx + c.dx) + d.dx);
const double 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;
}
// Slightly modified Gaffer code

int main()
{
// Gaffer
State state;
state.x = 100;
state.v = 0;

body b;
b.position.x = state.x;
b.velocity.x = 0;

double t = 0;
double dt = 0.1f;

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

cout << "Distance: " << fabs(state.x - b.position.x) << ' ' << fabs(state.v - b.velocity.x) << endl;

t += dt;
}

return 0;
}


### #3taby

Posted 14 June 2012 - 06:29 PM

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?

The object mercury has a position 3-vector and a velocity 3-vector. The acceleration() function returns an acceleration 3-vector, and takes in a position and velocity (useful for handling gravitation, drag, etc in one place).

// 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>

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 &pos, vector_3 &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;
}

// Original 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.5f, a);
Derivative c = evaluate(state, dt*0.5f, b);
Derivative d = evaluate(state, dt, c);

const double dxdt = 1.0f/6.0f * (a.dx + 2.0f*(b.dx + c.dx) + d.dx);
const double 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;
}
// Original Gaffer code

int main()
{
// Gaffer
State state;
state.x = 100;
state.v = 0;

body b;
b.position.x = state.x;
b.velocity.x = 0;

double t = 0;
double dt = 0.1f;

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

cout << "Distance: " << fabs(state.x - b.position.x) << ' ' << fabs(state.v - b.velocity.x) << endl;

t += dt;
}

return 0;
}


PARTNERS