{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Demo - Poisson's equation 1D\n",
    "=======================\n",
    "\n",
    "In this demo we will solve Poisson's equation\n",
    "    \n",
    "\\begin{align}\n",
    "\\label{eq:poisson}\n",
    "\\nabla^2 u(x) &= f(x), \\quad \\forall \\, x \\in [-1, 1]\\\\\n",
    "u(\\pm 1) &= 0, \n",
    "\\end{align}\n",
    "\n",
    "where $u(x)$ is the solution and $f(x)$ is some given function of $x$.\n",
    "\n",
    "We want to solve this equation with the spectral Galerkin method, using a basis composed of either Chebyshev $T_k(x)$ or Legendre $L_k(x)$ polynomials. Using $P_k$ to refer to either one, Shen's composite basis is then given as \n",
    "\n",
    "$$\n",
    "V^N = \\text{span}\\{P_k - P_{k+2}\\, | \\, k=0, 1, \\ldots, N-3\\},\n",
    "$$\n",
    "\n",
    "where all basis functions satisfy the homogeneous boundary conditions.\n",
    "\n",
    "For the spectral Galerkin method we will also need the weighted inner product\n",
    "\n",
    "$$\n",
    " (u, v)_w = \\int_{-1}^1 u v w \\, {dx},\n",
    "$$\n",
    "\n",
    "where $w(x)$ is a weight associated with the chosen basis, and $v$ and $u$ are test and trial functions, respectively. For Legendre the weight is simply $w(x)=1$, whereas for Chebyshev it is $w(x)=1/\\sqrt{1-x^2}$. Quadrature is used to approximate the integral\n",
    "\n",
    "$$\n",
    "\\int_{-1}^1 u v w \\, {dx} \\approx \\sum_{i=0}^{N-1} u(x_i) v(x_i) \\omega_i,\n",
    "$$\n",
    "\n",
    "where $\\{\\omega_i\\}_{i=0}^{N-1}$ are the quadrature weights associated with the chosen basis and quadrature rule. The associated quadrature points are denoted as $\\{x_i\\}_{i=0}^{N-1}$. For Chebyshev we can choose between `Chebyshev-Gauss` or `Chebyshev-Gauss-Lobatto`, whereas for Legendre the choices are `Legendre-Gauss` or `Legendre-Gauss-Lobatto`. \n",
    "\n",
    "With the weighted inner product in place we can create variational problems from the original PDE by multiplying with a test function $v$ and integrating over the domain. For a Legendre basis we can use integration by parts and formulate the variational problem: \n",
    "\n",
    "Find $u \\in V^N$ such that\n",
    "\n",
    "$$     (\\nabla u, \\nabla v) = -(f, v), \\quad \\forall \\, v \\in V^N.$$\n",
    "\n",
    "For a Chebyshev basis the integration by parts is complicated due to the non-constant weight and the variational problem used is instead: \n",
    "\n",
    "Find $u \\in V^N$ such that\n",
    "\n",
    "$$     (\\nabla^2 u, v)_w = (f, v)_w, \\quad \\forall \\, v \\in V^N.$$\n",
    "\n",
    "We now break the problem down to linear algebra. With any choice of basis or quadrature rule we use $\\phi_k(x)$ to represent the test function $v$ and thus\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "v(x) &= \\phi_k(x), \\\\\n",
    "u(x) &= \\sum_{j=0}^{N-3} \\hat{u}_j \\phi_j(x),\n",
    "\\end{align}\n",
    "$$\n",
    "where $\\hat{\\mathbf{u}}=\\{\\hat{u}_j\\}_{j=0}^{N-3}$ are the unknown expansion coefficients, also called the degrees of freedom.\n",
    "\n",
    "Insert into the variational problem for Legendre and we get the linear algebra system to solve for $\\hat{\\mathbf{u}}$\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "(\\nabla \\sum_{j=0}^{N-3} \\hat{u}_j \\phi_j, \\nabla \\phi_k) &= -(f, \\phi_k), \\\\\n",
    "\\sum_{j=0}^{N-3} \\underbrace{(\\nabla \\phi_j, \\nabla \\phi_k)}_{a_{kj}} \\hat{u}_j  &= -\\underbrace{(f, \\phi_k)}_{\\tilde{f}_k}, \\\\\n",
    "A \\hat{\\textbf{u}} &= -\\tilde{\\textbf{f}},\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where $A = (a_{kj})_{0 \\ge k, j \\ge N-3}$ is the stiffness matrix and $\\tilde{\\textbf{f}} = \\{\\tilde{f}_k\\}_{k=0}^{N-3}$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation with shenfun\n",
    "\n",
    "The given problem may be easily solved with a few lines of code using the shenfun Python module. The high-level code matches closely the mathematics and the stiffness matrix is assembled simply as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from shenfun import *\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "N = 100\n",
    "V = FunctionSpace(N, 'Legendre', quad='LG', bc=(0, 0))\n",
    "v = TestFunction(V)\n",
    "u = TrialFunction(V)\n",
    "A = inner(grad(u), grad(v))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using a manufactured solution that satisfies the boundary conditions we can create just about any corresponding right hand side $f(x)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy\n",
    "x = sympy.symbols('x')\n",
    "ue = (1-x**2)*(sympy.cos(4*x)*sympy.sin(6*x))\n",
    "fe = ue.diff(x, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `fe` is the right hand side that corresponds to the exact solution `ue`. We now want to use `fe` to compute a numerical solution $u$ that can be compared directly with the given `ue`. First, to compute the inner product $(f, v)$, we need to evaluate `fe` on the quadrature mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fl = sympy.lambdify(x, fe, 'numpy')\n",
    "ul = sympy.lambdify(x, ue, 'numpy')\n",
    "fj = Array(V, buffer=fl(V.mesh()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`fj` holds the analytical `fe` on the nodes of the quadrature mesh.\n",
    "Assemble right hand side $\\tilde{\\textbf{f}} = -(f, v)_w$ using the shenfun function `inner`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_tilde = inner(-fj, v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All that remains is to solve the linear algebra system \n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "A \\hat{\\textbf{u}} &= \\tilde{\\textbf{f}} \\\\\n",
    "\\hat{\\textbf{u}} &= A^{-1} \\tilde{\\textbf{f}} \n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u_hat = Function(V)\n",
    "u_hat = A/f_tilde\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get solution in real space, i.e., evaluate $u(x_i) = \\sum_{j=0}^{N-3} \\hat{u}_j \\phi_j(x_i)$ for all quadrature points $\\{x_i\\}_{i=0}^{N-1}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uj = u_hat.backward()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = V.mesh()\n",
    "plt.plot(X, uj)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}