シンプレクティック数値解法書き下してみる - なぜか数学者にはワイン好きが多い

私はここ10年ほど，各種幾何的数値解析法の普及に努めているつもりなのですが，まだまだな気がしています．
物理学の分野では，ずいぶんと使われるようになった気がします．しかし，肝心の計算機科学の分野では，まだまだです．例えば物理シミュレーションの技術が必須の日本VR学会などに行っても，名前すら知られていません．

その原因の一つは，工学系の方は公式としてオイラー法やらRunge-Kutta法を使うことは得意でも，幾何的数値解析法の場合，公式の応用として対称合成法や対称分解法などの説明になると，教科書がとたんに鈴木・吉田の公式やらベイカー・キャンベル・ハウスドルフの公式やらリー代数やら，抽象的になって理解が難しいからだと思っています．

そこで，むちゃくちゃ具体的に，対称合成法の典型例を書き下してみようと思います．
扱うハミルトニアンは，もちろん，幾何的数値解析法の基本，調和振動子です．

$\frac{d p}{dt}=-\frac{\partial H}{\partial q}=-p$
$\frac{d q}{dt}=\frac{\partial H}{\partial p}=q$

この常微分方程式を，VR学会の方々がよくやるように，1次のオイラー法で差分化してみましょう．

$p'=p-\tau q$
$q'=q+\tau p$

$\tau$ は，いわゆる「刻み幅」， $p'$ は数値計算結果による $p$ の次のステップの値とします．
この式が「1次」と言われるのは，次の理由からです．
調和振動子は解析的に解く事ができ，解は初期値を元に以下のようになります．

$p=p_0 Cos(t)-q_0 Sin(t)$
$q=p_0 Sin(t)+q_0 Cos(t)$

この正しい解を，「三角関数は難しくて良く分からない」という理由で， $t_0=0$ の周りで $t=\tau$ についてテイラー展開してみます．

$p=p_0-\tau q_0-\frac{\tau^2}{2}p+\frac{\tau^3}{6}q+O(\tau^4)$
$q=q_0+\tau q_0-\frac{\tau^2}{2}q-\frac{\tau^3}{6}p+O(\tau^4)$

上の「1次のオイラー法の差分化」と見比べると， $\tau$ の1次の項まではぴったりと合っていることが分かります．
なので，「1次」の手法です．

さて，上の「1次のオイラー法」に対して，「1次のシンプレクティック・オイラー法」は，次のようなものです．

$p'=p-\tau q$
$q'=q+\tau p'$

びみょー．．．に違いますね． $p'$ が2回出てきておりますので，中学校で習った連立方程式の手法により，右辺の $p'$ を消すことができます．

$p'=p-\tau q$
$q'=q+\tau (p-\tau q)=q+\tau p-\tau^2 q$

さらに，上の「1次のオイラー法」に対して，「1次のシンプレクティック・オイラー法の随伴スキーム」というものがあります．

$p'=p-\tau q'$
$q'=q+\tau p$

またまたビミョー．．．に違いますね．今度は $q'$ が2回出てきています．やはり連立させて消すことができます．

$p'=p-\tau (q+\tau p) = p-\tau q - \tau^2 p$
$q'=q+\tau p$

ここで，シンプレクティック・オイラー法と，シンプレクティック・オイラー法の随伴スキームを合成させることを考えます．
これが対称合成法なのですが，次のステップを表す記号「'」がごっちゃになっておりますので，随伴スキームの方は「''」に置き換えることにします．ステップを一つずらすことになるので，p→p'，p'→p''のようにずらしてみます．
すると，連立させて重複変数を消した結果の式は，次のようになります．

$p''=p'-\tau (q'+\tau p') = p'-\tau q' - \tau^2 p'$
$q''=q'+\tau p'$

さて，このシンプレクティック・オイラー法の随伴スキームを，元のシンプレクティック・オイラー法のスキームと連立させて， $p'$ と $q'$ を消してみましょう．

$p''=p'-\tau (q'+\tau p') = p'-\tau q' - \tau^2 p'=(p-\tau q)-\tau (q+\tau p-\tau^2 q)-\tau^2 (p-\tau q)=p-2 \tau q-2 \tau^2 p + 2 \tau^3 q$
$q''=q'+\tau p' = q+\tau p-\tau^2 q + \tau (p - \tau q) = q + 2 \tau p - 2 \tau^2 q$

だいぶ式がややこしくなりましたが，これで最後ですので頑張って下さい．
最後に， $\tau$ を $1/2 \tau$ にします．これがベイカー・キャンベル・ハウスドルフ・ポアンカレの公式です．

ホントにややこしく見えてしまうので書きませんが，書き下してみて下さい．見事に解析解の $\tau^2$ の項までは一致することが分かると思います．
よって，2種類の1次の公式を組み合わせると，2次の公式が出来てしまいました．
ただのオイラー法よりも精度が良い上にシンプレクティック性という物理学的な保存量が厳密に保存されます．

次は，TCI(全保存型差分法）の具体例を紹介してみたいと思っています．