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.