From: Toby D. Young Date: Fri, 17 Jul 2009 10:15:51 +0000 (+0000) Subject: Tweaked documentation of step-36 X-Git-Tag: v8.0.0~7449 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c78e34ab04ebcc972acd0a8f63318080f3686f79;p=dealii.git Tweaked documentation of step-36 git-svn-id: https://svn.dealii.org/trunk@19114 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-36/doc/intro.dox b/deal.II/examples/step-36/doc/intro.dox index 685cdbb8d0..1df6e81ec8 100644 --- a/deal.II/examples/step-36/doc/intro.dox +++ b/deal.II/examples/step-36/doc/intro.dox @@ -3,95 +3,93 @@ This program was contributed by Toby D. Young and Wolfgang Bangerth. - +

Preamble

The problem we want to solve in this example is an eigenspectrum -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. +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, + L \Psi &= \varepsilon \Psi \qquad &&\text{in}\ \Omega\quad, \\ - \Psi &= 0 &&\text{on}\ \partial\Omega, + \Psi &= 0 &&\text{on}\ \partial\Omega\quad, @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 -SLEPc library; -SLEPc itself builds on the PETSc library for -linear algebra contents. - - - +Under suitable conditions, the 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 numerically finding all eigenvalues +(eigenfunctions) of such eigenvalue problems is a formidable +challange. In fact, if the set $\cal I$ is infinite, the callange is +insurmountable. Most of the time however we are really only +interested in a small subset of these values (functions); and +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 +SLEPc +library; SLEPc itself builds on the PETSc library +for linear algebra contents. + +

Introduction

-The basic equation of stationary quantum mechanics is the Schrö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ö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[ - \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] -In matrix and vector notation, this equation then reads: -@f[ - A \tilde{\Psi} = \varepsilon_h M \tilde{\Psi} \quad, -@f] -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_{h,\ell}$, with -associated eigenfunctions $\tilde{\Psi}_\ell=\sum_j \phi_j\tilde{\psi}_j$. - -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. - +The basic equation of stationary quantum mechanics is the +Schrödinger equation which models the motion of particles in an +external potential $V(\mathbf x)$ as if it were governed a wave +function $\Psi(\mathbf x)$ that satisfies a relation of the +(non-dimensionalized ) form +@f{align*} [-\Delta + V(\mathbf x)] +\Psi(\mathbf x) &= \varepsilon \Psi(\mathbf x) \qquad &&\text{in}\ +\Omega\quad, \\ \Psi &= 0 &&\text{on}\ \partial\Omega\quad. +@f} +As a consequence, 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. The +Copenhagen interpretation of quantum mechanics posits that, 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 eigenvalues/eigenfunctions, 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[ \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\quad. @f] In +matrix and vector notation, this equation then reads: @f[ A +\tilde{\Psi} = \varepsilon_h M \tilde{\Psi} \quad, @f] 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_{h,\ell}$, with associated eigenfunctions +$\Psi_\ell=\sum_j \phi_j\tilde{\psi}_j$. + +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.

Implementation details

diff --git a/deal.II/examples/step-36/doc/results.dox b/deal.II/examples/step-36/doc/results.dox index 60bc483cce..1a43689c6c 100644 --- a/deal.II/examples/step-36/doc/results.dox +++ b/deal.II/examples/step-36/doc/results.dox @@ -28,11 +28,10 @@ examples/step-36> make run 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: +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: @@ -62,8 +61,6 @@ corresponding eigenfunctions would look like this:
- -

Possibilities for extensions

It is always worth playing a few games in the playground! So here goes @@ -71,24 +68,25 @@ with a few suggestions: diff --git a/deal.II/examples/step-36/step-36.cc b/deal.II/examples/step-36/step-36.cc index f855166e11..33568dd10f 100644 --- a/deal.II/examples/step-36/step-36.cc +++ b/deal.II/examples/step-36/step-36.cc @@ -90,10 +90,10 @@ class EigenvalueProblem // a mass matrix for the right // hand side. We also need not // just one solution function, - // but a whole set of those for + // but a whole set of these for // the eigenfunctions we want to // compute, along with the - // corresponding eigenvectors: + // corresponding eigenvalues: PETScWrappers::SparseMatrix stiffness_matrix, mass_matrix; std::vector eigenfunctions; std::vector eigenvalues; @@ -122,7 +122,7 @@ class EigenvalueProblem // @sect4{EigenvalueProblem::EigenvalueProblem} // First up, the constructor. The - // main, new part is handling the + // main new part is handling the // run-time input parameters. We need // to declare their existence first, // and then read their values from @@ -258,22 +258,25 @@ void EigenvalueProblem::make_grid_and_dofs () // \int_K \nabla\varphi_i(\mathbf x) // \cdot \nabla\varphi_j(\mathbf x) + // V(\mathbf x)\varphi_i(\mathbf - // x)\varphi_j(\mathbf x)$. The - // function should be immediately - // familiar if you've seen previous - // tutorial programs. The only thing - // new would be setting up an object - // that described the potential - // $V(\mathbf x)$ using the - // expression that we got from the - // input file. We then need to - // evaluate this object at the - // quadrature points on each cell. If - // you've seen how to evaluate - // function objects (see, for example - // the coefficient in step-5), the - // code here will also look rather - // familiar. + // x)\varphi_j(\mathbf x)$ and + // $M^K_{ij} = \int_K + // \varphi_i(\mathbf + // x)\varphi_j(\mathbf x)$ + // respectively. This function should + // be immediately familiar if you've + // seen previous tutorial + // programs. The only thing new would + // be setting up an object that + // described the potential $V(\mathbf + // x)$ using the expression that we + // got from the input file. We then + // need to evaluate this object at + // the quadrature points on each + // cell. If you've seen how to + // evaluate function objects (see, + // for example the coefficient in + // step-5), the code here will also + // look rather familiar. template void EigenvalueProblem::assemble_system () { @@ -467,6 +470,21 @@ void EigenvalueProblem::solve () // visualization. It works as in many // of the other tutorial programs. // + // The whole collection of functions + // is then output as a single VTK + // file. +template +void EigenvalueProblem::output_results () const +{ + DataOut data_out; + + data_out.attach_dof_handler (dof_handler); + + for (unsigned int i=0; i::solve () // space. The result we also attach // to the DataOut object for // visualization. - // - // The whole collection of functions - // is then output as a single VTK - // file. -template -void EigenvalueProblem::output_results () const -{ - DataOut data_out; - - data_out.attach_dof_handler (dof_handler); - - for (unsigned int i=0; i projected_potential (dof_handler.n_dofs()); { FunctionParser potential; @@ -604,8 +605,10 @@ int main (int argc, char **argv) return 1; } - // ...or show that we are happy by - // exiting nicely: + // If no exceptions are thrown, + // then we can tell the program to + // stop monkeying around and exit + // nicely: std::cout << std::endl << "Job done." << std::endl;