6810: 1094 Activities 2
Online handouts:
gnuplot minimanual; gnuplot example; listings of several codes.
Your goals (in order of priority; you will
finish some in PS#1):
 Look at a roundoff error demonstration with quadratic equations
 Discuss a specialcase problem to help understand summing
in different orders
 Debug sample codes; make some plots with gnuplot
 Assess different methods of calculating spherical Bessel functions,
including with GSL
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 your working directory (e.g., for cygwin go to
the shared disk with cd /cygdrive/z, see Activities 1 for
creating one),
download session02.zip from the 6810 homepage and unzip it
in your directory.
RoundOff Errors with Quadratic Equations
This section assumes you have read the background notes for Session 2,
which describe roundoff errors and the quadratic equation.
 Review the discussion about roundoff errors and the quadratic
equation in the Session 2 notes with your partner.
Make sure you
go over the pseudocode for the test case and
understand it.
 Look at the code quadratic_equation_1a.cpp in an editor.
It represents the pseudocode as it has just been
typed in, mistakes and all.
Note that it is not indented consistently. Why might we care
about formatting?

Now try compiling with "make f make_quadratic_equation_1a".
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 below 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 to turn them on).
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:
 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
(of x1, x1p, x2, x2p)
are the most accurate? Explain why the others are not as accurate. (Hint: check the Session 2 notes.)
 Now look at quadratic_equation_2.cpp, which is the second
pass at the code. Check the comments for a description of changes.
List below up to three C++ lines from the code that you don't understand
completely (it may be most of them or it may be none!):
 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 from the Session 1 codes).
 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 the Smith1094 printer using
lpr P //PHYPCS/smith1094 quadratic_eq.ps
in cygwin or
lpr P smith1094 quadratic_eq.ps on Linux,
and hand it in at the end (all titles are not
necessary here, but will be in the future!).
 Is the plot qualitatively and
quantitatively consistent with the analysis in the notes? Explain
in a few sentences. [Hints: What does a straight line mean? What does the slope tell you?]
Summing in Different Orders (see Session 2 notes)
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/(N1) + ... + 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).
For the first three questions, predict the answer before running
the code.
 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.")
 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" here for a program that would compare
these two ways of summing.
Then look at the file order_of_summation1a.cpp in an editor;
does it look like an
implementation of your pseudocode?
 In a terminal
window, use the makefile make_order_of_summation1a to try and compile it
("make f make_order_of_summation1a").
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 order_of_summation1a.cpp") and check that
it is now nicely formatted.
 Challenge: Predict the results 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).
Bessel 1: Making Another Plot with Gnuplot
The next three tasks (Bessel 13) build on the discussion of spherical
Bessel functions in the Session 2 notes.
Skim this 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!)
Follow the instructions in the makefile to convert it to make_bessel.
Run the makefile using "make f make_bessel".
 If you run bessel.cpp,
you'll generate a data file "bessel.dat".
By looking at the code, figure out what the columns
in the file mean. What Bessel function is being output?
What are the columns?
 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: Error Assessment (complete in PS#1)
 Modify the output statements in bessel.cpp to increase
the number of digits in the output (look back at previous
codes for clues how to do this).
 Modify the code bessel.cpp (give it another name)
to compare upward and downward recursion methods by adding
a new column to the output that calculates the relative
difference between the results for each x. We'll define the
relative difference of a and b as ab/(a+b). Google to
find out how to take an absolute value in C++.
 Make an appropriate plot of the error vs. x and interpret the
different regions of the graph. Try with and without log scales.
Bessel 3: Using the GSL Scientific Library (complete in PS#1)
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?)
6810: 1094 Activities 2.
Last modified: 07:50 am, January 13, 2017.
furnstahl.1@osu.edu