]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Finish writing documentation.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 14 Jul 2009 21:04:59 +0000 (21:04 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 14 Jul 2009 21:04:59 +0000 (21:04 +0000)
git-svn-id: https://svn.dealii.org/trunk@19081 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-36/doc/intro.dox
deal.II/examples/step-36/doc/results.dox

index 5f665d5ac36701b40e19488f26ff4427ff640028..685cdbb8d09bea7804740691a72ced15fcf973a9 100644 (file)
@@ -7,66 +7,91 @@ Bangerth.  </i>
 <h1>Preamble</h1>
 
 The problem we want to solve in this example is an eigenspectrum
-problem which no longer contains a single solution, but a set of
-solutions for the various eigenfunctions we want to compute. The
-problem of finding all eigenvalues (eigenfunctions) of a generalized
-eigenvalue problem is a formidable challange; however, most of the
-time we are really only interested in a small subset of these values
-(functions). Fortunately, the interface to the SLEPc library allows us
-to select which portion of the eigenspectrum and how many solutions we
-want to solve for.
-
-To do this, everything here is built on top of classes provided by
-deal.II that wrap around the linear algebra implementation of the <a
-href="http://www.grycap.upv.es/slepc/" target="_top">SLEPc</a>
-library; which links to some the underlying features of the <a
-href="http://www.mcs.anl.gov/petsc/" target="_top">PETSc</a> library.
+problem. Eigenvalue problems appear in a wide context of problems, for example
+in the computation of electromagnetic standing waves in cavities, vibration
+modes of drum membranes, or oscillations of lakes and estuaries. One of the
+most enigmatic applications is probably the computation of stationary or
+quasi-static wave functions in quantum mechanics. The latter application is
+what we would like to investigate here, though the general techniques outlined
+in this program are of course equally applicable to the other applications
+above. 
+
+Eigenspectrum problems have the general form
+@f{align*}
+       L \Psi &= \varepsilon \Psi \qquad &&\text{in}\ \Omega,
+       \\
+       \Psi &= 0 &&\text{on}\ \partial\Omega,
+@f}
+where the Dirichlet boundary condition on $\Psi=\Psi(\mathbf x)$ could also be
+replaced by Neumann or Robin conditions; $L$ is an operator that generally
+also contains differential operators.
+
+Under suitable conditions, above equations have a set of solutions
+$\Psi_\ell,\varepsilon_\ell$, $\ell\in {\cal I}$, where $\cal I$ can be a finite or
+infinite set. In either case, let us note that there is no longer just a
+single solution, but a set of solutions (the various eigenfunctions and
+corresponding eigenvalues) that we want to compute. The problem of finding all
+eigenvalues (eigenfunctions) of such eigenvalue problems is a formidable
+challange; however, most of the time we are really only interested in a small
+subset of these values (functions). Fortunately, the interface to the SLEPc
+library that we will use for this tutorial program allows us to select which
+portion of the eigenspectrum and how many solutions we want to solve for.
+
+In this program, the eigenspectrum solvers we use are classes
+provided by deal.II that wrap around the linear algebra implementation of the
+<a href="http://www.grycap.upv.es/slepc/" target="_top">SLEPc</a> library;
+SLEPc itself builds on the <a
+href="http://www.mcs.anl.gov/petsc/" target="_top">PETSc</a> library for
+linear algebra contents.
+
 
 <a name="Introduction"></a>
 <h1>Introduction</h1>
 
-Start with a differential operator equation
-@f[
-       L(x,y) \Psi = \varepsilon \Psi \quad,
-@f]
-where $L$ is a differential operator. In the usual way in finite
-elements we multiply both sides with a test function and then replace
-$\Psi$ and the test function with discrete shape functions.
-@f[
-       \sum_j (\phi_i, L\phi_j) \tilde{\psi}_j =
-       \varepsilon \sum_j (\phi_i, \phi_j) \tilde{\psi}_j \quad,
-@f]
-where $\tilde{\psi}_j$, the reduced wavefunction, are the expansion
-coefficients of the approximation
-
+The basic equation of stationary quantum mechanics is the Schr&ouml;dinger
+equation. The Copenhagen interpretation of quantum mechanics posits that the
+motion of particles in an external potential $V(\mathbf x)$ is governed a wave
+function $\Psi(\mathbf x)$ that satisfies this Schr&ouml;dinger
+equation of the (non-dimensionalized) form
+@f{align*}
+  [-\Delta + V(\mathbf x)] \Psi(\mathbf x) &= \varepsilon \Psi(\mathbf x)
+  \qquad &&\text{in}\ \Omega,
+  \\
+  \Psi &= 0 &&\text{on}\ \partial\Omega.
+@f}
+As a consequence of this eigenvalue problem, this particle can only exist in a
+certain number of eigenstates that correspond to the energy eigenvalues
+$\varepsilon_\ell$ admitted as solutions of this equation, and if a particle has
+energy $\varepsilon_\ell$ then the probability of finding it at location $\mathbf
+x$ is proportional to $|\Psi_\ell(\mathbf x)|^2$ where $\Psi_\ell$ is the
+eigenfunction that corresponds to this eigenvalue.
+
+In order to numerically find solutions to this equation, i.e. a set of pairs
+of eigenvalue/eigenfunction, we use the usual finite element approach of
+multiplying the equation from the left with testfunctions, integrating by
+parts, and searching for solutions in finite dimensional spaces by
+approximating $\Psi(\mathbf x)\approx\Psi_h(\mathbf x)=\sum_{j}\phi_j(\mathbf
+x)\tilde\psi_j$, where $\tilde\psi$ is a vector of expansion coefficients. We
+then immediately arrive at the following equation that discretizes the
+continuous eigenvalue problem:
 @f[
-       \psi_i = \sum_j \phi_j \tilde{\psi}_j \quad.
+       \sum_j [(\nabla\phi_i, \nabla\phi_j)+(V(\mathbf x)\phi_i,\phi_j)]
+          \tilde{\psi}_j =
+       \varepsilon_h \sum_j (\phi_i, \phi_j) \tilde{\psi}_j.
 @f]
-It can easily be shown that the consequence of this, is that we want
-to solve the generalized eigenvalue problem
+In matrix and vector notation, this equation then reads:
 @f[
-       A \tilde{\Phi} = \varepsilon M \tilde{\Phi} \quad,
+       A \tilde{\Psi} = \varepsilon_h M \tilde{\Psi} \quad,
 @f]
-where $A$ si the stiffness matrix arising from the differential
+where $A$ is the stiffness matrix arising from the differential
 operator $L$, and $M$ is the mass matrix. The solution to the
-eigenvalue problem is an eigenspectrum $\varepsilon_j$, with
-associated eigenfunctions $\tilde{\Phi}=\sum_j \tilde{\phi}_j$.
+eigenvalue problem is an eigenspectrum $\varepsilon_{h,\ell}$, with
+associated eigenfunctions $\tilde{\Psi}_\ell=\sum_j \phi_j\tilde{\psi}_j$.
 
-To solve an actual physical problem we next define the differential
-operator $L$ in a way that mimicks the wave-equation of Schr\"odinger:
-@f[
-       L({\mathbf x}) = K({\mathbf x}) 
-       + V({\mathbf x})\quad,
-@f]
-and let $K({\mathbf x}) = \partial^2_{\mathbf x}$, and
-@f[
-       V({\mathbf x}) = \quad,
-@f]
+It is this form of the eigenvalue problem that involves both matrices $A$ and
+$M$ that we will solve in the current tutorial program. We will want to solve
+it for the lowermost few eigenvalue/eigenfunction pairs.
 
-In the parlance of Schroedinger eigenvalue problems, the solutions
-sought are the wavefunctions $\phi_j$. In finite elements the
-solutions are $\tilde{\phi}_j$, known as a reduced (or tilde-basis)
-wavefunctions.
 
 <h3>Implementation details</h3>
 
index 347bc65158c667ec94cb9d973872e8a1101a8254..60bc483ccec2e5a57a7bd18a07af32e7d3e8a2d9 100644 (file)
 <h1>Results</h1>
 
-<h1>Possibilities for extensions</h1>
+<h2>Running the problem</h2>
+
+The problem's input is parameterized by an input file <code>step-36.prm</code>
+which could, for example, contain the following text:
+
+@code
+set Global mesh refinement steps         = 5
+set Number of eigenvalues/eigenfunctions = 5
+set Potential                            = 0
+@endcode
+
+Here, the potential is zero inside the domain, and we know that the
+eigenvalues are given by $\lambda_{(mn)}=\frac{\pi}{4}(m^2+n^2)$ where
+$m,n\in{\mathbb N^+}$. Eigenfunctions are sines and cosines with $m$ and $n$
+periods in $x$ and $y$ directions. This matches the output our program
+generates:
+@code
+examples/step-36> make run
+============================ Running step-36
+   Number of active cells:       1024
+   Number of degrees of freedom: 1089
+
+   Eigenvalue 0 : 4.93877
+   Eigenvalue 1 : 12.3707
+   Eigenvalue 2 : 12.3707
+   Eigenvalue 3 : 19.8027
+   Eigenvalue 4 : 24.837
+
+Job done.
+@endcode
+These eigenvalues are exactly the ones that correspond to pairs $(m,n)=(1,1)$,
+$(1,2)$ and $(2,1)$, $(2,2)$, and $(3,1)$. A visualization of the
+corresponding eigenfunctions would look like this:
+
+<TABLE WIDTH="100%">
+<tr>
+<td>
+  @image html step-36.default.eigenfunction.0.png
+</td>
+<td>
+  @image html step-36.default.eigenfunction.1.png
+</td>
+</tr>
+
+<tr>
+<td>
+  @image html step-36.default.eigenfunction.2.png
+</td>
+<td>
+  @image html step-36.default.eigenfunction.3.png
+</td>
+</tr>
+
+<tr>
+<td>
+  @image html step-36.default.eigenfunction.4.png
+</td>
+<td>
+</td>
+</tr>
+</table>
+
+
+
+<h2>Possibilities for extensions</h2>
 
 It is always worth playing a few games in the playground! So here goes
 with a few suggestions:
 
 <ul> 
 
-<li> What happens to the world of the particle in a box if the box
-is bigger than we have considered in this application?</li>
+<li> The potential used above (called the <i>infinite well</i> because it is a
+flat potential surrounded by infinitely high walls, see below) is interesting
+because it allows for analytically known solutions. Apart from that, it is
+rather boring, however. That said, it is trivial to play around with the
+potential by just setting it to something different in the input file. For
+example, let us assume that we wanted to work with the following potential in
+2d: 
+@f[
+  V(x,y) = \left\{ 
+       \begin{array}{ll}
+         -100 & \text{if}\ \sqrt{x^2+y^2}<\frac 34 \ \text{and} 
+                         \ xy>0,
+         \\
+         -5 & \text{if}\ \sqrt{x^2+y^2}<\frac 34 \ \text{and} 
+                         \ xy\le 0,
+         \\
+         0 & \text{otherwise.}
+      \end{array} \right.
+@f]
+In other words, the potential is -100 in two sectors of a circle of radius
+0.75, -5 in the other two sectors, and zero outside the circle. We can achieve
+this by using the following in the input file:
+@code
+set Potential = if (x^2 + y^2 < 0.75^2, if (x*y > 0, -100, -5), 0)
+@endcode
+If in addition we also increase the mesh refinement by one level, we get the
+following results: 
+@code
+examples/step-36> make run
+============================ Running step-36
+   Number of active cells:       4096
+   Number of degrees of freedom: 4225
+
+   Eigenvalue 0 : -74.2562
+   Eigenvalue 1 : -72.7322
+   Eigenvalue 2 : -42.7406
+   Eigenvalue 3 : -42.2232
+   Eigenvalue 4 : -37.0744
+@endcode
+
+The output file also contains an interpolated version of the potential, which
+looks like this (note that as expected the lowest few eigenmodes have
+probability densities $|\Psi(\mathbf x)|^2$ that are significant only where the
+potential is the lowest, i.e. in the top right and bottom left sector of inner
+circle of the potential):
+
+@image html step-36.mod.potential.png
+
+The first five eigenfunctions are now like this:
+
+<TABLE WIDTH="100%">
+<tr>
+<td>
+  @image html step-36.mod.eigenfunction.0.png
+</td>
+<td>
+  @image html step-36.mod.eigenfunction.1.png
+</td>
+</tr>
+
+<tr>
+<td>
+  @image html step-36.mod.eigenfunction.2.png
+</td>
+<td>
+  @image html step-36.mod.eigenfunction.3.png
+</td>
+</tr>
+
+<tr>
+<td>
+  @image html step-36.mod.eigenfunction.4.png
+</td>
+<td>
+</td>
+</tr>
+</table>
+
+
+<li> In our derivation of the problem we have assumed that the particle is
+confined to a domain $\Omega$ and that at the boundary of this domain its
+probability $|\Psi|^2$ of being is zero. This is equivalent to solving the
+eigenvalue problem on all of ${\mathbb R}^d$ and assuming that
+the energy potential is finite only inside a region $\Omega$ and infinite
+outside. It is relatively easy to show that $|\Psi(\mathbf x)|^2$ at all
+locations $\mathbf x$ where $V(\mathbf x)=\infty$. So the question is what
+happens if our potential is not of this form, i.e. there is no bounded domain
+outside of which the potential is infinite? In that case, it may be worth to
+just consider a very large domain at the boundary of which $V(\mathbf x)$ is
+at least very large, if not infinite. Play around with a few cases like this
+and explore how spectrum and eigenfunction change as we make the computational
+region larger and larger.
 
 <li> What happens if we investigate the simple harmonic oscillator
-problem? This is a hyper-ballic symmetric function and as such it
-may be worth trying to change the geometry of the mesh to reflect the
-underlying potential which governs quantum mechanical problems.</li>
+problem $V(\mathbf x)=c|\mathbf x|^2$? This potential is exactly of the form
+discussed in the previous paragraph and one may want to use a large spherical
+domain to approximate the whole-space problem.
 
 <li> What happens if the particle in the box has <i>internal</i>
 degrees of freedom? For example, if the particle were a spin-$1/2$

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.