• Advertisement
Sign in to follow this  

Runge-Kutta integration

This topic is 4306 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

Consider the standard equations of motion: dx/dt = v, and dv/dt = a, where x is the position, v is velocity, a is acceleration. I'm interested in using the fourth-order Runge-Kutta (RK4) method to solve (i.e. integrate) these equations. My question is, within a single time-step, do I apply the RK4 method to each equation separately, or does the result of an application of RK4 to the first equation (dx/dt = v) affect the application of RK4 to the second equation (dv/dt = a)? I hope that makes sense.

Share this post


Link to post
Share on other sites
Advertisement
You apply each one separately. When solving higher-order differential equations, you treat them as a system of first-order differential equations X'=F(X). So each element of the vector X (in your case, x and v) is updated simultaneously, and none of the newly calculated values affects any other math until the following timestep.

Share this post


Link to post
Share on other sites
Example RK4 implementation:


struct RBState {
RigidBody * rb; // Pointer to rigid body associated with this state.
// State Variables at time t:
Vec3 center; // x : Center of rotation and center of mass.
Quat4 rot; // R : Rotational position.
Vec3 linM; // P : Linear momentum.
Vec3 rotM; // L : Rotional momentum

RBState() : rb(0) {}
RBState(RigidBody * irb,const Vec3 & c,const Quat4 & r,const Vec3 & lm,const Vec3 & rm) : rb(irb), center(c),rot(r),linM(lm),rotM(rm) {}

RBState derivative(flt currentTime);
};

inline RBState operator * (const RBState & s,flt n) {
return RBState(s.rb,s.center*n,s.rot*n,s.linM*n,s.rotM*n);
}
inline RBState operator * (flt n,const RBState & s) {
return RBState(s.rb,s.center*n,s.rot*n,s.linM*n,s.rotM*n);
}
inline RBState operator + (const RBState & a,const RBState & b) {
return RBState(a.rb,a.center+b.center,a.rot+b.rot,a.linM+b.linM,a.rotM+b.rotM);
}

inline void RigidBody::updateDerived(RBState & s) {
mrot = Mat3(s.rot.normalize());
rotVel = scaleLocal(s.rotM,rotMassIdealInv,mrot);
linVel = s.linM*linMassInv;
}// RigidBody::updateDerived

inline RBState RBState::derivative(flt currentTime) {
rb->clearForceAccum();
rb->computeForces(currentTime);
rb->updateDerived(*this);
Quat4 dRot = .5f*mulQuat4BwIs0(rot,rb->rotVel);
return RBState(rb,rb->linVel,dRot,rb->linForce,rb->rotForce);
} // RBState::derivative

#endif

class GODEI {
protected:
flt h;
flt time;
public:
GODEI() : h(0.f) {}
virtual void setTime(flt deltaTime,flt currentTime=0) {
h = deltaTime;
time = currentTime;
} // setTime
virtual void update(void)=0;
};

template <class T>
class TGODEI : public GODEI {
protected:
T * current;
public:
TGODEI() : current(0) {}
void setCurrent(T * icurrent) {
current = icurrent;
} // setCurrent
void setCurrentAndTime(T * icurrent,flt deltaTime,flt currentTime=0) {
setCurrent(icurrent);
setTime(deltaTime,currentTime);
} // setCurrentAndTime
};

template <class T>
class EulerODEI : public TGODEI<T> {
public:
void update(void) {
*current = *current + h*current->derivative(time);
} // update
};

template <class T>
class MidpointODEI : public TGODEI<T> {
protected:
flt h2;
public:
MidpointODEI() {}
void setTime(flt deltaTime,flt currentTime=0) {
h = deltaTime;
h2 = h*.5f;
time = currentTime;
} // setStepSize
void update(void) {
T k = *current + h2*current->derivative(time);
*current = *current + h*k.derivative(time+h2);
} // update
};

template <class T>
class RK4ODEI : public TGODEI<T> {
flt h2,h6;
T k1,k2,k3,k4;
public:
RK4ODEI() {}
void setTime(flt deltaTime,flt currentTime=0) {
h = deltaTime;
h2 = deltaTime*.5f;
h6 = deltaTime*(1.f/6.f);
time = currentTime;
} // setStepSize
void update(void) {
flt time2 = time+h2;
k1 = current->derivative(time);
k2 = (*current + k1*h2).derivative(time2);
k3 = (*current + k2*h2).derivative(time2);
k4 = (*current + k3*h).derivative(time);
*current = *current + h6*(k1 + 2.f*(k2 + k3) + k4);
} // update
};





See Baraff's paper here for physics system derivations and implementation details.

You can probably use Newton-Stormer-Verlet for a game-physics type of application (instead of Runge-Kutta). See this thread for more info.

Share this post


Link to post
Share on other sites
Sign in to follow this  

  • Advertisement