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