Schrödinger equation

Go back to Examples.

2D infinite square well

We are solving a Schrödinger equation for a quantum particle that is confined to a 2-dimensional square box, $0<x<L$ and $0<y<L$. Within this square the particle behaves as a free particle, but the walls are impenetrable, so the wave function is zero at the walls. This is described by the infinite potential

\[ V(x)= \begin{cases} 0, &0<x<L,\: 0<y<L, \\ \infty, &{\text{otherwise.}} \end{cases} \]

With this potential Schrödinger equation simplifies to \[ -{\hbar^2 \over 2m} \nabla^2 \Psi(x,y,t) = i\hbar {\partial \over \partial t}\Psi(x,y,t), \] where $m$ is the mass of the particle, $\hbar$ is reduced Planck constant and $\Psi(x,y,t)$ is the wave function. We set $\hbar / 2m = 1$ and $L=1$.

The solution $\Psi(x,y,t)$ must therefore satisfy \[ \nabla^2 \Psi = -i{\partial \over \partial t}\Psi \quad \text{in} \quad \Omega, \] with Dirichlet boundary condition \[ \Psi = 0 \quad \text{on} \quad \partial \Omega, \] where $\Omega=[0, 1]\times[0, 1]$ denotes the square box domain and $\partial \Omega$ is the boundary of the domain. In order to get the unique solution for the time propagation of the wave function, the initial state $\Psi(x, y,t= 0)$ has to be selected. We choose \[ \Psi(x,y,t=0) = \sin{(\pi x)} \sin{(\pi y)}. \]

Once we prepare the domain and operators, we can turn the Schrödinger equation into code. As we are solving the time propagation implicitly, we will have to solve a matrix equation for every step. For this reason a space matrix $M$ is constructed. Next the equations are implemented, $psi$ denotes the matrix solution for each step and $rhs$ the right-hand side of the matrix equation, which is in our case equivalent to the previous state $psi$.

SparseMatrix<complex<double>, RowMajor> M(N, N);

// Set initial state.
for(int i : domain.interior()) {
    double x = domain.pos(i, 0);
    double y = domain.pos(i, 1);
    rhs(i) = sin(PI * x) * sin(PI * y);
}

// Set equation on interior.
for(int i : domain.interior()) {
    op.value(i) + (-1.0_i) * dt * op.lap(i) = rhs(i);
}

// Set equation on boundary.
for (int i: domain.boundary()){
    op.value(i) = 0;
}

At each time step the matrix equation is solved using an Eigen linear equation solver, in this case the BICGStab algorithm, and $rhs$ is updated.

// Solve matrix system.
psi = solver.solve(rhs);

// Update rhs.
rhs = psi;

Time propagation of the real part of the solution $\Psi(x,y,t)$ from times 0 to $T=0.63$ is presented below. Time step used was $10^{-5}$.

As analytical solutions for infinite potential well are known this enables comparison between analytical and numerical solution computed with Medusa. On the graph bellow the error of the numerical solution is presented as maximal absolute error $e_{\infty}$, meaning the largest absolute difference between numeric and its analytical solution in node at specific time. Expected linear convergence of error on a log - log scale can be observed.

The complete exampe is avalible at infinite_well_2D.cpp.

3D infinite cubic well

Now we are dealing with a problem very similar to the previous one, the only difference is, that in this case we have a quantum particle in a 3-dimensional box. The equation and the boundary condition

\[ \nabla^2 \Psi = -i{\partial \over \partial t}\Psi \quad \text{in} \quad \Omega, \\ \Psi = 0 \quad \text{on} \quad \partial \Omega, \]

remain unchanged, merely the domain changes to unit cube, $\Omega=[0, 1]\times[0, 1]\times [0, 1]$. The only step left is to select an initial state of the wave function, we choose \[ \Psi(x,y,t=0) = \sin{(\pi x)} \sin{(\pi y) \sin(\pi z)}. \] The plot of the real and imaginary part of the solution $\Psi(x, y, z)$ at time $T=0.1$ is shown bellow.

The complete exampe is avalible at infinite_well_3D.cpp.

Go back to Examples.