Skip to main content
gamedev.net
🔒 Locked 🎮 Unity

Help me implement RK4 here

Started by xissburg Oct 16, 2005 at 12:24 PM 11 replies 4.9k views
Original Post
xissburg
xissburg
Hello... I was reading some articles about RK4 numerical integrator but I didnt understand it too much :( ...I learnt RK2 but RK4 is really strange for me :P So I want your help to implement RK4 in a simple particle system as an example. Look at this struct: 

struct PARTICLE
{
    VECTOR3 position;
    VECTOR3 speed;
    VECTOR3 force;
    float   invMass;
}Particle;
As you can see, my struct has the position, the speed, the inverse mass(for faster calculations:) and the force(the resultant force applied on the particle) of the particle. I want to integrate the position through the velocity and the velocity through the acceleration(force*invMass). Using RK2 I got this: 

.
.
.
//dt is the timestep
particle.speed += particle.force * particle.invMass * 0.5f * dt;
particle.position += particle.speed * dt;
particle.speed += particle.force * particle.invMass * 0.5f * dt;
//then update forces and so on...
.
.
.
Using RK2 it works correctly but I wanna know how the code would be when using RK4 . Please show me the way!! Thanks....
.
Zipster
Zipster
This contains a good series of articles on game physics and includes a section on writing an RK4 integrator. The basics of the itegrator is that you sample the slope of the curve between t and dt at four points, and use the weighted average as the final velocity for your object.
xissburg
xissburg
huh...I have already read this article but it didnt help me :D :D :D

I'm a little donkey on math/physics :P ....cant you show me the code for this case?!?
.
Zipster
Zipster
Well, the article contains all the code you need for the integrator. And I'm afraid that I wouldn't do a much better job explaining all the various issues than the article already does.

Is there anything specific you understand? You say you implemented an RK2 integrator, RK4 simply takes two more sample points and averages them in.
xissburg
xissburg
wow....I've read that again but I didn't learn it yet :( .. its complex D:

For example I dont understand that 'float acceleration(const State &state, float t)' function. Is the 't' the time that has passed when the acceleration begun acting on the object? I'm a lot confused about!!! please help :P

.
Zipster
Zipster
Ah, I see. Time 't' is the time since the acceleration you're interested in began. In the case of the article, they don't even use 't' in their calculations since they're modeling a spring-damper system, where the time since the start of some event doesn't matter, because the acceleration depends only on the properties of the input object. However, if you want to model simple gravity the function would look like this:
float acceleration(const State &state, float t){   return -10.0f;}

Where that value is anything that most accurately simulates gravity in your simulation. As you can see, time doesn't even matter since acceleration due to gravity doesn't vary as a function of time. However, let's say you want to model rocket thrusters, and that the force exerted by the thrusters does change as time goes on (let's say they're being turned off). And let's say that when the rocket first starts it exerts 10 units of acceleration, and that acceleration decreases linearly to 0 units in 10 seconds. Then your acceleration would look like this:
float acceleration(const State &state, float t){   float accel = t - 10;    return (accel > 0.0f) ? accel : 0.0f;}

What would you pass in for time 't'? The time since the rocket was started. When 't' is zero, you get the full acceleration, and when 't' is ten you get zero acceleration. When 't' is greater than ten you get a negative acceleration, but the conditional check will always make sure a non-negative acceleration is returned.

So to answer your question, time 't' depends exactly on where the acceleration is coming from. In the case of gravity it's irrelevant, but in the case of acceleration that varies as a function of time, it can be important, and you were correct in assuming it starts when the acceleration starts.
Koshmaar
Koshmaar
Sorry for hijacking thread xissburg, but since you already started thread about RK4 from that tutorial, I have one little question to someone who has good knowledge about various integration methods: wouldn't it be better to use Verlet? It's much simplier to write (one line), faster to execute and even more stable. So, is there any special reason against doing the same as Gaffer in his articles (fixed time movement etc), but only changing the integrator to Verlet, instead of RK4?
John Schultz
John Schultz
Quote:
Original post by Koshmaar
Sorry for hijacking thread xissburg, but since you already started thread about RK4 from that tutorial, I have one little question to someone who has good knowledge about various integration methods: wouldn't it be better to use Verlet? It's much simplier to write (one line), faster to execute and even more stable. So, is there any special reason against doing the same as Gaffer in his articles (fixed time movement etc), but only changing the integrator to Verlet, instead of RK4?


You can take larger timesteps with RK4 vs. Verlet where the force changes over time (as with a spring). However, since RK4 requires four derivative evaluations per step, it may not buy much advantage over running Verlet with a dt scaled by 1/4 (resulting in the same amount of derivs per unit time).

I currently prefer using Verlet (or "Euler-Verlet") with clamped springs and rigid Verlet positional constraints over RK4. RK4 will still blow up given high-stiffness elements: forces/accels/velocities must still be clamped to keep the code 100% stable. RK4 might be more useful in situations where the weighted averaging (Taylor series) of the deriv samples provides more stability than running Verlet/Euler-Verlet at 1/4 the time step. While not an issue for games (perhaps desireable in some cases), RK4 loses energy and smooths/dampens the simulation (acts as a low pass filter).

xissburg,

Here is my implementation of Euler, RK2, and RK4 (two implementations):

The RK4 variant is fairly short/simple to code:

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


Full source:

// godei.h// Generalized ODE integrator.// Created 3/23/00 John Schultz#ifndef GODEI_H#define GODEI_H#ifndef DYNARRAY_H#include "dynarray.h"#endif#if 0// Example code using the TGODEI class.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#if 0  void update(void) {    flt time2 = time+h2;    k1 = h*current->derivative(time);    k2 = h*(*current + k1*.5f).derivative(time2);    k3 = h*(*current + k2*.5f).derivative(time2);    k4 = h*(*current + k3).derivative(time);    *current += (1.f/6.f)*(k1 + 2.f*(k2 + k3) + k4);  } // update#endif  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};struct StateArrayFLT {  DynArray<flt,false> states;  flt * pos;  flt deltaTime; // deltaTime (in seconds) for this state.  flt time; // Time (in seconds) for this state.  void setMemSize(int size) {    states.setSizeMaxCount(size/sizeof(flt));  } // setMemSize  void startIter(void) { // Must call startIter before read or writing arrays (iteration).    pos = &states[0];  } // startIter  void write(const void * mem,int size) { // Size in bytes (write to state array)    memcpy(pos,mem,size);    pos = (flt *)((char *)pos + size);  } // write  void read(void * mem,int size) { // Size in bytes (read from state array)    memcpy(mem,pos,size);    pos = (flt *)((char *)pos + size);  } // read  void * getPosPtr(int incSize) {    void * p = (void *)pos;    pos = (flt *)((char *)pos + incSize);    return p;  } // getPosPtr  void * getPosPtr(void) {    return (void *)pos;  } // getPosPtr  void incPosPtr(int size) {    pos = (flt *)((char *)pos + size);  } // incPosPtr  void print(void) {    lprintf("StateArrayFLT [%d]:\n",states.maxCount);    for (int i=0; i < states.maxCount; i++) {      lprintf("[%d] %6.4f\n",i,states);    } // for  } // print};// Get/Set rules: object.getState(currentState): get object's state and store in currentState.//                object.setState(newState):     set object's state to newState.struct ObjStateFLT {//  virtual int getType(void) { return TRigid; } // Object type (Rigid, spring, etc.).  virtual int getType(void)=0; // Object type (Rigid, spring, etc.).  virtual int stateSize(void)=0; // In bytes.  virtual void setState(StateArrayFLT & newState)=0;     // Set object to state in newState  virtual void getState(StateArrayFLT & currentState)=0; // Store object's state in currentState  virtual void updateDerived(StateArrayFLT & currentState)=0;  virtual void clearForceAccumulators(StateArrayFLT & currentState)=0;  virtual void computeForces(StateArrayFLT & currentState)=0;  virtual void derivative(StateArrayFLT & deltaState,StateArrayFLT & currentState)=0;};typedef ObjStateFLT * ObjStatePtr;typedef DynArray<ObjStatePtr> ObjStatePtrArray;class ODEIARRAY {public:  enum {TEMP1=5,TEMP2=6,NUMSTATEARRAYS=5+2}; // 0-4 are needed for RK4, 5,6 are temps for currentState and newState.  StateArrayFLT k[NUMSTATEARRAYS];  StateArrayFLT * currentState, * newState;  int getStateMemSize(ObjStatePtrArray & objects) {    int memSize = 0;    for (int i=0; i < objects.maxCount; i++) {      memSize += objects->stateSize();    } // for    return memSize;  } // getStateMemSize  void initMem(ObjStatePtrArray & objects) { // Objects MUST be exact multiples of FLT!    int memSizeNeeded = getStateMemSize(objects)/sizeof(k[0].states[0]);    int currMemSize   = k[0].states.maxCount/sizeof(k[0].states[0]);    if (memSizeNeeded > currMemSize) {      for (int i=0; i < NUMSTATEARRAYS; i++) {        k.states.setSizeMaxCount(memSizeNeeded);      } // for      currentState = &k[TEMP1];      newState     = &k[TEMP2];    } // if  } // initMem  void initMem(StateArrayFLT & currentState,ObjStatePtrArray & objects) {    int memSizeNeeded = getStateMemSize(objects)/sizeof(k[0].states[0]);    int currMemSize = currentState.states.maxCount/sizeof(currentState.states[0]);    if (memSizeNeeded > currMemSize) {      currentState.states.setSizeMaxCount(memSizeNeeded);    } // if  } // initMem  void getState(StateArrayFLT & currentState,ObjStatePtrArray & objects) {    currentState.startIter();    for (int i=0; i < objects.maxCount; i++) {      objects->getState(currentState);    } // for  } // getState  void setState(ObjStatePtrArray & objects,StateArrayFLT & newState) {    newState.startIter();    for (int i=0; i < objects.maxCount; i++) {      objects->setState(newState);    } // for  } // setState  void getStateCurrent(ObjStatePtrArray & objects) {    getState(*currentState,objects);  } // getStateCurrent  void setStateNew(ObjStatePtrArray & objects) {    setState(objects,*newState);  } // setStateNew  void swapCurrentNewStates(void) {    StateArrayFLT * temp = currentState;    currentState = newState;    newState = temp;  } // swapCurrentNewStates// *** Should virtualize these: make derived ODEPhysics... ***// *** Only derivative here. Will inline directly - no perf loss. ***  void updateDerived(ObjStatePtrArray & objects,StateArrayFLT & currentState) {    currentState.startIter();    for (int i=0; i < objects.maxCount; i++) {      objects->updateDerived(currentState);    } // for  } // updateDerived  void clearForceAccumulators(ObjStatePtrArray & objects,StateArrayFLT & currentState) {    currentState.startIter();    for (int i=0; i < objects.maxCount; i++) {      objects->clearForceAccumulators(currentState);    } // for  } // clearForceAccumulators  void computeForces(ObjStatePtrArray & objects,StateArrayFLT & currentState) {    currentState.startIter();    for (int i=0; i < objects.maxCount; i++) {      objects->computeForces(currentState);    } // for  } // computeForces  void derivative(StateArrayFLT & deltaState,StateArrayFLT & currentState,ObjStatePtrArray & objects,flt deltaTime,flt time) {    currentState.deltaTime = time;    currentState.time = time;    updateDerived(objects,currentState);    clearForceAccumulators(objects,currentState);    computeForces(objects,currentState);    deltaState.startIter();    currentState.startIter();    for (int i=0; i < objects.maxCount; i++) {      objects->derivative(deltaState,currentState);    } // for  } // derivative// **************************************************************  void addScaled(StateArrayFLT & dest,StateArrayFLT & added,StateArrayFLT & scaled,flt scale) {    for (int i=0; i < dest.states.maxCount; i++) {      dest.states = added.states + scaled.states*scale;    } // for  } // addScaled  void updateEuler(ObjStatePtrArray & objects,flt deltaTime,flt currentTime) {//    *current = *current + h*current->derivative(time);    derivative(k[1],*currentState,objects,deltaTime,currentTime); // k[1] = current.derivative(currentTime)    addScaled(*newState,*currentState,k[1],deltaTime);  // newState = current + deltaTime*k[1]    newState->time = currentTime;  } // updateEuler  void updateMidpoint(ObjStatePtrArray & objects,flt deltaTime,flt currentTime) {//    T k = *current + h2*current->derivative(time);//    *current = *current + h*k.derivative(time+h2);    flt h2 = deltaTime*.5f;    derivative(k[1],*currentState,objects,h2,currentTime); // k[1] = current.derivative(currentTime)    addScaled(k[2],*currentState,k[1],h2);              // k[2] = current + h2*k[1]    derivative(k[1],k[2],objects,deltaTime,currentTime+h2);       // k[1] = k[2].derivative(currentTime+h2)    addScaled(*newState,*currentState,k[1],deltaTime);  // newState = current + deltaTime*k[1]    newState->time = currentTime;  } // updateMidpoint  void updateRK4(ObjStatePtrArray & objects,flt deltaTime,flt currentTime) {//    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);    flt h2 = deltaTime*.5f;    flt h6 = deltaTime*(1.f/6.f);    flt time2 = currentTime+h2;    derivative(k[1],*currentState,objects,h2,currentTime); // k[1] = current.derivative(currentTime)    addScaled(k[0],*currentState,k[1],h2);              // k[0] = current + h2*k[1]    derivative(k[2],k[0],objects,h2,time2);             // k[2] = k[0].derivative(time2)    addScaled(k[0],*currentState,k[2],h2);              // k[0] = current + h2*k[2]    derivative(k[3],k[0],objects,deltaTime,time2);      // k[3] = k[0].derivative(time2)    addScaled(k[0],*currentState,k[3],deltaTime);       // k[0] = current + deltaTime*k[3]    derivative(k[4],k[0],objects,h6,currentTime);       // k[4] = k[0].derivative(currentTime)    for (int i=0; i < newState->states.maxCount; i++) {      newState->states = currentState->states                          + h6*(k[1].states                          + 2.f*(k[2].states + k[3].states)                          + k[4].states);    } // for    newState->time = currentTime;  } // updateRK4};// ****************************************************************************// ******************************** Test Code *********************************// ****************************************************************************struct TestState {  flt pos,vel,force;  TestState() : pos(0), vel(0), force(0) {}  TestState(flt ipos,flt ivel,flt iforce) : pos(ipos), vel(ivel), force(iforce) {}};class TestObjectState : public ObjStateFLT {public:  TestState state;  int getType(void) { return 0; } //TRigid; } // Object type (Rigid, spring, etc.).  int stateSize(void) {    return sizeof(state);  } // stateSize  void updateDerived(StateArrayFLT & currentState) {    currentState.incPosPtr(sizeof(state));  } // updateDerived  void derivative(TestState & delta,TestState & currentState) {    delta.pos   = currentState.vel;    delta.vel   = currentState.force;    delta.force = 0.f;  } // derivative  void derivative(StateArrayFLT & deltaState,StateArrayFLT & currentState) {    TestState * delta   = (TestState *)deltaState.getPosPtr(sizeof(*delta));    TestState * current = (TestState *)currentState.getPosPtr(sizeof(*current));    derivative(*delta,*current);  } // derivative  void getState(StateArrayFLT & currentState) { // Get local state from object and write to state array.    currentState.write(&state,sizeof(state));   // Write current state to state array.  } // getState  void setState(StateArrayFLT & newState) {     // Load new values into object's local state.    newState.read(&state,sizeof(state));        // Read new state from state array and store in local state.  } // setState  void clearForceAccumulators(StateArrayFLT & currentState) {    TestState * current = (TestState *)currentState.getPosPtr(sizeof(*current));    current->force = 0.f;  } // clearForceAccumulators  void computeForces(StateArrayFLT & currentState) {    TestState * current = (TestState *)currentState.getPosPtr(sizeof(*current));    current->force += -9.81f;  } // computeForces};class TestODEIARRAY {public:  ODEIARRAY godei;  ObjStatePtrArray objects;  enum {NUMOBJECTS=3};  TestObjectState testObjects[NUMOBJECTS];  TestODEIARRAY() {    init();  }  void init(void) {    objects.setSizeMaxCount(NUMOBJECTS);    for (int i=0; i < NUMOBJECTS; i++) {      objects = &testObjects    } // for    godei.initMem(objects);  } // init  void update(flt deltaTime,flt time) {    godei.getStateCurrent(objects);//    godei.updateEuler(objects,deltaTime,time);//    godei.updateMidpoint(objects,deltaTime,time);    godei.updateRK4(objects,deltaTime,time);    godei.setStateNew(objects);    print();  } // update  void print(void) {    for (int i=0; i < NUMOBJECTS; i++) {      TestObjectState & testObject = testObjects;      lprintf("Object [%d]: pos %6.3f vel %6.3f\n",i,testObject.state.pos,testObject.state.vel);    } // for  } // print  void test(void) {    flt steps = 10.f;    flt time = 10.f;    flt h = time/steps;    lprintf("\nStarting TestODEI:\n");    for (int i=0; i < steps; i++) {      update(h,flt(i)/flt(steps));      print();    } // for    lprintf("Stopping TestODEI.\n");  } // test};#endif// godei.h
Amazing Curves Racing
Zipster
Zipster
Quote:
Original post by Anonymous Poster
instead of
float accel = t - 10;
it should be
float accel = 10 - t;
I think

Depends on what your coordinate system is. The idea what that you have a linearly changing acceleration as a function of time. In this case the rocket was thrusting to toward the ground which is kinda of silly [smile]

Quote:
Sorry for hijacking thread xissburg, but since you already started thread about RK4 from that tutorial, I have one little question to someone who has good knowledge about various integration methods: wouldn't it be better to use Verlet? It's much simplier to write (one line), faster to execute and even more stable. So, is there any special reason against doing the same as Gaffer in his articles (fixed time movement etc), but only changing the integrator to Verlet, instead of RK4?

From my understanding, Verlet is faster and more stable but less accurate since you're not dealing with the true velocities. RK4 is stable and more accurate however slower since it requires 4 derivative evaluations.
Koshmaar
Koshmaar

Hmmm, general information about my background: I'm going to implement standard Gamasutra-Jacobsen Verlet particle system, with handling collisions by projection and solving constraints through infinitely stiff springs (or sth like that) and using all other tricks he have presented there.

Quote:
You can take larger timesteps with RK4 vs. Verlet where the force changes over time (as with a spring). However, since RK4 requires four derivative evaluations per step, it may not buy much advantage over running Verlet with a dt scaled by 1/4 (resulting in the same amount of derivs per unit time).


That's quite strange, I remember that somewhere I read about different integration methods, and there was written that explicit Euler is accurate only to 1 order, RK4 to X, while Verlet is accurate to X+1. I also thought that it means: Verlet is the most accurate integration method.


Quote:
I currently prefer using Verlet (or "Euler-Verlet") with clamped springs and rigid Verlet positional constraints over RK4. RK4 will still blow up given high-stiffness elements: forces/accels/velocities must still be clamped to keep the code 100% stable. RK4 might be more useful in situations where the weighted averaging (Taylor series) of the deriv samples provides more stability than running Verlet/Euler-Verlet at 1/4 the time step. While not an issue for games (perhaps desireable in some cases), RK4 loses energy and smooths/dampens the simulation (acts as a low pass filter).


Hmmm, what articles about physics stuff (Verlet, integration methods, physics in games etc.) would you recommend to read? Apart from Jacobsen, Glenn Fiedler and GameDev Lonelysocks's article about fixed time verlet vs. normal one? :-)
John Schultz
John Schultz
Quote:
Original post by Zipster
From my understanding, Verlet is faster and more stable but less accurate since you're not dealing with the true velocities. RK4 is stable and more accurate however slower since it requires 4 derivative evaluations.


With the same timestep, "Velocity-less (VL)" Verlet is (at least 4x) faster but less stable than RK4 for springs (Jacobsen-style rigid constraints are very fast and 100% stable with VL Verlet; not so with springs).

Updating velocity before position is sometimes called semi-implicit Euler; I use a variant I call Euler-Verlet (as dt is premultiplied into the velocity in the same way as VL Verlet). Verlet is considered very accurate, and is called a symplectic integrator, as it preserves area under the integration curve. RK4 on the other hand, is more stable, but less accurate, performing a low-pass filter on the resulting motion (and losing energy over time).

Keep in mind there are many Verlet variants (each with different behavior), as well as many Runge-Kutta methods.

Koshmaar,

If you look through my previous physics-related posts, you'll find many links to many papers (sorry I don't have time to dig them all up), as well as more info on Euler-Verlet (including source code examples).
Amazing Curves Racing
Koshmaar
Koshmaar
Quote:
Verlet is considered very accurate, and is called a symplectic integrator, as it preserves area under the integration curve. RK4 on the other hand, is more stable, but less accurate, performing a low-pass filter on the resulting motion (and losing energy over time).


Oh, now I understand! Accuracy != stability :-) while before I thought that accuracy <=> stability.

Quote:
If you look through my previous physics-related posts, you'll find many links to many papers (sorry I don't have time to dig them all up), as well as more info on Euler-Verlet (including source code examples).


Ok, but be warned that I'm going to look very closely at what you were doing in the past... [wink]

Topic Locked

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

Sign in to reply to this topic.