next up previous
Next: 3 ケプラー問題 Up: 物理シミュレーションの手法と結果の検証 Previous: 1 調和振動

2 シンプレクティック法

学生B:ところで、僕が書いた最初の間違ったプログラムですが、どうしてあれのほうが正しいオイラー法よりもいいんですか?

赤木: そうねえ、まあ偉そうにいうと、あれは「一次のシンプレクティッ ク法」というものにたまたまなっていたからうまくいくのね。ハミルトン力学 系でハミルトニアンが分離可能な時、まあ、普通は要するに運動方程式が

\begin{displaymath}
\frac{d^2 \mathbf{x}}{dt^2} = -\nabla \Phi(\mathbf{x})
\end{displaymath} (9)

と書ける場合ね、そういう時に、一次のシンプレクティック法っていうのは、 例えば
$\displaystyle \mathbf{v}_{i+1}$ $\textstyle =$ $\displaystyle \mathbf{v}_i - \Delta t \nabla \Phi(\mathbf{x}_i)$ (10)
$\displaystyle \mathbf{x}_{i+1}$ $\textstyle =$ $\displaystyle \mathbf{x}_i + \Delta t \mathbf{v}_{i+1}$ (11)

というふうに、最初に速度を更新して、更新した速度を使って位置を更新する というふうに位置と速度で違うやり方で計算する方法なわけ。多変数だったら、 もちろんまず速度の全要素を更新してから位置を更新しないといけないわけね。

これがどうしてうまくいくかってことだけど、理屈はいろいろあるけど、詳しくは吉田さんの解説[2]でもみて。とり あえずさっきのオイラー法みたいに行列で書いてみて。

学生B:あ、はい、、、って、$x_{i+1}$ の式に $v_{i+1}$ があるから行列にならないですけど、、、

赤木: え、でも、$v_{i+1}$$x_i$$v_i$ で書けてるんだから、それを代入すればいいでしょ。

学生B:あ、なるほど、そうすると、

\begin{displaymath}
\mathbf{x}_{i+1} = \pmatrix{
1-\Delta t^2 & \Delta t \cr
-\Delta t & 1 } \mathbf{x}_i
\end{displaymath} (12)

です。左上の要素に $-\Delta t^2$ がつくのが違いますね。

赤木: そうね。形を見ればわかるけど、これは回転+定数倍にはなって ないわけね。 まあ、固有値をまずは求めてみましょ。

学生B:えーと、固有方程式が

\begin{displaymath}
(1-\Delta t^2 -\lambda)(1-\lambda) + t^2 = 0
\end{displaymath} (13)

で、整理すれば
\begin{displaymath}
\lambda^2 - (2-\Delta t^2)\lambda + 1 = 0
\end{displaymath} (14)

あれ、複素解ですよ。 $\Delta t < 2$ なら複素解で、絶対値が 1 ですね。

赤木: そういうこと。だから、固有ベクトルの空間で考えると厳密に円 軌道をとるわけね。もとの空間ではどうなるかしら?

学生B:えーと、なんか楕円ですよね。

赤木: そういうことなわけ。つまり、時間刻みを有限にとっているのに、数値解は何周回っても必ずある楕円の上にのるのよ。

学生B:えーと、それで $x$の誤差が 3次になることが説明できます?

赤木: できるわよ。でも、それあんまり難しくないし、あとで考えてみ て。調和振動だけで話が終わっちゃうと書くほうは楽でいいけど、さすがにそ れだけですますわけにもいかないから話を進めましょう。次は簡単な非線型問 題がいいんだけど、どんなのにしましょ?

学生B:ええと、ケプラー問題とかはどうでしょう?非線型だけど解析解があ るから、とりあえずプログラム書くほうとしては安心できます、、、

赤木: そうね、最終的にはもちろん解析解がないものを考えるんだけど、 調和振動からいきなりなにもわからない世界にいくのもつらいかしらね。でも、 3変数か、平面問題としても2変数よね。プログラムどんなふうに書くのがいい かな、、、まあ、とりあえずおまかせするわ。

学生B:はい、、、



Jun Makino
平成14年10月4日