表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;
}