Numerical integration methods

Shawn Halayka · 2012-06-15T18:19:13

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

Math and Physics Programming

Started by taby June 14, 2012 11:38 PM

13 comments, last by taby 11 years, 10 months ago

arkane7

213

June 15, 2012 03:51 PM

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.

phestermcs

123

June 15, 2012 03:58 PM

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.

Promit

13,404

June 15, 2012 04:00 PM

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.

SlimDX | Ventspace Blog | Twitter | Diverse teams make better games. I am currently hiring capable C++ engine developers in Baltimore, MD.

arkane7

213

June 15, 2012 04:13 PM

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.

jjd

2,140

June 15, 2012 04:55 PM

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
--linkedin
--irc.freenode.net#gdnet

taby

1,504

Author

June 15, 2012 06: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.

Numerical integration methods

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Numerical integration methods

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Reticulating splines