2 数値計算法

2.1 オイラー法

常微分方程式を数値計算で解く方法として、もっとも単純ですが、最も精度の 悪い方法です。よっぽどのことが無い限り、この方法で微分方程式を計算して はいけません。ただし、常微分方程式を数値計算することのイメージはつかみ やすいでで述べておきます。

もう一度、初期条件を含めて数値計算により解くべき方程式を示します。

$$\frac{dy}{dx}=f(x,y) \hspace{10mm} 初期条件 \hspace{3mm} y(a)=b$$ (6)

この微分方程式の解を$y=y(x)$とすると、$x_i$のまわりのテイラー展開は、

$$y_{i+1}=y(x_i+\Delta x) =y(x_i)+\frac{dy}{dx}\Bigm\vert _{x=x_i}\... ...x_i}\Delta x^2+ \frac{1}{6}\frac{d^3y}{dx^3}\Bigm\vert _{x=x_i}\Delta x^3+\dots$$ (7)

です。この式の右辺第2項は、式(6) から計算できます。したがって、テイラー展開は、次のように書き表すことが 出来ます。

$$y_{i+1}=y_i+f(x_i,y_i)\Delta x+O(\Delta x^2)$$ (8)

オイラー法での数値計算では、計算の刻み幅$\Delta x$は十分に小さいとして、

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

を計算します。式(5)と全く同じです。このとき計算の精 度は1次と言います²。

オイラー法をまとめると、以下に示すように微分方程式は差分方程式に変換さ れます。

$$\left\{ \begin{aligned}&\frac{dy}{dx}=f(x,y)\\ &y(a)=b \end{ali... ...y_i}{\Delta x}=f(x_i,y_i)\\ &x_0=a\\ &y_0=b \end{aligned} \right.$$

オイラー法での数値計算は、次の漸化式に従い計算します。解である $(x_0,y_0), (x_1,y_1), (x_2,y_2), \cdots$が同じ手続きで計算できます。 実際にプログラムを行うときは、forやwhileを用いて繰り 返し計算を行い、結果の$x_i$と$y_i$は、配列x[i] や y[i]に格納します。

$$\left\{ \begin{aligned}&x_0=a\\ &y_0=b\\ &x_{i+1}=x_i+\Delta x\\ &y_{i+1}=y_i+f(x_i,y_i)\Delta x\\ \end{aligned} \right.$$

この方法の計算のイメージは、図4の通りです。明らか に、出発点の導関数のみ利用しているために精度が悪いことが分かります。式 も対称でないため、逆から計算すると元に戻りません。

図 4: オイラー法。ある区間での$y$の変化$\Delta y$は、計算の始めの 点の傾きに区間の幅$\Delta x$を乗じて、求めている。

2.2 2次のルンゲクッタ法

2次のルンゲ・クッタと呼ばれる方法は、いろいろあります。ここでは、教科 書に載っているホイン法と中点法を示します。オイラーは1次の精度でしたが、 これらは2次の精度があります。

2.2.1 ホイン法

先に示したように、オイラー法の精度は1次です。それに対して、2次のルンゲ・ クッタ法の精度は2次です。教科書に書いてある方法で公式を導き出し、その 精度を示します。今まで刻み幅を$\Delta x$と記述していましたが、教科書に 合わせて$h$と表現します。

2次の精度ということは、テイラー展開より

$$y(x_0+h)=y(x_0)+y^\prime(x_0)h+\frac{1}{2}y^{\prime \prime}(x_0)h^2+O(h^3)$$ (12)

となっていることを意味します。即ち、計算アルゴリズムが、

$$\Delta y=y^\prime(x_0)h+\frac{1}{2}y^{\prime \prime}(x_0)h^2+O(h^3)$$ (13)

となっていることを言います。最終目標が分かったので、教科書と同じ方法で 計算公式を導きましょう。

式(13)から分かるように、$y$の増分$\Delta y$を計算する ためには、1階微分と2階微分の2項を満たす式が必要です。そうすると少なく とも、2点の値が必要となります。2点として、計算区間の両端の導関数の値を 使います。この導関数は問題として与えられているので、計算は簡単です。そ うして、区間の増分を $\alpha, \beta$をパラメーターとした和で表すことに します。即ち、以下の通りです。

$$\Delta y=h\{\alpha y^\prime(x_0)+\beta y^\prime(x_0+h)\}$$ (14)

この $\alpha, \beta$を上手に選ぶことにより、式(13)と同 一にできます。この式を$x_0$の回りでテイラー展開すると

$$\Delta y=(\alpha+\beta)y^\prime(x_0)h+\beta y^{\prime\prime}(x_0)h^2+O(h^3)$$ (15)

となります。これを、式(13)と比較すると、

$$\left\{ \begin{aligned}\alpha &=\frac{1}{2}\\ \beta &=\frac{1}{2} \end{aligned} \right.$$

とすれば良いことが分かります。これで、必要な式は求まりました。まとめる と、式(6)を数値計算で近似解を求める には次式を使うことになります。

$$\left\{ \begin{aligned}k_1&=hf(x_n,y_n)\\ k_2&=hf(x_n+h,y_n+k_1)\\ y_{n+1}&=y_n+\frac{1}{2}(k_1+k_2) \end{aligned} \right.$$

この式は、図5のようになります。オイラー法の図 4との比較でも、精度が良いことが分かります。

図 5: ホイン法(教科書の方法)。ある区間での$y$の変化$\Delta y$は、 計算の始めと終わりの点付近の平均傾きに区間の幅$\Delta x$を乗じて、求 めている。

よく見ると、この式(17)は、本当に2次の精度があるのでしょう か?。$\alpha$や$\beta$のパラメーターを計算したときの$x+h$の導関数は $y^\prime(x+h)$を使いました。一方、式(17)では、 $f(x_n+h,y_n+k_1)$を使っています。ほんのちょっと違いますので、式 (17)が2次の精度を持っているか、検証してみましょう。この式を 変形することで、精度を確認しましょう。紙面の都合上、精度の確認は2段階 で行います。まず初めは、少なくとも2次の精度があることを確認します。そ の後、3次の精度はないことを示します。

$$\begin{aligned}y_{n+1}&=y_n+\frac{1}{2}(k_1+k_2)\\ &=y_n+\frac{... ...1}{2}\frac{d^2y}{dx^2}\Bigm\vert _{x=x_n}h^2+O(h^3) \end{aligned}$$

この結果は、まさに式(7)と同じ形をしており、少なくとも 2次の精度があることが確認できます。式(18)を変形 するときに、次式を用いたので注意が必要です。

$$\begin{aligned}\frac{d^2y}{dx^2}&=\frac{d}{dx}\left(\frac{dy}{d... ...rtial f}{\partial x}+\frac{\partial f}{\partial y}f \end{aligned}$$

次に3次の精度がないことを示します。テイラー展開の3次の項は、係数は無視 すると、

$$\begin{aligned}\frac{d^3y}{dx^3}&=\frac{d}{dx}\left(\frac{d^2y}... ...l}{\partial x}+f\frac{\partial}{\partial y}\right)f \end{aligned}$$

となります。一方、ホイン法の2次公式の$h^3$の項、即ち式()の右辺のテイラー展開の2次の項は、以下の通りです。

$$\begin{aligned}\frac{d^2 f(x_n+h,y_n+hf(x_n,y_n))}{dh^2}&=\frac... ...{\partial x}+f\frac{\partial}{\partial y}\right)^2f \end{aligned}$$

明らかに、テイラー展開の3次の項である式()とホイン法の3次の項の式(21)は異なってい ます。したがって、ホイン法は3次の精度がないことが分かります。少なくと も2次の精度があって、3次の精度がないことが示され、ホイン法は2次の精度 であることが証明されました。

2.2.2 中点法

これも、ホイン法と同じ2次の精度です。ホイン法は区間の両端の点の導関数 を使いましたが、中点方は始点と中点を使います。2点ありますので2次の精度 になります。ホイン法の式(14)に対応するものは、

$$\Delta y=h\{\alpha y^\prime(x_0)+\beta y^\prime(x_0+\frac{h}{2})\}$$ (22)

となります。これを$x_0$の回りでテイラー展開すると、

$$\Delta y=(\alpha+\beta)y^\prime(x_0)h+\frac{\beta}{2} y^{\prime\prime}(x_0)h^2+O(h^3)$$ (23)

となります。これを、式(13)と比較すると、

$$\left\{ \begin{aligned}\alpha &=0\\ \beta &=1 \end{aligned} \right.$$

となります。したがって、中点法の公式は、

$$\left\{ \begin{aligned}k_1&=hf(x_n,y_n)\\ k_2&=hf(x_n+\frac{h}{2},y_n+\frac{k_1}{2})\\ y_{n+1}&=y_n+k_2 \end{aligned} \right.$$

となります。この公式は、図6のようになります。これが2次の 精度であることの証明は、式(18)と同じ手順で、以下の 通りです。

$$\begin{aligned}y_{n+1}&=y_n+k_2\\ &=y_n+hf(x_n+\frac{h}{2},y_n+... ...1}{2}\frac{d^2y}{dx^2}\Bigm\vert _{x=x_n}h^2+O(h^3) \end{aligned}$$

これで少なくとも2次の精度があることが分かります。一方、3次の精度がない ことは、以下の通り明らかである。式(21)と比べ て、微小変位$h$は、 $\frac{1}{2}$異なるだけですので、計算結果は、 と直ちに分かります。

$$\begin{aligned}\frac{d^2 f(x_n+\frac{h}{2},y_n+\frac{hf(x_n,y_n... ...{\partial x}+f\frac{\partial}{\partial y}\right)^2f \end{aligned}$$

これもまた、式(20)と異なりますので、3次 の精度がないことが分かります。

図 6: 中点法。ある区間での$y$の変化$\Delta y$は、中点付近の傾きに 区間の幅$\Delta x$を乗じて、求めている。

2.3 4次のルンゲ・クッタ法

今まで示したオイラー法や2次のルンゲ・クッタ法のように、パラメーターを 増やして誤差項の次数を上げていく方法で、最良の方法と言われるのが4次の ルンゲ・クッタ法です。パラメーターを増やして、5, 6, 7, $\cdots$と誤差 項を小さくすることは可能ですが、同じ計算量であれば4次のルンゲ・クッタ の刻み幅を小さくするほうが精度が良いです。私は、5次以上のルンゲ・クッ タの公式は見たことがありません。

ということで、皆さんが常微分方程式を計算する必要が生じたときは、何はと もあれ4次のルンゲ・クッタで計算してください。「この問題を解く場合、4次 のルンゲクッタだな」と一言いって、プログラムを書き始めると、出来るなと 思われること間違いなしです。間違っても「2次のルンゲ・クッタ$\cdots$」 と言ってはいけません。「4次の方が$\cdots$」と言う輩が出てきます。普通 の科学に携わる人にとって、4次のルンゲ・クッタは常微分方程式の最初で最 後の解法なのです。

ただし、4次のルンゲ・クッタ法よりも精度の良い方法がないわけでは有りま せん。より高精度な方法として、Bulirsch-Store法や予測子・修正法などがあ ります。進んだ勉強をしたいときに、学習してください。

4次のルンゲ・クッタの公式は、式(28)に示す通りです。そし て、これは図7のように表せます。

2次の場合と同じ手順で、公式を導き、そして4次の精度であることが証明でき ると思います。しかし、計算は明らかに大変なので、腕力のある人はトライし てみてください。

$$\left\{ \begin{aligned}k_1&=hf(x_n,y_n)\\ k_2&=hf(x_n+\frac{h}{... ...y_{n+1}&=y_n+\frac{1}{6}(k_1+2k_2+2k_3+k_4) \end{aligned} \right.$$

図 7: 4次のルンゲ・クッタ法。ある区間での$y$の変化$\Delta y$は、区 間内の4点の傾きのある種の加重平均に幅$\Delta x$を乗じて、求めている。

ホームページ: Yamamoto's laboratory
著者: 山本昌志