Medusa: Coordinate Free Mehless Method implementation - User contributions [en]

Complex-valued problems

2019-07-19T14:29:30Z

APirnat:

Go back to [[Medusa#Examples|Examples]].

Our library natively supports complex numbers. To demonstrate this a tutorial for solving Poisson's equation is presented below. For more complicated complex-valued problems see [[Electromagnetic scattering]] and [[Schrödinger equation]].

== 2D Complex Poisson's equation ==
We are solving 2D complex Poisson equation, on a unit square with Dirichlet boundary conditions

\[
\begin{align*}
i \Delta u &= f &&\text{in } \Omega, \\
u &= 0 &&\text{on } \partial \Omega,
\end{align*}
\]
where $u(x,y)$ is the solution to the problem and $\Omega = [0, 1] \times [0, 1]$ the square domain. We will consider $f(x,y) = 2\pi^2\sin(\pi x)\sin(\pi y)$, as it makes for
a simple solution $u(x,y) = i\sin(\pi x)\sin(\pi y)$.

We write the equations for our problem by directly translating the mathematical formulation above into the code.

<syntaxhighlight lang="cpp">
for (int i : domain.interior()) {
double x = domain.pos(i, 0);
double y = domain.pos(i, 1);
1.0_i * op.lap(i) = 2*PI*PI*sin(PI * x)*sin(PI * y);
}
for (int i : domain.boundary()) {
op.value(i) = 0.0;
}
</syntaxhighlight>

The solution $u(x,y)$ is plotted below. In our case solution consists only of the imaginary part, the real part is equal to zero.

[[File:complex_poisson_2D_real.png|500px]]
[[File:complex_poisson_2D_imag.png|500px]]

Go back to [[Medusa#Examples|Examples]].

Complex-valued problems

2019-07-19T14:23:11Z

APirnat:

Go back to [[Medusa#Examples|Examples]].

Our library natively supports complex numbers. To demonstrate this a tutorial for solving Poisson's equation is presented below. For more complicated complex-valued problems see [[Electromagnetic scattering]] and [[Schrödinger equation]].

== 2D Complex Poisson's equation ==
We are solving 2D complex Poisson equation, on a unit square with Dirichlet boundary conditions

\[
\begin{align*}
i \Delta u &= f &&\text{in } \Omega, \\
u &= 0 &&\text{on } \partial \Omega,
\end{align*}
\]
where $u(x,y)$ is the solution to the problem and $\Omega = [0, 1] \times [0, 1]$ denotes the square domain. We will consider $f(x,y) = 2\pi^2\sin(\pi x)\sin(\pi y)$, as it makes for
a simple solution $u(x,y) = i\sin(\pi x)\sin(\pi y)$.

We write the equations for our problem by directly translating the mathematical formulation above into the code.

<syntaxhighlight lang="cpp">
for (int i : domain.interior()) {
double x = domain.pos(i, 0);
double y = domain.pos(i, 1);
1.0_i * op.lap(i) = 2*PI*PI*sin(PI * x)*sin(PI * y);
}
for (int i : domain.boundary()) {
op.value(i) = 0.0;
}
</syntaxhighlight>

The solution $u(x,y)$ is plotted below. In our case solution consists only of the imaginary part, the real part is equal to zero.

[[File:complex_poisson_2D_real.png|500px]]
[[File:complex_poisson_2D_imag.png|500px]]

Go back to [[Medusa#Examples|Examples]].

Schrödinger equation

2019-07-19T14:19:31Z

APirnat:

Go back to [[Medusa#Examples|Examples]].

__FORCETOC__

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

<syntaxhighlight lang="cpp">
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;
}
</syntaxhighlight>

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.

<syntaxhighlight lang="cpp">
// Solve matrix system.
psi = solver.solve(rhs);

// Update rhs.
rhs = psi;
</syntaxhighlight>

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

[[File:video (1).mp4|400px]]

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.

[[File:infinite_well_2D_error.png|400px]]
[[File:infinite_well_2D_error_fit.png|400px]]

The complete exampe is avalible at [https://gitlab.com/e62Lab/medusa/blob/dev/examples/quantum_mechanics/infinite_well_2D.cpp 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.

[[File:infinite_well_3D_rsol.png|400px]]
[[File:infinite_well_3D_csol.png|400px]]

The complete exampe is avalible at [https://gitlab.com/e62Lab/medusa/blob/dev/examples/quantum_mechanics/infinite_well_3D.cpp infinite_well_3D.cpp].

Go back to [[Medusa#Examples|Examples]].

Complex-valued problems

2019-07-19T14:08:16Z

APirnat:

Go back to [[Medusa#Examples|Examples]].

Our library natively supports complex numbers. To demonstrate this a tutorial for solving Poisson's equation is presented below. For more complicated complex-valued problems see [[Electromagnetic scattering]] and [[Schrödinger equation]].

== 2D Complex Poisson's equation ==
We are solving 2D complex Poisson equation, on a unit square with Dirichlet boundary conditions

\[
\begin{align*}
i \Delta u &= f &&\text{in } \Omega, \\
u &= 0 &&\text{on } \partial \Omega,
\end{align*}
\]
where $u(x,y)$ is the solution to the problem and $\Omega = [0, 1] \times [0, 1]$ denotes the square domain. We will consider $f(x,y) = 2\pi^2\sin(\pi x)\sin(\pi y)$, as it makes for
a simple solution $u(x,y) = i\sin(\pi x)\sin(\pi y)$.

We write the equations for our problem by directly translating the mathematical formulation above into code.

<syntaxhighlight lang="cpp">
for (int i : domain.interior()) {
double x = domain.pos(i, 0);
double y = domain.pos(i, 1);
1.0_i * op.lap(i) = 2*PI*PI*sin(PI * x)*sin(PI * y);
}
for (int i : domain.boundary()) {
op.value(i) = 0.0;
}
</syntaxhighlight>

The solution $u(x,y)$ is plotted below. In our case solution consist only of the imaginary part, the real part is equal to zero.

[[File:complex_poisson_2D_real.png|500px]]
[[File:complex_poisson_2D_imag.png|500px]]

Go back to [[Medusa#Examples|Examples]].

Complex-valued problems

2019-07-19T14:07:36Z

APirnat:

Go back to [[Medusa#Examples|Examples]].

Our library natively supports complex numbers. To demonstrate this a tutorial for solving Poisson's equation is presented below. For more complicated complex-valued problems see [[Electromagnetic scattering]] and [[Schrödinger equation]].

== 2D Complex Poisson's equation ==
We are solving 2D complex Poisson equation, on a unit square with Dirichlet boundary conditions

\[
\begin{align*}
i \Delta u &= f &&\text{in } \Omega, \\
u &= 0 &&\text{on } \partial \Omega,
\end{align*}
\]
where $u(x,y)$ is the solution to the problem and $\Omega = [0, 1] \times [0, 1]$ denotes the square domain. We will consider $f(x,y) = 2\pi^2\sin(\pi x)\sin(\pi y)$, as it makes for
a simple solution $u(x,y) = i\sin(\pi x)\sin(\pi y)$.

We write the equations for our problem by directly translating the mathematical formulation above into code.

<syntaxhighlight lang="cpp">
for (int i : domain.interior()) {
double x = domain.pos(i, 0);
double y = domain.pos(i, 1);
1.0_i * op.lap(i) = -2*PI*PI*sin(PI * x)*sin(PI * y);
}
for (int i : domain.boundary()) {
op.value(i) = 0.0;
}
</syntaxhighlight>

The solution $u(x,y)$ is plotted below. In our case solution consist only of the imaginary part, the real part is equal to zero.

[[File:complex_poisson_2D_real.png|500px]]
[[File:complex_poisson_2D_imag.png|500px]]

Go back to [[Medusa#Examples|Examples]].

Complex-valued problems

2019-07-19T14:04:33Z

APirnat:

Go back to [[Medusa#Examples|Examples]].

Our library natively supports complex numbers. To demonstrate this a simple example of Poisson's equation is presented below. For more complicated complex-valued problems see [[Electromagnetic scattering]] and [[Schrödinger equation]].

== 2D Complex Poisson's equation ==
We are solving 2D complex Poisson equation, on a unit square with Dirichlet boundary conditions

\[
\begin{align*}
i \Delta u &= f &&\text{in } \Omega, \\
u &= 0 &&\text{on } \partial \Omega,
\end{align*}
\]
where $u(x,y)$ is the solution to the problem and $\Omega = [0, 1] \times [0, 1]$ denotes the square domain. We will consider $f(x,y) = 2\pi^2\sin(\pi x)\sin(\pi y)$, as it makes for
a simple solution $u(x,y) = i\sin(\pi x)\sin(\pi y)$.

We write the equations for our problem by directly translating the mathematical formulation above into code.

<syntaxhighlight lang="cpp">
for (int i : domain.interior()) {
double x = domain.pos(i, 0);
double y = domain.pos(i, 1);
1.0_i * op.lap(i) = -2*PI*PI*sin(PI * x)*sin(PI * y);
}
for (int i : domain.boundary()) {
op.value(i) = 0.0;
}
</syntaxhighlight>

The solution $u(x,y)$ is plotted below. In our case solution consist only of the imaginary part, the real part is equal to zero.

[[File:complex_poisson_2D_real.png|500px]]
[[File:complex_poisson_2D_imag.png|500px]]

Go back to [[Medusa#Examples|Examples]].

File:Complex poisson 2D imag.png

2019-07-19T14:03:20Z

APirnat: File uploaded with MsUpload

File uploaded with MsUpload

File:Complex poisson 2D real.png

2019-07-19T14:03:17Z

APirnat: File uploaded with MsUpload

File uploaded with MsUpload

Schrödinger equation

2019-07-15T14:35:30Z

APirnat:

Go back to [[Medusa#Examples|Examples]].

__FORCETOC__

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

<syntaxhighlight lang="cpp">
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;
}
</syntaxhighlight>

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.

<syntaxhighlight lang="cpp">
// Solve matrix system.
psi = solver.solve(rhs);

// Update rhs.
rhs = psi;
</syntaxhighlight>

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

[[File:video (1).mp4|400px]]

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.

[[File:infinite_well_2D_error.png|400px]]
[[File:infinite_well_2D_error_fit.png|400px]]

== 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.

[[File:infinite_well_3D_rsol.png|400px]]
[[File:infinite_well_3D_csol.png|400px]]

Go back to [[Medusa#Examples|Examples]].

Schrödinger equation

2019-07-15T14:21:31Z

APirnat:

Go back to [[Medusa#Examples|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$.

<syntaxhighlight lang="cpp">
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;
}
</syntaxhighlight>

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.

<syntaxhighlight lang="cpp">
// Solve matrix system.
psi = solver.solve(rhs);

// Update rhs.
rhs = psi;
</syntaxhighlight>

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

[[File:video (1).mp4|400px]]

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.

[[File:infinite_well_2D_error.png|400px]]
[[File:infinite_well_2D_error_fit.png|400px]]

== 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, only the domain changes to unit cube, $\Omega=[0, 1]\times[0, 1]\times [0, 1]$.
For initial state $\Psi(x, y, z,t= 0)$ 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 presented bellow.

[[File:infinite_well_3D_rsol.png|400px]]
[[File:infinite_well_3D_csol.png|400px]]

Go back to [[Medusa#Examples|Examples]].

Schrödinger equation

2019-07-15T14:19:00Z

APirnat:

Go back to [[Medusa#Examples|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$.

<syntaxhighlight lang="cpp">
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;
}
</syntaxhighlight>

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.

<syntaxhighlight lang="cpp">
// Solve matrix system.
psi = solver.solve(rhs);

// Update rhs.
rhs = psi;
</syntaxhighlight>

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

[[File:video (1).mp4|400px]]

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.

[[File:infinite_well_2D_error.png|400px]]
[[File:infinite_well_2D_error_fit.png|400px]]

== 3D infinite cubic well ==
Now we are dealing with a problem very similar to the previus one, the only diffrenece is, that in this case we have a quantom particle in a 3-dimensional box. 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, only the domain changes to unit cube, $\Omega=[0, 1]\times[0, 1]\times [0, 1]$.
For initial state $\Psi(x, y, z,t= 0)$ we choose
\[
\Psi(x,y,t=0) = \sin{(\pi x)} \sin{(\pi y) \sin(\piz)}.
\]
The plot of the real and imaginary part of the solution $\Psi(x, y, z)$ at time $T=0.1$ is plotted bellow.

[[File:infinite_well_3D_rsol.png|400px]]
[[File:infinite_well_3D_csol.png|400px]]

Go back to [[Medusa#Examples|Examples]].

File:Infinite well 3D csol.png

2019-07-15T14:18:33Z

APirnat: File uploaded with MsUpload

File uploaded with MsUpload

File:Infinite well 3D rsol.png

2019-07-15T14:18:31Z

APirnat: File uploaded with MsUpload

File uploaded with MsUpload

Schrödinger equation

2019-07-15T13:17:39Z

APirnat:

File:Infinite well 2D error fit.png

2019-07-15T13:15:54Z

APirnat: File uploaded with MsUpload

File uploaded with MsUpload

File:Infinite well 2D error.png

2019-07-15T13:15:53Z

APirnat: File uploaded with MsUpload

File uploaded with MsUpload

File:2D xerr convergence monomi 1.pdf

2019-07-15T13:07:05Z

APirnat: File uploaded with MsUpload

File uploaded with MsUpload

File:2D fit monomials 1.pdf

2019-07-15T13:07:03Z

APirnat: File uploaded with MsUpload

File uploaded with MsUpload

Schrödinger equation

2019-06-17T20:33:35Z

APirnat:

Go back to [[Medusa#Examples|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$ .

<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);
}

// set equation on boundary
for (int i: domain.boundary()){
op.value(i) = 0;
}
</syntaxhighlight>

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.

<syntaxhighlight lang="c++" line>
// solve matrix system
psi = solver.solve(rhs);

// update rhs
rhs = psi;
</syntaxhighlight>

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

[[File:video (1).mp4|400px]]

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.

[[File:infinite_well_2D_convergence.png|400px]]
[[File:infinite_well_2D_convergence_fit.png|400px]]

Go back to [[Medusa#Examples|Examples]].

Schrödinger equation

2019-06-17T20:32:14Z

APirnat:

Go back to [[Medusa#Examples|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$ .

<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);
}

// set equation on boundary
for (int i: domain.boundary()){
op.value(i) = 0;
}
</syntaxhighlight>

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.

<syntaxhighlight lang="c++" line>
// solve matrix system
psi = solver.solve(rhs);

// update rhs
rhs = psi;
</syntaxhighlight>

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

[[File:video (1).mp4|400px]]

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.

[[File:infinite_well_2D_convergence.png|400px]]
[[File:infinite_well_2D_convergence_fit.png|400px]]

Go back to [[Medusa#Examples|Examples]].

File:Infinite well 2D convergence fit.png

2019-06-17T20:31:04Z

APirnat: File uploaded with MsUpload

File uploaded with MsUpload

File:Infinite well 2D convergence.png

2019-06-17T20:30:02Z

APirnat: File uploaded with MsUpload

File uploaded with MsUpload

Schrödinger equation

2019-06-17T20:22:53Z

APirnat:

Go back to [[Medusa#Examples|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$ .

<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);
}

// set equation on boundary
for (int i: domain.boundary()){
op.value(i) = 0;
}
</syntaxhighlight>

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.

<syntaxhighlight lang="c++" line>
// solve matrix system
psi = solver.solve(rhs);

// update rhs
rhs = psi;
</syntaxhighlight>

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

[[File:video (1).mp4|400px]]

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.

[[File:infinite_well_2D_convergence_fit.pdf|400px]]
[[File:infinite_well_2D_xerr_convergence.pdf|400px]]

Go back to [[Medusa#Examples|Examples]].

File:Infinite well 2D xerr convergence.pdf

2019-06-17T20:21:29Z

APirnat: File uploaded with MsUpload

File uploaded with MsUpload

File:Infinite well 2D convergence fit.pdf

2019-06-17T20:21:28Z

APirnat: File uploaded with MsUpload

File uploaded with MsUpload

Schrödinger equation

2019-06-17T09:44:43Z

APirnat:

Schrödinger equation

2019-06-17T09:39:21Z

APirnat:

Schrödinger equation

2019-06-17T09:37:50Z

APirnat:

Schrödinger equation

2019-06-17T09:31:29Z

APirnat:

Schrödinger equation

2019-06-17T09:29:41Z

APirnat:

File:Video (1).mp4

2019-06-17T09:27:47Z

APirnat:

Schrödinger equation

2019-06-17T09:27:26Z

APirnat:

Schrödinger equation

2019-06-17T09:25:33Z

APirnat:

Schrödinger equation

2019-05-21T09:20:29Z

APirnat:

Solution procedure is still compiling ... so please wait for results :)

= Introduction =
The quantum world is governed by the [https://en.wikipedia.org/wiki/Schr%C3%B6dinger_equation Schrödinger equation]

\[{\displaystyle {\hat {H}}|\psi (t)\rangle =i\hbar {\frac {\partial }{\partial t}}|\psi (t)\rangle } \]

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].

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

\[\hat T = - \frac{\hbar^2}{2m} \nabla^2 .\]

The final version of the single particle Schrödinger equation can be written as

\[\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) \]

Quantum state function is a complex function, so it is usually split into the real part and imaginary part

\[ u, v \in C(\mathbb R)\colon \psi = u + i v , \]

which for a real $V$ yields a system of two real equations

\[\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) , \]
\[\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.

= Particle in a box =

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}}}.\]

Schrödinger equation

2019-05-21T09:17:41Z

APirnat:

Solution procedure is still compiling ... so please wait for results :)

= Introduction =
The quantum world is governed by the [https://en.wikipedia.org/wiki/Schr%C3%B6dinger_equation Schrödinger equation]

\[{\displaystyle {\hat {H}}|\psi (t)\rangle =i\hbar {\frac {\partial }{\partial t}}|\psi (t)\rangle } \]

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].

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

\[\hat T = - \frac{\hbar^2}{2m} \nabla^2 .\]

The final version of the single particle Schrödinger equation can be written as

\[\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) \]

Quantum state function is a complex function, so it is usually split into the real part and imaginary part

\[ u, v \in C(\mathbb R)\colon \psi = u + i v , \]

which for a real $V$ yields a system of two real equations

\[\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) , \]
\[\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.

= Particle in a box =

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 similar to Harmonic oscilator

\[\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}}}.\]

Schrödinger equation

2019-05-21T09:14:01Z

APirnat: Created page with "Solution procedure is still compiling ... so please wait for results :) = Introduction = The quantum world is governed by the [https://en.wikipedia.org/wiki/Schr%C3%B6dinger..."

Solution procedure is still compiling ... so please wait for results :)

= Introduction =
The quantum world is governed by the [https://en.wikipedia.org/wiki/Schr%C3%B6dinger_equation Schrödinger equation]

\[{\displaystyle {\hat {H}}|\psi (t)\rangle =i\hbar {\frac {\partial }{\partial t}}|\psi (t)\rangle } \]

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].

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

\[\hat T = - \frac{\hbar^2}{2m} \nabla^2 .\]

The final version of the single particle Schrödinger equation can be written as

\[\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) \]

Quantum state function is a complex function, so it is usually split into the real part and imaginary part

\[ u, v \in C(\mathbb R)\colon \psi = u + i v , \]

which for a real $V$ yields a system of two real equations

\[\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) , \]
\[\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.

= Particle in a box =

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 similar to Harmonic oscilator

\[\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}}}.\]

Quantum Mechanics

2019-05-21T09:12:59Z

APirnat: Blanked the page