Difference between revisions of "Schrödinger equation"

From Medusa: Coordinate Free Mehless Method implementation
Jump to: navigation, search
Line 1: Line 1:
Solution procedure is still compiling ... so please wait for results :)
+
Go back to [[Medusa#Examples|Examples]].
  
= Introduction =
+
== 2D infinite square well ==
The quantum world is governed by the [https://en.wikipedia.org/wiki/Schr%C3%B6dinger_equation Schrödinger equation]
+
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
  
\[{\displaystyle {\hat {H}}|\psi (t)\rangle =i\hbar {\frac {\partial }{\partial t}}|\psi (t)\rangle } \]
+
\[
 +
V(x)=
 +
\begin{cases}
 +
0, &0<x<L,\: 0<y<L,
 +
\\
 +
\infty, &{\text{otherwise.}}
 +
\end{cases}
 +
\]
  
where $\hat H$ is the [https://en.wikipedia.org/wiki/Hamiltonian_(quantum_mechanics) Hamiltonian], $|\psi (t)\rangle$ is the [https://en.wikipedia.org/wiki/Wave_function quantum state function] and $\hbar$ is the reduced [https://en.wikipedia.org/wiki/Planck_constant Planck constant].
+
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 Hamiltonian consists of kinetic energy $\hat T$ and potential energy $\hat V$. As in classical mechanics, potential energy is a function of time and space, whereas the kinetic energy differs from the classical world and is calculated as
+
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,\, 0)$ has to be selected. We choose
 +
\[
 +
\Psi(x,y;t=0) = \sin{(\pi x)} \sin{(\pi y)}.
 +
\]
  
\[\hat T = - \frac{\hbar^2}{2m} \nabla^2 .\]
+
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$ .  
  
The final version of the single particle Schrödinger equation can be written as
+
<syntaxhighlight lang="c++" line>
 +
Eigen::SparseMatrix<std::complex<double>, Eigen::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.0i) * dt * op.lap(i) = rhs(i);
 +
}
  
\[\left(- \frac{\hbar^2}{2m} \nabla^2 + V(t, \mathbf r)\right) \psi(t, \mathbf r) = i\hbar {\frac {\partial }{\partial t}}\psi(t, \mathbf r) \]
+
// set equation on boundary
 +
for (int i: domain.boundary()){
 +
    op.value(i) = 0;
 +
}
 +
</syntaxhighlight>
  
Quantum state function is a complex function, so it is usually split into the real part and imaginary part
+
At each time step the matrix equation is solved using an Eigen linear equation solver, in this case, the BICGStab algorithm.
  
\[ u, v \in C(\mathbb R)\colon \psi = u + i v , \]
+
<syntaxhighlight lang="c++" line>
 +
// solve matrix system
 +
psi = solver.solve(rhs);
  
which for a real $V$ yields a system of two real equations
+
// update rhs
 +
rhs = psi;
 +
</syntaxhighlight>
  
\[\left(- \frac{\hbar^2}{2m} \nabla^2 + V(t, \mathbf r)\right) u(t, \mathbf r) = -\hbar {\frac {\partial }{\partial t}} v(t, \mathbf r) , \]
+
The animation of the real part of solution $\Re \Psi(x,y; t)$ is presented below.
\[\left(- \frac{\hbar^2}{2m} \nabla^2 + V(t, \mathbf r)\right) v(t, \mathbf r) = \hbar {\frac {\partial }{\partial t}} u(t, \mathbf r) , \]
 
  
which may be easier to handle.
+
[[File:video.avi|400px]]
  
= Particle in a box =
+
Go back to [[Medusa#Examples|Examples]].
 
 
By selecting the potential $V(t, \mathbf r)$ and the initial state $\psi(0, \mathbf r)$ we get a unique solution for time propagation of the quantum state function. A theoretical one dimensional potential
 
 
 
\[\displaystyle V(x)={\begin{cases}0,&0<x<L,\\\infty ,&{\text{otherwise,}}\end{cases}}\]
 
 
 
is known as an infinite potential well. Its time independent eigenfunctions are
 
 
 
\[\sqrt{\frac{2}{L}}\psi_n(x) = \sin\left(k_n x \right), \qquad n = 1,2,3,...\]
 
 
 
where $k_n = \frac{\pi n}{L}$. With a time dependency
 
 
 
\[\psi_n(t, x) = \mathrm e ^ {-i \omega_n t} \psi_n(x),\]
 
 
 
where $\omega_n$ and $k_n$ are connected through dispersion relation through energy $E_n$
 
 
 
\[{\displaystyle E_{n}=\hbar \omega _{n}={\frac {n^{2}\pi ^{2}\hbar ^{2}}{2mL^{2}}}={\frac {\hbar ^{2} k_n^2}{2m}}}.\]
 

Revision as of 11:25, 17 June 2019

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,\, 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 // set initial state
 3 for(int i : domain.interior()) {
 4     double x = domain.pos(i, 0);
 5     double y = domain.pos(i, 1);
 6     rhs(i) = sin(PI * x) * sin(PI * y);
 7 }
 8 
 9 // set equation on interior
10 for(int i : domain.interior()) {
11     op.value(i) + (-1.0i) * dt * op.lap(i) = rhs(i);
12 }
13 
14 // set equation on boundary
15 for (int i: domain.boundary()){
16     op.value(i) = 0;
17 }

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

1 // solve matrix system
2 psi = solver.solve(rhs);
3 
4 // update rhs
5 rhs = psi;

The animation of the real part of solution $\Re \Psi(x,y; t)$ is presented below.

400px

Go back to Examples.