Computational Physics - Department of Physics

116 5 Numerical Integration

voidadaptive_integration(doublea,doubleb,doubleIntegral,intn,intsteps,double
(func)(double))
if( steps > maxrecursions){
cout <<'Too many recursive steps, the function varies too much'<< endl;
break;
}
doublec = (a+b)*0.5;
// the whole integral
doubleI0 = trapezoidal_rule(a, b,n, func);
// the left half
doubleI1L = trapezoidal_rule(a, c,n, func);
// the right half
doubleI1R = trapezoidal_rule(c, b,n, func);
if(fabs(I1L+I1R-I0) < tolerance ) integral = I0;
else
{
adaptive_integration(a, c, integral,intn, ++steps, func)
adaptive_integration(c, b, integral,intn, ++steps, func)
}
}
// end function adaptive_integration

The variablesintegralandstepsshould be initialized to zero by the function that calls the
adaptive procedure.

5.3 Gaussian Quadrature

The methods we have presented hitherto are taylored to problems where the mesh pointsxi
are equidistantly spaced,xidiffering fromxi+ 1 by the steph. These methods are well suited
to cases where the integrand may vary strongly over a certainregion or if we integrate over
the solution of a differential equation.
If however our integrand varies only slowly over a large interval, then the methods we
have discussed may only slowly converge towards a chosen precision^1. As an example,

I=

∫b

1

x−^2 f(x)dx,

may converge very slowly to a given precision ifbis large and/orf(x)varies slowly as function
ofxat large values. One can obviously rewrite such an integral by changing variables tot= 1 /x
resulting in

I=

∫ 1

b−^1

f(t−^1 )dt,

which has a small integration range and hopefully the numberof mesh points needed is not
that large.
However, there are cases where no trick may help and where thetime expenditure in
evaluating an integral is of importance. For such cases we would like to recommend methods
based on Gaussian quadrature. Here one can catch at least twobirds with a stone, namely,
increased precision and fewer integration points. But it isimportant that the integrand varies
smoothly over the interval, else we have to revert to splitting the interval into many small
subintervals and the gain achieved may be lost.
The basic idea behind all integration methods is to approximate the integral

(^1) You could e.g., impose that the integral should not change asfunction of increasing mesh points beyond the
sixth digit.