780.20: 1094 Session 10
Handouts: Pang excerpt, pendulum power spectrum graphs,
printouts of ode_test_class.cpp and classes, multifit_test.cpp and
In this session,
we'll take a look at nonlinear least-squares
fitting with GSL.
Your goals for today:
- Take a quick look at some pendulum power spectra.
- Try out and extend the adaptive ode code with classes.
- Try out nonlinear minimization on a demonstration problem.
- Try out nonlinear least-squares fitting on a demonstration
Please work in pairs (more or less).
The instructors will bounce around 1094 and answer questions.
Pendulum Power Spectra [Session 8 Revisited]
Look at the handout with the three rows of plots showing the
time dependence of a pendulum on the left and the corresponding
power spectrum on the right plotted in terms of frequency = 1/period.
- How can you most accurately determine the single frequency in the
first row from the graph on the left?
- Identify two of the frequencies in the second row.
- What is the characteristic of the power spectrum in the
third row that is consistent with a chaotic signal?
Start-Up File for Tcsh Users [skip if you use bash]
If you use tcsh rather than bash (type echo $SHELL to check),
then you want to include a .cshrc.more file in your
home directory rather than .bashrc.
- Edit cshrc.more (no initial period) and look at
the syntax. Keep what
you want and copy to your home directory: cp cshrc.more
~/.cshrc.more (note the period in the final name).
- The command source ~/.cshrc will activate the changes.
- Try the aliases! (Including the Session 8 Rsync exercise
if you skipped it.)
Multidimensional Minimization with GSL Routines
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" (there
is also a handout copy).
Identify the steps in the minimization process. What classes might
you introduce for this code?
- 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
What is it that is found to that accuracy? (E.g., is it the positions
of the minimum or something else?)
- [Bonus] Modify the function minimized so it is a function of three
variables (e.g., x, y, z), namely 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
Take a look at that now to see if
you have any immediate questions.
Nonlinear Least-Squares Fitting with GSL Routines
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/handout documentation and this
example from GSL (only slightly modified from what you will 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 ti
(in practice, ti is taken to be 0, 1, 2, ...). Noise is
added to each yi in the form a random number distributed like
a gaussian with a given standard deviation sigmai.
We'll revisit the GSL functions that generate this random distribution in the
next session; for now just accept that it works.
- 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 sigmai play?
- What is the "Jacobian" for? (Remember to check the online
- 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.] Attach a plot.
- 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?)
- [Bonus: come back to this.]
Modify the function by adding a fourth parameter and try to fit
(and plot!) again. For example, you could change "A" to "A cos(omega t)"
with omega another parameter. Your choice!
GSL Adaptive ODE Solver Revisited
In Session 9, we took a quick look at ode_test.cpp, which
implemented the GSL adaptive differential equation solver on the
Van der Pol oscillator. Here we look at a rewrite that uses classes.
- The details of ode_test_class.cpp and the Ode
and Rhs classes are described in detail in the Session 10
notes. Look at the printout while you go through the notes.
What questions do you have?
- Add two additional calls to evolve_and_print so that
all three initial conditions from Session 9,
[x0=1.0, v0=0.0], [x0=0.1, v0=0.0], and [x0=-1.5, v0=2.0], are
generated with the same run. Did the three output files
get generated? (You can try plotting to check correctness.)
- Add another instance of the Rhs_VdP class called
vdp_rhs_2 with mu=3 and generate results for the same
initial conditions. Plot these. Is it still an isolated
- [Bonus. Come back to this if you finish everything else.]
Add another class Rhs_Pendulum that implements the
pendulum differential equation with natural frequency omega0 and
a driving force with amplitude f_ext and frequency omega_ext.
(You will want to modify or replace evolve_and_print.)
780.20: 1094 Session 10.
Last modified: 09:00 am, February 16, 2009.