5 微分方程式の差分化

差分化
オイラー法
修正オイラー法
ルンゲ-クッタ法(4次のルンゲ-クッタ法)
など有名な方法がある。
オイラー法は、計算精度が悪い場合があり学習目的以外にはあまり 使われない方法。4次のルンゲ-クッタ法は、計算精度が高い方法である。 しかしながら、計算手続きが多く高速な計算を要求される場合はあまり 使われない。その他に、粒子系のシミュレーションでは、 ベルレ法が良く使われる。

オイラー法
1階常微分方程式

$$\frac{d y}{d x} = f(x,y(x))$$ (28)

の差分化
初期値 $y(x_{0})=y_0$ ,
$x_{n}=x_{0} + n \Delta x$ における近似解$y_{n}$ ($x$ を等間隔に分割.分割幅$\Delta x$ )

$$\frac{d y}{d x} = f(x,y(x))$$ (29)

(29)を区間 $[x_n,x_n+1]$ で形式的に積分すると、

$$\int_{x_n}^{x_{n+1}}\frac{dy}{dx}\cdot dx = \int_{x_n}^{x_{n+1}}f(x, y(x))dx$$ (30)

$$y_{n+1}-y_{n} = \int_{x_n}^{x_{n+1}}f(x, y(x))dx$$ (31)

被積分関数$f$ が$y(x)$ に依存しているので、$y_{n+1}$ が決まったわけではない。

オイラー法では、被積分関数$f$ 中の$y$ を$y=y_{n}$ に近似して積分を実行する。 被積分関数で計算できる値を使う必要があるので$x=x_{n}$ で近似する。

$$y_{n+1}-y_{n} = \int_{x_n}^{x_{n+1}}f(x_n, y_n)dx$$ (32)

$$y_{n+1}-y_{n} = \Delta x f(x_n, y_n)$$ (33)

$$y_{n+1} = y_{n} + \Delta x f(x_n, y_n)$$ (34)

右辺は全て既知の量、左辺を計算する漸化式となっている。このように $y_{n+1}$ が既知の量によりすぐに計算できる方法を陽解法と呼ぶ。 これに対して、陰解法があり、それは$y_{n+1}$ が関数$f$ の中に 含まれていて一般に非線形方程式を解く方法のもの。

式(34)は、テーラー展開の1次まで正しい。テーラー展開を1次まで あらわに書くと

$$y_{n+1} = y_{n} + \Delta x y^{\prime} (x_n) + o((\Delta x)^{2})$$ (35)

であり、今 $f(y_n,x_n)=y^{\prime} (x_n)$ である。差分化による誤差は、 $o((\Delta x)^{2})$ と分かる。つまり$\Delta x$ の1次の精度で正しいと いうことになる。

修正オイラー法
修正オイラー法では、差分化による誤差を$\Delta x$ の2次の精度まで 正しいものを作りたいとする。関数$y$ のテーラー展開の2次では、

$$y_{n+1} = y_{n} + \Delta x y^{\prime} (x_n) + \frac{1}{2}(\Delta x)^{2} y^{\prime\prime} (x_n) + o((\Delta x)^{3})$$ (36)

となるので、 $y^{\prime\prime} (x_n)$ を求めておけばよいことになる。 しかし関数 $f(x, y(x))$ の$x$ に対する微分を求めて置く必要があり、 数値計算上不便なことが多い。(解析的に微分が求められないかも しれないし、数値計算の微分ではさらに誤差を生じる。) 例えば、微分は

$$\frac{d^{2} y}{d x^{2}} = \frac{d}{d x} f(x,y(x))$$ (37)

$$= \frac{\partial f(x, y)}{\partial x} - \frac{\partial f(x, y)}{\partial y} \frac{d y(x)}{d x}$$ (38)

$$= f_{x} + f_{y} f$$ (39)

と計算できるが、あらかじめ$f_{x}$ や$f_{y}$ が計算できていることが求められる。 $\partial$ は偏微分を示す記号であり、全微分では$d$ を微分記号として 用いる。

そこで、微分ではなく、$x$ が異なった値での$y$ の値を用いて、書くことが できるとする。つまり

$$y_{n+1} = y_{n} + \Delta x \{ A f(x_n, y_n) + B f( x_n+ \alpha \Delta x, y_n + \beta \Delta x f(x_n, y_n)) \}$$ (40)

とおく。これを$\Delta x$ について2次まで展開して、テーラー展開の形式と 同じになるように係数 $A, B, \alpha, \beta$ を決定する。 式(40)を展開すると

$$y_{n+1} = y_{n} + \Delta x (A+B) f(x_n, y_n) + (\Delta x)^{2} B ( \alpha f_{x}(x_n, y_b) + \beta f(x_n, y_n) f_{y}(x_n, y_n)$$ (41)

ここで、式(36)と比較して

$$A + B = 1$$ (42)

$$\alpha = \beta$$ (43)

$$B\alpha = \frac{1}{2}$$ (44)

と求められる。この条件を満たす

$$A = B = \frac{1}{2}$$ (45)

$$\alpha = \beta = 1$$ (46)

がよく使われる。このとき、

$$y_{n+1} = y_{n} + \frac{1}{2} \Delta x f(y_n,x_n) + f(x_{n}+ \Delta x , y_{n} + \Delta x f(y_{n},x_{n} ))$$ (47)

これを計算コードに埋め込みやすいように書くと

$$k_{1} = f(x_{n}, y_{n})$$ (48)

$$k_{2} = f(x_{n}+ \Delta x , y_{n} + \Delta x k_{1} )$$ (49)

$$y_{n+1} = y_{n} + \frac{1}{2} \Delta x (k_{1} + k_{2})$$ (50)

これを修正オイラー法と呼んだり、2次のルンゲ-クッタ法と呼んだりする。 計算誤差は、 $o((\Delta x)^{3})$ ということになる。

4次のルンゲ-クッタ(Runge-Kutta)法

ルンゲ-クッタ法と呼んだ場合、通常4次のルンゲ-クッタ法を指すことが多い。 4次のテーラー展開まで正しくなるように、係数を決めるので計算誤差は、 $o((\Delta x)^{5})$ である。結果だけを示す。

$$k_{1} = f(x_{n}, y_{n})$$ (51)

$$k_{2} = f(x_{n}+ \frac{1}{2} \Delta x, y_{n} + \frac{1}{2} \Delta x k_{1} )$$ (52)

$$k_{3} = f(x_{n}+ \frac{1}{2} \Delta x, y_{n} + \frac{1}{2} \Delta x k_{2} )$$ (53)

$$k_{4} = f(x_{n}+ \Delta x, y_{n} + \Delta x k_{3} )$$ (54)

$$y_{n+1} = y_{n} + \frac{1}{6} \Delta x (k_{1} + 2 k_{2} + 2 k_{3} + k_{4})$$ (55)

$k_{1}, k_{2}, k_{3}, k_{4}$ の相互関係から、これらの式を上から順に 計算しなければならないことと、関数$f$ を4回計算しなければならず 必要な計算量が多いことが不便な点である。 一方、$\Delta x$ を小さく採ることで計算誤差は小さくかなり正確な計算値を 得ることができる。例えば、 $\Delta x=0.01$ にとれば、 $(\Delta x)^{5}=10^{-10}$ であるので、9桁程度は正確に計算できることが 予想される。

問：ルンゲ-クッタ法で関数$f$ が$y$ に依存しないとすると、シンプソンの 積分公式を得ることを示せ。(シンプソンの積分公式は、１年次の講義「計算科学」 の配布物を参照)

運動方程式(2階の常微分方程式)

運動方程式は、2階の微分方程式で、次ぎのように1階ずつに分ける。

$$\frac{dv}{dt} = \frac{1}{m} f(t,v,x)$$ (56)

$$\frac{dx}{dt} = v$$ (57)

これらを連立して解くことになる。 一般的に

$$\frac{dv}{dt} = f_{1} (t,v,x)$$ (58)

$$\frac{dx}{dt} = f_{2} (t,v,x)$$ (59)

と書くと、ルンゲ-クッタ法では、

$$k_{1} = \Delta t f_{1}(t_{n}, v_{n}, x_{n})$$ (60)

$$\ell_{1} = \Delta t f_{2}(t_{n}, v_{n}, x_{n} )$$ (61)

$$k_{2} = \Delta t f_{1}(t_{n}+ \frac{1}{2} \Delta t,\ v_{n}+ \frac{1}{2} k_{1},\ x_{n}+ \frac{1}{2} \ell_{1})$$ (62)

$$\ell_{2} = \Delta t f_{2}(t_{n}+ \frac{1}{2} \Delta t,\ v_{n}+ \frac{1}{2} k_{1},\ x_{n}+ \frac{1}{2} \ell_{1})$$ (63)

$$k_{3} = \Delta t f_{1}(t_{n}+ \frac{1}{2} \Delta t,\ v_{n}+ \frac{1}{2} k_{2},\ x_{n}+ \frac{1}{2} \ell_{2})$$ (64)

$$\ell_{3} = \Delta t f_{2}(t_{n}+ \frac{1}{2} \Delta t,\ v_{n}+ \frac{1}{2} k_{2},\ x_{n}+ \frac{1}{2} \ell_{2})$$ (65)

$$k_{4} = \Delta t f_{1}(t_{n}+ \Delta t,\ v_{n}+ k_{3},\ x_{n}+ \ell_{3})$$ (66)

$$\ell_{4} = \Delta t f_{2}(t_{n}+ \Delta t,\ v_{n}+ k_{3},\ x_{n}+ \ell_{3})$$ (67)

$$v_{n+1} = v_{n} + \frac{1}{6} (k_{1} + 2k_{2} + 2k_{3} + k_{4})$$ (68)

$$x_{n+1} = x_{n} + \frac{1}{6} (\ell_{1} + 2\ell_{2} + 2\ell_{3} + \ell_{4})$$ (69)

と書くことができる。

ベルレ法

次のようにベルレ法の漸化式を導出する。 位置座標 $x(t+\Delta t)$ および $x(t-\Delta t)$ をテーラー展開する。

$$x(t+\Delta t) = x(t)+ x^{\prime}(t) \Delta t + \frac{1}{2} x^{\prime\prime}(t)(\D... ...)^{2} +\frac{1}{3!} x^{\prime\prime\prime}(t)(\Delta t)^{3} + o((\Delta t)^{4})$$ (70)

$$x(t-\Delta t) = x(t)- x^{\prime}(t) \Delta t + \frac{1}{2} x^{\prime\prime}(t)(\D... ...)^{2} -\frac{1}{3!} x^{\prime\prime\prime}(t)(\Delta t)^{3} + o((\Delta t)^{4})$$ (71)

この両辺を加えたもの、差し引いたものを作り、式変形すると

$$x^{\prime\prime}(t) = \frac{1}{(\Delta t)^{2}} \left\{ -2x(t) + x(t+\Delta t) + x(t-\Delta t) \right\} +o((\Delta t)^{2})$$ (72)

$$x^{\prime}(t) = \frac{1}{2 \Delta t} \left\{ x(t+\Delta t)-x(t-\Delta t)\right\} + o((\Delta t)^{2})$$ (73)

運動方程式より $x^{\prime\prime}(t)=\frac{f(x)}{m}, x^{\prime}(t)=v(t)$ なので

$$x(t+\Delta t) = 2x(t) - x(t-\Delta t) + \frac{(\Delta t)^{2}}{m} f(x) +o((\Delta t)^{4})$$ (74)

$$v(t) = \frac{1}{2 \Delta t} \left\{ x(t+\Delta t)-x(t-\Delta t)\right\} + o((\Delta t)^{2})$$ (75)

式(74)より位置座標の計算誤差が$\Delta t$ の4次であることが 分かる。座標$x$ と速度$v$ の時刻がずれていることが短所である。この短所を 改良した速度ベルレ法と呼ばれる方法がある。

速度ベルレ法

$$x(t+\Delta t) = x(t) + \Delta t v(t) + \frac{(\Delta t)^{2}}{2m} f(t) +o((\Delta t)^{3})$$ (76)

$$v(t+\Delta t) = v(t) + \frac{\Delta t}{2m} \left\{ f(t+\Delta t) + f(t)\right\} + o((\Delta t)^{3})$$ (77)

問：速度ベルレ法を導出せよ。

問：全エネルギー(運動エネルギー＋ポテンシャルエネルギー)の誤差を、 ベルレ法および速度ベルレ法について評価せよ。