在一个单自由度质量-弹簧系统中实现半隐式后向欧拉方法

Question

在一个单自由度质量-弹簧系统中实现半隐式后向欧拉方法

algorithmsimulationinterpolation

6

我有一个简单的弹簧系统，其中两个点通过弹簧连接。其中一个点固定在天花板上，因此我想使用数值方法计算第二个点的位置。所以，基本上我得到了第二个点的位置和速度，并想知道这两个值在一个时间步长后如何更新。

下列力会影响该点：

重力力，由-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 = ma。

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

- Etan

已删除spring标签，因为它与Java Spring框架相关。 - Sean Patrick Floyd

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Jonathan Dursi · Accepted Answer

我可以帮助您翻译。以下是需要翻译的内容：

因此，最好先考虑完全隐式方法，然后再转向半隐式方法。

隐式欧拉公式为（我们将其称为方程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)

您好，以下是您需要翻译的内容：

$\left ( \begin{array}{cc} 1 & -\Delta t \ \frac{k}{m} & 1 - \frac{d}{m}\end{array} \right ) \left ( \begin{array}{c} x^\mathrm{new} \ v^\mathrm{new} \end{array} \right ) = \left ( \begin{array}{c} x^\mathrm{old} \ v^\mathrm{old} - g \Delta t\end{array} \right )$

在这个矩阵中，您已知所有元素和右侧的数值，您只需要解出向量(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)ⁿ)/f'((x,v)ⁿ)

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

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

而f'(newPos, newVel)是对应矩阵的Jacobian矩阵

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

进行半隐式计算的过程与标准计算类似，但稍微容易一些 - 方程(1)中并非所有右手项都是新量。通常的做法是：

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

例如，速度取决于位置的旧时间值，而位置取决于速度的新时间值。（这与“跳跃式”积分非常相似。）使用这组略微不同的方程组应该很容易完成上述步骤。基本上，在上面的矩阵中，k/m项消失了。