00-sympy.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import sympy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Symbolic algebra\n",
    "\n",
    "Here we'll demonstrate a bit about the computer algebra system [SymPy](https://www.sympy.org).\n",
    "SymPy does symbolic math, much like Maple or Mathematica, entirely in Python.\n",
    "The core data type in SymPy is a *symbol*, from which you can build more complicated expressions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = sympy.symbols('x')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = sympy.exp(x) + sympy.sin(x)\n",
    "f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we've created a symbolic representation of the function $x \\cdot \\sin x$, we'll want some way to actually plug in values of $x$ and evaluate it.\n",
    "The simplest way to do this is to use the method `subs` of the expression `f`, which will substitute in numeric values for `x`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f.subs(x, 1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SymPy also has symbolic representations of transcendental numbers like $e$ and $\\pi$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f.subs(x, sympy.pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we used a numerical approximation of $\\pi$, which we can get from the package numpy, we'll get an inexact answer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "f.subs(x, np.pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What this shows is that SymPy is smart enough to recognize identities like $\\sin(k \\cdot \\pi) = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `subs` method works perfectly fine, but if we need to evaluate this function repeatedly on a lot of numeric arguments it's quite slow.\n",
    "So there's another function call `lambdify` that will take in a symbolic function and return a numeric one.\n",
    "The lambdified version is much faster to call over and over again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F = sympy.lambdify(x, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F(np.pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you know you'll be using numpy for doing big array operations, you can tell sympy that you want a lambdified function that'll work with numpy too:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F = sympy.lambdify(x, f, modules='numpy')\n",
    "X = np.array([0., 1/6, 1/4, 1/3, 1/2, 2/3, 3/4, 1.]) * np.pi\n",
    "F(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Symbolic differentiation\n",
    "\n",
    "One of the most useful features of these kinds of tools is that they can calculate symbolic derivatives."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sympy.diff(f, x)\n",
    "g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we could only calculate derivatives numerically by the method of finite differences, we could get really bad results because of floating-point underflow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "G = sympy.lambdify(x, g, modules='numpy')\n",
    "G(0.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(F(1e-16) - F(0.)) / 1e-16"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observe how if we put in a symbolic argument we get a symbolic result, and if we put in a numeric argument we get a numeric result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g.subs(x, sympy.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g.subs(x, np.pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Symbolic integration\n",
    "\n",
    "Differentiation is a \"nice\" problem in that, if you have an symbolic expression for an elementary function, its derivative is elementary and there is an algorithm that will terminate in finite time to come up with a symbolic expression for its derivative.\n",
    "Integration is not nice.\n",
    "There are elementary functions with no elementary anti-derivative, for example $f(x) = e^{-x^2}$.\n",
    "Nonetheless, SymPy and other symbolic algebra systems usually have a rich set of heuristics for calculating anti-derivatives."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sympy.integrate(f, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I tried to stump the algorithm but they have so many special functions built in that it's kind of hard!\n",
    "For example, the anti-derivative of $e^{-x^2}$ isn't elementary, but it's important enough that people decided they'd call it the *error function* (this comes from statistics)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h = sympy.exp(-x**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sympy.integrate(h, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to actually evaluate the integral over specified bounds rather than get an indefinite integral, you can pass those integration bounds.\n",
    "The variable `sympy.oo` stands for $\\infty$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.float64(sympy.integrate(h, (x, 0., sympy.oo)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Several variables\n",
    "\n",
    "So far we've only looked at functions of a single variable $x$.\n",
    "But we could just as easily use more symbols to look at functions of several variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = sympy.symbols('x y')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = x**2 - y**2\n",
    "f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sympy.diff(f, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sympy.diff(f, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can calculate higher-order derivatives and mixed derivatives too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sympy.diff(f, x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sympy.diff(f, x, x) + sympy.diff(f, y, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a *very* handy trick when you want to test whether your PDE solver is working right.\n",
    "For example, let's say you wrote some code to solve the Poisson equation\n",
    "\n",
    "$$-\\nabla\\cdot k\\nabla u = f$$\n",
    "\n",
    "where $k$ is a conductivity coefficient and $f$ is a right-hand side.\n",
    "You could painstakingly pick some $k$ and $f$ so that this problem has an exact solution.\n",
    "Or you could take any random old $k$, any random old $u$, and then *define* a right-hand side $f$ through the last equation.\n",
    "You've already picked what the true honest $u$ will be, so you know exactly what to compare the result of your numerical solver with.\n",
    "This is called the **method of manufactured solutions** and it's the standard way that the professionals make sure the PDE solvers do what they're supposed to.\n",
    "For complex problems, doing all that algebra by hand is tedious and error-prone, but with a symbolic algebra system it's super easy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finite element analysis\n",
    "\n",
    "Filling matrices and vectors for the finite element method is really error-prone.\n",
    "You can use a tool like SymPy to automate a lot of that.\n",
    "For example, let's say we wanted to solve the boundary value problem\n",
    "\n",
    "$$-u'' = f.$$\n",
    "\n",
    "Eventually we're going to have to evaluate a bunch of integrals that look like\n",
    "\n",
    "$$A_{ij} = \\int_0^1 \\phi_i'\\cdot \\phi_j'\\, dx$$\n",
    "\n",
    "where $\\phi_i$ are the shape functions.\n",
    "For linear shape functions this is easy to do by hand, but for quadratic or higher shape functions it starts to get awful.\n",
    "But SymPy can save us from doing all that!\n",
    "\n",
    "First, let's make a list of the quadratic basis functions on the interval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = sympy.symbols('x')\n",
    "shape_functions = [\n",
    "    (1 - x) * (1 - 2 * x),\n",
    "    4 * x * (1 - x),\n",
    "    x * (2 * x - 1)\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make sure we tabulated them right, we can plot them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "fig, axes = plt.subplots()\n",
    "X = np.linspace(0., 1., 51)\n",
    "for index, function in enumerate(shape_functions):\n",
    "    f = sympy.lambdify(x, function, modules='numpy')\n",
    "    axes.plot(X, f(X), label=index)\n",
    "axes.legend(loc='upper right')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can evaluate all the entries of the element stiffness matrix by iterating over the shape functions twice!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.zeros((3, 3))\n",
    "for i, ϕ_i in enumerate(shape_functions):\n",
    "    for j, ϕ_j in enumerate(shape_functions):\n",
    "        integrand = sympy.diff(ϕ_i) * sympy.diff(ϕ_j)\n",
    "        A[i, j] = sympy.integrate(integrand, (x, 0., 1.))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conclusion\n",
    "\n",
    "SymPy and symbolic algebra systems generally are really nice tools for automating away tedious and error-prone math.\n",
    "The example above shows how you might use symbolic algebra to automate one of the most important (and most error-prone) parts of finite element analysis: calculating local stiffness matrices.\n",
    "The software package FEniCS was the first to take this idea to its logical conclusion: scientists like you should be able to specify symbolically what problem you want to solve, and the software will automatically fill matrices and vectors for you.\n",
    "**Combining symbolic and numerical analysis is the biggest current trend in computational science today.**\n",
    "In the following examples, I'll use Firedrake, an offshoot of FEniCS, to show how that works."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "firedrake",
   "language": "python",
   "name": "firedrake"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}