薛定谔方程随时间不能正确演变？

Question

我正在使用Crank-Nicolson方案在python中编写代码来演化时间相关的薛定谔方程。我不知道如何处理这种潜力所以我环顾四周，找到了this question的方法，我已经从其他几个来源验证了这一点。据他们说，对于谐振子电位，C-N方案给出了

                                      AΨn+1=A∗Ψn     

其中A的主对角线上的元素是dj = 1 + [（iΔt）/（2m（Δx）^ 2）] + [（iΔt（xj）^ 2）/ 4]以及上下对角线上的元素是a =-iΔt/ [4m（Δx）^ 2]

我理解它的方式，我应该以矩阵Ψn的形式给出一个初始条件（我选择了一个相干态），我需要计算矩阵Ψn+ 1，这是时间Δt后的波函数。为了获得给定步骤的Ψn+ 1，我将矩阵A反转并将其与矩阵A *相乘，然后将结果乘以Ψn。然后，得到的矩阵变为Ψn，用于下一步骤。

但是当我这样做时，我得到的动画不正确。波包应该在边界之间振荡，但在我的动画中，它几乎没有从其初始平均值移动。我只是不明白我做错了什么。我对这个问题的理解是错的吗？或者这是我的代码中的缺陷？请帮忙！我在下面发布了我的代码和我的动画here的视频。我很抱歉代码和问题的长度，但它让我发疯，不知道我的错误是什么。

import numpy as np
import matplotlib.pyplot as plt
L = 30.0
x0 = -5.0
sig = 0.5
dx = 0.5
dt = 0.02
k = 1.0
w=2
K=w**2
a=np.power(K,0.25)
xs = np.arange(-L,L,dx)
nn = len(xs)

mu = k*dt/(dx)**2
dd = 1.0+mu
ee = 1.0-mu
ti = 0.0
tf = 100.0
t = ti
V=np.zeros(len(xs))
u=np.zeros(nn,dtype="complex") 
V=K*(xs)**2/2            #harmonic oscillator potential
u=(np.sqrt(a)/1.33)*np.exp(-(a*(xs - x0))**2)+0j    #initial condition for wave function
u[0]=0.0          #boundary condition
u[-1] = 0.0      #boundary condition

A = np.zeros((nn-2,nn-2),dtype="complex")     #define A
for i in range(nn-3):
    A[i,i] = 1+1j*(mu/2+w*dt*xs[i]**2/4)
    A[i,i+1] = -1j*mu/4.
    A[i+1,i] = -1j*mu/4.
A[nn-3,nn-3] = 1+1j*mu/2+1j*dt*xs[nn-3]**2/4 

B = np.zeros((nn-2,nn-2),dtype="complex")    #define A*
for i in range(nn-3):
B[i,i] = 1-1j*mu/2-1j*w*dt*xs[i]**2/4
    B[i,i+1] = 1j*mu/4.
    B[i+1,i] = 1j*mu/4.
    B[nn-3,nn-3] = 1-1j*(mu/2)-1j*dt*xs[nn-3]**2/4

X = np.linalg.inv(A)    #take inverse of A
plt.ion()
l, = plt.plot(xs,np.abs(u),lw=2,color='blue')   #plot initial wave function
T=np.matmul(X,B)                                #multiply A inverse with A*

while t<tf:
    u[1:-1]=np.matmul(T,u[1:-1]) #updating u but leaving the boundary conditions unchanged
    l.set_ydata((abs(u)))              #update plot with new u
    t += dt
    plt.pause(0.00001)

Answer 1

经过大量的修补，它降低了我的步长。这对我有用 - 我缩小了步长，程序也运行了。如果有人遇到与我相同的问题，我建议玩步长。如果其余代码没问题，这是唯一可能的错误区域。