如果你的离散方程正确地代表了你的问题，你应该有正确的特征值，尽管在数值上你可以在解决特征值问题时得到一个小的虚部。

例如，让我们考虑圆柱井的轴对称情况。薛定谔方程可以写成

$$\frac{\partial^2 \psi}{\partial r^2} + \frac{1}{r}\frac{\partial \psi}{\partial r} = E \psi\, .$$

使用朴素的有限差分，我们可以使用以下近似值

$$\begin{align} &\frac{\partial^2 \psi(r_i)}{\partial r^2} \approx \frac{\psi_{i+1} - 2\psi_i + \psi_{i-1}}{\Delta r^2}\, ,\\ &\frac{\partial \psi(r_i)}{\partial r} \approx \frac{\psi_{i+1} - \psi_{i-1}}{2\Delta}\, , \end{align}$$

其中第二个近似导致一个斜对称矩阵。离散问题将是

$$[D_{2r} + R_\text{inv} D_{r}] \{\Psi\} = E \{\Psi\}\, ,$$

存在$D_{2r}$二阶导数矩阵， $D_{r}$一阶导数矩阵和$R_\text{inv}$一个对角矩阵$\frac{1}{r_i + \Delta r/2}$作为条目。

下面的代码片段解决了这个问题

from __future__ import division, print_function
from scipy.special import j0, jn_zeros
import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import eigs
import matplotlib.pyplot as plt

# Parameters
r_max = 1.0
dr = 0.01
r = np.arange(0, r_max, dr)
npts = r.shape[0]
nvals = 10

# Matrices
D2 = diags([1, -2, 1], [-1, 0, 1], shape=(npts - 1, npts - 1))/dr**2
D1 = diags([-0.5, 0, 0.5], [-1, 0, 1], shape=(npts - 1, npts - 1))/dr
r_new = r[0:npts-1] + dr/2
R_inv = diags(1/r_new, 0)
A = D2 + R_inv.dot(D1)

# Finite differences
vals, vecs = eigs(-A, k=nvals, which="SM") 
vecs = np.vstack((vecs, np.zeros(vecs.shape[1])))

# Analytic results
vals_anal = jn_zeros(0, nvals)
vecs_anal = j0(np.outer(r, vals_anal)/r_max)


## Plots

# Eigenvalues
plt.figure()
plt.plot(vals_anal**2)
plt.plot(vals, "ok")
plt.legend(["Analytic eigenvalues",
            "Numeric eigenvalues"])
plt.xlabel(r"$n$")
plt.ylabel(r"$E_n$")


# Eigenvectors
plt.figure()
plt.subplot(121)
plt.plot(r, np.abs(vecs[:, 0:4]/vecs[0, 0:4])**2)
plt.title("Finite differences")
plt.xlabel(r"$r$")
plt.ylabel(r"$|\psi|^2$")

plt.subplot(122)
plt.plot(r, np.abs(vecs_anal[:, 0:4])**2)
plt.title("Analytic solution")
plt.xlabel(r"$r$")
plt.ylabel(r"$|\psi|^2$")

plt.tight_layout()
plt.show()

具有以下特征值

和这些特征向量