6810: 1094 Session 3

Online handouts: Gnuplot fitting example; C++ formatting; listings of codes

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

• Finish some leftover tasks from Session 2
• Look at a quick demo about comparing floating-point numbers
• Use a makefile to compile a project with multiple .cpp files and a header (.h) file
• Duplicate the figure in the notes 3d with a log-log plot, then fit the slopes with Gnuplot
• Practice coding an algorithm: Add a new subroutine with the 3/8 rule and analyze [PS#2]
• Add an appropriate integration routine from the GSL and analyze [PS#2] the error

Please discuss and compare answers with others. The instructors will bounce around 1094 to ask and answer questions.

You should make or change to a 6810 sub-directory, download session03.zip from the 6810 homepage and unzip it.

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. Make sure to finish the entire section. Be careful in part 5: There is one bug (a typo) that cannot be caught by the compiler.
2. Bessel 1 (you'll do Bessel 2-3 in PS#1b).

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 == (not just one =), 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. Compile, link, and run the program (there is a makefile). Why doesn't it give the answer you want? Add statements to print out x1 and x2 to check your response. What did you add?
   
   
   
   
3. Suggest a better way to compare two numbers based on the idea that the relative inaccuracy of any number can be as large as a specified precision eps (which may be greater than the machine precision).
   
   
   
   
   

Numerical Derivatives: Pass 1

The Session 3 notes have a short introduction to numerical differentiation (see the Hjorth-Jensen notes online for more detail). These are among the simplest algorithms for us to derive and to verify the theoretical approximation and round-off errors (the errors for most real-world algorithms are noisy).

1. Look at the file derivative_test_simple.cpp while reviewing the discussion in the notes. (Also look at the last section in the notes, which describes pointers to functions, which are used here.) What part(s) of the code do you not understand? What is being printed?
   
   
   
   
2. Use the makefile to compile, link and run the code, generating the file derivative_test_simple.dat. There is a missing #include statement you need to add (compare with previous codes). Add a statement to print a header line to the output file with a label for each column (and start the line with # so that it doesn't screw up gnuplot). Show an instructor your output file.
   
   
3. Make a graph with gnuplot with two plots: the logarithm of the relative error for forward difference vs. the logarithm of h (which is Delta-x) and the analogous plot for central difference. Print the graph and attach.
   
   
4. Are the slopes in each region consistent with the analysis of errors in the notes? Which is the better algorithm? Explain. What value of h is optimal for each algorithm? If you switched to single precision, would the slopes of the lines change? What would the graph look like?
   
   
   
   
   
   
   

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:

• integ_test.cpp, which has the main program and the function to be integrated
• integ_routines.cpp, which has the integration functions themselves
• integ_routines.h, which has function prototypes (a prototype tells the compiler the function return type and the type of all its arguments)

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. Compare the trapezoid_rule and Simpsons_rule functions in integ_routines.cpp to equations (3.15)-(3.19) and the table in the Session 3 class notes. Can you see how the algorithms are implemented?
   What is the advantage of having the routines in a separate file from the main program if you want to test that they work on known integrals or if you later improve the algorithm?
   
   
   
   
2. Create the executable integ_test using the makefile make_integ_test and run it to generate the integ.dat output file. Does the output file makes sense?
   
   
   
   
3. Change integ_test.cpp to output relative (rather than absolute) errors. Note that only the files that have changed since the last compilation are recompiled each time!
4. Use Gnuplot to reproduce the figure in Section 3d of the notes. Briefly explain what you can learn from the plot (remember slopes!). Consider all regions of the graphs.
   
   
   
   
   
   
5. Bonus: Change the loop in integ_test.cpp so that the points on the log-log plot are evenly spaced. What did you change?
   
   
   
   
   
   

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. What are the slopes of the trapezoid and Simpson's rule plots in the regions where they are linear?
   
   
   
   
3. Are the slopes consistent with the analysis in the text? Now try to fit the round-off error region and interpret the "slope".
   
   
   
   
   

Extra: Coding an Algorithm (complete in PS#2)

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_rule, which implements the 3/8 rule from the table in Section 3d of the notes.

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_rule) 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!]
   
   
   
   

Extra: GSL Scientific Library Yet Again (complete in PS#2)

1. Go to the web page with the GSL manual (it is linked from the 6810 web page). Find an appropriate integration routine for the test integral we've been working on. Which one did you choose, and why?
   
   
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.