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

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

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

隐式欧拉将有（让我们称这些方程（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 求解器来完成此操作（在这个简单的情况下，可以手动进行高斯消除）。但既然你提到了雅可比行列式，你可能想用牛顿-拉夫森迭代或类似的方法来解决这个问题。

在这种情况下，您本质上是要求解方程的零点

((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)^{n)/f'((x,v)n)}

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

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 项消失了。