780.20: 1094 Session 2

Handouts: Landau-Paez chapter 3 excerpt; gnuplot mini-manual; gnuplot example; printouts of several codes.

Your goals for this session (in order of priority; you won't finish them all!):

Please work in pairs and briefly answer questions asked here on this sheet. The instructors will bounce around 1094 and answer questions (don't be shy to ask about anything!).

You should go to a directory on the shared disk (see Session 1 for creating one), download session02.zip from the 780 homepage and unzip it in your directory.

Leftover Tasks from Session 1

You should finish or correct the following subtasks from Session 1, which we will build on. If you are correcting something, check it with an instructor sometime during the class today.

  1. Dev-C++: Make sure that you've correctly modified area.cpp to calculate the volume of a sphere and verified that it works.
  2. Overflows, Underflows, and Machine Precision: Complete part 2, modifying "precision.cpp" so that the result for double precision is correct (10e-12 is incorrect).
  3. Setting g++ Compile Options: Make sure to do this.
  4. Using the GSL Scientific Library: Complete parts 1, 2, and 3.

Round-Off Errors with Quadratic Equations

This section assumes you have read the background notes for Session 2, which describe round-off errors and the quadratic equation.

  1. Review the discussion about round-off errors and the quadratic equation in Session 2 notes with your partner. Make sure you go over the pseudo-code for the test case and understand it.
  2. Look at the code quadratic_equation_1a.cpp in Dev-C++ and/or on the handout. It represents the pseudo-code as it has just been typed in, mistakes and all. Note that it is not indented consistently. Under Tools->Editor Options, change the tab size at the bottom right to 2. Highlight everything between the main () { line and the final } (not including these) and then hit Tab to indent. The convention we'll use is 2 spaces for each sub-level.
  3. Now try compiling it in Dev-C++. You should get a bunch of warnings and errors. Each one comes with the line number of the error. Take each one in turn and mark on the handout what you think the error is. Each warning is a clue to a mistake in the code (such as a typo)! [One hint: "setprecision" is defined in the header file "iomanip".] If you get stuck, quadratic_equation_1.cpp is the corrected version (without the "a" after the "1"). List the errors here:

  4. Correct these errors and see that it compiles and runs. Try a=1, b=2, c=0.1 and see that the roots differ. Which two are the most accurate? (Hint: check the Session 2 notes.)

  5. Now look at quadratic_equation_2.cpp, which is the second pass at the code. Check the comments for a description of changes. (This code is printed on the same handout as 1 and 1a.) Mark every C++ line on the printout that you don't understand completely (it may be most of them!).
  6. Compile and run quadratic_equation_2.cpp, and enter a=1, b=2, c=.1 again. The output is now more detailed and an output file named "quadratic_eq.dat" is created. We'll discuss outputting to a file more later; for now, this is an example you can follow to do it on your own (there were also examples in Session 1).
  7. Follow the "Plotting Data from a File with Gnuplot" handout to duplicate the plot on the back. Try some of the extra commands, including output to a postscript file. Print out the file on Smith1011s and hand it in at the end.
  8. Is the plot qualitatively and quantitatively consistent with the analysis in the notes? Explain in a few sentences.

Summing in Different Orders (cf. Landau-Paez 3.4-3 Handout)

Here we consider another example of subtractive cancellation and how computer math is not the same as formal math. Imagine we need to sum the numbers 1/n from n=1 to n=N, where N is a large number. We're interested in whether it makes any difference whether we sum this as:
   1 + 1/2 + 1/3 + ... + 1/N (summing up)
or as:
   1/N + 1/(N-1) + ... + 1/2 + 1 (summing down)
In formal mathematics, of course, the answer are identical, but not on a computer!

To help think about the sum up vs. sum down problem, it's useful to think first about a simpler version. Suppose you want to numerically add 107 small numbers a = 5x10-8 to 1 (so the exact answer is 1.5). For the first three questions, predict the answer before running the code.

  1. If you add 1 + a + a + ... in single precision, what do you expect to get for the total? (Hint: recall the discussion of "machine precision.")

  2. If, instead, you add a + a + ... + 1 in single precision, do you expect a different answer? Which will be closer to the exact answer?

  3. If you repeat the exercise in double precision, what do you expect will happen?

  4. Write a brief "pseudocode" here for a program that would compare these two ways of summing. Then look at the file demo1a.cpp in Dev-C++ (also on a handout); does it look like an implementation of your pseudocode?

  5. Try to compile demo1a.cpp. Identify and fix the errors, until it compiles and runs correctly (verify the results against the comments in the file). Correct your answers to 1,2,3 if necessary (and understand them!). How should you indent it (see demo1.cpp)?
  6. Challenge: Predict the results for 108 additions (instead of 107):

    Now change the code, and see if your prediction is correct (or that you can explain it after the fact).

Bessel 1: Making Another Plot with Gnuplot

The next three tasks (Bessel 1-3) build on the discussion in sections 3.6 to 3.9 of the Landau-Paez book. These sections are copied on one of our handouts. Skim these sections before going further.

  1. Compile and run bessel.cpp. You'll generate a data file "bessel.dat". By looking at the code (and chapter 3), figure out what the columns in the file mean. What Bessel function is being output? What are the columns?

  2. Make a plot of the third column as a function of the first column using gnuplot. (Use the example in the handout as a guide. You might also find the "GNUPLOT Manual and Tutorial" useful.) Try plotting BOTH the second and third columns vs. the first column on the same plot (i.e., two curves). Print out your plot.

Bessel 2: Landau-Paez Section 3.9 Assessment (page 40 on handout)

  1. Modify the output statements in bessel.cpp to change the number of digits in the output (look back at previous codes for clues).
  2. Modify the code bessel.cpp (give it another name) to do part 2 of 3.9.
  3. To answer part 3 of 3.9, make an appropriate plot of the error vs. x and interpret the graph.

Bessel 3: Using the GSL Scientific Library (continued)

Make sure you have completed the GSL task from Session 1 before doing this.

  1. Modify bessel.cpp to also calculate the spherical Bessel function using the GSL routine gsl_sf_bessel_jl for comparison. The GSL documentation gives the usage of this function as:
        Function: double gsl_sf_bessel_jl (int l, double x)
    "This routine computes the regular spherical Bessel function of order l, j_l(x), for l >= 0
    and x >= 0."
    How does the accuracy of the downward recursion routine in bessel.cpp compare to that of the GSL function? (How can you tell?)

Are Round-Off Errors Random? (Advanced)

In the book "Computational Physics", Landau and Paez say that the round-off errors in single-precision are distributed approximately randomly. Your task is to test whether this is true and what the distribution actually looks like. More specifically, if we generate a large set of zc = z (1 + epsilon) numbers from some calculation (e.g., by taking the square root of a bunch of numbers), how are the different epsilons distributed?

  1. Devise a scheme to test for the distribution of round-off errors.
  2. Carry out your scheme with a C++ program.
  3. Make appropriate plots to explore your results.

780.20: 1094 Session 2. Last modified: .