demo_eigenvalue.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A simple eigenvalue solver\n",
    "==========================\n",
    "\n",
    "We recommend that you are familiar with the demo for the Poisson\n",
    "equation before looking at this demo.\n",
    "\n",
    "If you want a more complex problem, we suggest that you look at the\n",
    "other eigenvalue demo.\n",
    "\n",
    "Implementation\n",
    "--------------\n",
    "\n",
    "This demo is implemented in a single Python file, demo\\_eigenvalue.py\n",
    "and reads in an external mesh box\\_with\\_dent.xml.gz.\n",
    "\n",
    "The eigensolver functionality in DOLFIN relies on the library SLEPc\n",
    "which in turn relies on the linear algebra library PETSc. Therefore,\n",
    "both PETSc and SLEPc are required for this demo. We can test whether\n",
    "PETSc and SLEPc are available, and exit if not, as follows :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dolfin import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Specifying matplotlib as backend:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "parameters[\"plotting_backend\"]=\"matplotlib\"\n",
    "# Test for PETSc and SLEPc\n",
    "\n",
    "if not has_linear_algebra_backend(\"PETSc\"):\n",
    "    print(\"DOLFIN has not been configured with PETSc. Exiting.\")\n",
    "    exit()\n",
    "\n",
    "if not has_slepc():\n",
    "    print(\"DOLFIN has not been configured with SLEPc. Exiting.\")\n",
    "    exit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we need to construct the matrix $A$. This will be done in three\n",
    "steps: defining the mesh and the function space associated with it;\n",
    "constructing the variational form defining the matrix; and then\n",
    "assembling this form. The code is shown below :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define mesh, function space\n",
    "mesh = Mesh(\"box_with_dent.xml.gz\")\n",
    "V = FunctionSpace(mesh, \"Lagrange\", 1)\n",
    "\n",
    "# Define basis and bilinear form\n",
    "u = TrialFunction(V)\n",
    "v = TestFunction(V)\n",
    "a = dot(grad(u), grad(v))*dx\n",
    "\n",
    "# Assemble stiffness form\n",
    "A = PETScMatrix()\n",
    "assemble(a, tensor=A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we (in this example) first define the matrix `A` as a\n",
    ":pyPETScMatrix &lt;dolfin.cpp.PETScMatrix&gt; and then assemble the form\n",
    "into it. This is an easy way to ensure that the matrix has the right\n",
    "type.\n",
    "\n",
    "In order to solve the eigenproblem, we need to define an eigensolver. To\n",
    "solve a standard eigenvalue problem, the eigensolver is initialized with\n",
    "a single argument, namely the matrix `A`. :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create eigensolver\n",
    "eigensolver = SLEPcEigenSolver(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we ready solve the eigenproblem by calling the :pysolve\n",
    "&lt;dolfin.cpp.SLEPcEigenSolver.solve&gt; method of the eigensolver.\n",
    "Note that eigenvalue problems tend to be computationally intensive and\n",
    "may hence take a while. :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute all eigenvalues of A x = \\lambda x\n",
    "print(\"Computing eigenvalues. This can take a minute.\")\n",
    "eigensolver.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is kept by the eigensolver, but can fortunately be extracted.\n",
    "Here, we have computed all eigenvalues, and they will be sorted by\n",
    "largest magnitude. We can extract the real and complex part (`r` and\n",
    "`c`) of the largest eigenvalue and the real and complex part of the\n",
    "corresponding eigenvector (`ru` and `cu`) by asking for the first\n",
    "eigenpair as follows :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract largest (first) eigenpair\n",
    "r, c, rx, cx = eigensolver.get_eigenpair(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we want to examine the results. The eigenvalue can easily be\n",
    "printed. But, the real part of eigenvector is probably most easily\n",
    "visualized by constructing the corresponding eigenfunction. This can be\n",
    "done by creating a :pyFunction\n",
    "&lt;dolfin.functions.function.Function&gt; in the function space `V` and\n",
    "the associating eigenvector `rx` with the Function. Then the\n",
    "eigenfunction can be manipulated as any other :pyFunction\n",
    "&lt;dolfin.functions.function.Function&gt;, and in particular plotted :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Largest eigenvalue: \", r)\n",
    "\n",
    "# Initialize function and assign eigenvector\n",
    "u = Function(V)\n",
    "u.vector()[:] = rx\n",
    "\n",
    "# Plot eigenfunction\n",
    "plt.figure()\n",
    "plot(u)\n",
    "interactive()"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 1
}