Python implementation of the Numerov method applied to the one-dimensional quantum harmonic oscillator, solving the time-independent Schrodinger equation.
The companion Java implementation is available at
numerov-schrodinger-java.
The time-independent Schrodinger equation for a particle in a harmonic potential can be rewritten in dimensionless form as a second-order ODE where the second derivative of the wavefunction equals a position-dependent coefficient times the wavefunction itself.
The Numerov method solves this type of equation with sixth-order accuracy by expressing each new wavefunction value in terms of the two previous ones, weighted by the local coefficient evaluated at those points. This makes the method both highly accurate and simple to implement.
Eigenvalues are found by scanning for sign changes of the wavefunction at the right boundary as a function of the trial energy, then halving the search step each time a sign change is detected until the desired tolerance is reached.
| n | Exact energy | Numerical energy | Error (%) |
|---|---|---|---|
| 0 | 0.5 | 0.499999998 | 4.8e-07 |
| 1 | 1.5 | 1.499999998 | 1.6e-07 |
| 5 | 5.5 | 5.499999998 | 4.3e-08 |
| 10 | 10.5 | 10.499999998 | 2.3e-08 |
Errors are 2 to 4 orders of magnitude smaller than those obtained with the original implementation using a narrower grid, without any increase in the number of grid points.
Running the script generates three PNG files in the current directory:
| File | Description |
|---|---|
wavefunctions.png |
Wavefunctions for n = 0 to 5, Numerov vs Analytical |
penetration.png |
Probability of finding the particle in the classically forbidden region |
energy_error.png |
Percentage error per quantum state on a log scale |
numpy>=1.24
matplotlib>=3.7
scipy>=1.11
git clone https://github.com/YuriCastroDev/numerov-schrodinger-python.git
cd numerov-schrodinger-python
pip install -r requirements.txt
python numerov.pyAll parameters are defined at the top of numerov.py:
| Parameter | Default | Description |
|---|---|---|
X_MIN / X_MAX |
-8 / +8 | Grid boundaries |
N_POINTS |
10 000 | Number of grid points |
TOLERANCE |
5e-9 | Eigenvalue search tolerance |
N_MAX |
10 | Highest quantum number computed |
Increasing N_MAX beyond 15 may require wider grid boundaries to maintain accuracy,
since the wavefunctions of higher states extend further from the origin.
Two subtle issues were identified and corrected during development:
Boundary divergence — The Numerov solution grows exponentially near the grid edges due to floating-point contamination from the unphysical growing solution that is always present in numerical integration. Both boundaries are scanned inward until the growth stops, and those regions are zeroed before normalization.
Spurious integration panel — Computing the penetration probability with a single boolean mask over the forbidden region causes the trapezoidal rule to treat the gap between the two forbidden segments as a single wide panel, significantly inflating the result. The fix is to integrate the left and right tails separately and sum them.
| Language | Computation time | Speedup |
|---|---|---|
| Python | ~10.2 s | 1x |
| Java | ~0.17 s | ~60x |
The Numerov loop is inherently sequential since each step depends on the two previous values. This prevents NumPy vectorization and makes the CPython interpreter overhead dominant.
Java implementation: numerov-schrodinger-java