2 常微分方程式

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

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

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

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

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

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

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

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

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

2.2 プログラム(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.3 高階の常微分方程式

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

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

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

この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*}

この式を用いて,式(4)を書き直すと

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


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


no counter