2 常微分方程式の数値計算法

数値計算により、近似解を求める微分方程式は、

である。これは問題として与えられ、この式に従う$x$と$y$の関係を計算する。

2.1 オイラー法

常微分方程式の解を$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$$ (8)

である。この式の右辺第2項は、式(7) から計算できる。したがって、テイラー展開は、次のように書き表わせる。

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

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

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

を計算する。この場合、計算の精度は1次と言う²。

オイラー法では、次の漸化式に従い数値計算を進める。解である $(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.$$

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

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

2.2 2次のルンゲクッタ法(ホイン法)

2次のルンゲ・クッタと呼ばれる方法は、いろいろある。ここでは、ホイン法をしめす。

先に示したように、オイラー法の精度は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.$$

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

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

この式は、図2のようになる。オイラー法の図 1との比較でも、精度が良いことが分かる。

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

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

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

ということで、皆さんが常微分方程式を計算する必要が生じたときは、何はと もあれ4次のルンゲ・クッタで計算すべきである。普通の科学に携わる人にとっ て、4次のルンゲ・クッタは常微分方程式の最初で最後の解法である。

4次のルンゲ・クッタの公式は、式(18)に示す通りである。そし て、これは図3のように表すことができる。

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

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

2.4 プログラム(4次のルンゲ・クッタ法)

実際の微分方程式

を4次のルンゲ・クッタ法で計算するプログラムを示す。計算結果は、配列「x[]」と 「y[]」に格納される。実際にプログラムでは、この結果を利用してグラフにしたりする のであるが、ここでは計算のみとする。

	#include <stdio.h>
	#include <math.h>
	#define IMAX 100001
	double func(double x, double y);

	/*================================================================*/
	/*       main function                                            */
	/*================================================================*/
	int main(void){
	  double x[IMAX], y[IMAX];
	  double final_x, h;
	  double k1, k2, k3, k4;
	  int ncal, i;


	  /*--- set initial condition and cal range ---*/

	  x[0]=0.0;
	  y[0]=0.0;

	  final_x=10.0;
	  ncal=10000; 

	  /* --- size of calculation step --- */

	  h=(final_x-x[0])/ncal;

 	
	  /* --- 4th Runge Kutta Calculation --- */

	  for(i=0; i < ncal; i++){
	    k1=h*func(x[i],y[i]);
	    k2=h*func(x[i]+h/2.0, y[i]+k1/2.0);
	    k3=h*func(x[i]+h/2.0, y[i]+k2/2.0);
	    k4=h*func(x[i]+h, y[i]+k3);

	    x[i+1]=x[i]+h;
	    y[i+1]=y[i]+1.0/6.0*(k1+2.0*k2+2.0*k3+k4);
	  }


	  return 0;
	}


	/*================================================================*/
	/*       define function                                          */
	/*================================================================*/
	double func(double x, double y){
	  double dydx;

	  dydx=sin(x)*cos(x)-y*cos(x);

	  return(dydx);
	}

2.5 高階の常微分方程式

2.5.1 1階の連立微分方程式に変換

ここまで示した方法は、1階の常微分方程式しか取り扱えないので不便である。そこで、 高階の常微分方程式を1階の連立微分方程式に直す方法を示す。要するに、高階の常 微分方程式を連立1階常微分方程式に直し、4次のルンゲ・クッタ法を適用すれば良いので ある。例えば、次のような3次の常微分方程式があったする。

$$y^{\prime\prime\prime}(x)=f(x,y,y^{\prime},y^{\prime\prime})$$ (20)

この3階常微分方程式を次に示す式を用いて変換する。

$$\left\{ \begin{aligned}y_0(x)&=y(x)\\ y_1(x)&=y^{\prime}(x)\\ y_2(x)&=y^{\prime\prime}(x) \end{aligned} \right.$$

この式を用いて、式(20)を書き直すと

$$\left\{ \begin{aligned}y_0^{\prime}(x)&=y_1(x)\\ y_1^{\prime}(x)&=y_2(x)\\ y_2^{\prime}(x)&=f(x,y_0,y_1,y_2) \end{aligned} \right.$$

となる。これで、3階の常微分方程式が3元の1階の連立常微分方程式に変換できた。2階で あろうが4階$\cdots$でも同じ方法で連立微分方程式に還元できる。

2.5.2 練習問題

以下の高次常微分方程式を連立1階微分方程式に書き換えなさい。

$$(1) \hspace{2mm} y^{\prime\prime}+3y^{\prime}+5y=0 \hspace{30mm} (2) \hspace{2mm} y^{\prime\prime}+6y^{\prime}+y=0$$

$$(3) 5y^{\prime\prime}+2xy^{\prime}+3y=0 (4) y^{\prime\prime\prime}+y^{\prime}+xy=0$$

$$(5) 5y^{\prime\prime}+y^{\prime}+y=\sin(\omega x) (6) xy^{\prime\prime}+y^{\prime}+y=e^{x}$$

$$(7) 5y^{\prime\prime}y^{\prime}+y^{\prime}+y=0 (8) y^{\prime\prime}y^{\prime}+x^2y^{\prime}y+y=0$$

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