Matrix Numerov Method

Matrix Numerov improves the centered second derivative by relating neighboring second derivatives to neighboring wavefunction values. On a uniform grid with smooth coefficients and consistent boundaries, it gives a compact fourth-order eigenvalue discretization. The familiar sixth-order local Numerov identity should not be confused with sixth-order global eigenvalue accuracy.

Required background. Finite differences derives the second-order kinetic matrix; boundary conditions fixes the endpoint problem.

The compact derivative identity

Taylor expansion around $x_j$ yields

$$y_{j-1}-2y_j+y_{j+1} ={}\frac{h^2}{12} (y''_{j-1}+10y''_j+y''_{j+1}) +O(h^6).$$

Let $A$ be the usual second-difference matrix divided by $h^2$ on the interior Dirichlet grid, and let

$$B=I+\frac{h^2}{12}A.$$

Then $B$ has diagonal $10/12$ and adjacent diagonals $1/12$. Dividing the identity by $h^2$ gives $Ay=By''+O(h^4)$, which explains the fourth-order differential approximation. The argument needs sufficient smoothness and boundary closures consistent with the identity.

For $-\alpha y''+Vy=Ey$, with $V$ the diagonal matrix of sampled potential values,

$$-\alpha Ay+BVy=EBy.$$

The endpoints are zero and the potential is assumed finite there, so endpoint contributions from $(V-E)y$ vanish. Singular endpoints require a separate analysis; inserting an infinite potential value at a zero endpoint is not a valid matrix operation.

A Hermitian matrix formulation

Multiplying by $B^{-1}$ gives

$$H_Ny=Ey, \qquad H_N=-\alpha B^{-1}A+V.$$

For this uniform Dirichlet construction, $A$ and $B$ are real symmetric and commute because $B$ is a polynomial in $A$. Thus $B^{-1}A$ is symmetric and $H_N$ is Hermitian. The eigenvalues of $B$ lie between $2/3$ and one, so its inversion is well conditioned.

However, the left matrix $-\alpha A+BV$ in the earlier generalized equation is generally not symmetric: $B$ and a nonconstant diagonal $V$ need not commute. Passing that pair to a generalized Hermitian eigensolver silently violates its hypotheses. Use the Hermitian $H_N$ action or a solver appropriate to the actual matrix pair.

Do not explicitly form a dense inverse for a large problem. Factor the tridiagonal $B$ once and apply $H_Nv$ by computing $Av$, solving $Bz=Av$, and adding $Vv$. The inverse is dense, but its action costs only a tridiagonal solve. If an iterative inner solve is used instead, its tolerance becomes another error source.

Discrete square-well dispersion

For $V=0$, the sine modes diagonalize both $A$ and $B$. With $\theta_n=n\pi/(N+1)$,

$$E_n^{(N)}= \frac{4\alpha\sin^2(\theta_n/2)} {h^2[1-\tfrac13\sin^2(\theta_n/2)]}.$$

Writing $k_n=n\pi/L$ and holding $n$ fixed gives

$$E_n^{(N)}=\alpha k_n^2 -\frac{\alpha k_n^6h^4}{240}+O(h^6).$$

The $h^2$ error cancels. Halving the spacing therefore reduces the leading error by about sixteen in the resolved regime. This provides an independent analytic test of the matrix implementation and its claimed order.

For $L=\alpha=1$ and $h=0.1$, the ground energy is approximately $9.869202265$, compared with $\pi^2\approx9.869604401$. Its error is about $4.02\times10^{-4}$, much smaller than the three-point finite-difference error at the same spacing, but still nonzero. The highest grid modes do not share the low-mode asymptotic accuracy.

Smoothness, interfaces, and radial equations

A discontinuity in $V$ generally invalidates the smooth Taylor expansion across that interface. Splitting the domain and enforcing matching conditions can retain high accuracy within smooth pieces. A Coulomb or centrifugal singularity at the origin likewise needs a boundary treatment informed by the regular radial solution, rather than blindly reusing the nonsingular endpoint stencil.

For radial calculations, apply the method to the reduced wavefunction with its correct norm and near-origin behavior. Vary the first grid location, domain size, and spacing separately. A visually smooth numerical wavefunction does not establish fourth-order convergence; compare successive spacings quantitatively.

The construction here relies on the commuting uniform-grid matrices $A$ and $B$. On a nonuniform grid or with modified endpoint rows, rederive the discrete metric and symmetry instead of assuming the same proof holds. Coordinate transformations offers a systematic route to mapped problems.

Exercises

Why can the Numerov identity have a remainder $O(h^6)$ while the resulting derivative approximation is fourth order?

Solution

The identity relates a second difference to $h^2$ times an average of second derivatives. Dividing by $h^2$ turns its $O(h^6)$ remainder into an $O(h^4)$ derivative error. Global eigenvalue order also requires stability, smoothness, and compatible boundary treatment.

References

• LeVeque, R. J., Finite Difference Methods for Ordinary and Partial Differential Equations, SIAM, 2007 — consistency and boundary-dependent convergence analysis.
• Pillai, M., J. Goglio, and T. G. Walker, “Matrix Numerov Method for Solving Schrödinger’s Equation,” American Journal of Physics 80, 1017–1019 (2012), author copy — matrix kinetic representation and numerical examples.