{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A primer on the Finite Element Method\n",
    "\n",
    "This notebook provides a demonstration of the Finite Element Method, a technique to obtain a numerical solution of a partial differential equation. This is by no means a complete introduction, but serves more as a numerical example of how the technique works in practice. For more detail on the mathematics behind it, there are numerous books and resources online; for open-source implementations, you can check out [FEniCS](https://fenicsproject.org/), [DUNE](https://www.dune-project.org/modules/dune-fem/), [deal.II](www.dealii.org) and many others.\n",
    "\n",
    "## Description\n",
    "We are looking for the solution to the partial differential equation \n",
    "\n",
    "$$-\\Delta u=f$$ \n",
    "\n",
    "on a domain $\\Omega$. Here, $\\Delta$ is the Laplace operator and $f$ is a known function. We assume that on the boundary of our domain $\\partial \\Omega$, $u$ takes value $0$.\n",
    "\n",
    "To find a numerical solution, we apply the finite element method (FEM). The starting point of this method is that instead of solving the equation in the form above, we solve a _weak form_, obtained by multiplying the equation with a _test function_ $\\phi$ and integrate over the domain $\\Omega$. After integration by parts, this results in the following problem:\n",
    "Find $u$ such that for all test functions $\\phi$ we have\n",
    "\n",
    "$$\\int_\\Omega\\nabla u\\cdot \\nabla\\phi = \\int_\\Omega f \\phi.$$\n",
    "\n",
    "The mentioned test functions belong to a set of functions with a specific property (in this case, a Sobolev space). If we find a function $u$ such that our equation is satisfied for all those test functions, we have found a _weak solution_.\n",
    "\n",
    "One of the differences between a weak and a regular (strong) solution, is that the weak solution does not require the same amount of 'smoothness'. In our case, the weak solution does not need to have a second derivative (while the strong one does). However, it's good to know a weak solution is a generalization of a strong solution, meaning that if a strong solution exists, it will be equal to the weak one.\n",
    "\n",
    "In this notebook, we apply FEM to the one-dimensional example (where $\\Delta u$ reduces to $u''$), and investigate how the grid influences the exactness of the solution.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We solve the Laplace equation on the domain $\\Omega=[0,1]$. We leave the right hand side open to experiment with different functions. As boundary conditions, we prescribe $u(0) = u(1) = 0$.\n",
    "To obtain a FEM approximation, we require the following steps.\n",
    "1. Discretize the domain\n",
    "2. Choose basis functions $\\phi$\n",
    "3. Compute stiffness matrix and mass vector\n",
    "4. Solve the linear system\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Discretize the domain\n",
    "We discretize the domain by partitioning it into a grid with grid points $x_0,...,x_{n}$.\n",
    "The patch between subsequent grid points is commonly referred to as an _element_. This puts the 'Element' in Finite Element Method. (The 'Finite' stems from the fact that we cannot choose infinitely small elements; they need to have a _finite_ size). For now, we keep the size of an element constant ($h$). In two dimensions an element commonly takes a triangular (or quadrilateral) shape, in three dimensions a tetrahedon. We are going to compute an approximation of the solution for each element."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Number of elements\n",
    "n = 10\n",
    "# Element size\n",
    "h = 1./n\n",
    "# Nodes\n",
    "xs = np.linspace(0,1,n+1)\n",
    "# Right hand side function (expressed as an anonymous function)\n",
    "f = lambda x: 6*x - 2\n",
    "# Using this function, we know the solution must be u=x**2 - x**3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Choose basis functions\n",
    "Now that the domain $[0,1]$ is partioned in $n$ elements, we can fix our test functions.\n",
    "Instead of looking at the entire set of test functions, we look at combinations of _basis functions_. Much like basis vectors in linear algebra, other test functions can be expressed as linear combinations of these basis functions.\n",
    "\n",
    "By fixing a finite sized set of basis functions $\\phi_i$, we can expand the weak form to a system of equations: for all $i=1,...,n-1$, find $u$ such that\n",
    "\n",
    "$$\\int_\\Omega\\nabla u\\cdot \\nabla\\phi_i = \\int f \\phi_i$$\n",
    "\n",
    "is satisfied.\n",
    "\n",
    "In this example, we choose piecewise linear basis functions $\\phi_i$, $1\\leq i \\leq n-1$, defined as\n",
    "\n",
    "$$ \\phi_i = \\begin{cases}0 & \\mbox{ if }x < x_{i-1}\\\\ \\frac{x-x_{i-1}}{x_{i}-x_{i-1}}& \\mbox{ if }x_{i-1} < x < x_{i}\\\\ \\frac{x_{i+1}-x}{x_{i+1}-x_{i}}& \\mbox{ if }x_{i} < x < x_{i+1}\\\\ 0 & \\mbox{ if }x > x_{i+1}\\end{cases}$$\n",
    "\n",
    "In other words: within each element, the functions are linear. On their corresponding grid node they have value 1, on every other node they have value 0. For higher precision, people often choose higher-degree polynomials, like quadratic or qubic functions.\n",
    "\n",
    "For $n=10$, they look like this ($i=4$ is highlighted). Note that we won't actually use the `phi`-variable in the FEM formulation, it just serves as an illustration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEWCAYAAABBvWFzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztvXmYFfWV//8+QLessi8KsqgooAhCyybIKjZbV/KYReNo\nxkSN+Zp5MjOZSeI3k5htJpNkEjNGHWMSf06+STQ+iVpXdmRHZEcBFRUBFVxYRDbZOb8/zv3Qty+3\n7q2691Ofut2c1/P009213Hv6dN069XnX57yLmBmKoiiKUiqNkg5AURRFaRhoQVEURVGsoAVFURRF\nsYIWFEVRFMUKWlAURVEUK2hBURRFUaygBUVRsiCiWUT0xSL2+zQRvUtEh4jo6jhiC3jfW4horqv3\nU5QgSPtQlPoKEW0H0BnAKQAnACwHcDczv5tQPG8B+Gdm9mN8j54AtgGoYOaTcb2PohSDjlCU+s40\nZm4J4AIAHwL4dYKx9ADwSoLvryiJogVFaRAw81EAfwXQzywjoilEtJ6IDqSlqO9nrGtKRH8kor1E\n9DERrSaizul1i4jojvTPlxLRYiLaT0R7iOgv2e9NROcR0SEAjQG8nB6pgIiYiC7N2O5xIvpx+ucx\nRLSDiL5BRLuI6H0iuj1j22ZE9Asiejv93suIqBmAJelNPk5La8OJ6O+JaFnGviPSf8/+9PcRGesW\nEdGPiOgFIjpIRHOJqEOhnChKGLSgKA0CImoO4PMAVmQsPgzgNgBtAEwB8FUi+lR63RcBtAZwEYD2\nAO4GcCTHS/8IwFwAbQF0Q44REDMfS4+SAGAAM18SMuwu6Ri6AvgygIeIqG163X8BGAxgBIB2AL4J\n4DSA69Lr2zBzS2Z+MfMFiagdgBkAHkj/Xb8EMIOI2mds9gUAtwPoBKASwL+kl4fNiaLkRAuKUt95\nlog+BrAfwPUAfm5WMPMiZt7IzKeZeQOAJwCMTq8+ATlpXsrMp5h5LTMfyPH6JyBS1oXMfJSZl+XY\nplhOAPghM59g5pkADgG4nIgaAfgSgK8z8850fMuZ+ViI15wC4E1m/n/MfJKZnwCwGcC0jG3+P2Z+\ng5mPAHgKwMCMeMLkRFFyogVFqe98ipnbAGgK4GsAFhNRFwAgoqFEtJCIdhPRfsgVd4f0fv8PwBwA\nTxLRe0T0MyKqyPH63wRAAFYR0StE9CWLse/NurH+CYCW6RibAniriNe8EMDbWcvehoyCDB/keE8g\nfE4UJSdaUJQGQfqK+mnIjK+R6cV/BpACcBEztwbwCKQ4ID0q+AEz94PISlMh8lj2637AzHcy84UA\nvgLg4cz7IgX4BEDzjN+7hNxvD4CjAHJJZ4WmZb4HGVFl0h3AzkJvGjYnihKEFhSlQUCCB7nX8Vp6\ncSsAHzHzUSIaArl3YLYfS0T9iagxgAMQued0jtf9LBF1S/+6D3JCP2u7AF4C8AUiakxE1aiV2/LC\nzKcBPAbgl0R0YXr/4UR0HoDd6fe/OGD3mQAuI6IvEFETIvo8ZKLC9ELvGzYnihKEFhSlvvNceobV\nAQD/DuCLzGym7v4fAD8kooMAvge5X2DoApkVdgBSgBZDJJ9srgGwMv0eKch9ja0hY/s65N7FxwBu\nAfBshL/rXwBsBLAawEcAfgqgETN/Avk7X0jPxBqWuRMz74WMLL4BYC9EspvKzHtCvGfYnChKTrSx\nUVEURbGCjlAURVEUKzgrKER0UXrGzavp2TJfz7ENEdEDRLSFiDYQ0aCMddVE9Hp63bddxa0oiqKE\nw+UI5SSAb6RnkAwDcA8R9cvaZhKA3umvuwD8DwCkbxI+lF7fD8DNOfZVFEVREsRZQWHm95l5Xfrn\ng5Cbfl2zNvMA/IGFFQDaENEFAIYA2MLMW5n5OIAn09sqiqIoZUKTJN6UxDH1agArs1Z1BZDpFLsj\nvSzX8qEBr30XZHSDFi1aDO7Tp4+VmBVFUc4F1q5du4eZOxazr/OCQkQtAfwNwD/GYevAzI8CeBQA\nqqqqeM2aNbbfQlEUpcFCRNlOC6FxWlDSNg5/A/CndFdzNjshxnSGbullFQHLFUVRlDLB5SwvAvB7\nAK8x8y8DNksBuC0922sYgP3M/D6kuas3EfUiokoAN6W3VRRFUcoElyOUawHcCmAjEb2UXvZ/IT5D\nYOZHILYRkwFsgfgg3Z5ed5KIvgYxrmsM4LGMbmhFURSlDHBWUNK231RgGwZwT8C6mZCCoyiKopQh\n2imvKIqiWEELiqIoimIFLSiKoiiKFbSgKIqiKFbQgqIoiqJYQQuKoiiKYgUtKIqiKIoVtKAoiqIo\nVtCCoiiKolhBC4qiKIpiBS0oiqIoihW0oCiKoihW0IKiKIqiWEELiqIoimIFLSiKoiiKFbSgKIqi\nKFbQgqIoiqJYQQuKoiiKYgVnjwAmoscATAWwi5mvzLH+XwHckhFXXwAdmfkjItoO4CCAUwBOMnOV\nm6gVRVGUsLgcoTwOoDpoJTP/nJkHMvNAAPcCWMzMH2VsMja9XouJoihKGeKsoDDzEgAfFdxQuBnA\nEzGGoyiKolim7O6hEFFzyEjmbxmLGcDzRLSWiO5KJjJFURQlH87uoURgGoAXsuSukcy8k4g6AZhH\nRJvTI56zSBecuwCge/fu8UerKIqiACjDEQqAm5AldzHzzvT3XQCeATAkaGdmfpSZq5i5qmPHjrEG\nqiiKotRSVgWFiFoDGA3Az1jWgohamZ8BTASwKZkIFUVRlCBcTht+AsAYAB2IaAeA+wBUAAAzP5Le\n7NMA5jLz4YxdOwN4hohMvH9m5tmu4lYURVHC4aygMPPNIbZ5HDK9OHPZVgAD4olKURRFsUVZSV6K\noihK/UULiqIoimIFLSiKoiiKFbSgKIqiKFbQgqIoiqJYQQuKoiiKYgUtKIqiKIoVtKAoiqIoVtCC\noiiKolhBC4qiKIpiBS0oiqIoihW0oCiKoihW0IKiKIqiWEELiqIoimIFLSiKks3Ro0BFBUAEvP12\n0tEoSr1BC4qiZLNxI3DyJNC2LdCjR9LRKEq9QQuKomSzfr18v/rqZONQlHqGFhRFyUYLiqIUhRYU\nRcnGFJSBA5ONQ1HqGc4KChE9RkS7iGhTwPoxRLSfiF5Kf30vY101Eb1ORFuI6NuuYlbOIZ5+Grj+\nerlvsnKlLLvjDmDAAODRR5ONTVHqCU0cvtfjAB4E8Ic82yxl5qmZC4ioMYCHAFwPYAeA1USUYuZX\n4wpUOYf45BPgi18E/vpX+b1pU/neqBHADGzYAHzlK8C2bcBPfpJcnIpSD3A2QmHmJQA+KmLXIQC2\nMPNWZj4O4EkAntXglHOT06eBm26SYnLxxcDixbWjkSFDgI8+Am69VX7/2c+AN95ILlZFqQeU2z2U\nEUS0gYhmEdEV6WVdAbybsc2O9LKcENFdRLSGiNbs3r07zliV+s5vfgM89xzQogUwfz5w3XXASy/J\nuoEDZfkjjwCtW0vxeeaZZONVlDKnnArKOgDdmfkqAL8G8GwxL8LMjzJzFTNXdezY0WqASgPixAng\nhz+Un3/8Y6BnT/k5s6AAQPPmwLBh8vNbbzkNUVHqG2VTUJj5ADMfSv88E0AFEXUAsBPARRmbdksv\nU5TimTkT+OADoLJSbr4bXn5ZvmfO8GrZUr6fPp37tf70J+mqN/dfFOUcpWwKChF1ISJK/zwEEtte\nAKsB9CaiXkRUCeAmAKnkIlUaBPPmyffRo2sLxo4dwN69ckO+f//abT/8UL5363b26xw6BHzzmyKP\nKco5jrNZXkT0BIAxADoQ0Q4A9wGoAABmfgTAZwB8lYhOAjgC4CZmZgAniehrAOYAaAzgMWZ+xVXc\nSgPFSFtDh569rHdvkboA4NQpmekF5G50/PGPgVatgLFja2eKKco5irOCwsw3F1j/IGRaca51MwHM\njCMu5Rzl3fQ8j86da5dl3z8BgKVLgQMHRBobPbrua7z5JnD//XKz/qmn4o1XUeoBZSN5KYpTRF2V\nYmHIVVAeeEC+33gj0KZN3df4x3+UkcnkyfHFqSj1CC0oyrnJJZfI9wULapdlF5Snn5bRR2UlcN99\ndfefMQOYOxf45S/jj1VR6glaUJRzkxtvlO/z5wMPPwzs3w9s3SrLrrgC+N3vgFtukd9/+Uvg8str\n9z1+HPinfwLuvhvo189t3IpSxmhBUc5N7rgDGD5cfr7nHuCCC8RqBQAuuwy4807pVbn/flmfyf33\nA3v2AD/4gduYFaXM0YKinJtUVopk9S//Alx0EXDsmCwnks74L3wBWLdO7pNksn+/zOy68065/7J9\nu3wdOiQFaft2YNcux3+MopQHxOaqrAFSVVXFa9asSToMpT7w938P/O//At/5jhSMILZvB3r1yv9a\nU6YA06fbjE5RnEFEa5m5qph9XboNK0r5kqtDPhedOuX29HrgAWDZMpk+3KWL/fgUpR6gBUVRTpwA\nXk0/DWHAgPzbNm8OfOpTZy9/9lnpsM+1TlHOEfQeiqJs3iwzt1q0qJ1OrChKZLSgKIqRu/r3l1FG\nMTz+OHD0qLWQFKU+ogVFUf7u72SG1osvJh2JotRrtKAoiqIoVtCCoiiKolhBC4qiKIpiBS0oiqIo\nihW0oCiKoihW0IKiKIqiWEELiqIoimIFZwWFiB4jol1EtClg/S1EtIGINhLRciIakLFue3r5S0Sk\nbo+KoihliMsRyuMAqvOs3wZgNDP3B/AjAI9mrR/LzAOLdcFUFEVR4sWZOSQzLyGinnnWL8/4dQWA\nbnHHpCiKotijXO+hfBnArIzfGcDzRLSWiO7KtyMR3UVEa4hoze7du2MNUlEURaml7OzriWgspKCM\nzFg8kpl3ElEnAPOIaDMzL8m1PzM/irRcVlVV1XCfHqYoilJmlNUIhYiuAvA7AB4z7zXLmXln+vsu\nAM8AGJJMhIqiKEoQZVNQiKg7gKcB3MrMb2Qsb0FErczPACYCyDlTTFEURUkOZ5IXET0BYAyADkS0\nA8B9ACoAgJkfAfA9AO0BPExEAHAyPaOrM4Bn0suaAPgzM892FbeiKIoSDpezvG4usP4OAHfkWL4V\nQIHnsiqKoihJUzaSl6IoilK/0YKiKIqiWEELiqIoimIFLSiKoiiKFbSgKIqiKFbQgqIoiqJYQQuK\noiiKYgUtKIqiKIoVtKAoiqIoVtCCoiiKolhBC4qiKIpiBS0oiqIoihUiF5S0nXzjOIJRFEVR6i8F\nCwoRNSKiLxDRDCLaBWAzgPeJ6FUi+jkRXRp/mIqiKEq5E2aEshDAJQDuBdCFmS9i5k6QR/SuAPBT\nIvq7GGNUFEVR6gFhnocygZlPZC9k5o8A/A3A34iownpkiqIoSr2i4AjFFBMi+m9KPzYxaBtFURTl\n3CXKTfmDAFLp57qDiG4gohfiCUtRFEWpb4QuKMz8bwCeALAoXUj+GcC3w+5PRI8R0S4i2hSwnojo\nASLaQkQbiGhQxrpqIno9vS70eybFz34GLFxYd9nChbL8XIujHGIojyDKI4xyiKFc4iiHGMopDisw\nc6gvAOMhN+gXAXgdwOVh903vfx2AQQA2BayfDGAWAAIwDMDK9PLGAN4CcDGASgAvA+gX5j0HDx7M\nSbBgAXOHDvI91+/nUhzlEEN5BFEeYZRDDOUSRznEUE5xGACs4Qjn9syvKAVhAYCR6Z/7A3gJwLhI\nbwb0zFNQfgPg5ozfXwdwAYDhAOZkLL8XwL1h3i+pgsIsB8P55zNPm5bswTF9OnPTpszf+EZycSxY\nwNyqFfOnP51gLhYsYG7fnnnIEPme0D/E5OLGG5P9f5RBKnjBAuaWLZk/97lkc9GuHfPQocl+Tk0u\nbr452TiYHRWUs3aUk/3yiPvkKyjTTcFK/z4fQBWAzwD4XcbyWwE8mOc97gKwBsCa7t2728xzJE6f\nZm7bVjJ8772JhcG//73EADB/97vJxHDqlHxYAObvfCeZGJiZubpagpg8ObEQTpyQAg8w/9u/JRYG\njx8vMXhecjEcOcJcUZHsscnMfN11EsONNyYXw4EDzI0aJZ8L5tIKSpjGxqCZXe+nZbDAbZKAmR9l\n5ipmrurYsWNicTz+OLBvn/z88MNna6Su+P3v5Xv37sD//E8ycTzyCHDokPz84IMJ5WLhQmD+fPl5\n3rzE/iG//jVw9Kj8nFQuFi4EFi+Wn2fNSu7YvP9+4ER6fmiSuXghPbVo+vTkcvGLXwCnT8vPDz2U\nXBwlU6jiQO6Z/AOA7lnLKwGMA/C/AP4+TPXCOSJ5LVjA3Ly5XG00bcpcU5PMMHbmTImhslKufp55\nJhmdulkzef+KioTkDaPxNGokwTRpkojWY3LRpImEcuutyfw/2rWT46JZMzk+k7pv0LSpHJsA8513\nJp+Lli2Ty8V559Xm4mtfq7/3UMLM8qoGcArAE0T0XtpyZRuANwHcDOBXzPx4yZUNSAG4LT3baxiA\n/SyjoNUAehNRLyKqBHBTetuyZfVqGREMGQJMmgSsXQv85S+y3CVPPinff/QjufrZvx946im3caxe\nDXTuDFx3HTBuHLB+fQK5WL0auPtuScL3vw+cPAncc4/zf8iqVUDr1sDEicC11wIbNiTz//jSl+Tn\n739fRkv//M/uj81Vq4BmzQDPAwYNAjZtSiYXX/iC/HzffTKKvvde97lYsQKoqABuvhno0wd47TX3\nubBGlOoDoAIyamgTtXJBphy/D+AEgB0AvgzgbgB3p9cTgIcgM7o2AqjK2HcygDfS674T9j2TGqHs\n3ClXGv/+78yPPy4/r17tPo7bb2du3Zr52DHmrl3lprhrtmyRv//++5kfflh+fvVV93Hw5z/P3KkT\n8/HjMjq55RbnIWzYIH//b37D/F//JT9v2+Y8DJ42jblHD+ajR+Wq/M473cewYoX8/X/8I/MPfsBM\nxPzBB+7jGD+euW9f5oMHZZTw9a+7j2HBAsnF008zf+tbMoLdt899HAa4vCkP4IcA/gPA5wBcVuwb\nu/hKqqA88ohkduNG5t27Rd5wfQP25Enmjh2Zb7pJfv/qV0WG++QTt3H84heSi7feYn73Xfn5Jz9x\nGwMfOyZT7r70Jfn9i19kbtNGiotDfvQj+ft37mR+4w35+b//22kIfOiQSE3/8A/y+2c+w9yli0yc\ncMm99zI3bsy8dy/z+vWSi9/+1m0M+/bJyftb35LfJ09m7tlTJtS45Otfl2J28CDzCy9ILv70J7cx\nZFJKQYlkX5+WnH4C4FcA9gP4NBH9trixUcMllQIuvhi44gqgQweRN1KORboVK4Ddu0VSAOT7J58A\nCxa4jSOVAq68UvLRrRsweLD7XGDxYuDAgbrJ+PhjYNkyp2GkUiKDXngh0Ls30Lev+1w8/7zIXJmp\n+OAD9/JKKiUyaLt2wIABQI8e7nMxa5aon5m52L5d5DdXMMvfPWEC0LIlMHQo0KlTAp8RS4QuKET0\ndYhk9SaAxQB6M/NPmfnOuIKrjxw6JJOJPA8wc988T/TybdvcxeH7ostOmiS/jxkDtGoly12xdy+w\ndGntBxaQn1esAD780F0c8H0R7CdMkN8nTgSaNnWajPfek5N2di4WL5ba5grfB9q0kZM5AEyeDDRu\n7Pa4eOst4JVXanNBBNTUyOS7w4fdxeH7cn9v6FD5fdq02uWu2LRJzgsmF40bSxyzZgHHj7uLwxZh\npg3/NxF9EcDXAfRl5m6Qrvd+RPTDuAOsb8yZAxw7dvaJA3B71eH7UkRat5bfzztPistzz9VOT4yb\nGTPkvbJzwSxxOMFcAk6cCDRvLstatJDi4vuy3gHmf5+di5MngZkznYSAU6dkauzkyXKxAcgI4brr\n3J5EzXtl5+LoUSkqLjh+XE7a06YBjdJnwQsukOLiOhdEtcUMkFwcOAAsWuQuDltEeR5KBwDLiWgd\ngJ9DbpDfTERtY4yv3uH78iG99traZZdeCvTr5+5A3bwZeOONuh9YoFbeWLXKTRy+L/LO4MG1y/r3\nB3r2dPihXb8eePfd3MnYvh3YuNFJGL4PXHKJHAeGIUOALl3c5eLFF+vKoAbPA159FdiyxU0cvg9c\ndZUcB4brrpORk6tcLFpUVwU1eB6wZg2wc6ebOHxfiliXLrXLJkyQax+Xhc0WYezrn2Xm70EepuUB\nmADgcQAnAbQDsICI3oozyPrCyZNyVT51KtAk60kzngcsWQJ89FH8cZgDsaam7vJJkyQuFwfq0aMy\nWqupqb0CBORqzPNEy3cib5hLwKlT6y6fNk2WO0jGwYNy7ypTBgUkL0beOHYs9jDOyKDV1XWXm+PE\nxXGxZ4/cuso+kVdUyMhp+nQZScWN78tJe/z4ustdqgk7d0rxys5Fs2YyoE6lnA2grRHlpvw9AP4I\n4BcQk8crAWxk5qsB9I0htnrHsmVSMLJP5IAsO3XKjbzh+8DVVwMXXVR3edu27uSN+fOlYATl4uhR\nYO7c+OOA7wMjRgDZrglGPHeQjNmzRWIJysXBg/HLG8zyp44dC5x/ft11vXrJyNHFcWFk0KBc7NkD\nLF8ebwyZKmizZnXX9e0rI0kXuTBFKygXO3YA69bFH4dNotjXvwlgKIC/AmgKYAOAT6fX1cPbR/ZJ\npeRexQ03nL3OyBtxX/l8+KHc9M6+6jF4njROvflmvHGkUjJrZdy4s9eNGiXyRuxXgW+/Dbz8cv5k\nrF0rn9wYSaXOlkEN48fLlXLcuXj9dfmf50vFCy/ICT1OUimga9e6Mqhh0iQZqcSdi/Xr5V+eKxdm\nBL1ggUhiceL7Iof3zXE5PnWqjGDr22yvSNOGmfk4M89g5v9g5l8z8764AqtvmCvA8ePlRJqNK3lj\n+nSJJd+JA4j3QD19Wl6/uloKbDYVFcCUKQ7kjVx3wjMxy2OcIXDiRLAMCsgV8g03xC9vBMmgBs+T\n/9uMGfHFkCmD5nL/O/98GUHFPVfC9+XzmK2CGjxP/m9z5sQXw4EDuWVQQ8eOMrCub/dRIhUUJZhX\nXgG2bg0+dwGy7tCheI3ffF/m9A8YkHu9WRfngbp6tdz8L5SL2OUN3xcvi8suy72+Tx9pCIkxGcuW\niUlooVzELW/4vowKunXLvX7wYBk5xHlcGBm0UC7efFMmlsSF78tosUOH3OtHjADat483F3PmSNEq\nlIuXX5a5I/UFLSiWMAdf5vS/bMaPlxmrcR2ohw/LtMugK0BD3PKG78t8+smTg7eprgYqK2P80H78\nsTR55PvEOtA3fF9GaRMnBm8zZYpcMceVi0IyKFDbCzJnDnDkSDxx+L70Qo0ZE7xN3BMEtm/Pr4IC\nMpKcOlVGa8YN2Ta+LwVtxIjgbZJoNygVLSiWMNP/LrggeJumTWvljTh6QebNq9sFHYSRN6ZPtx8D\nILkwXdBBtGol91dikzdmzqzbBh2E0Tdmz7YegpFBTRd0EB06ACNHxncSfe65/DKowbgpGJd/m5w+\nLXFMmpRbBjUYN4W4cpHvRngmxkxh6VL7MWTKoI0bB29n3BTqk+ylBcUCpgu60EEKyDbvvSf3gm3j\n+9LIaLqgg7j6avngxnGgbtkiPQ1hc7Fli0wSsI7vi4fFkCH5txs+XM7oMSRj40a5Ig6bi7jcFIwM\n2r9//u3idFNYtUpk0LC5WLlStreN78tJunfv/Ntdf70UvjhysXSpFKuwuVi8uPbZSuWOFhQLFLr3\nm4mRN2wPY3N1QQdh5I25c+3LG1FyYT5Q1of0mW3Q+S4BAVk/daqMaCzrG+bvyieDGuKaH3D4sPT8\nBN38zeS880SKjMNNIZUqLIMajJuC7RH0vn2FVVBDy5YysoxjsoTvi1qRTwY1eJ58tmfNshtDXGhB\nsUAqdXYXdBBxyRsvvij3RMJ8WID45A3flyvhXr0Kb9u1K1BVFcNV4KJF0twRJRkx6BthZFBDXG4K\nYWVQg+fJPRfbbgq+D4weLb1QhbjqKhlR2c7FrFlyco6SC9tmCplmkC1aFN5+6FBpmaovspcWlBI5\nePBsM8hCeJ4cpDbljWwzyEKMGSPTNG0eqEFd0PnwvBjkDdMGbcwgC3H99dbNIoO6oPNhzCJtyhvG\nDHLUqHDbx2EWaWTQsLmIy03B96UXrJAKaojDTMHIoGFz4dpNoVS0oJTInDmisEQ9cQD2DtR8XdBB\nVFbaN4vMZQZZCOtmkfnaoINo0UKKisUZAlGkP4ORN2y5KRgZdMqUwjKooW1bGUnYPIkW6oHJhTGL\ntOWmcOzY2WaQhejSxb6ZQi4zyEJ4nhs3BRtoQSkR35c56/mm/2VzySXyrBRbB+rmzfm7oIMw8sbK\nlXbi8P3gLuggrrxS5DFrH9p164LboPNRUyOd9Rs2WAkjXxd0ENdcY9cscvlyGTVGOZED9t0UcplB\nFsK4KdjKRVQV1GDbTMH3gWHDRMYKi3FTqA+ylxaUEijUBZ0PzxPJ3oZZZDFXgIBds8gjR/J3QQeR\nKW8cOlR6HGfaoKdMibafRX2jUBd0ELbljSAzyELY7AXZs0d6nqKeyDPdFE6eLD2OIDPIQtjsBdmx\nQ4pT1Fy4clOwgdOCQkTVRPQ6EW0hom/nWP+vRPRS+msTEZ0ionbpdduJaGN63RqXcQdhuqCjnsgB\nu2aRqRQwaFBwF3QQbdrYkzcWLJCb/MXm4tgxS8/CSKVym0EWonNnuXS0cOYwXdDF5MK4KZQqbxgZ\ndNy48DKooWdPGVHYOIkWI4MaamrkIW0vvlhaDEYFveEGuVUWhT59ZKRpIxdhe2ByUVMj9+XK3SzS\nWUEhosYAHgIwCUA/yLNU6syLYuafM/NAZh4I4F4Ai5k58xp+bHp9lau482HMIMNM/8vmmmtk9k+p\nJ/MwXdD58Lza56eUgumCHjs2+r6jRol2X3JhC9MGnQ9L+kYxMqjBlpvC5s1yM7yUVNhwU/B9udAZ\nNCj6vtXVMlIpNRfr1snJuJhc2DRTSKWk/6VPn+j7GrPIcpe9XI5QhgDYwsxb0+7ET0KerxLEzQCe\ncBJZEYTtgg7CyBuzZ5cmb4Ttgg7CRi+I6YIOMoMsRJMmluSNYu6EZ2JB3yhFBgXquimUIm8UK4Ma\nbLgpFCuDGs4/346bQrEqqMGGmUKxMqihQwfxH9OCUktXAO9m/L4jvewsiKg5gGoAf8tYzACeJ6K1\nRHRX0Juom17OAAAgAElEQVQQ0V1EtIaI1uzevdtC2LnJfhZ0MRh5Y8GC4l/DdEFfdVVx+/foAQwc\nWNqBarqgS83F3r0lmkWGbYMOwhhJlpAM0wVdai527izNTcH3pcena85PWGGMhFrKcTF/vsigpeai\nVDcF35feryAzyEKMGFG6mcLs2YXNIAvhefG5KdiiXG/KTwPwQpbcNTIthU0CcA8R5TQYYeZHmbmK\nmas6RtXRI1DM9L9sxo0rTd6I0gWdD8+TE3mx9TeMGWQhbrihRLPIKG3Q+aipETvo/fuL2j1KF3QQ\nkyeXJm988IHM3Ct2dALYcVMwMujo0cXHUeoEgW3b5CRcymFhw0zB9+W23vDhxcdRH8wiXRaUnQAy\nnyHYLb0sFzchS+5i5p3p77sAPAOR0BIj17Ogo9K0qchExZpFzp0brQs6iFLljShd0EG0aiX3D4qW\nN2bOjNYGHUQJ+kamDBqmCzqIUt0USpVBDcZN4fnno+8b1gyyEKW6KZSqghqMmcKSJdH3PXFCDs9C\nZpCFMG4KWlCE1QB6E1EvIqqEFI2zUkNErQGMBuBnLGtBRK3MzwAmAtjkJOocmC7oUq4ADTU1wPvv\nFydvpFLRuqCDGDhQHhdczIf2zTdFjrCVi7feKlLeSKVkplbYNugghg+XS8kiPrUbN0ori41clOKm\nkErJTK1CZpCFGD1aCn0xJ7BVq2TCSKkncqA0s8hUSk7Cl1xSWgzGTKGYXCxZEt4MshDlbhbprKAw\n80kAXwMwB8BrAJ5i5leI6G4iujtj008DmMvMmaYLnQEsI6KXAawCMIOZ7fuNh8R0ddv4sEyZUpzV\nRRQzyEJkyhuffBJtX1tXgEAJ8kYxbdBBlKBv2JBBDcXKG7ZkUEBGFsW6Kfi+TEooRQY1FGucaUsF\nBWTEOWFCcSNoI4Nef33pcdh2U7AOMzfYr8GDB3McVFczX3op8+nTdl5v9GjmK6+Mts+SJcwA81/+\nYieGuXPl9VKpaPuNGsV81VV2YmBmvuYa5qFDI+40e7YE/9xzdoJ49ll5veefj7Tb4MHMw4bZCYGZ\n+YormMeOjbbP009L6AsW2InhT3+S11u+PNp+ffsyjx9vJ4bTp5l79mSeMiXafn/8o8S+YoWdOH77\nW3m9l14Kv8/p08w9ejBPm2YnhlOnmLt0Yf7sZ+28Xi4ArOEiz7nlelO+bDl4sLTpf7nwPJk1tnVr\n+H2K7YIOYvTo6GaRxXZB58OYRb7/foSdim2DDqIIs8hiu6Dz4Xkil0RxU/B9uZc1cqSdGCZPju6m\nYGRQW7ko1k3B96XX65pr7MRRjJnChg0ig9rKRbmbRWpBicjs2dHNIAsR1SyylC7oICor5eTx3HMy\npA7D9OnFd0EHEVneyGyDDmsGWYjmzSObRdqU/gxR3RROnrQngxqKcVMotQcmF54nJ9CwZpE2VVCD\nMVOImgsiUVFtYdoNFi6095q20IISkVK6oIO4+GIxSQx7oL72Wmld0EF4HrBrV3izSGMGWUwXdBBX\nXCH5CP2hXbu2+DbofHge8M470nkfAt8vvgs6iKhuCsuXSy9PHKmI4qbg+8CAAdLjZIuRI6O5KSxc\nKCfdOHKxbh3w7ruFtwUk3uHDo5lBFsKYRZbjbC8tKBEI+yzoYqipEW+wvXsLbxvlSYBRqK4OL28c\nOSJXi8V2QQdhJgjMnx9S3kilSmuDDmLqVAkmxKf2wAE5gdnORVQ3hVRKRpq2ZFBDFDeF3bulsNk+\nkVdUyMhrxoxwbgqplNxIHzfObhwmF2FG0O++K8XH5kgNsOemEAdaUCKwbFnpXdBBRJm94ftiER/V\nDLIQbdrIg7fCFBQbXdBBRJI3fF88KYptgw6ic2e5tAyRDBtd0EGElTcyZdBWrezG0KOHjDjCHBel\nmEEWIqybQilmkIXo00dGomFyEYcMarDhphAHWlAiYKMLOoiqqnDyhumCjuMgBeR1X39dvvJhuqDH\njLEfw8iRQLt2IT60Ntqg8xFS3/B9qWc2ZVBDWDeFuGRQQ1g3Bd+Xnqarr7YfQ3V1ODeFuFRQoHaC\nQBgzhVRKnHxsyqCGKVPK0yxSC0pIbHVBB9GokQyNZ8+W7vcgbHVBBxFG3rDVBR1EaLPIOC8BM183\nTzJsdUEHEdZNwZxYbMughjBuCnHJoIZWrcKZRRorINsqqCGMmcL+/VJ04jo0S3VTiAstKCGJ+izo\nYvA8aUzLZxbp+3a6oIPo3l2uLvMdqCtX2uuCDsLzZLrsCy/k2cj3pQ360kvjCeLyywuaRZou6Dhz\nUVMDvPdefnmjVDPIQlx9dWE3heefL/6ZOGHxPHFTePXV4G2MGWT79vHEYMwU8uUiThnUUIqbQlxo\nQQmJzS7oIMaNEyv8oAP10CF7XdD5MPLGrl2515su6EmT4ouhoFnkRx/J2TzOTywgr79oUaC+YbML\nOohCbgrvvx+vDAqEc1PwfZnGHocMajCfv6BcbNsmJ9k4cxHGTMGYQQ4bFl8c5WgWqQUlJMYM0ub0\nv2zOO09OpEFWF/Pmyc3qOK8AAXl95mB5I5UCrruuNDPIQrRsKdMjA2eyzJolsxjiTkYefcPc/I1L\nBjW0by9X3EEnDvN/iru21tSIrDV//tnrMmXQysr4YjBmkUG5KOWpiFGoqZFrjFxmkXHLoIZLLpEB\nejnJXlpQQhBHF3QQnidXnGtyPOTY9+2YQRZi4ECRvnIdqLa7oPORV97wfbF6LtUMshDDhgXqG7a7\noPORT97wfaBXL+llipMxY4LdFFaulBGtq1wEuSn4vvQylWoGWYh8ZgqLF0uxcZWLqG4KcaIFJQQ2\nzSALESRvmC7oKVPsdUEHYeSNefPOljdMXC5yEShvxNEGHUTjxvI+OfQNFzKoIchNwZUMCsjIw5hF\nZrspuJBBDUFuCq5UUEBGpEFmCqmUmDbEKYMays0sUgtKCOLogg6iXTsZgWSfOOLqgg7C80TeyH4W\nRhxd0EFceKEMQM4qKHG1QQfheXLJuXhxncW+LwOYOGVQQ5Cbwty5Ul9dpiKXm4LvywimTZv4Y7jy\nShmRZefC1iNxwpLLTMHMBr3+eulmj5uobgpxowWlAKU+C7oYPA945RWRewy+H08XdBCjRwOtW9c9\nUOPqgs6H58mzNd57L2Oh78sloi0zyEJMmCCXnBnJMF3QLnNRUyOPGM6UN2ybQRZi0qSz3RTeeEOs\nWeK+b2EwvSDZbgrGDLKqyk0cxkwhMxcvvyxFxtVxEdVNIfZ4kg6g3HEx/S+b7OeCxNkFHYSxusiU\nN+IwgyzEWfLG6dPxtUEHkcMsMu4WmFwYeWPGDPk9UwZt0sRNDLncFOIwgyyEcVOYM0d+P3ZMPqs1\nNfGroIZcZgpxmEEWopzMIrWgFMB0QZfyLOioGHnDnLRee01GKy4/sIC83+7dwIoV8nsqJXYvcXRB\nB9Gvn+TjzKyedetkuOI6GZ4nw5K0vpFKiQx6+eXuQjBuCiYXy5fLaMVlUQMk9ZluCqmUTORwIYMa\njFmkyYVRQZP4jKxfX2umYMwgO3VyF0NYNwUXaEHJg6vpf7nwPJE39u5N5goQEHmjokLeP+4u6CDO\nkjd8Px4zyEJk6BuZXdAuc5HppnDsWK0MesMN7mIA6ropJCGDAme7KRgV1LYZZCEye0HeeUeKi+tc\nZJpFRn2ypm20oORh6dL4u6CDMFYXM2bE3wUdROvWtfKG6YJOKhdn5A3TBm3bDLIQnTqJWZfvJyKD\nGoy8MX++pGL8eHcyqKFHDxmR+H4yMqjBuCksXSon0+pqdyqooU+fWjOFJGRQg+cVdlNwgdOCQkTV\nRPQ6EW0hom/nWD+GiPYT0Uvpr++F3TcOXHRBBzF4sMxyevLJ+Lug8+F5ctP1D3+Ivws6iGuvldlv\nL/zRQRt0PtL6hv/nw+jY0a0Mahg7Vq7EH388GRnUYNwUnnpKLFkGDnQfg3FT+O1v5WSa5Gdk0SLg\n6adFAnUpgxoKuSk4o9hnB0f9AtAYwFsALgZQCeBlAP2ythkDYHox++b6KuWZ8rafBV0Md9/NXFkp\nz7HesCGZGN55R96/ZUvmz38+mRiYmW+7jfne5vdLMFu2JBPE5s18HE24dbOjfPvtyYTAzHzjjcyt\nWkkqdu5MJoZ16+T9KyuZ77knmRiYmSdNYm7ThrlxY+a9e5OJYdkyyUXjxszf/GYyMTAzjx7NfOWV\npb8O6skz5YcA2MLMW5n5OIAnAYS9pihl36Jw2QUdhOfJ44Y7d46/CzqIiy6SIb3Lto9ceB5w/Sc+\nDvd00AYdxOWXY3G3v8P+I+clnouDB2XCwoUXJhPDwIGiOtp+HHZUPE9k6UGDZBSbBMOGiezosgcm\nF54HbNoEbN2aXAwuC0pXAJkPltiRXpbNCCLaQESziOiKiPuCiO4iojVEtGZ3oYc35CGJ6X/ZXHON\nfO/Y0e3N32w6dpTvgwcnF8PEqo8wCkvxYseENJ40qU53oCmOYMI1BR6GESPm/2D+L0lAVPv+5jhN\ngquuku9J5qJxYymuRG5nQGYT5cmacVFuN+XXAejOzFcB+DWAZ6O+ADM/ysxVzFzVsYSjzGUXdBCm\nMXvnzmRnb5imwqVLk4uh5ZKZaIJTeHinl9hjT5kBf2cVrsc8tFgyK5kgIE8OBeS4SIpTp2q9tHIZ\nJLrCdOzv2JFcDMePy+McmAs8biFmLrlEfMySvI/isqDsBHBRxu/d0svOwMwHmPlQ+ueZACqIqEOY\nfW2yY4f7LuhcmKmQ+/YBq1cnE8Mbb4ghYdu2Cd/w830cbn0Bnn3vGrzySjIhvPwy8M6H58FrtTDR\nZPi+OBBv2ZKcvLFypUhNWQYCzvF9uejbsCHLTcEhixfLDMgwT5OMG9NukJRZpMuCshpAbyLqRUSV\nAG4CUGdwRkRdiETcIaIh6fj2htnXJklO/zOYLujJk5OdvZHZA5PLLNIJpg162jQwGiWaizMy6MyZ\ncmnqGDNlONtNwTWZZpC5zCJdYKYMm1xkm0W6wphBTphQ+GmScZPtpuAaZwWFmU8C+BqAOQBeA/AU\nM79CRHcT0d3pzT4DYBMRvQzgAQA3pSce5Nw3rlh9P75nQYflhRfkA/O5z8mzR5I8cQwcCNx6qzya\neN68BIJYsAA4dAgtbvYwdGiyuRg+HOh801gxecsyi3TBnDlSX2+7LbdZpCuMGeRnP1vXTcElM2bI\nyfOOO8RNIYlcGFukiROBG2+UjvmXXnIfh8G4KSQ2Uip2elh9+Cpm2vDHHzNXVDD/679G3tUq//RP\nMiXzwAHmX/1KpiW++abbGD78kJmI+b77mI8fZ27dmpOZLvuVrzC3aMF85Aj/x38kM1327bflfX/6\nU2Y+fJi5WbNE5sveeitzu3bMJ04wf+c7zI0aMe/Z4zaGzZslF7/+tXxemjRJ5vNy443MF17IfOpU\n3c+LS8z06cceq/t5SZKMj0tRoJ5MG64XmC7opBrGgNonARozyKRmb8yYIbHU1NSaRU6f7ljeMI8C\nTJtBJiVvmPerqYGYRU6cmOdxkvFw8qT8T4wZpHFTcP0sjMynIho3BdfH5tGjZ1TQM5Y0x4+LPZBL\njAw6ZYqYKQwfnvwjeWtqgMOHkzGL1IKShXkWdBJd0IZXX5UuaHMPp1cvoH9/98NY35c+FDMV0vNE\n3njxRYdBrF1bpw26Xz+ZzZJELurIoMYs0qG+YWRQc1wYN4UkcnH11fJUT0DiyTSLdMHChXLSNLkY\nOVL6UJLIxYgRtWaQnid+Xu+84zaOTJI0i9SCkkGSZpCZ5DKD9DyZLrpnj5sYPvnkbDPITLNIZ/i+\n/DPSZpCZZpEHD7oJYf9+sdaoM0kj18MwYibbDDLTLPLoUTcx7Np1thlkEhMEfB9o2bLWDNKYRc6Y\nISM5F7z9tlxPZOYi0ywyKZo2FV+zJMwitaBksGSJu2dB58P3pVkssws60yzSBc8/Lw7Dmbk4/3zx\nknL6YTFmkO3bn1lkHATMszDiZtasHGaQHTueMYt0gbn5O368nEgNRt5YsMBJGJg+vVYGNXTvLiMW\nVwXFPBKnuho477za5cYs0vTpxE2uR4MbL69ymD78/vvAmjVu31cLSga+7+5Z0EG89548oTC7qA0e\nLG7Drg5U35cCMnp03eXGLHLzZgdBbN0qXhJZyRgxQuqLy1x07CiNrnXwPLlEffvt2GN45RVJR/Zx\nMW6cFBiXueje/WwzSM8TKfTDD+OPYc0aOVlm5+KGG6TAuMyFcRvOxJhFfvyxmzhykZRZpBaUNOYK\n0NWzoIPIddUDiLpSUyNX5UeOxBvDqVMSx+TJIrFk4lTeMG+SlYwmTURxmjFDRg5xcvy4yKDTpuWQ\nQR3qGyYV06bVXX7eee7kjU8+kWnjuZ6J43nyGZo+Pd4YgFoVdPLkustbtpQRnItekI8/ziGDpvE8\nkd1mJWemgHbtgFGj3EtvWlDSuH4WdBCplMypv+KKs9d5nnyo45Y3Vq6Um++5ctGtm4yWnBSUVEqa\nLS6++KxVnicOAnFbXSxZIi0nOY8Lc5fewac2lQKGDMltBul5wAcfxC9v5JJBDQMGyMjFxQkslZLe\nrFxmkJ4nzg5xuynMmiVFI1cuhg6Vm/RJz/ZKwixSC0qazOl/SZHZBZ3LDHLMGJlGHPfJ3HRBV1fn\nXl9TI41sscob2W3QWVx/vRt5w8igEyYEbOBA3wiSQQ2u3BR8X6YJZ8ugQO0IOm43BaOCBk3rN2au\nLnLRqZMU+WwaN5Y4EjJTOEMSkyW0oKRJpdJd0AmaQZou6KATh5E3nnsuXnnDdEG3aZN7vZE3Yu0F\nMW3QAclo2TJ+q4tQMqgDfSNIBjW0axe/m0KmDFpRkXsbz5MRTJxuCgEq6BkuvFBO8nHm4vhx+Xfn\nlEFRG19CZgpnuPhi924KWlAg7QTlYgbZrp1MagrCyBurVsUTg+knyJeLq66Sx8DGeqD6vnhIVFUF\nbmLkjU2b4gnhpZfk2Mh7XAwZIpeqMSbD9+Xk0K9f8DY1NSLzvPVWPDGsWCEyaL6G39GjZQQT53Hh\n+9KT1atX8DaeJ2aqcZlFLlqURwZNM2FC8saZQK1Z5N69bt5PCwrKxwwysws6iLjljVw9MNmYXpDn\nn5cpq9YxbdA1NdJsEcC0afG2goR6Jk7jxhLIrFmx6BsHD4oM6nn5n4ljjt04c1FRIb1IQcTtprB3\nr5wcC31O454rkUrJiDVQBkViZgpn4brdQAsK5MOS1LOgDcuW1e2CDqJtW7kSjPPEkdkFHYTnxWgW\nuWBB3TboALp0Qaxmkdld0IEYfWPRIusxzJkT7qmIcbspGBm0dev828XppjBjhpwcC+UiTjcFY4s0\ncaKMQPJhzBTWr7cfR1iMm4KrCQLnfEHJ2QWdAKlU3S7ofHge8NprwJtv2o1h1y45EYTJxahRco8l\nlhNYKlW3DToPniezm2w/bOqdd87ugg5kwgS5JI3hU5tKiQx67bWFtzVuCrbljddfl96jMLmorpaR\nShwnsFRKTo6FnhxqRtALFth3U1i/PoQMmsaYKST6BEXHbgrnfEExXdBJm0Hm6oIOIi6zyFxd0EHE\nJm+YNmjTpVaAuMwiMw0QC9KsWSz6hpFBp07NL4Ma4pI3wsigBmMWaftCI1MFDfM4bGMWadtNwffl\nJB1mNqhjM4VAXLopnPMFJZUK6IJ2SFAXdBA9e8qNcdsHalAXdBCeJ95iy5dbDCKoDTqAvn2BSy+N\nJxeRZNAY9I2wMqghLrNI3wcGDRKj0DDE4aYQUgU9w7XXxmMWaWTQsE8Xd2imEIhLN4VzuqAYM8h8\n0/9cENQFnQ/Pk6Y+W2aR+bqggzDyhtUDNcsMshBxyBv5uqADmTJFLl0tJsP3ZZA2cWK47TPdFGzJ\nGx9+KDJolBF8HP0Pvi89WGPHhts+DjeF7dulATrKcVEOZpEu3RTO6YKyeHH5mEEGdUEHYeQNW1YX\n8+YFd0EHcf75cvVjtRfE9+UGTa426ACMWeTs2XZCyNcFHYhlfSOqDGrwPLmSnz/fShhnZNAoubjo\nIhnR2CooQWaQhTBuCrbMIgv1A+XCmCkkLXu5clM4pwtKwS5oB7z3nsyZj1rUBg0SGxRbB2q+Luh8\neB6wZYsleeOtt0T/i5iMESOADh3s5qJTJ5lBFgnPk0vY7dtLjmHTJumxiXpcjB1r103B96XnaMCA\naPt5nj03hdWr5WQYNRcTJ9p1U/B9kVh79462n+fJxWuSZpGu3BScFhQiqiai14loCxF9O8f6W4ho\nAxFtJKLlRDQgY9329PKXiKjkOpv5LOhyNIMshJE35s4t3Szy1Cm5Es3XBR2EkemsHKiR7oTXYqwu\nbMgbYbqgAzH/RAszBEwqosiggF03hWJkUINNs8hUKrcZZCGMm4KNuRIffyxFoRg1o1zMIuN2UwAc\nFhQiagzgIQCTAPQDcDMRZff+bgMwmpn7A/gRgEez1o9l5oHMHNw+HRLTBZ3k7C4gXBd0EDU18qEv\nVd4I0wUdhFWzSN8PNIMsRE2NfOhLlTcWL5aWkqKOi9695RLWQjJ8X0ZIF1wQfV8jb6xeXVoM8+bJ\nvZhiTqI23RR8X06GbdtG37emxo6bwsyZUhSKOS4cmCmEIm43BcDtCGUIgC3MvJWZjwN4EkCdQ5WZ\nlzPzvvSvKwB0iyuYUF3QMRO2CzoIW2aRYbqg8+F54lD8wQclBBG2DTqAiRPlSXU2clGSDFpTU7K+\nYWTQYi92Jk2yI28YGfS666Lvm2kWWYqbQpEq6BlsjaB9X3z+IsugiN1MITRxuykAbgtKVwDvZvy+\nI70siC8DyBwkMoDniWgtEd1VajCpVMgu6BgJ2wUdxHnnycmjVHkjbBd0EFbMIsO2QQfQokXpZpGZ\nXdBFy6BG35g5s8gXKN0KyIa8UYoMarDhplDIDLIQF1xQupvCsWO1MmgeJ6C8xGimEJq43RSAMr0p\nT0RjIQXlWxmLRzLzQIhkdg8R5bxuIqK7iGgNEa3ZvXt3ztd/5x1pFyiH2V1hu6CD8Dy58blyZXH7\nb94cvgs6iP79pTempKmRvh+uDToPnif3wzduLG7/KF3QgQwdKpeyJXxqfV+sQ4qRQQ2eB7z6qkyY\nKIYXXwx+Jk5YrruudDcF3xf5rGfP4l+jVDeFxYtFTSglF8ZMIWnZKy43BYPLgrITQGZrVLf0sjoQ\n0VUAfgfAY+YzfzYz70x/3wXgGYiEdhbM/CgzVzFzVceA7qNyMIM8cSJaF3QQkyfL/sUeqFG6oIMo\n2Szy6FEZrhUwgyxEqWaRpgu6JBm0UaNafePYsci7HzwoPTXFyqCGUuWNUmVQoHQ3hT175ORX6ue0\n1F4Q35diMH588THEZKYQmbjNIl0WlNUAehNRLyKqBHATgDr/YiLqDuBpALcy8xsZy1sQUSvzM4CJ\nAIq+zRb0LGiXvPCCzJEv9cPSpo1M9S32w5JKReuCDsLIG3PnFrFz1DboADp3FseDUnIRpQs6EM+T\nylDEwzBKlUENxk2hlFyMHSu9RqVg3BSKMYucObMkFfQMxk2hmFwYGfSGGwqbQRbC84AdO8rDLDKu\nkZKzgsLMJwF8DcAcAK8BeIqZXyGiu4no7vRm3wPQHsDDWdODOwNYRkQvA1gFYAYzF9XGZrqgy2F2\nV5Qu6HzU1BRnFllMF3QQI0eWIG/4vszxDNsGnYeamuLkjbfflpl/Vo6L8eOL1jd8H2jfXgpbqRh5\nI6qbgg0Z1FCKm4LvA127ygVPKZgJAsW4KaxbJ0XAxnERg5lCZOJwU6gDMzfYr8GDB3M2f/4zM8D8\nwgtnrXLG6dPMvXoxT55s5/W2b5e/6ec/j7bf734n+61fbyeOW25hbt+e+cSJCDudOsXcpQvzZz5j\nJYZXX5W/6eGHo+33wAOy3+uvWwmD+VOfYu7WTf7ZITl+nLltW+bbbrMTwurV8jc9/ni0/f7zP2W/\nd96xE8fEicy9e0dKBR85wtyiBfNXv2onhsWL5W966qlo+333u8yNGjHv3m0njpEjmQcMsPNaxTJr\nluRi+vTc6wGs4SLPuWV5Uz5OUqkiu6AtUmwXdBCmkznqlU+xXdBBeJ7c7ItkFllsG3QAffpIO0gx\nubAqgxp9Y9260LssW2ZHBjUMHixX+MXkwoYMavA8GT1HcVOYP9+KCnqGESNk5FdMLq69VpwYbGDR\nTKFobLspZHJOFZTjx+uvGWQhPE9O5AET287i8OHiu6CDqK6WZ7pEOlCNGWTUNugAMs0iDxwIt08p\nXdCBTJ0aWd+wKYMCdeWNsG4KH34oja42c1GMWaQxgxwzxk4MxZhFbtsGbNhgNxflZBZpw00hm3Oq\noJgu6HKYLlxsF3QQUc0iS+mCDqJVqyLMIk0bdAQzyEJ4npw0wppFmi5oq8dFhw5yaRvyLGqsgCZM\niGYGWQjPi+am8Nxz0c0gCxHVTeH0aYlj0qRoZpCF8Dy5eFi6NNz2xdoi5cOimUJJ2HJTyOacKijl\nYAa5c6fcNLZd1K6+Wj64Ya98Uqniu6Dz4XnS3fzaayE23rJFmiUsJ2P4cDmfR8lFLDKo58klbgh9\nY9Mm2cz2cWHcFKLkokcPmSFmkyhuCpZV0DMYs8iwufB96QW69FK7cRizyH37Cm8bF3GZRZ4zBSXK\ns6DjxFz12J5lFsUs0kYXdBCRrC6KNIMsRBSzyEwzyBJaYHIT4dGaccigQDQ3BSODltoDk4uamvBm\nkUYFLaUHJhdR3BT27ZOTfhyzQWtq5DOYpFlk27bxmEWeMwXFShe0BWx0QQdh5I3nn8+/nY0u6CC6\ndgWqqkIeqL4vbfa9elmPI6y8sWhRjDJoBH3DyKBdutgPo6ZGrvhXrcq/nZFB4ziJRjGL9H3prSrG\nDLIQYd0UZs6Uk34cx4UFMwUrlOqmkItzpqCkUha6oEvEVhd0EGPGSCNaoQPVRhd0Poy88f77eTay\n1QYdwPXXhzOLNF3QscmgIfSNuGRQQ1h5w/ell8i2DAqEd1OISQU9Q1g3Bd+X4j4kpx9HaZRopmCN\nOMcCwS8AAAvFSURBVMwiz5mCEvVZ0HEwe7adLuggKitr5Y0gqwtz89dGF3QQoR4LUqIZZCFatJCi\nkk/ecCKDep78M/KYRcZtBdS2rVzx5ztxxCmDGsK4KZRqBlmILl0Km0XaMIMshDFTSNIsslQ3hVyc\nEwXFdEGXg9xlqws6CM8Ddu0KNovcvFl6AuLMxZVXioqV90A1bdAlmEEWwvPkf79hQ+71pgs61uNi\nyBA5i+U5g/m+3Pjt2ze+MDwvv5vC8uUyaIwzF6NGFXZT8H3pi+rRI744PA9Yu1b+97lYtAg4dCje\nXJRgpmCVYt0UgjgnCkpDMoMsxKRJ8vpBJ/OY7oPXIVPeOHQoxwaZZpBxaH9ppk6Vl8+Xi0aNxBIj\nNoy+MXt2Tn3jwIF4ZVBDof6HVEpGJtXV8cVQUSG5DjKL3LNHPO7i/pwWGkHbMIMsRLNm4g/W0Mwi\nz4mCYrqgoz4L2ibLlslN4rg9xIxZZNCVj+mC7hbbo8uEmho5f+aUN+bPl9kDMSfDmEXmy4UTGTSP\nvjFnjlxsxH0SzeemYGTQcePik0ENNTXBbgpGBY37M9Knj4wIg3JhzCCbNo03jpoauX8WwUzBOoMG\nFeemEESDLyixdEEXge0u6Hx4Xq3BXyZxdEEHMWqUaPc5D1TTBm3BDLIQQfLG9u1igeHkuBg3LlDf\nMDLo8OHxh1FTIyOAbHnDyKAuDFPzmUX6vlzolGoGWYh8bgpr18pJ3sVxUYSZgnWKcVPIR4MvKLNm\nxdAFHZG4uqCDCLK6iKMLOogmTWrljZMnM1aYNujqartt0AEEST1OZdAAfcOVDGoIclOw8UycsJx/\nfm43hSNHnKigZwhyUzDPxIlVBk0T0UwhNqK6KeSjwReUUp4FbYuNG+Ppgg6iRw9g4MCzD1RjBmm7\nCzoIzwM++kiuis+walU8bdABGMPHXLno29ehDOp5cum7du2ZRUuXygja1XFhpM5cuRg8OH4Z1OB5\nMj04003BqKCucjFihJzQc+Vi5Eh7ZpCFMGYK27a5eb9cRHVTyEeDLijM8U//C4Pvy1WX7S7ofBiz\nyF275PfDh+Umedw3fzO54QaZylznQLVsBhkGzwMWLgT275ffTRe001Frjodh+L7o9C5kUCC3m8IH\nH8iMQJe5yDWC9n0ZvdgygyyEcVOYObPWTWHbNrn4c5mLcjGLDOumUIgGXVAOHiwPM8hUKr4u6CA8\nTwqqmb0RhxlkIVq1kpkydeSNVCq+NugAjLwxZ478PmtWfF3QgXToIJe+6TOHufk7YYL0zLgiW96Y\nPt2dDGowbgrmJJppBllZ6S6ObDeFJGaDXnqpuGYkWVCAWrPIQm4KhWjQBeXjj+Of/lcI0wXt+gmR\nAwfK8yzMVaDpgh41ym0cNTViFvnqq6htg3acjGHDZCZXZi46d46nCzovGfqGaxnUYOSNzFz07CkO\nOC6pqak1i1y1SiaMuP6MZLspGDPISy5xG0dNTfJmkZMm2TGLbPAFJWkzyKR6YDLljYMH4++CDqKO\nvBF3G3QAmfLGoUMJyqAZZpFGBnVtBZTppnDwoMigrm6EZ2JG0M89J4dFkyZOVVAAdc0i9+4FlixJ\nRs0IYaYQO2HcFMLQoAuKi/n9hXDRBR2E54lW/uCD8XdBB3HhhcA112QUlKuukktixxh548EH5USa\nyHFh9A3fj9UMshCeJyOChx5yL4Ma+veXw8D3a80g27RxH4dxU/jNbxKQQdOEMFNwgnFTKIlinx1c\nzBeAagCvA9gC4Ns51hOAB9LrNwAYFHbf3O83mHfvZu7Zdh9X4FhxD2Augko6xpd13sf79zNXVDB/\n4xvMl3Xex5XkLgYTx/nYx1deKXHs3+8+jko6xhe22MftsZtPN2rE/N3vOo+BiLljR+ZmzZj792du\n3lx+J3IWAjMzN6dDPKr5an63UXcGmH/yE+YrOu/i5nTIaQx9Ou7iJk2Yr7ySuU0b5r6d3MZg4mjX\n7DBXVMjzzR94IJlcXNZxNxMx9+vHfMEFzP0cx2DiaNv0E27ZkvnoUVnmOhcVOMZdz/+Y03c7i36m\nvMti0hjAWwAuBlAJ4GUA/bK2mQxgVrqwDAOwMuy+ub5athzMPdvuY+A092y7z+5/IA+XdZb3vOj8\nfQzwmRgu6+wuhsw4zsc+njix9neXcZzJBbYyA9y37XvOY+jcWY70igrmxo35zEmsc2dnITCznCSA\n09wf6xlgvrT9XgZO8xWddzmPoXmTo9y4MfP5533iPIbMOOQExnxZhz2J5aJp42PcqBFz26bJ52L2\n7NrfXcZhzlEVjU6WVFCIHRnJENFwAN9n5hvSv98LAMz8k4xtfgNgETM/kf79dQBjAPQstG/u96xi\nwPIzLhVLOBbty4YEjZvOglBe8SSJ5qKWKjCvLeoD6qA/9wxdAbyb8fsOANnthrm26RpyXwAAEd0F\n4C75rT2AKgDr1ubaNn4GDa49UJOKoVziKIcYMq2N1yb2/yAMGswgEBicUC7KIYZyiaMcYiifOAYN\nBt4uem+XBcUJzPwogEcBgIjWMO+pSjikskBywed8LjQPtRDRmtOaCwCai0yIaE2x+7osKDsBXJTx\ne7f0sjDbVITYV1EURUkQl9OGVwPoTUS9iKgSwE0AsvtDUwBuI2EYgP3M/H7IfRVFUZQEcTZCYeaT\nRPQ1AHMgs7YeY+ZXiOju9PpHAMyEzPTaAuATALfn2zfE2z5q/y+pt2guBM1DLZqLWjQXtRSdC2ez\nvBRFUZSGTYPulFcURVHcoQVFURRFsUK9LyhEVE1ErxPRFiL6do71REQPpNdvIKKYHzCaHCFycUs6\nBxuJaDkRDUgiThcUykXGdtcQ0Uki+ozL+FwSJhdENIaIXiKiV4hosesYXRHiM9KaiJ4jopfTubg9\niTjjhogeI6JdRLQpYH1x581iW+zL4Qsl2Lk0tK+QuRgBoG3650nnci4ytlsAmQzymaTjTvC4aAPg\nVQDd0793SjruBHPxfwH8NP1zRwAfAahMOvYYcnEdgEEANgWsL+q8Wd9HKEMAbGHmrcx8HMCTALL9\nQj0Af2BhBYA2RHSB60AdUDAXzLycmc1TF1ZA+nkaImGOCwD4BwB/A7DLZXCOCZOLLwB4mpnfAQBm\nbqj5CJMLBtCKiAhAS0hBOek2zPhh5iWQvy2Ios6b9b2gBFm1RN2mIRD17/wy5AqkIVIwF0TUFcCn\nAfyPw7iSIMxxcRmAtkS0iIjWEtFtzqJzS5hcPAigL4D3AGwE8HVmLvHBuPWSos6bDc56RSkMEY2F\nFJSRSceSIL8C8C1mPk2uny5VfjQBMBjAeADNALxIRCuY+Y1kw0qEGwC8BGAcgEsAzCOipcx8INmw\n6gf1vaCUYufS0Aj1dxLRVQB+B2ASM+91FJtrwuSiCsCT6WLSAcBkIjrJzM+6CdEZYXKxA8BeZj4M\n4DARLQEwAEBDKyhhcnE7gP9kuZGwhYi2AegDoMSnrdc7ijpv1nfJqxQ7l4ZGwVwQUXcATwO4tYFf\nfRbMBTP3YuaezNwTwF8B/J8GWEyAcJ8RH8BIImpCRM0hTt6lPruvHAmTi3cgIzUQUWcAlwPY6jTK\n8qCo82a9HqFwCXYuDY2QufgexNP/4fSV+UlugA6rIXNxThAmF8z8GhHNhjwl9TSA3zFzzumk9ZmQ\nx8WPADxORBshM5y+xcx7Egs6JojoCcizpjoQ0Q4A90FMeEs6b6r1iqIoimKF+i55KYqiKGWCFhRF\nURTFClpQFEVRFCtoQVEURVGsoAVFURRFsYIWFEVRFMUKWlAURVEUK2hBURSHENFCIro+/fOPiejX\nScekKLao153yilIPuQ/AD4moE4CrAdQkHI+iWEM75RXFMeknIrYEMIaZDyYdj6LYQiUvRXEIEfUH\ncAGA41pMlIaGFhRFcUT6iXd/gjwN7xARVScckqJYRQuKojggbQv/NIBvMPNrEFfb+5KNSlHsovdQ\nFEVRFCvoCEVRFEWxghYURVEUxQpaUBRFURQraEFRFEVRrKAFRVEURbGCFhRFURTFClpQFEVRFCv8\n/xJGvTSaW3f6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x109a6d4a8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in range(1,n):\n",
    "    phi=(np.zeros(xs.shape))\n",
    "    phi[i] = 1\n",
    "    if i==4:\n",
    "        plt.plot(xs,phi,'r-x')\n",
    "    else:\n",
    "        plt.plot(xs,phi,'b-x')\n",
    "plt.title(r'Basis functions ')\n",
    "plt.text(4./n, 1.4, r'$\\phi_4$', fontdict={'color':'red','size':24})\n",
    "plt.xlabel(r'$x$')\n",
    "plt.axis([0,1,0,2])\n",
    "plt.ylabel(r'$\\phi_i(x)$')\n",
    "plt.show()\n",
    "# First and last basis function are omitted, since we know the value of u in the first and last node."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute the stiffness matrix and mass vector\n",
    "Recall the weak form:\n",
    "\n",
    "$$\\int_\\Omega\\nabla u\\cdot \\nabla\\phi = \\int_\\Omega f \\phi.$$\n",
    "\n",
    "We still have a fully continuous solution $u$. To discretize $u$, we use the assumption that we can express it in the very same basis functions: we will be looking for solutions of the form \n",
    "\n",
    "$$u = \\sum_{j=1}^{n-1}\\vec{u}_j\\phi_j.$$\n",
    "\n",
    "Since we establised basis functions earlier, we are essentally looking for the coefficient vector $\\vec{u}$.\n",
    "Substituting this in our system of equations gives us:\n",
    "For all $i$, find $u_i$ such that \n",
    "\n",
    "$$\\int_\\Omega \\sum_{j=1}^{n-1} u_j \\nabla \\phi_j\\nabla\\phi_i = \\int_\\Omega f \\phi_i.$$\n",
    "\n",
    "In the system above, all quantities except for $(u_i)_i$ are known; this can be expressed as a matrix-vector system:\n",
    "\n",
    "$$A\\vec{u} = \\vec{b}$$\n",
    "\n",
    "where matrix $A$ is defined as \n",
    "\n",
    "$$A_{ij} = \\int_\\Omega \\nabla \\phi_i\\nabla\\phi_j$$\n",
    "\n",
    "and vector $\\vec{b}$ is defined as \n",
    "\n",
    "$$\\vec{b}_i= \\int_\\Omega f\\phi_i.$$\n",
    "\n",
    "$A$ is called the _stiffness matrix_, $\\vec{b}$ is called the _load vector_. These terms refer to a classical case in engineering, where the finite element method is used to compute the strain on a beam for a given load.\n",
    "\n",
    "Because the basis functions (and their derivatives) only overlap with their direct neighbours, the stiffness matrix is _sparse_, it contains many zero entries. This is convenient, since many solution algorithms can take advantage of that fact to find a solution much faster than in case of full matrices. We won't make use of that here.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stiffness matrix\n",
    "We picked piecewise linear basis functions, so we evaluate piecewise constant gradients to compute $A$.\n",
    "In our one-dimensional case, we obtain the following stiffness matrix:\n",
    "\n",
    "$$A_{ij} = \\int \\nabla \\phi_i\\nabla\\phi_j= \\int_0^1 \\phi_i'(x)\\phi_j'(x)dx= \\begin{cases}2/h & \\mbox{ if } i=j\\\\-1/h& \\mbox{ if } i = j\\pm 1\\\\ 0 & \\mbox{ else }\\end{cases}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Stiffness matrix:\n",
      " [[ 20. -10.   0.   0.   0.   0.   0.   0.   0.]\n",
      " [-10.  20. -10.   0.   0.   0.   0.   0.   0.]\n",
      " [  0. -10.  20. -10.   0.   0.   0.   0.   0.]\n",
      " [  0.   0. -10.  20. -10.   0.   0.   0.   0.]\n",
      " [  0.   0.   0. -10.  20. -10.   0.   0.   0.]\n",
      " [  0.   0.   0.   0. -10.  20. -10.   0.   0.]\n",
      " [  0.   0.   0.   0.   0. -10.  20. -10.   0.]\n",
      " [  0.   0.   0.   0.   0.   0. -10.  20. -10.]\n",
      " [  0.   0.   0.   0.   0.   0.   0. -10.  20.]]\n"
     ]
    }
   ],
   "source": [
    "# Construct the tri-diagonal matrix\n",
    "A = (np.diag(2*np.ones(n-1)) + np.diag(-np.ones(n-2),1) + np.diag(-np.ones(n-2),-1))/h\n",
    "print(\"Stiffness matrix:\\n\",A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load vector\n",
    "While in this case it is possible to evaluate the expression for $\\vec{b}$ exactly, for general $f$ or more complicated basis functions this is not the case. For our model problem, we choose a rough but consistent approximation, which happens to be exact in case $f$ is a linear function:\n",
    "\n",
    "$$ \\vec{b}_i= \\int f\\phi_i= \\int_0^1 f(x) \\phi_i(x)dx \\approx hf\\left(x_i\\right).$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Element centers: [ 0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9]\n",
      "Load vector: [-0.14 -0.08 -0.02  0.04  0.1   0.16  0.22  0.28  0.34]\n"
     ]
    }
   ],
   "source": [
    "# Construct the right hand side vector\n",
    "# Centers of the elements\n",
    "centers = xs[1:-1]\n",
    "print(\"Element centers:\",centers)\n",
    "b = h*f(centers)\n",
    "print(\"Load vector:\",b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve the system\n",
    "Now that we have expressed our partial differential equation as a linear system, all that remains is to solve this system. Our system is quite small, but we note that for higher precision, the number of elements in the stiffness matrix increases drastically. In real life examples, the discretized domain can easily exceed a hundred thousand elements. For such large systems, people often have to resort to different solution methods. A popular class of algorithms is of the _iterative_ kind, where one repeats a sequence of steps, each time coming closer to the actual solution of the system. \n",
    "\n",
    "We have no such problems and we can solve our system simply by using Numpy's `linalg.solve` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finite element solution:  [ 0.009  0.032  0.063  0.096  0.125  0.144  0.147  0.128  0.081]\n"
     ]
    }
   ],
   "source": [
    "# Solve linear system\n",
    "u = np.linalg.solve(A,b)\n",
    "print(\"Finite element solution: \",u)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The vector $\\vec{u}$ contains the coefficients of the basis functions, which in this case correspond to the value of the approximate solution at any interior grid point $x_i$. We can plot the approximation and compare it with the exact solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4FOX6xvHvk04SWiBAIIGEGjpCCE2CNBVU2sFDFbEh\ngqJHUfFnr0fPwS7iQVREpQmIiBQVxUJNQHqR0AMBQi8hZbPv748NXpFilpDd2c0+n+vKRXZ3Zvce\nhTuTd2beEWMMSimlfIef1QGUUkq5lxa/Ukr5GC1+pZTyMVr8SinlY7T4lVLKx2jxK6WUj9HiV0op\nH6PFr5RSPkaLXymlfEyA1QEupWLFiiY2NtbqGEop5TVWr159xBgT6cyyHln8sbGxpKSkWB1DKaW8\nhojscXZZHepRSikfo8WvlFI+RotfKaV8jEeO8SulVGFyc3NJS0sjKyvL6ihuFRISQnR0NIGBgUV+\nDy1+pZRXSktLo3Tp0sTGxiIiVsdxC2MMR48eJS0tjbi4uCK/jw71KKW8UlZWFhUqVPCZ0gcQESpU\nqHDVv+Vo8SulvJYvlf55xbHNOtSjVElmDJzYA4c2Q84ZsGWBLdvxlZcNthwIrwQV60LFOhBeGXyw\nTH2NFr9SJcmJvXDg9/yvtZC+Fs4dd3794DJQoZbjB0F0S+zxPfArU9l1eb2cv78/jRs3/vPxnDlz\n2L17Nz179vzLGPzYsWPp0qULIsKgQYP4/PPPAbDZbERFRdGqVSvmzZvnttxa/Ep5u5P7YcOXsH4G\nHN7keM4vACo1gPo9oGozqNIESpUH/yAICHZ8+Qdj/AI4cWgvR/duJvPAFjiynZBTO4k49CMV10/H\nfPsom0o1p2an2ynVuCeUKmfttnqYUqVKsXbt2r88t3v3btq3b3/JIg8LC2Pjxo2cO3eOUqVK8f33\n31OtWjV3xf2TFr9S3ij7NGyeC+unw65fAAPRLeHGVyEmESo1hMAQALJy80g7nsnejEz2HTvH3mNH\n2Hssk335X2dz8gABGgANqBgeREyFUBLDDtEm8ydqHlxIqfmjMAtHI3Wvh6YDIP4mHRIqou7du/Pt\nt9/St29fpk6dyoABA/j111/dmkGLXylvcmIf/DoW1k0H2zkoH4dJeowjNXuxy1Rh77FM9m7OJO23\nrY7vj2Vy+HT2X94iJNCPmPKhVI8IpXXNClSPCCUmwvE4unwpwoIL1kIflqVm8MTnM+hmX8o/9yQT\ntHUeVGsB178ENdq6d/sv4/lvNrH5wKlifc8GVcvw7C0N/3aZc+fO0axZMwDi4uL46quvAPj111//\nfB5g1qxZ1KpVC4D+/fvzwgsvcPPNN7N+/XruvPNOLX6l1MXOZOwl+6exlN86BWNgTUQ3FgZ25qcz\nNUj7MYuc73YBuwDHjnhUmRBiIkJJqhtJ9fxSj4koRUxEKJHhwVd0Zkjb2pFUvX8od33akBeODuSz\nhJ202jUePukG8TdDl+ccB4Z90KWGeoDLDvUANGnShN27dzN16lS6d+/u6oiXpMWvlAeznz7M71Of\npeH+LwnGzrS8Drxn68Xpo1WoHhFKvSqhdG1Qhej8cq8eEUrVciEEB/gXa47YimHMHtGOB6b+Tr9V\nwrA2XzCm3I/4LX0LxrWChDugwxgId2pW4GJX2J65p+nRowejR49myZIlHD161O2fr8WvlCfKzSJ3\nyX+xL3uPZvZsUsrdQGbrR2gcU48FEaGUDS365fpFVbZUIB/fnsBL325hwrLd/FEviXfvHUjpFa9D\nyieOg8s93oGGvd2ezdvceeedlCtXjsaNG7NkyRK3f75TF3CJyI0isk1EUkVkzCVejxeR5SKSLSKj\nL/G6v4j8LiLuO19JKW+1dwW299sRuHQsi2zXMKvNVyQ+NI2ObRJpHF3WktI/L8Dfj+d6NOTl3o34\nbfsRen+6nT2tX4CRKyEyHr4cCgvGOK4P8GHnx/jPf82cOfMvr0dHRzNq1CiL0oEYY/5+ARF/4A+g\nK5AGJAMDjDGbCyxTCagB9AKOG2PGXvAeDwMJQBljzM2FhUpISDB6Ixblc7LPwOLnMas+JJ2KPGu/\nh/4DhtK5vmeeR79sxxHu+3wNfgLjB7egdfXS8P0zsHI8RCfCrZOgrOtOVdyyZQv169d32ft7sktt\nu4isNsYkOLO+M3v8iUCqMWanMSYHmAb0LLiAMeawMSYZyL1wZRGJBm4CJjoTSCmflLoY3m+DWfUh\nn9lv4PaQtxk9YoTHlj5A21oV+XpkOyLCghg8cSXT1hyEbq86Cv/wZvhfe9jxo9Ux1SU4U/zVgH0F\nHqflP+est4DHAPsVrKOUb8g+A3NGwud9OJrjR9/sZ5hX9SGm3d+FelVKW52uUOcP+ratXZExszfw\nwjebyavfC4YtgbBK8FkfWPIa2PWfvydx6SRtInIzcNgYs9qJZYeJSIqIpGRkZLgyllKe4UgqTOyC\nWTeFheUH0vb4c9Rs3pnP725FhfBgq9M57fxB36FtY/l46S7u+jSZU+GxcM9iaNIPlrwCc+4De57V\nUVU+Z4p/PxBT4HF0/nPOaAf0EJHdOIaIOonI55da0BgzwRiTYIxJiIy05pQwpdxm81yYcB15Zw7x\nVPgL3HfwZkZ3b8p/+jYhKMD7Js298KBvn/eXsec00PsD6PgUrJ+m5e9BnPkblgzUEZE4EQkC+gNz\nnXlzY8wTxphoY0xs/no/GmMGFzmtUt4uzwbfPQUzbuNs2Vr0zP03c07WYeKQBO5Jqun10wwPalWD\nyXclknE6m57jlrJi1zHo8Ch0esoxvcRXw7X8PUChxW+MsQH3A4uALcAMY8wmERkuIsMBRKSKiKQB\nDwNPiUiaiJRxZXClvM6Zw/BZL1j2Lrvi+tMq/RFOBFZi9oh2Hn0Q90qdP+hb4fxB31V7IelR6PwM\nbJgBX93r+AGoLOPUBVzGmPnA/Aue+6DA9wdxDAH93XssAZZccUKlSoIDa2Fqf8y5Eyys/Rz3baxL\nYmwE4wc396rxfGcVvNJ3zOwN/HHoDE/e9DD+CCx+Howdek8Af+++hvTCaZn79+/PmDEXXepUJGvX\nruXAgQMumdbBu/+rK+UNUhfDjCHYQ8rxStQ7TNwYxq0tonm5d2OvHM93VsErfT9euosdGWd4d+AD\nlBE/+OFZx01i+nzo1eV/ubl6isPatWtJSUlxSfGX3L91SnmCddNhyj/JLVOdof4v81FqGE92r++1\nB3GvVMGDvktT8w/61r8Hur4Im2bDN6McPwBKkJMnT1KvXj22bdsGwIABA/jwww8BuO+++0hISKBh\nw4Y8++yzf66TnJxM27Ztadq0KYmJiZw8eZJnnnmG6dOn06xZM6ZPn16sGb33R61SnswYWPo2/PAs\np6Pa0CPjPjJyQ/jo9mZ0ii854/nOGtSqBnEVw7jv8zX0HLeUDwYPonWHM/Dza1ClMbS+7+o+YMEY\nOLiheMKeV6Wx44K0v1FwWmaAJ554gn79+vHee+8xdOhQHnzwQY4fP84999wDwMsvv0xERAR5eXl0\n7tyZ9evXEx8fT79+/Zg+fTotW7bk1KlThIaG8sILL5CSksJ7771XvNuFFr9Sxc9uh0VPwMoP2F+t\nGzfsGUj5MuHMvqcldSt7/kVZrnL+oO9dnyYzeOJKXuo5kP7xm2DRk455fmp1tDriFbvcUE/Xrl35\n8ssvGTlyJOvWrfvz+RkzZjBhwgRsNhvp6els3rwZESEqKoqWLVsCUKaM68+L0eJXqjjZsmH2MNg8\nh9VRA+i74yZaxlYssQdxr9RfDvp+tYndrR/i8YqpyJdDYdhPEFGzaG9cyJ65u9ntdrZs2UJoaCjH\njx8nOjqaXbt2MXbsWJKTkylfvjxDhw4lKyvLknwlf5BRKXfJzYJpA2HzHGZWGM4/dt1C3xbVve5K\nXFcreKXvBysO81jAExiAqQMdt5QsAd58803q16/PlClTuOOOO8jNzeXUqVOEhYVRtmxZDh06xIIF\nCwCoV68e6enpJCcnA3D69GlsNhulS5fm9GnX/PfQ4leqOORmwfRBkPoDb4eN4tEDST51EPdKnT/o\n+0rvxny1J4jH/R7GHNmWf4GX98zrc36M//zXmDFj2LZtGxMnTuT111+nffv2JCUl8dJLL9G0aVOu\nueYa4uPjGThwIO3atQMgKCiI6dOn88ADD9C0aVO6du1KVlYWHTt2ZPPmzS45uFvotMxW0GmZlVcp\nUPov+d/HNFtH3hngmwdxi2LZjiOM+GINg808RptP4bon4LrCz4XXaZldOy2zUupyCpT+k3nDWBR8\nPbNHtNXSvwJta1Vkzoh2LAjrzey89rDk37DlG6tjlWha/EoVVW4WZvpgSP2Bx3LvYXu1PswZ0c6n\nz9wpqtiKYcweeS3fxo5hrb0mWTOHYzueZnWsEkuLX6miyM0ib9ogJPV7Hs+9B665TQ/iXqWypQL5\n39C2/NLoVey2XLZMGIo97+/H+z1xqNrVimObtfiVulJ5uWRNvQ3/HT8wJvce6tw4gtf+oQdxi0OA\nvx+j/nkDv9d7kMbnktn87buXXTYkJISjR4/6VPkbYzh69CghISFX9T56cFepK2EMx6bcTcT2mbxo\n7qbdgMd0PN8F8vLyWPdKB+LzUin14AqkfOxFy+Tm5pKWlmbZufBWCQkJITo6msDAwL88fyUHd/UC\nLqWuQOqUh6m9fSYfBfan393P6Xi+i/j7+5PeYSx1F/fkxNR7KT98Afj99TeqwMBA4uLiLEro3fR3\nU6WctGzyM9Te/jELQ2+h16i3tfRd7Pp2ibwXeAflD6/AJH9odZwSRYtfKSesn/c+bXe+zZrSHen0\n0CdUKH11Y6yqcIH+fkR3Hs6SvKbYv3sWju6wOlKJocWvVCEOrPqKBslPsi6wGQ1HTiEoKLDwlVSx\nuLVlDK+HjCTL7gdzRuhtG4uJFr9SfyMz9TcqzB/GHxJL5WEzCQ4JtTqSTwkO8Kd3h0Sezr4N9q2A\nFe9bHalEcKr4ReRGEdkmIqkictG11CISLyLLRSRbREYXeD5GRH4Skc0isklEHizO8Eq5kjmyHTOl\nP/tNBc7cOo0qkZFWR/JJAxKr83NIZ1aXags/vgTH91gdyesVWvwi4g+MA7oBDYABItLggsWOAaOA\nsRc8bwMeMcY0AFoDIy+xrlKeJ/MYpz7qTVYerGo7gcSG9axO5LNKBflzd1It7j8+ALsBvn/G6khe\nz5k9/kQg1Riz0xiTA0wDehZcwBhz2BiTDORe8Hy6MWZN/vengS1AtWJJrpSr2HI4Oak/IZnpTKr+\nMv2ub291Ip93W5saZJaqwrwy/WDzHNi91OpIXs2Z4q8G7CvwOI0ilLeIxALXACuvdF2l3MYYzs5+\ngLKHV/JG6CiGDx6EiFidyueFBwdwZ7s4Hku/jtzwqrDwcT3QexXccnBXRMKBWcBDxphTl1lmmIik\niEhKRkaGO2IpdZHcX94kbPM0PjD/oN+djxAWrNc4eoqhbWMJCA7js/C7HPfX/f1zqyN5LWeKfz8Q\nU+BxdP5zThGRQByl/4UxZvblljPGTDDGJBhjEiL1IJqywua5BP70PN/ktabmrS9TMzLc6kSqgLKh\ngQxpU4MX98RzrkpL+PFFyLrkfqQqhDPFnwzUEZE4EQkC+gNznXlzcfyO/BGwxRjzRtFjKuVi+9dg\nm3kPa+y1SW3zGtc3irI6kbqEu66NIyQggPGl7oazGfDLf62O5JUKLX5jjA24H1iE4+DsDGPMJhEZ\nLiLDAUSkioikAQ8DT4lImoiUAdoBtwGdRGRt/ld3l22NUkVxKp3cL/pxMC+cj6JfZtSNTa1OpC6j\nQngwg1pVZ9y2Mpyp3x9WjNcreotAZ+dUvs2WTe5H3chN38Tw4Fd5e9QgyocFWZ1K/Y1Dp7Jo/5+f\nuL1xKZ7cMQjikmDAVKtjWU5vvaiUk+zzHyUwfTVj8obz2JA+WvpeoHKZEPolxDBpfSanWj4I2+bD\njp+sjuVVtPiV70r5BL81nzLO1oOkXvfQqFpZqxMpJ93boSbGwDtnu0C5GrDwCcizWR3La2jxK9+0\nbxX2bx/l57wmHGzxCH1bRFudSF2B6PKh9Glejc9SDnEy6VnI2ALrp1sdy2to8Svfc/ogtqmDSbNH\nMLHykzx9SxOrE6kiGHFdbXLz7LyfXh+imjrO8NG9fqdo8SvfYsshb9pt5Gae4NGAx/jvbdfpvXK9\nVGzFMHo0rcpnK/dyuvVoOL4LNsywOpZX0L/xyqeYBY/jv38Vj9nu5V+DelOlrN5QxZuN7FibzJw8\nJhysC1Wa6F6/k7T4le/4/Qtk9cd8YLuFpjfcQeuaFaxOpK5Sncql6daoCpOW7eFsm9FwbCds+NLq\nWB5Pi1/5hkObyJv3MMvsDdlU/0HuulZv0l1SjOxYm9PZNj45Uh+qNNa9fido8auSL+sUuVNv41he\nCG+XfZzXbr1GZ9wsQRpVK0vn+Ep8tHQ3WW1Hw7EdsHGm1bE8mha/KtmMIW/uKPxO7OIx8xCv3t6F\n0CCdcbOkGdmpNsczc5l8vCFU1r3+wmjxqxLNrPoQ/81fMTb3nwzqN5C4imFWR1Iu0Lx6ea6tXZEJ\nv+4h59rRcDQVNs6yOpbH0uJXJdf+1dgXPsHivGsITHqILg0qW51IudD9nWpz5Ew2U081gUoNHXv9\nerOWS9LiVyXTueNkT72Ng/ZyzKr+FA92jbc6kXKx1jUrkBgbwQe/7CK3/WNwdDtsvOwtQHyaFr8q\neex2smfei9+ZQzwf8iivDErC308P5vqC+zvVJv1kFrMym0GlBvDza7rXfwla/KrEyVv2LsE7FvFq\n3mAevH0A5UJ1xk1f0b5ORZpGl+X9n3eRd36vf9NXVsfyOFr8qmRJWw2Ln2d+XiINe42mYVWdcdOX\niAj3d6rD3mOZfJ3TAirUhuXjwAPvO2IlLX5VcmSd5OyUIaTby7O22fP0aRFT+DqqxOlSvxLxVUoz\nbslO7InD4cAaSEu2OpZH0eJXJYMxnJr5AMFnD/B+xBOM7tna6kTKIiLCA53qsCPjLN8FdoTgsrDi\nfatjeRSnil9EbhSRbSKSKiJjLvF6vIgsF5FsERl9JesqVRzOrfqUMqlf8z///jx4x2CdcdPH3dio\nCrUiw3jrlwOY5kNg81w4mWZ1LI9R6L8OEfEHxgHdgAbAABFpcMFix4BRwNgirKvUVbEf2orfwsdY\nZm9Iy0EvULmMzrjp6/z9hJEda7P14Gl+i+gDGFj1odWxPIYzu0WJQKoxZqcxJgeYBvQsuIAx5rAx\nJhnIvdJ1lboquVkcmzyY0/Zg9nR4k8RakVYnUh6iR9OqxESUYuzKTEz9W2D1JMg5a3Usj+BM8VcD\n9hV4nJb/nDOcXldEholIioikZGRkOPn2ytftn/EIFc9uZ2bMk/TvlGh1HOVBAvz9GHFdbdalnWRd\n1QGQdQLWTbM6lkfwmIFQY8wEY0yCMSYhMlL32lThjiTPotr2z5kV3IshQ+7WGTfVRf7RPJqosiG8\nsqEsRDWDlR+A3W51LMs5U/z7gYLnxUXnP+eMq1lXqcvKOrqX4Pmj2ERNWtz5ls64qS4pKMCP4R1q\nsWrPcVJrDYEjf8DOH62OZTlnij8ZqCMicSISBPQH5jr5/lezrlKXZOx57PvodvztuZzqPp7YyuWt\njqQ8WL+WMVQMD+bFXfUgvDKsGG91JMsVWvzGGBtwP7AI2ALMMMZsEpHhIjIcQESqiEga8DDwlIik\niUiZy63rqo1RvmHN9Jeok7mGpbVH0yZRz9dXfy8k0J9hSXH8vOMkB+oOhtQfIGOb1bEsJcYDL2VO\nSEgwKSkpVsdQHmj72qXU+OoWNoS25prR8/Dz95jDVMqDnc22ce1rP5JUTXj7wCC4ZhDc/KbVsYqV\niKw2xiQ4s6z+q1FeIyvzNEFzh3FCylL77o+19JXTwoIDuOvaOL7ensPxWr1g7VTIPGZ1LMvovxzl\nNTZ8Mooa9jQOdXqTshWqWB1HeZkhbWMpHRLAuHNdwHYO1nxqdSTLaPErr7B5yQxaZsxmaaUBNE7q\nZXUc5YXKhAQytG0sE7eHkVm1Naz+1GdP7dTiVx7vVEYaVZY8QqpfHM2HvmF1HOXF7mgXR2iQP7NM\nJzi+C/YstTqSJbT4lWczhv2T7iDUnMPWawKlQkOtTqS8WERYELe1rsEru+tiDyoDayZbHckSWvzK\no238+g3qn13FsloPEd9Ep2RQV++u9nHY/UNYEd4JtsyFcyesjuR2WvzKYx3dvZFaa19ldWAL2g98\nwuo4qoSoVDqEAYnVee1QS7BlwYYvrY7kdlr8yiMZWw6nptxBlgkiYuCHBAb4Wx1JlSDDkmqymTjS\nQ2rD759ZHcfttPiVR9o47Wnicv5gbdPniIurZXUcVcJULVeKvi1imHj2WkhfB+nrrY7kVlr8yuOk\nb/qF+tsn8GtoFzr0utvqOKqEuq9Dbb7Ka4dNgnxur1+LX3mUvKzTMPteDkkF6gx9Hz8/nWpZuUb1\nCqFc16wui/ISsK+fAblZVkdyGy1+5VG2fPoglW3p7Gw3liqVKlsdR5VwI66rzRTbdfhlnYCt86yO\n4zZa/Mpj7Fk+i0bps1gc8U+u7aJ36FSuV7tSOOUbdmafqYQteZLVcdxGi195hKwThyjz3b/YTg0S\nhr6ud9NSbjOyU11m2JII2PsrHNtldRy30OJX1jOG3ZPuIdR+luPdxlG+bGmrEykfUj+qDIdq/QM7\nQnaKbxzk1eJXlkv94SPiT/zMj1WHkdiqvdVxlA8a1LUtv+Q1IXf152DPszqOy2nxK0udObSLKkuf\nZr1ffTrc/pzVcZSPahpTjg2VexCefYjsbd9bHcflnCp+EblRRLaJSKqIjLnE6yIi7+S/vl5Emhd4\n7V8isklENorIVBEJKc4NUF7Mbid98p1g7Pj1+R+hIcFWJ1I+rE23wRw1pTnw04dWR3G5QotfRPyB\ncUA3oAEwQEQaXLBYN6BO/tcwYHz+utWAUUCCMaYR4I/jhutKsXXuWOqcXcPSWg/TqFFTq+MoH5dQ\nqworwztT7fASss+U7LtzObPHnwikGmN2GmNygGnAhefa9QQmG4cVQDkRicp/LQAoJSIBQChwoJiy\nKy92fM9GYtf+l1WBLek4YLTVcZQCoFr7IQRhY/XCkn2Q15nirwbsK/A4Lf+5QpcxxuwHxgJ7gXTg\npDHmu6LHVSWBseVw4os7OWeCqTDwfwQF6gRsyjM0SezIsaCqNDxWsmvKpQd3RaQ8jt8G4oCqQJiI\nDL7MssNEJEVEUjIyMlwZS1ls04zniMvZxurGT1NLJ2BTHkT8/IhoNYCy6cvgTMntIWeKfz8QU+Bx\ndP5zzizTBdhljMkwxuQCs4G2l/oQY8wEY0yCMSYhMjLS2fzKyxzcupx62z7g11Id6dRnmNVxlLpY\no75g7LB5jtVJXMaZ4k8G6ohInIgE4Tg4O/eCZeYCQ/LP7mmNY0gnHccQT2sRCRXHpZidgS3FmF95\nkbzsTGxfDuMYZak5RCdgUx6qcgOIrA8bZlqdxGUKLX5jjA24H1iEo7RnGGM2ichwERmev9h8YCeQ\nCnwIjMhfdyUwE1gDbMj/vAnFvRHKO2z8/FGi8/byR+tXqRZV1eo4Sl1e43/AvhVwYl/hy3ohMcZY\nneEiCQkJJiUlxeoYqhjtWb2ImLn9+KXsLXT412Sdi0d5tmM74Z1roOsL0O5Bq9M4RURWG2MSnFlW\nr9xVLpd99gTB397PfqlM46Fva+krzxdRE6q1KLHDPVr8yuW2fjqKSnkZHOz0JhUiIqyOo5RzGv0D\nDq6HI9utTlLstPiVS23/dSZND3/NL5EDaJnU3eo4SjmvYR9AYOMsq5MUOy1+5TJnjh8mYvFodkh1\nEu4Ya3Ucpa5MmSiIvdYx3OOBx0Kvhha/cpkdk+6ltDnFuZvHEx4WZnUcpa5coz5wdLtjyKcE0eJX\nLrFh0cc0Pfkjy2PuplGLa62Oo1TRNOgFfgElbrhHi18Vu2MH91B9+dNs9a9L6yEvWh1HqaILjYBa\nnWDjbLDbrU5TbLT4VbEydjsHJt9NkMkhsO8EgoN0jn3l5Rr1hZP7IG2V1UmKjRa/Kla/z3mLRpmr\nSKn7ELXqX2N1HKWuXnx3CAgpUef0a/GrYnNw9xbi173K+qBmtO1/0Y3alPJOwaWh7g2OSdvybFan\nKRZa/KpY2G02Tky5Cxv+VBg4EX9/nWNflSANe8PZDNi30uokxUKLXxWLlGkvEp+zic3NnqJabB2r\n4yhVvGp3Af9g2Pqt1UmKhRa/ump7Nq+i2fb3WB16La163md1HKWKX3BpqHkdbJ1XIi7m0uJXVyUn\nOwvbrHs5LWHUuH0C4qd/pVQJFX8TnNgDhzZZneSq6b9SdVXWfDaGWnk72dP2FSpWvvBWzEqVIPW6\nAVIihnu0+FWRbUv5gZb7JrGqbDeaX3/JWykrVXKEV4KYVrD1G6uTXDUtflUkmWdOEP7tSA75RRJ/\n5/tWx1HKPeJvgoMb4Pgeq5NcFS1+VSQbPhlFlP0Qx7q+Q5myOse+8hHxNzn+3Dbf2hxXyaniF5Eb\nRWSbiKSKyEVX5uTfZP2d/NfXi0jzAq+VE5GZIrJVRLaISJvi3ADlfsnfTaPV0a9ZGTWQRm27WR1H\nKfepUMtxI3YvH+cvtPhFxB8YB3QDGgADRKTBBYt1A+rkfw0Dxhd47W1goTEmHmiK44btykvt37+P\nuGWPsds/lhZDX7c6jlLuF38T7FkKmcesTlJkzuzxJwKpxpidxpgcYBrQ84JlegKTjcMKoJyIRIlI\nWSAJ+AjAGJNjjDlRjPmVG+Xa8tg7eRhlzFmC//kRQSGlrI6klPvF3wTGDn8stDpJkTlT/NWAfQUe\np+U/58wycUAG8ImI/C4iE0VE78jhpRZNeZs22cvY0fhBouolWB1HKWtUvQZKV/Xq4R5XH9wNAJoD\n440x1wBngUvO3iUiw0QkRURSMjIyXBxLXamVv68jacd/2BXWlPp9nrQ6jlLWEXHs9acuhpxMq9MU\niTPFvx+IKfA4Ov85Z5ZJA9KMMednNpqJ4wfBRYwxE4wxCcaYhMjISGeyKzc5fCoT/7n3ESAQdfsk\n8NMJ2JTyB/cVAAASyklEQVSPi78JbOdg509WJykSZ4o/GagjInEiEgT0B+ZesMxcYEj+2T2tgZPG\nmHRjzEFgn4jUy1+uM7C5uMIr17PbDT989DQJZhOnOr5ESKWaVkdSynqx10JwWa8d7gkobAFjjE1E\n7gcWAf7Ax8aYTSIyPP/1D4D5QHcgFcgE7ijwFg8AX+T/0Nh5wWvKw3357Xz6nviEfVW6EJN0l9Vx\nlPIM/oGOOfq3LXDM0e9faJV6FKfSGmPm4yj3gs99UOB7A4y8zLprAT0S6IVWpx6gecqjnAsoR/SQ\nCY6xTaWUQ/xNsGEG7Fvh+A3Ai+iVu+qSTmTmsHPqaOrIfgL6/g8Jq2B1JKU8S+3OXjtHvxa/uogx\nhsmfTeTWvG/JaHgXYfW7Wh1JKc/jxXP0a/Gri8z45XcGHHiVY2G1iez1itVxlPJc8TfBib1waKPV\nSa6IFr/6i037T1Bh8aOUk0zKDZ4EgSFWR1LKc9W90fHn9u+szXGFtPjVn85m21g4+T908Ush57qn\n8ItqbHUkpTxb6cpQpQls/8HqJFdEi1/96Z0ZC7kvayInqrQjLGmU1XGU8g51usK+lZB10uokTtPi\nVwDMSd7JLdufRAJDKDdwIui9c5VyTu0uYPJg5xKrkzhN/3Urdmac4cS8p2jkt5ugPu9DmapWR1LK\ne0QnOq7iTfWe4R4tfh+Xbctj0qcTGSrfcrbpUPwb3Gx1JKW8i38A1OzgGOf3ktM6tfh93DtfL+OB\n069zpmxdwm5+1eo4SnmnOl3h9AE47B1TkWnx+7DvN6WTuPb/KOeXRfigyRCoN1ZRqkhqd3H86SXD\nPVr8PurAiXOsm/kKHfzXw43/hkr1rY6klPcqUxUqNYTt31udxCla/D7Ilmfn7c9m8KD9C87WvJHA\nRJ11U6mrVqcL7F0B2aetTlIoLX4fNP67ddyb8Qq5pSIJ6zteZ91UqjjU7gr2XNj1i9VJCqXF72OW\nbc+g2rKniPU7RGj/jyA0wupISpUMMa0gKNwrhnu0+H3IkTPZ/DDtTfr4/4at/eNeN4e4Uh4tIMgx\nW2fqYo8/rVOL30fY7YY3v5jDo7YPOVPtWoI6PmZ1JKVKntqd4eReOPKH1Un+lha/j5i0ZCN37H8O\nE1yG8AGf6A3TlXKF2vn3rvDw0zqdKn4RuVFEtolIqoiMucTrIiLv5L++XkSaX/C6v4j8LiLziiu4\nct7avccp/9MT1PQ7SKn+H0N4JasjKVUylYuByHiPH+cvtPhFxB8YB3QDGgADRKTBBYt1A+rkfw0D\nxl/w+oPAlqtOq67YqaxcFnw+lt7+v5LTbjRSs4PVkZQq2Wp3gT1LIees1Ukuy5k9/kQg1Riz0xiT\nA0wDel6wTE9gsnFYAZQTkSgAEYkGbgImFmNu5QRjDO9MnctD2RM4FdWWkM4X/bKmlCputbtAXg7s\n/s3qJJflTPFXA/YVeJyW/5yzy7wFPAbYi5hRFdHM5dvov+spTFBpygycpOP6SrlDjbYQGOrRwz0u\nPbgrIjcDh40xq51YdpiIpIhISkZGhitj+YQ/Dp4ieOEj1PRLJ6T/x447BSmlXC8gGOKSIPV7jz2t\n05ni3w/EFHgcnf+cM8u0A3qIyG4cQ0SdROTzS32IMWaCMSbBGJMQGRnpZHx1Kedy8vhu0ov08PuN\nzLaP4VfrOqsjKeVbaneB47vh2E6rk1ySM8WfDNQRkTgRCQL6A3MvWGYuMCT/7J7WwEljTLox5glj\nTLQxJjZ/vR+NMYOLcwPUxT6dPpV7z33E0WqdCO+i4/pKud352To9dLin0OI3xtiA+4FFOM7MmWGM\n2SQiw0VkeP5i84GdQCrwITDCRXlVIb5buY7eqU9yulQUFQZ/ordQVMoKEXFQPs5j5+0JcGYhY8x8\nHOVe8LkPCnxvgJGFvMcSYMkVJ1RO23v4JBXm30tZOYf/kHlQqpzVkZTyXXFJsGkO5Nkcd+nyILo7\nWELk2Oys/fgBWsgWzt7wJoFVG1sdSSnfVrMDZJ+Eg+usTnIRLf4SYv7Ud+mR9TW7aw+hQptBVsdR\nSsW2d/zpgcM9WvwlwKqVv3J96svsDmtC7IA3rI6jlALH1CiVGsDOn61OchEtfi93+FA6UQvuJssv\nlCp3TQP/QKsjKaXOi+vguCuXLdvqJH+hxe/FcnOyOfxRfyqbI2T2nkRIxIUXVCulLBWXBLZzkJZi\ndZK/0OL3Url5dpa9dzeNctay9poXiG5yndWRlFIXqtEWxA92edZwjxa/F7Ll2Zn1wXN0ODWXdTWG\nktjrb8+kVUpZpVQ5iGrmcQd4tfi9jC3PzvuffETfw++yp2ISTW/Xg7lKebSaHSAt2aOmadbi9yK2\nPDuvfP4tQ/Y9y8nwmtS4Z4rOuKmUp4tLArsN9i63OsmftPi9RJ7d8NS03xi041GCgwKpcPcsCC5t\ndSylVGFiWoN/kEed1ulZ1xGrS8qzGx6bsZoeW58kNiAD/0FzoXys1bGUUs4ICoXoRI8a59c9fg+X\nZzc8+uVaWm58kQ7+6/G/5U2IbWd1LKXUlYhLgvR1cO641UkALX6PZrcbHp+1nuob3qF/wBJIegya\nD7E6llLqSsUlAcZjbseoxe+h7HbDmNnrCVg7mYcCZkOzwdDx/6yOpZQqimotIDDMY4Z7dIzfA9nt\nhv/7agNH1sxlYtAnjps63PIWiFgdTSlVFAFBUKONxxS/7vF7GLvd8OScDWxN+Yn/hbyLRDWGWz/V\nOXiU8nZxSZCxFU4fsjqJFr8nMcbw9NcbWZ68iilhbxBQNgoZ9CUEh1sdTSl1teI6OP70gL1+LX4P\nYYzhma838d3K9XxV5g1KBfkjg2c7pnZVSnm/Ko0hpJxHzNvjVPGLyI0isk1EUkXkort3599k/Z38\n19eLSPP852NE5CcR2Swim0TkweLegJLAGMNzczcxb8UG5pcfSzn7cWTgDKhQy+poSqni4ucPsdd6\nxx6/iPgD44BuQANggIg0uGCxbkCd/K9hwPj8523AI8aYBkBrYOQl1vVpxhie/2YzXy3fxILyr1Mx\n9wAycDpEJ1gdTSlV3OI6wIk9cHy3pTGc2eNPBFKNMTuNMTnANKDnBcv0BCYbhxVAORGJMsakG2PW\nABhjTgNbAJ00Pp8xhhfnbWHmss0siHiLyjl7kH5f5J/zq5QqcWp6xji/M8VfDdhX4HEaF5d3ocuI\nSCxwDbDyUh8iIsNEJEVEUjIyMpyI5d2MMbz87RamLd3C/ArvUDXrD+TWT6FOF6ujKaVcpWJdCK/s\nFcV/1UQkHJgFPGSMOXWpZYwxE4wxCcaYhMjISHfEsowxhn8v2Mpnv21jXsX3iMnciPxjIsR3tzqa\nUsqVRKBGO9izDIyxLIYzxb8fiCnwODr/OaeWEZFAHKX/hTFmdtGjlgzGGF5duJVJv2zjm8gPiDvz\nO9LrA2jY2+poSil3qNEWTu2HE3sti+BM8ScDdUQkTkSCgP7A3AuWmQsMyT+7pzVw0hiTLiICfARs\nMcb4/B1DjDH8Z9E2Jv+8mXmR46h7egVyy9vQtJ/V0ZRS7lK9jePPvSssi1Bo8RtjbMD9wCIcB2dn\nGGM2ichwERmev9h8YCeQCnwIjMh/vh1wG9BJRNbmf/nkeIYxhte/+4MvlqxnQcSb1DmTAre8Ay1u\ntzqaUsqdKjWA4LKwd5llEZyaq8cYMx9HuRd87oMC3xvgohu/GmN+A3SCGeDNH7Yz7afVLCz3OlFZ\ne5C+n0DDXlbHUkq5m58fVG8Fe6y7I5deuesGb/3wB7N/XMbCMi8RlXcAGTRDS18pX1a9DRzZBmeP\nWvLxWvwu9vYP25m3+Cfmhb1EBb8zyJCvoVYnq2MppaxUo63jT4vuw6vF70LvLt7O4sUL+Dr0ZcqG\n+CF3LICYRKtjKaWsVvUa8A+2rPh1Pn4XGfdTKn8snsSskAkElI5ChsyBiJpWx1JKeYKAYMe0LHus\nOcCre/wuMP6n7bD4ed4Neo+AmBbI3Yu19JVSf1W9jeM+vDln3f7RWvzF7OPF66jz4z2MDJiLvfnt\nyJC5EF6yr0RWShVB9TZg8iAt2e0frcVfjKYvXEL7nwfQ0X8ded3+i98tbztuuaaUUheKSQTxs+S0\nTh3jLyYL5kzhht9HExDgD4Pm4F+rg9WRlFKeLKQMVG5kyYVcusd/tXKz2PjJKG74fQRngysRPOIX\nLX2llHNqtIW0FMjLdevHavFfjYMb4MOONNrzKcvL30Klf/1KYMU4q1MppbxF9TaQm+k4yOtGOtRT\nFPY8WPYu/PgShEbAwC9pV/d6q1MppbzNnxO2LXfrXfd0j/9KHd8Nk26GH56Fet3gvuWgpa+UKorS\nlR2nerv5AK/u8Tsr9xwsfw9+e8txJL73/6BJP8eNFZRSqqiqt4Vt88Fud0zg5ga6x18Yux3WTYN3\nWziGdmpeB/cthab9tfSVUlevRhs4dwyO/OG2j9Q9/r+z+zdY9CSkr4WoZtBnAsRea3UqpVRJ8uc4\n/zKoFO+Wj9Tiv5AxsH8N/PYGbJ0HZapB7wnQ+Fa3/RqmlPIhETUdN2DfsxwS7nTLR2rxn5eTCRtn\nQvJEx6lVQaWh09PQZiQElrI6nVKqpBKB6q3deitGp3ZhReRGEdkmIqkiMuYSr4uIvJP/+noRae7s\nupY7sh0WjIE34mHuA2DLge5j4eHNkDRaS18p5XrV28LJvXAyzS0fV+gev4j4A+OArkAakCwic40x\nmwss1g2ok//VChgPtHJyXffKOumYCnXXr7D7F8dFWH6B0KAHtLzbMd6mB22VUu5UI3+cf89yaHKr\nyz/OmaGeRCDVGLMTQESmAT2BguXdE5icf+/dFSJSTkSigFgn1i1+djucOw5nM+DsYThz2FHwu35x\nHKg1dsdNEGISoctz0GwQhFdyaSSllLqsyo0guIzjAK+HFH81YF+Bx2k49uoLW6aak+sWD2NgQgc4\nfRDOHnFMd1qQX6Djyrj2oyEuCaJbQmCIS6IopdQV8fN37Ii66UIujzm4KyLDgGEA1atXL8obQGQ8\nRDWFsEgIqwRhFR178mGRUK46BIUVc2qllCom9XvAgTVuuZDLmeLfD8QUeByd/5wzywQ6sS4AxpgJ\nwASAhIQE40Sui/WZUKTVlFLKci1ud3y5gTM/VpKBOiISJyJBQH9g7gXLzAWG5J/d0xo4aYxJd3Jd\npZRSblToHr8xxiYi9wOLAH/gY2PMJhEZnv/6B8B8oDuQCmQCd/zdui7ZEqWUUk4Rx4k4niUhIcGk\npKRYHUMppbyGiKw2xjg1t7POQaCUUj5Gi18ppXyMFr9SSvkYLX6llPIxWvxKKeVjPPKsHhHJAPYU\ncfWKwJFijOMNdJtLPl/bXtBtvlI1jDGRzizokcV/NUQkxdlTmkoK3eaSz9e2F3SbXUmHepRSysdo\n8SullI8picXvizO16TaXfL62vaDb7DIlboxfKaXU3yuJe/xKKaX+hlcW/9Xc/N1bObHNg/K3dYOI\nLBORplbkLE6FbXOB5VqKiE1E+roznys4s80icp2IrBWRTSLys7szFjcn/m6XFZFvRGRd/jbfYUXO\n4iIiH4vIYRHZeJnXXd9fxhiv+sIxvfMOoCYQBKwDGlywTHdgASBAa2Cl1bndsM1tgfL533fzhW0u\nsNyPOKYG72t1bjf8fy6H457V1fMfV7I6txu2+f+A1/K/jwSOAUFWZ7+KbU4CmgMbL/O6y/vLG/f4\n/7z5uzEmBzh/A/eC/rz5uzFmBXD+5u/eqtBtNsYsM8Ycz3+4AsfdzryZM/+fAR4AZgGH3RnORZzZ\n5oHAbGPMXgBjjLdvtzPbbIDSIiJAOI7it7k3ZvExxvyCYxsux+X95Y3Ff7kbu1/pMt7kSrfnLhx7\nDN6s0G0WkWpAb2C8G3O5kjP/n+sC5UVkiYisFpEhbkvnGs5s83tAfeAAsAF40Bhjd088S7i8vzzm\nZuuqeIhIRxzFf63VWdzgLeBxY4zdsTPoEwKAFkBnoBSwXERWGGP+sDaWS90ArAU6AbWA70XkV2PM\nKWtjeS9vLP6rufm7t3Jqe0SkCTAR6GaMOeqmbK7izDYnANPyS78i0F1EbMaYOe6JWOyc2eY04Kgx\n5ixwVkR+AZoC3lr8zmzzHcCrxjEAnioiu4B4YJV7Irqdy/vLG4d6rubm796q0G0WkerAbOC2ErL3\nV+g2G2PijDGxxphYYCYwwotLH5z7u/01cK2IBIhIKNAK2OLmnMXJmW3ei+M3HESkMlAP2OnWlO7l\n8v7yuj1+cxU3f/dWTm7zM0AF4P38PWCb8eIJrpzc5hLFmW02xmwRkYXAesAOTDTGXPK0QG/g5P/n\nF4FJIrIBx5kujxtjvHbWThGZClwHVBSRNOBZIBDc11965a5SSvkYbxzqUUopdRW0+JVSysdo8Sul\nlI/R4ldKKR+jxa+UUj5Gi18ppXyMFr9SSvkYLX6llPIx/w/uxzTX6cIagwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d15fdd8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(centers,u)\n",
    "# As a comparison we plot the exact solution in the same graph\n",
    "x = np.linspace(0,1)\n",
    "u_exact = x**2 - x**3\n",
    "plt.plot(x,u_exact)\n",
    "l = plt.legend(['FEM','Exact'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that even for low values of $n$, like $n=10$, the solution is reproduced exactly in each grid node, showing that this approximation by linear basis functions is, in some sense, the best approximation possible. The reason for this is that the load vector $\\vec{b}$ is perfectly approximated; would $f$ have a higher degree, this exact approximation would not be guaranteed. \n",
    "\n",
    "In the general case, two ways to increase FEM accuracy are by\n",
    "\n",
    "* Decreasing the element size\n",
    "* Increasing the degree of the basis functions\n",
    "\n",
    "You can try the first one yourself by adjusting the code. Just increase `n` to 100 or 1000 and you start seeing improvements. Increasing the polynomial order requires more tweaks: the stiffness matrix will get different coefficients. Normally, \n",
    "\n",
    "The FEM-implementations mentioned before have many high-degree basis functions readily available.\n",
    "\n",
    "In the next notebook, we explore how to express and quantify the error between the approximation and the exact solution as a function of $n$ (or $h$, to be precise)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}