*Handouts: excerpt on simulated annealing *

In this session, we'll first consider a toy problem comparing multidimensional minimization with standard downhill routines for local minima to an application of Monte Carlo techniques to find the global minimum: simulated annealing. Then we'll consider the bound-state Schrodinger equation in momentum space, which takes the form of an integral equation.

Standard optimization methods are very good in
finding local minima near where the minimization was started, but not
good at finding the *global* minimum. Finding the global minimum
of a function (such as the energy) is often (but not always) the goal.
One strategy using conventional minimizers is to start the minimization
at many different places in the parameter space (perhaps chosen at
random) and keep the best minimum found.

Here we'll adapt our Monte Carlo algorithm for generating a canonical
Boltzmann distribution of configurations at a temperature T to mimic how
physical systems find their ground states (i.e., the energy minimum at
T=0). At high temperature (compared to characteristic energy
spacings), the equilibrium distribution will include many states. If
the system is cooled *slowly*, then it will have the opportunity
to settle into the lowest energy state as T goes to zero. This is
called *annealing*. If the system is cooled quickly, it can get
stuck in a state that is not the minimum ("quenching"); this is
analogous to the routines we looked at in Session 10, which rapidly go
"downhill" to find local minima.

The strategy here is to simulate the annealing process by treating the function to be minimized as an energy (it might actually be an energy!), introducing an artificial temperature T, generating a sequence of states in a canonical distribution via the Metropolis algorithm. Then we lower the temperature according to a "schedule" and let the system settle into (hopefully!) a configuration that minimizes the energy.

We first consider the simple example of a damped sine wave in one
dimension, which has
many local minima.

`f(x) = e ^{-(x-1)2}sin(8x)`

(This example is from the GSL documentation.) Suppose we start with an initial guess of 15.5; let's see how "downhill" minimization methods compare to simulated annealing.

- Plot the function in Mathematica (see test_function.nb)
so we know what to expect ahead of
time (remember that in general we won't be able to
do this!). Use FindMinimum starting near 15.5 to find the local
minimum. Then find the global minimum (using the graph to identify
a good starting value).
*Answers:*

- The code multimin_sa_compare.cpp (compile and link with
make_multimin_sa_compare) applies the GSL routines we used in a
previous session to find the minimum and then applies a simple
simulated annealing algorithm. Run the program and look at the code to
figure out where the simulated annealing control parameters (starting
and finishing temperatures, maximum step size, cooling rate, # of
iterations at each temperature) are set.
*How well does each method work?*

- Now adjust the simulated annealing control parameters to give it a
chance! The step size should be large enough (or small enough)
so that successes are neither assured nor happen too infrequently
(50% would be ideal!).
The initial temperature should be considerably larger
than the largest typical delta_E normally encountered.
You shouldn't cool too rapidly. Make these adjustments and rerun the
code.
*What did you do to get it to work best?*

- Just to check that the Metropolis algorithm is doing something,
once the code is working, change the line with "random < exp()"
to "random > exp()" and verify that it no longer works.
*What happened?*

- Bonus: Switch the code to use classes defined in previous sessions.

As discussed in the background notes, the bound-state Schrodinger equation in momentum space can be solved as a matrix eigenvalue problem. We have three versions of a code to solve it: two with the matrix in symmetric form, which use GSL routines and LAPACK routines, respectively, and one that uses a LAPACK routine for non-symmetric matrices. We also dynamically allocate arrays using "new" and "delete".

The delta-shell potential is a convenient choice for a test potential because we can easily find the bound-state eigenvalue (numerically) and it takes a simple form in momentum space. It has only one bound state; if we only print out the negative-energy (E<0) eigenvalues then we should find just one. See the background notes for additional details. [Warning: It is not a well-behaved potential in momentum space, so you will be challenged to find good solutions.]

- What are the units of lambda in the delta-shell potential?
[Note: hbar=1 here.]

- For mu=0.5, b=2.0, and lambda=-1.0, use Mathematica
to find the correct eigenvalue
for the (single) bound state in the delta-shell potential by
using the appropriate equation in the background notes. [Hint: FindRoot]
*What is the answer?*

- Look at deltashell_boundstates_gsl.cpp.
Note how the non-gsl arrays
are dynamically allocated and removed using "new" and "delete".
*If we had an int called*`size`and wanted to create an array called`my_array`with dimension`size`, what command would we give?

- Compile and link deltashell_boundstates_gsl.cpp using the corresponding makefile.
Run with lambda=-1.0 but specify on your own
the number of gauss points to use and the upper bound for the
integrations of momentum. The latter may be much larger than you
would guess because the delta-shell potential is actually rather
poorly behaved in momentum space!
*Based on how the potential falls off with k, what do you think is a reasonable momentum cutoff? (What should you do to discover how the potential falls off?)*

- Now look at deltashell_boundstates_lapack.cpp. This is essentially
the same code but it uses a Fortran LAPACK function instead of a GSL function
to diagonalize the matrix. (Note: You have to be on a Linux machine for
this to work.) We'll also use the Intel compiler and the MKL
library. Load them with:
`module load intel-12.1.3-64`

Look at the corresponding makefile to see how the code is compiled and linked.*Can you reproduce the results from the GSL code?*

- Run one of the codes
with different numbers of gauss points and different momentum
cutoffs to try to improve the predicted eigenvalue.
Estimate the precision of your best value for the bound-state
energy.
*How could you improve the numerical accuracy of the calculation?*(For example, would it help to change "job" when calling "gauss"? [Hint, hint.])

- What does the momentum-space wave-function look like (it is
printed out to a file)?
- Does it fall off at large k? Does it oscillate? Is it well behaved
at the origin? Can you account for its characteristics by some simple
arguments?

- The analytic wave function in coordinate space is given in the background notes. Try to calculate the momentum space wave function from this (as discussed in the notes) using Mathematica.
*Is the momentum wave function normalized? What is the normalization condition?*

- Does it fall off at large k? Does it oscillate? Is it well behaved
at the origin? Can you account for its characteristics by some simple
arguments?
- Bonus: Verify that, regardless of the potential's strength, there is only a single bound state.
- Bonus: Modify the code to calculate the l=1 bound state.
- Bonus: Make a version of the code with classes that can use
*either*GSL or LAPACK.

furnstahl.1@osu.edu