Runge-Kutta-Fehlberg法とは

前回にRunge-Kutta法による数値解析を投稿しました。また、前々回においては、Euler法による数値解析を投稿しました。これらの方法のでは、シミュレーション時間を固定幅の刻み幅で分割し、分割した単位時間をそれぞれの方式の演算区間としていました。対して、Runge-Kutta-Fehlberg法では、シミュレーション時間の刻み幅を可変にすることが出来ます。これにより、数値の変化が基準より少ない場合には、刻み幅を大きくし、数値の変化が基準値より大きい場合には、刻み幅を小さくする制御として使用することが出来ます。従って、Runge-Kutta法よりも精度が高くなります。今回においても対象とする回路は、直流電源のRC回路とし、入力を電源電圧、出力をコンデンサ電圧とします。次にC言語によるソースコードを示します。

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

double func(double t, double Vc) {

	return 100000 - 1000 * Vc;

}

int main(void) {

	//変数設定//

	long double Vc, t, step, R, C, end, r;
	long double k1, k2, k3, k4, k5, k6;
	int FLAG = 0;
	long double step_MIN, step_MAX, tol, e, Vc_o5, Vc_o4, diff;
	int error;
	FILE *fi;

	//変数初期化//
	R = 10;
	C = 0.0001;
	r = 0.0;
	step = ((R*C) / 100.0) * 100;
	//step = 0.00001; //step数は関数に応じて任意に設定する
	diff = 0.0;
	Vc = 0.0;
	Vc_o4 = 0.0;
	Vc_o5 = 0.0;
	step_MIN = 0.00001;
	step_MAX = 0.005;
	//tol = 2.059999; //偏差の基準は関数に応じて任意に決定する
	tol = 2.06;
	e = 0;
	t = 0.0;
	end = 0.01;

	//textファイルの作成//
	if ((error = fopen_s(&fi, "Runge_Kutta_Fehlberg_data.txt", "w")) != 0) {

		printf("\nThis file can't opened \n\n");

		return 0;

	}

	//初期値の書き込み//
	fprintf(fi, "%f %f\n", t, Vc);

	//Runge-Kutta-Fehlberg法の演算//
	while (FLAG != 1) {
		
	/*	k1からk6までを求める演算では、コメントアウトしている式が正しい式なのですが、
		分数のまま演算すると後で0になってしまうので、わざわざ小数にしています*/

		k1 = step*func(t, Vc); printf("k1 = %lf\n", k1);
		k2 = step*func(t + step*0.25, Vc + (k1*0.25)); printf("k2 = %lf\n", k2);
		//k3 = step*func(t + step*(3/8), Vc + ( (k1*(3/32)) + (k2*(9/32)))); 
		k3 = step*func(t + step*(0.375), Vc + ((k1*(0.09375)) + (k2*(0.28125))));
		//k4 = step*func(t + step*(12/13), Vc +  (((1932/2197)*k1) - ((7200/2197)*k2) + ((7296/2197) * k3) )); 
		k4 = step*func(t + step*(0.92307), Vc + (((0.87938)*k1) - ((3.27719)*k2) + ((3.320892) * k3))); printf("k4 = %lf\n", k4);
		//k5 = step*func(t + step, Vc +  (((439 / 216)*k1) - (8 * k2) + ((3680 / 513)*k3) - ((845 / 4104)*k4)));
		k5 = step*func(t + step, Vc + (((2.032407)*k1) - (8 * k2) + ((7.173489)*k3) - ((0.20589)*k4))); printf("k5 = %lf\n", k5);
		//k6 = step*func(t + step*0.5, Vc + (-((8 / 27)*k1) + (2 * k2) - ((3544 / 2565)*k3) + ((1859 / 4104)*k4) - ((11/40)*k5)));
		k6 = step*func(t + step*0.5, Vc + (-((0.29629)*k1) + (2 * k2) - ((1.38167)*k3) + ((0.45330)*k4) - ((0.275)*k5))); printf("k6 = %lf\n", k6);

		Vc_o5 = ((0.1185 * k1) + (0.51898 * k3) + (0.50613 * k4) - (0.18 * k5) + (0.03636 * k6));
		//Vc_o5 = (((16 / 135) * k1) + ((6656 / 12825) * k3) + ((28561 / 56430) * k4) - ((9 / 50) * k5) + ((2 / 55) * k6));
		//Vc_o4 = Vc + (((25 / 216)*k1) + ((1408 / 2565)*k3) + ((2197 / 4104)*k4) - (1 / 5)*k5);
		Vc_o4 = ((0.11574 * k1) + (0.54892 * k3) + (0.53533 * k4) - (0.2)*k5);

		diff = Vc_o5 - Vc_o4;

		if (diff < 0) { 

			diff = diff*(-1);

		}

		r = (diff / step); printf("r = %lf\n", r);

		if (r <= tol) {

			t = t + step;
			//Vc = Vc + (25/216)*k1 + (1408/2565)*k3 + (2197/4104)*k4 - (1/5)*k5;
			Vc = Vc + (0.11574) * k1 + (0.54892)*k3 + (0.535722)*k4 - (0.2)*k5;
			fprintf(fi, "%f %f\n", t, Vc);

		}

		e = pow((tol / (2 * r)), 0.25);

		//次のif文の基準は、関数に応じて任意に設定する
		if (e <= 0.1) {

			step = step * 0.1;  //step幅を小さくする

		}
		else if (e >= 1.5) {

			//step = step * 100000;
			step = step * 2; //step幅を大きくする

		}
		else {

			step = e*step;

		}

		if (step > step_MAX) {

			step = step_MAX;

		}

		printf("t = %lf\n", t);

		if (t >= end) {

			FLAG = 1;

		}
		else if ((t + step) > end) {

			step = 0.01 - t;

			if (step < step_MIN) {

				FLAG = 1;
				printf("\nMINIMUM step exceeded !!\n");

			}
		}
	}

	//textファイルの保存//
	fclose(fi);

	return 0;
}

上記のソースコードを実行すると、時間領域におけるRC回路のコンデンサ電圧の値がtextファイルとして出てきます。
この出てきたファイルをgunuplotで実行させるとRC回路のコンデンサ電圧におけるEuler法を用いた数値解をプロットする事が出来ます。

以下に、簡単なGunuplotのソースコードを示します。


se grid
se xrange [0:0.01]
se yrange [0:110]
se xlabel "s[sec]"
se ylabel "Vc[V]"
plot "Euler_data.txt" using 1:2 with lines