Skip to main content
gamedev.net
🔒 Locked

Integration function implementation

Started by _V_ Jan 5, 2007 at 1:48 AM 4 replies 900+ views
Original Post
_V_
_V_
Hi, I am after some help/advice on how to write a reasonably reusable integration function. The calculations for Euler integration are fairly simple: y(t + dt) = y(t) + dt * f(t, y(t)), where dt is the step size and f() is the derivative evaluated at t, y(t). To move a particle, the following steps are performed: position = position + velocity * dt velocity = velocity + acceleration * dt What I have so far is a function euler() which takes the following parameters: t, dt, y(t) and a function pointer to the function f(t, y(t)) which evaluates the slope at that point. 

euler(t, dt, y(t), f(t, y(t)))
{
    return y(t) + f(t, y(t)) * dt;
}

position = euler(t, dt, position, velocity());
velocity = euler(t, dt, velocity, acceleration());
In this simple case lets say we have a particle with initial a=10, v=0, p=0. The function acceleration simply returns the constant 10. Function velocity() return the current velocity. This works fine. But let's say we make things slightly more complex: 

acceleration(t, y(t))
{
    return -0.5 * y(t);
}
Now acceleration is no longer constant. My question is, what is the correct way to formulate the velocity function so that the proper velocity is returned? I would like the same function interface to work for other integrators so that when I decide to implement RK4 at some point I will only need to replace all instances of euler() with rk4(). A lot of examples that I've come across have all of the variables integration crammed into the one function. Thanks for your attention. [Edited by - _V_ on January 5, 2007 4:53:23 AM]
nmi
nmi
Your numerical integrator needs the first derivate (function f) as function pointer, so it can call it multiple times. For instance, an RK4 could be implemented as:
template <typename T>T rk4(double t, double dt, T y, T (*f)(double, T)){  T k1 = f(t, y);  T k2 = f(t + dt/2, y + k1 * dt/2);  T k3 = f(t + dt/2, y + k2 * dt/2);  T k4 = f(t + dt, y + k3 * dt);  return y + (k1 + 2*k2 + 2*k3 + k4) * dt/6;}

Making the function a template allows to use any datatype (for instance a mathematical vector) for the calculation.

_V_
_V_
Yes, I did mention that I am passing f as a function pointer. Perhaps I should have written it using proper C/C++ code to make it clear, rather than the pseudocode I have above. I have actually used templates in my actual code, but left them out of the problem description as I didn't feel they were relevant there.


Thanks for the suggestion.
nmi
nmi
The acceleration of an object is a second derivative. You need to map it to a first derivative . This can be accomplished by packing the pair (position, velocity) into a state vector, then letting f calculate the first derivative of each resulting in (velocity, acceleration).

The formulate for the euler integrator then would look like this:
(newpos, newvel) = (pos, vel) + dt * (vel, acc)

So this would result in code like this:
// input is current time t and current state y// output is first derivative of y at time tvec2d f(double t, vec2d y){  vec2d result;  result[0] = y[1]; // copy velocity value  result[1] = acceleration; // constant acceleration in this example  return result;}


The slope function f maps (pos, vel) to (vel, acc) to produce the first derivative. The velocity can simply be copied. The acceleration, if a constant, can be set to that value.
Christian Weis
Christian Weis
I don't get it!
Maybe it's because of my poor math skills, but what's so wrong about for not just using your current velocity for your position integration?!
Your velocity is integrated by your non-constant acceleration function. So all looks fine for me!

cu,
Chris
_V_
_V_
Quote:
Original post by nmi
The acceleration of an object is a second derivative. You need to map it to a first derivative . This can be accomplished by packing the pair (position, velocity) into a state vector, then letting f calculate the first derivative of each resulting in (velocity, acceleration).

The formulate for the euler integrator then would look like this:
(newpos, newvel) = (pos, vel) + dt * (vel, acc)

So this would result in code like this:
// input is current time t and current state y// output is first derivative of y at time tvec2d f(double t, vec2d y){  vec2d result;  result[0] = y[1]; // copy velocity value  result[1] = acceleration; // constant acceleration in this example  return result;}


The slope function f maps (pos, vel) to (vel, acc) to produce the first derivative. The velocity can simply be copied. The acceleration, if a constant, can be set to that value.


That will most likely work, but it seems like a lot of effort to manage the state structure. Is there maybe another way this could be handled?

Topic Locked

This topic has been locked by a moderator. New replies are not allowed.

Sign in to reply to this topic.