Physics 780.20: Assignment #2

1. Follow-up/Completion of Session 3
   • 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.)