javascriptsimulationgame-enginephysicsrunge-kutta

Runge Kutta implementation is less accurate than Euler implementation

So I am trying to code a 2D Physics engine, and I am in need of a more accurate ODE solver so that my spring simulation does not become unstable. I want to make a stiff, or rigid spring that acts like a pendulum, which is simulated more accurately by make the stiffness as high as possible.

Currently, I have implemented the following RK4 method, but through testing, it appears that the Spring is more unstable and breaks at lower stiffness values compared to Euler's method:

const acc = this.getTotalAcc();

const k1 = acc.divide(state.preset.options.stepsPerFrame);
const k2 = acc.add(k1.divide(2)).divide(state.preset.options.stepsPerFrame);
const k3 = acc.add(k2.divide(2)).divide(state.preset.options.stepsPerFrame);
const k4 = acc.add(k3.divide(2)).divide(state.preset.options.stepsPerFrame);

const deltaVel = k1
  .add(k2.multiply(2))
  .add(k3.multiply(2))
  .add(k4)
  .divide(6);
this.vel = this.vel.add(deltaVel);

const deltaPos = this.vel.divide(state.preset.options.stepsPerFrame);
this.pos = this.pos.add(deltaPos);

I checked this a couple of times, and I don't believe I made any calculation mistakes, but I could be wrong.

For reference, here is my Euler's implementation:

this.vel = this.vel.add(acc.divide(state.preset.options.stepsPerFrame))
this.pos = this.pos.add(this.vel.divide(state.preset.options.stepsPerFrame));

Here is my spring interaction function in case you need to see how it works. this.end and this.start are the objects attached to each end of the spring:

interact() {
  const distance = this.end.pos.subtract(this.start.pos);
  const force = this.stiffness * (distance.magnitude() - this.length);

  // Distribute tensions between masses
  const startTension =
    (this.start.inverseMass /
      (this.start.inverseMass + this.end.inverseMass)) *
    force;
  const endTension =
    (this.end.inverseMass / (this.start.inverseMass + this.end.inverseMass)) *
    force;

  this.start.accs[this.name] = distance.unit().multiply(startTension);
  this.end.accs[this.name] = distance.unit().multiply(-endTension);
}
Solution

  • The k4 computation should not have division by 2

    You are using division by StepsPerFrame instead of multiplication by the step size. This needs also to apply to the final computation, after the division by 6.

    You should also compute your position updates with the RK4 method. For mechanical systems that are described by second order equations, this can be integrated into the velocity computation as described in https://scicomp.stackexchange.com/questions/34257/solving-coupled-odes-using-runge-kutta-method, general remarks in https://math.stackexchange.com/a/4177218/115115

    These are the structural errors regarding the parameters of the RK4 method. You have a deeper error in applying the method. It should be something like

    k2 = getTotalAcc( vel.add(k1.divide(2).divide(stepsPerFrame)) )

    or

    k2 = getTotalAcc( vel.add(k1.divide(2)) ).divide(stepsPerFrame)

    depending on if the k vectors are the derivatives or the increments (equal derivative times time step).

    The derivatives are computed in every stage of the RK4 step. What you do makes no sense, try to translate the function calls back into mathematical formulas to see this.

    Making the springs stiff literally increases the stiffness of the ODE system. You need to keep the step size small, or stepsPerFrame large enough, to remain inside the stability region of the method. This can not be directly checked, you would need to compare computation runs with different step sizes.

    Your Euler method is semi-explicit or symplectic Euler.