*Handouts: Pang excerpt, pendulum power spectrum graphs,
printouts of multifit_test.cpp and
multimin_test.cpp.*

In this session,
we'll finish up session 9 and take a look at nonlinear least-squares
fitting with GSL.

*
Your goals for today (and probably Assignment 4):
*

- Finish up as much of the Session 9 interpolation tasks as you can do in the first hour.
- Try out nonlinear minimization on a demonstration problem.
- Try out nonlinear least-squares fitting on a demonstration problem.
- Think about applying nonlinear minimization to a real physics problem (from the Pang excerpt). [We'll finish this in a later session.]

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

Budget your time so that you get through "GSL Differential Equation Solver" and "GSL Interpolation Routines" in the first hour. If you have time, start "Cubic Splining"; don't worry about finishing it, as we'll do the rest in Assignment 4.

The basic problem is to find a minimum of the scalar function f(xvec) [scalar means that the result of evaluating the function is just a number] where xvec = (x_0, x_1, ... , x_N) is an N-dimensional vector. This is a much harder problem for N > 1 than for N=1. The gsl library offers several algorithms; the best one to use depends on the problem to be solved. Therefore, it is useful to be able to switch between different choices to compare results. Gsl makes this easy. Here we'll try a sample problem. For more details, see the online gsl manual under "Multidimensional Minimization".

The routines used here find local minima only, and only one at a time. They have no way of determining whether a minima is a global one or not (we'll return to this point later).

- Look at the test code "multimin_test.cpp" and compare it to the online gsl documentation under "Multidimensional Minimization". Identify the steps in the minimization process.
- Compile and link the program with "make_multimin_test" and run it.
Your results may differ from what is in the comments of the code.
Adjust the program so that your answer is good to 10
^{-6}accuracy. What is it that is found to that accuracy? (E.g., is it the positions of the minimum or something else?)

- Modify the function minimized so it is a three-dimensional paraboloid with your choice of minimum. Change the code so that the dimension of xvec is a parameter. (Be sure to change ALL of the relevant functions.) Verify that the minimum is still found.
- What algorithm is used initially to do the minimization?
Modify the code to try a couple of the other algorithms.
Which is "best" for this problem?

- In Session 14, we'll return to try to
reproduce the results illustrated in Fig. 4.2 of the
Pang handout (stable structures of
Na
_{3}Cl_{2}^{+}and Na_{3}Cl_{3}). Take a look at that now to see if you have any immediate questions

Using the nonlinear least-squares fitting routines from gsl is similar to using the gsl minimizer (and involves a special case of minimization). Your main job is to figure out from the online documentation and this example from gsl (only slightly modified from what you would find in the documentation) how to use the routines. Here we briefly explore the sample problem, which also illustrates how to create a "noisy" (i.e., realistic) test case.

The model is that of a weighted exponential with a constant background:
` y = A e ^{-lambda t} + b `

You'll generate pseudo-data at a series of time steps t

- Take a look at the multifit_test.cpp
code (there is a printout), run it (compile
and link with make_multifit_test) and figure
out the flow of the program. What is the fitting function to be minimized?
What role does sigma
_{i}play?

- What is the "Jacobian" for? (Remember to check the online
documentation!)

- Modify the code so that you can plot both the initial data and the fit curve using gnuplot. Look up how to add error bars for the data in gnuplot. [Hint: Try "help set style", "help errors", and "help plot errorbars" for some clues.]
- The covariance matrix of the best-fit parameters can be used to
extract information about how good the fit is. How is it used in the
program to estimate the uncertainties in the fit parameters? How do these
uncertainties scale with the number of time-steps used in the fit?
How do they scale with the magnitude of the noise? (I.e., if you
change the "amplitude" of the gaussian noise, how do the errors in the
fit parameters change?)

- Modify the function by adding a fourth parameter and try to fit (and plot!) again. For example, you could change "b" to "b sin(omega t)" with omega another parameter, or "A" to "A cos(omega t)". Your choice!

furnstahl.1@osu.edu