Run and deconstruct a code comparing numerical derivative methods,
including Richardson extrapolation.

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.

Modify the integration code from Session 3 to use a different function,
then evaluate an integral with a singularity.
This is one of the Assignment 2 problems, but you can get a
head start here.

Please work in pairs (more or less).
The instructors will bounce around 1094 and answer questions.

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, which
is a generalization of the derivative_test_simple.cpp code from
Session 3. It uses "void pointers" to pass extra parameters (like alpha)
to the functions.

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. Add code to print column headings.

Edit the plot file to add the corresponding plots and fits for
the central derivative and extrapolated central derivative
approximations.
Print the postscript version of the plot and attach it.

In Session 3, you saw
the functions forward_diff and central_diff (notice that
the function name is passed as a pointer).
Explain the slope of the extrap_diff graph (see the discussion
in the Session 4 notes).

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. Describe
how you would get an even higher-order result.

What is happening on the left side of the graph
(smaller mesh sizes)? Why are the slopes the same?

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 test_function.
Your job is to modify the code so that two variables,
alpha and beta, are passed to the function alpha*xbeta.

Start by defining a structure with the two parameters
alpha and beta (see pointer_test.cpp for an example).

Modify test_function and
test_function_derivative for alpha*xbeta
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).

Test the numerical derivatives at x=2.
(You should be able to use the gnuplot plot
file with small modifications to generate a postscript file,
which you should attach.)

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 set-up. In
particular, the program in eigen_test.cpp creates a Hilbert matrix (as
described in the Session 4 background notes) of user-specified 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 zip file 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 the eigenvalues of a
4 by 4 Hilbert matrix.
Bring up Mathematica and generate the matrix Amat using
the following Mathematica code. (You could
also use HilbertMatrix after loading
LinearAlgebra`MatrixManipulation`; see the Help Browser for more
information.)

The last two commands present the matrix in the conventional way
and calculates the 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?). Note: you'll need the matrix dimension to be 100 or more.

More Complicated Integrals

Choose one of the integrals with singularities from the
"Integrals with Singularities or Discontinuities"
handout [Eq.(9) is a good choice!].

Determine an "exact" answer analytically or using Mathematica
(see instructors for assistance).

Modify your integration code from Session 3
so that you can analyze the integral with one
of the integration rules several ways, using an error plot to compare
(these may not all apply to a given integral):