Physics 6730 Assignment 4

In this assignment you will solve the radial Schroedinger equation for the hydrogen atom. See the hydrogen atom handout for a discussion of the equation.

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

Exercise 1 Preliminaries

1. Convert the radial Schroedinger equation, a linear second order differential equation, into two first order equations for y(0) (the wave function) and y(1) (the derivative of the wave function). Use the simplified units of the handout, Eq (3). We call rho r.

   The equation we are solving here has singular points at both endpoints of integration, namely r = 0 and r = infinity. The best strategy is always to integrate away from a singular point. So our shooting procedure involves integrating from both ends toward a common matching point, r = c. At c we must match the values of the function and its derivative resulting from the two integrations. At the starting point of both integrations, we select initial values of the function and its derivative. For each integration the initial value is zero and the derivative is arbitrarily set to 1.

2. On [a,c] we call the solution u₁(x) and its derivative du₁(x). On [c,b] we get the solution u₂(x) and its derivative du₂(x). Now suppose we multiply the solution on [a,c] by a constant K to get v₁(x) = u₁(x). Does it still satisfy the differential equation? What boundary condition would it satisfy at x = a?

3. Show that if u₁(x) and u₂(x) and their derivatives satisfy the matching condition

       u1(c) du2(c) - u2(c) du1(c) = 0
   

   then so do v₁(x) and u₂(x) and their derivatives.

   Note that the same conclusion applies to a change in the normalization of the solution u₂(x). In this way we can always adjust the normalization of our original solution on the right and left so the wavefunction and its derivative are continuous on the whole interval and the square integral of the wavefunction is 1.

Please put your answers in a file setup.txt

Exercise 2 Completing the Code

File	Purpose of procedure
bohr.cc	Driver (main) program. (provided)
bsrf.cc	Finds the zero. (provided)
wavef.cc	Your modified driver. (Incomplete).
hydrogen_derivs.cc	Evaluates the rhs of the ODEs. (You write)
param.h	Defines parameter classes. (provided)
Makefile	Rule for compiling the above. (provided)

Although most of the code is provided, it is important that you understand how it works. Here is a diagram showing the module names and hierarchy of subprogram calls for this code:

			 bohr
			  |
			  v
			 bsrf
			  |
			  v
			matchf
			  |
			  v
			wavef
			  |
			  v
		       NR::rkqs
			  |
			  v
		       NR::rk4
			  |
			  v
		   hydrogen_derivs

1. The main program bohr.cc

   This program sets the parameters of the integration, prompts for the angular momentum, and prompts for a range of eigenvalues in which to search for the eigenvalue. The eigenvalue must be bracketed. This routine calls bsrf to locate the zero of the function matchf (as a function of the energy) and reports the answer, if successful.
   
   

2. The subprogram bsrf.cc

   Adapted from Numerical Recipes. The design of this subprogram should be somewhat familiar to you from last quarter. It is a general purpose routine for finding the zeros of any function of a single variable. It locates the zero of a function using two methods: first the bisection method for a rough determination, and then the regula falsi method to polish. It calls your routine matchf.
   
   

3. The file wavef.cc contains two subprograms: wavef and matchf
   • The subprogram wavef

     This subprogram is built after the pattern of your program anhrkqc from the previous assignment. It integrates any second-order linear eigenvalue problem specified by your subprogram, in this case hydrogen_derivs, starting from a zero value and specified slope at the endpoint a and going to the endpoint b with initial stepsize dr and error scale eps. It returns the value of y and its derivative at the second endpoint. Please note that the notation for the endpoints a and b has local scope only. These are variables set in matchf, depeding on which half of the wave function is being constructed.
     
     

   • The subprogram matchf (incomplete)

     This subprogram calls wavef twice, once to do the integration from a to c and then to do the integration from b to c. It then evaluates the matching expression and returns the value. To complete the code, insert the correct formula for the matching expression in the location specified by comments in the code.

   
   
   
4. hydrogen_derivs

   The subprogram hydrogen_derivs is called by NR::rk4 and wavef. The purpose of hydrogen_derivs is the same as anhderivs in Assignment 3.

Exercise 3 Checking your understanding

1. Although we want to integrate from 0 to infinity, what range did we actually take? Why? How might this choice affect the accuracy of the eigenvalue?
2. What exactly (in terms of solutions to the differential equation) is the function whose zeros are being sought in bsrf?
3. Why did we choose the root finder bsrf, which uses bisection and regula falsi instead of the root finder nrbisect from last semester, which used the combination Newton Raphson and bisection?
4. What internally-set parameters should we consider adjusting to improve the accuracy of the solution?
5. What would happen if we changed the matching point from 1 to 2?
6. What if the matching point turns out to be precisely at a zero of the wavefunction? How does the code handle that situation?

Exercise 4 Running the code

1. Run your program and determine the first three eigenvalues for angular momentum 0 and angular momentum 1 and compare with what you expect (see handout). Put your results in a file results.txt.
2. To plot a wave function use the following command:

       bohr | sort -n > wavef.txt
   

   and hit Crtl-d after the first successful eigenvalue determination. Use gnuplot to plot the contents of wavef.txt. Do this for your third eigenvalue for angular momentum 2. Your result for this eigenvalue was not as good as the others. After comparing with a couple other plots, can you suggest a reason for the larger discrepancy? What should you do to improve the answer? Put your explanation and answers (but not the plot file) in results.txt.