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 Schrodinger 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 Schrodinger 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$ .

 1 Eigen::SparseMatrix<std::complex<double>, Eigen::RowMajor> M(N, N);
 2 
 3 // set initial state
 4 for(int i : domain.interior()) {
 5     double x = domain.pos(i, 0);
 6     double y = domain.pos(i, 1);
 7     rhs(i) = sin(PI * x) * sin(PI * y);
 8 }
 9 
10 // set equation on interior
11 for(int i : domain.interior()) {
12     op.value(i) + (-1.0i) * dt * op.lap(i) = rhs(i);
13 }
14 
15 // set equation on boundary
16 for (int i: domain.boundary()){
17     op.value(i) = 0;
18 }

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.

1 // solve matrix system
2 psi = solver.solve(rhs);
3 
4 // update rhs
5 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, 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.

Go back to Examples.