クォータニオンの微分方程式

クォータニオンの微分方程式
$\dfrac{d}{dt} \begin{bmatrix} q_0\\ q_1\\ q_2\\ q_3 \end{bmatrix} = \begin{bmatrix} 0& -w_x& -w_y& -w_z\\ w_x& 0& w_z& -w_y\\ w_y & -w_z & 0& w_x\\ w_z & w_y & -w_x & 0 \end{bmatrix} \begin{bmatrix} q_0\\ q_1\\ q_2\\ q_3 \end{bmatrix}$

Matlabで解析解が出るっぽい

>> syms w_x w_y w_z
>> syms q0(t) q1(t) q2(t) q3(t)
>> A=[0 -w_x -w_y -w_z;w_x 0 w_z -w_y; w_y -w_z 0 w_x; w_z w_y -w_x 0]
A =
[   0, -w_x, -w_y, -w_z]
[ w_x,    0,  w_z, -w_y]
[ w_y, -w_z,    0,  w_x]
[ w_z,  w_y, -w_x,    0]
>> odes = diff(Y)==A*Y
odes(t) =
 diff(q0(t), t) == - w_x*q1(t) - w_y*q2(t) - w_z*q3(t)
   diff(q1(t), t) == w_x*q0(t) - w_y*q3(t) + w_z*q2(t)
   diff(q2(t), t) == w_x*q3(t) + w_y*q0(t) - w_z*q1(t)
   diff(q3(t), t) == w_y*q1(t) - w_x*q2(t) + w_z*q0(t)
>> C=Y(0)==[1;0;0;0]
C =
 q0(0) == 1
 q1(0) == 0
 q2(0) == 0
 q3(0) == 0
>> [q0(t) q1(t) q2(t) q3(t)]= dsolve(odes,C)
q0(t) =
(w_y*exp(t*(- w_x^2 - w_y^2 - w_z^2)^(1/2))*((w_y*(- w_x^2 - w_y^2 - w_z^2)^(1/2))/(w_y^2 + w_z^2) + (w_x*w_z)/(w_y^2 + w_z^2)))/(2*(- w_x^2 - w_y^2 - w_z^2)^(1/2)) + (w_y*exp(-t*(- w_x^2 - w_y^2 - w_z^2)^(1/2))*((w_y*(- w_x^2 - w_y^2 - w_z^2)^(1/2))/(w_y^2 + w_z^2) - (w_x*w_z)/(w_y^2 + w_z^2)))/(2*(- w_x^2 - w_y^2 - w_z^2)^(1/2)) + (w_z*exp(t*(- w_x^2 - w_y^2 - w_z^2)^(1/2))*((w_z*(- w_x^2 - w_y^2 - w_z^2)^(1/2))/(w_y^2 + w_z^2) - (w_x*w_y)/(w_y^2 + w_z^2)))/(2*(- w_x^2 - w_y^2 - w_z^2)^(1/2)) + (w_z*exp(-t*(- w_x^2 - w_y^2 - w_z^2)^(1/2))*((w_z*(- w_x^2 - w_y^2 - w_z^2)^(1/2))/(w_y^2 + w_z^2) + (w_x*w_y)/(w_y^2 + w_z^2)))/(2*(- w_x^2 - w_y^2 - w_z^2)^(1/2))
q1(t) =
(w_z*exp(t*(- w_x^2 - w_y^2 - w_z^2)^(1/2))*((w_y*(- w_x^2 - w_y^2 - w_z^2)^(1/2))/(w_y^2 + w_z^2) + (w_x*w_z)/(w_y^2 + w_z^2)))/(2*(- w_x^2 - w_y^2 - w_z^2)^(1/2)) - (w_y*exp(-t*(- w_x^2 - w_y^2 - w_z^2)^(1/2))*((w_z*(- w_x^2 - w_y^2 - w_z^2)^(1/2))/(w_y^2 + w_z^2) + (w_x*w_y)/(w_y^2 + w_z^2)))/(2*(- w_x^2 - w_y^2 - w_z^2)^(1/2)) - (w_y*exp(t*(- w_x^2 - w_y^2 - w_z^2)^(1/2))*((w_z*(- w_x^2 - w_y^2 - w_z^2)^(1/2))/(w_y^2 + w_z^2) - (w_x*w_y)/(w_y^2 + w_z^2)))/(2*(- w_x^2 - w_y^2 - w_z^2)^(1/2)) + (w_z*exp(-t*(- w_x^2 - w_y^2 - w_z^2)^(1/2))*((w_y*(- w_x^2 - w_y^2 - w_z^2)^(1/2))/(w_y^2 + w_z^2) - (w_x*w_z)/(w_y^2 + w_z^2)))/(2*(- w_x^2 - w_y^2 - w_z^2)^(1/2))
q2(t) =
(w_y*exp(t*(- w_x^2 - w_y^2 - w_z^2)^(1/2)))/(2*(- w_x^2 - w_y^2 - w_z^2)^(1/2)) - (w_y*exp(-t*(- w_x^2 - w_y^2 - w_z^2)^(1/2)))/(2*(- w_x^2 - w_y^2 - w_z^2)^(1/2))
q3(t) =
(w_z*exp(t*(- w_x^2 - w_y^2 - w_z^2)^(1/2)))/(2*(- w_x^2 - w_y^2 - w_z^2)^(1/2)) - (w_z*exp(-t*(- w_x^2 - w_y^2 - w_z^2)^(1/2)))/(2*(- w_x^2 - w_y^2 - w_z^2)^(1/2))

たとえばq0について, 初期値q0(0)=1, q1(0)=0, q2(0)=0, q3(0)=0とすると、
f:id:seinzumtode:20200518165339p:plain
このまま使うことはなさそう