シンプレクティック形式と数値積分

　ハミルトン力学の位置座標と運動量座標からは、シンプレクティック形式を自然に定義することができます。そうすると、例えば、シンプレクティック形式を保存する写像によって、正準変換やシンプレクティック積分による数値積分を記述することができます。

準備


　n次元ベクトル場 X を次のように書きます。 x^i は i 番めの座標です。 X^i は i 番めの成分です。一般には \frac{\partial}{\partial x^i} は基底を表す単なる記号であり、偏微分ではありません。しかし、変換則が分かりやすいので、この形で書かれます。

\begin{align} X=\sum_{i=1}^{n} X^i \frac{\partial}{\partial x^i} \end{align} \tag{1}

　特に、多様体 M 上では、各点ごとに接ベクトル空間を定義できます。点 p \in M での接ベクトル空間を、 T_p (M) と書きます。


　また、微分 k 形式は一般に次のように書けます。 f_{l_1, \cdots, l_k}(x_1, \cdots, x_n) は関数です。

\begin{align} u=\sum_{l_1< \cdots < l_k} {f_{l_1, \cdots, l_k}(x_1, \cdots, x_n) dx^{l_1} \wedge \cdots \wedge dx^{l_k}} \end{align} \tag{2}

　関数 f の微分 df は、ベクトル場 X を実数に移す写像 df: T_p(X)\rightarrow \mathbb{R} とみなせます。 微分形式の微分が次のようになるので、外積の反対称性から d(du)=0 が成り立ちます。

\begin{align} du=\sum_{i=1}^n \sum_{l_1< \cdots < l_k} {\frac{\partial f_{l_1, \cdots, l_k}}{\partial x^i} dx^i \wedge dx^{l_1} \wedge \cdots \wedge dx^{l_k}} \end{align} \tag{3}

　次に、開集合 U から V への可微分同相写像を \Phi とします。 (y^1, \cdots , y^n)=\Phi (x^1 ,\cdots , x^n) です。座標変換もこの写像で表せます。 \Phi によるベクトル場の変換を次のように表します。

\begin{align} \Phi_{*} X = \Phi_{*} \left( \sum_{i=1}^{n} X^i \frac{\partial}{\partial x^i} \right) = \sum_{j=1}^{n} \left( \sum_{i=1}^n {X^I \frac{\partial y^j}{\partial x^i}} \right) \frac{\partial}{\partial y^j} \end{align} \tag{4}

　逆に、 V から U への微分形式の引き戻しを、次のように表します。ここでの f は先程とは異なり、 V 上の関数です。

\begin{align} \Phi^* (u) = \Phi^* \left( \sum_{l_1< \cdots < l_k} {f_{l_1, \cdots, l_k} dy^{l_1} \wedge \cdots \wedge dy^{l_k}} \right) \\ = \sum_{l_1< \cdots < l_k} {(f_{l_1, \cdots, l_k} \circ \Phi) d\Phi^{l_1} \wedge \cdots \wedge d\Phi^{l_k}} \end{align} \tag{5}


ハミルトンベクトル場


　ハミルトン力学系の話に移ります。ハミルトニアン H によるハミルトンベクトル場は次のように定義されます。 q^i は位置座標、p^i は運動量座標です。

\begin{align} X_H = \sum_{i=1}^n {\frac{\partial H}{\partial p^i} \frac{\partial}{\partial q^i} - \frac{\partial H}{\partial q^i} \frac{\partial}{\partial p^i}} \end{align} \tag{6}

　これは単なるベクトル場ですが、形式的にリウヴィル演算子 i\mathcal{L}_Hと同じ形になっています。ポアソン括弧を使うと i\mathcal{L}_H=\{\cdot, H\} です。相空間上で定義される関数 \gamma(q,p,t) について、次の式が成り立ちます。

\begin{align} \frac{d\gamma}{dt} = i\mathcal{L}_H \gamma + \frac{\partial \gamma}{\partial t} \end{align} \tag{7}

　例えば、 \gamma =p,q の場合は、ハミルトン方程式になります。

　
　また、ハミルトン方程式の解の各々 q(t), p(t) を積分曲線と言い、関数 G(p, q) のうち、任意の積分曲線上で \frac{dG}{dt}=0 となるものを第一積分または保存量と言います。式(7) から、 G(p,q) が第一積分である必要十分条件は i\mathcal{L}_HG=0 またはポアッソン括弧を使って \{G,H \}=0 です。

シンプレクティック形式


　次の微分2形式がシンプレクティック形式です。これは閉形式（ d\omega = 0）です。

\begin{align} \omega = \sum_{I=1}^n {dp^i \wedge dq^i} \end{align} \tag{8}

　写像 \Phi: U \rightarrow V による変換 (P, Q)=\Phi (p, q) を考えます。シンプレクティック形式が変換によって不変であることと、ハミルトンベクトル場が変換によって不変であることが等価であることを確認します。このことは、次のように表せます。 \Omega は (P,Q) 空間のシンプレクティック形式です。

\begin{align} \omega = \Phi^* \Omega \longleftrightarrow \Phi_* X_{H(p,q)} = X_{H(P.Q)} \end{align} \tag{9}

　道具として、微分形式に対する次のような線形作用素 i_X を考えます。

\begin{align} i_{\frac{\partial}{\partial x^{l_j}}} (dx^{l_1} \wedge \cdots \wedge dx^{l_n}) = (-1)^{j-1} dx^{l_1} \wedge \cdots \wedge dx^{l_{j-1}}\wedge dx^{l_{j+1}} \wedge \cdots \wedge dx^{l_n} \end{align} \tag{10}
\begin{align} i_{fX+gX'}(u)=f i_{X}(u) + g i_{X'}(u) \\ i_{X}(fu+gu') = f i_{X}(u) + g i_{X}(u') \end{align} \tag{11}

i_{\frac{\partial}{\partial x^i}} は、微分形式から dx^i の項を先頭に移した後に、この項を除く操作を表しています。微分形式 u が dx^i を含まない場合には、 i_{\frac{\partial}{\partial x^i}}(u)=0 となります。


　一般に、 \Phi : U \rightarrow V が定義された時、式(4)、(5)、(10)、(11) から、次の関係が成り立つことが分かります。 u は V 上の微分形式です。

\begin{align} \Phi^* i_{\Phi_* X}(u) = i_{X}(\Phi^* u) \end{align} \tag{12}

　 \Phi が相空間 U: (p,q) から V: (P,Q) への写像であったことから、式(6)、式(8) を参照すると、

\begin{align} i_{X_H}(\omega) = dH \end{align} \tag{13}

が成り立ちます。式(12) によって、 \omega = \Phi^* \Omega から、 \Phi^* i_{\Phi_* X_{H(p,q)}}(\Omega) = i_{H(p,q)}(\omega) となります。 \Phi^* dH(P,Q) を dH(p,q) と同一視することで、式(13) と一意性から、 \Phi_* X_{H(p,q)}=X_{H(P,Q)} となっていることが分かります。


正準変換

　変換 \Phi :U \rightarrow V によってシンプレクティック形式が不変であるとき、この変換は正準変換です。シンプレクティック形式が不変ならハミルトンベクトル場も不変であり、変換後の座標でもハミルトニアンが存在します。


　正準変換によって結びつく二つのシンプレクティック形式が等しいことから、

\begin{align} \sum_i ( dp^i \wedge dq^i -dP^i \wedge dQ^i)=0 \end{align} \tag{14}

ですが、これを

\begin{align} d\left( \sum_i ( p^i \wedge dq^i -P^i \wedge dQ^i)\right) =0 \end{align} \tag{15}

と書くことで、次のような関数 S(q,Q) が存在することが分かります（ U が単連結な空間の場合）。

\begin{align} dS = \sum_i ( p^i \wedge dq^i -P^i \wedge dQ^i) \end{align} \tag{16}

　 S(q,Q) を生成関数と呼びます。生成関数を用いて下記の式が得られます。

\begin{align} p^i = \frac{\partial S}{\partial q^i} \\ P^i = -\frac{\partial S}{\partial Q^i} \end{align} \tag{17}

　式(15) の書き方に応じて、4通りの生成関数 S(q,Q), S(q,P), S(p,Q), S(p,P)を定義できます。生成関数を用いて、正準変換 (P,Q)=\Phi(p,q) を構成することができます。


シンプレクティック数値積分法


　ハミルトン方程式の解を数値積分によって近似的に求めることを考えます。リウヴィル演算子を使って

\begin{align} \frac{d}{dt}\left( \begin{array}{c} p \\ q \end{array} \right) = i\mathcal{L}\left( \begin{array}{c} p \\ q \end{array} \right) \end{align} \tag{18}

ですから、 \exp = \sum_n {\frac{1}{n!} \frac{\partial ^n}{\partial t^n}} と定義してテイラー展開と比較することで、時刻 t から t+\tau までの積分を次のような時間並進演算子として書くことができます。

\begin{align} q(t+\tau) = \exp {(i\mathcal{L}_H\tau)}q(t) \end{align} \tag{19}

p についても同様です。 \exp {(i\mathcal{L}_H\tau)} を \tau が小さい場合の近似として具体的に書き下すことで、数値積分を行うことができます。


　ハミルトニアンを分解し、 H(p,q)=K(p)+U(q) とします。リウヴィル演算子の定義から、 i\mathcal{L}_H = i\mathcal{L}_K + i\mathcal{L}_U です。一粒子系なら i\mathcal{L}_f = \frac{\partial f}{\partial p}\frac{d}{dq} - \frac{\partial f}{\partial q}\frac{d}{dp} ですから、

\begin{align} i\mathcal{L}_K q = \frac{\partial K}{\partial p} \\ i\mathcal{L}_K p = 0 \\ i\mathcal{L}_U q = 0 \\ i\mathcal{L}_U p = -\frac{\partial U}{\partial q} \end{align} \tag{20}

これらを使って、次の式が得られます。

\begin{align} \exp{(i\mathcal{L}_K\tau)} q = q + \frac{\partial K}{\partial p} \tau \\ \exp{(i\mathcal{L}_K\tau)} p = p \\ \exp{(i\mathcal{L}_U\tau)} q = q \\ \exp{(i\mathcal{L}_U\tau)} p = p - \frac{\partial U}{\partial q} \tau \end{align} \tag{21}

式(21) から、 \exp{(i\mathcal{L}_K\tau)} や \exp{(i\mathcal{L}_U\tau)} による変換がシンプレクティック形式 dp \wedge dq を不変に保つことはすぐに分かります。これらをシンプレクティック写像と言います。シンプレクティック写像の積（合成写像）もシンプレクティック写像です。そのため、次のように時間並進演算子を近似することでシンプレクティック写像が得られます。

\begin{align} \exp{(i\mathcal{L}_H\tau)} \approx \exp{(i\mathcal{L}_K\tau)}\exp{(i\mathcal{L}_U\tau)} \end{align} \tag{22}

　式(22) に従って位置と運動量を更新するのが、一次のシンプレクティック数値積分法です。式(22) 右辺の二つの項は逆に並べても構いません。具体的なアルゴリズムは、式(21) を用いて次のように簡単に実行できます。

\begin{align} q(t+\tau)=q(t) + \frac{\partial K(p(t))}{\partial p} \tau \\ p(t+\tau)=p(t) - \frac{\partial U(q(t+\tau))}{\partial q} \tau \end{align} \tag{23}

このアルゴリズムは、最も単純な数値積分法であるオイラー法と同じ程度に簡単ですが、実際のハミルトニアン H(p,q) に近い近似的なハミルトニアン \tilde{H}(p,q) に厳密に従う時間発展を記述しています。


　式(22) の右辺から決まる下記の \tilde{H} が、一次のシンプレクティック数値積分法において厳密に保存される近似的なハミルトニアンです。

\begin{align} \exp{(i\mathcal{L}_\tilde{H}\tau)} = \exp{(i\mathcal{L}_K\tau)}\exp{(i\mathcal{L}_U\tau)} \end{align} \tag{23}

　真のハミルトニアン H と 近似ハミルトニアン \tilde{H} との誤差は、時間幅 \tau に対して、シンプレクティック法の次数程度までしか拡大しません。従って、シンプレクティック法で数値積分を行う限り、エネルギー誤差が発散することはありません。これが、シンプレクティック形式を不変に保たないアルゴリズムと比較した場合の大きな利点です。


高次のシンプレクティック数値積分法と近似精度


　一次のシンプレクティック数値積分法は式(22) の近似に依拠しますが、式(22) の両辺は一致しません。ベーカー・キャンベル・ハウスドルフの公式を使って、 \tau の次数に対する両辺の差を評価することができます。

\begin{align} \exp{(A)}\exp{(B)}=\exp{\left( A+B+\frac{1}{2}[A,B] + \cdots \right) } \end{align} \tag{24}

例えば、式(22) は

\begin{align} \exp{(i\mathcal{L}_K\tau)}\exp{(i\mathcal{L}_U\tau)}=\exp{\left( i\mathcal{L}_H\tau+ O(\tau^2) \right) } \end{align} \tag{25}

なので、一次精度です。


　より一般に、シンプレクティック形式を不変に保つように \exp{( i\mathcal{L}_H\tau)} を近似する方法として、次の形が考えられます。

\begin{align} \exp{(i\mathcal{L}_H\tau)} \approx \prod_j {\exp{(c_j (i\mathcal{L}_K\tau))}\exp{(d_j (i\mathcal{L}_U\tau))}} \end{align} \tag{26}

ただし、一次近似を成立させるために \sum_j {c_j} = 1, \sum_j {d_j} = 1 である必要があります。


　例えば、次の形で近似すると、二次のシンプレクティック数値積分法（リープ・フロッグ法）が得られます。

\begin{align} \exp{(i\mathcal{L}_H\tau)} \approx \exp{(\frac{i\mathcal{L}_K\tau}{2})}\exp{(i\mathcal{L}_U\tau)}\exp{(\frac{i\mathcal{L}_K\tau}{2})} \end{align} \tag{27}

　式(24) で評価してみると、ちょうど [i\mathcal{L}_K, i\mathcal{L}_U] \tau^2 の項が打ち消しあうことで、 \exp{(i\mathcal{L}_K\tau)}\exp{(i\mathcal{L}_U\tau)} = \exp{(i\mathcal{L}_H\tau) + O(\tau^3)} となり、二次精度の近似になっていることが分かります。さらに、式(27) は演算子として右から作用させても左から作用させても同じ、対称な形をしており、時間反転対称性を満たしています。

　さらに高次のシンプレクティック数値積分法も発見されています。例えば四次のシンプレクティック法は下記のように与えられています（Forest, 1989）。 \alpha = 1 - 2^{\frac{1}{3}} として、

\begin{align} \exp{(i\mathcal{L}_H\tau)} \approx \exp{\left( \frac{\tau}{2(1+\alpha)}i\mathcal{L}_K \right) } \exp{\left( \frac{\tau}{1+\alpha}i\mathcal{L}_U\right) } \exp{\left( \frac{\alpha \tau}{2(1+\alpha)}i\mathcal{L}_K\right) } \\ \exp{\left( \frac{(\alpha -1)\tau}{1+\alpha}i\mathcal{L}_U\right) } \exp{\left( \frac{\alpha \tau}{2(1+\alpha)}i\mathcal{L}_K\tau\right) } \exp{\left( \frac{\tau}{1+\alpha}i\mathcal{L}_U\right) } \exp{\left( \frac{\tau}{2(1+\alpha)}i\mathcal{L}_K\right) } \end{align} \tag{27}

　その他、シンプレクティック積分については、散逸系のような非ハミルトン系へ適用するための Split 法や、Partitioned-ルンゲクッタ法としての表現といった話題があります。