1

So, we are trying to solve the equation
    -y''(r)-y(r)=E y(r),
or, more specifically, the first excited state of this.

First, a few words on the excited states. The discrete states of a quantum mechanical system can be characterized by so-called quantum numbers, which then serve to determine the energy, momentum, angular momentum etc. of the system. Suppose a system is completely determined by its energy, or its energy quantum number. This applies well to the states in question here, i.e. the lowest energy states of the hydrogen atom. The energy levels of the system, in this case the electron, are of the form
f(k), where k is the energy quantum number. The range of possible values of k determines the range of the energy values.

Usually, in numerical calculations, we don't know the function f. From the outset of the problem, from the zero of the Coulomb potential, we can however deduce that E<0, so as k increases, the energy approaches zero from below. The ground state energy was found to be -1, so all (an infinite number of) excited states are in the interval -1<E<0. There is, a priori, no way of telling where the 'second' state is. The only way of finding out is to calculate the values of the function Δ(E) in sufficiently small steps, starting from -1 and see where its next zero lies.

What is the significance of the function Δ, defined as
    Δ(E)=-,
where (r) is the solution in the interval <r< and (r) gives the solution for <r<? Here is the point corresponding to the origin in our numerical calculations, and corresponds to infinity (actual zero and infinity are hard to handle in numerics). The solution of the above equation starts from the formulation of the boundary conditions, in this case they are y()= and y()=, where and are some constants, which are determined by matching of the partial solutions and the normalization (for the beginning part of the solution, these constants can be set to unity). Corresponding to these boundary conditions, the differential equation has solutions only for a certain set of values of E. A possible solution y of the equation, corresponding to some (guessed) value of E, has to satisfy the conditions that y and y' are continuous, i.e. ()=() and '()='(). Dividing these equations yields the condition that Δ(E)=0. What this all means is that if Δ(E)≠0, there are no solutions corresponding to this value of E.

Alright, let's get on with the solution. First, we'll input some definitions. The crucial ones here are the values of r demarcating the different regions. With an inappropriate choice of , the function Δ might behave in a way making it difficult to find the zeros. Again, there is no way of telling what is an appropriate choice, but as a rule, you should make large enough and not too large. In the second excercise, a better way of finding the zeros of Δ is presented. Note also that the constants and have been set equal to one in the boundary condition functions y0 and y1.

In[1]:=

Then, give some off-the-hat values to E and calculate the corresponding value of Δ. You can do this simply by hand, setting one value after the other, or devise some neat automatic method. For example, you could give E the values -0.9,-0.8,... ,-0.2 and discover that Δ(E) changes sign in the interval -0.3<E<-0.2.

In the following, diff1=Δ(e1) and diff2=Δ(e2).

In[6]:=

Now, since diff1 and diff2 have opposite signs, it has to be the case that Δ(E)=0 for some value of E for which e1<E<e2. We will search a more accurate approximation of the eigenvalue by dividing the interval [e1,e2] in half and calculate Δ(), from the sign of which we can then deduce the location of the eigenvalue with double precision, i.e. it is either in the interval [e1,] or [,e2]. This process of dividing the interval in half is then repeated until some predefined, desired accuracy is reached.

To decide which of the subintervals to examine next, we calculate the product diff*diff1, where diff stands for Δ(). If this product is negative, it means that diff and diff1 have different signs and the eigenvalue must lie between e1 and . If it is positive, the eigenvalue is in the other interval.

In[19]:=

This solution can be verified by comparison with the analytical solution, =, i.e. =-.

Then we can construct the wave function. A preliminary solution has been calculated above: y(r)={

In[23]:=

We'll start the search of the constant from the normalization condition. We'll denote y(r)={

1=|ψ(r)dV=dΩdr =4π dr

= dr=dr+dr=(dr+dr)

=

()=()

=

In[26]:=

Then we calculate the normalization coefficient. Note that the actual wave function is the function y(r)/r. That is plotted in the following.

In[30]:=

This can then be compared to the analytical solution. In Wikipedia, there is some visualization of the orbitals, and you can do a comparison if you wish. http://en.wikipedia.org/wiki/Hydrogen_atom#Visualizing_the_hydrogen_electron_orbitals

Created by Mathematica  (April 10, 2007)