使用python将Spring物理应用于四元数

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和{}我可以得到一个更快/更重的结果，总体上我对它的感觉很满意：

现在我想对quaternions执行相同的操作。因为我没有使用quaternions进行物理操作的经验，所以我走了一条幼稚的道路，应用了相同的函数，只做了一些修改来处理quaternion翻转和规范化。在

^{pr2}$

令我大吃一惊的是，它给了我非常好的结果！在

因为我是凭直觉去做的，我确信我的解决方案是有缺陷的（比如如果q和q_u是对立的），并且有一个正确/更好的方法来实现我想要的。在

问题：

模拟quaternions上的弹簧力的正确方法是什么（至少）考虑到被拖动对象的mass，弹簧stiffness和damping factor。在

实际的python代码将非常感谢，因为我在读博士论文时遇到了很大的困难：）。另外，对于quaternion操作，我通常指的是Christoph Gohlke's excellent library，但请随意在您的答案中使用其他任何操作。在

Answer 1

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

1. 四元数通过操作为您提供对象的方向

   orientation = vector_part(q_ r q_*)
   

   其中星表示共轭，r是一个固定的方向（比如说“沿z的单位向量”，它在你的系统中对所有物体都是唯一的）。q_u，r和q_x的乘法被假定为“四元数乘法”。

2. 例如，用这个方向定义能量函数

   energy = dot_product(orientation, unit_z)
   

3. 取能量相对于四元数的“减导数”，你将有力作用于系统。

您当前的代码做了一个“四元数空间中的阻尼振子”，这可能是您的问题的一个很好的解决方案，但它不是一个物体上的弹簧力：-）

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

Answer 2

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

你有点像谋杀四元数的概念，因为你的步长和位移都很小。根据您的应用程序，这可能实际上是可以的（游戏引擎经常利用这样的技巧来简化实时计算），但如果您的目标是准确性，或者您希望增加步长而不是得到不稳定的结果，则需要使用四元数。在

正如@z0r在评论中解释的那样，由于四元数通过乘法变换旋转，它们之间的“区别”就是乘法逆运算——基本上是四元数除法。在

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

现在，就像对于小的theta，theta =~ sin(theta)，这个x和减法的结果差别不大，只要差别很小。像这样滥用“小角度定理”在所有类型的模拟中经常使用，但重要的是要知道何时打破它们，以及它对模型的限制。在

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

^{pr2}$

组成单位旋转

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

同样，对于小角度（即dt与速度相比很小），和和与乘积非常接近。在

如果有必要的话，像以前一样正常化。在

q_ = unit_vector(q_)

我认为这应该行得通，但会比你以前的版本慢一点，可能会有非常相似的结果。在