# Runge-Kutta integration

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

## 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 on other sites
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 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::updateDerivedinline 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#endifclass 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.

##### Share on other sites
Thanks for setting me straight. That thread is pretty damned interesting.

1. 1
Rutin
46
2. 2
3. 3
4. 4
5. 5

• 13
• 10
• 12
• 10
• 13
• ### Forum Statistics

• Total Topics
632993
• Total Posts
3009762
• ### Who's Online (See full list)

There are no registered users currently online

×