Int. J. Adv. Appl. Math. and Mech.

2(2) (2014) 102 - 109 (ISSN: 2347-2529)

Journal homepage: www.ijaamm.com

International Journal of Advances in Applied Mathematics and Mechanics

Bound state eigenfunctions of an anharmonic oscillator in

one dimension: A Numerov method approach
Research Article

Balkrishna P. Shah ∗
Department of Physics, Faculty of Science, The M.S.University of Baroda, 390002, Vadodara, India.

Received 04 October 2014; accepted (in revised version) 24 November 2014

Abstract: Numerical methods are useful when analytical solutions are difficult to derive. An oscillator is described by the
second order differential equation. In this article the numerical method for the special case of differential equa-
tion is applied for the solution of the wave function of a harmonic oscillator quantum mechanically in classical
as well as non classical region. The Numerov’s method is applied for evaluating the wave function with initial
condition. The computational flow is explain schematically. The method is further applied to anharmonic oscil-
lator. The evaluated wave function is plotted along with harmonic potential. The possible solution leads to the
bound state energy eigen values. The inward and outward integrations are executed to achieve smooth matching
at boundaries. The exact bound state solutions are obtained by the numerical integration corresponds to n = 0
to n = 4. The numerical approach gives the symmetric and antisymmetric wave functions. Computational issues
concerned is energy convergence, which is controlled by the step size h , range of potential xm and lower bound
guessed eigenvalue as well. The present results obtained are also compared with that of Bahttacharya.

MSC: 34L20 • 65L12

Keywords: Numerov method • Anharmonic oscillator • Inward integration• Outward integration

c 2014 IJAAMM all rights reserved.

1. Numerov’s method
One class of differential equations which often occur in physics are second order differential equations which are
independent of the first derivative. An example is the Newton’s equation of motion of a particle under the influence
2
of an external force m dd t x2 = f (x , t ), where m is the mass of the particle and f (x , t ) is the external force. Another
example is Schrödinger equation of a quantum system. For a quantum particle of mass m in a potential V (x ), the
Schrödinger equation is

ħ2 2
h
− ∇ ψ(x ) + V (x )ψ(x ) = E ψ(x ), (1)
2m
where ψ(x ) is the wave function of the system and E is the energy eigenvalue. Numerov’s method is an elegant
scheme to solve differential equations which belong to this class[1,2]. The general form of this class of differential
equations is

d2y
= U (x ) + V (x )y . (2)
d x2
The Taylor series expansion of the solution y at (x − h ) about x is

d y (x ) h 2 d 2 y (x ) h 3 d 3 y (x ) h 4 d 4 y (x )
y (x − h ) = y (x ) − h + − + − ..., (3)
dx 2 d x2 6 d x3 24 d x 4

∗ Corresponding author.
E-mail address: [email protected]

102
 Balkrishna P. Shah / Int. J. Adv. Appl. Math. and Mech. 2(2) (2014) 102 - 109

and similarly, for the solution at (x + h )

d y (x ) h 2 d 2 y (x ) h 3 d 3 y (x ) h 4 d 4 y (x )
y (x + h ) = y (x ) + h + + + + ..., (4)
dx 2 d x2 6 d x3 24 d x 4
For simplicity of notation, denote the solution y (xn ) at the n th value of the independent variable xn as yn and add
Eq.(3) and (4), which gives
 2 2   
h d yn  h 4 d 4 yn 
yn+1 + yn −1 = 2yn + 2 + + O (h 6 ) (5)
2 d x 2  xn 24 d x 4  xn

The form of the equation hints the order of the error, which is O (h 6 ). For further calculations define d 2 y /d x 2 | xn =
Fn and after rearranging

2 h 4 d 2 F 
yn+1 = 2yn − yn−1 + h Fn + + O (h 6 ). (6)
12 d x 2  xn

Replace d 2 F /d x 2 in Eq.(6) with the central difference approximation

Fn +1 − 2Fn + Fn −1

d 2 F 
= Fn00 = , (7)
d x xn
2  h2

which is of order O (h 2 ). Then the solution at xn +1 is

h2
yn+1 = 2yn − yn−1 + (10Fn + Fn+1 + Fn −1 ) . (8)
12
From the general equation in Eq.(2) Fn = Un + Vn yn , using this expression for Fn+1 in the above equation

h2

yn+1 = 2yn − yn−1 + 10Fn + Un +1 + Vn +1 yn +1 + Fn−1 . (9)
12
The variables Vn+1 and Un +1 on the right hand side, from the definitions in Eq.(2), are functions of xn +1 only. After
rearranging the terms in the equation
2
2yn − yn −1 + h12 (Un +1 + 10Fn + Fn −1 )
yn+1 = 2
+ O (h 6 ). (10)
1 − h12 Vn +1

This is the general expression of the scheme. In practice, if the initial value of the independent variable is x0 , the
method is applicable starting from x2 . The solution at this point
2
2y1 − y0 + h12 (U2 + 10F1 + F0 )
y2 = 2
+ O (h 6 ). (11)
1 − h12 V2

However, in this equation y1 is an also unknown. To use the method y1 must be estimated with another scheme to
O (h 5 ) accuracy. The O (h 5 ) accuracy is essential to maintain the O (h 6 ) accuracy of the Numerov method. Evaluating
y1 from Taylor series expansion

h2 h 3 0 h 4 00
y1 = y0 + h y00 + F0 + F + F + O (h 5 ), (12)
2 3! 0 4! 0
here y00 is an intial condition. Further, expand F1 and F2 in the Taylor series about F0 up to O (h 2 ) as

h 2 00
F1 = F0 + h F00 + F + O (h 3 ), (13)
2 0
2h 2 00
F2 = F0 + 2h F00 + F + O (h 3 ). (14)
2 0
Evaluating 6F1 − F2 and rearranging terms in the resulting equation, we get

h2 h2 h 3 0 h 4 00
7F0 + 6F1 − F2 ≈ F0 + F + F . (15)
24 2! 3! 0 4! 0
Collectively, this is an approximation of F0 , F00 and F000 dependent terms in Eq.(12). Substituting this expression in
Eq.(12),

h2
y1 = y0 + h y00 + (7F0 + 6F1 − F2 ) + O (h 5 ) (16)
24
103
 Bound state eigenfunctions of an anharmonic oscillator in one dimension: A Numerov method approach

Noting the earlier definition F (x , y ) = U (x )+V (x )y , denote F1 = U1 +V1 y1 and F2 = U2 +V2 y2 . Substituting F1 , F2 and
rearranging terms the equation is
 
h2 h2 h2
1− V1 y1 + V2 y2 = y0 + h y00 + (7F0 + 6U1 − U2 ) (17)
4 24 24

In a similar way, rearranging the terms in Eq.(11)

   
5h 2 h2 h2
− 2+ V1 y1 + 1 − V2 y2 = −y0 + (F0 + 10U1 + U2 ) . (18)
6 12 12

These two equations Eq.(17) and (18), are linear algebraic equation with y1 and y2 unknowns. In short form, the
equation are

a 11 y1 + a 12 y2 = b1 ,
a 21 y1 + a 22 y2 = b2 ,

where the coefficients and constants are

h2 h2 h2
a 11 = 1 − V1 , a 12 = V2 , b1 = y0 + h y00 + (7F0 + 6U1 − U2 ),
4 24 24
h2 h2 h2
a 21 = −2 − 5V1 , a 22 = 1 − V2 , b2 = −y0 + (F0 + 10U1 + U2 ).
6 12 12

The solution for y1 with the required O (h 5 ) accuracy is

a 22 b1 − a 12 b2
y1 = . (19)
a 11 a 22 − a 12 a 21

Then, with this value of y1 the solution at the next point y2 is calculated from Eq.(17). The Numerov scheme is then
applicable to evaluate the solutions at the other points. An important advantage of the Numerov method is the local
error of O (h 6 ) with one evaluation of U and V per step. Whereas in Runge-Kutta algorithm six function evaluations
are needed to achieve the same local error and two previous values of the solution are needed to calculate a new
one.

Fig. 1. Sequential steps and direction of arrow show the flow of evaluation of the Numerov’s method for simultaneous 1st order
ODE. yn is calculated with the accuracy of the O (h 5 ) from yn−1 with step size h .

2. Eigenfunctions of simple harmonic oscillator

Solution of Schrodinger equation for simple harmonic oscillator potential using Numerov method. Consider realis-
tic case of potential which is continuous functions of position, the simple harmonic oscillator potential V (x ) ∝ x 2 .

104
 Balkrishna P. Shah / Int. J. Adv. Appl. Math. and Mech. 2(2) (2014) 102 - 109

The simple harmonic oscillator is very important in physics as it is the prototype for any system involving oscil-
lations about a point of stable equilibrium. The potential function must have a minimum at a position of stable
equilibrium, thus it can be well approximated by a parabola. Keeping the origin of the x-axis and energy axis at the
minimum, we can write
C 2
V (x ) = x , C = constant (20)
2
For an example a particle moving under the potential experiences a linear restoring force, with C being the force con-
stant of the corresponding linear restoring force. Another important feature is the total energy of the particle. Classi-
cally it can have any value, whereas quantum mechanics predicts only a discrete set of values. The time-independent
Schrodinger equation for simple harmonic oscillator must be solved to find out allowed energy values.

Fig. 2. Simple Harmonic potential well. Solution of Schrodinger equation (Eq. (1)) using Numerov’s method. The ground state
wave function in classical and non-classical region.

The eigenfunctions are expressed in terms of dimensionless variable[3]

h 1/2 x
u = (C m)1/4 /ħ
 

The proportionality constant depends on the properties of the oscillator. In classically allowed region i.e. En < V (x ),
the number of oscillations increases with increasing n , because there are n values of x for which ψ(x ) = 0. The values
of x are the locations of the nodes of ψ. Outside the classical regions the eigenfunctions decreases very rapidly. The
eigenfunction of lowest allowed energy is of even parity.
The time-independent Schrodinger equation for the potential is

ħ 2 d 2ψ C 2
h
− + x ψ=Eψ (21)
2m d x 2 2
1
The Force constant in terms of the classical oscillation frequency is ν = 2π C /m or ν2 2π2 m = C /2.
p

The Schrödinger equation becomes

h 2 d 2ψ d 2ψ 2πmν 2 2
 
2m E
 ‹
ħ
− + 2π2 mν2 x 2 ψ = E ψ or + − x ψ = 0.
2m d x 2 d x2 h2
ħ h
ħ
2πmν
Introducing parameters α = h
ħ h2
and β = 2m E /ħ

d 2ψ
+ (β − α2 x 2 )ψ = 0. (22)
d x2
Consider dimensionless variable for numerical calculations
 ‹ 1/2
2πm 1 C 1/2 (C m)1/4

p
u = αx = x= x
h 2π m
ħ h 1/2
ħ
We have
dψ dψ du du dψ p dψ d 2ψ d dψ du d dψ d 2ψ
 ‹  ‹
= = = α and = = =α
dx du dx dx du du d x2 d x dx dx du dx d u2
Finally the Schrodinger equation in terms of variable u becomes

d 2ψ d 2ψ β
 ‹
α + (β − αu 2 ) ψ = 0 or + − u2 ψ = 0 (23)
d u2 d u2 α

Notice that α and β has unit of length−2 . Thus dimensionless Schrodinger equation is

105
 Bound state eigenfunctions of an anharmonic oscillator in one dimension: A Numerov method approach

d 2ψ 2E
+ (ξ − u 2 )ψ = 0, ξ = β /α =
d u2 hν
or

d 2ψ
= (u 2 − ξ)ψ (24)
d u2
Theoretically

2E
= 2n + 1; n = 0, 1, 2 . . .
hν
where E is the eigenvalues of the simple harmonic oscillator. For n = 0 and given ν ξ = 1, so the good guess value for
ξ may be 0.94 and increment it by very small constant. Similarly for n = 1, E = (1+1/2)h ν ⇒ ξ = 3.0 i.e. starting guess
value can be 2.4. When ξ will be close to the true value the logarithmic derivatives of the wave function obtained
by integration from boundary points to inside using Numerov method match with the required accuracy at the
common point x = x c . In the present case because of the symmetric wave function, it could be origin i.e. center of
the well. This verifies both correct set of wave functions and a corresponding bound state energy.
Equation (24) is integrated by recurrence relation (9) of Numerov algorithm
2
2ψi − ψi −1 + h12 (Ui +1 + 10Fi + Fi −1 )
ψi +1 = 2
(25)
1 − h12 Vi +1

where U = 0, V = u 2 −ξ and F = (u 2 −ξ)ψ(u ). For starting the solution ψ0 , ψ1 are needed, where ψ0 is the initial con-
dition and ψ1 = ε, where ε is a small constant. It is typically sufficient to have about 5 significant figures of accuracy,
so that ε = 10−6 is a reasonable choice. But it reduces magnitude of ψ, so consider ε = 0.001 in the calculation.
It is very important to start the Numerov iterations deep inside a classically forbidden region, where the true solution
decreases. Then integrate out towards lower potential energy region. Because the truncation error accumulation in
forbidden region would not be dominant in the solution. The trick is to integrate the solution in the direction of
increasing magnitude.

2.1. Outward-Inward integration

Starting from the left side of the common matching point with step height +h . The acceptable behavior of the wave
function is ψ(x ) = 0 as |x | → ∞. The potential well is bounded here by the region [−u m , u m ], where its center is
origin. Since the true solution decays exponentially, the wave function can safely be set to zero at a point that is
deep enough in to the classically forbidden region. Assume the distance −2u m and 2u m corresponds to infinite
region where ψ(x ) → 0. Thus outward integration begins from x = −2u m where ψ(x ) = 0 (initial condition).

Fig. 3. Outward integration (left to right) by Numerov method from −∞ = −2u m to the matching point (origin u = 0), with
positive step height.


Starting from x = 2u m , the right side of the well, where ψ(x ) x =2u = 0 (initial condition) up to the matching point
m
with step height −h .

106
 Balkrishna P. Shah / Int. J. Adv. Appl. Math. and Mech. 2(2) (2014) 102 - 109

Fig. 4. inward integration (right to left) by Numerov method from +∞ = +2u m to the matching point (origin u = 0), with
nagative step height.

2.2. Matching Criteria

The acceptable solution of the Schrodinger equation, the wave function and its first derivative must be continuous
at all the points of region [−2u m , 2u m ], so as at the matching point u = 0, for selected ξ. The matching is confidently
done by comparing the logarithmic derivative of both the solutions[4]. If γ = ψ0 /ψ and if x = x c is the matching
point, then matching condition is
   
γl  x = γr  x or γl − γr  < δ (26)
c c

where γl and γr are the logarithmic derivative of solutions starting from left and right regions respectively. Increment
in ξ from the guess value must be small, as an example if |γl −γr | ≤ 0.0001, then required increment is ξ = ξ+0.0001.
In the case of antisymmetric wave function ψ1 = −ε and matching condition at origin x = x c would be
 
1  1 
=
γl  x =xc γr  x =xc

In addition to ψ, its fist derivative is also required at x = x c for matching. One can use the central difference formula
of order O (h 4 ).

− f2 + 8 f1 − 8 f−1 + f−2
f 0 (x0 ) =
 12h   
−ψ x +2h + 8 ψ x +h − 8 ψ x + 8 +ψ x
c −h c −2h
or ψ0 (x c ) = c c
(27)
12h
The integration must be continued up to two more points after the matching boundary. Note that direction in which
derivative is evaluated for both the wave functions from outward and inward integration must be the same.
It is obvious from the figure that width of the potential must cover the bound energy level, therefore xm should be
p
such that V (xm ) ≥ En . Thus, range of the dimensionless independent variable u is [−u m , u m ], where u m = α xm .
For the case n = 0 : V (xm ) > 12 h ν i.e. u m > ξ or u m > ξ
2
p

Fig. 5. Simple harmonic potential wave functions and associated eigenvalues for n = 0 (left) and n = 1 (right).

107
 Bound state eigenfunctions of an anharmonic oscillator in one dimension: A Numerov method approach

3. Eigenfunctions of Anharmonic Oscillator

First two eigenfunctions and eigenvalues for a harmonic oscillator are evaluated by solving Schrödinger equation
with the Numerov method as shown in fig.5. Now the numerical derivation of an anharmonic oscillator is followed
by the potential energy of the form V (x ) = C2 x 2 + D2 x 4 . The Schrodinger equation for the anharmonic potential is
followed from Eq.(22) will be of the form

d 2ψ D
+ (β − α2 x 2 − x 4 ) ψ = 0 (28)
d x2 2
p
Considering dimensionless variable u = αx

d 2ψ D u4
α + (β − αu 2 − )ψ = 0
d u2 2 α2
d 2ψ D
+ (ξ − u 2 − δu 4 )ψ = 0 ; δ =
d u2 2α3
or

d 2ψ
= (u 2 + δu 4 − ξ)ψ. (29)
d u2

3.1. Result and Discussion

The numerical procedure discussed in previous section is tested in sequential executions of corresponding code. It
is obvious that the particular value of δ is constant for the differential equation (29). The independent variable is u
and ξ is the parameter results from guessed energy eigenvalue. The eigenvalue selected through iterative process
which satisfy proper matching of logarithmic derivatives while developing eigenfunction about the origin. Thus
eigenvalues and corresponding eigenfunctions are executed from IDL 0.7 code based on the above Numerov ap-
proach prepared by the author. The graph of eigenfunctions are generated within the code. The plotted graphs has
the same y-axis to compare wave functions and related potential widths. The height of the wave functions is pro-
portionally reduce for perfect matching. The horizontal lines along y-axis represents energy eigenvalues which is
combined together with their respective wave functions.

Table 1. Ground state energies of quartic anharmonic oscillator evaluated by Numerov method and compared with that of
Bhattacharya

δ Numerov Method Ranjan Bhattacharya

0.1 1.0674000 1.06529
0.2 1.1192000 1.11829
1.0 1.3923000 1.3923
4.0 1.9031000 1.9031
10.0 2.4491000 2.4491
50.0 4.0039000 4.00399

V(x)

Fig. 6. Anharmonic potential and wave functions for n = 0 to n = 4 generated by Numerov method with δ = 0.10 .

108
 Balkrishna P. Shah / Int. J. Adv. Appl. Math. and Mech. 2(2) (2014) 102 - 109

The results of ground state energies of anharmonic oscillator are obtained by R. Bhattcharya solving polynomial
equation[5]. By the application of Numerov method the energy eigenvalues evaluated exactly by taking care of energy
convergence in the present work. The graphical representation of anharmonic potential along with corresponding
bound states and their wave functions within a single graph clarify distribution of wave function along the effective
range of potential and its nature. The input parameters have been playing major role in the computational task are
mentioned here.

1. The guessed eigenvalue should be the lower bound value to stop the divergence.

2. The search for eigenvalue is restricted to energy region where the wave function has same number of nodes
about the origin.

3. xm will be different for different states.

4. Accuracy of the solution depends on the step size h and choice of xm . These parameters must be varied till
energy convergence is achieved.

Acknowledgement
The author specially thanks to Dr. Dilip Angoom, Physical Research Laboratory, Ahmedabad for his fruitful guidance
during the present work.

References

[1] D.R.Hartree, The calculation of atomic structures, John Wiley & sons, New york, 1957.
[2] P. C. Chow, Computer solution to Schrödinger equation, Am. J. Phys. 40 (1972) 730-734.
[3] R. M. Eisberg, R. Resnick,Quantum Physics of atoms,molecules, solid, nuclei and particals, Wiley 2nd Edition,
India,1985
[4] D. G. Kanhere, Physics through computation-II, IAPT Physics Education (April 2007) 55-60.
[5] R. Bhattacharya, D. Roy, S. Bhowmick, Simple systematics in the ground state energies of quantum anharmonic
oscillators, Phy. Lett. A. 244 (1998) 9.

109