Physics 6810: Assignment #3

• The updated version of derivative_test.cpp is available here.
• 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 h⁶.
• 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*h² 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?