NUMERICAL INTEGRATION, Continued SPH 5421 notes.031
The trapezoid rule and Simpson's rule are both special cases of another
algorithm for performing numerical integration: Romberg Integration. For a
very sketchy description, see section 4.3 of Numerical Recipes. The basic idea
with Romberg integration is to fit a polynomial to the function values at
equally spaced points in the interval of integration, and then compute what the
integral of that polynomial would be. The result is a highly efficient method
for performing numerical integration. The program below implements Romberg
integration. It is based on the subroutine QROMB in Numerical Recipes.
This program will give (almost) exact answers for integrals of polynomials up to
to degree 2 * K - 1, where K is one of the parameters that you can change in the
early steps of the program. In the example below, K = 5.
Romberg integration is remarkably more efficient than either the trapezoid
rule or Simpson's rule. For the function shown in the program below, with
convergence criterion eps = 1e-6, the three algorithms all gave almost exactly
the same answers, but the numbers of function evaluations were:
Trapezoid rule : 4097 evaluations
Simpson's rule : 129 evaluations
Romberg : 17 evaluations
Clearly for practical purposes, you should use Romberg integration in
preference to the other two. The other two programs are essentially obsolete.
Unfortunately, the efficiency of the three programs is inversely related to their
ease of understanding: the trapezoid rule program is simple and straightforward,
the Simpson's rule program is just slightly more complicated, and the Romberg program
is nearly indecipherable. There is really no difference in ease of use:
in all three, the basic problem is to provide a module which defines the function
to be integrated.
Another collection of efficient integration algorithms is based on
Gauss quadrature. The main difference with Gauss quadrature is that the
algorithms are based on unequal spacing of the subdivisions of the interval of
estimation, with the exact spacings determined by the roots of sets of
orthogonal polynomials (Legendre polynomials). Gauss quadrature is described
in section 4.6 of Numerical Recipes, and will not be discussed further here.
==================================================================================
options linesize = 80 mprint ;
footnote "~john-c/5421/romberg.sas &sysdate &systime" ;
/* Proc iml program to perform numerical integration ... */
/* Uses Romberg algorithm (Numerical Recipes, p. 114-115). */
proc iml ;
/*********************************************************************/
start func(x, fvalue, evals) ;
fvalue = x**4 * log(x + sqrt(x*x + 1)) ;
evals = evals + 1 ;
finish func ;
/*********************************************************************/
/* Begin Romberg main program ... */
a = 0 ; b = 2 ;
eps = 1e-6 ;
jmax = 20 ;
jmaxp = jmax + 1 ;
nmax = 10 ;
s = j(20, 1, 0) ;
h = j(20, 1, 0) ;
c = j(20, 1, 0) ;
d = j(20, 1, 0) ;
k = 5 ;
km = k - 1 ;
evals = 0 ;
ss = 0 ;
normend = 0 ;
h[1] = 1 ;
do j = 1 to jmax ; /* Begin do loop A */
/*********************************************************************/
/* Start of trapezoid section.... */
/*********************************************************************/
if j = 1 then do ; /* Begin do loop B */
run func(a, fa, evals) ;
run func(b, fb, evals) ;
s[j] = 0.5 * (b - a) * (fa + fb) ;
sj = s[j] ;
it = 1 ;
goto jump1 ;
end ; /* End do loop B */
tnm = it ;
del = (b - a) / tnm ;
x = a + .5 * del ;
sum = 0 ;
do iit = 1 to it ; /* Begin do loop C */
run func(x, fx, evals) ;
sum = sum + fx ;
x = x + del ;
end ; /* End do loop C */
sj = s[j] ;
sjm1 = s[j - 1] ;
s[j] = (sjm1 + (b - a) * sum / tnm) / 2 ;
it = 2 * it ;
sj = s[j] ;
/*********************************************************************/
/* End of trapezoid section.... */
/*********************************************************************/
jump1:
if j >= k then do ; /* Begin if loop D */
/*********************************************************************/
/* POLINT (polynomial interpolation) section ... */
/*********************************************************************/
offset = j - k ;
x = 0 ;
ns = 1 ;
diff = abs(x - h[1 + offset]) ;
do jk = 1 to k ; /* Begin do loop E */
difftemp = abs(x - h[jk + offset]) ;
if (difftemp < diff) then do ; /* Begin if loop F */
ns = jk ;
diff = difftemp ;
end ; /* End if loop F */
c[jk] = s[jk + offset] ;
d[jk] = s[jk + offset] ;
end ; /* End do loop E */
y = s[ns + offset] ;
ns = ns - 1 ;
do m = 1 to k - 1 ; /* Begin do loop G */
do inm = 1 to k - m ; /* Begin do loop H */
ho = h[inm + offset] - x ;
hp = h[inm + m + offset] - x ;
w = c[inm + 1] - d[inm] ;
den = ho - hp ;
if (den = 0) then stop ;
den = w / den ;
d[inm] = hp * den ;
c[inm] = ho * den ;
dinm = d[inm] ; cinm = c[inm] ;
end ; /* End do loop H */
cns1 = c[ns + 1] ; dns = d[ns] ;
ns2 = 2 * ns ;
kmm = k - m ;
cns1 = c[ns + 1] ;
dns = d[ns] ;
if (ns2 < kmm) then dy = c[ns + 1] ;
if (ns2 >= kmm) then do ; /* Begin do loop I */
dy = d[ns] ;
ns = ns - 1 ;
end ; /* End do loop I */
y = y + dy ;
ss = y ;
dss = dy ;
end ; /* End do loop G */
/*********************************************************************/
/* End of POLINT section ... */
/*********************************************************************/
abdss = abs(dss) ;
abss = abs(ss) ;
eabss = eps * abss ;
if ((abdss < eabss) | ss = 0) then do ; /* Begin if loop J */
print " Integral, error estimate, function evals: " ss dss evals ;
normend = 1 ;
stop ;
end ; /* End if loop J */
end ; /* End if loop D */
s[j + 1] = s[j] ;
h[j + 1] = .25 * h[j] ;
end ; /* End do loop A */
if normend = 0 then print " ... Too many steps ..." j jmax evals ss dss ;
quit ;
==================================================================================
The SAS System 1
10:05 Tuesday, April 18, 2000
SS DSS EVALS
Integral, error estimate, function evals: 8.1533644 1.0718E-7 17
==================================================================================
The FORTRAN subroutine package IMSL has a version of the Romberg integration
algorithm which is somewhat shorter than the code above, and in fact is amazing
for its efficiency in implementing the algorithm. The program below,
romberg3.sas, is a SAS PROC IML version of this program. It has most of the
capabilities of the program above and runs somewhat faster. It does not print
an error estimate (DSS in the above program), and the parameter K does not
appear explicitly in the program. See the IMSL manual for a description and
explanation of the underlying algorithm.
==================================================================================
options linesize = 80 mprint ;
footnote "~john-c/5421/romberg3.sas &sysdate &systime" ;
/* Proc iml program to perform numerical integration ... */
/* Uses Romberg algorithm ... from IMSL. */
proc iml ;
/*********************************************************************/
start func(x, fvalue, evals) ;
fvalue = cos(x) ;
evals = evals + 1 ;
finish func ;
/*********************************************************************/
/* Begin Romberg main program ... */
xl = 0 ; xu = 20 ;
h = xu - xl ;
ndim = 20 ;
eps = 1e-6 ;
aux = j(200, 1, 0) ;
evals = 0 ;
delt2 = 0 ;
p = 1 ;
run func(xl, fxl, evals) ;
run func(xu, fxu, evals) ;
aux[1] = .5 * (fxl + fxu) ;
e = eps / abs(h) ;
hh = h ;
jj = 1 ;
do i = 2 to ndim ;
y = aux[1] ;
delti = delt2 ;
h0 = hh ;
hh = hh / 2 ;
p = p / 2 ;
x = xl + hh ;
sm = 0 ;
do j = 1 to jj ;
run func(x, fx, evals) ;
sm = sm + fx ;
x = x + h0 ;
end ;
im1 = i - 1 ;
aux[i] = .5 * aux[i - 1] + p * sm ;
q = 1 ;
ji = i - 1 ;
do j = 1 to ji ;
ii = i - j ;
q = 4 * q ;
aux[ii] = aux[ii + 1] + (aux[ii + 1] - aux[ii]) / (q - 1) ;
end ;
delt2 = abs(y - aux[1]) ;
if ((i - 5) < 0) then goto jump1;
if ((delt2 - e) <= 0) then do ;
y = h * aux[1] ;
print "Algorithm converged ... " ;
print "Integral estimate = " y " Number of fn evals = " evals ;
abort ;
end ;
if ((delt2 - delti) < 0) then goto jump1;
if ((delt2 - delti) >= 0) then do ;
print "Some kind of error occurred ... ";
print "Integral estimate = " y " Number of fn evals = " evals ;
abort ;
end ;
jump1:
jj = jj + jj ;
end ;
y = h * aux[i] ;
print "No convergence ... " ;
print "Integral estimate = " y " Number of fn evals = " evals ;
abort ;
quit ;
==================================================================================
PROBLEM 33
Use Romberg integration to find the mean and variance of a distribution
with the following PDF:
f(x) = 4/3 if 0 <= x <= .5,
f(x) = -(8/3) * (x - 1) if .5 < x <= 1.
==================================================================================
/home/walleye/john-c/5421/notes.031 Last update: May 2, 2000