Physics 6810: Assignment #2

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).
   • 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?
   • 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!