計算物理学 No.6

1 常微分方程式の数値解法 3

1.1 Runge-Kutta 法

常微分方程式の数値解法で、オイラー法から修正オイラー法にすると精度が劇的に改善した。より高精度な公式を使うと計算の能率はさらに上がる。この中で古典的なものの１つは Runge-Kutta 法である。 4 次の Runge-Kutta 法

k_{1} = hf(x_{n}, y_{n})
k_{2} = hf(x_{n} + \frac{h}{2}, y_{n} + \frac{1}{2}k_{1})
k_{3} = hf(x_{n} + \frac{h}{2}, y_{n} + \frac{1}{2}k_{2})
k_{4} = hf(x_{n} + h, y_{n} + k_{3})
y_{n+1} = y_{n} + \frac{1}{6}(k_{1} + 2k_{2} + 2k_{3} + k_{4}) + O(h^{5})  (1)

で、誤差は累積すると O(h^{4}) である。中間値として 4 回 f(x, y) を計算しなければならないが、精度の向上が著しいのでメリットがある。

Runge-Kutta 法の導出と、誤差の評価については「工学のための応用数値計算法入門」第 8 章 3 節 [1] など参照。

多元連立の常微分方程式に対しては、 k_{1}, k_{2}, k_{3}, k_{4} を配列に置き換えれば良い。手続きは修正オイラー法の連立常微分方程式の項も参考のこと。

1.2 その他

常微分方程式の数値解法には、 Runge-Kutta 法の他、予測子・修正子法、 Richardson 補外を応用した Bulirsch-Stoer 法などある。

参考文献 “Numerical Recipes in C,Fortran”

2 Fortran 文法について

2.1 配列

Fortran では、配列宣言は基本的には 1 から始まる。 0 または負の数から配列宣言したい時には

real*8 x(0:100,-10:20)

と言った風にする。

よくあるミスは、配列宣言の範囲外の変数を読み込み・書き出しするものである。この場合、コンパイルエラーなどにはならないのでしばしば気づき難いが、予想とは異なる出力結果となる。

3 レポート課題

No.5 での課題の他、

\frac{d^{2}y}{dx^{2}} + sin(y) = 0  (2)

（大振幅の振子の問題; 初期値は各自与えること）をオイラー法、修正オイラー法で扱ってみよ。刻み幅に対して誤差がどのように収束するかも調べること。

単に数字を出力結果として出すだけでなく、 gnuplot などでグラフの形にしてみよ。調和振動子や、振り子の課題では、 x, y(x) のグラフだけでなく、 y, dy/dx (位相空間) でも表示して、物理的意味の考察もしてみよ。振り子では、大振幅で回転運動になることに注意。

余裕があれば、 4 次の Runge-Kutta 法にも取り組んでみよ。

この他、各自物理的に興味ある常微分方程式も扱ってみよ。

レポート提出期限は 6 月 23 日まで。

参考文献

[1] 「工学のための応用数値計算法入門」篠崎、松下共著、コロナ社