Runge Kutta 实现不如 Euler 实现准确

Question

所以我正在尝试编写一个 2D 物理引擎，并且我需要一个更准确的 ODE 求解器，以便我的弹簧模拟不会变得不稳定。我想制作一个像钟摆一样坚硬或刚性的弹簧，通过使刚度尽可能高来模拟更准确。

目前，我已经实现了以下 RK4 方法，但通过测试，与欧拉方法相比，弹簧似乎更不稳定并且在较低刚度值下断裂：

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

我检查了几次，我不相信我犯了任何计算错误，但我可能是错的。

作为参考，这是我的欧拉实现：

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

这是我的 spring 交互函数，以防您需要了解它是如何工作的。

this.end

和

this.start

是连接到弹簧两端的物体：

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

Answer 1

k4 计算不应除以 2

您正在使用除以 StepsPerFrame 而不是乘以步长。这也需要应用于除以 6 后的最终计算。

您还应该使用 RK4 方法计算您的位置更新。对于由二阶方程描述的机械系统，可以将其集成到速度计算中，如 https://scicomp.stackexchange.com/questions/34257/solving- Coupled-odes-using-runge-kutta-method 中所述，一般性评论 https://math.stackexchange.com/a/4177218/115115

这些是关于 RK4 方法参数的结构错误。您在应用该方法时存在更深层次的错误。应该是这样的

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

或

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

取决于 k 个向量是导数还是增量（等于导数乘以时间步长）。

导数在 RK4 步骤的每个阶段进行计算。你所做的毫无意义，尝试将函数调用转换回数学公式来看看这一点。

---

使弹簧变硬实际上会增加 ODE 系统的刚度。您需要保持步长较小或stepsPerFrame足够大，以保持在方法的稳定区域内。这不能直接检查，您需要比较不同步长的计算运行。

---

您的欧拉方法是半显式或辛欧拉方法。