*Handouts:* Background notes for Session 2; Landau-Paez chapter 3
excerpt;
gnuplot mini-manual; gnuplot example; printouts of several codes.

*
Your goals for this session (in order of priority):
*

- Finish (some) leftover tasks from Session 1 (if any)
- Look at a round-off error demonstration with quadratic equations
- Discuss a special-case problem to help understand Landau-Paez 3.4 problem 3
- Debug and format (using "indent") sample codes
- Make some plots with gnuplot
- Assess different methods of calculating spherical Bessel functions (section 3.9 of Landau-Paez)
- Use the spherical Bessel function from the Gnu Scientific Library (GSL) in a code
- Are round-off errors random? (Advanced)
- Mathematica Scavenger Hunt

Please work in pairs. Dick and Daniel will bounce around 2082 and answer questions (don't be shy to ask about anything!).

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

It is not essential that you finish absolutely everything from Session 1, but you need to complete the following subtasks, which we will build on. Check them with Dick sometime during the class today (even if you did them last time).

**Linux Environment:**Successfully run "area", answer part 7, and modify area.cpp to calculate the volume of a sphere (part 8).**Overflows, Underflows, and Machine Precision:**Complete part 2, running "precision.cpp".**Using a Makefile:**Successfully complete part 3 (noting that "floats" should be "flows"!).**Using the GSL Scientific Library:**Complete parts 1 and 2.

[Note: Last year this was presented in "lecture" mode. As an experiment, we'll have you simply read about it and run the demonstration codes.]

- Read the discussion about round-off errors and the quadratic equation in the background notes for Session 2 (section 2b). Make sure you go over the pseudo-code for the test case and understand it.
- Look at the code quadratic_equation_1a.cpp in an editor 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. We can fix this with the "indent" tool. Just type "indent quadratic_equation_1a.cpp" at the command line and then look at the code again. It will have been indented automagically!
- Now try compiling with "make -f make_quadratic_equation_1a". You should get a bunch of warnings and errors. Each one comes with a number, which is the line number of the error. Take each one in turn and mark on the handout what you think the error is (you'll need to identify the line numbers in your editor; ask an instructor if you don't know how). 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").
- Correct these errors and see that it compiles and runs.
Try a=1, b=2, c=0.1 and see how the roots differ. Which are the best
ones (check the notes)?

- Now compile and run 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.) 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.
- Follow the "Plotting Data from a File with Gnuplot" handout to duplicate the plot on the back. (You don't need to print the file.)
- Is the plot consistent with the analysis in the notes? Explain.

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
10^{7}
small numbers a = 5x10^{-8} to 1 (so the exact answer is 1.5).

- If you add 1 + a + a + ... in single precision, what do you
expect to get for the total?

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

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

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

- In a terminal window, use the makefile make_demo1a to try and compile it ("make -f make_demo1a"). 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!).
- Run indent on the file ("indent demo1a.cpp") and check that it is now nicely formatted.
- Challenge: Predict the result for 10
^{8}additions (instead of 10^{7}):

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

The next four tasks (Bessel 1-4) 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. Read through the discussion before going further.

- Copy make_area to make_bessel and open make_bessel in an editor. (If you are in the session_02 subdirectory, you can use "cp ../session_01/make_area ." to create a new copy where you are or "cp ../session_01/make_area make_bessel" to make a copy directly; do you understand what all the "."'s mean? If not, ask!) Change all "area"'s to "bessel"'s. (Can you do this all in one step with the editor?) Run the makefile using "make -f make_bessel".
- If you run "bessel", 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?
- 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.

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

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

- 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?)

Mathematica doesn't have built-in spherical Bessel function routines but you can define them in terms of the more general Bessel functions.

- Using the Help Browser (with Master Index), find the name of the general Bessel function and how to use it.
- Define a function of x and l to calculate
spherical Bessel functions, using:

j_{l}(x) = (pi/2x)^{1/2}J_{l+1/2}(x) - How can you test that your definition is correct?
- Try to print 16 digits for l=10 and x=36.2. [Hint: Look up NumberForm and ScientificForm in the Help Browser.]
- How can you reveal the small x behavior of the spherical Bessel function?

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 z_{c} = 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?

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

Mathematica has A LOT of useful built-in functions, IF you can find them. The Help Browser (under the Help Menu) is your key; it's like the Google for Mathematica. Here we'll try to solve some problems to exercise your ability to use the Help Browser. Using the "Master Index" search is often the most efficient method.

- What is the atomic weight of uranium?

- How much is a boardfoot in MKS units?

- Is 15485863 a prime number? How many smaller primes are there?

furnstahl.1@osu.edu