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

Question

I'm trying to simulate the 2D Schrödinger equation using the explicit algorithm proposed by Askar and Cakmak (1977). I define a 100x100 grid with a complex function u+iv, null at the boundaries. The problem is, after just a few iterations the absolute value of the complex function explodes near the boundaries.

I post here the code so, if interested, you can check it:

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()

As you can see the code is extremely simple: I cannot see where this error is coming from (I don't think is a stability problem because alpha<1/2). Have you ever encountered anything similar in your past simulations?

Davide Fiocco · Accepted Answer · 2022-01-01T18:58:38.373

1

I'd try setting your dt to a smaller value (e.g. 0.001) and increase the number of integration steps (e.g fivefold). The wavefunction looks in shape also at Ntsteps=150 and well beyond when trying out your code with dt=0.001.

Checking integrals of the motion (e.g. kinetic energy here?) should also confirm that things are going OK (or not) for different choices of dt.

edited Jan 01 '22 at 18:58

answered Jan 01 '22 at 16:53

Davide Fiocco

5,350
5
35
72

I can't believe it was this kind of problem! I thought that with alfa=1/4 I should be safe: at the same time the fact that the explosion was at the boundaries made me suspect was another kind of misterious problem. I'm very grateful for your help Davide, you really made my day :) :) – Rob Tan Jan 01 '22 at 17:40

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

1 Answers1