全エネルギーを保存する時間積分法 -

全保存量を保つコマの離散化
http://www.riam.kyushu-u.ac.jp/fluid/meeting/16ME-S1/papers/Article_No_27.pdf

Stackel系の全ての保存量を保つ離散化
http://ci.nii.ac.jp/naid/110000167002

離散Kepler 運動の時間補正
http://www.riam.kyushu-u.ac.jp/fluid/meeting/16ME-S1/papers/Article_No_36.pdf

正則化法と全保存型差分法を用いた重力N体問題のシミュレーション
https://qir.kyushu-u.ac.jp/dspace/bitstream/2324/14286/1/Article_No_14.pdf

可積分な相対論的運動方程式に対する全保存型差分法とその時間対称化
http://ci.nii.ac.jp/Detail/detail.do?LOCALID=ART0008214646&lang=en

あたりを読んだ。いずれもTotally Conservative Integrator(全保存型差分法,TCI)の話。Hamilton力学系の全ての保存量を保存する差分法という、凄そうな話ではあるけれども使える系が限定されてるので実はそうでもない。で、上の文献読む限りでは、エネルギーを保存する差分法(Greenspanのスキーム)というものに、うまく帰着させるというのがミソなのだけど、肝心のエネルギー保存差分スキームが詳しく説明されておらず、70年代の古い文献を引いているだけでよく分からないので勝手に想像してみた。というわけで、以下のはわたしの想像であり本来のGreenspanのスキームとは違うかもしれません(エネルギーは保存すると思うけど)。ついでに、実はよく知られた方法で有名だったりするのかもしれませんが、多分そんなことはないんじゃないかと思います。

簡単のために、一番単純なPoisson構造
$\{q_i,p_j\}=\delta_{ij},\{q_i,q_j\}=\{p_i,p_j\}=0$
のケースを考える。まず通常のHamilton力学系で、エネルギーが保存する理由を思い出す。Poisson括弧が反対称だからというのが一番簡潔な答えだけども、もうちょっと遡って考えると、
$\frac{dH}{dt} = \sum(\frac{dp_i}{dt}\frac{\partial H}{\partial p_i}+\frac{dq_i}{dt}\frac{\partial H}{\partial q_i})$
$(\frac{dp_i}{dt},\frac{dq_i}{dt}) = (-\frac{\partial H}{\partial q_i},\frac{\partial H}{\partial p_i})$
という2つの関係式から、エネルギー保存が出る。前者は、単なる合成関数の微分。後者は、正準方程式。そして、第二式を使えば、第一式が0になることは明らか。

上の式を見ていれば、空間微分も差分化してみようという気分になる。空間微分を何らかの決まった仕方で差分化して、正準方程式を要請してやると、第一式の右辺は必ず0になる。つまり、エネルギー保存を妨げる障害となるのは、第一式自体が成立しないということが根源といえる。ところで、まだ時間差分の仕方も決めてなかったけど、時間差分は通常のEuler法の場合と同様に考える。つまり、物理量Fに関する時間微分を
$\100dpi \frac{F(\bold{p}_{i+1},\bold{q}_{i+1})-F(\bold{p}_{i},\bold{q}_{i})}{t_{i+1}-t_i}$
のように差分化する。このとき、考えるべき問題は、以下のように集約される。

問題：Euler法と同じ時間差分を採用した場合、「合成関数の微分に関する公式」を差分化後にも成立させるためには、どのような空間差分を導入すればよいか?

これに対する答えは例えば以下の通り。書くのが面倒なので、相空間が2次元の場合に限定して書くと、
$\100dpi \frac{H(p_{i+1},q_{i+1})-H(p_{i},q_{i})}{t_{i+1}-t_i}=\frac{p_{i+1}-p_i}{t_{i+1}-t_i}\frac{H(p_{i+1},q_{i+1})-H(p_i,q_{i+1})}{p_{i+1}-p_i}+\frac{q_{i+1}-q_i}{t_{i+1}-t_i}\frac{H(p_i,q_{i+1})-H(p_i,q_i)}{q_{i+1}-q_i}$
となる。この式が正しいことは、中学生でも分かる。厳密に言えば、pやqが時間によらず一定の時などは、0除算が発生するけれど、この空間差分の分母は分子を因数に持つので、そのような場合でも式は意味を持つ。空間差分幅は一定でなく、時間発展によって決まる。ついでに、普通、時間差分幅の取り方は一定にするけれども、別に一定にしなくても、エネルギー保存は成り立つ。あと、よく見ると、pに関する差分とqに関する差分がびみょ〜に違うことが分かる。pとqの差分化が対称でないので、
$\100dpi \frac{H(p_{i+1},q_{i+1})-H(p_{i},q_{i})}{t_{i+1}-t_i}=\frac{p_{i+1}-p_i}{t_{i+1}-t_i}\frac{H(p_{i+1},q_i)-H(p_i,q_i)}{p_{i+1}-p_i}+\frac{q_{i+1}-q_i}{t_{i+1}-t_i}\frac{H(p_{i+1},q_{i+1})-H(p_{i+1},q_i)}{q_{i+1}-q_i}$
という差分化も可能になる。相空間の次元があがった時の差分化の仕方は、明らか。

例として、一次元調和振動子を考える。つまり、
$\100dpi \{q,p\}=1$
$\100dpi H=\frac{1}{2m}p^2 + \frac{k}{2}q^2$
というPoisson構造とhamiltonianを考える。得られる正準方程式を上の手続きに基づいて差分化すると、空間差分の仕方によらず
$\100dpi q_{n+1}-q_n = \frac{\Delta t}{2m}(p_{n+1}+p_n)$
$\100dpi p_{n+1}-p_n = \frac{-k \Delta t}{2}(q_{n+1}+q_n)$
となる。面倒なので、時間差分幅は一定とした。陽的な形に整理すると、
$\100dpi p_{n+1} = \frac{1}{4m+k(\Delta t)^2}\{(4m-k(\Delta t)^2)p_n - 4\Delta t km q_n\}$
$\100dpi q_{n+1} = \frac{1}{4m+k(\Delta t)^2}\{4\Delta t p_n + (4m-k(\Delta t)^2)q_n\}$
となる。これは線形で、行列式を計算すると1になる。つまり、面積保存写像になっていて、symplectic。多分、こんなことは一般には成り立たないと思う。

上で述べたようなことは、もっと一般のPoisson構造、例えば相空間が奇数次元の場合とかSymplecticでないPoisson構造に対しても使える。Poisson構造が、
$\100dpi \{f,g\} = \sum J_{ij} \frac{\partial f}{\partial x_i}\frac{\partial g}{\partial x_j}$
と書けるとすると、
$\100dpi \frac{dH}{dt} = \sum \frac{dx_i}{dt}\frac{\partial H}{\partial x_i} = \sum J_{ij}\frac{\partial H}{\partial x_i}\frac{\partial H}{\partial x_j}$
が成り立つ。最初の等号は合成関数の微分則によって、二番目の等号は正準方程式による。そして、Jの反対称性(=Poisson括弧の反対称性)から、これは0となる。なので、先ほどと同じ手続きで、空間微分を差分化してやれば、エネルギーは保存する。

SymplecticでないPisson構造とHamiltonianの例として、Eulerのコマを考える。
$\100dpi \{X,Y\}=Z,\{Y,Z\}=X,\{Z,X\}=Y$
$\100dpi H=\frac{1}{2}(\frac{X^2}{A}+\frac{Y^2}{B}+\frac{Z^2}{C})$
がPoisson構造とHamiltonian。(A,B,C)は慣性モーメントで定数。上に書いた手続きで差分化すると、
$\100dpi \frac{X_{n+1}-X_n}{\Delta t} = \frac{Z_n}{B}(Y_n+Y_{n+1}) - \frac{Y_n}{C}(Z_n+Z_{n+1})$
$\100dpi \frac{Y_{n+1}-Y_n}{\Delta t} = \frac{X_n}{C}(Z_n+Z_{n+1}) - \frac{Z_n}{A}(X_n+X_{n+1})$
$\100dpi \frac{Z_{n+1}-Z_n}{\Delta t} = \frac{Y_n}{A}(X_n+X_{n+1}) - \frac{X_n}{B}(Y_n+Y_{n+1})$
となる。エネルギーが保存することは簡単に分かる。でもって若干複雑ではあるけれども、陽的に解くこともできる。

Eulerのコマには、もう一つ全角運動量という保存量があるけれども、差分化すると、これは保存しなくなる。これは単に
$\100dpi L=X^2+Y^2+Z^2$
というもの。式をじっと見ていると、
$\100dpi \frac{X_{n+1}-X_n}{\Delta t} = \frac{Z_n+Z_{n+1}}{B}(Y_n+Y_{n+1}) - \frac{Y_n+Y_{n+1}}{C}(Z_n+Z_{n+1})$
という形にしてやると(残りの2式は面倒なので省略)、LもHも保存することが分かる。空間微分を差分化しても、Poisson構造を決める係数Jについては何もしなかったけど、これについても修正するとよいのかもしれない(どういう原理に基づいて修正すればよいかは不明。例えばJacobi恒等式とか?)。ただまあ、こいつは、もはや陽的には解けない。一応、この式は、

Integrable discretizations of the Euler top
http://arxiv.org/abs/solv-int/9803016

で同じものが導かれているけれども、Eulerのコマの差分化としては、広田・木村による差分化の方がメジャーっぽい。こっちは陽的に解けるからだと思う。ただ、広田・木村による離散Euler topは、エネルギーと角運動量に対応する保存量は存在するものの、式の形が変わってしまってるので、数値計算とかに使うには微妙だし、導出の仕方も機械的に出てくるようなものでもない。

おまけ。世の中で普及している時間積分法と言えば、Symplectic時間積分法である(というのは若干嘘で、普及している実用的なアプリで4次のSymplectic積分とか使われてるのを見たことない。Verletでみんな満足らしい)。Symplectic積分は、誰が言い出したのかは知らないが、要するに、時間発展は正準変換であるべし！という積分法で、生のHamiltonianは保存しないけど、影のHamiltonianとか呼ばれる代替物は保存する。まあでも、厳密にHamiltonianを保存してほしいというのは自然な要請。

以下では、H=T+V型のhamiltonianを考える。Tが運動量、Vがポテンシャル。普通に物理で出てくるHamiltonianはこういう形をしている。でもって、物理量fに対して一階の線形微分作用素 $\100dpi D_f(g)=\{g,f\}$ を定義する。すると物理量fの時間発展は、形式的に $\100dpi exp(t D_H)f=exp(t(D_T+D_V))f$ で書ける。そして、これは
$\100dpi exp(t(D_T+D_V))=exp(t D_T)exp(t D_V) + O(t^2)$
$\100dpi exp(t(D_T+D_V))=exp(\frac{1}{2}t D_T)exp(t D_V)exp(\frac{1}{2}t D_T) + O(t^3)$
といった具合で、tの2次、3次、...のオーダーで近似できる。 $\100dpi exp(t D_T),exp(t D_V)$ は、それぞれTとVによる時間発展なので、例えば2次の場合は、時間tだけVをhamiltonianとする時間発展してその後時間tだけTをHamiltonianとする時間発展を合成すればよい。TやVによる時間発展は厳密に解けて、厳密に解けるってことは、正確に正準変換で解の時間発展を記述できる。というのが、Symplectic積分の要諦。

結局、Symplectic数値積分のやってることは、Hによる時間発展を厳密に解ける時間発展に分解しているだけ。さて、こいつがHを保存しないということは、Tによる時間発展やVによる時間発展がHを保存しないということ。それはまあ当然で、TとVは一般にPoisson可換ではないので、TやVによる時間発展がHを保存するわけがない。

そんなわけで、標準的なSymplectic積分法とエネルギー保存を両立させるのは、相当無理筋のように見える。多分、全然違う発想で、Symplecticな時間積分法を考えてやる必要があるんじゃないかと思う。