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