780.20: 2082 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:

Please work in pairs (more or less). Dick and Daniel will bounce around 2082 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.

  1. 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.)
  2. 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 (find "seascape" in the menus to flip the graph) and then print it (you can also print directly from gv):
       gv derivative_test_plt.ps &
       lpr -Psmith2076 derivative_test_plt.ps
    (For a color plot, "set term postscript color" in the plotfile and print on smith2097c.)
  3. 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?
  4. 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 higher-order result.
  5. 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*xbeta.

  1. Start by defining a structure with the two parameters alpha and beta (see pointer_test.cpp for an example).
  2. Modify funct and funct_deriv 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).
  3. 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).
  4. 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 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?).

  1. 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.
  2. 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]
    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.
  3. 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!).
  4. 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.
  5. 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?).

Leftover Tasks from Session 3

Come back to these after you complete the regular Session 4 tasks. 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!

  1. 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?
  2. GSL Scientific Library Again. Look at the qags_test.cpp example (you have a printout AND the code in Session 4).
  3. More Complicated Integrals. This is one of the Assignment 2 problems.

780.20: 2082 Session 4. Last modified: 08:15 am, January 17, 2005.