- 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).
If you're having some difficulty seeing how to code the three-eighths
or Milne rules, consider this rewrite of the heart of the
double x1, x2;
for (n=2; n<num_pts; n+=2) // step through the integral, 2 points at a time
x1 = min + interval * (n-1);
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?
If you are having trouble getting a gsl integration routine to work,
take a look at the
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);
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
double f (double x, void * )
double f = log(2.*x) / sqrt(x);
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 hard-coded 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
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 *)
which simply evaluates the other function. Then pass this
function to GSL and everyone is happy!