• Create Account

Numerical integration methods

Old topic!
Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.

15 replies to this topic

#1taby  Members   -  Reputation: 336

Like
3Likes
Like

Posted 14 June 2012 - 05:38 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;
}


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

#2thecoast47  Members   -  Reputation: 255

Like
2Likes
Like

Posted 14 June 2012 - 06:07 PM

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]

#3Washu  Senior Moderators   -  Reputation: 5189

Like
2Likes
Like

Posted 14 June 2012 - 09:53 PM

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.

In time the project grows, the ignorance of its devs it shows, with many a convoluted function, it plunges into deep compunction, the price of failure is high, Washu's mirth is nigh.
ScapeCode - Blog | SlimDX

#4Promit  Moderators   -  Reputation: 7190

Like
1Likes
Like

Posted 14 June 2012 - 10:05 PM

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.

#5taby  Members   -  Reputation: 336

Like
0Likes
Like

Posted 15 June 2012 - 01:21 AM

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, 15 June 2012 - 03:22 AM.

#6taby  Members   -  Reputation: 336

Like
0Likes
Like

Posted 15 June 2012 - 01:35 AM

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.

#7Álvaro  Crossbones+   -  Reputation: 13317

Like
2Likes
Like

Posted 15 June 2012 - 08:04 AM

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.

#8samoth  Crossbones+   -  Reputation: 4783

Like
1Likes
Like

Posted 15 June 2012 - 08:49 AM

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.

#9h4tt3n  Members   -  Reputation: 1051

Like
3Likes
Like

Posted 15 June 2012 - 09:36 AM

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, 15 June 2012 - 09:54 AM.

#10Promit  Moderators   -  Reputation: 7190

Like
2Likes
Like

Posted 15 June 2012 - 09:41 AM

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.

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, 15 June 2012 - 09:41 AM.

#11arkane7  Members   -  Reputation: 213

Like
1Likes
Like

Posted 15 June 2012 - 09:51 AM

I feel alvaro is correct. Anytime you are not going to alter a parameter through a reference, at least for small data types, it makes sense to pass by value.
References are actually masked pointers, they are basically the same thing just easier to use in general, so there should be some overhead of dereferencing etc.

If you are not going to alter the parameter state, and the parameter is fairly small, it is better to pass by value rather than reference. Copying an int or a double is trivial, but having to go into memory is much slower. Const reference should generally be reserved for very large class/struct objects where copying is a waste of memory.
Always improve, never quit.

#12phestermcs  Members   -  Reputation: 123

Like
1Likes
Like

Posted 15 June 2012 - 09:58 AM

The book "Game Physics" by David H. Eberly, includes free-to-use C source for several numerical integrators (and a great deal of other game related routines). The source can be found here for free http://www.geometrictools.com/LibMathematics/NumericalAnalysis/NumericalAnalysis.html.

Also, the book "Numerical Recipes in C++" has (not free to get, but free to use) source for several numerical integrators.

Both books do spend some time describing trade-offs, which basically come down to accuracy/stability vs. performance.

Edited by phestermcs, 15 June 2012 - 10:28 AM.

#13Promit  Moderators   -  Reputation: 7190

Like
1Likes
Like

Posted 15 June 2012 - 10:00 AM

Let's please not do the C references/pointers thing. This is a thread about numerical integration and I'll delete any further posts outside that.

#14arkane7  Members   -  Reputation: 213

Like
1Likes
Like

Posted 15 June 2012 - 10:13 AM

aren't there any more efficient integration methods besides sem-implicit euler?
i agree its very simple and can see why it could be used, but even the article (and people one here) claim it is inaccurate at times. Wouldn't that pose a problem, like a particle/character/whatever moving in a non-realistic way?
Always improve, never quit.

#15jjd  Crossbones+   -  Reputation: 2075

Like
1Likes
Like

Posted 15 June 2012 - 10:55 AM

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.

Sorry, just to clarify: Verlet is a symplectic integrator and what h4tt3n was talking about is 'velocity verlet' and how it handles velocity differences. Nonetheless, both are symplectic integrators, which means that under certain conditions they are very stable and accurate. However, for a lot of games those conditions are not met -- e.g. if you have friction in your physics, you do not get the same guarantee of stability because the system is not time-reversible.

Personally, I still go with semi-implicit Euler by default because I feel that it is the best 'bang-for-your-buck' integrator out there: super easy to code, debug, and use. If you are making games, then I think using something like RK4 is generally overkill.

-Josh

--www.physicaluncertainty.com
--irc.freenode.net#gdnet

#16taby  Members   -  Reputation: 336

Like
0Likes
Like

Posted 15 June 2012 - 12:19 PM

OK, so the symplectic (as close to energy conserving as possible) integrators aren't fool-proof when it comes to modeling dissipative (energy leaking) systems / forces. That kind of makes sense, if I'm seeing the big picture correctly.

And, unless your data type is *much* larger in size than the pointer type, don't bother using pass-by-reference because the time saved by skipping the copy is more than offset by the amount of time it takes to set up a reference and then dereference it.

And, you are limited to the number of posts that you can give a positive rating to per day.

And, even the cases where spin RPM is like 8000, a smaller time step with semi-Implicit Euler still does as good of a job as RK4 with a larger time step, and faster too.

Edited by taby, 15 June 2012 - 12:43 PM.

Old topic!
Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.

PARTNERS