ルンゲクッタ法のサンプルコード

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

int main(void){
    double f,f1,f2,dx,x,k1,k2,k3,k4,k;

    f1 = 1;
    dx = 0.01;

    for (int i=0; i<=100; i++) {
        x = i*dx;
        k1 = dx*x*f1;
        k2 = dx*(x + dx/2.)*(f1 + k1/2.);
        k3 = dx*(x + dx/2.)*(f1 + k2/2.);
        k4 = dx*(x + dx)*(f1 + k3);
        k = (k1 + 2*k2 + 2*k3 + k4)/6.;

        f2 = f1 + k;
        f = exp(x*x/2.);

        printf("x = %f, 計算値 -> f1 = %f, 理論値 -> f = %f \n", x, f1, f);

        f1 = f2;
    }
    return 0;
}

出力結果

    x = 0.000000, 計算値 -> f1 = 1.000000, 理論値 -> f = 1.000000 
    x = 0.010000, 計算値 -> f1 = 1.000050, 理論値 -> f = 1.000050 
    x = 0.020000, 計算値 -> f1 = 1.000200, 理論値 -> f = 1.000200
    x = 0.030000, 計算値 -> f1 = 1.000450, 理論値 -> f = 1.000450
    x = 0.040000, 計算値 -> f1 = 1.000800, 理論値 -> f = 1.000800
    x = 0.050000, 計算値 -> f1 = 1.001251, 理論値 -> f = 1.001251 
    x = 0.060000, 計算値 -> f1 = 1.001802, 理論値 -> f = 1.001802
    x = 0.070000, 計算値 -> f1 = 1.002453, 理論値 -> f = 1.002453
    x = 0.080000, 計算値 -> f1 = 1.003205, 理論値 -> f = 1.003205
    x = 0.090000, 計算値 -> f1 = 1.004058, 理論値 -> f = 1.004058
    x = 0.100000, 計算値 -> f1 = 1.005013, 理論値 -> f = 1.005013
    x = 0.110000, 計算値 -> f1 = 1.006068, 理論値 -> f = 1.006068
    x = 0.120000, 計算値 -> f1 = 1.007226, 理論値 -> f = 1.007226
    x = 0.130000, 計算値 -> f1 = 1.008486, 理論値 -> f = 1.008486
    x = 0.140000, 計算値 -> f1 = 1.009848, 理論値 -> f = 1.009848
    x = 0.150000, 計算値 -> f1 = 1.011314, 理論値 -> f = 1.011314
    x = 0.160000, 計算値 -> f1 = 1.012882, 理論値 -> f = 1.012882
    x = 0.170000, 計算値 -> f1 = 1.014555, 理論値 -> f = 1.014555
    x = 0.180000, 計算値 -> f1 = 1.016332, 理論値 -> f = 1.016332
    x = 0.190000, 計算値 -> f1 = 1.018214, 理論値 -> f = 1.018214
    x = 0.200000, 計算値 -> f1 = 1.020201, 理論値 -> f = 1.020201
    x = 0.210000, 計算値 -> f1 = 1.022295, 理論値 -> f = 1.022295 
    x = 0.220000, 計算値 -> f1 = 1.024495, 理論値 -> f = 1.024495
    x = 0.230000, 計算値 -> f1 = 1.026803, 理論値 -> f = 1.026803
    x = 0.240000, 計算値 -> f1 = 1.029219, 理論値 -> f = 1.029219
    x = 0.250000, 計算値 -> f1 = 1.031743, 理論値 -> f = 1.031743
    x = 0.260000, 計算値 -> f1 = 1.034378, 理論値 -> f = 1.034378
    x = 0.270000, 計算値 -> f1 = 1.037122, 理論値 -> f = 1.037122
    x = 0.280000, 計算値 -> f1 = 1.039978, 理論値 -> f = 1.039978
    x = 0.290000, 計算値 -> f1 = 1.042947, 理論値 -> f = 1.042947
    x = 0.300000, 計算値 -> f1 = 1.046028, 理論値 -> f = 1.046028
    x = 0.310000, 計算値 -> f1 = 1.049223, 理論値 -> f = 1.049223
    x = 0.320000, 計算値 -> f1 = 1.052533, 理論値 -> f = 1.052533
    x = 0.330000, 計算値 -> f1 = 1.055960, 理論値 -> f = 1.055960
    x = 0.340000, 計算値 -> f1 = 1.059503, 理論値 -> f = 1.059503
    x = 0.350000, 計算値 -> f1 = 1.063165, 理論値 -> f = 1.063165 
    x = 0.360000, 計算値 -> f1 = 1.066946, 理論値 -> f = 1.066946
    x = 0.370000, 計算値 -> f1 = 1.070847, 理論値 -> f = 1.070847
    x = 0.380000, 計算値 -> f1 = 1.074870, 理論値 -> f = 1.074870
    x = 0.390000, 計算値 -> f1 = 1.079017, 理論値 -> f = 1.079017
    x = 0.400000, 計算値 -> f1 = 1.083287, 理論値 -> f = 1.083287
    x = 0.410000, 計算値 -> f1 = 1.087683, 理論値 -> f = 1.087683
    x = 0.420000, 計算値 -> f1 = 1.092207, 理論値 -> f = 1.092207
    x = 0.430000, 計算値 -> f1 = 1.096858, 理論値 -> f = 1.096858
    x = 0.440000, 計算値 -> f1 = 1.101640, 理論値 -> f = 1.101640
    x = 0.450000, 計算値 -> f1 = 1.106553, 理論値 -> f = 1.106553
    x = 0.460000, 計算値 -> f1 = 1.111600, 理論値 -> f = 1.111600
    x = 0.470000, 計算値 -> f1 = 1.116781, 理論値 -> f = 1.116781
    x = 0.480000, 計算値 -> f1 = 1.122098, 理論値 -> f = 1.122098
    x = 0.490000, 計算値 -> f1 = 1.127553, 理論値 -> f = 1.127553
    x = 0.500000, 計算値 -> f1 = 1.133148, 理論値 -> f = 1.133148
    x = 0.510000, 計算値 -> f1 = 1.138885, 理論値 -> f = 1.138885
    x = 0.520000, 計算値 -> f1 = 1.144766, 理論値 -> f = 1.144766
    x = 0.530000, 計算値 -> f1 = 1.150792, 理論値 -> f = 1.150792
    x = 0.540000, 計算値 -> f1 = 1.156965, 理論値 -> f = 1.156965
    x = 0.550000, 計算値 -> f1 = 1.163287, 理論値 -> f = 1.163287
    x = 0.560000, 計算値 -> f1 = 1.169762, 理論値 -> f = 1.169762
    x = 0.570000, 計算値 -> f1 = 1.176389, 理論値 -> f = 1.176389
    x = 0.580000, 計算値 -> f1 = 1.183173, 理論値 -> f = 1.183173
    x = 0.590000, 計算値 -> f1 = 1.190115, 理論値 -> f = 1.190115
    x = 0.600000, 計算値 -> f1 = 1.197217, 理論値 -> f = 1.197217
    x = 0.610000, 計算値 -> f1 = 1.204482, 理論値 -> f = 1.204482
    x = 0.620000, 計算値 -> f1 = 1.211913, 理論値 -> f = 1.211913
    x = 0.630000, 計算値 -> f1 = 1.219511, 理論値 -> f = 1.219511
    x = 0.640000, 計算値 -> f1 = 1.227280, 理論値 -> f = 1.227280
    x = 0.650000, 計算値 -> f1 = 1.235221, 理論値 -> f = 1.235221
    x = 0.660000, 計算値 -> f1 = 1.243338, 理論値 -> f = 1.243338 
    x = 0.670000, 計算値 -> f1 = 1.251634, 理論値 -> f = 1.251634
    x = 0.680000, 計算値 -> f1 = 1.260111, 理論値 -> f = 1.260111
    x = 0.690000, 計算値 -> f1 = 1.268773, 理論値 -> f = 1.268773
    x = 0.700000, 計算値 -> f1 = 1.277621, 理論値 -> f = 1.277621
    x = 0.710000, 計算値 -> f1 = 1.286660, 理論値 -> f = 1.286660 
    x = 0.720000, 計算値 -> f1 = 1.295893, 理論値 -> f = 1.295893
    x = 0.730000, 計算値 -> f1 = 1.305322, 理論値 -> f = 1.305322
    x = 0.740000, 計算値 -> f1 = 1.314952, 理論値 -> f = 1.314952
    x = 0.750000, 計算値 -> f1 = 1.324785, 理論値 -> f = 1.324785
    x = 0.760000, 計算値 -> f1 = 1.334825, 理論値 -> f = 1.334825
    x = 0.770000, 計算値 -> f1 = 1.345075, 理論値 -> f = 1.345075
    x = 0.780000, 計算値 -> f1 = 1.355540, 理論値 -> f = 1.355540
    x = 0.790000, 計算値 -> f1 = 1.366223, 理論値 -> f = 1.366223
    x = 0.800000, 計算値 -> f1 = 1.377128, 理論値 -> f = 1.377128
    x = 0.810000, 計算値 -> f1 = 1.388258, 理論値 -> f = 1.388258
    x = 0.820000, 計算値 -> f1 = 1.399619, 理論値 -> f = 1.399619
    x = 0.830000, 計算値 -> f1 = 1.411214, 理論値 -> f = 1.411214
    x = 0.840000, 計算値 -> f1 = 1.423047, 理論値 -> f = 1.423047
    x = 0.850000, 計算値 -> f1 = 1.435122, 理論値 -> f = 1.435122
    x = 0.860000, 計算値 -> f1 = 1.447445, 理論値 -> f = 1.447445
    x = 0.870000, 計算値 -> f1 = 1.460020, 理論値 -> f = 1.460020
    x = 0.880000, 計算値 -> f1 = 1.472851, 理論値 -> f = 1.472851
    x = 0.890000, 計算値 -> f1 = 1.485944, 理論値 -> f = 1.485944
    x = 0.900000, 計算値 -> f1 = 1.499303, 理論値 -> f = 1.499303
    x = 0.910000, 計算値 -> f1 = 1.512933, 理論値 -> f = 1.512933
    x = 0.920000, 計算値 -> f1 = 1.526840, 理論値 -> f = 1.526840
    x = 0.930000, 計算値 -> f1 = 1.541028, 理論値 -> f = 1.541028
    x = 0.940000, 計算値 -> f1 = 1.555505, 理論値 -> f = 1.555505
    x = 0.950000, 計算値 -> f1 = 1.570274, 理論値 -> f = 1.570274
    x = 0.960000, 計算値 -> f1 = 1.585342, 理論値 -> f = 1.585342
    x = 0.970000, 計算値 -> f1 = 1.600714, 理論値 -> f = 1.600714
    x = 0.980000, 計算値 -> f1 = 1.616398, 理論値 -> f = 1.616398
    x = 0.990000, 計算値 -> f1 = 1.632398, 理論値 -> f = 1.632398
    x = 1.000000, 計算値 -> f1 = 1.648721, 理論値 -> f = 1.648721