函数在边界处发散：薛定谔 2D，显式方法

Question

我正在尝试使用the explicit algorithm proposed by Askar and Cakmak (1977) 模拟二维薛定谔方程。我用复杂函数 u+iv 定义了一个 100x100 网格，边界处为空。问题是，仅仅经过几次迭代，复函数的绝对值就会在边界附近爆炸。

我把代码贴在这里，有兴趣的可以去看看：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm                     
from mpl_toolkits.mplot3d import Axes3D

#Initialization+meshgrid

Ntsteps=30
dx=0.1
dt=0.005
alpha=dt/(2*dx**2)

x=np.arange(0,10,dx)
y=np.arange(0,10,dx)
X,Y=np.meshgrid(x,y)          

#Initial Gaussian wavepacket centered in (5,5)

vargaussx=1.
vargaussy=1.
kx=10
ky=10

upre=np.zeros((100,100))
ucopy=np.zeros((100,100))
u=(np.exp(-(X-5)**2/(2*vargaussx**2)-(Y-5)**2/(2*vargaussy**2))/(2*np.pi*(vargaussx*vargaussy)**2))*np.cos(kx*X+ky*Y)
vpre=np.zeros((100,100))
vcopy=np.zeros((100,100))
v=(np.exp(-(X-5)**2/(2*vargaussx**2)-(Y-5)**2/(2*vargaussy**2))/(2*np.pi*(vargaussx*vargaussy)**2))*np.sin(kx*X+ky*Y)

#For the simple scenario, null potential

V=np.zeros((100,100))

#Boundary conditions

u[0,:]=0
u[:,0]=0
u[99,:]=0
u[:,99]=0
v[0,:]=0
v[:,0]=0
v[99,:]=0
v[:,99]=0

#Evolution with Askar-Cakmak algorithm

for n in range(1,Ntsteps):

   upre=np.copy(ucopy)
   vpre=np.copy(vcopy)
   ucopy=np.copy(u)     
   vcopy=np.copy(v)

   #For the first iteration, simple Euler method: without this I cannot have the two steps backwards wavefunction at the second iteration
   #I use ucopy to make sure that for example u[i,j] is calculated not using the already modified version of u[i-1,j] and u[i,j-1]

   if(n==1):
       upre=np.copy(ucopy)
       vpre=np.copy(vcopy)

   for i in range(1,len(x)-1):
       for j in range(1,len(y)-1):
           u[i,j]=upre[i,j]+2*((4*alpha+V[i,j]*dt)*vcopy[i,j]-alpha*(vcopy[i+1,j]+vcopy[i-1,j]+vcopy[i,j+1]+vcopy[i,j-1]))
           v[i,j]=vpre[i,j]-2*((4*alpha+V[i,j]*dt)*ucopy[i,j]-alpha*(ucopy[i+1,j]+ucopy[i-1,j]+ucopy[i,j+1]+ucopy[i,j-1]))

#Calculate absolute value and plot

abspsi=np.sqrt(np.square(u)+np.square(v))

fig=plt.figure()
ax=fig.add_subplot(projection='3d')
surf=ax.plot_surface(X,Y,abspsi)        
plt.show()

如您所见，代码非常简单：我看不出这个错误来自哪里（我认为不是稳定性问题，因为 alpha

Answer 1

我会尝试将您的 dt 设置为较小的值（例如 0.001）并增加集成步骤的数量（例如五倍）。 在使用dt=0.001 尝试您的代码时，波函数的形状也在Ntsteps=150 上。

检查运动的积分（例如这里的动能？）还应该确认dt的不同选择情况是否正常。