代码之家 › 专栏 › 技术社区 › Etan bbum

半隐式向后欧拉在一自由度质量弹簧系统中的实现

interpolation simulation algorithm

Etan bbum · 技术社区 · 14 年前

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

以下力作用于该点:

重力,由 -g * m
弹簧力,由 k * (l - L)
阻尼力,由 -d * v

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

例如,应用显式Euler,可以导出以下内容:

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

不过,我现在想用 半隐式向后欧拉

1 回复 | 直到 14 年前

Jonathan Dursi 14 年前

因此,最容易看到的可能是,从首先考虑完全隐式方法,然后再考虑半隐式方法。

隐式Euler应该有(我们称之为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)

$\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 )$

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

也就是说,f(newPos,newVelocity)=(0,0)。你有一个前一个值作为开始猜测(oldPos,oldpvelocity)。现在你只想重复

^n+1个 =(x,v) ^n个 +f((x,v) ^n个 ^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项消失了。