角速度ベクトル、回転時間が与えられた時の回転行列

対角化が面倒そうだったので、加法定理を求める方針でやってみる事にした。
つまり、

(I+δt(x Sx+y Sy+z Sz))(a0(t)I+a1(t)Sx+a2(t)Sy+a3(t)Sz)
=a0(t+δt)I+a1(t+δt)Sx+a2(t+δt)Sy+a3(t+δt)Sz ... (1)

を満たし、かつa0(0)=1, a1(0)=a2(0)=a3(0)=0であるような関数(連続性とか微分可能性とかの細かな扱いは知らない。必要なら行列の指数関数の話と併用すれば論理の隙間は埋められるだろう)を求めれば良い。それが求まるとa0(t)I+a1(t)Sx+a2(t)Sy+a3(t)Szが、角速度ベクトル(x,y,z), つまり、回転軸(x,y,z)/|(x,y,z)|, 角速度|(x,y,z)|でtだけ回転した時の回転行列、つまり、回転軸(x,y,z)/|(x,y,z)|, 回転角|(x,y,z)|tの回転行列を表す。

(x,y,z)をベクトルωで表し、(a1(t),a2(t),a3(t))をベクトルα(t)で表す事にする(ベクトルを明示するうまい方法がないのでここでだけギリシャ文字はベクトルを表すことにする)。しばしばtをオミットする。

すると、(1)式から
da0/dt=-ω・α
dα/dt=a0ω+ω×α
が得られる。
左辺にそれぞれa0, a1, a2, a3を書けて足しあげると0となる。これはa0^2+a1^2+a2^2+a3^2が不変である事を表す。t=0ではa0=1で他は0だから、連続性も考慮して、
a0=sqrt(1-|α|^2)
と書ける事が分かる。後はαの各成分だけ求めれば良い。αの微分方程式は

dα/dt=sqrt(1-|α|^2)ω+ω×α ... (2)

であるから、左からωを外積すると右辺第一項が消えて
d(ω×α)/dt=ω×(ω×α)
を得る。右辺はω×αに直交するから、|ω×α|は保存する(まじめにやるならば両辺とω×αの内積を取って右辺が消滅する事をがんばって計算すればd(|ω×α|^2)/dt=0を示せる。超めんどい。)。その値はt=0を参照すれば、0である。とすれば、ω×αは零ベクトルである。(2)式はいまや

dα/dt=sqrt(1-|α|^2)ω ... (3)

となる。αの各成分は初期値が等しく、微分係数が定数にωの各成分を乗じた形になっているので、αは全ての時刻でωに比例する。α=f(t)ω, とおくとf(t)が満たす微分方程式は
f(0)=0,
df/dt=sqrt(1-f^2|ω|^2)
である。必要ならばさらにg(t)=f(t)|ω|と変数変換しても良いが、ともかく、解が
f(t)=|ω|^(-1)sin|ω|t
と求まる。結局、

exp(t(x Sx+y Sy+z Sz))=cos(sqrt(x^2+y^2+z^2)t)
   +sin(sqrt(x^2+y^2+z^2)t)(x Sx+y Sy+z Sz)/sqrt(x^2+y^2+z^2)

と求まる。特に、sqrt(x^2+y^2+z^2)=1ならばexp(t(x Sx+y Sy+z Sz))は軸(x,y,z)周りの角度tだけの回転行列を表すが、それは、

exp(t(x Sx+y Sy+z Sz))=cos(t)I
   +sin(t)(x Sx+y Sy+z Sz)

である。

で、記憶によればこの結果は間違っていて、cos(t/2), sin(t/2)みたいになるはずなんだけど・・・結局、そのうち対角化の方もやってみなきゃならなそう。あと、そろそろWPFのQuaternionクラスのコンストラクタに色々なパラメタを入れたり、それからRotateTransform3Dを作って行列を調べたりして、ここまでの話と照らし合わせて見よう。あっちは、右手系だからちゃんと整合せいると思うんだけど・・・。Quaternionのコンストラクタには軸と角度を与えるコンストラクタがあるのだが、正規化されていないベクトルを与えた時の挙動が少し気になる。場合によっては、これらの引数は角速度ベクトルと回転時間としての役割を持つかもしれないし、内部で強制的に正規化されるかもしれない。この辺り全くドキュメント化されていない。