使用python将弹簧物理应用于四元数

Question

我想在 python 中创建一个简单的物理系统，它的工作quaternions方式与velocity/position. 它的主要目标是模拟一个被拖动的对象，并随着时间的推移尝试赶上另一个对象。模拟使用 3 个变量：：k弹簧常数，d：阻尼因子，和m：拖动对象的质量。

使用经典的欧拉积分，我可以用以下方法解决位置：

import numpy as np
from numpy.core.umath_tests import inner1d

# init simulation values
dt = 1./30. # delta t @ 30fps
k = 0.5 # Spring constant
d = 0.1 # Damping
m = 0.1 # Mass

p_ = np.array([0,0,0]) # initial position of the dragged object
p  = np.array([1,2,0]) # position to catch up to, in real life this value could change over time
v  = np.array([0,0,0]) # velocity buffer (init speed is assumed to be 0)

# Euler Integration
n = 200 # loop over 200 times just to see the values converge
for i in xrange(n):
    x  = (p-p_)
    F  = (k*x - d*v) / m # spring force
    v  = v + F * dt # update velocity
    p_ = p_ + v * dt # update dragged position

    print p_ # values oscillate and eventually stabilize to p
    
    
    

这对位置很有用。通过改变k，d而且m我可以得到一个迅捷/重结果和总体来说，我很满意我的感觉：

现在我想对quaternions. 因为我没有用 做物理的经验quaternions，所以我走上了天真的道路，并应用了相同的函数，并进行了一些修改来处理quaternion翻转和归一化。

# quaternion expressed as x,y,z,w
q_ = np.array([0., 0., 0., 1.]) # initial quaternion of the dragged object
q  = np.array([0.382683, 0., 0., 0.92388]) # quaternion to catch up to (equivalent to xyz euler rotation of 45,0,0 degrees)
v  = np.array([0.,0.,0.,0.]) # angular velocity buffer (init speed is assumed to be 0)

# Euler Integration
n = 200
for i in xrange(n):
    
    # In a real life use case, q varies over time and q_ tries to catch up to it.
    # Test to see if q and q_ are flipped with a dot product (innder1d).
    # If so, negate q
    if inner1d(q,q_) < 0:
        q = np.negative(q)
    
    x  = (q-q_)
    F  = (k*x - d*v) / m # spring force
    v  = v + F * dt # update velocity
    q_ = q_ + v * dt # update dragged quaternion
    q_ /= inner1d(q_,q_) # normalize
    
    print q_ # values oscillate and eventually stabilize to q

令我惊讶的是，它给了我非常不错的结果！

因为我凭直觉行事，所以我确信我的解决方案是有缺陷的（比如如果 q 和 q_ 是对立的）并且有一种正确/更好的方法来实现我想要的。

题：

模拟弹簧力的正确方法是quaternions（至少）考虑mass拖动对象，弹簧stiffness和damping factor.

实际的 python 代码将不胜感激，因为我很难阅读博士论文:)。对于quaternion操作，我通常会参考Christoph Gohlke 的优秀库，但请随时在您的答案中使用任何其他库。

Answer 1

“更好”在这里实际上是非常主观的。

你有点逃避四元数的概念，因为你的步长和位移很小。根据您的应用程序，这实际上可能没问题（游戏引擎经常利用这样的技巧来使实时计算更容易）但是如果您的目标是准确性，或者您想增加步长而不是获得不稳定的结果，您将需要使用四元数，因为它们应该被使用。

正如@z0r 在评论中所解释的那样，由于四元数通过乘法变换旋转，因此它们之间的“差异”是乘法逆 - 基本上是四元数除法。

qinv = quaternion_inverse(q)  # Using Gohlke's package 
x = quaternion_multiply(q_, qinv)

现在，对于 small theta，theta =~ sin(theta)，x只要差异很小，这与减法的结果没有太大区别。像这样滥用“小角度定理”经常用于所有类型的模拟，但重要的是要知道何时打破它们以及它对模型施加了哪些限制。

加速度和速度仍然增加，所以我认为这仍然有效：

F  = (k*x - d*v) / m 
v  = v + F * dt 

组成单位旋转

q_ = quaternion_multiply(q_, unit_vector(v * dt)) # update dragged quaternion

再一次，对于小角度（即dt与速度相比很小），总和和乘积非常接近。

如有必要，然后像以前一样标准化。

q_ = unit_vector(q_)

我认为这应该可行，但会比您之前的版本慢一点，并且可能会产生非常相似的结果。

Answer 2

对于“模拟四元数上的弹簧力的正确方法是什么”这个问题，答案是正确写下势能。

1. 四元数通过运算给出物体的方向

   orientation = vector_part(q_ r q_*)
   

   其中星号表示共轭，r 是固定方向（例如“沿 z 的单位向量”，它在系统中对于所有对象必须是唯一的）。q_、r 和 q_* 的乘法被假定为“四元数乘法”。

2. 例如，定义具有该方向的能量函数

   energy = dot_product(orientation, unit_z)
   

3. 取能量相对于四元数的“减去导数”，您将获得施加到系统的力。

您当前的代码执行“四元数空间中的阻尼振荡器”，这可能是解决您问题的好方法，但它不是物体上的弹簧力:-)

PS：评论太长了，希望对你有帮助。PS2：我没有使用直接代码来解决上述问题，因为（i）我发现阅读上面的库文档并不容易，（ii）问题的第一部分是做数学/物理。