I was asked to use Romberg integration to evaluate the integral$$\int_0^1x^{-x}dx=\sum_{n=1}^\infty n^{-n}$$ and compare the result with the result I get from the sum. And I also need to estimate how many function evaluation Romberg integration will require to achieve 12 digit accuracy. I looked it up on Wikipedia but I still don't know how to estimate it using Romberg integration. And I guess that to achieve the accuracy, we could do that by looking into the error? Thanks!
This is not really an answer, but it is way too long for a comment. Here is an example of Romberg-type quadrature:
The base rule is a left hand Riemann sum on a uniform partition with subintervals of length $h=1$. The problem of interest is to integrate $e^x$ on $[0,1]$ with a tolerance of $10^{-1}$. So the final result should be somewhere between about $1.62$ and about $1.82$.
Our initial estimate is $1$, by implementing the base rule with $h=1$. We consider a refinement to $h=1/2$. This refinement would give an estimate of $\frac{1+e^{1/2}}{2} \approx 1.3$. This is too far from our initial estimate. The clever thing about Romberg integration is that instead of just using the value of the refinement, we can do a kind of weighted average of these two results to get an even better estimate. All we need to know is the order of our method (1 in this case). The idea is called Richardson extrapolation: if you have a function $T(h)$ representing the approximation of the desired value $I$ given by your method, then for a first order (convergent) method you have the expansion
$$T_1(h)=I + a_1 h + o(h)$$
where $h$ is the step size. Therefore
$$T_1(h/2)=I + a_1 h/2 + o(h)$$
Now $T_1(h)-2T_1(h/2)$ cancels out the linear terms, leaving $-I+o(h)$. Switching the sign gives $I+o(h)$, so we have a higher order method, namely $T_2(h):=2T_1(h/2)-T_1(h)$.
Returning to our example, this means we should try the estimate
$$2(e^0 \cdot 1/2 + e^{1/2} \cdot 1/2)-(e^0 \cdot 1)=e^0+e^{1/2}-e^0=e^{1/2} \approx 1.64.$$
$T_2$ is actually the midpoint rule, when simplified, but you don't need to know that (and in general such nice simplification doesn't occur).
To extrapolate further, we need to know the order of the midpoint rule, which is 2. So we look at Richardson extrapolation for order 2:
$$T_2(h)=I+a_2h^2+o(h^2) \\ T_2(h/2)=I+a_2h^2/4+o(h^2).$$
Therefore $T_2(h)-4T_2(h/2)=-3I+o(h^2)$, so $T_3(h):=\frac{4T_2(h/2)-T_2(h)}{3}$ is an even higher order method.
Now we need to calculate $T_3(1)$ in order to know whether to stop. We already took $T_2(1)$, it was $e^{1/2} \approx 1.6$. Now $T_2(1/2)$ is $e^{1/4} \cdot 1/2 + e^{3/4} \cdot 1/2 \approx 1.700$. So $T_3(1)$ is $\frac{2 e^{1/4} + 2 e^{3/4} - e^{1/2}}{3} \approx 1.718$. This is within $10^{-1}$ of our previous answer (about $1.64$), so we call this our final result.
I find it easier to work with this by just rederiving it through the Richardson extrapolation procedure rather than trying to learn any formulae.
By the way, this presentation is not adaptive. You get an adaptive variant by separately estimating the error on subintervals, and only using the higher order method on a given subinterval when the error estimate for that subinterval alone is too high.
Romberg integration involves successive subdivisions of the region of integration combined with extrapolation (from step $h$ to step $h/2$) to estimate the integral and get an approximation to the error.
It is usually done using a canned routine, but the formulas are simple enough that you should be able to implement it yourself.
What is your problem?
