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

関連記事