# Physics 6810: Assignment #2

Here are some hints, suggestions, and comments on the assignment.
1. Follow-up/Completion of Session 3
• 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 non-zero lower limit. The fix is easy: just replace "x = interval*(n-1)" by "x = min + interval*(n-1)" and so on. (Note that the Gaussian quadrature routine is ok.)
• If you're having some difficulty seeing how to code the three-eighths 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 * (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?
• 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 round-off 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 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
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!