Python, numpy and complex functions (PDE's)

Question

Update 4

I have almost given up on getting this right. This is the solution to the time-independent Schrodinger's equation, so the analytical solution is: $\psi(x,t) = \psi(x,0)e^{\frac{-iE t}{\hbar}}$. That means $|\psi(x,t)|^2 = |\psi(x,0)|^2$ for all t. However, I cannot make that work numerically, and I do not know why.

My latest failed attempt:

from numpy import *
from scitools.all import *

V0= 10.0 #eV
a= 1.0 #nm
c = 299792.0 # nm per ps
m= (0.510*10**6)/(c**2) #eV/c**2
hbar= 6.58*10**(-4) #eV ps

z = 1.456 # fra c
alpha = 16.2 # fra c

def Psi0(x):
    psi = x.astype(complex128)
    ka = (z/a)*sqrt((alpha/z)**2 -1.)
    A = exp(ka*a)*cos(z)/(a + 1./ka)
    B = 1./(a+ 1./ka)
    for i in range(len(x)):
        if x[i]<-a or x[i]>a:
            psi[i] = A*exp(-ka*abs(x[i]))
        else:
            psi[i] = B*cos((z/a)*x[i])
    return psi


x = linspace(-8.0,8.0,1001).astype(complex128)
t = linspace(0.0,0.1,2001).astype(complex128)

dx = x[2]-x[1]
dt = t[2]-t[1]
k1 = 1.j*hbar/(2.0*m)
k2 = 1.j/hbar

def V(x,a=1.0):
    V = zeros(len(x)).astype(complex128)
    for i in range(len(x)):
        if x[i] >= -a and x[i] <= a:
            V[i] = -V0
    return V

def d2(Psi, dx):
    D2P = zeros(len(Psi)).astype(complex128)
    for i in range(len(D2P)-2):
        D2P[i+1] = (Psi[i+2] - 2*Psi[i+1] + Psi[i])/(dx**2)
    return D2P

def sch(Psi,x,dx):
    return k1*d2(Psi,dx) - k2*V(x)*Psi

Psi = []
Psi.append(Psi0(x).astype(complex128))

a = abs(Psi[0])**2
plot(x.real, a, xlabel='x', ylabel='f')

''' Euler '''
'''for i in range(len(t)):
    b = dx*sum(abs(Psi[-1] + dt*sch(Psi[-1],x,dx))**2)
    Psi.append((Psi[-1] + dt*sch(Psi[-1],x,dx))/sqrt(b))
    print b
'''
'midpoint'
for i in range(len(t)):
    a = Psi[-1] + dt*sch(Psi[-1] + (dt/2.)*sch(Psi[-1],x,dx),x,dx)
    b = dx*sum(abs(a)**2)
    Psi.append(a/sqrt(b))
    print b

counter = 0
forhver = 10

for i in range(len(t)):
    if counter == forhver:
        a = abs(Psi[i])**2
        plot(x.real, a, axis=[-8, 2, -1, 3], xlabel='x', ylabel='f', legend='t=%4.2f' % i, savefig='obl8%04d.png' % i)
        counter = 0
    counter += 1

movie('obl8*.png')

Update 3:

I changed the units from meter per second to nanometers per picoseconds. It made the code work well. However, it just does not look right: http://sineofmadness.com/wp-content/uploads/2014/11/movie.gif

No matter which $\psi(x,0)$ I choose, and I can choose between five, it turns into that double bellcurve. I have used the midpoint method, and Euler's method. It gives more or less the same result. The animation above screams numerical error to me. Yet, I have now seen that result with several different methods.

Does this look right to anybody?

Original

I am struggling with a specific type of problem in quantum physics (one dimensional): Let's say I have a function $\psi(x,0)$, and I have to show how it evolves over time using the Schrodinger equation (by animating $|\psi(x,t)|^2$): $$ -\frac{\hbar^2}{2m}\frac{d^2\psi}{dx^2} + V(x)\psi = i \hbar\frac{d\psi}{dt}$$ So, what I did was to define functions for the potential and the second derivative, and use Euler's method. At first, the plot function kept throwing away the imaginary parts of the function, so I decided to store the values for $\psi(x,t)$ in a separate array. However, it still does not work. I get problems like every element in the arrays being Not a Number, or the functions remains constant(i.e. $\psi(x,t) = \psi(x,0)$ for all t after some $t_0$).

Answer 1

Caveat: I did not read beyond the statement

So, what I did was to define functions for the potential and the second derivative, and use Euler's method.

So there may be other issues with your code, but this one is already fundamental.

You are trying to solve a wave equation with Euler's method. This is numerically unstable. To see why, note that the second-derivative operator (and its centered finite-difference discretization) is symmetric and has only negative real eigenvalues. When you combine this with the imaginary unit multiplying the time-derivative, you find that the dynamics of this equation are wavelike -- the eigenvalues of the semi-discretization are purely imaginary. Meanwhile, the stability region of Euler's method:

$$\{z\in{\mathbb C : |1+z|\le 1}$$

contains none of the imaginary axis (except for the origin).

Instead of Euler's method, you should use a time discretization that is stable on the imaginary axis, like the midpoint method or the trapezoidal method. Some good methods and references can be found in answers to this question:

Are there simple ways to numerically solve the time-dependent Schödinger equation?

Answer 2

First, you should simplify the problem, so you can test the pieces individually. Get the code working first with V = 0. A simpler Psi0 would also help.

Second, you need to have a better understanding of exactly where things go wrong. Python includes the python debugger pdb. If you put the line: "import pdb;pdb.set_trace()" into your code, it will give you a prompt where you can look at the contents of variables, try out different commands, plot variables, and so on. The pdb documentation can help you with this. I would recommend doing this right before your time loop, to see if everything has been initialized correctly. Then, stick it inside the loop in order to identify exactly where things fail. Debugging one line is much easier than a whole code.

Third, keep at it! There were a bunch of other things I was going to say based on your first code, but the second one solved most of them.