Solving 1D Wave Function Evolution with Quantum Simulation

[Mephisto]

I am looking for a way to evolve a wave function in one dimension. I tried searching already but couldn't find anything. Basically I have:

V[x] gives the potential at x
Psi0[x] gives initial conditions on wave function for each x

I would like to calculate Psi[t,x], the wave function at any time t given those two arrays. How could I go about doing this?

 

[Ben Niehoff]

You'll need two arrays for Psi, for the real and imaginary parts. Omitting various constants, Schrödinger's equation is

$$i \partial_t \psi(x,t) = - \partial_x^2 \psi(x,t) + V(x) \psi(x,t)$$

A simple way to implement time-evolution, if you don't need to be extremely precise, is to simply discretize Psi and V and rewrite Schrödinger's equation as a finite difference equation. You can approximate as follows:

$$\partial_t \psi(x,t) \approx \frac{\psi(x,t+\Delta t) - \psi(x,t)}{\Delta t}$$

$$\partial_x^2 \psi(x,t) \approx \frac{\psi(x+\Delta x,t) - 2 \psi(x,t) + \psi(x - \Delta x,t)}{\Delta x^2}$$

You'll have to take some care at the endpoints (x min and x max) to keep errors from propagating due to the finite length of your Psi array. That is, you need some way of assuming what the value of Psi is beyond the limits of the array, so that you can take accurate finite differences.

Once you figure that out, merely plug these approximations into Schrödinger's equation, and you'll get an algebraic equation that you can solve for $\psi(x, t+\Delta t)$ in terms of your Psi array at the current instant of time.

Note: If you want to, you can even include a time-varying potential! Then you can simulate, for example, what happens to a particle in a potential when it is hit by a pulse of light. Have fun!

 

[Mephisto]

Thank you very much for that. Surprisingly, it is really hard to find anything about this on the internet!

I implemented the difference equations, and it looks like it could have worked. I don't actually know, because my solution almost always "blows up". I set up potential V=0 with a gaussian distribution initially in the Real part of Psi, and let it evolve. It nicely diffuses, but then starts to oscillate wildly and spikes form inside it, and make it blow up. Its weird. Maybe I got the equations wrong? Or the boundary conditions? Everywhere outside my range I set Psi=0...

Here is my code for the most crucial part. self.at returns real, imaginary parts of Psi at the passed parameter.

        dtdx = self.dt/(self.dx*self.dx)        
        for x in range(0, self.num):
            r1,i1 = self.at(x-1)
            r2,i2 = self.at(x)
            r3,i3 = self.at(x+1)
            
            newpsir[x] = r2 + dtdx*(-i3 + 2.0*i2 - i1) + i2*self.v[x]*self.dt
            newpsii[x] = i2 + dtdx*(r3 - 2.0*r2 + r1) - r2*self.v[x]*self.dt

Maybe I could try re-normalizing my Psi after every update or something? hmm

 

[signerror]

I am too sleepy to take a close look, but can I point out you can simplify by using Python's built in complex number type?

>>> x = 1 + 2j
>>> x**2
(-3+4j)

On another unsolicited note, Paul Falstad wrote some (much more sophisticated) Java applets on variations of this, and they are open source! For instance, he has 1D quantum mechanics applets (both wells, and periodic boundary conditions ("quantum crystal")), and 1D quantum mechanics with time-varying potentials ("quantum transitions"), as well as many 2D and 3D versions:

http://www.falstad.com/mathphysics.html

 

[Mephisto]

thank you for the suggestion signerror!
I implemented it using COMPLEX() in python, and got the same results, but MUCH slower. So I put it back to what it was before, handling imaginary and real separately. So that means the math is right, and its just that my timestep is too large, or the approximation not good enough or something? :(

 

[signerror]

It shouldn't be - there's nothing wrong with the complex type, it's just a struct {double real; double imag}. I just tired a simple example benchmark (square a million complex numbers), and the complex-type version was fully TWICE as fast - user 1.476s, vs 2.930s (uncompiled).

What is COMPLEX()? Is that a class constructor or function? Don't do that - function calls are hugely inefficient (unless they're compiled inline functions, which aren't technically functions anyway). Use complex literals: x = 2+3j, y = x*z, and so forth.