Numerov Algorithm with python

 

Numerov’s method

Numerov’s method¹ (also called Cowell’s method) is a numerical method to solve ordinary differential equations of second order in which the first-order term does not appear. It is a fourth-order linear multistep method. The method is implicit, but can be made explicit if the differential equation is linear. Numerov’s method was developed by the Russian astronomer Boris Vasil’evich Numerov.

Given the differential equation of this form

where

and are continuous function on the domain . To derive the Numerov’s method, we first discretize the interval using equally spaced grid, where the grid space is . We then proceed with the Taylor expansion of the function around the point

:

Similarly, we have

By summing the two equations, we have

To get an expression for

, remember that

. By repeating the same procedure we did above

If we substitute this into the proceeding equation, we get

or equivalently

This yields the Numerov’s method if we ignore the term of order

.

where

. It is easy to show that

With this relation, we can now find the value of

given two previous values and . All we have to do to find an approximate solution is to find the values and at the first two points and , then we use the Numerov’s method to find the value of all following points in the interval

.

Case 1. Hydrogen-like atom

The Schrödinger equation for hydrogen-like atom is ²

Here, we assume that there is only one electron and thus neglecting the Coulomb interaction between the electrons. By variable separation, the wavefunction can be written as

, where is the spherical harmonics and

follows the equation

Introducing the dimensionless variables

Å

where

is the Bohr radius and

is the Rydberg energy.

The equation for

then becomes

The boundary conditions are

1. From known energy eigenvalue to the wavefunctions

As is well known that the energy levels of hydrogen-like atoms are

Let’s first show that given the correct energy levels, we can obtain the corresponding wavefunctions.

CLICK TO SHOW THE CODE!

In the code above, we have defined two integration method to integrate the radial Schrödinger equation

• radial_wfc_scipy uses odeint implemented in scipy.³

• radial_wfc_numerov uses the Numerov’s method.

• We implemented forward/backward integration scheme in radial_wfc_scipy function, where forward integration means integrating from

and backward integraton starts from

.

Of course we can not integrate from

, we used a large enough cutoff and assumes

• .

CLICK TO SHOW THE CODE!

Figure 1. Left: forward and backward integration of the radial Schrödinger using scipy.integrate.odeint. The solid cyan line shows the exact solution. Right: the difference between the exact solution and the numerical results. Two different integration schemes are used, i.e. odeint in scipy and Numerov method.

We then try to get the radial wavefunction of hydrogen (

) with quantum number

. Some remarks here:

• It turns out that the forward integration of the equation is not stable, as can be seen in the left panel of Figure 1.

• It is better to use logarithmic mesh for radial variable rather than linear. Radial functions need smaller number of points in logarithmic mesh. For example, use np.logspace(1E-6, 2, 500) to generate the radial grid and pass to the radial_wfc_scipy function. Note that Numerov’s method only support linear mesh, as can be seen in the method section.

One can show that the Numerov’s method also works fine for higher excited states. For example, below we compare the radial wavefunction

for

obtained from Numerov method with the exact solutions. The script can be found here.

Figure 2. The lowest few radial wavefunctions

with

.

2. Unknown energy levels

If we don NOT know the correct energy levels, we can utilize the boundary condition to get the energy. Recall that the boundary condition is

where

is large enough so that . If corresponds to the correct energy level, when we integrate the radial wavefunction from backward using Numerov’s method, we should get . A deviation of the energy from will result in

. Therefore we can use the so called “shooting” method to get the correct energy. The basic procedure is

1. Start wite a guess energy, e.g.

. In practice, the guess energy should be smaller than the smallest potential energy. In our case, should be smaller than . With the guess energy , integrate the equation and get the value of the function at , which we will denote as . Meanwhile, Set another energy

.

Increase

by an amount and get a new energy . Integrate the equation to get the value at , which we will denote as

.

Go back to step 2 if

.

At this step, we should have the correct energy enclosed in

4. . Use root finding method, e.g. scipy.optimize.brentq to get the correct energy.

CLICK TO SHOW THE CODE!

References

1. Numerov’s method ↩

2. Solving Hydrogen Atom with Python ↩

3. scipy.integrate.odeint
   

 

http://staff.ustc.edu.cn/~zqj/posts/Numerov-Algorithm/