pypy.com/

Python、Unity、FX自動化などを勉強しています。あと、コーラと車も好きです。そこらへんについて、たまに記事を書きます。

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

こんにちは。

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

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

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

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

taiboman.com - taiboman リソースおよび情報

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

ルンゲクッタ法の考え方

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

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

考え方は以下の図を使って説明したいと思います。
f:id:SaidaTaisei:20200506192206p:plain
まず、ある点(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);
}

結果をプロットしてみると、以下のようになりました。

f:id:SaidaTaisei:20200506225828p:plain

yのプロットがほとんど見えませんが、これはほぼ重なっているためです。

実際にx=3.0の時のyの値を比較すると、
厳密解が360.06853で、
ルンゲクッタ法で出した解が360.06825
となっていて、ほぼ同じ値になっていることがわかります。

ということで、今回はルンゲクッタ法をやってみました。

また、数値計算系の記事は書くと思うので、興味のある方は見てくださるとうれしいです。