# Solving the time dependent Schrödinger equation with leapfrog integration in 1D

To my frustration I am struggling to implement leapfrog integration for the time dependent Schrödinger equation.

To the best of my knowledge this was first explicitly done in "A fast explicit algorithm for the time‐dependent Schrödinger equation" in 1991 (see https://aip.scitation.org/doi/abs/10.1063/1.168415).

I'm following the treatment in chapter 10 of "Computational Physics" (2nd edition) by Giordano and Nakanishi, which I think lacks some important details.

The idea is essentially to split the Schrödinger equation into real and imaginary parts, and to solve them separately in a staggered fashion, using centered finite differences - solve the real part at $$t$$, then solve the imaginary part at $$t+\Delta t$$, etc.

I wrote a Python script trying to do this for a stationary gaussian wave packet - I expect the packet to gently spread out, but instead it looks like there is something wrong with my code, as the wave function quickly starts looking pathological. See below:  Here is my Python code. Any suggestions as to what the problem might be are welcome.

import math
import numpy as np
import matplotlib.pyplot as plt

from scipy.ndimage import shift

dx = 0.0005                 # discretised spatial step
dt = 5e-7                   # discretised time step

x0 = 0.4                    # center of Gaussian wave packet
sigma = 10.0**(-3/2)        # width of Gaussian wave packet

def psi_0(x):
'''Initial wavefunction, psi_0(x, t=0), is a Gaussian wavepacket'''
psi0 = np.exp(-(x - x0)**2 / sigma**2)
# normalisation
s = math.sqrt((np.sum(psi0.real**2) + np.sum(psi0.imag**2))*dx)
return psi0/s

def potential():
'''For now we return a constant 0 potential'''
return 0

def psi_real_step(psi):
t1 = psi.real
t2 = -dt/(2*dx**2) * (shift(psi.imag, -1, cval=0) - 2.0*psi.imag + shift(psi.imag, 1, cval=0))
t3 = dt*potential()*psi.imag

return t1 + t2 + t3

def psi_imag_step(psi):
t1 = psi.imag
t2 = dt/(2*dx**2) * (shift(psi.real, -1, cval=0) - 2.0*psi.real + shift(psi.real, 1, cval=0))
t3 = -dt*potential()*psi.real

return t1 + t2 + t3

def time_step(psi):
psi_real = psi_real_step(psi)
psi_imag_old = psi.imag
psi = psi_real + psi.imag*1j
psi_imag = psi_imag_step(psi)

return psi_real, psi_imag, psi_imag_old

def main():
x_range = np.arange(0, 1, dx)

psi = psi_0(x_range)
psi_real = psi.real
psi_imag_old = psi.imag
print(f"psi_imag_old: {psi_imag_old}")

plt.plot(x_range, psi_real**2 + psi_imag_old**2, label="t=0")
plt.legend()
plt.show()

# take the initial dt/2 time step for Im(psi(x, t))
psi_imag = psi_imag_step(psi)
psi = psi_real + psi_imag*1j

for t in np.arange(0, 4e-3, dt):
print(psi)
psi_real, psi_imag, psi_imag_old = time_step(psi)
psi = psi_real + psi_imag*1j

plt.plot(x_range, (psi_real**2 + psi_imag*psi_imag_old), label=f"t={t}")
plt.legend()
plt.show()

if __name__ == '__main__':
main()



Edit I hope this addendum will help some other confused person in the future. The accepted answer of user helloworld922 is correct - I was using too large a time step.

As to WHY I was using those specific parameter values - they are the ones specified in the caption of figure 10.17 in Giordano and Nakanishi, which is the figure I was attempting to reproduce. From the text one might reasonably conclude that it was generated using leapfrog integration. In fact, if you read the errata for the book (available here) you will see that the Crank-Nicolson method was used instead, which has a different stability condition.

• I haven't done the analysis for the Schrodinger equation, but typically for wave equations the center difference is unconditionally unstable. Also, with explicit time stepping you are almost always restricted by the size of the timestep ("CFL" condition for hyperbolic equations). Have you tried using a smaller timestep or doing a Von-Neumann stability analysis to figure out the stability conditions? May 28, 2022 at 21:14

I read through the paper you linked and they give the stability condition for this method to be (eq. A6) $$\frac{-2}{\Delta t} \le V \le \frac{2}{\Delta t} - \frac{2}{m \Delta r^2}$$ This has to be true for V at every grid point. Assuming V is strictly non-negative, then we get $$\Delta t \le \frac{2}{\max(V) + \frac{2}{m \Delta r^2}}$$ When I calculate this using your parameters of $$V=0$$ and $$\Delta r = 0.0005$$ I get $$\Delta t \le 2.5 \cdot 10^{-7}$$, which is half of what you have in your code snippet.