[Previous: 1 後期中間試験の傾向と対策] [Up: 常微分方程式と連立方程式，補間法のまとめ] [Next: 3 連立一次方程式(消去法)]

2 常微分方程式

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

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

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

2.1 テイラー展開

数値計算は言うに及ばず科学技術全般の考察にテイラー展開(Taylor expansion)は，重要 な役割を果たす．電気の諸問題を考察する場合，いたるところにテイラー展開は顔を出す ので，十分理解しなくてはならない．まずは，$ x=a$の回りのテイラー展開は，

\begin{align*}\begin{aligned}f(x)&=f(a)+f^\prime(a)(x-a)+\frac{1}{2!}f^{\prime\p... ...ots\\ &=\sum_{n=0}^\infty\frac{1}{n!}f^{(n)}(a)(x-a)^n \end{aligned}\end{align*} (2)

と書ける．0の階乗は$ 0!=1$となることに注意． $ f^{(n)}(a)$は，関 数$ f(x)$を$ n$回微分したときの$ x=a$の値である．テイラー展開は，任意の関数 $ f(x)$は，無限のべき級数に展開できると言っている．その係数が， $ \frac{1}{n!}f^{(n)}(a)$である．次に，$ a=0$としてみる．この場合，

\begin{align*}\begin{aligned}f(x)&=f(0)+f^\prime(0)x+\frac{1}{2!}f^{\prime\prime... ...+\cdots\\ &=\sum_{n=0}^\infty\frac{1}{n!}f^{(n)}(0)x^n \end{aligned}\end{align*} (3)

となる．これがマクローリン展開(Maclaurin expansion)である．次に $ x-a=\Delta x$とする．そして，$ a=x_0$とすると

\begin{align*}\begin{aligned}f(x_0+\Delta x)&=f(x_0)+f^\prime(x_0)\Delta x +\fra... ... &=\sum_{n=0}^\infty\frac{1}{n!}f^{(n)}(x_0)\Delta x^n \end{aligned}\end{align*} (4)

となる．これは，しばしばお目にかかるパターン．

通常，我々がテイラー展開を使う場合，この式 (2)〜式(4) のいずれかである．

2.2 オイラー法

2.2.1 計算原理

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

$\displaystyle 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$ (5)

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

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

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

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

を計算する．この場合，計算の精度は1次と言う2．

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

\begin{equation*}\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.\end{equation*}

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

図 1: オイラー法．ある区間での$ y$の変化$ \Delta y$は，計算の始めの 点の傾きに区間の幅$ \Delta x$を乗じて，求めている．
\includegraphics[keepaspectratio, scale=0.7]{figure/diff_eq/Euler.eps}

2.2.2 プログラム例

次の微分方程式

  $\displaystyle \if 11 \frac{\mathrm{d}y}{\mathrm{d}x} \else \frac{\mathrm{d}^{1} y}{\mathrm{d}x^{1}}\fi =x^2\sin x+\sqrt{y}\cos y$   $\displaystyle y(0)=0$   (9)

の近似解をオイラー法で計算するプログラムをリスト1に示す．計算結果 は，ファイルに格納している．計算領域は $ 0\leqq x \leqq 2$，計算のステップ幅は $ \mathrm{d}x = 0.001$となっている． 
   1 #include <stdio.h>
   2 #include <math.h>
   3 #define NSTEPS 2000                // 計算ステップ数
   4 #define NDIM 30000               // 用意する配列の大きさ
   5 #define X_MAX 2.0                // 計算する最大 x
   6 
   7 double f(double x, double y);    // プロトタイプ宣言
   8 
   9 //============================================================
  10 // main 関数
  11 //============================================================
  12 int main(void){
  13   double x[NDIM], y[NDIM];      // 計算結果を入れる
  14   double dx;
  15   int i;
  16   FILE *out;                    // 計算結果を保存するファイル
  17   
  18 
  19   dx = X_MAX/NSTEPS;            // 計算のきざみ幅dxの計算
  20   x[0] = 0;
  21   y[0] = 0;
  22 
  23 
  24   //------- オイラー法の計算 --------
  25 
  26   for(i=0; i<NSTEPS; i++){
  27     x[i+1] = x[0]+(i+1)*dx;       //x[i+1]=x[i]+dxよりも精度が良い
  28     y[i+1] = y[i]+f(x[i],y[i])*dx;
  29   }
  30 
  31 
  32   //------- 計算結果をファイルに保存 --------
  33 
  34   out = fopen("result.txt","w");
  35   for(i=0; i<=NSTEPS; i++){
  36     fprintf(out,"%20.15f\t%20.15f\n", x[i],y[i]);
  37   }
  38 
  39 
  40   fclose(out);
  41 
  42   return 0;
  43 }
  44 
  45 //============================================================
  46 // dy/dx = f(x,y) の　f(x,y)を計算する関数
  47 //============================================================
  48 double f(double x, double y){
  49   return x*x*sin(x)+sqrt(y)*cos(y);
  50 }

2.3 2次のルンゲクッタ

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

2.3.1 ホイン法

先に示したように，オイラー法の精度は1次であるが，2次のルンゲ・ クッタ法では2次となる．今まで刻み幅を$ \Delta x$と記述していたが，簡単のため$ h$と 表現する．

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

$\displaystyle 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)$ (10)

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

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

となっている．

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

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

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

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

となる．これを，式(11)と比較すると，

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

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

\begin{equation*}\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.\end{equation*}

この式は，図2のようになる．オイラー法の図 1との比較でも，精度が良いことが分かる．
図 2: ホイン法．ある区間での$ y$の変化$ \Delta y$は， 計算の始めと終わりの点付近の平均傾きに区間の幅$ h$を乗じて，求 めている．
\includegraphics[keepaspectratio, scale=0.7]{figure/diff_eq/RK2_1.eps}

2.3.2 中点法

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

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

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

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

となる．これを，式(11)と比較すると，

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

となる必要がある．したがって，中点法の漸化式は，

\begin{equation*}\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.\end{equation*}

となる．この公式のイメージを，図3に示しておく．
図 3: 中点法．ある区間での$ y$の変化$ \Delta y$は，中点付近の傾きに 区間の幅$ h$を乗じて，求めている．
\includegraphics[keepaspectratio, scale=1.0]{figure/diff_eq/RK2_2.eps}

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

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

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

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

\begin{equation*}\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.\end{equation*}

図 4: 4次のルンゲ・クッタ法．ある区間での$ y$の変化$ \Delta y$は，区 間内の4点の傾きのある種の加重平均に幅$ h$を乗じて，求めている．
\includegraphics[keepaspectratio, scale=0.7]{figure/diff_eq/RK4.eps}

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

実際の微分方程式

\begin{equation*}\left\{ \\ \begin{aligned}&\frac{dy}{dx}=\sin x \cos x -y\cos x\\ &y=0\quad(初期条件 x=0の時) \end{aligned} \right.\end{equation*}

を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.6 高階の常微分方程式

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

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

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

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

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

この式を用いて，式(22)を書き直すと

\begin{equation*}\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.\end{equation*}

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

2.6.2 練習問題

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

  $\displaystyle (1)$ $\displaystyle \hspace{2mm}$ $\displaystyle y^{\prime\prime}+3y^{\prime}+5y=0$ $\displaystyle \hspace{30mm}$ $\displaystyle (2)$ $\displaystyle \hspace{2mm}$ $\displaystyle y^{\prime\prime}+6y^{\prime}+y=0$    
  $\displaystyle (3)$   $\displaystyle 5y^{\prime\prime}+2xy^{\prime}+3y=0$   $\displaystyle (4)$   $\displaystyle y^{\prime\prime\prime}+y^{\prime}+xy=0$    
  $\displaystyle (5)$   $\displaystyle 5y^{\prime\prime}+y^{\prime}+y=\sin(\omega x)$   $\displaystyle (6)$   $\displaystyle xy^{\prime\prime}+y^{\prime}+y=e^{x}$    
  $\displaystyle (7)$   $\displaystyle 5y^{\prime\prime}y^{\prime}+y^{\prime}+y=0$   $\displaystyle (8)$   $\displaystyle y^{\prime\prime}y^{\prime}+x^2y^{\prime}y+y=0$    


[Previous: 1 後期中間試験の傾向と対策] [Up: 常微分方程式と連立方程式，補間法のまとめ] [Next: 3 連立一次方程式(消去法)]
[Top:ページの先頭に戻る]

ホームページ: Yamamoto's laboratory
著者: 山本昌志
Yamamoto Masashi
平成19年11月27日

no counter