Subsections

• 3.1 簡単な例題
• 3.2 リープフロッグ公式
• 3.3 数値例
  • 3.3.1 調和振動子
  
  
• 3.4 非線形振動
• 3.5 シンプレクティック公式
• 3.6 陽解法
• 3.7 シンプレクティック公式の意味
• 3.8 シンプレクティック公式の問題点と対応
• 3.9 対称型公式とは
• 3.10 ハミルトン系用の陽的対称型RK公式
• 3.11 対称型線形多段法
• 3.12 エルミート型公式
• 3.13 対称型一段公式における時間刻みの変更

3 ハミルトン系向けの方法

3.1 簡単な例題

古典力学系のもっとも簡単な例といえば、一次元調和振動子、つまり、運動方 程式では

$$\frac{d^2x}{dt^2} = -kx$$ (17)

である。 これの解はもちろん単振動するものであり、それは固有値が純虚数であるから である。

さて、このような、固有値が純虚数である、いいか えれば中立安定な場合、安定性解析はちょっと厄介になる。というのは、通常 の意味で安定というのは、数値解がいつかは原点にいってしまうということ、 いいかえればエネルギーが保存しないということを意味するからである。とい うのは、安定であるとは、線型差分方程式の固有値の絶対値が全て 1 より小 さいということだからである。

つまり、安定領域のある公式では、ほとんどのもので虚軸上の一定範囲では数 値解が「安定」になってしまう。つまり、本当の解はいつまでも振動を続ける のに、数値解は減衰してしまうのである。もちろん、タイムステップを大きく とれば、安定領域から外れるのが普通であり、振動しながら 広がっていくような解になる。安定領域に虚軸を含まないような公式（例えば 前進オイラー公式）の場合では、かならず軌道が広がっていくことになる。つ まり、虚数軸が安定性限界になっていない限り、必ず軌道から一方的に外れて いってしまうのである。

安定領域とはそもそも何かという話をしておく。線形常微分方程式

$$\frac{dy}{dx} = kx$$ (18)

を刻み $h$ で積分した時にi、これは結局 $y_n = C \lambda^n$ という形の一 般解を持つ。ここで、すべての固有値 $\lambda$ の絶対値が1を超えないよう な、複素平面上での $kh$ の領域のことを安定性領域という。

$k$ が複素数の場合も考えるのは、元の微分方程式が連立方程式なら固有値が 実数でない場合もありえるからである。

つまるところ、線形系では安定性限界が虚数軸である（数値解に純虚数の固有 値がある）というのが、数値解が厳密に振動的であるための必要十分条件とい うことになる。この時、数値解はもとの系の性質をある意味で良く表している といえるであろう。

実は多変数の場合や非線形の場合でも本質的な事情は変わらない。ただし、こ ちらはまだいろいろ理論的にはっきりわかっていないことがいろいろある。が、 大雑把にいうと同じような事情が成り立っている。つまり、もとが周期解であ るときに、大抵の積分法では周期解から一方的にずれていく、つまり、保存量 であるはずのエネルギー等が保存しなくなる。

もちろん、エネルギーが厳密には保存しないということが問題であるかどうか というのも実は難しい問題である。というのは、保存しないといってももちろ ん数値解の精度では保存しているわけだから、それでかまわないのではないか とも考えられるからである。

3.2 リープフロッグ公式

ハミルトン系のシミュレーションの例には例えば以下のようなものがある

• プラズマシミュレーション
• （古典）分子動力学
• 銀河動力学

これらはいずれも、古典的な多体問題として定式化される。プラズマの場合は 磁場も入るのでちょっと面倒だが、それを無視すれば運動方程式が

$$m_i\frac{d^2\mbox{\boldmath$x$}_i}{dt^2} = \sum_{j\ne i} \mbox{\boldmath$f$}_{ij}$$ (19)

つまり、ある粒子 $i$ の加速度は他のすべての粒子からの力の合計というこ とになる。もちろん、古典力学系は他にもいろいろあるが、基本的なアプロー チは似たようなものである。

このような問題に対して、もっとも普通に使われる公式は以下の leapfrog といわれるものである。

$$v_{i+1/2} = v_{i-1/2} + \Delta t a(x_i)$$ (20)

$$x_{i+1} = x_{i} + \Delta t v_{i+1/2}$$ (21)

これでは速度と位置がずれた時間でしか定義されない が、出発用公式として

$$v_{1/2} = v + \Delta t a(x_0)/2$$ (22)

を使い、さらに終了用公式として

$$v_{i} = v_{i-1/2} + \Delta t a(x_i)/2$$ (23)

を使うことで最初と最後を合わせることが出来る。この形は、実は

$$x_{i+1} = x_{i} + \Delta t v_{i} + \Delta t^2 a(x_i)/2$$ (24)

$$v_{i+1} = v_{i} + \Delta t [a(x_i)+a(x_{i+1})]/2$$ (25)

と数学的に等価である（証明してみること）また、以下のようにも書ける

$$v_{i+1/2} = v_i + \Delta t a(x_i)/2$$ (26)

$$x_{i+1} = x_{i} + \Delta t v_{i+1/2} + \Delta t^2 a(x_i)/2$$ (27)

$$v_{i+1} = v_{i+1/2} + \Delta ta(x_{i+1})/2$$ (28)

さらにまた、速度を消去して$x_{i-1}$, $x_{i}$, $x_{i+1}$ の関係式の形で 書いてあることもあるかもしれない。なお、すべて違う名前がついていたり するが、結果は丸め誤差を別にすれば完全に同じである。時々、これらが違う ものであるかのように書いてある教科書があるので注意すること。

この公式は局所誤差が $O(\Delta t^3)$、大域誤差が$O(\Delta t^2)$である。

さて、局所誤差という観点からはこれは決して良い公式というわけではないが、 現実にはこの公式は非常に広く使われている。

これは、別にもっとよい方法を知らないからとかではなく、実はこの方法がい くつかの意味で非常に良い方法であるからである。いくつかの意味とは、例え ばエネルギーや角運動量のような保存量が非常によく保存するということであ る。これらの保存量については多くの場合に誤差がある程度以上増えない。

この、ある程度以上誤差が増えないというのは目覚しい性質である。通常の方 法では、エネルギーの誤差は時間に比例して増えていく。従って、長時間計算 をしようとすればそれだけ正確な計算をする必要がある。ところが、エネルギー の誤差は溜っていかないのならば、かならずしも精度を上げる必要はないとも 考えられる。

もちろん、エネルギーが保存していればそれだけで計算が正しいということに はならない。が、理論的には、いくつかの重要な結果が得られている。まず、

1. leapfrog は symplectic method のもっとも簡単なものの一つである。
2. leapfrog は symmetric method のもっとも簡単なものの一つである。

Symplectic method については、以下のようなことが知られている

1. symplectic method は、すくなくともある種のハミルトニアンに対して 使った場合に、それに近い別のハミルトン系に対する厳密解を与えることがあ る。
2. 周期解を持つハミルトン系に対して使った場合に、どんな量でも誤差が 最悪で時間に比例してしか増えない。
3. 時間刻を変えると上のようなことは成り立たなくなる

これに対し、 symmetric method については以下のようなことが知られている

1. 周期解を持つ時間対称な系に対して使った場合に、どんな量でも誤差が 最悪で時間に比例してしか増えない。
2. 時間刻を変えてもうまくいくようにすることも出来る。

以下、まず数値例でこれらの性質を確認し、それからなぜこのようなうまい性 質を持つのかについて直観的な説明を与えよう（厳密な説明／証明には力学系 の積分可能性に関するかなり深い知識が必要となる）。

3.3 数値例

3.3.1 調和振動子

能書きを聞いていても良くわからないので、実例を見てみよう。まず、簡単な 例として調和振動子

$$\frac{d^2x}{dt^2} = -x$$ (29)

を leap frog と 2階のルンゲクッタで解いた例を示す。初期条件は $(x,v) = (1,0)$ で、時間刻みは $1/4$ である。

軌道とエネルギーを図に示す。非常に特徴的なのは、 leap frog ではエネル ギーが周期的にしか変化しないのに対し、ルンゲクッタでは単調に増えている ことである。

Figure 2: 調和振動子の数値積分。左：軌道。右：エネルギー。破線は2次のル ンゲクッタ、実線は leapfrog の結果

ルンゲクッタでは単調に変化するというのは前に説明した通り （2階のルンゲクッタでは虚軸を安定領域に含まないため）である。 さて、これに対し、 leapfrog ではエネルギーが変化していないが、これはど ういうことなのであろう？

実は、この調和振動子の場合には、 leapfrog 公式は以下の量

$$H' = \frac{1}{2}(x^2 + v^2) - \frac{h^2}{4}x^2$$ (30)

を保存するということを確かめることが出来る。つまり、 $(x,v)$ で与えら れる位相平面の上で考えると、 leapfrog 公式の解は上の式で与えられる楕円 の上にのっているのである。このために、エネルギーの誤差がある値よりも大 きくなり得ないことになる。

3.4 非線形振動

では、非線形振動ではどうだろう？簡単な例として

$$\frac{d^2x}{dt^2} = -x^3$$ (31)

を leap frog と 2階のルンゲクッタで解いた例を示す。初期条件は $(x,v) = (1,0)$ で、時間刻みは $1/8$ である。

軌道とエネルギーを図に示す。調和振動の場合と同様に、 leap frog ではエネル ギーが周期的にしか変化しないのに対し、ルンゲクッタでは単調に増えている。

Figure 3: 非線形振動の数値積分。左：軌道。右：エネルギー。破線は2次のル ンゲクッタ、実線は leapfrog の結果

この場合も保存量があるので、頑張れば求まるかもしれない。

3.5 シンプレクティック公式

さて、 leapfrog 公式は、上で見たようにハミルトン系に対してエネルギー誤 差が有界に留まるという大きな特長がある。が、2次精度でしかない。もっと 高次の方法はないのだろうか、またあるとすればどういう原理で作れるのだろ うか。

高次の方法を構成する一つのアプローチが、シンプレクティック公式といわれ るものである。これはなにかというと、積分公式がシンプレクティック写像に なるように作るということである。

シンプレクティック写像とは、ようする に正準変換のことである。正準変換とはなにかというのは解析力学を思い出し てみて欲しいが、直観的には、力学系を不変に保つ、つまり変換まえの座標系で 求めた軌道を変換したものと、変換後の座標系での力学系の軌道が厳密に同じ になるようなものである。

ハミルトン力学系の解そのもの（ある時刻 $t$ での座標から、 $t+h$ での座 標に移す変換）もシンプレクティックである。まあ、だから、シンプレクティッ クになっているような積分公式は、そうでないものに比べてなんとなく力学系 の性質にあっているような気はする。

で、いいたかったことは何かっていうと、上の leapfrog 公式はこのシンプレ クティック性を満たしているということだった。詳しくは、「数理科学」 1995年6月号にのった吉田春夫さんによる解説記事でもみてもらうことにして、ここで は高次の公式にはどんなものがあるかという話をしておく。

3.6 陽解法

陽解法の組み立て方はいろいろある。一つは、 RK 系の公式で、係数をシン プレクティック性を満たすように決めるということである。これはここ15年で 無数に論文がでた。4次、あるいは6次の公式としては、吉田や鈴木による作用 素分解に基づく公式が良く知られている。

これらの方法の原理は、要するに上の leap frog をタイムステップを変えて いくつか組合わせるというものである。うまくタイムステップを組み合わせる と誤差の高次の項を消すことができるわけである。3段4次の公式、7段6次の公 式等が吉田によって導かれている。

なお、実は陽解法はハミルトニアンが $T(p) + V(q)$ の形の場合にしか使え ないが、大抵の問題はこう書けることはいうまでもないであろう。

また、 RK 系の公式を力任せに構成する試みもあり、4次から8次までの公式 が作られている。計算精度という観点からはこちらのほうが leapfrog を組合 わせるものよりもいいものもある。

3.7 シンプレクティック公式の意味

さて、シンプレクティックであるということと、「良い」ということの間には ちゃんとした理論的な関係があるのであろうか？一応あるということになって いる。つまり、あるハミルトニアン $H$ で表される系をシンプレクティック な $p$ 次の公式を使って積分した解は、別のハミルトニアン $H'$で与えられ るシステムの厳密解になっていて、$H$ と $H'$ の間に

$$H = H' + h^{p+1}H_p + O(h^{p+2})$$ (32)

という関係がある（$H'$を求める数列が収束すれば）ということがわかってい る。

なお、例えば上の調和振動子の場合には実際に数列が収束し、$H'$ が求まっ ている。が、これは極めて例外的な場合で、一般の場合には収束するかどうか も明らかではないようである。

収束するかどうか明らかではないのでは、使っていいことがあるという保証は ないではないかと思うかもしれない。理論的にはその通りなのであるが、実際 にはいろいろな問題に適用されて、従来の方法よりも高い精度が得られるとい うことが確認されている。

3.8 シンプレクティック公式の問題点と対応

さて、式(32) をみるとわかるように、シンプレクティック公式 に付随するハミルトニアン$H'$ には時間刻み $h$ が入っている。従って、 $h$ をふらふら変えると $H'$ も変わって、結果的に求まった数値解は一つの 力学系の軌道ではないなんか変なものになってしまう。ということは、可変時 間刻みにするとシンプレクティック公式はうまく働かないのではないかという ことが想像される。

じつはその通りで、例えばシンプレクティックで埋め込み型のルンゲクッタと いうものを作って、実際に時間刻みを変えてみた人がいる。その結果、普通の ルンゲクッタよりも良くならないということが発見された (1992年頃)。問題 によってはこれは致命的な欠陥となるので、さまざまな対応法が精力的に研究 されているのが現状である。が、あまり一般的にうまくいく方法というのは見 つかっていないようである。つまり、 万能な公式というのはシンプレクティック公式に関する限り知られていない。

もう一つのシンプレクティック公式の問題点は、「写像」であるということか ら一段法、具体的にはルンゲクッタ型の公式であるので、局所誤差に対する計 算量という観点からは線形多段階法に比べて必ず悪いということである。

これらを解決しようという研究はいろいろあるが、一つの方向が以下に述べる対称型の公式である。

3.9 対称型公式とは

さて、まず対称型公式とはなにかということだが、まず時間反転に対する対称性というものを定義しておこう。時間反 転に対する対称性とは、常微分方程式

$$\frac{dy}{dx} = f(x,y)$$ (33)

に対する一段法

$$y_{i+1} = F(x_i,y_i,f, \Delta x)$$ (34)

が、

$$y_i = F(x_{i+1},y_{i+1},-f, \Delta x)$$ (35)

を満たすこと、つまり、直観的には、ある微分方程式系があって、それを1ス テップ分数値解を求めたとする。で、そこから逆に戻ると厳密に元の値に戻る ということである。（ここでは丸め誤差はないものとする）

で、一段法の場合には、この意味で対称なものを対称型公式ということにする。

具体例で考えてみよう。例えば前進オイラー法は対称型ではない。これは、いっ た先で導関数を計算すれば、一般にもとのところでの導関数とは違うから当然 であろう。前進オイラー法が対称ではないのだから、その逆写像である後退オ イラー法も対称型ではないことになる。

では、対称型にはどういうものがあるかということだが、明らかに対称型であ るものとしては、台形公式

$$y_{i+1} = y_i + \frac{h}{2} (f_i + f_{i+1})$$ (36)

がある。これが上の対称性を満たしていることはいうまでもない。

さて、なぜこの型の公式を考えるかということであるが、 実験的には以下のようなことが知られている。

• ハミルトン系に対して対称型公式を使った場合、エネルギーや角運動量 などの保存量の誤差がある範囲内にとどまる
• 時間刻みを変えても上の性質を保つことが出来る

前のほうの性質はシンプレクティック公式と同じであるが、後の方の性質はシ ンプレクティック公式よりもある意味でよいものである。

ここではとりあえず対称型公式にはどんなものがあ るかという話をしておこう。

さて、一階の方程式に対する一段法という制限をつけた場合、つまり、ルンゲ クッタ法の場合、対称な公式はかならず陰的公式になる。逆に、陰的公式であ れば対称型のものを作るのは容易である。つまり、ルンゲクッタ法を与える行 列 $A$ とかベクトル $b$, $c$が対称になっていればいい。

一階の系に対する一段法では、陰的 RK のほかにうまい方法があるわけではな い。それでは、

• ハミルトン系専用の解法ではどうだろうか
• 多段法ではどうなっているのだろう
• それ以外の方法はないのだろうか

以下、順番に考えていくことにする。

3.10 ハミルトン系用の陽的対称型RK公式

実は、すでに述べたように、 leapfrog 公式は対称型でありしかも陽公 式である。で、 leapfrog の組合わせで作られる公式も、実は すべて対称型になっている。

3.11 対称型線形多段法

線形多段階法でも、対称型というものを考えることはできる。但し、これら の公式は、今のところ

• 一階の微分方程式用のものは次数が高いと不安定である。

• 二階の微分方程式用のものは、線形解析では安定であっても非線型振動 では数値不安定を起こす。

ということが知られており、実用になっていない。これは 国立天文台の福島さんのグループにより精力的に研究が進められている。

3.12 エルミート型公式

一般に対称型というわけではないが、4次の対称型公式を導く特別な方法とし てエルミート型（あるいはエルミート・オブレヒホフ型）と呼ばれる公式のク ラスについてここで触れておこう。

エルミート型公式は、線形多段階法のアダムス型公式の一つの一般化である。 どのような意味において一般化であるかというと、アダムス型と同じように補 間多項式を積分することで新しい時刻での値を求めるところは同じである。違 うところは、アダムス型公式では補間多項式を作るのに関数の値を使う（ラグ ランジュ・ニュートン補間）が、エルミート型公式では、関数の導関数の値も 使うのである（エルミート補間）。

もっとも単純な発想としては、テイラー法というものがある。つまり、微分方 程式

$$\frac{dy}{dx} = f(x,y)$$ (37)

というものがあったとしたら、

$$\frac{df}{dx} = \frac{\partial f}{\partial x} + \frac{\partial f}{\partial y} \frac{dy}{dx}$$ (38)

という関係式を使って $f$の全微分を求めることが出来、同様に高次の導関数 もどんどん計算していくことが出来る（微分が式で書ければ）。これらを使っ て直接テイラー級数を評価して近似解を求めるのである。

高次の導関数が簡単に求められる場合（例えば線形な系とか）には、この方法 はかなり強力である。線形多段階法やルンゲクッタに比べて誤差項の係数がずっ と小さく、タイムステップが大きくとれるからである。

が、多くの問題において、高階導関数の直接計算は現実的ではない。これは、 計算量が指数関数的に爆発するのが普通だからである。とはいえ、$f$の一階 導関数や二階導関数くらいならば指数関数的といってもたかがしれている。し かしながら、これではテイラー法では2次や3次の公式しか作れない。実用的な 公式としては、（leap frog のような特別に良い性質をもっているとかいうの でなければ）もうちょっと高次のものが必要である。

そこで考えられるのが、高階導関数を使ってさらにルンゲ・クッタ型とか線形 多段階法のような公式を作ることである。例によってもっとも簡単な場合とい うのを考えると、それは $y_i$ と $y_{i+1}$ のところで $f$ の一階導関数 の値を使うもの（高次導関数を使わないものであれば台形公式に相当）の一般 化ということになろう。

台形公式は、2次の陰的アダムス型公式と考えることができる。つまり、 $f$ を線 形補間してそれを積分しているわけである。これに対し、 $f$ と $f' = df/dx$ を使う補間というのはどういうものかということを考えてみると、こ れは 3次多項式を構成できることがわかる。具体的には、導関数の近似多項式 を

$$f(x_0 + h) = f_0 + a h + \frac{b}{2} h^2 + \frac{c}{3!} h^3$$ (39)

としたとき、

$$a = f_0'$$

$$b = { -6(f_0 - f_1) - 2\Delta x (2 f_0' + f_1') \over \Delta x^2}$$

$$c = { 12(f_0 - f_1) + 6\Delta x (f_0' + f_1') \over \Delta x^3}$$ (40)

という形で係数を求めることが出来る。さらに、上の近似多項式を積分してや れば、結局

$$x_1 = x_0 + {h \over 2}(f_0 + f_1) + {h^2 \over 12}(f_0'-f_1')$$ (41)

という形の公式が得られる。これは台形公式と同じく陰的公式である。対称型 であることは明らかであろう。

この公式はプログラムが簡単であるわりには精度が高く、また対称型であるというよ い性質をもつこともあり、良く使われるようになりつつあるようである。収束 させるための反復には普通の直接代入が使える。

3.13 対称型一段公式における時間刻みの変更

対称型一段公式には、対称性を保ったままで時間刻みを変更できるというシン プレクティック公式では（普通の意味では）なかった良い性質がある。以下で は、どのようにしてそれが可能であるかを簡単に述べておく。

今、任意の対称型一段法があって、その時間刻みを与える関数として $h(x,y)$ が与えられているとする。ここで $h(x,y)$ はなんでもいいが、これも一段法 である、つまり前のステップの値とかを必要としないものであるとする。

一般に、 $h(x,y)$ によって刻み幅を決めて求めていった数値解は、時間反転 に対する対称性を満たしていない。つまり、1ステップ積分してから、逆に戻 すと、タイムステップが変わってしまうためにもとのところに戻らない。この ことのために、例えば周期軌道の場合に誤差が周期的にならないということが 起こるわけである。

刻み幅を変えつつ時間反転に対する対称性を保つ一つの方法は（これ以外の方 法もあるかもしれないが、今のところ知られていない）、 $h(x,y)$ 自体に対 称性を持たせることである。つまり、一段法を

$$y_1 = y_0 + h(x_0, y_0)F(x_0,y_0)$$ (42)

という形に書いた時、

$$h(x_0, y_0) = h(x_1, y_1)$$ (43)

がなりたつことを保証するように $h$ を決めるのである。具体的には、上の 対称性が成り立っていないような時間刻みの式 $h$ があった時に、

$$h_s = \frac{h(x_0, y_0) + h(x_1, y_1)}{2}$$ (44)

によって対称化した時間刻みを作ればいいことになる。この場合時間刻み自体 が陰的に決まることになるが、とにかく対称性という要求は満たされるのであ る。