# How can I evaluate more accurate energy eigenvalues from Schrodinger equation using shooting method?

I am trying to use the "shooting method" for solving Schrodinger's equation for a reasonably arbitrary potential in 1D. But the eigenvalues so evaluated in the case of potentials that do not have hard boundaries are not much accurate compared to analytical results.

I thought that the problem could be solved, by making the spatial grid finer, but changing the spatial grid does not practically have any effect on the eigenvalues. I am not making the energy-grid finer, because the job of fining down to the correct eigenvalue is tackled by the bisection method from SciPy, and the wavefunction is evaluated by solving the concerned IVP by odeint from SciPy, these functions are accurate enough.

Finally, changing the 2nd boundary to make the wavefunction die out at deeper part of classically forbidden region also did not bring a practical improvement in the eigenvalue (Changes found only in 9th or 10th place of decimal but made wavefunctions of lower energy state divergent at endpoints to make things worse).

I cannot find what to modify to obtain more accurate eigenvalues. Boundary condition or stepsize? Did my implementation go wrong, or is it due to rounding errors or other "Python things"?

## Example: Morse potential

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import bisect

def V(x, xe=1.0, lam=6.0):
"""Morse potential definition"""
return lam**2*(np.exp(-2*(x- xe)) - 2*np.exp(-(x - xe)))

def func(y, x):
"""
Utility function for returning RHS of the differential equation.
"""
psi, phi = y # psi=eigenfunction, phi=spatial derivative of psi
return np.array([phi, -(E - V(x))*psi])

def ivp(f, initial1, initial2, X):
"""Solve an ivp with odeint"""
y0 = np.array([initial1, initial2])
return odeint(f, y0, X)[:, 0]

def psiboundval(E1):
"""
Find out value of eigenfunction at bound2 for energy E1
by solving ivp.
"""
global E;
E = E1
S = ivp(func, bval1, E1, X)
return S[(len(S)) - 1] - bval2

def shoot(Erange):
"""
Find out accurate eigenvalues from approximate ones in
Erange by bisect.
"""
global E
Y = np.array([psiboundval(E) for E in Erange])
eigval = np.array([bisect(psiboundval, Erange[i], Erange[i + 1])
for i in np.where(np.diff(np.signbit(Y)))[0]])
return eigval

#%% Solution
xe, lam = 1.0, 6.0 # parameters for potential
# Bval, Bval2 = wavefunction values at x = bound1, bound2
bound1, bound2, bval1, bval2 = 0, xe + 15, 0, 0
X = np.linspace(bound1, bound2, 1000) # region of integration
Erange = np.geomspace(-lam**2, -0.0001, 100) # region of Energy level searching
print("Numerical results:", np.round(shoot(Erange), 4))
print("Analytical results:",
[-(lam - n - 0.5)**2 for n in range(0, int(np.floor(lam - 0.5) + 1))])


Output

Numerical results: [-30.2483 -20.2432 -12.2361  -6.2318  -2.2343  -0.2438]
Analytical results: [-30.25, -20.25, -12.25, -6.25, -2.25, -0.25]


For higher energy states, accuracy is seen to decrease.It is desirable that accuracy is of at least up to 4th decimal place (if not more), for all states.

• First, the shooting method is a rather expensive approach to solving this problem. Just solve it as an eigenvalue problem to a two-point boundary value problem, rather than trying to treat it as an ODE. May 7 '20 at 14:08
• Second, I don't think I quite understand your question. Are you saying that if you have hard boundaries (infinite potential) then your solution is highly accurate whereas if you have a finite potential you do not? Can you state in mathematical terms how you truncate the domain if you have a finite potential? May 7 '20 at 14:09
• For the second,what you understand is exactly correct.There is no "mathematics" behind the truncation of domain with finite potential.Infact,I guessed i will have to extend the domain to get more accuracy for higher energy terms,so I extended the domain just in hope I will get better eigenvalue,but unfortunately nothing worked,and that is the problem exactly,and the extension is done purely based on guess. The truncation is exactly what will be needed for optimal point,I thought. May 7 '20 at 15:16
• There are several questions in this site related to this topic as well. For example: scicomp.stackexchange.com/a/19679/9667 May 7 '20 at 23:15
• My conclusion after several years of numerical work on the Schrödinger equation is: Don't use the shooting method, and even more so if you're out for accuracy. The shooting method might have its application for non-linear problems or exotic boundary conditions, but certainly not for the 1D linear Schrödinger equation. Use a diagonalization procedure instead. It will give you ten digits more accuracy than in the example in Nicoguaro's answer. And if that's not enough (it usually is), you can improve on these results. May 9 '20 at 16:17

Your problem was the lower limit of integration. It should have been $$-x_e$$ instead of 0, since $$x_e$$ is the equilibrium point for the potential and not the minimum distance.

After correcting that, you get the following

#%% Solution
xe, lam = 1.0, 6.0 # parameters for potential
xmax = 10
# Bval, Bval2 = wavefunction values at x = bound1, bound2
bound1, bound2, bval1, bval2 = -xe, xmax, 0, 0
X = np.linspace(bound1, bound2, xmax) # region of integration
Erange = np.geomspace(-lam**2, -0.01, 100) # region of Energy level searching
print("Numerical results:", np.round(shoot(Erange)[:6], 6))
print("Analytical results:",
[-(lam - n - 0.5)**2 for n in range(6)])


with results

Numerical results: [-30.25     -20.25     -12.25      -6.25      -2.25      -0.240849]
Analytical results: [-30.25, -20.25, -12.25, -6.25, -2.25, -0.25]


I also tried the method for a simple harmonic oscillator and returned the expected eigenvalues.

As mentioned before, I don't think that this method is the best one for the task. You should try a domain discretization method, such as the Finite Element Method or Finite Differences or a variational method. The latter is commonly used in quantum chemistry codes with Gaussian bases.

The two main drawbacks that I see for this method is that you should have a (finely enough discretized) range for your eigenvalues, that might not be known beforehand. Also, I don't see how to generalize the method to higher-dimensional problems.

A scenario where I think it might be useful is when doing a perturbation analysis. In that case a set of (approximated) eigenvalues is available and the same goes for the eigenfunctions.

• Comments are not for extended discussion; this conversation has been moved to chat. Jul 26 at 18:37