函数在边界发散：Schrödinger 2D，显式方法

Question

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

我post这里的代码是这样的，有兴趣的可以看看：

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<1/2）。你在过去的模拟中有没有遇到过类似的事情？

Answer 1

我会尝试将您的 dt 设置为较小的值（例如 0.001）并增加积分步骤的数量（例如五倍）。当使用 dt=0.001.

尝试您的代码时，波函数的形状也在 Ntsteps=150 和更远的地方

检查运动的积分（例如这里的动能？）还应该确认对于 dt 的不同选择，事情进展顺利（或不顺利）。

函数在边界发散：Schrödinger 2D，显式方法

Function diverging at boundaries: Schrödinger 2D, explicit method

python

waveform

numerical-stability