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... ...\\ &x_{i+1}=x_i+\Delta x \\ &x_0=a\\ &y_0=b \end{aligned} \right.$$

これれから、オイラー法での数値計算の漸化式

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

となる。初期値$(x_0,y_0)$が決まれば、 $(x_1,y_1), (x_2,y_2), \cdots$ が同 じ手続きで、芋づる式に計算できるのである。この芋づる式がコンピューターの得意なところで る。通常、初期値$(a,b)$は問題で与えられる。

実際にプログラムを行うときは、forやwhileを用いて繰り返し計算を行う(芋 づる式の部分)。そして、計算結果の$x_i$と$y_i$は、配列x[i] やy[i]に格納 する。

	x[0]=a;
	y[0]=b;

	while(計算終了条件){
	    delta_x や delta_y の計算
	  x[i+1]=x[i]+delta_x;
	  y[i+1]=y[i]+f(x[i],y[i])*delta_x;
	}

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

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

[練習1]

以下の微分方程式をオイラー法で計算して見よ。最初は刻み 幅を2として、$x$の範囲[0,10]で計算せよ。次に、刻み幅を その半分にして見よ。

$$\frac{\mathrm{d}y}{\mathrm{d}x}=2x$$ (12)

初期条件は、$x=0$の時、$y=0$とする。

2.2 2次のルンゲクッタ法

2次のルンゲ・クッタと呼ばれる方法は、いろいろある。ここでは、代表的なホイン法と 中点法を示す。オイラー法は1次の精度であったが、これらは2次の精度になる。

2.2.1 ホイン法

2.2.1.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)$$ (13)

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

$$\begin{aligned}\Delta y &=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) \end{aligned}$$

となっている必要がある。

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

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

とあらわすのである。この $\alpha, \beta$を上手に選ぶことにより、式()と同一にできる。

この式を$x_0$の回りでテイラー展開すると

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

となる。これを、式(14)と比べると、 $\alpha+\beta=1,\quad\beta=1/2$になるので

$$\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の示すところである。オイラー法では、区間の 平均の傾きを出発点だけで決めていたが、ホイン法は両端で決めているのである。これに より、計算精度が向上するのである。

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

2.2.1.2 精度の検証

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

まずは、少なくとも2次の精度があることを確認である。漸化式は、

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

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

次に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}$$

となる⁴。

一方、ホイン法の3次の精度を表すのは、式(19)の右辺のテ イラー展開の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次の項である式(20)とホイ ン法の3次の項の式(21)は異なっている。したがって、ホイ ン法は3次の精度がないことが分かる。少なくとも2次の精度があって、3次の精度がない ことが示されわけで、、ホイン法は2次の精度であることが証明されたことになる。

2.2.2 中点法

2.2.2.1 漸化式

これも、ホイン法と同じ2次の精度である。ホイン法は区間の両端の点の導関数 を使ったが、中点法は出発点と中点で漸化式を作る。先ほど同様、2点を使うので、2次の精度 にすることができる。ホイン法の式(15)に対応するものは、

$$\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)

となる。これを、式(14)と比較すると、

$$\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.2.2.2 精度の検証

式(19)と同じ手順でを用いることにより、中点法が2次の精 度であることが証明できる。漸化式をテーラー展開すると、

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

が導かれる。ホイン法の場合と同様、これは、式(7)の2次の部分まで等 しいので、少なくとも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
著者: 山本昌志