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$を乗じて、求めている。

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)

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

$$\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$を上手に選ぶことにより、式(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)と比較すると、 $\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のようになります。オイラー法の図 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次の精度はないことを示します。

まずは、少なくとも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}$$

となります³。

一方、ホイン法の2次公式の$h^3$の項、即ち式(18)の右辺の テイラー展開の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次の項である式(19)とホイ ン法の3次の項の式(20)は異なっています。したがって、ホ イン法は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})\}$$ (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次の精度である ことの証明は、式(18)と同じ手順で、以下の通りです。

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

これで少なくとも2次の精度があることが分かります。一方、3次の精度がない ことは、以下の通り明らかである。式(20)と比べ て、微小変位$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}$$

これもまた、式(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法や予測子・修正法などがあります。進んだ勉強を したいときに、学習してください。

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
著者: 山本昌志