6810: 1094 Activities 4

Online handouts: GSL eigensystems documentation; singular integrals eigen_test.cpp and derivative_test.cpp, pointer_test.cpp, and qags_test listings.

Your (always optimistic!) goals:

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 the source of error 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*x^beta.

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*x^beta and its derivative, getting alpha and beta from the passed params_ptr (again, see pointer_test.cpp for an example in 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. Use the gnuplot plot file with small modifications to generate a postscript file, which you should attach.
Did you find the same slopes with your new function? Why are the slopes the same but the intercepts different?

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 task 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.)
```
MyHilbertMatrix[n_] := Table[1/(i + j - 1), {i, 1, n}, {j, 1, n}]
Amat = MyHilbertMatrix[4]
MatrixForm[Amat]
Eigenvalues[Amat] // N
```
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?). Describe your method. (For a plus: attach a graph showing the scaling, with a fit.) [Note: you'll need the matrix dimension to be 200 or more.]

Extra: 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):
1. Brute force (unmodified)
2. Subtracting the singularity
3. Changing variables
What do you find?

6810: 1094 Session 4. Last modified: .

furnstahl.1@osu.edu