09a_ODE_BVP.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "<table>\n",
    " <tr align=left><td><img align=left src=\"./images/CC-BY.png\">\n",
    " <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli</td>\n",
    "</table>\n",
    "\n",
    "Note:  The presentation below largely follows part I in \"Finite Difference Methods for Ordinary and Partial Differential Equations\" by LeVeque (SIAM, 2007)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "init_cell": true,
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "%matplotlib inline\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Solving Boundary Value Problems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The Problem\n",
    "\n",
    "We want to solve an ODE (PDE) that instead of having initial conditions is contained to an interval and has values at the edges of the interval.  This naturally comes about when we consider spatial problems.  One of the simplest cases for this is the Poisson problem in one-dimension\n",
    "\n",
    "$$\n",
    "    u_{xx} = f(x)\n",
    "$$\n",
    "\n",
    "where we will use the short-hand\n",
    "\n",
    "$$\n",
    "    u_{xx} = \\frac{\\text{d}^2 u}{\\text{d} x^2} \\quad \\text{or} \\quad \\frac{\\partial^2 u}{\\partial x^2}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Note that due to the order of the derivative we require two conditions to solve this.  The simplest case where we are on the domain $x \\in [a,b]$ is to have conditions such that we require $u(a) = u_a$ and $u(b) = u_b$ and are commonly termed boundary value problems (BVP).  If these conditions are both at one end of the domain then we can actually phrase the ODE (PDE) again as an initial value problem (IVP).  So what do we need to do to solve these types of problems?  We will consider several approaches to this problem:\n",
    "\n",
    "1. Recast our BVP to an IVP and use our standard methods for ODEs in a *shooting method*\n",
    "1. Use finite differences to represent the unknowns as a linear system and solve the resulting system.\n",
    "1. Use finite differences and Newton's method to solve non-linear BVP's\n",
    "1. Use Galerkin finite elements to find the best-fit piecewise linear function that satisfies the ODE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The Shooting Method\n",
    "\n",
    "The shooting method takes the approach that we want to use our ability to solve IVP problems and so tries to rewrite the problem as a root finding problem for the higher order initial condition that we are not given.  This is best illustrated by an example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Consider the problem\n",
    "$$ \n",
    "    u_{xx} + \\sin u = 0\n",
    "$$\n",
    "\n",
    "with\n",
    "$$\n",
    "    x \\in [0, 2] \\quad \\text{and} \\quad u(0) = 0.0, \\quad u(2.0) = \\frac{\\pi}{2}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can rewrite this problem as a system of two first order ODEs by defining\n",
    "\n",
    "$$\n",
    "    \\mathbf{v} = \\begin{bmatrix} u \\\\ u_x \\end{bmatrix} = \\begin{bmatrix} v_1 \\\\ v_2 \\end{bmatrix} \n",
    "$$\n",
    "\n",
    "then setting\n",
    "\n",
    "$$\n",
    "    \\frac{d}{dx}\\mathbf{v} = \\mathbf{f}(\\mathbf{v}) =  \\begin{bmatrix} v_2 \\\\ -\\sin v_1 \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We know that we want $v_1(0) = 0$ but what do we use for $v_2(0)$?  Making an initial guess at $v_2(0)$ and solving the associated IVP ODE we can then find out what these initial conditions produces on the right boundary of the problem.  Using a root-finding approach (or minimization routine) we can write this procedure as\n",
    "$$\n",
    "    \\min_{v_2(0)} \\left | \\pi / 2 - v_1(2) \\right |\n",
    "$$\n",
    "where the parameter we vary is $v_2(0)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Basic Shooting Method solving u_xx -sin(u) = 0\n",
    "from scipy.integrate import solve_ivp\n",
    "\n",
    "# Algorithm parameters\n",
    "TOLERANCE = 1e-8\n",
    "MAX_ITERATIONS = 100\n",
    "\n",
    "# Problem Statement\n",
    "a = 0.0\n",
    "b = 2.0\n",
    "N = 100\n",
    "x = numpy.linspace(a, b, N)\n",
    "u_a = 0.0\n",
    "u_b = numpy.pi / 2.0\n",
    "\n",
    "# RHS function\n",
    "def f(x, u):\n",
    "    return numpy.array([u[1], -numpy.sin(u[0])])\n",
    "\n",
    "# Initial guess\n",
    "# Slope at RHS\n",
    "u_prime_rhs = 1.0\n",
    "# Initial step size\n",
    "du_prime = 0.5\n",
    "\n",
    "# Plotting\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "\n",
    "# Main loop\n",
    "success = False\n",
    "u = numpy.empty((2, N))\n",
    "u0 = numpy.empty(2)\n",
    "convergence = numpy.zeros(MAX_ITERATIONS)\n",
    "for n in range(MAX_ITERATIONS):\n",
    "    \n",
    "    # Initial condition\n",
    "    u0 = [u_a, u_prime_rhs]\n",
    "\n",
    "    # Compute solution - note that we are only producing the intermediate values\n",
    "    # for demonstration purposes\n",
    "    \n",
    "    sol = solve_ivp(f, [a, b], u0, dense_output=True)\n",
    "    u_end = sol.y[0,-1]\n",
    "\n",
    "    # Stopping Criteria \n",
    "    # this should really be a bracketing scheme\n",
    "    convergence[n] = numpy.abs(u_end - u_b)\n",
    "    if convergence[n] < TOLERANCE:\n",
    "        success = True\n",
    "        break\n",
    "    else:\n",
    "        if u_end < u_b:\n",
    "            u_prime_rhs += du_prime\n",
    "        else:\n",
    "            u_prime_rhs -= du_prime\n",
    "        du_prime *= 0.5\n",
    "\n",
    "    axes.plot(x, sol.sol(x)[0], 'b')\n",
    "    axes.plot(b, u_b, 'ro')\n",
    "\n",
    "axes.set_title(\"Shooting Method Iterations\")\n",
    "axes.set_xlabel(\"$x$\")\n",
    "axes.set_ylabel(\"$u(x)$\")\n",
    "axes.grid()\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "n_range = numpy.arange(n)\n",
    "axes.semilogy(n_range, convergence[:n])\n",
    "axes.set_title(\"Convergence of Shooting Method\")\n",
    "axes.set_xlabel(\"step\")\n",
    "axes.set_ylabel(\"$|u(b) - U(b)|$\")\n",
    "axes.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The tricky part of this procedure is coming up with the search criteria, i.e. coming up with the decision of how to change $v_2(0)$ with respect to the position of $v_2(2)$ compared to what we want $u(2)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In general any minimization or rootfinding routine can be used in a shooting method.  These approaches are generally very effective at approaching non-linear BVPs where the next method we will discuss is too expensive to perform."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Exercise\n",
    "\n",
    "rewrite the previous example to combine solve_ivp with your favorite rootfinder e.g. root_scalar using brent to produce an efficient method for solving ODE BVP's using the shooting method.  As a start consider the interface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def shoot_bvp(f, x, u_a, u_b, i_a, i_b, rtol=1.e-3, atol=1.e-6):\n",
    "    \"\"\"\n",
    "    Solve the two-point boundary value problem on the interval x\\in [a,b], using a shooting method that combines \n",
    "        scipy.integrate.solve_ivp and scipy.optimize.root_scalar and allows a range of boundary conditions\n",
    "        \n",
    "        \n",
    "    parameters:\n",
    "    -----------\n",
    "    f: calleable \n",
    "        vector value function for righthand side of the ODE with interface f(t,u). returns ndarray of length 2\n",
    "    x: numpy array\n",
    "        coordinates array for solution  on interval [a,b] with x[0] = a, x[-1] = b\n",
    "    u_a:  numpy array (length 2)\n",
    "        provide initial boundary conditions  [u, u' ] at x=a\n",
    "    u_b:  numpy array (length 2)\n",
    "        target boundary condition at x = b\n",
    "    i_a: integer\n",
    "        index of known boundary condition at x = a.  i.e.\n",
    "        if dirichlet conditions : i_a = 0 and u(a) is known  \n",
    "        if neumann conditions   : i_a = 1 and u'(a) is known\n",
    "        the complementary index is adjusted to match the boundary condition at b\n",
    "    i_b: integer\n",
    "        index of known boundary condition at x = b. i.e.\n",
    "        if dirichlet conditions : i_b = 0 and u(b) is known  \n",
    "        if neumann conditions   : i_b = 1 and u'(b) is known\n",
    "        the complementary index is ignored at b\n",
    "    rtol:  float\n",
    "        relative tolerance\n",
    "    atol:  float\n",
    "        absolute tolerance\n",
    "        \n",
    "    returns:\n",
    "    --------\n",
    "    u: solution u(x) for x (uses dense_output from solve_ivp to interpolate solution onto x)\n",
    "    \"\"\"\n",
    "    \n",
    "    def udiff(u_var,u_b,ib):\n",
    "        \"\"\"\n",
    "            Solves the IVP problem using solve_ivp with adjustable initial condition u_var \n",
    "            and returns the difference between sol.y[ib] - u_b[ib]\n",
    "            \n",
    "        \"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Finite Difference methods for linear problems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Formulation\n",
    "\n",
    "The second approach we will consider involves the formation of a system of equations to solve based on finite difference approximations.  Again let's consider an example problem where\n",
    "\n",
    "$$\n",
    "    u_{xx} = f(x)\n",
    "$$\n",
    "\n",
    "with the boundary conditions $u(a) = u_a$ and $u(b) = u_b$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We know from our finite difference discussion that the second order centered difference approximation for the second derivative for a function $u(x)$ on a uniform mesh with mesh spacing $\\Delta x$ is\n",
    "$$\n",
    "    u_{xx} \\approx \\frac{u(x_{i-1}) - 2 u(x_i) + u(x_{i+1})}{\\Delta x^2}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If we discretize the domain of the original BVP into $N$ points (not including the boundaries) such that\n",
    "\n",
    "$$\n",
    "    x_i = a + \\frac{b - a}{N+1} \\cdot i ~~~ \\text{where} ~~~ i = 1, \\ldots, N\n",
    "$$\n",
    "\n",
    "we can then write the finite difference approximation as a system of linear equations!  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If for instance we take $N = 5$ then\n",
    "$$\\begin{aligned}\n",
    "    (U_{xx})_1 &\\approx \\frac{U_a - 2 U_1 + U_2}{\\Delta x^2} \\\\\n",
    "    (U_{xx})_2 &\\approx \\frac{U_1 - 2 U_2 + U_3}{\\Delta x^2} \\\\\n",
    "    (U_{xx})_3 &\\approx \\frac{U_2 - 2 U_3 + U_4}{\\Delta x^2} \\\\\n",
    "    (U_{xx})_4 &\\approx \\frac{U_3 - 2 U_4 + U_5}{\\Delta x^2} \\\\\n",
    "    (U_{xx})_5 &\\approx \\frac{U_4 - 2 U_5 + U_b}{\\Delta x^2} \\\\\n",
    "\\end{aligned}$$\n",
    "where we have used $U_a = u(a)$ and $U_b = u(b)$ as the boundary conditions.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Using these approximations to the derivatives we can then write the ODE as\n",
    "$$\n",
    "    \\frac{1}{\\Delta x^2}\\begin{bmatrix}\n",
    "    -2 &  1 &    &    &    \\\\\n",
    "     1 & -2 &  1 &    &    \\\\\n",
    "       &  1 & -2 &  1 &    \\\\\n",
    "       &    &  1 & -2 &  1 \\\\\n",
    "       &    &    &  1 & -2 \\\\\n",
    "    \\end{bmatrix} \\begin{bmatrix}\n",
    "        U_1 \\\\ U_2 \\\\ U_3 \\\\ U_4 \\\\ U_5\n",
    "    \\end{bmatrix} = \n",
    "    \\begin{bmatrix}\n",
    "        f(x_1) \\\\ f(x_2) \\\\ f(x_3) \\\\ f(x_4) \\\\ f(x_5) \\\\\n",
    "    \\end{bmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Note that our previous example used for the shooting method is difficult in the current context as the unknown function $u$ is in the function $f$ so that we would need to actual solve a non-linear system of equations.  This is still possible in this context using an approach such as a Newton solver.  We will address that problem after we consider the linear problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Boundary Conditions\n",
    "\n",
    "This does not include the boundary conditions though.  We can add these values easily for Dirichlet boundary conditions by sending the values we know to the $b$ vector:\n",
    "$$\\begin{aligned}\n",
    "    \\frac{U_a - 2 U_1 + U_2}{\\Delta x^2} = f(x_1) &\\Rightarrow& \\frac{- 2 U_1 + U_2}{\\Delta x^2} = f(x_1) - \\frac{U_a}{\\Delta x^2} \\\\\n",
    "    \\frac{U_4 - 2 U_5 + U_b}{\\Delta x^2} = f(x_5) &\\Rightarrow& \\frac{U_4 - 2 U_5}{\\Delta x^2} = f(x_5) - \\frac{U_b}{\\Delta x^2}\n",
    "\\end{aligned}$$\n",
    "so that final system looks like\n",
    "$$\n",
    "    \\frac{1}{\\Delta x^2} \\begin{bmatrix}\n",
    "    -2 &  1 &    &    &    \\\\\n",
    "     1 & -2 &  1 &    &    \\\\\n",
    "       &  1 & -2 &  1 &    \\\\\n",
    "       &    &  1 & -2 &  1 \\\\\n",
    "       &    &    &  1 & -2 \\\\\n",
    "    \\end{bmatrix} \\begin{bmatrix}\n",
    "        U_1 \\\\ U_2 \\\\ U_3 \\\\ U_4 \\\\ U_5\n",
    "    \\end{bmatrix} = \n",
    "    \\begin{bmatrix}\n",
    "        f(x_1) - \\frac{U_a}{\\Delta x^2} \\\\ f(x_2) \\\\ f(x_3) \\\\ f(x_4) \\\\ f(x_5) - \\frac{U_b}{\\Delta x^2} \\\\\n",
    "    \\end{bmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example\n",
    "\n",
    "Here we will  solve the BVP\n",
    "$$\n",
    "    u_{xx} = e^x, \\quad x \\in [0, 1] \\quad \\text{with} \\quad u(0) = 0.0, \\text{ and } u(1) = 3\n",
    "$$\n",
    "via the construction of a linear system of equations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This problem can actually be solved relatively easily analytically by simply integrating twice which introduces two constants of integration\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "    u_{xx} &= e^x \\\\\n",
    "    u_x &= A + e^x \\\\\n",
    "    u &= Ax + B + e^x\\\\\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "Note: this solution is just two indefinite integrals of the right-hand-side $f(x)$ plus an arbitrary line $Ax + B$ which is in the kernel of the operator $d^2/dx^2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The specific solution is then just found by substituting in the boundary conditions to solve for $A$ and $B$\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "u(0) &= B + 1 = 0 \\Rightarrow B = -1 \\\\\n",
    "u(1) &= A - 1 + e^{1} = 3 \\Rightarrow A = 4 - e\\\\ \n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "With solution\n",
    "\n",
    "$$\n",
    "    u(x) = (4 - e) x - 1 + e^x\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Problem setup\n",
    "a = 0.0\n",
    "b = 1.0\n",
    "u_a = 0.0\n",
    "u_b = 3.0\n",
    "f = lambda x: numpy.exp(x)\n",
    "u_true = lambda x: (4.0 - numpy.exp(1.0)) * x - 1.0 + numpy.exp(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Discretization\n",
    "N = 10\n",
    "x_bc = numpy.linspace(a, b, N + 2)\n",
    "x = x_bc[1:-1]\n",
    "delta_x = (b - a) / (N + 1)\n",
    "\n",
    "# Construct matrix A\n",
    "A = numpy.zeros((N, N))\n",
    "diagonal = numpy.ones(N) / delta_x**2\n",
    "A += numpy.diag(diagonal * -2.0, 0)\n",
    "A += numpy.diag(diagonal[:-1], 1)\n",
    "A += numpy.diag(diagonal[:-1], -1)\n",
    "\n",
    "# Construct RHS with Dirichlet BC's\n",
    "b = f(x)\n",
    "b[0] -= u_a / delta_x**2\n",
    "b[-1] -= u_b / delta_x**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Solve system\n",
    "U = numpy.empty(N + 2)\n",
    "# insert Boundary conditions\n",
    "\n",
    "U[0] = u_a\n",
    "U[-1] = u_b\n",
    "\n",
    "#solve\n",
    "U[1:-1] = numpy.linalg.solve(A, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Plot result\n",
    "fig = plt.figure(figsize=(16,6))\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "axes.plot(x_bc, U, 'o', label=\"Computed\")\n",
    "axes.plot(x_bc, u_true(x_bc), 'k', label=\"True\")\n",
    "axes.set_title(\"Solution to $u_{xx} = e^x$\")\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"u(x)\")\n",
    "axes.grid()\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "err = numpy.abs(U - u_true(x_bc))\n",
    "axes.plot(x_bc, err, 'ro-')\n",
    "axes.set_title(\"Absolute Error\")\n",
    "axes.set_xlabel('x')\n",
    "axes.set_ylabel('err')\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If we instead have Neumann boundary conditions it is no longer clear how to handle the boundary conditions using the above approach.  Instead a **ghost cell** approach is often used.  These **ghost cells** are added unknowns that represent the boundary values that we actually know.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For instance, if we had the BVP\n",
    "$$\n",
    "    u_{xx} = e^x, \\quad x \\in [-1, 1] \\quad \\text{with} \\quad u(-1) = 3, \\text{ and } u_x(1) = -5\n",
    "$$\n",
    "then we could keep the boundary values in the vector of unknowns so that now\n",
    "$$\n",
    "    U = \\begin{bmatrix} U_0 \\\\ U_1 \\\\ \\vdots \\\\ U_N \\\\ U_{N+1} \\end{bmatrix}\n",
    "$$\n",
    "where here $U_0$ and $U_{N+1}$ are actually the boundary points.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The matrix $A$ is then modified to have the appropriate relationships.  In this case the left boundary condition leads to\n",
    "$$\n",
    "    A = \\begin{bmatrix}\n",
    "  1 &    &    &    &    &    \\\\\n",
    "  \\frac{1}{\\Delta x^2} & \\frac{-2}{\\Delta x^2} &  \\frac{1}{\\Delta x^2} &    &    &    \\\\\n",
    "    &  \\frac{1}{\\Delta x^2} & \\frac{-2}{\\Delta x^2} &  \\frac{1}{\\Delta x^2} &    &    \\\\\n",
    "    & & \\ddots & \\ddots & \\ddots\n",
    "    \\end{bmatrix} \\quad \\text{and} \\quad b = \\begin{bmatrix}\n",
    "        u(a) \\\\ f(x_1) \\\\ f(x_2) \\\\ \\vdots\n",
    "    \\end{bmatrix}\n",
    "$$\n",
    "which multiplied out simply gives\n",
    "$$\n",
    "    U_0 = u(-1) = 3.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For the right boundary condition we can use the second order backward finite difference approximation for the first derivative\n",
    "$$\n",
    "    u_x(b) \\approx \\frac{3 U_{N+1} - 4 U_{N} + U_{N - 1}}{2.0 \\Delta x} = -5\n",
    "$$\n",
    "which can be incorporated into the matrix $A$ and vector $b$ as\n",
    "$$\n",
    "    A =  \\begin{bmatrix}\n",
    "     \\ddots & \\ddots & \\ddots &    &    \\\\\n",
    "            & \\frac{1}{\\Delta x^2} &     \\frac{-2}{\\Delta x^2}&  \\frac{1}{\\Delta x^2} &    \\\\\n",
    "            &        &      \\frac{1}{\\Delta x^2} & \\frac{-2}{\\Delta x^2} &  \\frac{1}{\\Delta x^2} \\\\\n",
    "            &        &      \\frac{1}{2 \\Delta x} &  \\frac{-4}{2 \\Delta x} &  \\frac{3}{2 \\Delta x} \\\\\n",
    "    \\end{bmatrix} ~~~~ \\text{and} ~~~~ b = \\begin{bmatrix}\n",
    "        \\vdots \\\\ f(x_N) \\\\ u_x(b)\n",
    "    \\end{bmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "All together the new system looks like\n",
    "$$\n",
    "    \\begin{bmatrix}\n",
    "     1 &    &    &    &    &    \\\\\n",
    "     \\frac{1}{\\Delta x^2} & \\frac{-2}{\\Delta x^2} &  \\frac{1}{\\Delta x^2} &    &    &    \\\\\n",
    "       &  \\ddots & \\ddots &  \\ddots &    \\\\\n",
    "       &    & \\frac{1}{\\Delta x^2} & \\frac{-2}{\\Delta x^2} &  \\frac{1}{\\Delta x^2} \\\\\n",
    "            &        &      \\frac{1}{2 \\Delta x} &  \\frac{-4}{2 \\Delta x} &  \\frac{3}{2 \\Delta x} \\\\\n",
    "    \\end{bmatrix} \\begin{bmatrix}\n",
    "        U_0 \\\\ U_1 \\\\ \\vdots \\\\ U_N \\\\ U_{N+1}\n",
    "    \\end{bmatrix} = \n",
    "    \\begin{bmatrix}\n",
    "        u(a) \\\\ f(x_1) \\\\ \\vdots \\\\ f(x_N) \\\\ u_x(b)\n",
    "    \\end{bmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example\n",
    "\n",
    "Want to solve the BVP\n",
    "$$\n",
    "    u_{xx} = e^x, \\quad x \\in [-1, 1] \\quad \\text{with} \\quad u(-1) = 3.0, \\text{ and } u_x(1) = -5.0\n",
    "$$\n",
    "via the construction of a linear system of equations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "General Solution\n",
    "$$\n",
    "    u(x) = A x + B + e^x \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Apply Boundary conditions\n",
    "\\begin{align*}\n",
    "    u(-1) &= -A + B + e^{-1} = 3\\\\\n",
    "    u_x(1) &= A + e^1 = -5\\\\\n",
    "\\end{align*}\n",
    "\n",
    "Which implies\n",
    "\\begin{align*}\n",
    "    A &= -5 - e \\\\\n",
    "    B & = A - e^{-1} + 3\\\\\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Which gives the Specific solution\n",
    "\n",
    "$$\n",
    "    u(x) = -(5 + e) x -(2 + e + e^{-1}) + e^{x}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Problem setup\n",
    "a = -1.0\n",
    "b = 1.0\n",
    "u_a = 3.0 \n",
    "u_x_b = -5.0\n",
    "f = lambda x: numpy.exp(x)\n",
    "u_true = lambda x: -(5.0 + numpy.exp(1.0)) * x - (2.0 + numpy.exp(1.0) + numpy.exp(-1.0)) + numpy.exp(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Discretization\n",
    "N = 20\n",
    "x_bc = numpy.linspace(a, b, N + 2)\n",
    "x = x_bc[1:-1]\n",
    "delta_x = (b - a) / (N + 1)\n",
    "\n",
    "# Construct matrix A\n",
    "A = numpy.zeros((N + 2, N + 2))\n",
    "diagonal = numpy.ones(N + 2) / delta_x**2\n",
    "A += numpy.diag(diagonal * -2.0, 0)\n",
    "A += numpy.diag(diagonal[:-1], 1)\n",
    "A += numpy.diag(diagonal[:-1], -1)\n",
    "\n",
    "# Construct RHS\n",
    "b = f(x_bc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Boundary conditions\n",
    "A[0, 0] = 1.0\n",
    "A[0, 1] = 0.0\n",
    "A[-1, -1] = 3.0 / (2.0 * delta_x)\n",
    "A[-1, -2] = -4.0 / (2.0 * delta_x)\n",
    "A[-1, -3] = 1.0 / (2.0 * delta_x)\n",
    "\n",
    "b[0] = u_a\n",
    "b[-1] = u_x_b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Solve system\n",
    "U = numpy.empty(N + 2)\n",
    "U = numpy.linalg.solve(A, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Plot result\n",
    "fig = plt.figure(figsize=(16,6))\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "axes.plot(x_bc, U, 'o', label=\"Computed\")\n",
    "axes.plot(x_bc, u_true(x_bc), 'k', label=\"True\")\n",
    "axes.set_title(\"Solution to $u_{xx} = e^x$\")\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"u(x)\")\n",
    "axes.grid()\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "err = numpy.abs(U - u_true(x_bc))\n",
    "axes.plot(x_bc, err, 'ro-')\n",
    "axes.set_title(\"Absolute Error\")\n",
    "axes.set_xlabel('x')\n",
    "axes.set_ylabel('err')\n",
    "axes.grid()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Generalizations and Sparse Matrices for $A u = f$\n",
    "\n",
    "Most finite difference approximations to linear two-point boundary value problems will  lead to large sparse linear systems of the form $A \\mathbf{u} = \\mathbf{f}$ where $\\mathbf{u}$ is a discrete sampling of the continous function $u$ at a vector of coordinates $\\mathbf{x}$ and $\\mathbf{f} = f(\\mathbf{x})$ is the discrete sampling of the RHS, and $A$ is a finite difference approximation to some Linear differential operator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Using the tools we have already developed for general finite difference stencils (e.g. with the routine fdcoeffV), we can write some useful utility routines for assembly of general sparse operators.  Here we will use a combination of `fdcoeffV` and `scipy.sparse` to assemble general sparse matrices for the second derivative operator on an arbitrary mesh with a choice of consistent dirichlet or Neumann BC's"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from fdcoeffV import fdcoeffV\n",
    "from scipy.sparse import lil_matrix, identity\n",
    "from scipy.sparse.linalg import spsolve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def D2(x,bcs=['dirichlet', 'dirichlet']):\n",
    "    \"\"\"  \n",
    "        Assemble a general sparse second-order finite-difference approximation to d/dx^2 on a possibly irregular mesh\n",
    "        First and last rows are set by string bcs\n",
    "        \n",
    "        parameters:\n",
    "        -----------\n",
    "        x: numpy.array\n",
    "            mesh coordinates\n",
    "        bcs: list of strings for boundary conditions e.g [left_string, right_string] where\n",
    "            the strings can be either 'dirichlet' or 'neumann'\n",
    "    \"\"\"\n",
    "    N = len(x)\n",
    "    A = lil_matrix((N,N))\n",
    "    if bcs[0] == 'dirichlet':\n",
    "        A[0,0] = 1.\n",
    "    elif bcs[0] == 'neumann':\n",
    "        A[0,0:3] = fdcoeffV(1,x[0],x[:3])\n",
    "    else:\n",
    "        raise ValueError('no known BC type for left boundary {}'.format(bcs[0]))\n",
    "        \n",
    "    if bcs[1] == 'dirichlet':\n",
    "        A[-1,-1] = 1.\n",
    "    elif bcs[1] == 'neumann':\n",
    "        A[-1,-3:] = fdcoeffV(1,x[-1],x[-3:])\n",
    "    else:\n",
    "        raise ValueError('no known BC type for right boundary {}'.format(bcs[1]))\n",
    "        \n",
    "    for i in range(1,N-1):\n",
    "        A[i, i-1:i+2] = fdcoeffV(2, x[i], x[i-1:i+2] )        \n",
    "    return A.tocsr()\n",
    "\n",
    "def RHS(x, f, bvalues):\n",
    "    \"\"\" Set the rhs vector\n",
    "    \n",
    "        parameters\n",
    "        ----------\n",
    "        x: numpy.array\n",
    "            mesh coordinates\n",
    "        f: callable\n",
    "            rhs function for interior points called on f(x[1:-2])\n",
    "        bvalues:  numpy.array (len 2)\n",
    "            values for boundary conditions (either dirichlet or neumann)        \n",
    "    \"\"\"\n",
    "    \n",
    "    N = len(x)\n",
    "    rhs = numpy.empty(N)\n",
    "    rhs[[0, N-1]] = bvalues\n",
    "    rhs[1:-1] = f(x[1:-1])\n",
    "    \n",
    "    return rhs\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# check our routines for a small mesh\n",
    "x = numpy.linspace(0,1,5)\n",
    "print('x = {}'.format(x))\n",
    "print()\n",
    "A = D2(x,['dirichlet','neumann'])\n",
    "print('A = \\n{}'.format(A.todense()))\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6,6))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.spy(A)\n",
    "axes.grid()\n",
    "axes.set_title('Visualization of Sparse Matrices')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "b = RHS(x, numpy.exp, [1.,1.])\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Now let's solve our previous problem again\n",
    "\n",
    "$$\n",
    "    u_{xx} = e^x \\quad u(0) = 0.,\\quad u(1) = 3.0\n",
    "$$\n",
    "\n",
    "on a regular mesh with 100 points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10\n",
    "a = 0.\n",
    "b = 1.\n",
    "x = numpy.linspace(a, b, N)\n",
    "\n",
    "A = D2(x)\n",
    "f = RHS(x, numpy.exp, [0., 3. ])\n",
    "\n",
    "u = spsolve(A,f)\n",
    "u_true = lambda x: (4.0 - numpy.exp(1.0)) * x - 1.0 + numpy.exp(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Plot result\n",
    "fig = plt.figure(figsize=(16,6))\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "axes.plot(x, u, 'o', label=\"Computed\")\n",
    "axes.plot(x, u_true(x), 'k', label=\"True\")\n",
    "axes.set_title(\"Solution to $u_{xx} = e^x$\")\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"u(x)\")\n",
    "axes.grid()\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "err = numpy.abs(u - u_true(x))\n",
    "axes.plot(x, err, 'ro-')\n",
    "axes.set_title(\"Absolute Error\")\n",
    "axes.set_xlabel('x')\n",
    "axes.set_ylabel('err')\n",
    "axes.grid()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### And Check convergence\n",
    "\n",
    "From our discussion of finite differences, our current discrete operator should have errors that scale as $O(\\Delta x^2)$.  We can easily check that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "N = [ 2**n for n in range(3,9) ]\n",
    "rel_err = numpy.empty(len(N))\n",
    "delta_x = numpy.empty(len(N))\n",
    "\n",
    "for i,n in enumerate(N):\n",
    "    x = numpy.linspace(a, b, n)\n",
    "    A = D2(x)\n",
    "    f = RHS(x, numpy.exp, [0., 3. ])\n",
    "    u = spsolve(A,f)\n",
    "    rel_err[i] = numpy.linalg.norm(u - u_true(x))/numpy.linalg.norm(u_true(x))\n",
    "    delta_x[i] = x[1] - x[0]\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# calculate best fit slope and Plot result\n",
    "p = numpy.polyfit(numpy.log(delta_x),numpy.log(rel_err),1)\n",
    "dx = numpy.logspace(numpy.log10(delta_x[0]),numpy.log10(delta_x[-1]),100)\n",
    "err = numpy.exp(p[1])*dx**p[0]\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.loglog(delta_x, rel_err, 'o',label='error', markersize=10 )\n",
    "axes.loglog(dx,err,'r--',label=\"p={}\".format(p[0]))\n",
    "axes.set_title(\"Convergence\")\n",
    "axes.set_xlabel(\"$\\Delta x$\")\n",
    "axes.set_ylabel(\"$||u - u_{true}||/||u_{true}||$\")\n",
    "axes.legend(loc='best')\n",
    "axes.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### An exercise:\n",
    "\n",
    "Use these tools to try and  solve a more interesting two-point linear  BVP\n",
    "\n",
    "$$\n",
    "    u_{xx} + 4u = x , \\quad\\quad u(0) = 0,  u(2\\pi) = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Finite Difference methods for non-linear problems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "Suppose we wanted to solve our original non-linear 2-pt BVP \n",
    "$$\n",
    "    u_{xx} + \\sin{u} = 0 , \\quad\\quad u(0) = 0,  u(2) = \\frac{\\pi}{2}\n",
    "$$\n",
    "\n",
    "Using finite difference methods...How would you go about setting up and solving this problem?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If we first replace our continuous function $u$ with its discrete approximation \n",
    "\n",
    "$$\n",
    "    \\mathbf{u} = u(\\mathbf{x}) \n",
    "$$ \n",
    "\n",
    "where $\\mathbf{x}$ is a vector of coordinates in $\\mathbb{R}^n$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "we can approximate the second derivative operator using finite differences\n",
    "$$\n",
    "    u_{xx} \\approx D\\mathbf{u}\n",
    "$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "then our initial ODE becomes the discrete **Non-linear** problem\n",
    "\n",
    "$$\n",
    "    D\\mathbf{u} + \\sin(\\mathbf{u}) = \\mathbf{0} \\quad u_0 = 0,\\, u_N = \\frac{\\pi}{2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or more succinctly\n",
    "\n",
    "$$\n",
    "    \\mathbf{F}(\\mathbf{u}) = \\mathbf{0}\n",
    "$$\n",
    "\n",
    "which can be attacked with a Newton solver, Given an appropriate Jacobian. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We have already discussed that for a linear problem with\n",
    "\n",
    "$$ \n",
    "    \\mathbf{F}(\\mathbf{u}) =  A\\mathbf{u} -\\mathbf{b} \n",
    "$$\n",
    "\n",
    " each equation in the residual can be written as\n",
    "\n",
    "$$\n",
    "    F_i(\\mathbf{u}) = \\sum_j a_{ij}u_j - b_i\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Therefore the Jacobian is simply\n",
    "$$\n",
    "    J_{ij} = \\frac{\\partial F_i}{\\partial u_j} =a_{ij}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or\n",
    "\n",
    "$$\n",
    "    J = A\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For our non-linear residual\n",
    "\n",
    "$$\n",
    "  \\mathbf{F}(\\mathbf{u}) =   D\\mathbf{u} + \\sin(\\mathbf{u}) = \\mathbf{0} \n",
    "$$\n",
    "\n",
    "each equation can be written\n",
    "\n",
    "$$\n",
    "    F_i(\\mathbf{u}) = \\sum_j D_{ij}u_j  +\\sin(u_i) = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or\n",
    "$$\n",
    "    F_i(\\mathbf{u}) = \\sum_j \\left[ D_{ij}u_j  +\\delta_{ij}\\sin(u_j)\\right] = 0\n",
    "$$\n",
    "\n",
    "where $\\delta_{ij}$ is the kronecker delta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "which suggests that the Jacobian is just\n",
    "\n",
    "$$\n",
    "    J_{ij}(\\mathbf{u})=  D_{ij}  +\\delta_{ij}\\cos(u_j) = 0\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "$$\n",
    "    J(\\mathbf{u}) = D +\\mathrm{diag}(\\cos(\\mathbf{u}))\n",
    "$$\n",
    "\n",
    "which we should be able to hand off to a Newton solver. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Question?\n",
    "\n",
    "How to handle Boundary conditions?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Let's code this up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We'll need to describe a sparse diagonal matrix and a sparse solver\n",
    "from scipy.sparse import spdiags\n",
    "from scipy.sparse.linalg import spsolve\n",
    "\n",
    "# and a Newton Solver for sparse matrices\n",
    "def newton(F,J,x0,tol=1.e-6,MAX_ITS=100,verbose=True):\n",
    "    \"\"\" Solve F(x) = 0 using Newton's method until ||F(x)|| < tol or reaches Max iterations\n",
    "    \n",
    "    Params:\n",
    "    -------\n",
    "        F: calleable: \n",
    "            Function returning residual F(x)\n",
    "        J: calleable\n",
    "            Function returning Jacobian J(x)\n",
    "        tol: float\n",
    "            stopping criteria for ||F|| < tol\n",
    "        MAX_ITS: integer\n",
    "            maximum number of iterations\n",
    "        verbose: bool\n",
    "            if true spits out norm of the residual at each iteration\n",
    "            \n",
    "    Returns:\n",
    "    --------\n",
    "        x: list\n",
    "            list of points for each Newton iteration (just used for plotting intermediate results)\n",
    "            the solution is x[-1]\n",
    "    Raises:\n",
    "    -----------\n",
    "        ValueError if number of its exceeds MAX_ITS\n",
    "        \n",
    "    \"\"\"\n",
    "    x = [ x0 ]\n",
    "    for k in range(MAX_ITS+1):\n",
    "        xk = x[k]\n",
    "        res = numpy.linalg.norm(F(xk))\n",
    "        if verbose:\n",
    "            print('k = {}, ||F|| = {}'.format(k,res))\n",
    "        if res < tol: \n",
    "            return numpy.array(x), k\n",
    "        delta = spsolve(J(xk),-F(xk))\n",
    "        x.append( xk + delta)\n",
    "        \n",
    "        \n",
    "                \n",
    "    raise ValueError('Maximum number of iterations exceeded {}'.format(MAX_ITS))\n",
    "    \n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# test on a linear problem\n",
    "\n",
    "N = 20\n",
    "x = numpy.linspace(0.,1.,N)\n",
    "A = D2(x)\n",
    "# set the residual with 0 dirichlet condtions\n",
    "def F(u):\n",
    "    f =  A.dot(u) - numpy.exp(x)\n",
    "    # set residual to 0 on dirichlet BCs\n",
    "    f[[0,-1]] = 0\n",
    "    return f\n",
    "\n",
    "# set the jacobian\n",
    "def J(u):\n",
    "    return A\n",
    "    \n",
    "\n",
    "# initial guess use a straight line with the correct boundary conditions\n",
    "\n",
    "u0 = numpy.linspace(0.,3.,N)\n",
    "sol,its = newton(F,J,u0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "u_true = lambda x: (4.0 - numpy.exp(1.0)) * x - 1.0 + numpy.exp(x)\n",
    "\n",
    "fig = plt.figure(figsize=(16,6))\n",
    "axes = fig.add_subplot(1,2,1)\n",
    "for k in range(its+1):\n",
    "    axes.plot(x,sol[k],label='$k={}$'.format(k))\n",
    "axes.plot(x,u_true(x),'k--',label='$u_{true}$')\n",
    "axes.grid()\n",
    "axes.legend(loc='best')\n",
    "axes.set_xlabel('x')\n",
    "axes.set_ylabel('u')\n",
    "\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "err = numpy.abs(sol[-1] - u_true(x))\n",
    "axes.plot(x, err, 'ro-')\n",
    "axes.set_title(\"Absolute Error\")\n",
    "axes.set_xlabel('x')\n",
    "axes.set_ylabel('err')\n",
    "axes.grid()\n",
    "plt.show()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# test on a non-linear problem\n",
    "\n",
    "N = 20\n",
    "x = numpy.linspace(0.,2.,N)\n",
    "A = D2(x)\n",
    "# set the residual with 0 dirichlet condtions\n",
    "def F(u):\n",
    "    f =  A.dot(u) + numpy.sin(u)\n",
    "    # set residual to 0 on dirichlet BCs\n",
    "    f[[0,-1]] = 0\n",
    "    return f\n",
    "\n",
    "# set the jacobian\n",
    "def J(u):\n",
    "    n = len(u)\n",
    "    return A + spdiags(numpy.cos(u),0,n,n)\n",
    "    \n",
    "\n",
    "# initial guess use a straight line with the correct boundary conditions\n",
    "\n",
    "u0 = numpy.linspace(0.,numpy.pi/2,N)\n",
    "sol,its = newton(F,J,u0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "for k in range(its+1):\n",
    "    axes.plot(x,sol[k],label='$k={}$'.format(k))\n",
    "axes.grid()\n",
    "axes.legend(loc='best')\n",
    "axes.set_xlabel('x')\n",
    "axes.set_ylabel('u')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Introduction to Finite Element methods\n",
    "\n",
    "In the previous problem we used Finite Difference methods to transform linear and non-linear continous ODE's into discrete linear or non-linear problems of form $A\\mathbf{u} = \\mathbf{f}$ or $\\mathbf{F}(\\mathbf{u}) = \\mathbf{0}$ which we can solve with standard linear and non-linear solvers.\n",
    "\n",
    "Finite Element Methods (FEM),  will, in the end,  accomplish the same thing,  mapping from a continous, infinite dimensional problem to a finite dimensional problem, but in a mathematically more satisfying approach"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Outline\n",
    "\n",
    "* Finite Dimensional Function Spaces (piecewise linear functions)\n",
    "* Bases for Function spaces (the hat functions)\n",
    "* Global Basis vs Element Basis\n",
    "* Galerkin Projection as a least-squares problem\n",
    "    * Review of Linear Least Squares\n",
    "    * The Projection problem\n",
    "    * The Mass Matrix and Load vector\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Piecewise Polynomial Function Spaces\n",
    "\n",
    "In our discussion of interpolation, we introduced the idea of a piecewise polynomial interpolation and illustrated it with a **Piecewise Linear**  \"connect-the-dots\" function defined on a set of nodes on the interval $x\\in[0,L]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "nodes = numpy.array([ 0., 2.,  4.,  6.,  7.,  9., 10.        ])\n",
    "f = numpy.array( [0.5 ,  1.5, 2.0, 3., 2.75 , 1.,  0.5])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "L=10\n",
    "#x = numpy.linspace(0,L,N)\n",
    "#x[1:-1] += numpy.random.rand(N-2)\n",
    "#f = .1*(x*(L - x) + numpy.random.rand(N))\n",
    "#nodes = numpy.array([ 0., 2.07005551,  4.30576527,  5.70618816,  6.72033573,  9.30961848, 10.        ])\n",
    "#f = numpy.array( [0.5595035,  1.72355142, 2.52531147, 2.49691638, 2.24858267, 0.73643982,  0.02774723])\n",
    "N = len(nodes)\n",
    "\n",
    "fig = plt.figure(figsize=(10,5))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.plot(nodes,f,'bo-')\n",
    "axes.set_xlabel('x')\n",
    "axes.set_title('Piecewise linear function with {} Nodes'.format(N))\n",
    "axes.grid()\n",
    "\n",
    "#%precision 3\n",
    "#print('nodes  = {}'.format(nodes))\n",
    "#print('values = {}'.format(f))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### Points: \n",
    "\n",
    "* The function is continuous $f(x)\\in C^0$ and defined for all $x\\in[0,L]$\n",
    "* However, it is defined by only 7 points\n",
    "* **All** piecewise linear functions on this mesh form a **Function Space** that is isomorphic to  $\\mathbb{R}^7$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true
   },
   "outputs": [],
   "source": [
    "# define the hat functions \n",
    "def hats(x, nodes):\n",
    "    \"\"\" mainly for pedagogy, calculate the hat functions given a set of nodes \"\"\"\n",
    "    N = len(nodes)\n",
    "    l0 = lambda x, x0, x1: (x - x1)/(x0 - x1)\n",
    "    l1 = lambda x, x0, x1: (x - x0)/(x1 - x0)\n",
    "\n",
    "    phi = numpy.zeros((N, len(x)))\n",
    "    conds = [x <= nodes[1], x >= nodes[1]]\n",
    "    phi[0] = numpy.piecewise(x, conds, [ lambda x: l0(x, nodes[0], nodes[1]),  lambda x: 0.])\n",
    "    for i in range(1,N-1):\n",
    "        conds = [x < nodes[i-1], (x >= nodes[i-1]) & (x <= nodes[i]), (x >= nodes[i]) & (x <= nodes[i+1]), x >= nodes[i+1]]\n",
    "        funcs = [ lambda x: 0., lambda x: l1(x,nodes[i-1],nodes[i]), lambda x: l0(x,nodes[i],nodes[i+1]), lambda x: 0.]\n",
    "        phi[i] = numpy.piecewise(x,conds,funcs)\n",
    "                       \n",
    "    conds = [x <= nodes[-2], x >= nodes[-2]]\n",
    "    phi[-1] = numpy.piecewise(x, conds, [   lambda x: 0., lambda x: l1(x, nodes[-2], nodes[-1])])\n",
    "    return phi\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The Hat functions:  A basis for piecewise linear functions\n",
    "\n",
    "As any vector in $\\mathbb{R}^7$ can  be describe uniquely given 7 linearly independent **basis vectors**,  any piecewise-linear function in this function space can be described uniquely  given a linearly independent set of **basis functions**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "A useful basis for a piecewise linear function are the global \"hat\" functions which are piecewise linear functions defined such that $\\phi_i =1$ at node $i$ and is 0 at all other nodes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "This figure shows two hat functions, $\\phi_2$ and $\\phi_3$ which control the functions at points 2 and 3 respectively"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0.,L,1000)\n",
    "phi = hats(x,nodes)\n",
    "\n",
    "fig = plt.figure(figsize=(10,5))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.plot(nodes,f,'bo-')\n",
    "axes.plot(x,phi[2],'r',label='$\\phi_2$')\n",
    "axes.plot(x,phi[3],'c--',label='$\\phi_3$')\n",
    "#axes.plot(x,phi[5],'g--',label='$\\phi_5$')\n",
    "axes.set_xlabel('x')\n",
    "axes.legend()\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "All 7 Hat Functions for this particular mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0.,L,1000)\n",
    "phi = hats(x,nodes)\n",
    "\n",
    "fig = plt.figure(figsize=(10,5))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.plot(nodes,f,'bo-')\n",
    "for i in range(phi.shape[0]):\n",
    "   axes.plot(x,phi[i],label='$\\phi_{}$'.format(i))\n",
    "axes.set_xlabel('x')\n",
    "axes.legend()\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Given a basis for the space of all piecewise linear functions on a given mesh.  We can now write any function in this space as a linear combination of the basis functions\n",
    "i.e. \n",
    "$$\n",
    "    \\tilde{f}(x) = \\sum_j^N w_j \\phi_j(x)\n",
    "$$\n",
    "\n",
    "where $w_j$ is a weight.  It should be clear that due to the properties of the hat functions.  \n",
    "\n",
    "$$\n",
    "    \\tilde{f}(x_i) = \\sum_j^N w_j \\phi_j(x_i) = \\sum_j^N w_j \\delta_{ij} = w_i\n",
    "$$\n",
    "and the weights are just the value of the function at node $j$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "w = f.copy()\n",
    "print(w)\n",
    "phi = hats(x,nodes)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "w_star = w.copy()\n",
    "w_star[3] = 2*w[3]\n",
    "func = w_star.dot(phi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10,5))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.plot(nodes,f,'bo--',label='$f(x)$')\n",
    "axes.plot(x, func,'g-',label='$f^*(x)$')\n",
    "axes.plot(x,phi[2],'r',label='$\\phi_2$')\n",
    "axes.plot(x,phi[3],'c--',label='$\\phi_2$')\n",
    "#axes.plot(x,phi[5],'g--',label='$\\phi_5$')\n",
    "axes.set_xlabel('x')\n",
    "axes.legend()\n",
    "axes.grid()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The element view\n",
    "\n",
    "This may all feel a little familiar and may become clearer if we consider just the element between nodes 2 and 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10,5))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.plot(nodes,f,'bo--',label='$f(x)$')\n",
    "#axes.plot(x, func,'g-',label='$f^*(x)$')\n",
    "axes.plot(x,phi[2],'r',label='$\\phi_2$')\n",
    "axes.plot(x,phi[3],'c--',label='$\\phi_3$')\n",
    "axes.plot(nodes[2]*numpy.ones(2),[0.,f[2]],'g',linewidth=3)\n",
    "axes.plot(nodes[3]*numpy.ones(2),[0.,f[3]],'g',linewidth=3)\n",
    "#axes.plot(x,phi[5],'g--',label='$\\phi_5$')\n",
    "axes.set_xlabel('x')\n",
    "axes.legend()\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "For a single element it should be come clear, that the hat functions are just the lagrange polynomials for linear interpolation within element 2. And we can define local basis functions for each element $N_1$ and $N_2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Interpolation of a function $f(x)$ onto $P_1$\n",
    "\n",
    "Given a set of elements and a local interpolant,  the hat functions form a basis for the space $\\cal V$ of all Piecewise linear functions on the mesh.  If the interpolant is Linear and continuous we call this space $P_1$ and all functions $\\tilde f$ in $P_1$ can be written as\n",
    "\n",
    "$$\n",
    "    \\tilde{f}(x) = \\sum_j w_j\\phi_j(x)\n",
    "$$\n",
    "\n",
    "where $\\phi_j(x)$ is the $j$th basis (hat) function.\n",
    "\n",
    "To *interpolate* a continuous function $f(x)$ onto $P_1$, we simply require that \n",
    "$$\n",
    "    \\tilde{f}(x_i) = f(x_i)\n",
    "$$\n",
    "\n",
    "i.e the interpolant agrees with the function at each node.  It's not hard to see that in this case\n",
    "\n",
    "$$\n",
    "    w_i = f(x_i)\n",
    "$$ \n",
    "and the weights are just the value of the function at the nodes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example\n",
    "\n",
    "let $f(x) =\\sin{\\pi x}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0,L)\n",
    "f = lambda x: numpy.sin(x)\n",
    "w = f(nodes)\n",
    "phi = hats(x,nodes)\n",
    "ftilde = w.dot(phi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10,5))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.plot(x,f(x),'b-',label='$f(x)$')\n",
    "axes.plot(x,ftilde,'r--',label='$\\\\tilde{f}(x)$')\n",
    "axes.plot(nodes,w,'ro')\n",
    "axes.set_xlabel('x')\n",
    "axes.legend()\n",
    "axes.grid()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Galerkin Projection of  a function $f(x)$ onto $P_1$\n",
    "\n",
    "Given a function and its interpolant (or any other function in  $P_1$), we can consider their difference\n",
    "\n",
    "$$\n",
    "    r(x) = f(x) - \\tilde{f}(x) \n",
    "$$\n",
    "\n",
    "which is a function of $x$ we call the residual.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10,5))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.plot(x,f(x),'b-',label='$f(x)$')\n",
    "axes.plot(x,ftilde,'r--',label='$\\\\tilde{f}(x)$')\n",
    "axes.plot(x,f(x)-ftilde,'g:', label='$r(x)$')\n",
    "axes.plot(nodes,w,'ro')\n",
    "axes.set_xlabel('x')\n",
    "axes.legend()\n",
    "axes.grid()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The **projection** of $f$ onto $P_1$ is a different function $f_h$ defined as the function in $P_1$ that minimizes the residual in a least squares sense, i.e. \n",
    "\n",
    "find $f_h\\in P_1$ that minimizes the $L_2$ norm\n",
    "\n",
    "$$\n",
    "    || f(x) - f_h(x) ||^2 = \\int_\\Omega ( f(x) - f_h(x))^2 dx\n",
    "$$ \n",
    "\n",
    "where $\\Omega$ is the domain, in this case $x\\in[0,L]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This problem is reminiscent of the Linear least squares problem for an overdetermined problem \n",
    "\n",
    "$$\n",
    "    A\\mathbf{x} = \\mathbf{b}\n",
    "$$ \n",
    "\n",
    "for $A\\in\\mathbb{R}^{m\\times n}$, $\\mathbf{x}\\in\\mathbb{R}^n$ $\\mathbf{b}\\in\\mathbb{R}^m$ with $m>n$ \n",
    "\n",
    "where we want to find $\\hat{\\mathbf{x}}$ that minimizes  the Euclidean norm of the residual $\\mathbf{r} = A\\mathbf{\\hat{x}} - \\mathbf{b}$, where $A\\mathbf{\\hat{x}}=\\mathbf{p}$ is the *projection* of $\\mathbf{b}$ onto the column space of $A$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We will review that problem first, as if you understand classical linear least-squares it is easy to understand Galerkin projection onto a finite-dimensional function space."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Review:  Orthogonal projection and Least-Squares Problems\n",
    "\n",
    "Consider the overdetermined problem $A\\mathbf{x}=\\mathbf{b}$ where $A\\in\\mathbb{R}^{3\\times2}$ and $\\mathbf{b}\\in\\mathbb{R}^3$ i.e. \n",
    "\n",
    "$$\n",
    "    \\begin{bmatrix} | & | \\\\\n",
    "                    \\mathbf{a}_1 & \\mathbf{a}_2\\\\\n",
    "                     | & | \\\\\n",
    "     \\end{bmatrix} \\begin{bmatrix} x_1 \\\\ x_2\\\\ \\end{bmatrix} \n",
    "         = \\begin{bmatrix} |   \\\\\n",
    "                    \\mathbf{b} \\\\\n",
    "                     |  \\\\ \n",
    "                     \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "and $\\mathbf{a}_1$, $\\mathbf{a}_2$ are linearly independent vectors that span a two-dimensional subspace of $\\mathbb{R}^3$.         "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Geometrically** this problem looks like\n",
    "\n",
    "<img align=center, src=\"./images/least_squares_geometry.jpg\" alt=\"Drawing\" width=600/>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "If $\\mathbf{b}\\notin C(A)$, then there is clearly no solution to $A\\mathbf{x}=\\mathbf{b}$.  However, we can find the point $\\mathbf{p}=A\\hat{\\mathbf{x}}\\in C(A)$ that minimizes the length of the the error $\\mathbf{e}=\\mathbf{b}-\\mathbf{p}$.  While we could resort to calculus to find the values of $\\hat{\\mathbf{x}}$ that minimizes $||\\mathbf{e}||_2$.  It should be clear from the figure that the shortest error (in the $\\ell_2$ norm) is the one that is perpendicular to every vector in $C(A)$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To be orthogonal to every vector in $C(A)$ only requires that the error is orthogonal to each basis vector. i.e. \n",
    "\n",
    "$$\n",
    "     \\mathbf{a}_i^T\\mathbf{e} =0, \\quad  i=1,2\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Which is equivalent to saying\n",
    "\n",
    "$$\n",
    "    A^T\\mathbf{e} = \\begin{bmatrix} - & \\mathbf{a}_1^T & - \\\\ \n",
    "                                    - & \\mathbf{a}_2^T & - \\\\ \\end{bmatrix}\\mathbf{e} =  \\mathbf{0}\n",
    "$$\n",
    "(or that the error must be in the left-Null Space $N(A^T)$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "This implies\n",
    "\\begin{align}\n",
    "    \\mathbf{0} &= A^T\\mathbf{e} \\\\\n",
    "    &= A^T(\\mathbf{b}-\\mathbf{p})\\\\\n",
    "    &= A^T(\\mathbf{b} - A\\hat{\\mathbf{x}})\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or we just need to solve the \"Normal Equations\"  \n",
    "$$ A^T A\\hat{\\mathbf{x}} = A^T\\mathbf{b}$$\n",
    "\n",
    "for the least-squares solution $\\hat{\\mathbf{x}}$  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "and the **projection** $\\mathbf{p}$ of $\\mathbf{b}$ onto $C(A)$ is just\n",
    "\n",
    "$$\n",
    "    \\mathbf{p}= A\\hat{\\mathbf{x}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "While this problem can be solved using calculus,  the better approach is to use linear algebra to find a residual that is orthogonal to $C(A)$, i.e. is orthogonal to all the columns of $A$.  In essence, the smallest residual will be the one that cannot be constructed out of any of the basis vectors of $A$.   (see Lecture Notes `11_LA_QR.ipynb` for more details)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The Galerkin projection of a function onto a finite-dimensional functions space can be found by essentially the same process if we extend the definition of orthogonality and the dot product. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "For vectors in $\\mathbb{R}^n$ the inner product is just the dot product i.e.\n",
    "\n",
    "$$\n",
    "    <\\mathbf{x},\\mathbf{y}> = \\mathbf{x}^T\\mathbf{y}\n",
    "$$  \n",
    "\n",
    "with Euclidean norm\n",
    "\n",
    "$$\n",
    "    ||\\mathbf{x}||^2 = <\\mathbf{x},\\mathbf{x}>\n",
    "$$ \n",
    "\n",
    "Moreover we consider two vectors $\\mathbf{x}$ and $\\mathbf{y}$ to be **orthogonal** when \n",
    "\n",
    "$$<\\mathbf{x},\\mathbf{y}> = 0$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Likewise we can define the $L_2$ inner product for functions as\n",
    "\n",
    "$$\n",
    "    <f, g> = \\int_\\Omega f(x)g(x)dx\n",
    "$$\n",
    "\n",
    "with norm\n",
    "\n",
    "$$\n",
    "    ||f||^2_{L_2} = <f,f>\n",
    "$$\n",
    "\n",
    "and orthogonality of functions whenever $<f,g>=0$\n",
    "\n",
    "Then it becomes straightforward to construct the problem to solve for the Galerkin projection of $f(x)$ onto $P_1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Given a finite dimensional function space $\\cal V$ (e.g. $P_1$) with basis functions \n",
    "\n",
    "$$\n",
    "    \\phi_0(x),\\phi_2(x),\\ldots,\\phi_N(x)\n",
    "$$\n",
    "\n",
    "Then any function in $\\cal V$ can be written as\n",
    "$$\n",
    "    f_h(x) = \\sum_{j=0}^N w_j\\phi_j(x)\n",
    "$$\n",
    "\n",
    "if we also let $f(x)$ be any continuous function (not necessarily in $\\cal V$) then we can define the residual as\n",
    "\n",
    "$$\n",
    "    r(x) = f_h(x) - f(x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "to find the smallest residual, we seek $f_h(x)$ such that $r(x)$ is orthogonal to all the basis functions, i.e.\n",
    "\n",
    "find $\\mathbf{w}$ such that\n",
    "\n",
    "$$\n",
    "    <\\phi_i, r(x)> = 0 \\quad\\forall\\, i=0,\\ldots N\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or\n",
    "\n",
    "$$\n",
    "    \\int_\\Omega \\phi_i\\left[\\sum_j w_j\\phi_j(x)- f(x)\\right]dx = 0, \\quad \\forall\\, i=0,\\ldots N\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or rearranging,   solve for the $w_j$ such that\n",
    "\n",
    "$$\n",
    "    \\int_\\Omega \\phi_i\\sum_j w_j\\phi_j(x)dx = \\int_\\Omega f\\phi_i dx,  \\quad \\forall\\, i=0,\\ldots N\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now the first term of\n",
    "$$\n",
    "    \\int_\\Omega \\phi_i\\sum_j w_j\\phi_j dx = \\int_\\Omega f\\phi_i dx,  \\quad \\forall\\, i=0,\\ldots N\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "can be re-arranged to swap the sum and the integral to get\n",
    "\n",
    "$$\n",
    "    \\sum_j \\left[ \\int_\\Omega \\phi_i\\phi_j dx\\right]w_j = \\int_\\Omega f\\phi_i dx,  \\quad \\forall\\, i=0,\\ldots N\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or \n",
    "$$\n",
    "    \\sum_j <\\phi_i,\\phi_j> w_j = <f,\\phi_i>  \\quad \\forall\\, i=0,\\ldots N\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or letting\n",
    "$$\n",
    "    M_{ij} = <\\phi_i,\\phi_j>,\\quad f_i = <f, \\phi_i>\n",
    "$$\n",
    "\n",
    "the problem reduces to a linear algebra problem $M\\mathbf{w} = \\mathbf{f}$ where $M$ is the **Mass Matrix** and $\\mathbf{f}$ is the **load vector**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Assembling the Mass Matrix and load vector\n",
    "\n",
    "if $\\phi_i$ are the hat functions,  it's not hard to see that the Mass Matrix should be tridiagonal as the only hat functions that should interact around node $i$ should be $\\phi_{i-1}$, $\\phi_i$, $\\phi_{i+1}$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0,L,500)\n",
    "phi = hats(x,nodes)\n",
    "\n",
    "fig = plt.figure(figsize=(10,5))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.plot(x,phi[2],'g--',label='$\\phi_{i-1}$')\n",
    "axes.plot(x,phi[3],'r-',label='$\\phi_i$')\n",
    "axes.plot(x,phi[4],'b--',label='$\\phi_{i+1}$')\n",
    "\n",
    "\n",
    "#axes.plot(x,phi[5],'g--',label='$\\phi_5$')\n",
    "axes.set_xlabel('x')\n",
    "axes.legend()\n",
    "axes.grid()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "i.e.  in the above picture, it should be clear that $<\\phi_{i-1}, \\phi_{i+1} > = 0$.  So that for the $i$th row of $M$,  the only three values should be \n",
    "\n",
    "$$<\\phi_{i-1},\\phi_i>,\\, <\\phi_{i},\\phi_i>,\\, <\\phi_{i},\\phi_{i+1}>$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Assembling the Element tensor\n",
    "\n",
    "It is possible to calculuate the terms of the global mass matrix directly from the hat functions.  However, this can be cumbersome due to the piecewise nature of the hat functions $\\phi_i(x)$.  The standard approach is to loop over elements, and just calculate the contribution from each element.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0,L,500)\n",
    "phi = hats(x,nodes)\n",
    "\n",
    "fig = plt.figure(figsize=(10,5))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.plot(x,phi[2],'g--',label='$\\phi_{2}$')\n",
    "axes.plot(x,phi[3],'r-',label='$\\phi_3$')\n",
    "axes.plot(x,phi[4],'b--',label='$\\phi_4$')\n",
    "for node in nodes:\n",
    "    axes.plot(node*numpy.ones(2),[0.,1],'k:',linewidth=2)\n",
    "for i in range(len(nodes)-1):\n",
    "    axes.text(0.5*(nodes[i+1]+nodes[i])-.06, 0.1, '$e_{}$'.format(i),size=16)\n",
    "    \n",
    "axes.plot(nodes[2]*numpy.ones(2),[0.,1],'m-.',linewidth=3)\n",
    "axes.plot(nodes[3]*numpy.ones(2),[0.,1],'m-.',linewidth=3)\n",
    "\n",
    "\n",
    "#axes.plot(x,phi[5],'g--',label='$\\phi_5$')\n",
    "axes.set_xlabel('x')\n",
    "axes.legend(loc='best')\n",
    "#axes.grid()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "For example.  Element $e_2$,  will contribute to 4 entries in the global Mass matrix.\n",
    "\n",
    "For row $i=2$ (for $\\phi_2$), it will contribute part of $<\\phi_2,\\phi_2>$ and all of $<\\phi_2,\\phi_3>$.\n",
    "\n",
    "For row $i=3$ (for $\\phi_3$), it will contribute all of $<\\phi_2,\\phi_3>$ and part of $<\\phi_3,\\phi_3>$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Assembling the Element tensor\n",
    "\n",
    "For each element we can define local basis Functions $N_1$ and $N_2$ which are just the degree 1 Lagrange polynomials for the linear interpolant over element $i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0,L,500)\n",
    "phi = hats(x,nodes)\n",
    "\n",
    "N1 = lambda t: (1 - t)/2.\n",
    "N2 = lambda t: (1 + t)/2.\n",
    "\n",
    "# affine transformation of element \n",
    "xt = lambda t, coords: coords[0] + (coords[1] - coords[0]) * (t + 1)/2\n",
    "\n",
    "t = numpy.linspace(-1,1)\n",
    "\n",
    "coords = nodes[2:4]\n",
    "\n",
    "fig = plt.figure(figsize=(10,5))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.plot(xt(t,coords),N1(t),'g--',label='$N_1$')\n",
    "axes.plot(xt(t,coords),N2(t),'r-',label='$N_2$')\n",
    "axes.plot(x,numpy.zeros(x.size),'k')\n",
    "for node in nodes:\n",
    "    axes.plot(node*numpy.ones(2),[0.,1],'k:',linewidth=2)\n",
    "for i in range(len(nodes)-1):\n",
    "    axes.text(0.5*(nodes[i+1]+nodes[i])-.06, 0.1, '$e_{}$'.format(i),size=16)\n",
    "    \n",
    "axes.plot(nodes[2]*numpy.ones(2),[0.,1],'m-.',linewidth=3)\n",
    "axes.plot(nodes[3]*numpy.ones(2),[0.,1],'m-.',linewidth=3)\n",
    "\n",
    "\n",
    "#axes.plot(x,phi[5],'g--',label='$\\phi_5$')\n",
    "axes.set_xlabel('x')\n",
    "axes.legend(loc='best')\n",
    "#axes.grid()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Within Element $e_2$ we can construct an *element Mass Matrix* and *element load vector*\n",
    "\n",
    "$$\n",
    "    M^e_2 = \\begin{bmatrix} <N_1,N_1> & <N_1, N_2> \\\\ <N_1, N_2> & <N_2, N_2> \\\\ \\end{bmatrix}\\quad\n",
    "    f^e_2 = \\begin{bmatrix} <N_1, f> \\\\  <N_2, f> \\\\ \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "where each of the definite integrals in the inner product can be evaluated either analytically or usually using quadrature"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Assembling the Global tensors\n",
    "\n",
    "Given each element tensor,  they can then be *added* to the global Mass Matrix and load vector, to assemble the whole system.  In particular, for this problem,  the element mass matrices form overlapping diagonal blocks in $M$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "Me = numpy.ones((2,2))\n",
    "fe = numpy.ones(2)\n",
    "\n",
    "N = 7\n",
    "\n",
    "prop_cycle = plt.rcParams['axes.prop_cycle']\n",
    "colors = prop_cycle.by_key()['color']\n",
    "markers = ['s', 'o', '^', 'P', '*', 'v', 'd']\n",
    "\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(13,6))\n",
    "axes = fig.add_subplot(1,2,1)\n",
    "for i in range(N-1):\n",
    "    M = lil_matrix((N,N))\n",
    "    M[i:i+2,i:i+2] = Me\n",
    "    axes.spy(M,marker=markers[i],color=colors[i])\n",
    "    axes.text(i+.4,i+.6,'$M^e_{}$'.format(i))\n",
    "axes.grid()\n",
    "\n",
    "axes = fig.add_subplot(1,2,2)\n",
    "for i in range(N-1):\n",
    "    f = lil_matrix((N,1))\n",
    "    f[i:i+2] = fe\n",
    "    axes.spy(f,marker=markers[i],color=colors[i])\n",
    "    axes.text(-.1,i+.6,'$f^e_{}$'.format(i))\n",
    "\n",
    "\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Here is some code that assemble both element tensors (using GL quadrature for the integrals in the inner products) and the global tensors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def element_tensors(f, coords):\n",
    "    \"\"\" returns element mass matrix and load vector for P1 elements using 2-point Gauss-Legendre quadrature rules\n",
    "    \n",
    "    parameters\n",
    "    ----------\n",
    "    f: calleable\n",
    "        RHS function\n",
    "    coords: ndarray\n",
    "        x coordinates of element\n",
    "    \"\"\"\n",
    "    # element basis functions (lagrange polynomials on [-1,1])\n",
    "    N1 = lambda t: (1 - t)/2.\n",
    "    N2 = lambda t: (1 + t)/2.\n",
    "    \n",
    "    # 2-pt GL points and weights\n",
    "    tQ = numpy.array([-1., 1.])/numpy.sqrt(3)\n",
    "    wQ = numpy.ones(2)\n",
    "    \n",
    "    # element width\n",
    "    h = coords[1] - coords[0]\n",
    "    \n",
    "    # affine transformation of element \n",
    "    x = lambda t: coords[0] + h * (t + 1)/2\n",
    "    \n",
    "    # element mass matrix\n",
    "    Me = numpy.empty((2,2))\n",
    "    Me[0, 0] = wQ.dot(N1(tQ)*N1(tQ))\n",
    "    Me[0, 1] = wQ.dot(N1(tQ)*N2(tQ))\n",
    "    Me[1, 0] = Me[0, 1]\n",
    "    Me[1, 1] = wQ.dot(N2(tQ)*N2(tQ))\n",
    "    Me *= h/2.\n",
    "    \n",
    "    # load vector\n",
    "    \n",
    "    fe = numpy.empty(2)\n",
    "    fe[0] = wQ.dot(f(x(tQ))*N1(tQ))\n",
    "    fe[1] = wQ.dot(f(x(tQ))*N2(tQ))\n",
    "    fe *= h/2.\n",
    "    \n",
    "    return Me, fe\n",
    "\n",
    "\n",
    "def assemble_global(func,nodes):\n",
    "    \"\"\" Assemble the global Mass Matrix and load vector given a P1 mesh with coordinates at nodes\n",
    "    \n",
    "    parameters\n",
    "    -----------\n",
    "    func: calleable\n",
    "        RHS function\n",
    "    nodes: ndarray\n",
    "        x coordinates of mesh\n",
    "    \"\"\"\n",
    "    # number of nodes and elements\n",
    "    N = len(nodes)\n",
    "    Nel = N - 1\n",
    "    \n",
    "    # allocate sparse matrix and vector \n",
    "    M = lil_matrix((N,N))\n",
    "    f = numpy.zeros(N)\n",
    "    \n",
    "    # loop over elements and assemble global matrix\n",
    "    for k in range(Nel):\n",
    "        Me, fe = element_tensors(func,nodes[k:k+2])\n",
    "        M[k:k+2,k:k+2] += Me\n",
    "        f[k:k+2] += fe\n",
    "        \n",
    "    return M.tocsr(), f\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Solve the projection problem\n",
    "\n",
    "again on the same mesh use $f = \\sin(x)$ but find the projection $f_h$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "N=11\n",
    "x = numpy.linspace(0., L, 100)\n",
    "nodes = numpy.linspace(0.,L, N)\n",
    "phi = hats(x, nodes)\n",
    "\n",
    "func = lambda x: numpy.sin(x)\n",
    "M, f = assemble_global(func,nodes)\n",
    "w = spsolve(M,f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10,5))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.plot(x,func(x),'b-',label='$f(x)$')\n",
    "axes.plot(x,w.dot(phi),'r--',label='$f_h(x)$')\n",
    "axes.plot(x,func(nodes).dot(phi),'g:',label='$\\\\tilde{f}(x)$')\n",
    "axes.plot(nodes,w,'ro')\n",
    "axes.set_xlabel('x')\n",
    "axes.legend()\n",
    "axes.grid()\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Convergence Analysis\n",
    "\n",
    "Beyond just finding the least squares solution,  there are some important convergence results that can also be shown. (e.g. see Larson and Bengzon, \"The Finite Element Method\" for more details)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### Theorem 1: Best approximation\n",
    "\n",
    "The $L^2$-projection $f_h$, satisfies the best approximation result\n",
    "\n",
    "$$\n",
    "    || f - f_h ||_{L^2} \\leq || f - v ||_{L^2},\\ \\ \\forall v \\in V_h\n",
    "$$\n",
    "\n",
    "where $V_h$ is the finite element function space with basis $\\phi_1, \\phi_2,\\ldots,\\phi_N$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### Theorem 2:  Convergence\n",
    "\n",
    "for $V_h=P_1$, the space of all piecewise linear polynomial spaces, then an *a priori* error estimate can be found such that\n",
    "\n",
    "$$\n",
    "    || f - f_h ||_{L^2} \\leq C h^2|| f''|_{L^2}\n",
    "$$\n",
    "\n",
    "where $h$ is the width of the largest element.  Therefore, as $h\\rightarrow 0$, the projection converges to the true function in the $L^2$ norm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Convergence Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "Ns = [ (1+2**n) for n in range(3,9) ]\n",
    "rel_err = numpy.zeros(len(Ns))\n",
    "delta_x = numpy.zeros(len(Ns))\n",
    "\n",
    "for i,N in enumerate(Ns):\n",
    "    mesh = numpy.linspace(0, L, N)\n",
    "    phi = hats(x,mesh)\n",
    "    M, f = assemble_global(func,mesh)\n",
    "    u = spsolve(M,f)\n",
    "    rel_err[i] = numpy.linalg.norm(u.dot(phi) - func(x))/numpy.linalg.norm(u_true(x))\n",
    "    delta_x[i] = mesh[1] - mesh[0]\n",
    "    print('n: {}, rel_err = {}'.format(N,rel_err[i]))\n",
    "\n",
    "# calculate best fit slope and Plot result\n",
    "p = numpy.polyfit(numpy.log(delta_x),numpy.log(rel_err),1)\n",
    "dx = numpy.logspace(numpy.log10(delta_x[0]),numpy.log10(delta_x[-1]),100)\n",
    "err = numpy.exp(p[1])*dx**p[0]\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.loglog(delta_x, rel_err, 'o',label='error', markersize=10 )\n",
    "axes.loglog(dx,err,'r--',label=\"p={}\".format(p[0]))\n",
    "axes.set_title(\"Convergence, FEM method\")\n",
    "axes.set_xlabel(\"$\\Delta x$\")\n",
    "axes.set_ylabel(\"$||u - u_{true}||/||u_{true}||$\")\n",
    "axes.legend(loc='best')\n",
    "axes.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6,6))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.spy(M)\n",
    "axes.grid()\n",
    "axes.set_title('Visualization of Mass Matrix')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Solving 2-pt boundary value problems using Galerkin Finite Elements\n",
    "\n",
    "If you understand the Galerkin projection problem, extensions to solving two-point boundary value problems such as \n",
    "\n",
    "$$\n",
    "    u_{xx} - u + f(x) = 0,\\quad x\\in[0,L] \n",
    "$$\n",
    "\n",
    "with appropriate boundary conditions are relatively straightforward and follow much of the same workflow. \n",
    "\n",
    "Here we will work out a few of the details to solve this exact problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Given the 2-pt Boundary value problem\n",
    "\n",
    "$$\n",
    "    u_{xx} - u + f(x) = 0\n",
    "$$\n",
    "\n",
    "on the domain $x\\in[0,L]$ with dirichlet boundary conditions\n",
    "\n",
    "$$\n",
    "    u(0)=u_0,\\quad u(L)=u_L\n",
    "$$\n",
    "\n",
    "we want to find the best Piecewise Linear function in $P_1$ that satisfies the ODE in a \"weak sense\".  I.e. it should be pretty clear that there is no function in $P_1$ that actually satisfies the ODE (because $u_{xx}$ is not well defined in $P_1$).  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "However,  we can still find a function $u_h\\in P_1$ that minimizes the residual\n",
    "\n",
    "$$\n",
    "    r(x) = (u_h)_{xx} - u_h + f(x)\n",
    "$$\n",
    "\n",
    "By requiring that the residual is orthogonal to all the basis functions of $P_1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or\n",
    "$$\n",
    "    \\int_0^L \\phi_i(x) r(x) dx = 0,\\quad \\forall\\, i=0,\\ldots,N\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "expanding the definition of the residual gives the \"weak form\" of the equations\n",
    "\n",
    "$$\n",
    "    \\int_0^L \\phi_i \\frac{d^2 u_h}{dx^2} dx -  \\int_0^L \\phi_i u_h dx  +  \\int_0^L \\phi_i f dx = 0,\\quad \\forall\\, i=0,\\ldots,N\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Note:**  The last two terms are basically the same as in the projection problem.  The only new addition is the first term"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Integration by parts\n",
    "\n",
    "A standard trick in Finite elements is to take terms involving higher order derivatives (in this case $u_{xx}$) and integrating by parts to move one derivative off the trial function and put it on the basis functions (aka test functions).  i.e with no loss of generality\n",
    "\n",
    "$$\n",
    "\\int_0^L \\phi_i \\frac{d^2 u_h}{dx^2} dx  = - \\int_0^L \\frac{d \\phi_i}{dx}\\frac{du_h}{dx}dx + \\left.\\phi_i\\frac{du_h}{dx}\\right|_0^L\n",
    "$$\n",
    "\n",
    "where the last term accomodates Neumann Boundary conditions on the edge of the domain (if they exist) if $\\frac{d u}{dx}$ is known"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Weak Form of the ODE\n",
    "\n",
    "The full weak form of the ODE  can be written,\n",
    "\n",
    "find $u_h\\in P1$ such that \n",
    "\n",
    "$$\n",
    "    \\int_0^L \\frac{d \\phi_i}{dx}\\frac{du_h}{dx}dx +  \\int_0^L \\phi_i u_h dx  =  \\int_0^L \\phi_i f dx + \\left.\\phi_i\\frac{du_h}{dx}\\right|_0^L,\\quad \\forall\\, i=0,\\ldots,N\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Following the same procedure as in the projection problem and substituting in \n",
    "\n",
    "$$\n",
    "    u_h = \\sum_{j=0}^N w_j\\phi_j(x)\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "transforms the weak form into the Linear Algebraic problem\n",
    "\n",
    "$$\n",
    "    (A + M)\\mathbf{w} = \\mathbf{f} \n",
    "$$\n",
    "\n",
    "where $M$ and $\\mathbf{f}$ are the mass matrix and load vector (possibly modified for boundary Conditions terms) as before "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "and $A$ is the \"stiffness\" matrix with entries\n",
    "\n",
    "$$\n",
    "    A_{ij} = <\\frac{d \\phi_i}{dx},\\frac{d \\phi_j}{dx}>\n",
    "$$\n",
    "\n",
    "With possible modifications for Dirichelt BC's"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### More code\n",
    "\n",
    "And here is a bit of code for assembling all of the appropriate tensors again using GL quadrature and assembling over elements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def element_tensors(f, coords):\n",
    "    \"\"\" returns element stiffness and mass matrix and load vector for P1 elements using 2-point Gauss-Legendre quadrature rules\n",
    "    \n",
    "    parameters\n",
    "    ----------\n",
    "    f: calleable\n",
    "        RHS function\n",
    "    coords: ndarray\n",
    "        x coordinates of element\n",
    "    \"\"\"\n",
    "    # element basis functions (lagrange polynomials on [-1,1])\n",
    "    N1 = lambda t: (1 - t)/2.\n",
    "    N2 = lambda t: (1 + t)/2.\n",
    "    \n",
    "    # Gradient of basis functions \n",
    "    dN1dt = lambda t: -1./2.*numpy.ones(t.shape)\n",
    "    dN2dt = lambda t: 1./2.*numpy.ones(t.shape)\n",
    "    \n",
    "    # 2-pt GL points and weights\n",
    "    tQ = numpy.array([-1., 1.])/numpy.sqrt(3)\n",
    "    wQ = numpy.ones(2)\n",
    "    \n",
    "    # element width\n",
    "    h = coords[1] - coords[0]\n",
    "    \n",
    "    # affine transformation of element \n",
    "    x = lambda t: coords[0] + h * (t + 1)/2\n",
    "    \n",
    "    # element mass matrix\n",
    "    Me = numpy.empty((2,2))\n",
    "    Me[0, 0] = wQ.dot(N1(tQ)*N1(tQ))\n",
    "    Me[0, 1] = wQ.dot(N1(tQ)*N2(tQ))\n",
    "    Me[1, 0] = Me[0, 1]\n",
    "    Me[1, 1] = wQ.dot(N2(tQ)*N2(tQ))\n",
    "    Me *= h/2.\n",
    "    \n",
    "    # element Stiffness matrix\n",
    "    Ae = numpy.empty((2,2))\n",
    "    Ae[0, 0] = wQ.dot(dN1dt(tQ)*dN1dt(tQ))\n",
    "    Ae[0, 1] = wQ.dot(dN1dt(tQ)*dN2dt(tQ))\n",
    "    Ae[1, 0] = Ae[0, 1]\n",
    "    Ae[1, 1] = wQ.dot(dN2dt(tQ)*dN2dt(tQ))\n",
    "    Ae *= 2./h\n",
    "    \n",
    "    # load vector\n",
    "    \n",
    "    fe = numpy.empty(2)\n",
    "    fe[0] = wQ.dot(f(x(tQ))*N1(tQ))\n",
    "    fe[1] = wQ.dot(f(x(tQ))*N2(tQ))\n",
    "    fe *= h/2.\n",
    "    \n",
    "    return Me, Ae,  fe\n",
    "\n",
    "\n",
    "def assemble_global(func,nodes):\n",
    "    \"\"\" Assemble the global  Matrix (A+M) and load vector given a P1 mesh with coordinates at nodes\n",
    "    This version assumes Homogeneous dirichlet BC's\n",
    "    \n",
    "    parameters\n",
    "    -----------\n",
    "    func: calleable\n",
    "        RHS function\n",
    "    nodes: ndarray\n",
    "        x coordinates of mesh\n",
    "    \"\"\"\n",
    "    # number of nodes and elements\n",
    "    N = len(nodes)\n",
    "    Nel = N - 1\n",
    "    \n",
    "    # allocate sparse matrix and vector \n",
    "    M = lil_matrix((N,N))\n",
    "    f = numpy.zeros(N)\n",
    "    \n",
    "    # loop over elements and assemble global matrix A+M\n",
    "    for k in range(Nel):\n",
    "        Me, Ae, fe = element_tensors(func,nodes[k:k+2])\n",
    "        M[k:k+2,k:k+2] += Me + Ae\n",
    "        f[k:k+2] += fe\n",
    "        \n",
    "    # fix Boundary conditions\n",
    "    M[0,0:3] = [ 1.,0.,0.]\n",
    "    M[-1,-3:] = [0., 0., 1.]\n",
    "    f[[0,-1]] = 0.\n",
    "        \n",
    "    return M.tocsr(), f\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Test with a \"manufactured solution\"\n",
    "\n",
    "While it is not difficult to solve this problem analytically for various simple source terms, another approach that is very powerful for checking numerical methods is to create a \"manufactured solution\",  i.e. choose a known functional form for $u$ that satisfies the boundary conditions and then calculate the appropriate forcing function $f(x)$ with that should produce you chosen solution.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example\n",
    "\n",
    "consider the problem\n",
    "\n",
    "$$\n",
    "    u_{xx} - u + f(x) = 0\n",
    "$$\n",
    "\n",
    "with homogeneous dirichlet BC's on the interval $x\\in[0,L]$ (i.e. $u(0) = u(L) = 0$).  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "If we assume the true solution is \n",
    "\n",
    "$$u(x) = \\sin\\left(\\frac{n\\pi x}{L}\\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "then this solution would be satisfied if \n",
    "\n",
    "$$ f(x) = u - u_{xx} = \\left[1 + \\left(\\frac{n\\pi}{L}\\right)^2\\right]\\sin\\left(\\frac{n\\pi x}{L}\\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# set the manufactured solution\n",
    "\n",
    "L=1.\n",
    "n=3\n",
    "u_true = lambda x: numpy.sin(n*numpy.pi*x/L)\n",
    "func = lambda x : (1+ (n*numpy.pi/L)**2)*u_true(x)\n",
    "\n",
    "#Assemble matrices and RHS vectors for a uniform mesh with N points\n",
    "N = 11\n",
    "mesh = numpy.linspace(0,L,N)\n",
    "\n",
    "\n",
    "A, f = assemble_global(func,mesh)\n",
    "\n",
    "# and solve\n",
    "u = spsolve(A,f)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0,L,513)\n",
    "phi = hats(x,mesh)\n",
    "\n",
    "err = u_true(x) - u.dot(phi)\n",
    "rel_err=numpy.linalg.norm(err)/numpy.linalg.norm(u_true(x))\n",
    "fig = plt.figure(figsize=(10,5))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.plot(x,u_true(x),'b-',label='$u_{true}(x)$')\n",
    "axes.plot(mesh,u,'r-o',label='$u_h(x)$')\n",
    "axes.plot(x,err,'g:',label='error')\n",
    "axes.set_xlabel('x')\n",
    "axes.legend()\n",
    "axes.grid()\n",
    "axes.set_title('$||e|| = {}$'.format(rel_err))\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### And the always necessary convergence plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "Ns = [ (1+2**n) for n in range(3,9) ]\n",
    "rel_err = numpy.zeros(len(Ns))\n",
    "delta_x = numpy.zeros(len(Ns))\n",
    "\n",
    "for i,N in enumerate(Ns):\n",
    "    mesh = numpy.linspace(0, L, N)\n",
    "    phi = hats(x,mesh)\n",
    "    A, f = assemble_global(func,mesh)\n",
    "    u = spsolve(A,f)\n",
    "    rel_err[i] = numpy.linalg.norm(u.dot(phi) - u_true(x))/numpy.linalg.norm(u_true(x))\n",
    "    delta_x[i] = mesh[1] - mesh[0]\n",
    "    #print('n: {}, rel_err = {}'.format(N,rel_err[i]))\n",
    "\n",
    "# calculate best fit slope and Plot result\n",
    "p = numpy.polyfit(numpy.log(delta_x),numpy.log(rel_err),1)\n",
    "dx = numpy.logspace(numpy.log10(delta_x[0]),numpy.log10(delta_x[-1]),100)\n",
    "err = numpy.exp(p[1])*dx**p[0]\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.loglog(delta_x, rel_err, 'o',label='error', markersize=10 )\n",
    "axes.loglog(dx,err,'r--',label=\"p={}\".format(p[0]))\n",
    "axes.set_title(\"Convergence, FEM method\")\n",
    "axes.set_xlabel(\"$\\Delta x$\")\n",
    "axes.set_ylabel(\"$||u - u_{true}||/||u_{true}||$\")\n",
    "axes.legend(loc='best')\n",
    "axes.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Themes and variations\n",
    "\n",
    "These notes provide only the briefest of introductions to Finite Element Methods which are a vast and mathematically rich field of numerical methods and analysis.  However, they demonstrate how the ideas of interpolation, quadrature, function approximation and linear algebra can be combined to solve ODE's and PDE's.\n",
    "\n",
    "Some critical issues that still need to be explored include\n",
    "\n",
    "* Error and Stability Analysis:  this really requires a solid understanding of functional analysis\n",
    "* Higher Dimensional FEM methods\n",
    "* Higher Order FEM methods\n",
    "* More exotic function spaces than lagrange elements (Discontinous Galerkin, H-div elements)\n",
    "* Better software packages for solving FEM problems\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.7"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}