Back to the original topic:
I used a C implementation of the Romberg algorithm to find the integral of
f(x) = 1 / sqrt(1 - x) over the interval 0 - 1
It runs for a relatively long time (several seconds on my 3 GHz Intel i5) and it came up with this result:
result = 2.000052194143781214563660
function evaluations = 1073741825 (1 billion of function evaluations!)
I set an accuracy of 0.0001 with max steps = 100. I also used the interval from 0 to 0.99999999999 or else it would evaluate to NAN.
If i set an accuracy of 0.01 I get this:
result = 2.005444550190635499831160
function evaluations = 16777217
Code:
#include <stdio.h>
#include <math.h>
void
dump_row(size_t i, double *R){
printf("R[%2zu] = ", i);
for(size_t j = 0; j <= i; ++j){
printf("%.12f ", R[j]);
}
printf("\n");
}
double
romberg(double (*f/* function to integrate */)(double), double /*lower limit*/ a, double /*upper limit*/ b, size_t max_steps, double /*desired accuracy*/ acc){
double R1[max_steps], R2[max_steps]; //buffers
double *Rp = &R1[0], *Rc = &R2[0]; //Rp is previous row, Rc is current row
double h = (b-a); //step size
Rp[0] = (f(a) + f(b))*h*.5; //first trapezoidal step
dump_row(0, Rp);
for(size_t i = 1; i < max_steps; ++i){
h /= 2.;
double c = 0;
size_t ep = 1 << (i-1); //2^(n-1)
for(size_t j = 1; j <= ep; ++j){
c += f(a+(2*j-1)*h);
}
Rc[0] = h*c + .5*Rp[0]; //R(i,0)
for(size_t j = 1; j <= i; ++j){
double n_k = pow(4, j);
Rc[j] = (n_k*Rc[j-1] - Rp[j-1])/(n_k-1); //compute R(i,j)
}
//Dump ith column of R, R[i,i] is the best estimate so far
dump_row(i, Rc);
if(i > 1 && fabs(Rp[i-1]-Rc[i]) < acc){
return Rc[i-1];
}
//swap Rn and Rc as we only need the last row
double *rt = Rp;
Rp = Rc;
Rc = rt;
}
return Rp[max_steps-1]; //return our best guess
}
int functionCounter = 0;
double myfunc(double x) {
functionCounter++;
double l = 1 / sqrt( 1 - x );
return l;
}
int main() {
printf("result = %.24f\nfunction evaluations = %d\n", romberg(&myfunc, 0, 0.99999999999, 100, 0.0001), functionCounter);
}
The code source is from wikipedia:
https://en.wikipedia.org/wiki/Romberg%27s_method