External forces and RK4

Recommended Posts

JimmyDeemo    156
So i started this [url="http://www.gamedev.net/topic/597235-getting-my-head-around-rk4/"]thread[/url] a while ago, and the replies dried up. I don't think i made it clear as to my problem so i wanted to re word it so that i might get some more responses.

Firstly let me point out that i know RK4 is over kill, but it's mainly just so that i can learn and practise my coding techniques, so while i understand there are good intentions to pointing out other systems lets try and keep this about RK4. So i have implemented a system similar to that detailed on [url="http://gafferongames.com/"]http://gafferongames.com[/url] but i am still confused as to where external forces fit into the system.

From my understanding the way RK4 works is that it needs to calculate the force of a derivative implicitly from the current state, rather than explicitly applying forces once per update. So if that's the case how do i say apply some thrust to a ship for example?

Let me break it down furthur i have the position and momentum of the RigidBody prior to any calculation, i also have the velocity of the current derivative step. What i want to be able to do at this point is apply things like collision input and 'push' forces, and the end result should be to calculate the forces for this current derivative step. I also have the current time step also.

Here is some code for you, hopefully its clear enough, please just ask if anything is confusing.

[code]

/**
* Integrate a state forward by a delta time.
*
* This function uses the RK4 integration method to update the state by a given
* time.
*
* @param state Reference to the state that needs integrating.
* @param deltaTime Time to move the state on by.
*/
void CPhysController::Integrate(CRidgidBody::CRBState &state, float deltaTime)
{
//RK4 uses 4 derivatives of the current state at different points within the delttime.
//Each evaluation uses the previous derivative to calculate next.
CRidgidBody::CRBDerivative a = Evaluate(state, deltaTime, CRidgidBody::CRBDerivative() );
CRidgidBody::CRBDerivative b = Evaluate(state, deltaTime*0.5f, a);
CRidgidBody::CRBDerivative c = Evaluate(state, deltaTime*0.5f, b);
CRidgidBody::CRBDerivative d = Evaluate(state, deltaTime, c);

//Using a wighted sum that comes from the Taylor Series expansion, the best overall derivative
//is calculated. This can then increment the state.
state.UpdatePosition( 1.0f/6.0f * deltaTime * (a.mVelocity + 2.0f*(b.mVelocity + c.mVelocity) + d.mVelocity) );
state.UpdateMomentum( 1.0f/6.0f * deltaTime * (a.mForce + 2.0f*(b.mForce + c.mForce) + d.mForce) );
}

CRidgidBody::CRBDerivative CPhysController::Evaluate(CRidgidBody::CRBState &state, float deltaTime,
const CRidgidBody::CRBDerivative &prevDerivative)
{
//First move the state on base on this time step value, using euler integration.
state.UpdatePosition( prevDerivative.mVelocity * deltaTime );
state.UpdateMomentum( prevDerivative.mForce * deltaTime );

//Prep our derivative
CRidgidBody::CRBDerivative derivative = CRidgidBody::CRBDerivative();
derivative.mVelocity = state.GetVelocity();
Resolve(state, derivative);
return derivative;
}

void CPhysController::Resolve(CRidgidBody::CRBState &currentState, CRidgidBody::CRBDerivative &derivative)
{
//???
}
[/code]

Any help understanding this would be great. I did think that perhaps i should be coming at the problem from a different angle, perhaps someone here could enlighten me.

Share on other sites
johnstanp    267
I haven't taken a lookt at your code: it is not necessary, since I'll give a general answer.
You can consider that the total force applied to the rigid body is constant between t and t+dt: this is what I do for my update function using RK4.

Share on other sites
johnstanp    267
An example of how you could implement it:

[source lang="cpp"]
#include <iostream>

struct Vector3;
struct State;
struct RigidBody;

struct Vector3
{
float x, y, z;

Vector3();
Vector3( float x, float y, float z );

Vector3 & operator+=( Vector3 const & );
Vector3 operator+( Vector3 const & )const;
Vector3 & operator-=( Vector3 const & );
Vector3 operator-( Vector3 const & )const;
Vector3 & operator*=( float );
Vector3 operator*( float )const;
Vector3 & operator/=( float );
Vector3 operator/( float )const;
Vector3 cross( Vector3 const & )const;
float dot( Vector3 const & )const;
};

struct State
{
Vector3 position_;
Vector3 linear_momentum_;

State & operator+=( State const & );
State operator+( State const & )const;
State & operator-=( State const & );
State operator-( State const & )const;
State & operator*=( float );
State operator*( float )const;
State & operator/=( float );
State operator/( float )const;
};

struct RigidBody
{
float mass_;
State state_;
Vector3 forceAcc_;

RigidBody();
RigidBody( float mass, State const & state );

void apply_force( Vector3 const & force, Vector3 const & location );
State state_derivative( State const & state );
void update( double dt );
};

State RigidBody::state_derivative( State const & state )
{
State derivative;
derivative.position_ = state.linear_momentum_ / mass_;
derivative.linear_momentum_ = forceAcc_;

return derivative;
}

//the value of forceAcc doesn't change during the call to update
//forceAcc(t) is constant between t and t+dt
void RigidBody::update( double dt )
{
// k1 = dxdt( x, t )
State k1 = state_derivative( state_ );
double hdt = dt / 2.0;

//x1 = x + k1 * dt / 2;
State x1 = state_ + k1 * hdt;

// k2 = dxdt( x1, t + dt / 2 );
State k2( state_derivative( x1 ) );

// x2 = x + k2 * dt / 2;
State x2 = state_ + k2 * hdt;

// k3 = dxdt( x2, t + dt / 2 );
State k3( state_derivative( x2 ) );

// x3 = x + k3 * dt;
State x3 = state_ + k3 * hdt;

// k4 = dxdt( x3, t + dt );
State k4( state_derivative( x3 ) );

// x4 = x + k4 * dt;
State x4 = state_ + k4 * dt;

//xf = x + ( k1 + k4 ) * dt / 6 + ( k2 + k3 ) * dt / 3
State u_state = state_;
u_state += ( k1 + k4 ) * ( dt / 6.0 );
u_state += ( k2 + k3 ) * ( dt / 3.0 );

state_ = u_state;
}

int main( int argc, char ** argv )
{

return 0;
}
[/source]

Share on other sites
JimmyDeemo    156
[quote name='johnstanp' timestamp='1305822367' post='4813056']
I haven't taken a lookt at your code: it is not necessary, since I'll give a general answer.
You can consider that the total force applied to the rigid body is constant between t and t+dt: this is what I do for my update function using RK4.
[/quote]

Thanks for your reply, as well as the code but i am afraid i find you answer quite cryptic. So let me just see if i am on the right track here.

What you are saying is that for a given state of a Rigid body object, forces would be accumulated in the cycle (e.g. +thrust, -drag, etc.), when it comes to each of the four derivatives it can be considered that the forces are the same and are not influenced by the time. Is that right? So therefore the momentum is just the same as the total forces acting on the object.

If i have understood that right, then that's great and I'm pretty sure i can sort my code for that. But i just wondered why this was the case?

Share on other sites
johnstanp    267
[quote name='JimmyDeemo' timestamp='1305832100' post='4813131']
[quote name='johnstanp' timestamp='1305822367' post='4813056']
I haven't taken a lookt at your code: it is not necessary, since I'll give a general answer.
You can consider that the total force applied to the rigid body is constant between t and t+dt: this is what I do for my update function using RK4.
[/quote]

Thanks for your reply, as well as the code but i am afraid i find you answer quite cryptic. So let me just see if i am on the right track here.

What you are saying is that for a given state of a Rigid body object, forces would be accumulated in the cycle (e.g. +thrust, -drag, etc.), when it comes to each of the four derivatives it can be considered that the forces are the same and are not influenced by the time. Is that right? So therefore the momentum is just the same as the total forces acting on the object.

If i have understood that right, then that's great and I'm pretty sure i can sort my code for that. But i just wondered why this was the case?
[/quote]
Yes, you accumulate the forces before the update and keep the accumulator constant for the RK4 integrator (update function).
The state of the rigid body is constant untill you assign the result of the integration step to it i.e untill this->state_ = u_state;

But no, the momentum is not the same...Since x1 != x2 and x2 != x3 and x3 != x4 if forceAcc_ != 0. Do the computing yourself, I mean by hand, that'll be a lot more clear. This is how I understood it. Take a pen, a paper and a calculator and do it by hand.
I could have used a function that computes the force for a given state of the rigid body and an instant (time), meaning that I would have known beforehand the forces that would have acted upon the body, but you would have been distracted by the details of the implementation (the use of a function pointer or a functor). If you understand this first part, then I will post a version where the force actually varies during t and t+dt. But you have to understand that unless the forces acting upon the body are known beforehand, you are the only person to actually decide how they will vary between two updates. That'll be clearer after but for now, do the previous computations with a pen and a paper.

Share on other sites
JimmyDeemo    156
[quote name='johnstanp' timestamp='1305833799' post='4813152']
Yes, you accumulate the forces before the update and keep the accumulator constant for the RK4 integrator (update function).
The state of the rigid body is constant untill you assign the result of the integration step to it i.e untill this->state_ = u_state;

But no, the momentum is not the same...Since x1 != x2 and x2 != x3 and x3 != x4 if forceAcc_ != 0. Do the computing yourself, I mean by hand, that'll be a lot more clear. This is how I understood it. Take a pen, a paper and a calculator and do it by hand.
I could have used a function that computes the force for a given state of the rigid body and an instant (time), meaning that I would have known beforehand the forces that would have acted upon the body, but you would have been distracted by the details of the implementation (the use of a function pointer or a functor). If you understand this first part, then I will post a version where the force actually varies during t and t+dt. But you have to understand that unless the forces acting upon the body are known beforehand, you are the only person to actually decide how they will vary between two updates. That'll be clearer after but for now, do the previous computations with a pen and a paper.
[/quote]

Apologies, what i meant was 'In each derivative calculation, the momentum is just the same as the forces, because time isn't taken into account.'. I will take some time an do as you suggested. The slight hurdle i have got to get over is that your method looks a little different to mine (and the articles i followed), but i am sure its all there. Once i get my head around it i will post again.

Share on other sites
johnstanp    267
[quote name='JimmyDeemo' timestamp='1305837080' post='4813178']
Apologies, what i meant was 'In each derivative calculation, the momentum is just the same as the forces, because time isn't taken into account.'. I will take some time an do as you suggested. The slight hurdle i have got to get over is that your method looks a little different to mine (and the articles i followed), but i am sure its all there. Once i get my head around it i will post again.
[/quote]
Keep in mind that the force is the derivative of the linear momentum. We're computing the derivative of a state: that's why the derivative of the state is multiplied by a time to yield a delta state...
The force is the derivative of the linear momentum with respect to time: it has the dimension of a linear momentum divided by time. So yes, the derivative of the linear momentum stay constant, but the linear momentum doesn't stay constant: just like the slope of a line is constant, but not the values of y ( y = m*x + b, with m a constant =>
dy/dt = m and dy/dt is constant)
You may know perfectly well those notions, but using the right terms remove any possibility of confusion.

P.S: State x4 isn't needed: so that line should be commented.