我试图使问题尽可能详细。我有一个为薛定谔方程编写的非常简单的求解器，但具有虚时间，它基本上将其转换为扩散方程（带有潜在项）。该方法在此页面上有很好的记录，我基本上几乎完全按照这些步骤进行操作。这是链接：
http://jkwiens.com/2010/01/02/finite-difference-heat-equation-using-numpy/
以及相应的 Python 代码：
http://jkwiens.com/heat-equation-using-finite-difference/
这是我的代码：
import scipy as sp
from scipy import integrate, sparse, linalg
import scipy.sparse.linalg
import pylab as pl
nx = 8000
dx = 0.0025
dt = 0.00002
niter = 20
nonlin = 0.0
gridx = sp.zeros(nx)
igridx = sp.array(range(nx))
psi = sp.zeros(nx)
pot = sp.zeros(nx)
depth = 0.01
# Set up grid, potential, and initial state
gridx = dx*(igridx - nx/2)
pot = depth*gridx*gridx
psi = sp.pi**(-1/4)*sp.exp(-0.5*gridx*gridx)
# Normalize Psi
#psi /= sp.integrate.simps(psi*psi, dx=dx)
# Plot parameters
xlimit = [gridx[0], gridx[-1]]
ylimit = [0, 2*psi[nx/2]]
# Set up diagonal coefficients
Adiag = sp.empty(nx)
Asup = sp.empty(nx)
Asub = sp.empty(nx)
bdiag = sp.empty(nx)
bsup = sp.empty(nx)
bsub = sp.empty(nx)
Adiag.fill(1 - dt/dx**2)
Asup.fill(dt/(2*dx**2))
Asub.fill(dt/(2*dx**2))
bdiag.fill(1 + dt/dx**2)
bsup.fill(-dt/(2*dx**2))
bsub.fill(-dt/(2*dx**2))
# Construct tridiagonal matrix
A = sp.sparse.spdiags([Adiag, Asup, Asub], [0, 1, -1], nx, nx)
b = sp.sparse.spdiags([bdiag, bsup, bsub], [0, 1, -1], nx, nx)
# Loop through time
for t in range(0, niter) :
# Calculate effect of potential and nonlinearity
psi *= sp.exp(-dt*(pot + nonlin*psi*psi))
# Calculate spacial derivatives
psi = sp.sparse.linalg.bicg(A, b*psi)[0]
# Normalize Psi
psi /= sp.integrate.simps(psi*psi, dx=dx)
# Output figures
pl.plot(gridx, psi)
pl.plot(gridx, psi*psi)
pl.plot(gridx, pot)
pl.xlim(xlimit)
pl.ylim(ylimit)
pl.savefig('outputla/fig' + str(t))
pl.clf()
该nonlin变量只是非线性的一个额外术语（稍后我将求解非线性薛定谔方程），但如果我将其设置为零，那么它就无关紧要了。
我现在的问题是空间运算符没有做任何事情（它应该扩散得相当快），最终整个事情都爆炸了（典型的不稳定性......线条变得锯齿状，最终情节看起来很疯狂）。
我已经检查了我能想到的一切，我完全不知道为什么会这样。我已经改变了线性求解器，密集矩阵而不是稀疏矩阵，系数A和b是正确的（我认为..它们与博客链接不同，因为 Schro. eq. 有一个i术语，所以会发生一些负正号交换......这个在下面引用的论文中提到）。
另外，这是一篇描述计算的论文（但我不使用递归关系，只是一个运算符）：http ://arxiv.org/abs/0904.3131