Simulating Schrödinger Equation in 1 Dimension

The time dependent schrodingers equation is:

i ℏ \frac{\partial ψ (x, t)}{\partial t} = - \frac{ℏ^{2}}{2 m} \frac{\partial^{2} ψ (x, t)}{\partial x^{2}} + V (x) ψ (x, t)

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

\frac{\partial ψ (x, t)}{\partial t} = \frac{- \frac{ℏ^{2}}{2 m} \frac{\partial^{2} ψ (x, t)}{\partial x^{2}} + V (x) ψ (x, t)}{i ℏ}

and $\frac{1}{i} = - i$ so

\frac{\partial ψ (x, t)}{\partial t} = i \frac{ℏ}{2 m} \frac{\partial^{2} ψ (x, t)}{\partial x^{2}} + \frac{- i}{ℏ} V (x) ψ (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) = e^{i (k x - ω t)}

$k = \frac{τ}{λ}$ is the wave vector, the radians per distance, and $λ$ is the wave length of the wave.

$ω = τ v$ is the angular velocity.

By definition of $h$ , the momentum of a wave is $p = \frac{h}{λ} = ℏ k$

Then:

\frac{\partial ψ (x, t)}{\partial x} = i k e^{i (k x - ω t)} = i k ψ (x, t)

\frac{\partial^{2} ψ (x, t)}{\partial x^{2}} = i^{2} k^{2} e^{i (k x - ω t)} = - k^{2} ψ (x, t)

\frac{\partial^{2} ψ (x, t)}{\partial x^{2}} = - \frac{p^{2}}{ℏ^{2}} ψ (x, t)

p^{2} ψ (x, t) = - ℏ^{2} \frac{\partial^{2} ψ (x, t)}{\partial x^{2}}

Calculating kinetic energy is $K = \frac{1}{2} m v^{2} = \frac{p^{2}}{2 m}$ , so:

2 m K (x, t) ψ (x, t) = - ℏ^{2} \frac{\partial^{2} ψ (x, t)}{\partial x^{2}}

K (x, t) ψ (x, t) = - \frac{ℏ^{2}}{2 m} \frac{\partial^{2} ψ (x, t)}{\partial x^{2}}

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

H = K + V

H (x, t) ψ (x, t) = K (x, t) ψ (x, t) + V (x) ψ (x, t)

H (x, t) ψ (x, t) = - \frac{ℏ^{2}}{2 m} \frac{\partial^{2} ψ (x, t)}{\partial x^{2}} + V (x) ψ (x, t)

And that is the same as the original equation so

H (x, t) ψ (x, t) = - \frac{ℏ^{2}}{2 m} \frac{\partial^{2} ψ (x, t)}{\partial x^{2}} + V (x) ψ (x, t) = i ℏ \frac{\partial ψ (x, t)}{\partial t}

i ℏ \frac{\partial ψ (x, t)}{\partial t} = H (x, t) ψ (x, t)

\frac{\partial ψ (x, t)}{\partial t} = - \frac{i}{ℏ} H (x, t) ψ (x, t)

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

\frac{\partial ψ (x, t)}{\partial t} = - \frac{i}{ℏ} \hat{H} ψ (x, t)

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

\frac{\partial ψ (x, t)}{\partial t} = - \frac{i}{ℏ} \hat{H} ψ (x, t)

ψ (x, t) = e^{- \frac{i}{ℏ} \hat{H} t} ψ (x, 0)

The Hamiltonian operator is equal to $\hat{H} = - \frac{ℏ^{2}}{2 m} \frac{\partial^{2}}{\partial x^{2}} + V (x)$ .

The second derivative can be calculated by:

\lim_{h \to 0} \frac{\frac{d f (x + h)}{d x} - \frac{d f (x)}{d x}}{h}

\lim_{h \to 0} \frac{\frac{f (x + h) - f (x)}{h} - \frac{f (x) - f (x - h)}{h}}{h}

\lim_{h \to 0} \frac{(f (x + h) - f (x)) - (f (x) - f (x - h))}{h^{2}}

\lim_{h \to 0} \frac{f (x + h) - 2 f (x) + f (x - h)}{h^{2}}

When converted to discrete form:

\frac{f (x + Δ x) - 2 f (x) + f (x - Δ x)}{Δ x^{2}}

And for each position of $- \frac{ℏ^{2}}{2 m} \frac{\partial^{2}}{\partial x^{2}}$ we can make a matrix like:

- \frac{ℏ^{2}}{2 m} \frac{1}{Δ x^{2}} [\begin{array}{ccccc} - 1 & 1 \\ 1 & - 2 & 1 \\ ⋱ \\ 1 & - 2 & 1 \\ 1 & - 1 \end{array}]

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

[\begin{array}{ccccc} V (0) \\ V (1) \\ ⋱ \\ V (n - 1) \\ V (n) \end{array}]

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

\hat{H} = - \frac{ℏ^{2}}{2 m} \frac{1}{Δ x^{2}} [\begin{array}{ccccc} - 1 & 1 \\ 1 & - 2 & 1 \\ ⋱ \\ 1 & - 2 & 1 \\ 1 & - 1 \end{array}] + [\begin{array}{ccccc} V (0) \\ V (1) \\ ⋱ \\ V (n - 1) \\ V (n) \end{array}]

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

ψ_{t + 1} = e^{- \frac{i}{ℏ} \hat{H} Δ t} ψ_{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 $ψ$ making this algorithm $\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.