計算物理学:No.5

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

常微分方程式の数値解法で、単純なオイラー法よりも精度が上がる方法を紹介する。一つは修正オイラー法、もう一つは 2 次の Runge-Kutta 法である。 1 変数の常微分方程式については簡単に拡張できるが、多変数の常微分方程式に関してはやや注意がいる。

1.1 修正 Euler 法

y(x + h) の厳密な値は、積分を用いて

y(x + h) = y(x) + ∫_{x}^{x+h} f(x, y(x))dx  (1)

と表すことができる。オイラー法は、この積分を

∫_{x}^{x+h} f(x, y(x))dx ≈ hf(x, y(x))  (2)

で近似すること (図 1-a) に相当するが、離散近似による誤差をもっと小さくするためには、台形公式を用いて、

∫_{x}^{x+h} f(x, y(x))dx ≈ \frac{h}{2}(f(x, y(x)) + f(x + h, y(x + h)))  (3)

で近似する (図 1-b)。

しかし、微分方程式の計算の場合には数値積分と違い、被積分関数の中に未知関数があり、 (3) の右辺には左辺で予測すべき量が含まれている。

そこで、まずオイラー法で y(x + h) の近似値

ỹ = y(x) + hf(x, y(x))  (4)

を計算し、それを (3) の右辺に代入する

y(x + h) = y(x) + \frac{h}{2} [f(x, y(x)) + f(x + h, ỹ)]  (5)

これが修正オイラー法である。

図 1: (a) オイラー法 (b) 修正オイラー法 (c) 2 次のルンゲクッタ法に相当。

1.2 連立常微分方程式

一般には (多変数、高階の）常微分方程式は m 元連立 1 階常微分方程式

\frac{dy^{i}(x)}{dx} = f^{i}(x, y^{1}, y^{2}, …, y^{m}),  i = 1, …, m  (6)

になるが、その場合も y, f をベクトル（配列）の形にすれば良い。この場合修正オイラー法の手続きは次のようになる。パラメーター x_{n} = a + nh に対して、求めるべき値の y^{i}(x_{n}) の近似値を y_{n}^{i} で表す。

\underline{ステップ 0} x, y^{i} の組が与えられた時、

f^{i}(x, y^{1}, y^{2}, …, y^{m})  (7)

を i = 1, ⋯ m まで 1 次元配列として出すサブルーチンを作成する (これは問題ごとに異なる)。

\underline{ステップ 1}

ỹ^{i} = y_{n}^{i} + hf^{i}(x_{n}, y_{n}^{i})  (8)

を i = 1, ⋯ m まで 1 次元配列として出す。 \underline{ステップ 2}

y_{n+1}^{i} = y_{n}^{i} + \frac{h}{2}[f^{i}(x_{n}, y_{n}^{i}) + f^{i}(x_{n+1}, ỹ^{i})]  (9)

で、次のステップの y^{i} を出す。 \underline{ステップ 3} x_{n+1} と、求めた y_{n+1}^{i} を新たにステップ 1,2 での x_{n}, y_{n}^{i} に置き換える。すると次に y_{n+2}^{i}

が求まる。以下順次これを繰り返す。

1.3 修正オイラー法の誤差

オイラー法の誤差は刻み幅 h に対し O(h) 程度なので精度を上げようとすると計算量が非常に増える。これに対し、修正オイラー法の誤差は O(h^{2}) なのでずっと効率が良い。

y(x + h) の x の回りでのテイラー展開を考えると、

y(x + h) = y(x) + hf(x) + \frac{1}{2}h^{2}y′′(x) + O （h^{3}）  (10)

逆に y(x) の x + h の回りでのテイラー展開を考えると、

y(x) = y(x + h) − hf(x + h) + \frac{1}{2}h^{2}y′′(x + h) + O （h^{3}）  (11)

この 2 つの式から、

y(x + h) = y(x) + \frac{h}{2}[y′(x) + y′(x + h)] + O(h^{3})
= y(x) + \frac{h}{2}[f(x) + f(x + h)] + O(h^{3})  (12)

なので、修正オイラー法の常微分方程式全体としての誤差は、 O(h^{2}) に押えられる。

修正オイラー法はやや複雑であるが、 100 倍精度を上げるのに、 10 倍の計算量で済むので効果的である。さらにこの結果と収束加速法の組合せも考えられる。

1.4 2 次の Runge-Kutta 法

離散近似による誤差を小さくするには、台形公式の他、中点公式を用いて

∫_{x}^{x+h} f(x, y(x))dx ≈ hf (x + \frac{h}{2}, y(x + \frac{h}{2}))  (13)

で近似することも考えられる (図 1-c)。しかし、 y(x + h/2) の評価が必要なのでオイラー法で y(x + h/2) の近似値

ỹ = y(x) + \frac{h}{2}f(x, y(x))  (14)

を計算し、それを (13) の右辺に代入する

y(x + h) ≈ hf (x + \frac{h}{2}, ỹ)  (15)

これを 2 次のルンゲクッタ法という。誤差は O(h^{2}) である。

2 課題

前回の課題を修正オイラー法で扱ってみよ。刻み幅を半分にすると誤差が約 1/4 になることを確認せよ。余裕があれば 2 次のルンゲクッタ法でも扱ってみよ。