| |
- The basis of Romberg's method is the (composed) trapezoidal rule with equidistant
nodes
T(h)= (h/2)(f(a)+2*∑i=1,..,n-1 f(a+i*h)+f(b)), h=(b-a)/n .
Using the intermediate value theorem it follows that this converges to the
true integral as n -> ∞ for any Riemann integrable integrand f.
If f is two times continuously differentiable on [a,b] then
the quadrature error takes the form
-(1/12)(b-a)*h2*f"(xs) with xs somewhere in [a,b].
- Romberg's method is based on a refined error formula (which can be obtained from the
Euler McLaurin expansion) and which reads
T(h) = I + ∑i=1,..,m pi h2i + O(h2m+2) (**)
.
for f 2m+2 times continuously differentiable on [a,b].
Here I denotes the true integral value and
pi= B2i/(2i!)(f(2i-1)(b)-f(2i-1)(a))
with the Bernoulli numbers B2i . For these holds
2/(2*pi)2i <= |B2i|/(2i!) <= 4/(2*pi)2i .
This looks as if the error of T(h) for a periodic f with period b-a would be zero,
but attention:
(**) is not a power series and the constant in the O(..)-term depends strongly on m.
So this is the case only for special f and sufficiently small h.
Specifically, in f can be expanded in a powerseries at the midpoint of the interval and the
convergence radius is larger than (b-a)/2 and there exists a constant C such that
max x in [a,b] | f(k)(x) | ≤ Ck for all k
and
h < 2 π / C
then (**) can be extended as a infinite series. Consequently, if
f is in addition periodic with period b-a then T(h) will be the exact integral (in exact arithmetic).
- Using (**) and considering T(h) and T(2h), then it is easy to see that
T(h)+(1/3)*(T(h)-T(2h)) = I + O(h4).
Indeed this is the composed Simpson's rule with equidistant nodes.
If one reads T(h) as a polynomial in h2 plus an error term, then one sees
that the unique polynomial of degree <= k, which interpolates
the pairs of numbers (hj2,T(hj)), j=i-k,..,i
and evaluated at the point h=0 delivers an approximate integral value with an error
of the order hi-k2k+2 and exactness order 2k+1.
This is a special case of Richardson's extrapolation and sometimes is called ''extrapolation
to the stepsize zero''.
Interpolation and evaluation can be combined in Neville's algorithm, delivering the simple recursion
Ti,0 = T((b-a)/2i), i=0,1,2,..
Ti,k = Ti,k-1+( Ti,k-1- Ti-1,k-1)/(22k-1), k=1,2,.. i=k,k+1,...
and Ti,i denotes the value of the polynomial mentioned above.
- It is not necessary to begin this recursion with i=0, but any value h=(b-a)/N may serve
as the initial one. Then one computes first the values
T(h), T(h/2), T(h/4), ...
and then the triangular array of approximate integral values Ti,k.
- The values in column k=0,1,... have the order
2k+2 and the total scheme converges columnwise and diagonalwise to the true integral value,
this being true for any Riemann integrable function.
- If f is analytic in an ellipse of the complex plane containing the interval [a,b] in its
interior, then the convergence in the diagonals is even superlinear and the difference to its left neighbor is
an asymptotically correct error estimate for each value Ti,k.
- It is not necessary to use only stepsize halving for using this technique, other stepsize sequences are
possible, but the strong theoretical statements above are proven for the continued halving only.
- If f is not smooth, then this technique looses its advantages, however.
|
| |
- You get back the triangular array of the Ti,k, the so called Romberg tableau.
The first column contains values of the trapezoidal rule and the following columns the extrapolated
values.
If indeed f is 2k+2 times continuously differentiable, then the
quotients
(Ti,k-Ti-1,k)/(Ti-1,k-Ti-2,k)*4k+1,
these are the so called ''control coefficients'', should all be <=1. This can be seen as a
''regularity test'' for f. But this holds for exact computation only. Due to roundoff in the
function evaluation and the summation there might dominate the roundoff errors in the differences
due to cancellation of significant figures, leading to unsensical values in the table of these
control coefficients.
- If a control coefficient cannot be evaluated at all (division by zero) then it is set to zero.
- You also get a table of the quadrature errors if they are known, i.e. in the predefined
testcases and in your own case, if you know the true value.
- finally, you get plots of the integrand and of the quadrature errors for
the different columns. These error plots are in double-log format.
Hence if anything is as in theory (and the roundoff does not begin to
dominate) then you should be able to verify the corresponding orders of
convergence.
The errorplots are suppressed if there are too less data.
We add 1.0e-16 to all errors. Hence if an error occurs as zero, you will not
see the plot since it coincides with the upper boundary of the plot.
- The program tries to compute the true integral with a relative error
tol, but uses 20 halvings of the stepsize at most, i.e. 1048577 function
evaluations at most.
If you see less than columns 0 to 7 then this means that higher
columns have been discarded since the control coefficients couldn't be
evaluated safely due to roundoff accumulation.
|