2 数値計算法

2.1 オイラー法

常微分方程式を数値計算で解く方法の中でもっとも単純ではあるが，最も精度の悪い方法 である．よっぽどのことが無い限り，この方法で微分方程式を計算してはならない．ただ し，常微分方程式を数値計算することのイメージがつかみやすいので，学習する価値はあ る．

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

$$\if 11 \frac{\mathrm{d}y}{\mathrm{d}x} \else \frac{\mathrm{d}^{1} y}{\mathrm{d}x^{1}}\fi =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)+ \if 11 \frac{\mathrm{d}y}{\mathr... ...\frac{\mathrm{d}^{3} y}{\mathrm{d}x^{3}}\fi \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}& \if 11 \frac{\mathrm{d}y}{\mathrm{d}x}... ...\\ &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;
	dx = きざみ幅の計算

	for(計算終了条件){
	  x[i+1]=x[i]+dx;
	  y[i+1]=y[i]+f(x[i],y[i])*dx;
	}

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

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

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)$$ (12)

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

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

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

式(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$を上手に選ぶことにより，式()と同一にできる．

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

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

となる．これを，式(13)と比べると， $\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 精度の検証

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

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

$$\begin{aligned}y_{n+1}&=y_n+\frac{1}{2}(k_1+k_2)\\ &=y_n+\frac{... ...ac{\mathrm{d}^{2} y}{\mathrm{d}x^{2}}\fi h^2+O(h^3) \end{aligned}$$

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

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

$$\begin{aligned}\if 13 \frac{\mathrm{d}y}{\mathrm{d}x} \else \fr... ...l}{\partial x}+f\frac{\partial}{\partial y}\right)f \end{aligned}$$

となる⁴．

一方，ホイン法の3次の精度を表すのは，式(18)の右辺のテ イラー展開の2次の項である．これは，

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

となる．

明らかに，テイラー展開の3次の項である式(19)とホイ ン法の3次の項の式(20)は異なっている．したがって，ホイ ン法は3次の精度がないことが分かる．少なくとも2次の精度があって，3次の精度がない ことが示されわけで，ホイン法は2次の精度であることが証明されたことになる．

2.2.2 中点法

2.2.2.1 漸化式

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

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

である．これを$x_0$の回りでテイラー展開すると，

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

となる．これを，式(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.2.2.2 精度の検証

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

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

が導かれる．ホイン法の場合と同様，これは，式(7)の2次の部分まで等 しいので，少なくとも2次の精度があることが分かる．一方，3次の精度がない ことは，以下の通り明らかである．式(20)と比べ て，微小変位$h$は， $\frac{1}{2}$異なるだけですので，計算結果は，

$$\begin{aligned}\if 12 \frac{\mathrm{d}f(x_n+\frac{h}{2},y_n+\fr... ... \frac{\partial^{1} }{\partial y^{1}}\fi \right)^2f \end{aligned}$$

と直ちに導くことができる．これは，式(19)と異なりま すので，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法や予測子・修正法などがある．進んだ勉強を したいときに，学習するのがよいだろう．例えば文献 [1]等に詳しくかかれてい る．

4次のルンゲ・クッタの公式は，式(27)に示す通りである．そして，これ のイメージは図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
著者: 山本昌志