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)))]) 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))])
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.