Runge-Kutta法とは
Runge-Kutta法とは、常微分方程式における数値解析の1種である。Euler法では、シミュレーション時間を刻み幅で分割し、刻み幅を用いて数値解を算出していた。この刻み幅において、Euler法では演算を1回行い、その結果を次の数値解としていたが、Runge-Kutta法では、任意の回数だけ刻み幅の単位時間に演算を行い、次の数値解を算出している。従って、一般的にEuler法よりも精度は高くなる。しかし、次元を大きくしていくに従って、計算処理量が多くなってしまう。よって、Runge-Kutta法の演算次元は4次までとしていることが多い。今回、対象とする常微分方程式は、RC回路のコンデンサ電圧における回路方程式とします。この回路では、直流電源を用いるとします。
Cによるソースコード
#include<stdio.h>
#include<math.h>
double func(double t, double Vc) {
return 100000 - 1000 * Vc;
}
int main(void) {
double Vc, t, step, R, C, end;
double k1, k2, k3, k4;
int error;
FILE *fi;
R = 10;
C = 0.0001;
step = ((R*C) / 100.0)*100;
Vc = 0.0;
end = 0.01;
if ((error = fopen_s(&fi, "Runge_Kutta_data_100.txt", "w")) != 0) {
printf("\nThis file can't opened \n\n");
return 0;
}
for (t = 0.0; t <= end; t = t + step) {
fprintf(fi, "%f %f\n", t, Vc);
//Runge-Kutta法による演算
k1 = step*func(t, Vc);
k2 = step*func(t + step*0.5, Vc + k1*0.5);
k3 = step*func(t + step*0.5, Vc + k2*0.5);
k4 = step*func(t + step, Vc + k3);
Vc = Vc + (k1 + 2.0*k2 + 2.0*k3 + k4) / 6.0;
}
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