Runge-Kutta-Fehlberg

rk-f.c

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

double f(double t, double x, double y, double ix) {
    const double k = 1.0;
    const double a = 1.0;
    const double c = 2.0;
    const double D = 0.5;
    return k * (a*x / (1 + a*x + c*ix)) - D * y;
}

int main(void) {
    double A[6][6] = {
        {    0.0       ,    0.0       ,    0.0       ,    0.0       ,    0.0       , 0.0},
        {    1.0/   4.0,    0.0       ,    0.0       ,    0.0       ,    0.0       , 0.0},
        {    3.0/  32.0,    9.0/  32.0,    0.0       ,    0.0       ,    0.0       , 0.0},
        { 1932.0/2197.0,-7200.0/2197.0, 7296.0/2197.0,    0.0       ,    0.0       , 0.0},
        {  439.0/ 216.0,-   8.0       , 3680.0/ 513.0,- 845.0/4104.0,    0.0       , 0.0},
        {-   8.0/  27.0,    2.0       ,-3544.0/2565.0, 1859.0/4104.0,-  11.0/  40.0, 0.0}
    };
    double cc[]  = { 0.0, 1.0/4.0, 3.0/8.0, 12.0/13.0, 1.0, 1.0/2.0 };
    double b4[] = { 25.0/216.0, 0.0, 1408.0/ 2565.0,  2197.0/ 4104.0, -1.0/ 5.0, 0.0    };
    double b5[] = { 16.0/135.0, 0.0, 6656.0/12825.0, 28561.0/56430.0, -9.0/50.0, 2.0/55 };

    const double thr = 1e-5;
    double t = 0.0, h = 1;
    double  a = 0.2,  b = 0.2,  c = 0.2;
    double  ia= 0.2, ib = 0.2,  ic= 0.1;

    while(t <= 300) {
        double  _da[6],  _db[6],  _dc[6];
        double _dia[6], _dib[6], _dic[6];

        for(int i=0; i<6; i++) {
            double  _a = 0,  _b = 0,  _c = 0;
            double _ia = 0, _ib = 0, _ic = 0;
            for(int j=0; j<i; j++) {
                _a  += h * A[i][j] * _da[j];
                _b  += h * A[i][j] * _db[j];
                _c  += h * A[i][j] * _dc[j];
                _ia += h * A[i][j] * _dia[j];
                _ib += h * A[i][j] * _dib[j];
                _ic += h * A[i][j] * _dic[j];
            }
            _a +=a, _ia += ia;
            _b +=b, _ib += ib;
            _c +=c, _ic += ic;

            _da[i]  = f(t + h * cc[i], _a,  _a, _ia);
            _db[i]  = f(t + h * cc[i], _b,  _b, _ib);
            _dc[i]  = f(t + h * cc[i], _c,  _c, _ic);
            _dia[i] = f(t + h * cc[i], _c, _ia, 0.0);
            _dib[i] = f(t + h * cc[i], _a, _ib, 0.0);
            _dic[i] = f(t + h * cc[i], _b, _ic, 0.0);
        }

        double da4 = 0, da5 = 0, db4 = 0, db5 = 0, dc4 = 0, dc5 = 0;
        double dia4 = 0, dia5 = 0, dib4 = 0, dib5 = 0, dic4 = 0, dic5 = 0;
        for(int i=0; i<6; i++) {
            da4  += h * b4[i] * _da[i];  da5  += h * b5[i] * _da[i];
            db4  += h * b4[i] * _db[i];  db5  += h * b5[i] * _db[i];
            dc4  += h * b4[i] * _dc[i];  dc5  += h * b5[i] * _dc[i];
            dia4 += h * b4[i] * _dia[i]; dia5 += h * b5[i] * _dia[i];
            dib4 += h * b4[i] * _dib[i]; dib5 += h * b5[i] * _dib[i];
            dic4 += h * b4[i] * _dic[i]; dic5 += h * b5[i] * _dic[i];
        }
        if(
                fabs(da4 - da5) > thr ||
                fabs(db4 - db5) > thr ||
                fabs(dc4 - dc5) > thr ||
                fabs(dia4 - dia5) > thr ||
                fabs(dib4 - dib5) > thr ||
                fabs(dic4 - dic5) > thr
          ) {
            h /= 2.0;
            continue;
        }

        t += h;
        a  += da5;
        b  += db5;
        c  += dc5;
        ia += dia5;
        ib += dib5;
        ic += dic5;

        printf("%.25f %.25f %.25f %.25f %.25f %.25f %.25f %.25f\n", t, h, a, b, c, ia, ib, ic);

        h *= 0.8 * pow( thr / fabs(da4 - da5), 0.2 ); 
    }

    return 0;
}