Computer Science/Mathematics 455 Lecture Notes
Lecture # 28
Numerical Integration – Romberg Integration
Need to compute
∫b
a
I =
f(x)dx
where f(x) has no known antiderivative.
Let xk = a + kh,k = 0,…,n where h = (b − a)/n. The composite
trapezoid rule is
n−1
T0(h) = h[f(x0) + f(xn)] + h ∑ f(xk)
2 k=0
To halve the interval size, we need to evaluate the function f at the
midpoints
We obtain the formula
xk+1/2 = [xk + xk+1]/2.
n−1 n−1
T0(h/2) = h[f(x0) + f(xn)] + h ∑ f(xk) + h ∑ f(xk+1/2)
4
= 1T0(h)+h∑f(xk+1/2)
n−1 2 2 k=0
2 k=0 2 k=0
Thus with the same function evaluations for T0(h/2), we can compute
both T0(h) and T0(h/2). This can be very useful. In fact, for the same func-
tion evaluations that we need to compute T0(h/2k) we can obtain T0(h), T0(h/2), . . . , T0(h/2k). That information may be used to obtain a more accurate approximation.
If I = ∫ b f (x)dx then assuming a sufficient number of continuous deriva- a
tives (2n + 2) we have
T0(h) = I + a2h2 + a4h4 + ··· + a2nh2n + O(h2n+2).
Likewise, for h/2 we have
T0(h/2) = I + a2h2/4 + a4h4/16 + · · · + a2nh2n/4n + O(h2n+2).
1
If we let then
where
T1(h) = (4T0(h/2) − T0(h))/3
T1(h) = I + a′4h4 + a′6h6 + ··· + a′2nh2n + O(h2n+2)
a′4 = −a4/4, a′2k = −(1 − 1/4k−1)a2k/3.
Thus all of these coefficients become smaller in magnitude.
T1(h) is the composite Simpson’s rule. There is a formula for it, but this
one is easier to use. We do not need to stop here.
If we have T0(h),T0(h/2), and T0(h/4) we can get T1(h) and T1(h/2).
Since
T1(h/2) = I + a′4h4/16 + a′6h6/64 + · · · + a′2nh2n/4n + O(h2n+2).
We can do the same trick again and eliminate the h4 term. Thus obtaining T (h) = (16T (h/2) − T (h))/15 = I + a(2)h6 + O(h8).
2116
Now we have beginnings of an algorithm to obtain higher order approxi- mations.
For instance, suppose we have Tk−1(h) from T0(h), T0(h/2), . . . , T0(h/2k−1) that satisfies
Tk−1(h) = I + a(k−1)h2k + a(k−1)h2k+2 + O(h2k+4). 2k 2k+2
So that once we generate T0(h/2k) followed by T1(h/2k−1), · · · Tk−1(h/2), the last of which satisfies
Tk−1(h/2) = I + a(k−1)h2k/4k + a(k−1)h2k+2/4k+1 + O(h2k+4). 2k 2k+2
Therefore we can compute
Tk(h) = (4kTk−1(h/2) − Tk−1(h))/(4k − 1) = I + a(k) h2k+2 + O(h2k+4). (1)
This technique of using several low order approximations to obtain a high order approximation is called Richardson extrapolation. When applied to the Trapezoid rule, it is called Romberg integration.
The basic idea is to set up a table as follows. 2
2k+2
T0(h) T0 (h/2) T0 (h/4) T0 (h/8) T0 (h/16)
T1(h) T1 (h/2) T1 (h/4) T1 (h/8)
T2(h) T3(h) T4(h) T2 (h/2) T3 (h/2)
T2 (h/4)
The rows (going to the right), the columns (going down), and the diag- onals (up and to the right), all converge to the integral. However, the rows converge the most quickly.
In a code, you do not need to keep this entire table around. Instead you can get away with T0(h/2k), T1(h/2k−1), . . . , Tk(h) as needed.
The method comes with its own error estimate. It is easily verified that
Tk−1(h) − Tk(h) = a(k−1)h2k + (a(k−1) − a(k) )h2k+2 + O(h2k+4) 2k 2k+2 2k+2
= Tk−1(h) − I − a(k) h2k+2 + O(h2k+4) 2k+2
= Tk−1(h) − I + O(h2k+4)
Thus Tk−1(h) − Tk(h) captures the first order term of the error in Tk−1(h)
and a good stopping criterion is
|Tk−1(h) − Tk(h)| ≤ ε
for some tolerence ε.
Incidentally, both the trapezoid method and Simpson’s rule also have er-
ror estimates that may be gleaned from extrapolation. For Simpson’s rule, suppose we have computed both T1(h) and T1(2h) (the second can be ob- tained at no additional cost).
T1(2h) − T1(h) = −15a′4h4 + O(h4) = 15(I − T1(h)) + O(h6). So a good error estimate is
Err = (T1(2h) − T1(h))/15. If |Err| < epsilon we generally accept
T2(2h) = T1(h) + Err
Verify for yourself that a good error estimate for the trapezoid rule is
as the integral.
Err = (T0(2h) − T0(h))/3 = I − T0(h) + O(h4).
3