Simulating Schrödinger Equation in 1 Dimension

The time dependent schrodingers equation is:

iℏ∂ψ(x,t)∂t=−ℏ22m∂2ψ(x,t)∂x2+V(x)ψ(x,t)i \hbar \dfrac{\partial \psi (x, t)}{\partial t} = -\dfrac{\hbar^2}{2m} \dfrac{\partial^2 \psi(x, t)}{\partial x ^2} + V(x)\psi(x, t)

Divide both sides by iℏi\hbar and then we have a way to do it iteratevely

∂ψ(x,t)∂t=−ℏ22m∂2ψ(x,t)∂x2+V(x)ψ(x,t)iℏ\dfrac{\partial \psi (x, t)}{\partial t} = \dfrac{-\dfrac{\hbar^2}{2m} \dfrac{\partial^2 \psi(x, t)}{\partial x ^2} + V(x)\psi(x, t)}{i \hbar}

and 1i=−i\dfrac{1}{i} = -i so

∂ψ(x,t)∂t=iℏ2m∂2ψ(x,t)∂x2+−iℏV(x)ψ(x,t)\dfrac{\partial \psi (x, t)}{\partial t} = i\dfrac{\hbar}{2m} \dfrac{\partial^2 \psi(x, t)}{\partial x ^2} + \dfrac{-i}{\hbar}V(x)\psi(x, t)

If this is implemented as simple as possible with Euler's method, it will usually be inaccurate and cause the wave to blow up. Simulating the equation will require a more accurate method.

First convert the equation to use the Hamiltonian, which is equal to the sum of the kinetic energy and potential energy.

A classical wave will look like

ψ(x,t)=ei(kx−ωt)\psi(x, t) = e^{i(kx - \omega t)}

k=τλk = \dfrac{\tau}{\lambda} is the wave vector, the radians per distance, and λ\lambda is the wave length of the wave.

ω=τv\omega = \tau v is the angular velocity.

By definition of hh, the momentum of a wave is p=hλ=ℏkp = \dfrac{h}{\lambda} = \hbar k

Then:

∂ψ(x,t)∂x=ikei(kx−ωt)=ikψ(x,t)\dfrac{\partial \psi(x, t)}{\partial x} = ike^{i(kx - \omega t)} = ik\psi(x, t)
∂2ψ(x,t)∂x2=i2k2ei(kx−ωt)=−k2ψ(x,t)\dfrac{\partial ^2 \psi(x, t)}{\partial x ^2} = i^2k^2e^{i(kx - \omega t)} = -k^2\psi(x, t)
∂2ψ(x,t)∂x2=−p2ℏ2ψ(x,t)\dfrac{\partial ^2 \psi(x, t)}{\partial x ^2} = -\dfrac{p^2}{\hbar^2}\psi(x, t)
p2ψ(x,t)=−ℏ2∂2ψ(x,t)∂x2p^2 \psi(x, t) = -\hbar^2 \dfrac{\partial ^2 \psi(x, t)}{\partial x ^2}

Calculating kinetic energy is K=12mv2=p22mK = \frac{1}{2}mv^2 = \dfrac{p^2}{2m}, so:

2mK(x,t)ψ(x,t)=−ℏ2∂2ψ(x,t)∂x22mK(x, t)\psi(x, t) = -\hbar^2 \dfrac{\partial ^2 \psi(x, t)}{\partial x ^2}
K(x,t)ψ(x,t)=−ℏ22m∂2ψ(x,t)∂x2K(x, t)\psi(x, t) = -\dfrac{\hbar^2}{2m} \dfrac{\partial ^2 \psi(x, t)}{\partial x ^2}

And the Hamiltonian is H=K+VH = K + V.

H=K+VH = K + V
H(x,t)ψ(x,t)=K(x,t)ψ(x,t)+V(x)ψ(x,t)H(x, t)\psi(x, t) = K(x, t)\psi(x, t) + V(x)\psi(x, t)
H(x,t)ψ(x,t)=−ℏ22m∂2ψ(x,t)∂x2+V(x)ψ(x,t)H(x, t)\psi(x, t) = -\dfrac{\hbar^2}{2m} \dfrac{\partial ^2 \psi(x, t)}{\partial x ^2} + V(x)\psi(x, t)

And that is the same as the original equation so

H(x,t)ψ(x,t)=−ℏ22m∂2ψ(x,t)∂x2+V(x)ψ(x,t)=iℏ∂ψ(x,t)∂tH(x, t)\psi(x, t) = -\dfrac{\hbar^2}{2m} \dfrac{\partial ^2 \psi(x, t)}{\partial x ^2} + V(x)\psi(x, t) = i \hbar \dfrac{\partial \psi (x, t)}{\partial t}
iℏ∂ψ(x,t)∂t=H(x,t)ψ(x,t)i \hbar \dfrac{\partial \psi (x, t)}{\partial t} = H(x, t)\psi(x, t)
∂ψ(x,t)∂t=−iℏH(x,t)ψ(x,t)\dfrac{\partial \psi (x, t)}{\partial t} = -\dfrac{i}{\hbar}H(x, t)\psi(x, t)

And we can replace the Hamiltonian at a point with the Hamiltonian operator H^\hat{H}

∂ψ(x,t)∂t=−iℏH^ψ(x,t)\dfrac{\partial \psi (x, t)}{\partial t} = -\dfrac{i}{\hbar}\hat{H}\psi(x, t)

For the computer simulation we will represent ψ\psi as a vector, so if we can make the H^\hat{H} a matrix then we can approximate the differential equation using matrix exponentiation.

∂ψ(x,t)∂t=−iℏH^ψ(x,t)\dfrac{\partial \psi (x, t)}{\partial t} = -\dfrac{i}{\hbar}\hat{H}\psi(x, t)
ψ(x,t)=e−iℏH^tψ(x,0)\psi (x, t) = e^{-\frac{i}{\hbar}\hat{H}t}\psi(x, 0)

The Hamiltonian operator is equal to H^=−ℏ22m∂2∂x2+V(x)\hat{H} = -\dfrac{\hbar^2}{2m} \dfrac{\partial ^2}{\partial x ^2} + V(x).

The second derivative can be calculated by:

lim⁡h→0df(x+h)dx−df(x)dxh\lim_{h \rightarrow 0}\dfrac{\dfrac{df(x+h)}{dx} - \dfrac{df(x)}{dx}}{h}
lim⁡h→0f(x+h)−f(x)h−f(x)−f(x−h)hh\lim_{h \rightarrow 0}\dfrac{\dfrac{f(x+h)-f(x)}{h} - \dfrac{f(x)-f(x-h)}{h}}{h}
lim⁡h→0(f(x+h)−f(x))−(f(x)−f(x−h))h2\lim_{h \rightarrow 0}\dfrac{(f(x+h)-f(x)) - (f(x)-f(x-h))}{h^2}
lim⁡h→0f(x+h)−2f(x)+f(x−h)h2\lim_{h \rightarrow 0}\dfrac{f(x+h)-2f(x)+f(x-h)}{h^2}

When converted to discrete form:

f(x+Δx)−2f(x)+f(x−Δx)Δx2\dfrac{f(x+\Delta x)-2f(x)+f(x-\Delta x)}{\Delta x^2}

And for each position of −ℏ22m∂2∂x2-\dfrac{\hbar^2}{2m} \dfrac{\partial ^2}{\partial x ^2} we can make a matrix like:

−ℏ22m1Δx2[−211−21⋱1−211−2]-\dfrac{\hbar^2}{2m}\dfrac{1}{\Delta x^2}\begin{bmatrix}-2 & 1 & & & \\1 & -2 & 1 & & \\ & & \ddots & & \\ & & 1 & -2 & 1 \\ & & & 1 & -2\end{bmatrix}

Then for the V(x)V(x) part, it is just a value at each point:

[V(0)V(1)⋱V(n−1)V(n)]\begin{bmatrix}V(0)& & & & \\ &V(1)& & & \\ & & \ddots & & \\ & & &V(n-1)& \\ & & & &V(n)\end{bmatrix}

So H^\hat{H} is just the sum of those two matricies

H^=−ℏ22m1Δx2[−211−21⋱1−211−2]+[V(0)V(1)⋱V(n−1)V(n)]\hat{H} = -\dfrac{\hbar^2}{2m}\dfrac{1}{\Delta x^2}\begin{bmatrix}-2 & 1 & & & \\1 & -2 & 1 & & \\ & & \ddots & & \\ & & 1 & -2 & 1 \\ & & & 1 & -2\end{bmatrix} + \begin{bmatrix}V(0)& & & & \\ &V(1)& & & \\ & & \ddots & & \\ & & &V(n-1)& \\ & & & &V(n)\end{bmatrix}

And then each discrete step of the simulation will go through a certain change in time

ψt+1=e−iℏH^Δtψt\psi_{t+1} = e^{-\frac{i}{\hbar}\hat{H}\Delta t}\psi_t

This method is more accurate but more computationally complex, since the matrix will have to be the same size as the size of vector ψ\psi making this algorithm ∈O(n2)\in O(n^2) of the number of elements in the vector.

I used this link for the explanation of the derivation and this link for explaining how to simulate it using the Hamiltonian.