在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

因此,最简单的方法是首先考虑完全隐式方法,然后转向半隐式方法.

隐含的欧拉会(让我们称之为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)^{n + 1} =(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))

通过半隐式的过程是相同的,但更容易一点 - 并非所有的方程(1)中的RHS术语都是新的数量.它通常采用的方式是

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

例如,速度取决于位置的旧时间值,以及速度的新时间值上的位置.(这与"越级"整合非常相似.)你应该能够通过这个稍微不同的方程组很容易地完成上述步骤.基本上,上面矩阵中的k/m项就会消失.