4元数の計算ルールに親しもう。
このページはマス旅の一部です。
ハミルトンさんの人物像や4元数の歴史などの周辺はすっとばして、
初回である今回は、計算中心で、
複素数の拡張版としての4元数の特徴をつかんでいきましょう。
１．4元数をベクトル・複素数につなげよう
＜ハミルトン流のベクトル積＝外積ー内積＞
R_3のベクトルの正規直交基底ベクトルをi,j,kとすると、
i^2=j^2=k^2=1,i・j=j・k=k・i=0であり、
コトバの意味から自明な結果だね。
一方で、
4元数が作る空間の4基底ベクトルをi,j,kとすると、
i^2=j^2=k^2=-1,
ij=-ji=k,jk=-kj=i,ki=-ik=j, ijk=-1
と基底の積については
交代性とサイクリック性があるという奇怪な定義がされる。
なぜだろうか？
これは、ベクトルの積についてのハミルトン流の定義があることから出てくる。
ハミルトンの4元数の基底ベクトルの積は外積ー内積で定義されている。
通常の内積a・b=|a||b|cosθから、a・a= |a|^2 、θ=π/2のとき、a・b=0
外積＝a×b(a,bに垂直なベクトルで、aからbに向かう右ねじの位置にある。）
外積は交代性がある。
ij=i×j -i・j =i×jから、直交基底どうしの積は外積と同じ性質を持ち、右手系なら、ij=i×j=kとなる。
また、i^2=0- |i|^2=-1となる。さらに、ijk=kk=-1もすぐわかる。
これで、R_3の正規直交基底ベクトルについて、ベクトルのハミルトン流の積＝外積ー内積
というルールの追加をするだけで、
4元数の基底ベクトルの性質がすべて導けたね。
＜4元数の演算と複素数の演算に似ている＞
4元数＝スカラー部＋R_3ベクトル部
という定義をする。
スカラー部をq0, ベクトル部qvの基底をi,j,kとしたときの成分が(q1,q2,q3)ならば、
4元数q=(q0,q1,q2,q3)
と４つの実数で表記できるから、4元数という。
「表記」について、
クォータ二オンはチルダを上につけて表すが、誤解がなければ何もつけない。
クォータというと4分の1というイメージが強いのでただ要素が４つあるという感じで
4元数という名称を使うことにする。
対象が4元数とわかるときはそれすら書かない。
q=q0+qv
や（q0,qv)とかくと面白い。4元数を2要素のようにかけるからだ。
複素数は共役は虚部の符号を反対にしたように、
共役*はベクトル部の符号を反転する。
複素数のノルムが共役との積の平方根だったように、4元数のノルムも同じだ。
複素数の逆数が共役をノルム2乗で割ったように、4元数の逆数も同じだ。
複素数の和差の実部は実部の和差、和差の虚部は虚部の和差でした。
4元数の場合も同様で、和差は実部どうし、ベクトルどうしになります。
すっとばしてしまったが、大切な4元数の積はどうだろう。
p=p0+pv,q=q0+qvの積だね。積＝外積ー内積
pq=(p0+pv)(q0+qv)=p0q0+q0pv+p0qv+pvqv
=p0q0+q0pv+p0qv+pv×qv-pv・qv
ここで、pv,qvの成分p1,p2,p3,q1,q2,q3と基底i,j,kにまで落とし込むと複雑な公式ができる。
それよりも、外積はベクトルで、内積は実数になることを考えて、実部とベクトル部に
わけてみましょう。
(p0+pv)(q0+qv)＝(p0q0-pv・qv)+(q0pv+p0qv+pv×qv)
とかけば、きれいに見えますね。
２．4元数を行列につなげよう
＜回転を行列にする＞
オイラー等式e^{ix}= cosｘ＋isinｘを使うと、
ｎを単位ベクトルにして、e^{nθ/2}でθ/2回転を表す4元数が作れる。
θ回転の4元数をq=cosθ/2+n sinθ/2＝q0+qv
とする。
なぜθ回転なのに、θ/2を入れるかというと、qとq*をセットにして使うからだ。
n^2が、i^2のときのように-1になり、qの共役はベクトル部だけ符号が反転することから、qのノルムは１になる。
位置ベクトルのθ回転はr'=qrq*で、
座標系のθ回転はr'=q*rqで
実現できる。
r・nの値をrnとすると、
位置ベクトルのθ回転でr'=qrq*=(cosθ/2+n sinθ/2)r(cosθ/2-n sinθ/2)
=(rn)n+(r-(rn)n)cosθ＋n×r sinθ
n=(0,0,1)とr=(x,y,z)とすれば、行列R=
{{cos θ, -sinθ, 0}, {sinθ, cosθ,0 }, { 0,0,1}}によって、r'=Rr
r'=(q0+qv)r(q0-qv)を展開して、qの成分でRを表すこともできるね。
同様に座標系の回転でr'=q*rq=(cosθ/2-n sinθ/2)r(cosθ/2+n sinθ/2)
=(rn)n+(r-(rn)n)cosθ＋r×n sinθ
n=(0,0,1)としてｚ軸回りの行列Sは
{{cos θ, sinθ, 0}, {-sinθ, cosθ,0 }, { 0,0,1}}によって、r'=S r
r'=(q0-qv)r(q0+qv)を展開して、qの成分でSを表すこともできるね。
RとSの関係は視点の変更でしかないから、転置行列になる。tR=S
しかも、RとSは直交行列なので、RS=I（単位行列）
また、どちらももとのベクトルに対して行列をかけるだけだから、
連続的な回転は、行列の積で実現できる。
具体的な設定で３D回転を4元数でやってよう。
#位置ベクトルkをj（ｙ軸）を軸にπ/2回転する場合
n=(0,j,0),r=(x,y,z)
q=cosπ/4+nsinπ/4=1/√2(1+j)
k'=1/√2(1+j) k 1/√2(1-j)=1/2(1+j)(k-kj)=1/2(k+jk-kj-jkj)
=1/2(k+jk+jk-k)=jk=i
k(z軸上）が90度回転して、i(ｘ軸上に移動する）
#位置ベクトル3j(y軸上）を直線y=xを軸にπ回転する場合
n=(1/√2,1/√2,0),r=(x,y,z)
q=cosπ/2+1/√2(i+j)sinπ/2=1/√2(i+j)
(3j)'=1/√2(i+j) 3j 1/√2(-i-j)=-3/2(i+j)(ji+jj)
=-3/2(i+j)(-k-1)=3/2(ik + i +jk +j)=-3/2(ki+kj-i-j)
=-3/2(j-i-i-j)=3i(x軸上）
基底ベクトルの積のルールから、回転の行先がわかるね。
予想どうりの結果になりましたね。
＜回転の合成＞
剛体Pの進行方向の方向ベクトルをp=[[px],[py],[pz]]とする。
z軸を中心にPをγ（ヨー角）回転するとき、cosをc、sinをsとかくと、
Rz={{cγ, -sγ, 0}, {sγ, cγ,0 }, { 0,0,1}}だった。
同様にして、
y軸を中心にPをβ（ピッチ角）回転するとき、
Ry={{cβ,0, sβ}, {0, 1,0 }, { -sβ,0,cβ}}となり、
x軸(ベクトルp)をα（ロール角）回転して姿勢変化するとき、
Rx={{1,0,0}, {0,cα,-sα }, { 0,sα,cα}}となる。
回転の合成によって、
G=Rz◦Ry◦Rx=
[[cγcβ, cγsβsα-sγcα, cγsβcα+sγsα],
[sγcβ, sγsβsα+cγcα, sγsβcα-cγsα],
[-sβ, cβsα , cβcα ]]
となる。
また、4元数のまま積を求めることもできる。
成分で表すと
rz=(c(γ/2),0, 0, s(γ/2)), ry=(c(β/2),0,s(β/2),0),rx=(c(α/2),s(α/2),0,0)
G=(c(γ/2),0, 0, s(γ/2))(c(β/2),0,s(β/2),0)(c(α/2),s(α/2),0,0)
= (
c(γ/2) c(β/2) c(α/2 ) + s(γ/2)s(β/2)s(α/2),
c(γ/2)c(β/2)s(α/2)-s(γ/2)s(β/2)c(α/2),
c(γ/2)s(β/2)c(α/2)+s(γ/2)c(β/2)s(α/2),
s(γ/2)c(β/2)c(α/2)-c(γ/2)s(β/2)s(α/2))
4元数を成分として表した方が、cとsの双対ペア表現が見られるね。
3次行列SO3（特殊直交3次行列）の方が実用的ではあります。
どちらで計算しても、位置ベクトルｐは
p = G (0, 1, 0, 0) G∗ = (0, cγcβ, sγcβ,−sβ)
つまり、pの成分にはαが出てきません。
当然ですが、剛体Pの向きは姿勢変更角αの影響を受けてないことがわかりますね。
３．振り返り
＜振り返り＞
こんなに4元数の計算がうまくできた理由は何だろう。
振り返ってみよう。
4元数の積が外積ー内積で定義されるという、単純明快なルールがあったから。
３つの基底i,j,kの積の交代性と2乗がー1になることはそのルールから導けたから。
共役がベクトル部の符号反転であることは複素数と似ていて違和感なく使えたから。
なによりも大きかったのは、
単位ベクトルnを軸にθ回転がq=e^{n θ/2}=cos θ/2 + n sin θ/2とかけることだった。
しかし、肝心のこの事実は天下りで使い、検証すらしていなかった。
r'=qrq*で2回かけているからθ/2でいいだろうくらいにしていた。
軸にする単位ベクトルをn=(ux,uy,uz)とするとき、
4元数q=(cosθ/2,ux sinθ/2,uy sinθ/2,uz sin θ/2)=(q0,q1,q2,q3)とおけば、
4元数r=(0,x,y,z)として、r'=qrq*を展開整理するとロドリゲスのθ回転の公式と
一致する。ロドリゲスの公式自体は、ベクトルの外積と内積から導くことができる。
課題：4元数成分の数式が回転後の３D座標になることを確認しよう。
タイトルは「目でみるロドリゲスの公式」
θ = Slider(0, 2π, 0.01) #回転角（これだけアニメーションにします。）
α = Slider(0, 2π, 0.01) #回転軸の傾き
β = Slider(-π/2, π/2, 0.01)
n = (cos(β)*cos(α), cos(β)*sin(α), sin(β)) #単位回転軸ベクトルn
R = (1, 2, 3) #回転させたい点
Vector((0,0,0), R) #位置ベクトルR
#4元数成分(q0, q1, q2, q3)
q0 = cos(θ / 2)
q1 = x(n) * sin(θ / 2)
q2 = y(n) * sin(θ / 2)
q3 = z(n) * sin(θ / 2)
x_p = (q0^2 + q1^2 - q2^2 - q3^2)*x(R) + 2*(q1*q2 - q0*q3)*y(R) + 2*(q1*q3 + q0*q2)*z(R)
y_p = 2*(q1*q2 + q0*q3)*x(R) + (q0^2 - q1^2 + q2^2 - q3^2)*y(R) + 2*(q2*q3 - q0*q1)*z(R)
z_p = 2*(q1*q3 + q0*q2)*x(R) + 2*(q2*q3 + q0*q1)*y(R) + (q0^2 - q1^2 - q2^2 + q3^2)*z(R)
R' = (x_p, y_p, z_p)
Vector((0,0,0), R') #回転後のベクトルR'
Circle((R n) n, Distance(R, (R n) n), n)#R'の軌道円（設定は点線）
見た目は、点がきれいに動くだけのことだけれども、
やっている内容はすごいですね。
4元成分を関数式扱いにして、そのまま３Dの座標にうまく入れ込んでいるのですからね。
詳しい計算は省きますが、4元数の積qrq^* を代数的に計算すると
たとえば、iを軸に回転したとき、r=xi+yj+zkをサンドイッチして、線形にわけると、
xiのサンドは何もしないと同じになるのですが、yj+zkのサンドは左からのかけ算と右からのかけ算が同じ方向として重なりあい、θ回転を生みます。
「挟み込むことで、軸方向を変えずに垂直面だけを綺麗に2倍（θ/2 + θ/2 = θ）回転させるフィルターとして機能している」という事実が、ロドリゲス公式との一致から一目瞭然になりますね。