表2の出力に用いたCプログラム(文字コードはUTF-8)
/* Romberg quadrature */ #include <stdio.h> #include <math.h> #define NUMAX 4 double func(double x); int main(void) { int k,nu,pow4k; double s,w,x,xa,xb,h,integral; double t[NUMAX+1][NUMAX+1]; xa = 0.0; xb = 1.0; integral = exp(1.0) -1.0; nu = 0; h = xb - xa; s = 0.5 * h * (func(xa) + func(xb)); t[nu][0] = s; for(nu = 1;nu <= NUMAX;nu++){ h *= 0.5; x = xa + h; w = 0.0; while(x < xb){ w += func(x); x += 2 * h; } s = 0.5 * s + h * w; t[nu][0] =s; pow4k = 1; for(k = 1;k <= nu;k++){ pow4k *= 4; t[nu-k][k] = t[nu-k+1][k-1]+(t[nu-k+1][k-1]-t[nu-k][k-1])/(pow4k-1.0); } } for(nu = 0;nu <= NUMAX;nu++){ printf("%2d ",nu); for(k = 0;k <= nu;k++){ printf("%18.15f", t[nu-k][k]); } printf("\n"); } for(nu = 0;nu <= NUMAX;nu++){ printf("%2d ",nu); for(k = 0;k <= nu;k++){ printf("%10.3e", t[nu-k][k]-integral); } printf("\n"); } return 0; } double func(double x) { double y; y = exp(x); return y; }