c言語で微分方程式を解く（ルンゲクッタ法）

こんにちは。

今日はc言語を使って、ルンゲクッタ法で微分方程式を解きたいと思います。

過去二回、それぞれオイラー法と修正オイラー法で解いてみましたが、ルンゲクッタ法はそれよりも精度の良い方法となります。

前回の修正オイラー法でも、それなりの精度での計算が可能でしたが、さらに精度が必要な場合はルンゲクッタ法を使うことになります。

オイラー法と修正オイラー法の記事は以下から見ることができるので、興味のある方は是非ご覧ください。

c言語で微分方程式を解く（オイラー法） - pypy.com/

c言語で微分方程式を解く（修正オイラー法） - pypy.com/

ルンゲクッタ法の考え方
実際にc言語でルンゲクッタ法のコードを書いてみる

ルンゲクッタ法の考え方

まず、以下の式を考えます。

$\frac{dy}{dx}=f(x,y)$

考え方は以下の図を使って説明したいと思います。

まず、ある点（xi,yi）の傾きを考え、図のk1を求めます。
$k_1=hf(x_i,y_i)$

次にk2を求めます。
これは、緑線のちょうど真ん中の点を使って求めます。
$k_2=hf(x_i+\frac{h}{2},y_i+\frac{k_1}{2})$

同様の考え方で、k3を求めます。
$k_3=hf(x_i+\frac{h}{2},y_i+\frac{k_2}{2})$

さらに、k4を求めます。
これは、赤の矢印の傾きを用いて考えます。
$k_4=hf(x_i+h,y_i+k_3)$

このk1、k2、k3、k4を用いて $x_{i+1}=x_{i}+h$ の時のy座標を以下のように推定します。
$y_{i+1}=y_{i}+\frac{1}{6}(k_1+2k_2+2k_3+k_4)$

これを繰り返すことで、微分方程式を解くことができます。

実際にc言語でルンゲクッタ法のコードを書いてみる

それでは、実際にコードを書いてみます。

今回は、前回の修正オイラー法と同じように、
\begin{align}　
\frac{dy}{dx} = xy　\tag{1}
\end{align}
を初期値 $(x_0,y_0)=(0.0,4.0)$
分割数100で解きたいと思います。

なお、(1)式の厳密解は
\begin{align}　
y = 4e^ \frac{x^2}{2}　
\end{align}
です。

厳密解との比較もしたいと思います。

#include <stdio.h>
#include <math.h>

#define N 100 //分割数

double f(double x, double y);
double y_true(double x);

int main(void) {
	//初期値
	double x = 0.0;
	double y = 4.0;
	double h = (3.0 - 0.0) / (double)N; //刻み幅

	printf("x       y       y_true\n");
	for (double i = 0; i < N; i++) {
		double k1 = h * f(x, y);
		double k2 = h * f(x + h/2.0, y + k1/2.0);
		double k3 = h * f(x + h / 2.0, y + k2 / 2.0);
		double k4 = h * f(x + h, y + k3);
		y = y + (k1 + 2.0*k2 + 2.0*k3 + k4) / 6.0;
		x = x + h;
		printf("%5.5lf, %5.5lf, %5.5lf\n", x, y, y_true(x));
	}
}

double f(double x, double y) {
	return x * y;
}
double y_true(double x) {
	return 4.0*exp(x*x / 2.0);
}