Physics 6810: Assignment #3

Note that extrap_diff calls central_diff with two different h's
and then extrap_diff2 (which you are to write) calls extrap_diff twice
(and NOT central_diff twice). Your error plot should show a slope of
h6.

When writing an adaptive code (first bonus problem),
you can assume that the
round-off error part of the error is given by the machine precision
divided by h (as described in Chapter 8). In practice this may
overestimate the actual error.

If you know that the relative error scales in the "approximation
error region" as b*h2 for the central difference formula.
How can you determine b from evaluating the derivative at two
different h values (without knowing the exact answer)?

Do not confuse the eigenvectors from the eigen_basis.cpp code
with the ones from the eigen_tridiagonal.cpp code. In the latter,
the eigenvector was itself the wave function on a discrete coordinate
mesh. In the former, we have expanded the wave function as a
sum over coefficients times harmonic oscillator wave functions.
The eigenvectors are those coefficients. So to reconstruct the
corresponding wave function u(r), for each value of r you need
to find the sum of the coefficient from the eigenvector times the
corresponding harmonic oscillator basis function. This requires
understanding the code well enough to find where those basis functions
are generated and deciding what values of r to use.

So you will need an additional loop over r.

Look where the Hamiltonian integrand is calculated to
find the harmonic oscillator functions.

Be sure to note the "common problems" in the next bullet.

Here are some common problems with the eigen_basis problem:

The eigenvector coefficients are labeled from j=0 to
j=dimension-1, but the corresponding n_j values for the
harmonic oscillator wave functions [called using
ho_radial(n_j,l,b_ho,r)] has n_j = j+1. So j=0 has n_0=1,
and so on.

We're working with "s-waves" here, so l=0.

You have the choice of r under your control.

You only need to get the first eigenvalue (i=0 in the code),
since we want the ground state.

The overall sign of the radial function is arbitrary (in
quantum mechanics, only the square of the wave function is
physical), so if your graph for some dimension size comes out
negative, you are free to flip it to positive (e.g., in
gnuplot, with: plot "output_file.dat" using ($1):(-($2)).

If you do the Coulomb potential (recommended), the exact
answer is u(r) = 2.*r*exp(-r), which you can plot on gnuplot
for comparison. (In gnuplot, "plot 2*x*exp(-x)" will do it.)

The eigen_basis
results are most interesting is you pick a value for the
harmonic oscillator parameter b that
is not optimal. Or, if possible, compare results with
two different values of b.

Assign variables to the minimum, maximum, and increment values for
the radial value r. For example, rmin=0., rmax=10., delta_r=0.05
are reasonable choices.

Plot the results for different dimension bases on the same gnuplot
plot. This is made much easier by using a plot file that uses a
different datafile for each basis size, rather than trying to make one
datafile with all of the results. (You can rename the output
files from your code each time you change the basis size.)

When commenting on the nature of the convergence, you might note
what part of the wave function is reproduced first (e.g., when is the
tail reproduced?). If you chose a different value for b, how might the
convergence change? (What is most important to get a fit with
dimension size 1?)

For the second bonus problem,
think about how you measure the goodness of fit when
you fit a straight line to some data points. Can you do something
analogous when you have a function of r like a wave function?