2017-06-22 60 views
2

我想創建一個簡單的物理系統在Python中的作品quaternionsvelocity/position類似的方式。它的主要目標是模擬一個被拖動的對象,並隨着時間的推移嘗試追趕另一個對象。模擬使用3個變量:k:彈簧常數,d:衰減係數,以及m:被拖動物體的質量。春季物理應用於四元數使用python

採用經典的歐拉積分,我可以解決的立場:

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 

這對於位置的偉大工程。通過改變kdm我可以得到一個迅捷/重結果和總體來說,我很滿意我的感覺:

A visual output my physics simulation

現在我想做同樣的事情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 Intration 
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 

而令我非常驚喜的是,它給了我非常不錯的結果!

quaternion version

因爲我有我的直覺去了,我相信我的解決方法是有缺陷的(例如,如果Q和Q_相向),並且有一個正確的/更好的方法來達到我想要的。

問題:

什麼是模擬在quaternions彈簧力的正確途徑,考慮到(至少)拖動對象的mass,彈簧stiffnessdamping factor

實際的Python代碼將不勝感激,因爲我讀博士論文非常困難:)。對於quaternion操作,我通常會參考Christoph Gohlke's excellent library,但請隨時在答案中使用其他人。

+1

我懷疑你在使用numpy四元數的時候比我們大多數人都領先。 – hpaulj

+0

我不知道如何正確查找'v',但'x'應該可能與[multiplicative inverse](https://stackoverflow.com/a/22167097/320036)一起找到,並且'q_'應該可能請使用[Slerp](https://en.wikipedia.org/wiki/Slerp#Quaternion_Slerp)。 – z0r

回答

0

在「模擬四元數彈簧力的正確方法是什麼」這個問題上,答案是正確地寫下潛在的能量。

  1. 四元數爲您提供了經由操作

    orientation = vector_part(q_ r q_*) 
    

    其中星號表示共軛和r是固定的方位(例如說「沿着z單位矢量」,它必須與對象的方向在所有對象的系統中都是唯一的)。 q_,r和q_ *的乘法被假定爲「四元數乘法」。

  2. 定義能量函數與這個方向上,例如

    energy = dot_product(orientation, unit_z) 
    
  3. 以「減衍生」的能量相對於四元數,你將不得不申請到系統的力量。

您當前的代碼做了這可能是一個很好的解決你的問題「在四元數空間阻尼振盪」,但它不是一個對象上的彈簧力:-)

PS:太長作爲評論,我希望它有幫助。 PS2:我沒有爲上述問題使用直接代碼,因爲(i)在上面的庫文檔中閱讀並不容易,(ii)問題的第一部分是執行數學/物理。

2

「更好」在這裏實際上會很主觀。

因爲你的步長和位移很小,所以你有點害怕四元數的概念。根據您的應用程序,這可能確實沒問題(遊戲引擎通常會利用這樣的技巧來簡化實時計算),但如果您的目標是準確的,或者您想要增加步長並且不會得到不穩定的結果,那麼您會需要使用四元數。

As @ z0r在評論中解釋過,由於四元數通過乘法變換旋轉,所以它們之間的「差異」是乘法逆 - 基本上是四元數除法。

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

現在,多爲小thetatheta =~ 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_) 

我認爲應該可以工作,但會比以前的版本慢一些,並且可能會有非常相似的結果。