ルンゲクッタ法の説明と実装例[C言語]

3日かけてルンゲクッタ法を理解する所までにきました．

とは言うものの，この記事では1階常微分方程式についての4次ルンゲクッタについて説明しています．

前回，前々回の記事である

・オイラー法と修正オイラー法の説明と実装例
・ニュートン・コーツの公式とシンプソンの公式の導出と実装

を読むと今回の記事で説明するルンゲクッタ法の導出部分が理解しやすいかと思います．

1.(4次)ルンゲクッタ法の簡単な説明

詳しい説明はしたところであまり興味がない人が多い気がするので，詳しい説明は最後にすっ飛ばして，実装するにおいて最低限の部分だけ示しておきます．

次の初期値がわかっている1階常微分方程式を解く事を考えます．
$$\begin{eqnarray*} f(x_0) &=& y_0 \\ \frac{df(x)}{dx} &=& g(x,f(x)) \end{eqnarray*}$$ この時，4次ルンゲクッタ法では次のようにして計算します． $$\int_{x_k}^{x_{k+1}} g(x,f(x)) dx \approx \frac{1}{6}(k_1+2k_2+2k_3+k_4 )\Delta x$$ ただしk1,...,k4は $$\begin{eqnarray*} k_1 &=& g(x_k,f(x_k)) \\ k_2 &=& g(x_{k+\frac{1}{2}},f(x_k) + k_1 \frac{\Delta x}{2}) \\ k_3 &=& g(x_{k+\frac{1}{2}},f(x_k) + k_2 \frac{\Delta x}{2}) \\ k_4 &=& g(x_{k+1},f(x_k)+k_3 \Delta x) \end{eqnarray*}$$ として計算しています．
この辺の導出は一番最後でしているので気になる人は見てってくださいな．

2.実装例

何を実装しようか悩んで適当にWikipediaを漁っていたらロジスティック方程式というモデルがあったので，これのシミュレーションをしてみたいと思います．

Wikipediaによれば，生物の個体数がどのように変化していくかを説明するモデルの一つであるとのこと．

求める常微分方程式はこれ． $$\frac{dN(t)}{dt} = r \frac{K-N(t)}{K} N(t)$$ rは一匹の個体の増加率，N(t)は時間tにおける生物の個体数，Kは個体数のしきい値で，最終的に安定するNの値．

例えば水槽なら最大でも20匹ぐらいしか入らない大きさならK=20みたいな値になると思う．

ちなみに厳密解は普通の変数分離法で多少計算が面倒だけど求まる．初期値をN_0とすると次のようになる． $$N(t) = \frac{N_0 K e^{rt}}{K - N_0 + N_0 e^{rt}}$$ 厳密解はさておき，この常微分方程式にルンゲクッタ法を適用する．

まず見通しを良くするために，1で示した式のg(x,f(x))に当たる部分g(t,N(t))を定義しておく． $$\begin{eqnarray*} \frac{dN(t)}{dt} &=& r \frac{K-N(t)}{K} N(t) \\ &=& g(t,N(t)) \end{eqnarray*}$$ そうすると$\Delta t$を微小な値として，$t_k = k \Delta t$とすると，漸化式はこのようになる． $$\int_{t_k}^{t_{k+1}} g(t,N(t)) dt \approx \frac{1}{6}(k_1+2k_2+2k_3+k_4 )\Delta t$$ ただしk1,...,k4は次の式で計算する． $$\begin{eqnarray*} k_1 &=& g(t_k,N(t_k)) \\ k_2 &=& g(t_{k+\frac{1}{2}},N(t_k) + k_1 \frac{\Delta t}{2}) \\ k_3 &=& g(t_{k+\frac{1}{2}},N(t_k) + k_2 \frac{\Delta t}{2}) \\ k_4 &=& g(t_{k+1},N(t_k)+k_3 \Delta t) \end{eqnarray*}$$ それではC言語による実装．関数g,を定義することで，数式とほぼ変わらないプログラムが書けます．実装の係数は適当なのでアレですが．

	/* 4次ルンゲクッタ法によるロジスティック方程式の解の計算 */
	#include <stdio.h>
	#define K 1000.0 // 個体数のしきい値
	#define r 0.9 // 1個体の増加率
	#define N_LOOP 10000 // ループ回数
	// 積分する関数g(t,N)
	double g(double t,double N_t){
	return r*(K-N_t)*N_t/K;
	}
	int main(void){
	int k; // ループカウンタ
	double dt = 1e-3; // 刻み幅Δt
	double N = 1.0; // 個体数t（初期値を代入しておく)
	double k1,k2,k3,k4; // 計算に使う諸変数
	double a,m,b; // t_k,t_k+0.5,t_k+1を格納する
	for( k = 0 ; k < N_LOOP ; k++ ){
	// t_k,t_k+0.5,t_k+1を計算
	a = k*dt;
	m = (k+0.5)*dt;
	b = (k+1.0)*dt;
	// k1,...,k4の値を計算
	k1 = g(a,N);
	k2 = g(m,N+k1*dt*0.5);
	k3 = g(m,N+k2*dt*0.5);
	k4 = g(b,N+k3*dt);
	// 積分して代入
	N = N + (dt/6.0)*(k1+2.0*k2+2.0*k3+k4);
	// 値を表示
	printf("%lf,%e\n",b,N);
	}
	return 0;
	}

view raw RK4_KogisticEquation.c hosted with ❤ by GitHub

シミュレーション結果．

このシミュレーションだとほとんど違いがわかりませんね．誤差をプロットするとこうなります．

オーダーが1e-5といったところでしょうか．誤差が振動しているのがわかります．シミュレーションの精度としては十分なのでしょう．

4次のルンゲクッタ法では4次の項まで級数展開が一致しているとかで，誤差は$O(\Delta t^5)$に比例して小さくなるとのことです．

流石にルンゲクッタ法でシミュレーションすれば精度どうなの？とか言われても大丈夫そう……？まだその辺はよく知らないけれど．

3.ルンゲクッタ法の説明

手元の本には書かれていないのではっきりしたことはわからないのですが軽く調べたところによると，修正オイラー法は2次ルンゲクッタ法とも呼ぶらしいです．

なので2次ルンゲクッタについては前回の記事を参照してください．ここでは4次ルンゲクッタについて説明します．

今回はとりあえず1階微分方程式についての4次ルンゲクッタについて説明するので，次の常微分方程式について考えます．
$$f(x_0) = y_0 \\ \frac{df(x)}{dx} = g(x,f(x))$$ ここで　$x_k = x_0 + k \Delta x$　と置くと　$f(x_{k+1})$　と　$f(x_{k})$　の漸化式を得ます． $$\begin{eqnarray*} \int_{x_k}^{x_{k+1}} \frac{df(x)}{dx} dx &=& f(x_{k+1})-f(x_k) \\ &=& \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \\ \\ ∴f(x_{k+1}) &=& f(x_k) + \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \end{eqnarray*}$$

この $\int_{x_k}^{x_{k+1}} g(x,f(x)) dx$ の部分をどう計算するかによってオイラー法，修正オイラー法，ルンゲクッタ法が分類できます．

オイラー法では矩形積分，修正オイラー法では台形積分，ルンゲクッタ法ではシンプソンの公式を使います．

シンプソンの公式は次のようになります． $$\int_{x_k}^{x_{k+1}}h(x) dx \approx \frac{1}{6}(h(x_k)+4h(x_k+\frac{\Delta x}{2}) + h(x_{k+1}) )\Delta x$$ これを使って積分値を計算すれば良いわけです．$x_{k+\frac{1}{2}} = x_k+\frac{1}{2}\Delta x \\ $ として $$\begin{eqnarray*} g_k &=& g(x_k,f(x_k)) \\ g_{k+\frac{1}{2}} &=& g(x_{k+\frac{1}{2}},f(x_{k+\frac{1}{2}})) \\ g_{k+1} &=& g(x_{k+1},f(x_{k+1})) \\ \end{eqnarray*}$$ という値を考えると， $\int_{x_k}^{x_{k+1}} g(x,f(x)) dx$ は次のように近似できます． $$\int_{x_k}^{x_{k+1}} g(x,f(x)) dx \approx \frac{1}{6}(g_k + 4g_{k+\frac{1}{2}} + g_{k+1} )\Delta x$$
しかしこの式には未知の値$f(x_{k+\frac{1}{2}})$や$f(x_{k+1})$を使っています．これらの値を矩形積分で求めて計算することもできますが，良い精度が得られません．

ルンゲクッタ法ではこれらの値を二段階の計算で計算します．
$$\begin{eqnarray*} k_1 &=& g(x_k,f(x_k)) \\ k_2 &=& g(x_{k+\frac{1}{2}},f(x_k) + k_1 \frac{\Delta x}{2}) \\ k_3 &=& g(x_{k+\frac{1}{2}},f(x_k) + k_2 \frac{\Delta x}{2}) \\ k_4 &=& g(x_{k+1},f(x_k)+k_3 \Delta x) \end{eqnarray*}$$ $g_{k+\frac{1}{2}}$に当たる部分をk2,k3という形で計算しています．これらを先ほどのシンプソンの公式に代入することで以下のような近似を行うものを4次のルンゲクッタ法と呼ぶそうです． $$\int_{x_k}^{x_{k+1}} g(x,f(x)) dx \approx \frac{1}{6}(k_1+2k_2+2k_3+k_4 )\Delta x$$

今日はここまで．