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

Question

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

使用经典的euler积分，我可以求解以下位置：

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_是对立的)，并且有一种正确/更好的方法来实现我想要的东西。

问题：

在考虑(至少)被拖动对象的mass、弹簧stiffness和damping factor的情况下，模拟quaternions上的弹力的正确方法是什么。

实际的python代码将非常感谢，因为我有很大的困难阅读PhD论文:)。此外，对于quaternion操作，我通常指的是Christoph Gohlke's excellent library，但请在您的答案中随时使用其他任何其他操作。

Answer 1

“更好”实际上在这里是相当主观的。

你可以忽略四元数的概念，因为你的步长和位移都很小。根据您的应用程序，这实际上可能是可以的(游戏引擎经常利用这样的技巧来使实时计算更容易)，但如果您的目标是准确性，或者您希望增加步长而不得到不稳定的结果，则需要使用四元数，因为它们是应该使用的。

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

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

现在，就像小的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. 四元数通过操作给出对象的方向

方向= vector_part(q_ r q_*)

其中星号表示共轭，r表示固定方向(例如，“沿z的单位向量”，它在系统中对于所有对象都必须是唯一的)。例如，假设q_、r和q_*的乘积是具有该取向的能量函数的四元数

能量=dot_product(unit_z，)

• 取能量相对于四元数的“减去导数”，你就有了作用于系统的力。

您当前的代码做了一个“四元数空间中的阻尼振荡器”，这可能是一个很好的问题解决方案，但它不是对象上的弹力:-)

附言:评论太长了，我希望它能有所帮助。PS2:我没有使用直接的代码来解决上面的问题，因为(i)我发现只读上面的库文档并不容易，(ii)问题的第一部分是做数学/物理。