2.5.2 Runge-Ketta 法

runge.c

/* runge.c */

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

int main()
{
    double a = 0.0, b = 1.0;
    double x0;
    double t, x, h, f(double, double);
    double k1,k2,k3,k4;
    int i, N;

    printf(" x0="); scanf("%lf", &x0);
    printf(" N="); scanf("%d", &N);

    h = (b - a) / N;

    t = a; x = x0;
    printf("%g %g\n", t, x);

    for (i = 0; i < N; i++) {
        k1 = h * f(t,x);
        k2 = h * f(t + h / 2, x + k1 / 2);
        k3 = h * f(t + h / 2, x + k2 / 2);
        k4 = h * f(t + h, x + k3);
        t = t + h;
        x = x + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
        printf("%g %g\n", t, x);
    }
    printf("%g %20.15e\n", t, x);
    return 0;
}

double f(double t, double x)
{
    return x;
}

実行結果

mathpc00% ./runge 
 x0=1
 N=10
0 1
0.1 1.10517
0.2 1.2214
0.3 1.34986
0.4 1.49182
0.5 1.64872
0.6 1.82212
0.7 2.01375
0.8 2.22554
0.9 2.4596
1 2.71828
1 2.718279744135166e+00
mathpc00%