在 1-DOF 质量弹簧系统中实现半隐式后向 Euler

Question

我有一个简单的（质量）弹簧系统，它有两个与弹簧相连的点。一个点固定在天花板上，所以我想用数值方法计算第二个点的位置。所以，基本上我得到了第二个点的位置和它的速度，并且想知道这两个值在一个时间步之后是如何更新的。

以下力作用于该点：

• 引力，由-g * m 给出
• 弹簧力，由k * (l - L) 给出，k 是刚度，l 是当前长度，L 是初始长度
• 阻尼力，由-d * v 给出

总结一下，这导致

• F = -g * m + k * (l - L)
• Fd = -d * v

应用例如显式欧拉，可以得出以下结论：

• newPos = oldPos + dt * oldVelocity
• newVelocity = oldVelocity + dt * (F + Fd) / m，使用F = m * a。

但是，我现在想使用 半隐式后向欧拉，但无法准确确定从哪里导出雅可比矩阵等。

Answer 1

因此，从首先考虑完全隐式方法，然后再考虑半隐式方法，这可能是最容易理解的。

隐式欧拉将有（我们称这些 eqn (1)）：

newPos = oldPos + dt * newVelocity
newVelocity = oldVelocity + dt * (-g * m + k*(newPos - L) - d*newVelocity)/m

现在让我们只测量相对于 L 的位置，这样我们就可以摆脱 -kL 项。重新排列我们最终得到

(newPos, newVelocity) - dt * (newVelocity, k/m newPos - d/m newVelocity) = (oldPos, oldVelocity - g*dt)

并将其放入矩阵形式

((1,-dt),(k/m, 1 - d/m)).(newPos, newVelocity) = (oldPos, oldVelocity -g*dt)

您知道矩阵中的所有内容以及 RHS 上的所有内容，您只需要求解向量（newPos、newVelocity）。您可以使用任何 Ax=b 求解器来执行此操作（在这个简单的情况下，手动高斯消除工作）。但是由于您提到雅可比行列式，您可能希望通过 Newton-Raphson 迭代或类似的方法来解决这个问题。

在这种情况下，您实际上是在寻找解方程的零点

((1,-dt),(k/m, 1-d/m)).(newPos, newVelocity) - (oldPos, oldVelocity -g*dt) = 0

也就是说，f(newPos, newVelocity) = (0,0)。您有一个以前的值可用作起始猜测，（oldPos，oldVelocity）。现在你只想迭代

(x,v)ⁿ⁺¹ = (x,v)ⁿ + f((x,v)ⁿ)/ f'((x,v)ⁿ)

直到你得到一个足够好的答案。这里，

f(newPos,newVel) = ((1,-dt),(k/m, 1-d/m)).(newPos, newVelocity) - (oldPos, oldVelocity -g*dt)

而f'(newPos, newVel)是矩阵对应的雅可比行列式

((1,-dt),(k/m, 1-d/m))

完成半隐式的过程是相同的，但更容易一些 - 并非 eqns (1) 中的所有 RHS 项都是新量。通常的做法是

newPos = oldPos + dt * newVelocity
newVelocity = oldVelocity + dt * (-g * m + k*oldPos - d*newVelocity)/m

例如，速度取决于位置的旧时间值，而位置取决于速度的新时间值。 （这与“跨越式”集成非常相似。）使用这组稍有不同的方程组，您应该能够非常轻松地完成上述步骤。基本上，上面矩阵中的 k/m 项消失了。