 The Milne integration rule is given in the table in the
Session 3 notes.
 Consider printing your results for the different rules into
different files. This allows you to use different numbers of points
for each rule (remember that the rules only work optimally for certain
values, like odd numbers for Simpson's rule, 3*n+4 for 3/8'th rule,
and so on).

The integration routines used in Session 3 (trapezoid and simpson), which
were moved into integ_routines.cpp, are not completely
general.
In particular, they assume that the lower limit of the integral
("min") is zero. You can see this by noting how x is calculated
in the "for loop" in each routine. So you'll get the wrong answer
for any integral that has a nonzero lower limit. The fix is easy:
just replace "x = interval*(n1)" by "x = min + interval*(n1)" and
so on. (Note that the Gaussian quadrature routine is ok.)

If you're having some difficulty seeing how to code the threeeighths
or Milne rules, consider this rewrite of the heart of the
simpson routine:
double x1, x2;
for (n=2; n<num_pts; n+=2) // step through the integral, 2 points at a time
{
x1 = min + interval * (n1);
x2 = x1 + interval;
sum += 4. * f(x1) + 2. * f(x2);
}
sum += f(min)  f(max); // note the "" sign
sum *= interval/3.
So now there is one routine rather than two, which covers each
"elementary interval" and then fixes the endpoints. Note
the minus sign in the latter; what does it do?

To apply the "Empirical Error Analysis" as described in
the Session 4 background notes to find
the approximation error, it is
not necessary to compare the result for N to 2N. You can compare
N to (2N+3) or (4N+1) (for the Milne integration rule, for example) or any
higher number of points for which the approximation
error is significantly smaller
than the error for N points. The point is that Equation (4.14) in
the notes holds
because the exact answers cancel and the error is dominated by the
error for N points. For the roundoff error region, the opposite is
true: the result with the greater number of points should have the
dominant error.

If you are having trouble getting a gsl integration routine to work,
take a look at the
example
for the routine gsl_integration_qags. You need to define a
gsl_integration_workspace and a gsl_function as in the example
(if you use the qng routine instead of qag or qags, you
won't need the workspace definition). If you are unfamiliar
with the use of C++ structures, please ask about one of the instructors
what is going on!
 When you define your function for a gsl routine, you need to
specify that it has two arguments: x and the "params"
pointer. In the example mentioned above, this is the function:
double f (double x, void * params)
{
double alpha = *(double *) params;
double f = log(alpha*x) / sqrt(x);
return f;
}
The purpose of "params" is to pass any arbitrary number of parameters
of arbitrary type to the function. So in the example, the integrand
depends on a parameter "alpha", which is obtained from params.
If you don't need to pass any parameters, you can define your function
like this:
double f (double x, void * )
{
double f = log(2.*x) / sqrt(x);
return f;
}
and it will work fine. (Note: The standard list of compiler
warnings will complain that "params" is not being used if you include
the word "params"
in the argument list of f. You tell the compiler it doesn't
appear in this case simply by leaving it out of the list, but keeping
the "void *" part, as shown in this example.)
 If you don't pass any parameters to your integrand function
(that is, any parameters are hardcoded into the function itself),
you might decide to use our standard routines that have only
an argument x:
double my_integrand (double x)
{
// stuff here
return my_integrand;
}
which causes problems when you try to call GSL routines in the
same program. That is because GSL insists that the function
have two arguments, a double and a void pointer, even if you don't
use the latter. You could write two independent functions with the
same integrand, but this could lead to errors; you really only want
to enter the integrand in one place. Here is a simple solution:
define my_integrand as above, with the actual integrand
function. Then define a function to pass to GSL:
double my_gsl_integrand (double x, void *)
{
return (my_integrand(x));
}
which simply evaluates the other function. Then pass this
function to GSL and everyone is happy!