rk-f.c

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

double f_ddy(double x, double y, double dy) {
    return -3.0 * dy - 2.0 * y + cos(x);
}
double f_dy(double x, double y, double dy) {
    return dy;
}

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 c[]  = { 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-14;
    double x = 0.0, y = 0.0, dy = 0.0, h = 100;
    while(x <= 2.0 * M_PI) {
        double _ddy[6], _dy[6];
        for(int i=0; i<6; i++) {
            double __dy = 0, __y = 0;
            for(int j=0; j<i; j++) {
                __dy += h * A[i][j] * _ddy[j];
                __y += h * A[i][j] * _dy[j];
            }
            __dy += dy; __y += y;

            _ddy[i] = f_ddy(x + h * c[i], __y, __dy);
            _dy[i] = f_dy(x + h * c[i], __y, __dy);
        }

        double dy4 = 0, dy5 = 0, ddy4 = 0, ddy5 = 0;
        for(int i=0; i<6; i++) {
            dy4 += h * b4[i] * _dy[i];
            ddy4 += h * b4[i] * _ddy[i];
            dy5 += h * b5[i] * _dy[i];
            ddy5 += h * b5[i] * _ddy[i];
        }
        if( fabs(dy4 - dy5) > thr || fabs(ddy4 - ddy5) > thr ) {
            h /= 2.0;
            continue;
        }

        x += h;
        y += dy5;
        dy+= ddy5;

        printf("%.25f %.25f %.25f %.25f\n", x, y, dy, h);

        h *= 0.8 * pow( thr / fabs(dy4 - dy5), 0.2); 
    }

    return 0;
}