780.20: 1094 Session 3

Handouts: Background notes for Session 3; Landau/Paez chapter 4 excerpts; Gnuplot fitting example; Integrals with singularities; C++ formatting

Your goals for today and next time (in order of priority):

Please work in pairs. The instructors will bounce around 1094 to ask and answer questions.

You should make or change to a 780 sub-directory, download session03.tarz from the 780 homepage and unpack it ("tar xfvz session03.tarz").

Leftover Tasks from Session 2

Here are the priority tasks to finish from Session 2. Don't worry about other items; when you complete these, move on to Session 3.

  1. Summing in Different Orders. Finish the entire section. Be careful in part 5: There is one bug that cannot be caught by the compiler (check against the answer in the comment to track it down).
  2. Bessel 1-3 Skip Bessel 4.
  3. Mathematica Scavenger Hunt. Do the first one only.

Comparing Floating Point Numbers

A common task in a computational problem is to check whether two floating point numbers are the same. In C++, the comparison operator is ==, so it might seem that a simple if statement would do the job. Here's an example of why this fails in general:

  1. Look at the file number_comparison.cpp in an editor. What does it do?

  2. Use make_number_comparison to compile and link it. Run the program. Why doesn't it give the answer you want? Add statements to print out x1 and x2 to check your response.

  3. Can you think of a better way to compare two numbers?

Makefiles for multiple project files (including header file)

Many of the example programs started as "all-in-one" C programs from the Landau/Paez text. One of these was integ.c, which we converted to C++ and then split up into:

There is also a function in gauss.cpp and there is make_integ_test to compile it all. In a subsequent step, we modified the codes so that the function is passed as an argument to the integration functions.

The idea is that the integration routines should be isolated in a file by themselves. The main program just invokes these routines. The header file conveys the prototype information about the integration routines to the main program and to any other functions that might call the routines. (Note: Later we'll consider going further and defining an integration class.)

  1. Look at the original file (integ.c) and see how it was converted and split up into integ_test.cpp, integ_routines.cpp, and integ_routines.h (the cpp files are printed out for you). Also look in the ORIGINALS subdirectory to see the change to passing the function for the integrand. Why do these changes improve the code?

  2. Change integ_test.cpp to output relative (rather than absolute) errors.
  3. Create the executable integ_test using the makefile make_integ_test and run it to generate the integ.dat output file. Note that only the files that have changed are recompiled by the makefile!
  4. Use Gnuplot to reproduce Figure 4.3 on page 59 of the Landau-Paez handout. Explain what you can learn from the plot to your partner (and then an instructor). Consider all regions of the graphs.

  5. Change the loop in integ_test.cpp so that the points on the log-log plot are evenly spaced. What did you do?

Finding the Approximation Error From a Log-Log Plot

From the plot in the previous section, we can estimate the approximation errors by eye. Now we want to actually fit lines to find the slope using Gnuplot. Use the handout as a guide.

  1. Modify the code so that it outputs the logarithm base 10 (log10) of N and the relative errors.
  2. Find the slopes of the trapezoid and Simpson's rule plots in the regions where they are linear.
  3. Are the results consistent with the analysis in the text? Can you fit the round-off error region?

Coding an Algorithm

This is primarily practice converting an algorithm to working code. Your goal is to add a new function to integ_routines.cpp, called three_eighth, which implements the 3/8 rule from Table 4.1.

  1. From the rule and the discussion in the text, write some pseudocode to explain how you would integrate a function with this rule.
  2. Implement your pseudocode in C++ by adding a new routine (called three_eighth) to integ_routines.cpp (be sure to add a prototype to integ_routines.h).
  3. Modify integ_test.cpp to output to a new file the results from the new routine, and add them to the previous plot. Explain the approximation error. [Warning: If you don't change the way that integ_test.cpp loops through the number of intervals, you will likely run into a subtle error that prevents you from getting the approximation predicted theoretically!]

GSL Scientific Library Yet Again

  1. Go to the web page with the GSL manual (it is linked from the 780.20 web page). Find an appropriate integration routine for the test integral we've been working on.
  2. Add another calculation to integ_test.cpp (with output to a file) using this GSL routine.
  3. Compare the accuracy of the GSL routine to the others on an error plot.

More Complicated Integrals

  1. Choose one of the integrals with singularities from the handout [Eq.(9) is a good choice!].
  2. Determine an "exact" answer analytically or using Mathematica (see instructors for assistance).
  3. Modify your code so that you can analyze the integral with one of the integration rules several ways, using an error plot to compare (these may not all apply to a given integral):
    1. Brute force (unmodified)
    2. Subtracting the singularity
    3. Changing variables
    What do you find?

780.20: 1094 Session 3. Last modified: 08:24 am, January 11, 2006.