解析力学の幾何学的側面

これはU-TOKYO AP Advent Calendar 2017の23日目の記事です。ふとみたら空いていたので、一応関係者ということで埋めてみます。

はじめに

力学とは、極論すれば運動方程式を解く学問です。では、運動方程式とはなんでしょうか？今、簡単のために一次元系を考えることにします。位置を$r$、速度を$v$とすると、運動方程式とは、それぞれの時間微分を現在の$r$や$v$で表現したものです¹。

\begin{align}
\frac{dr}{dt} &= f_r(r,v) \\ 
\frac{dv}{dt} &= f_v(r,v) 
\end{align}

時間に依存する変数の組について、時間微分がその変数で記述できている時、それを力学系と呼びます。従って運動方程式とは力学系の一種です。

つまり、運動方程式とは、系を記述する変数で張られた空間の点それぞれに、その時間微分を表すベクトルを対応付ける写像であると考えることができます。この写像は複数の関数から構成されますが、その関数は自由度の数だけ必要になります。今見たように、一次元一自由度系なら$f_r$と$f_v$と、二本の関数が必要になりますし、三次元の$3N$自由度系なら$6N$本の関数が必要になります。

しかし、解析力学の講義で学んだように、世の中の運動は単一のスカラー関数で支配されています。世の中にはただひとつのラグランジアンがあり、変分原理から必要な数だけの関数を作ることができます。さらに、ラグランジアンよりも、ハミルトニアンで考えた方がいろいろ見通しがよくなることも学んだでしょう。以下では、ルジャンドル変換により、　変数を$(r, v)$から$(q, p)$にとり、この運動が時間非依存なハミルトニアン$H(p,q)$に支配されていると考えましょう。

ハミルトンの運動方程式はこうかけます。

\begin{align}
\frac{dp}{dt} &= - \frac{\partial H}{\partial q} \\ 
\frac{dq}{dt} &=  \frac{\partial H}{\partial p}
\end{align}

たったひとつのスカラー関数$H$で運動が記述できています。多自由度系でも同様で、ハミルトニアンを自分に共役な変数で偏微分すれば、自分の時間微分が出てきます。

さて、ハミルトンの運動方程式は非常に対称性が良い形をしているため、いくつかすぐにわかる性質があります。よく知られているのは、ハミルトニアン$H$が時間不変量になることでしょう。

\begin{align}
\dot{H} &= \frac{\partial H}{\partial q} \dot{q} + \frac{\partial H}{\partial p} \dot{p} \\
&= \frac{\partial H}{\partial q} \frac{\partial H}{\partial p}
- \frac{\partial H}{\partial p}\frac{\partial H}{\partial q} \\
 &= 0
\end{align}

実はもう一つ重要な性質があります。それは、以下の等式が成り立つことです。

\frac{\partial \dot{q}}{\partial q} + \frac{\partial \dot{p}}{\partial p} = 0

多自由度系でも同様な式が成り立ちます。

\sum_i \frac{\partial \dot{q_i}}{\partial q_i} + \frac{\partial \dot{p_i}}{\partial p_i} = 0

既に見たように、運動方程式とは、位相空間にベクトル場を導入する写像です。このベクトル場は、位相空間における速度場だと思うことができます。速度場が定義された系に、ひとつトレーサーを置いてみましょう。速度場に添って流れていきます。この流線が運動方程式の軌跡となります。

さて、先程の等式は、速度場のベクトルを、各成分で偏微分して和をとったものです。これは実はベクトル場のdivergenceを計算していることになります。$(\dot{p}, \dot{q}) = \vec{z}$と表記すると、先程の式は

\frac{\partial \dot{q}}{\partial q} + \frac{\partial \dot{p}}{\partial p} = \mathrm{div} \vec{z} = 0

とかけます。これは、ハミルトンの運動方程式が位相空間に引き起こす流れ場のdivergenceがゼロ、すなわち、この流れが非圧縮であることを意味します。

流れ場が非圧縮である場合、流れの微小体積要素は時間発展により体積が保存されます。今、時刻$t$において$(p,q)$にいた系が、この運動方程式により時間発展した結果、時刻$t+h$に$(P,Q)$に移動したとすると、微小体積要素が保存するとは

dp dq = dP dQ

が成り立つことです。これは、$(p,q)\rightarrow (P,Q)$という写像が、シンプレクティック写像であることを意味します。

ある初期条件の集団の運動を考えます。この運動方程式により時間発展させた場合、その集団の形はどんどん変形していきますが、その体積は変わりません。従って、その密度(確率密度)も一定のままです。これは、この運動によりエントロピーが変化しないということを意味します²。

シンプレクティック変換とシンプレクティック積分

さて、ハミルトンの運動方程式による時間発展がシンプレクティック変換であることをもう少し詳しく見てみましょう。

ある変換により$(p,q)$を$(P,Q)$に変換するとは、$(P,Q)$を$(p,q)$の関数として表現するということです。この時、

dP dQ = \left|
\begin{matrix}
\partial_p P & \partial_p Q \\
\partial_q P & \partial_q Q 
\end{matrix}
\right| dp dq

ですから、この変換がシンプレクティック変換であるとは、ヤコビアン$\partial(P,Q)/\partial(p,q)$が1であることを意味します。

さて、いま運動方程式に従って系を時間発展させることを考えましょう。簡単な系を除いて、運動方程式は厳密には解けません。そこで、微小な時間ステップ$h$だけ時間が進んだ世界を近似的に計算し、それを必要なだけ繰り返すことにします。この計算手法を分子動力学法(Molecular dynamics method, MD)と言います。

簡単のため、以下の一次元調和振動子を考えましょう。

\begin{align}
\dot{p} &= -q \\
\dot{q} &= p 
\end{align}

「現在の時刻から時間$h$だけ進んだ世界を近似する」最も単純な手法は、一次のオイラー法でしょう。

\begin{align}
P &= p - q h \\
Q &= q + q h 
\end{align}

これは、例えば

\dot{p} = -q 

という式を$t$から$t+h$まで時間積分した式

P\equiv p(t+h) = p(t) + \int_t^{t+h}\dot{p} dt

において、$q$を定数だと思って積分を実行したことに対応します。

さて、よく知られているように、一次のオイラー法で数値積分するとあっというまにエネルギーが発散します。

euler.rb

h = 0.01
p = 1.0
q = 0.0
10000.times do
  p_old = p
  q_old = q
  p = p - q_old * h
  q = q + p_old * h
  puts "#{p} #{q}"
end

本来、円を描かなければならない起動が螺旋を描きながら大きくなっていることがわかります。ルンゲクッタや、予測子修正子法も普通に使うとエネルギーが保存せず、長時間安定して積分ができません。

さて、運動方程式を厳密に解くことはできません。そこで、軌道は正確性は妥協しつつ、元々の時間発展が持っていた「シンプレクティック性」は保持するような時間発展は作れないでしょうか？それを行うのがシンプレクティック積分です。

調和振動子の場合、一次のシンプレクティック積分は、

\begin{align}
P &= p - q h \\
Q &= q + P h 
\end{align}

と書けます。一次のオイラー法と異なり、$q$のアップデートに先程アップデートした$P$が使われています。こうするとエネルギーは保存し、時間発展は閉軌道を描きます。本来の円軌道とずれていることを表すため、時間刻みを大きくとってみましょう。

euler.rb

h = 0.5
p = 1.0
q = 0.0
10000.times do
  p = p - q * h
  q = q + p * h
  puts "#{p} #{q}"
end

軌道は楕円になり、円からずれているものの、閉軌道を描いていることがわかります。

さて、なぜシンプレクティック積分がエネルギーを保存するのか考えてみましょう。調和振動子の場合、系が線形なのでヤコビアンが行列になります。まず、オイラー法の場合の時間発展はこのように行列表示できます。

\left(
\begin{matrix}
P \\ Q
\end{matrix}
\right)
= \left(
\begin{matrix}
1 & -h \\
h & 1
\end{matrix}
\right)
\left(
\begin{matrix}
p \\ q
\end{matrix}
\right)

この場合、変数変換のヤコビアン(変換行列の行列式)は$1 + h^2$であり、1より大きいため、時間ステップの度に微小体積要素が大きくなっていきます。これがエネルギーが増加していく原因です。

さて、シンプレクティック積分の場合の変換行列はこうかけます。

\left(
\begin{matrix}
P \\ Q
\end{matrix}
\right)
= \left(
\begin{matrix}
1 & -h \\
h & 1 - h^2
\end{matrix}
\right)
\left(
\begin{matrix}
p \\ q
\end{matrix}
\right)

変換行列の行列式が1になったことがわかるかと思います。

影のハミルトニアン

さて、シンプレクティック積分を実行すると、エネルギーは揺らぐけれど保存します。実際、先程の調和振動子の場合、エネルギーをプロットするとこんな感じになります(時間刻み等を変えています)。

sym1e.rb

h = 0.1
p = 1.0
q = 0.0
t = 0
1000.times do
  puts "#{t} #{p**2+q**2}"
  p = p - q * h
  q = q + p * h
  t = t + h
end

エネルギーが、初期エネルギーである1.0のまわりを揺らいでいるけれど、そこから大きくはずれないことがわかります。これは、シンプレクティック積分が、ハミルトニアンに近い別の量を厳密に保存するからです。これを影のハミルトニアンと呼んだりします。これを求めてみましょう。

元々の保存量$p^2 + q^2$は、行列表示により

H=
p^2 + q^2
= 
\left(
\begin{matrix}
p & q
\end{matrix}
\right)
\left(
\begin{matrix}
1 & 0 \\
0 & 1
\end{matrix}
\right)
\left(
\begin{matrix}
p \\ q
\end{matrix}
\right)

と書けました³。

そこで、影のハミルトニアン$H'$が存在し、その形が

H'=
\left(
\begin{matrix}
p & q
\end{matrix}
\right)
X
\left(
\begin{matrix}
p \\ q
\end{matrix}
\right)

と書けることを仮定しましょう。時間発展行列

U = 
\left(
\begin{matrix}
1 & -h \\
h & 1 - h^2
\end{matrix}
\right)

を用いると、

\left(
\begin{matrix}
P \\ Q
\end{matrix}
\right)
= U
\left(
\begin{matrix}
p \\ q
\end{matrix}
\right)

と書けます。この時、$H'$が時間発展により厳密な保存量となるためには

H'(P,Q) = H'(p,q)

すなわち

\left(
\begin{matrix}
P & Q
\end{matrix}
\right)
X
\left(
\begin{matrix}
P \\ Q
\end{matrix}
\right)
= 
\left(
\begin{matrix}
p & q
\end{matrix}
\right)
X
\left(
\begin{matrix}
p \\ q
\end{matrix}
\right)

が満たされる必要があります。時間刻み$h$が0の時に$H'(p,q)=H(p,q)$となることを要請すると、

X = 
\left(
\begin{matrix}
1 & -h \\
0 & 1
\end{matrix}
\right)

と求めることができます。この時、影のハミルトニアンは

H'(p,q) = p^2 - h pq + q^2

となります。実際にこの量が保存されていることを確認してみましょう。

sym1e_2.rb

h = 0.1
p = 1.0
q = 0.0
t = 0
1000.times do
  puts "#{t} #{p**2+q**2} #{p**2 - h*p*q + q**2}"
  p = p - q * h
  q = q + p * h
  t = t + h
end

影のハミルトニアン$H'=p^2-hpq+q^2$が厳密に保存されていることがわかります。

さて、真のハミルトニアンと影のハミルトニアンの差が$O(h)$であることにも注意しましょう。これは数値積分法が一次のシンプレクティック積分だからです。

詳細は省きますが、時間発展演算子の鈴木トロッター分解により、高次のシンプレクティック積分法を構成することができます。

例えば二次のシンプレクティック積分であれば、

\begin{align}
p &\leftarrow p - q h /2 \\
q &\leftarrow q + p h \\
p &\leftarrow p - q h /2 
\end{align}

として構成できます。この時、一次の場合と同様にして影のハミルトニアンが

H'(p,q) = p^2 + \left(1 - \frac{h^2}{4}\right) q^2

であることがわかります。確認してみましょう。

sym2.rb

h = 0.5
p = 1.0
q = 0.0
t = 0
100.times do
  puts "#{t} #{p**2+q**2} #{p**2 + (1-h**2/4.0)*q**2}"
  p = p - q * h * 0.5
  q = q + p * h
  p = p - q * h * 0.5
  t = t + h
end

影のハミルトニアンが厳密に保存していることがわかります。また真のハミルトニアンとの差が$O(h^2)$であることもわかります。

余談ですが、影のハミルトニアンが楕円でいられる限界の時間刻みが存在することにも注意したいところです。一次、二次の手法ともに$h=2$を超えると楕円から双曲線になってしまうため、時間発展が破綻します。一般にシンプレクティック積分は、有限の時間刻みまでしか正しくありません(有限の収束半径を持つ)。

影のハミルトニアンの求め方についてはこちらの講義ノートも参照してください。

まとめ

要点を列挙します。

• 運動方程式とは位相空間の各点に速度場を導入する写像である
• 時間発展とはその速度場に従って流れるトレーサーの流線である
• ハミルトンの運動方程式が作る流れ場は非圧縮流れとなる
• 確率流が非圧縮であるから、時間発展によりエントロピーは増加しない(リュービルの定理)⁴
• ハミルトンの運動方程式に従うと、時間発展はシンプレクティック写像となる(ヤコビアンが1となる)
• シンプレクティック積分とは、時間発展がシンプレクティック写像となるように近似する手法である
• シンプレクティック積分では、影のハミルトニアンが保存する
• 影のハミルトニアンと真のハミルトニアンとのずれは、シンプレクティック積分の次数に依存する

なお、影のハミルトニアンを厳密に求めることができたのは、時間発展が線形だったからです。一般の非線形の系において影のハミルトニアンが厳密に求まった例を私は知りません。また、シンプレクティック積分の時間刻みは有限の収束半径を持ちますが、やはり調和振動子以外で厳密に収束半径を求められた例を知りません。

というわけで解析力学の幾何学的側面と、シンプレクティック積分についてまとめてみました。もうちょっと多様体とかそのあたりまで説明したかったのですが、時間的にここまでにしておきます。

追記：　多様体というか、そっち方面について続きを書きました→ 解析力学の幾何学的側面II

参考

• Velocity Verlet法とシンプレクティック積分 Velocity Verlet法がシンプレクティック積分になることの確認等
• Lie-Trotter公式の打切り誤差を調べる シンプレクティック積分の構築に使われるLie-Trotter公式の打切り誤差を確認
• 時間反転対称性とシンプレクティック積分
• Lie-Trotter公式における二次の対称分解 分子動力学法で、おそらく最も使われている二次のシンプレクティック積分(二次の対称分解)の精度の確認
• シンプレクティック変換と運動方程式の関係について 時間発展がシンプレクティックであることと、運動方程式が満たすべき条件の関係について

1. 今は簡単のために外力や拘束条件は考えないことにします。 ↩

2. インターネットでうかつにエントロピーについて言及するとエントロピー警察に怒られそうで怖いのですが・・・ ↩

3. 面倒なので1/2のファクターは除いてあります。 ↩

4. 雑な言い方をしています(といちいち逃げをうつ・・・)。 ↩