Animating wave packet with Python

[Avatrin]

Homework Statement

So, I start off with initial condition:

[tex]\psi(x,0) = e^{\frac{-(x-x_0)^2}{4a^2}}e^{ilx}[/tex]

This wavepacket is going to move from [itex]x_0 = -5pm [/itex] towards a potential barrier which is 1 MeV from x = 0 to x = 0.25 pm, and 0 everywhere else.

a = 1 pm
l = 2 pm^-1

Homework Equations

I can use any numerical method to solve this. So, Euler's method, Runge-Kutta methods or Crank-Nicolson's method are all methods I can use.

Also, we have Schrodingers equation:
[tex](-\frac{\hbar^2}{2m}\frac{d^2}{dx^2} + V(x))\Psi(x,t) = i\hbar\frac{d}{dt}\Psi(x,t)[/tex]

The Attempt at a Solution

I tried Euler's method:

[tex]\psi(x,t + \delta t) = \psi(x,t) + \frac{-i \delta t}{\hbar}\hat{H}\psi(x,t)[/tex]
Where:
[tex]\hat{H} = -\frac{\hbar^2}{2m}\frac{d^2}{dx^2} + V(x)[/tex]

I used vectors (linspace). So, V, the second derivative of Psi and Psi were all vectors. However, Python kept casting away the imaginary parts off the elements in the vectors. So, my Psi's end up being constant for all t:

[tex]\psi(x,t + \delta t) = \psi(x,t)[/tex]

I don't know how to get past this.

 

[soarce]

You need to split the equation into real and imaginary part and you will end up with two (real-number) equations.

 

[Avatrin]

But, the real and imaginaries change each other. The loop I am using is:
for i in range(10000):

psi = psi + 1j*dt*k*d2p(psi) - 1j*dt*V*psi
plot(x, psi, axis=[x[0], x[-1], 0, 10],
xlabel='x', ylabel='f', legend='t=%4.2f' % t,
savefig='tmp%04d.png' % i)​

psi starts off being the initial condition [itex] \psi(x,0)[/itex], and changes in the loop. D2p gives me the second derivative of the array I put into it. V(x) is the potential. k is just a real constant.

Even if I split this into two loops by first saving the values for psi them in a different array, and then plotting the element of that array:

for i in range(10000):

psi = psi + 1j*dt*k*d2p(psi) - 1j*dt*V*psi
y.append(psi)​

for i in y:

plot(x, i, axis=[x[0], x[-1], 0, 10],
xlabel='x', ylabel='f', legend='t=%4.2f' % t,
savefig='tmp%04d.png' % i)​

Even here, it does not work. I get a wave that does not move.

 

[soarce]

I mean to have two vectors, something like:
$$ \psi_r = \psi_r + {\rm Re} \left\{ i*dt*k*d2p(\psi) - i*dt*V*\psi \right\} $$
$$ \psi_i = \psi_i + {\rm Im} \left\{i*dt*k*d2p(\psi) - i*dt*V*\psi \right\} $$

unless Python has built-in complex number arithmetic.

LE: Also, check that you obtained correctly the dimensionless equation. If the dimensionless equation is not correct it may happen that the initial velocity which you put is too small to see any visible change in reasonable time window.

 

[Avatrin]

Well, I know longer get that error (I did not need to split it into two vectors). However, it still does not work. Here's my code. Hopefully somebody can tell me what's wrong (Some values suddenly turn into "nan+nanj".. Not a number? Why?):

from numpy import *
from scitools.all import *
import sys

def Psi0(x):
x0 = -5.0 #pm
a = 1.0 #pm
l = 2.0 #pm**-1

A = (1/(4*a**2))**(1/4.0)
K1 = exp(-(x-x0)**2/(4*a**2))
K2 = exp(1j*l*x)
return A*K1*K2

m = 0.511*10**6
hbar = 6.582*10**(-22)

x = linspace(-8.0,2.0,1001).astype(complex128)
t = linspace(0.0,0.01,1001).astype(complex128)

dx = x[2]-x[1]
dt = t[2]-t[1]def V(x,L=0.25):
V = zeros(len(x)).astype(complex128)
for i in range(len(x)):
if x >= 0 and x <= L:
V = 1.0
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)/dx**2

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

k1 = (dt*1.j*hbar/(2.0*m))
k2 = dt*1.j/hbar
V = V(x)
a = Psi[-1] + k1*d2(Psi[-1],dx) - k2*V*Psi[-1]
#plot(x.real,a.real)

for i in range(len(t)):
Psi.append(Psi[-1] + k1*d2(Psi[-1],dx) - k2*V*Psi[-1])

for i in range(len(Psi)):
print Psi[500]"""
for i in counter:
y2 = abs(Psi)**2

plot(x.real, y.real, axis=[-8, 2, -1, 1], xlabel='x', ylabel='f', legend='t=%4.2f' % i, savefig='tmp%04d.png' % i)

movie('tmp*.png')
"""

 

[BvU]

Don't speak Python. Perhaps you can learn something useful looking how these guys produced this movie

 

[soarce]

First of all I would get rid of the dimensional equatian by changing the spatial variable. For instance, divide the entire equation by ##\hbar## and then introduce the new spatial variable ##x_1 = x/2a## and then rewrite the time dependent equation as

[tex]\left(-\frac{\hbar^2}{8ma^2}\frac{d^2}{dx_1^2} + V(x_1)\right)\Psi(x_1,t)=i\hbar\frac{d}{dt}\Psi(x_1,t)[/tex]

as one can see the factor multiplying the second derivative the the dimension of Joule, i.e. energy. One must rewrite also the initial value ##\Psi(x,0)## in terms of ##x_1##. Then you divide the whole equation by ##\hbar^2/8ma^2##:

$$
\left[ -\frac{d^2}{dx_1^2} + V_1(x_1) \right] \Psi(x_1,t) = i \frac{8ma^2}{\hbar}\frac{d}{dt}\Psi(x_1,t) \equiv i \frac{d}{d\tau}\Psi(x_1,\tau)
$$
The last change of variable regard the time t, and we introduced ##\tau=\hbar t /8ma^2##. Pay attention to ##V_1## because it is also expressend in terms of ##\hbar^2/8ma^2##, i.e. is normalized (or scaled).

Although this procedure may look cumbersome and useless is in fact very useful:
(i) We remove as much as we can the dependence on the physical parametrs, in this way instead of varying, let's say, two physical parameters we discover that it is enough to vary their ratio.
(ii) handling dimensionless parameters you don't have to worry about the units and their multiples and submutliples. If you need to recover dimensional parameters you can do it at the end of the simulation process by reversing the changes of variables.
(iii) the dimensional parameters appearing in fron tof the spatial and temporal derivative helps you establish a quick correspondence between your simualtion parameters, spatial/temporal range and step, with the simulation of other person and moreover helps you decide whether you are in a regime where the numerical method is unstable. Usually the analysis of stable-unstable regime is done using dimensionless parameters.
(iv) removing the dependence on the physical parameters we are able to debugg much easier the simulation program.Now, back to simualtion program. If Python support complex number arithmetic, including vetors and arrays, is fine. x and t doesn't have to be complex numbers, only the wavefunction needs to be defined as complex number. Inside the for loop shouldn't be i-1 the index of previous step?

 

[Avatrin]

Now, back to simualtion program. If Python support complex number arithmetic, including vetors and arrays, is fine. x and t doesn't have to be complex numbers, only the wavefunction needs to be defined as complex number. Inside the for loop shouldn't be i-1 the index of previous step?

Python does not support complex numbers, but I am using a module, Numpy, that allows me to use complex numbers (that is also the module that gives me vectors that I can multiply element-wise).
The way I have written the program, I start off with an empty list. I add the array for [itex]\psi(x,0)[/itex] as an element to the list (which is going to contain arrays for each [itex]\psi(x,t_i)[/itex]. I did this to stop Python from casting away all the imaginary values). That is why I use Psi[-1], it being the last array, to create the next element/array which I am going to append to the array of arrays.

I am now going to try your suggestion. It may help. I noticed that the kinetic energy operator added very little to the sum while the potential energy operator added a lot. Thus, the parts of the function that were inside the potential barrier started dominating. So, maybe using a dimensionless version of the differential equation will help.