Physics 6730 Assignment 7

For each of the following exercises create the specified files with your answer(s). Submit your homework using the course submit utility.

Exercise 1 Fitting a single peak

The program ~p6730/examples/a07/spectrum1 fits synthetic spectral data to a Gaussian peak plus a constant background. 
    y(w) =  a1 + a2 * exp[-(w - a3)2/(2*a42))
The data could represent the result of the measurement of an absorption spectrum. A typical research objective would be to determine the amplitude a2, frequency a3 and "peak width" or standard deviation a4.

The program uses the Levenberg-Marquardt (LM) method to fit the data. This methods are described in "Numerical Recipes" in section 15.5. Please read this section to learn more about what the algorithms are doing.

The program source files and Makefile are also found in the directory ~p6730/examples/a07. You are to compile the code and run it.

To run the program spectrum1, use the command 

  spectrum1
You are prompted for the following input values:
  1. Names of input and output files.

    The input file for the first fit is peak1. The output file can be whatever you choose. The output file is designed for input to the axis plot program. It includes a comparison of the input data points and the fitted points.

  2. Initial parameter values (a1, a2, a3, a4).

The goal of this exercise is to study the sensitivity of the Levenberg-Marquardt method to the choice of input parameters. So try fitting the data in peak1 with the following starting parameters:

Attempt #1 

 	a1 = 0.1
	a2 = 2
	a3 = 130
        a4 = 2
Then run it with the following choice

Attempt #2 

 	a1 = 0.2
	a2 = 1
	a3 = 175
        a4 = 2
Copy your standard output to spectrum1.txt to hand in.

Then answer the following questions in the same file.

  1. What is the confidence level of the fit resulting from the first attempt? The second? (Use the course utility conf_level.)
  2. What is the relative error (fractional error) in the amplitude parameter a2 for the first attempt? For the second?
  3. What pair of parameters are the most strongly correlated in the fit from the first attempt? What about the second?
  4. Compare the plots. What went wrong in the second fit? Note: You are welcome to experiment with this fit by, for example, varying input parameters, etc. to see how well the code finds the minimum.

Exercise 2 Understanding the code

At this point it might be a good idea to consult Numerical Recipes. Please put answers to these questions in code.txt.

  1. Make a map of the code, showing the sequence of subroutine calls. (For an example of a map, see assignment 4, exercise 2.)
  2. What is the purpose of subroutine gaussj? Explain in terms of the notation of the lecture notes "Curve Fitting".
  3. What is the purpose of the array lista?
  4. What stopping condition is used?

Exercise 3 Fitting two peaks

We would like to be able to fit two peaks instead of just one, using the formula 
    y(w) =  a1 + a2 * exp[-(w - a3)2/(2*a42)) 
               + a5 * exp[-(w - a6)2/(2*a72))
To do this rewrite the main program spectrum1.cc. (Call your new version spectrum2.cc) and the subprogram peak1.cc (Call it peak2.cc).

Exercise 4 Running your code

Examine the data in the file peak2 and try fitting it with two peaks, one around 140 and the other around 150. You should be able to get a good fit. You may use my code in ~p6730/examples/107/spectrum2 to verify that yours is working properly.

In a file spectrum2.txt give the standard output from your fit to two peaks. Follow that with the answers to the following questions:

  1. What is the confidence level of your fit?
  2. What are the relative errors in the amplitudes of the two peaks?

Now for a little review of our in-class exercise: In the file peak3 you will find a synthetic spectrum generated from two peaks: 

       amplitude   freq     width   
          1        140       8
          .6       147       6
If you look at the plot, it looks like one peak, however. We will pretend we don't know whether there are supposed to be two peaks or one, so we fit this data using two models listed below. However, to speed things up, you may use the above information to help extract two peaks from the data. 
Model #1  Two peaks
Model #2  One peak
Include the standard output from your two fits in the same answer file spectrum2.txt. Then answer these questions in the same file:
  1. Do the best fit parameters agree with the values used in creating the synthetic data? (Use the errors on the parameters.)
  2. Compare the confidence levels of the two fits. Are they good fits?
  3. If I hadn't told you the data were generated from a model with two peaks, which model would you choose? Why?
This page is maintained by:
Carleton DeTar Mail Form
Last modified 7 March 2004