780.20: 1094 Session 4
Handouts: Session 4 background notes;
eigen_test.cpp and derivative_test.cpp, pointer_test.cpp, and qags_test
printouts;
GNUPLOT plot files;
GSL eigensystems documentation.
Your (always optimistic!) goals for today:
 Finish some leftover Session 3 tasks
(the rest become part of Assignment #2).
 Run and deconstruct a code comparing numerical derivative methods

Modify the derivative code
to apply to a different function with more than one passed
parameters (structures and pointers!).

Run a sample program to find eigenvalues and eigenvectors of
a matrix
and determine how the running time scales with the size of the
matrix.
Please work in pairs (more or less).
The instructors will bounce around 1094 and answer questions.
Leftover Tasks from Session 3
Work on these at most for the first 40 minutes of the period.
Note that any unfinished tasks from Session 3 are covered in Assignment
2. Or, if you get through them in class, you'll have a headstart
on your homework!
 Finding the Approximation Error from a LogLog Plot.
If you have not reviewed your error plot with Prof. Furnstahl,
please try to do so today.
 Coding an algorithm. Note the discussion in
the Session 4 background notes on approximation errors. Were you
careful about how many points were in each 3/8's interval?
 GSL Scientific Library Yet Again. Look at the
qags_test.cpp example (you have a printout AND the code in Session
4).
 More Complicated Integrals.
This is one of the Assignment 2 problems, but you can get a
head start here.
Numerical Derivatives and Richardson Extrapolation
Take a look at the derivative_test.cpp handout. The derivative_test.cpp
code evaluates the numerical derivative of a defined function (funct)
four ways: forward derivative, central derivative, extrapolated
central derivative, and with a GSL routine.
The function is defined according to the GSL conventions.
 Compile and link derivative_test.cpp (using make_derivative_test).
Run it to generate derivative_test.dat and then use
derivative_test.plt in gnuplot (see the handout on using a plot
file) to verify the figure on the back of the handout.
(You'll have to look at the code to identify the columns in
derivative_test.dat.)

Edit the plot file to add the corresponding plots and fits for
the central derivative and extrapolated central derivative
approximations.
Look at the postscript version of your plot with gv and then
print it (you can also print directly from gv):
gv derivative_test_plt.ps &
lpr Psmith1011s derivative_test_plt.ps
(For a color plot, "set term postscript color" in the plotfile and print
on smith2097c.)
 Look at the functions forward_diff and central_diff (notice that
the function is passed as a pointer).
Derive the expected approximation error in each case
(without peeking at notes or handouts, if possible)
by making
a Taylor expansion about x of each of the expressions returned by these
functions. Are the errors consistent with your graph?
 Can you account for the slope of the extrap_diff graph?
This method is an example of "Richardson extrapolation," in which
you use calculations at two different mesh sizes to derive a much
better estimate than either one individually. Discuss with your
partner how you would get an even higherorder result.
 What is happening on the left side of the graph
(smaller mesh sizes)? Can you explain the slopes?
How do you determine the optimal mesh size from
the graph for each approach?
Pointer Games
This exercise is practice in writing or modifying code based on
examples.
Take a look at the pointer_test.cpp handout. It gives
examples of how to pass several types of variables to a function using
the void pointer named params_ptr. In derivative_test.cpp,
a double named alpha is passed to the function funct.
Your job is to modify the code so that two variables,
alpha and beta, are passed to the function alpha*x^{beta}.
 Start by defining a structure with the two parameters
alpha and beta (see pointer_test.cpp for an example).
 Modify funct and funct_deriv for alpha*x^{beta}
and its derivative, getting alpha and beta from the passed
params_ptr (again, see pointer_test.cpp for an example in f_struct
or f_osu_parameters).
 Modify the main program in derivative_test.cpp to load alpha and
beta into your structure (see the main code in pointer_test.cpp for
an example). Simply choose values for alpha and beta (2. and 3./2., for
example).
 See how well the numerical derivatives work at x=2.
(You should be able to use the gnuplot plot
file with small modifications.)
Linear Algebra with GSL Routines
The GSL library has many functions defined to set up and manipulate
vectors and matrices. To do so, it defines various structures such as
"gsl_matrix" and "gsl_vector", and functions such as "gsl_matrix_set" to
set the value of an element in the matrix.
It is all a bit intimidating at first, so we'll
take a look at a basic example to see the general setup. In
particular, the program in eigen_test.cpp creates a Hilbert matrix (as
described in the Session 4 background notes) of userspecified dimension
and then calls a routine to find its eigenvalues and eigenvectors.
An important issue with numerical computations, and linear algebra in
particular, is how the computation time scales with the size of the
problem.
The program in eigen_test.cpp includes two calls to the "clock" function,
before and after the eigenvalue routine, to time how long the
routine takes to run.
Your job is to figure out how the time scales with the size of the
matrix (e.g., does it go like a power of the dimension of the matrix?).
 The session 4 tarball from the webpage
should contain eigen_test.cpp and make_eigen_test.
Compile and link eigen_test using make_eigen_test.
 You should always verify with test cases that a program is
working.
We'll use Mathematica to generate eigenvalues of a
4 by 4 Hilbert matrix.
Bring up Mathematica and generate the matrix Amat using Table, since a
matrix in Mathematica is just a nested list. (You could
also use HilbertMatrix after loading
LinearAlgebra`MatrixManipulation`; see the Help Browser for more
information):
MyHilbertMatrix[n_] := Table[1/(i + j  1), {i, 1, n}, {j, 1, n}]
Amat = MyHilbertMatrix[4]
MatrixForm[Amat]
Eigenvalues[Amat]
The last two commands present the matrix in the conventional way
and supposedly finds the eigenvalues.
But you get a cryptic result. Figure out how to get numerical
eigenvalues.
 Run eigen_test with a dimension 4 matrix
(i.e., 4x4) and compare the answers to the ones given by Mathematica.
Try to trace through the code on the printout
to identify what the different
GSL function calls do (you are NOT expected to
understand the calls in detail at this point!).
 Edit eigen_test.cpp and
comment out the section that prints to the screen, so that the
only output is the time the routine takes to run.
 Figure out a way to determine how the execution time scales with the
size of the matrix (e.g., is it a power law? If so, what is the
exponent?).
780.20: 1094 Session 4.
Last modified: 09:50 am, January 18, 2006.
furnstahl.1@osu.edu