これはrogy Advent Calendar 2015の17日目の記事です.

こんにちは、この度初めてrogyブログの記事を書かせて頂く事になりました、rogy15のくろすけ(@AstroBoy2030)と申します。
幸福の物理に配属されてますが宗教色は(今のところ)ないです。

はじめに

自分はrogy希少種である理学部の人間なので今回の記事もどちらかと言うと理学系のお話です。
しかし、理学系と言えば工学部やその他の一般の方からしたら象牙の塔の住人の巣窟と思われがちです(個人的な偏見)。
なので工学部の方にも興味を少しでも持っていただけるような話題としてFPU問題を取り上げることにしました。

Fermi-Pasta-Ulamの問題

本題のFermi-Pasta-Ulamの問題ですが、ざっくり言うとFPU問題とは「何個かある格子点を弦のように繋いで両端を固定させたときに隣り合う格子点同士の引っ張り合い方を複雑化(非線形化)させて、ある1つのモードにエネルギーを与えて動かすとモード毎のエネルギーがどのように分配されていくか」といったある種のエネルギー分配問題です。
ここでいう「モード」は弦の基準振動、
要するに山が何個かあるsin波の形のシンプルな振動の事を指します。
この問題は現代的コンピュータが登場した1950年代に研究されていたもので、当時のロスアラモス研究所が保有していた高性能計算機MANIACというコンピュータで数値解析されていたものです。
60年前は研究所で保有する大型コンピュータで行われていた実験が今や個人で持てるノートPCで出来るようになったわけです。
さてFPU問題の数式モデルを作ってみます。
まずは普通の弦を下の図のように$N+1$個の格子点の集合体に捉えたとして、
波動方程式を作るときと同じ要領で格子点$P_{i}$の変位$x_{i}$($i=1,2,...,N-1$)に関する運動方程式を立ててみます。
なお$P_{0}$と$P_{N}$は端点なので固定されていて運動方程式を考える必要はありません。
FPU01

このと
き簡単のため、格子点の質量と復元力の係数を両方$1$としています。
\[\frac{{\rm d}^{2} x_{i}}{{\rm d} t^{2}} = (x_{i+1}-x_{i})-(x_{i}-x_{i-1})\]
次に運動方程式に非線形の項を挿入します。
ここでは非線形項として隣合う格子点同士の相互作用として$1$次の項に続く$2$次の項を採用し、その係数を$\alpha$とします。
\[\frac{{\rm d}^{2} x_{i}}{{\rm d} t^{2}} = (x_{i+1}-x_{i})+(x_{i}-x_{i-1})+\alpha\left[(x_{i+1}-x_{i})^{2}-(x_{i}-x_{i-1})^{2}\right]\]
これが今回取り扱うFPUの運動方程式となります。
因みに元の論文で取り扱われている非線形項はこの$2$次の項に限らず、$3$次の項であったり折れ線グラフで表される複雑な$1$次関数だったりします。
要するに非線形項であれば良いのです。

Runge-Kutta法

まずは微分方程式を数値計算で解く方法として(多分)メジャーなものであるRunge-Kutta法という手法を使います。
Runge-Kutta法は一階の微分方程式を解く方法なので先程の運動方程式を格子点の速度$v_{i}$を使って下の式のように分解します。

\[ \frac{{\rm d} v_{i}}{{\rm d} t} = \left( x_{i+1} - x_{i} \right) - \left( x_{i} - x_{i-1} \right) + \alpha \left[ \left( x_{i+1}-x_{i}\right)^{2} - \left( x_{i}-x_{i-1} \right)^{2} \right] \\
   \frac{{\rm d} x_{i}}{{\rm d} t} = v_{i}\]
この2つの運動方程式を$i=1,2,..,N-1$に関して解けば格子点の位置が求まります。
Runge-Kutta法の手法は何種類かありますが効率よく計算できる4次の手法が使われることが多いので、ここでも4次の方法で計算します。
以下、細かい原理は省いた全体的な計算の流れです(沿え字$i$は省略)。
  1. 時刻$t_{0}$における値$v_{0}$、$x_{0}$を用いて$k_{(1)}=\displaystyle\frac{{\rm d} v}{{\rm d} t}{(v_{0},x_{0},t_{0})}$、$j_{(1)}=\displaystyle\frac{{\rm d} x}{{\rm d} t}{(v_{0},x_{0},t_{0})}$を計算する。
  2. $k_{(2)}=\displaystyle\frac{{\rm d} v}{{\rm d} t}{(v_{0}+k_{(1)}\cdot\frac{dt}{2},x_{0}+j_{(1)}\cdot\frac{dt}{2},t_{0}+\frac{dt}{2})}$、$j_{(2)}=\displaystyle\frac{{\rm d} x}{{\rm d} t}{(v_{0}+k_{(1)}\cdot\frac{dt}{2},x_{0}+j_{(1)}\cdot\frac{dt}{2},t_{0}+\frac{dt}{2})}$を計算
  3. $k_{(3)}=\displaystyle\frac{{\rm d} v}{{\rm d} t}{(v_{0}+k_{(2)}\cdot\frac{dt}{2},x_{0}+j_{(2)}\cdot\frac{dt}{2},t_{0}+\frac{dt}{2})}$、$j_{(3)}=\displaystyle\frac{{\rm d} x}{{\rm d} t}{(v_{0}+k_{(2)}\cdot\frac{dt}{2},x_{0}+j_{(2)}\cdot\frac{dt}{2},t_{0}+\frac{dt}{2})}$を計算
  4. $k_{(4)}=\displaystyle\frac{{\rm d} v}{{\rm d} t}{(v_{0}+k_{(3)}\cdot dt,x_{0}+j_{(3)}\cdot dt,t_{0}+dt)}$、$j_{(4)}=\displaystyle\frac{{\rm d} x}{{\rm d} t}{(v_{0}+k_{(3)}\cdot dt,x_{0}+j_{(3)}\cdot dt,t_{0}+dt)}$を計算
  5. 時間dt後の$x$、$v$の値$x_{dt}$、$v_{dt}$を$v_{dt}=v_{0}+\displaystyle\frac{dt}{6}\cdot(k_{1}+2k_{2}+2k_{3}+k_{4})$、$x_{dt}=x_{0}+\displaystyle\frac{dt}{6}\cdot(j_{1}+2j_{2}+2j_{3}+j_{4})$と計算する
  6. $(v_{dt},x_{dt})\to(v_{0},x_{0})$として手順1.に戻る
以上の操作で微分方程式の解が数値的に解かれます。
イメージとしては、次の時刻に進むとき現在時刻の周りで再計算して誤差を小さくしている感じです。

Runge-Kutta法の欠点

ある程度精度良く微分方程式の解を計算出来るRunge-Kutta法ですが、FPUのシミュレーションをする上でのある欠点を抱えています。
幾ら細かく計算したとしても、連続的な関数を離散的に計算している以上ある程度の誤差は生じます。
ですが、この誤差がシミュレーションを実行していく中で蓄積されていき、ついには計算結果が発散してしまう事が起きてしまうのです。
これはエネルギー保存則がある物理のシミュレーションでは非常にまずい問題であり、実際にFPUシミュレータを作っていて簡単にエネルギーが発散してしまいシミュレーションにならない事も多々ありました(というか長めにシミュレーションしてると大体発散しました)。

Symplectic積分

Runge-Kutta法だとエネルギーが発散し易いので別の手法を使います。
ここで先程立てたFPUのモデルに関してHamiltonianという新しい関数を導入します。
このHamiltonian(以下$H$)は運動量$p_{i}$と座標$q_{i}$の関数($i=0,1,...,N$)として表された運動エネルギー$T(p,q)$とポテンシャルエネルギー$V(p,q)$で
\[H=T+V\]
と定義されます。

※注意※

このHamitonianの定義は一般的なものではありません。
ここではあたかも$H$が全エネルギー$E=T+V$そのものに見えてますが、必ずしもHamiltonianが全エネルギーEを表したものとは限らないのです。
ですが、ここではHamiltonianが全エネルギーに一致する場合で話を進めます。

さて、FPUのモデルをHamiltonianで書き直してみます。
ある格子$P_{i}$に対応する運動量と座標をそれぞれ$p_{i}$、$q_{i}$とします。
ここではHamiltonianが全エネルギーの式と一致するとしているので、格子毎の運動エネルギーと弾性エネルギーを足し合わされば$H$を計算できます。
ここでは簡単のため格子点の質量$m$とバネ係数を両方$1$としてあります。
\[H=\sum_{i=0}^{N}\frac{p_{i}^{2}}{2}+\sum_{i=1}^{N}\frac{1}{2}(q_{i}-q_{i-1})^{2}+\sum_{i=1}^{N}\frac{\alpha}{3}(q_{i}-q_{i-1})^{3}\]
さて、Hamlitonianを導入したところで次に下の式で表される正準方程式なる式を使います。
\[\frac{\partial H}{\partial p_{i}}=\frac{{\rm d} q_{i}}{{\rm d} t} \\
  \frac{\partial H}{\partial q_{i}}=-\frac{{\rm d} p_{i}}{{\rm d} t}\]
この2式はHamiltonianを用いて表された運動方程式で、この方程式の解が格子点の運動を決定します。
そこで正準方程式を使って格子点の運動を数値的に解析する事を考えます。
単純に正準方程式の両辺を時間$t$で積分すると$p$、$q$の数値解を得ることが出来ますが、あまり精度が良くないのでより効率的な手法を用います。
先程の式をベクトルの形で書き直します($D$:Hamiltonianに依存する演算子)。
\[ D \left ( \begin{array}{ccc} p_{i} \\ q_{i} \end{array} \right ) = \displaystyle\frac{\rm d}{{\rm d} t}\left ( \begin{array}{ccc} p_{i} \\ q_{i} \end{array} \right ) \]
さて、上の式の通りベクトル$\left (\begin{array}{ccc} p_{i}\\q_{i}\end{array}\right)$に$D$を作用させると$\left (\begin{array}{ccc} p_{i}\\q_{i}\end{array}\right)$を時間微分したものが得られるので、$D$を$n$回だけ作用させると$n$回時間微分したものが得られるわけです。
ここで$\left (\begin{array}{ccc} p_{(\Delta t)}\\q_{(\Delta t)}\end{array}\right)$を$t=0$周りでTaylor展開すると
\[
\left (\begin{array}{ccc} p_{(\Delta t)}\\q_{(\Delta t)}\end{array}\right)=
\left (\begin{array}{ccc} p_{(0)}\\q_{(0)}\end{array}\right)+
\displaystyle\frac{1}{1!}\frac{\rm d}{{\rm d} t}\Delta t\left(\begin{array}{ccc} p_{(0)}\\q_{(0)}\end{array}\right)+\cdots+
\displaystyle\frac{1}{m!}\frac{{\rm d}^{m}}{{\rm d}^{m} t}(\Delta t)^{m}\left(\begin{array}{ccc} p_{(0)}\\q_{(0)}\end{array}\right)+\cdots
\]
ここで$D^{m} \left ( \begin{array}{ccc} p \\ q \end{array} \right ) = \displaystyle\frac{{\rm d}^{m}}{{\rm d}^{m} t}\left ( \begin{array}{ccc} p \\ q \end{array} \right )$であるので$D^{0}={\rm i.d.}$と約束すると
\[
\left (\begin{array}{ccc} p_{(\Delta t)}\\q_{(\Delta t)}\end{array}\right)=
\frac{1}{0!}D^{0}(\Delta t)^{0}\left (\begin{array}{ccc} p_{(0)}\\q_{(0)}\end{array}\right)+
\displaystyle\frac{1}{1!}\left(\Delta t D\right)^{1}\left(\begin{array}{ccc} p_{(0)}\\q_{(0)}\end{array}\right)+\cdots+
\displaystyle\frac{1}{m!}\left(\Delta t D\right)^{m}\left(\begin{array}{ccc} p_{(0)}\\q_{(0)}\end{array}\right)+\cdots
\]
と表すことができるので、結果的に
\[
\left (\begin{array}{ccc} p_{(\Delta t)}\\q_{(\Delta t)}\end{array}\right)={\rm Exp}\left(\Delta t D\right)\left(\begin{array}{ccc} p_{(0)}\\q_{(0)}\end{array}\right)\\
\left({\rm Exp}\left(\Delta t D\right)=\frac{1}{0!}(\Delta t D)^{0}+\frac{1}{1!}(\Delta t D)^{1}+\cdots+\frac{1}{m!}(\Delta t D)^{m}\cdots\right)
\]
と書けます。
上記の${\rm Exp}(\Delta t D)$の部分を$p$による要素と$q$による要素に分解してそれぞれの厳密な積分を行う事で$\Delta t$後の$p_{(\Delta t)}$と$q_{(\Delta t)}$が得られます。
${\rm Exp}(\Delta t D)$を$p$による要素と$q$による要素に分解するとき、これらの要素による作用子が可換であると仮定すると先程の正準方程式を単純に両辺時間積分した形となりますが、実は可換であるとは限らず一般には次式で作用子${\rm Exp}(\Delta t D)$が近似的に展開されます。
\[{\rm Exp}(\Delta t D)\simeq\prod_{i=1}^{N} {\rm Exp}(u_{i} \Delta t D_{p,i}) {\rm Exp}(k_{i} \Delta t D_{q,i})\]
$u_{i}$と$k_{i}$はそれぞれ総和が$1$となる$\Delta t D_{p,i}$と$\Delta t D_{q,i}$の係数です。
このとき係数は一意に定まらないので計算しやすいものを選ぶ余地があります。
ここで$N\to\infty$の極限をとると誤差が無い厳密解となります。
今回作成したシミュレータでは$N=2$として係数$\left(u_{1},u_{2},k_{1},k_{2}\right)=\left(\displaystyle\frac{1}{2},\displaystyle\frac{1}{2},1,0\right)$で計算しています(これは2次のSymplectic積分と呼ばれる)。

(外見だけ)完成した(気持ち悪い)シミュレータ

以下は山が1つのsin波を初期状態としてシミュレーションした結果です

この通りです。
見た目キモすぎてポです。
FPU-SimplecticGraph
これがFPUシミュレータで計算された系のエネルギーを散布図化したものです。
青が運動エネルギー、橙がポテンシャルエネルギー、灰色が全エネルギーを表しています。
みた通り気持ち悪いほど全エネルギーの値がぶれてません。
FPU-RungeKuttaGraph

因みに同じ条件を使ってRungeKutta法で計算した系のエネルギーも散布図化してみました。
気持ち悪いぶれ方をしてます。

関連資料とか 

  1. E.Fermi, J.Pasta, S.Ulam and M. Tsingou (1955), "Studies of nonlinear probems I", Los Alamos preprint LA-1940 → 本文(ミラノ大学)
  2. Joseph Ford (1992), "THE FERMI-PASTA-ULAM PROBLEM:TURNS DISCOVERY", School of Physics, Georgia Institute of Technology, Altanta, Ca 30332, USA → 本文(CBPF)
  3. Haruo Yoshida (1990), "Construction of higher order symplectic integrators", Physics Letters A, Volume 150, number 5,6,7 → 本文(CERN)  
  4. Symplectic数値積分法
1番はFPU問題の最初の論文で1955年に出されたものです。
そして2番はFPU問題のモデルをHamiltonianで書き換え、更なる解析を行ったものです(最期まで読んでませんが大体そんな感じかと)。
今回用いたHamiltonianもこの論文から引用したものを使わせて頂きました。
3番は高次のSymplectic積分に関する論文です。
前述の通りSymplectic積分では非可換な操作を行うため、作用子の係数は簡単に決められるものではありませんが、作用子の係数にフラクタル的構造を見出して漸化式を立て係数を求める事が出来るそうです。
今回は2次のSymplectic積分を適用しましたが、この論文では8次のSymplectic積分の係数まで求められています。  
4番はSymplectic積分を実装する際に参考にしたサイトです(このサイトをChromeで開くと文字化けしました。Firefoxで開くと何ともなかったので文字コードとかの問題かもしれません)

まとめ

Symplectic積分は†光†
あとモードのエネルギーがなんちゃらかんちゃらとかほざいておきながらモード別のエネルギーをプロットしてないのでFFT実装してスペクトル展開的なのをやりたいです(粉蜜柑)。