{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Accuracy of Newton-Cotes"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import numpy.linalg as la\n",
        "import matplotlib.pyplot as pt"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "A function to make Vandermonde matrices:\n",
        "    \n",
        "(Note that the ordering of this matrix matches the convention in our class but *disagrees* with `np.vander`.)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 34,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "def vander(nodes, ncolumns=None):\n",
        "    if ncolumns is None:\n",
        "        ncolumns = len(nodes)\n",
        "    result = np.empty((len(nodes), ncolumns))\n",
        "    for i in range(ncolumns):\n",
        "        result[:, i] = nodes**i\n",
        "    return result"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Fix a set of nodes:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 129,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "# nodes = [0.5] # Midpoint\n",
        "# nodes = [0]\n",
        "#nodes = [0, 1] # Trapezoidal\n",
        "nodes = [0, 0.5, 1] # Simpson's\n",
        "#nodes = [0, 1/3, 1]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Find the weights for the Newton-Cotes rule for the given nodes on $[0,1]$:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 130,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[ 0.16666667  0.66666667  0.16666667]\n"
          ]
        }
      ],
      "source": [
        "(a, b) = (0, 1)\n",
        "nodes = np.array(nodes)\n",
        "n = len(nodes)\n",
        "\n",
        "degs = np.arange(n)\n",
        "rhs = 1/(degs+1)*(b**(degs+1 - a**(degs+1)))\n",
        "weights = la.solve(vander(nodes).T, rhs)\n",
        "print(weights)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here is a function and its definite integral from $0$ to $x$:\n",
        "\n",
        "$$\\text{int_f}(x)=\\int_0^x f(\\xi)d\\xi$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 131,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fdeg = 9\n",
        "def f(x):\n",
        "    return sum((-x)**i for i in range(fdeg + 1))\n",
        "def int_f(x):\n",
        "    return sum(\n",
        "        (-1)**i*1/(i+1)*(\n",
        "            (x)**(i+1)-0**(i+1)\n",
        "            )\n",
        "        for i in range(fdeg + 1))\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plotted:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 132,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<matplotlib.legend.Legend at 0x7fe907aaa978>"
            ]
          },
          "execution_count": 132,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd0VHX6x/H3d9ILhCQklEAIvYQogqIoSBMIHWmKIrgI\ndrB394egqCgKrLiKDXZXVwQLrF0BQ1FBakIJSm+BdEhPJjPf3x8JDGKAQGbmTnle5+QcLrm58/Cc\n8OTmc+/9jtJaI4QQwjOZjC5ACCGE48iQF0IIDyZDXgghPJgMeSGE8GAy5IUQwoPJkBdCCA92wSGv\nlHpfKZWulEo5zz7/UErtVkptVUp1sG+JQgghLlV1zuQXAP3O9UmlVH+guda6JXAX8LadahNCCFFD\nFxzyWuu1QO55dhkK/Lty3/VAmFKqnn3KE0IIURP2yORjgMNnbB+t/DshhBAGs8eQV1X8nayVIIQQ\nLsDXDsc4AjQ+Y7sRkFbVjkopGf5CCHEJtNZVnVBfUHXP5BVVn7ED/A8YB6CUugY4obVOP9eBXnjh\nKM2a5dKhQwHbt5vRWnvlx9SpUw2vwVU+pBfSC+nF+T9qojq3UP4X+AVopZQ6pJT6m1LqLqXUnQBa\n62+A/UqpPcB84N7zHa9t27q88koILVvm0KWLlb//vQCzuUb/Brd04MABo0twGdILG+mFjfTCPi4Y\n12itb6nGPvdfzIv6+/tx222NueaaLN59V7FkiYn//MePq67yu5jDCCGEuABDn3ht2bIuL74YxtVX\n53DDDZopU/IpLjayIue5/fbbjS7BZUgvbKQXNtIL+1A1zXsu6sWU0p99Voqfn/9fPnf4cC7vvWcm\nPb02b7+tGDw4wGl1CSGEK1NKoR184dXhGjcOZ+rUuowYkc2ECRYGDSogLc1zb8ZJSkoyugSXIb2w\nkV7YnOpFXFwcSimv+IiLi7N7H11myAOYTCb69Ilh9myN2XyCtm3NzJpVhMVidGVCCKMcPHjQ8Ltb\nnPVx8OBBu/fPZeKaqiQnH2fhwgACAwP44AM/OneWC7NCeJvKqMLoMpziXP9Wj4hrqnL55fWZOTOU\nK67Ipk8fK3feWUB+vtFVCSGE+3DpIQ8Vt1uOHt2YV14pJiWlgKZNy3j33WLc/Qe7ZK820gsb6YWN\n9MI+XH7In9KwYR2efroekyZl8cILZXTsWMTGjeVGlyWE8HJ//PEHHTt2JCwsjHnz5hldzl+4dCZ/\nLiUlpSxdms7XX9dj6FAzr78eSmSknYoUQrgUV8/kJ06cSFhYGK+99lqNj+V1mfy5BAYGcPPNsbz6\najEHDpykRYsyZs8ulrtwhBBOd/DgQeLj440u45zccsifUr9+HR55pCEPPpjFP/9ZQnx8MUlJ7rEQ\njuSNNtILG+mFjTv0onfv3vz000/cd9991K5dmz179hhd0l+49ZCHil9jOnZsyMyZoXTtmsXw4eUM\nHFjA3r1Wo0sTQni4FStW0K1bN958803y8vJo0aKF0SX9hdsP+VP8/PwYOrQxc+aUo1QuHTpYuP/+\nQk6eNLqyqvXo0cPoElyG9MJGemFzMb1Qyj4fnshjhvwp4eG1uOuuxsyYcZJNm/Jo2rQiry+XG3GE\n8Fha2+fDE3nckD+ladO6PPlkPe6/P5P584tp1aqEL78sM7qs09whb3QW6YWN9MJGemEfHjvkoWIt\nnKuuiuHFF2uRmJjBpElmrr++kC1b5DYcIYR38Oghf4qfnx/9+8cye7aiXr1sevYsZ8SIQvbvN+73\nM8lebaQXNtILG3fphXLxMN8tH4aqqYyME3z2WT5r1tRnzJgyZswIISrK0JKEEOfg6g9D2ZM8DGUn\n0dF1uOeexrz44gn++COXFi3MPPVUMQUFzqtB8kYb6YWN9MJGemEfXjnkT4mLi+Lhh2N4+ulMVqw4\nQZMmZcyaVUKZ61yfFUKIGvHKuKYqVquVzZvT+PTTQE6erMULL8C4cQGYvPrHoBDGk7imZnGNDPmz\nWCzlrFmTxuef18bXN4AXXvBh9Gh/j31QQghXJ0NeMnm78vHxpUePWGbNCqZv3wwee6yEdu2K+fxz\ns10flpC80UZ6YSO9sJFe2IcM+XPw9/cnMbEJr7/uT/fuGTzwQAkJCcV8+WW5xz4ZJ4TwPBLXVFNx\ncRE//pjBN9/UJSrKh5de8iMx0VdiHCEcTOIayeSdqqiokB9+yOSbb6KIiVG89JI/ffr4Gl2WEB5L\nhrxk8k4VHBzCsGFxzJ4N7dtnMHZsGVddVcQ331xcjCN5o430wkZ6YSO9sA8Z8pcoJCSEESPimDtX\nEx+fycSJJVx2WTGffmrGKkvZCyFchMQ1dlJcXMTy5Rl8/304AQEB/P3vPtxyix8+PkZXJoR7c/e4\nJi8vjx49ehAaGsqKFSvw8/M7574S17iwoKBgBg+OY9Ysf/r2Pc60aQU0a1bK/PllmN3jHQmFEA4w\nc+ZMevfuzYkTJ9i9e7fTX1+GvJ0FBgbRr18cM2cGM2LEMebMySM2tozZs0soKbHtJ3mjjfTCRnph\n4wm90FqzcOFCxo8fz+rVq2nXrp3Ta5Ah7yD+/gH07BnHiy+GMX78cT74II+YGDPPPltMTo7R1Qkh\nnGH9+vX4+PjQvn176tSpY0gNksk7icViITk5je++82HHjmjGjCnj6aeDadLE6MqEcG3unMm/9NJL\nJCcns2jRomrtL5m8G/Px8aFjx8Y8+WR9pk7NYN++LBISzNx4YxFbtrjnN7AQ4vxWrVrFVVddZWgN\nMuSdzGQy0apVQ3r02Msrr+QCGdxwQylduxbx7bcWr1wywROyV3uRXti4ey+sViu//PILV1xxhaF1\nyJA3iFKKmJhoJkyIY86cYpo2TWfSpCLati1hwYIyWdNeCDe3ZcsWCgsLufzyyw2to1qZvFIqEZhD\nxQ+F97XWM8/6fGPgX0Cdyn2e0lp/W8VxvDaTr46SkiJWrcpgxYpQsrLCuPtuC1OmBMpbEwqv5q6Z\n/JtvvsmMGTNIS0ur9tcYkskrpUzAPKAfEA+MUUq1OWu3Z4FPtNYdgTHAPy+lGG8XGBhMv34Vd+Tc\ne+8xVq7Molmzcm65pYjkZPf7JhfCm61fv574+Hijy6hWXNMZ2K21Pqi1NgOLgKFn7WMFalf+uQ5w\n1H4leqZt25LO+TlfXz+uvDKWxx9vyIwZmRQUZNCzZxnXXlvM0qUWj1s2wd2zV3uSXti4ey/WrVtH\n+/btT2/PmTOHZ555hnfeecepdVRnyMcAh8/YPlL5d2eaBtymlDoMfAVMtk953s1kMtG0aQMmTYpj\nzpxCWrc+xsMP59OkSSmvvlpKfr7RFQohqpKbm8uePXtISEgAKpY2+OSTTxg6dChdu3Z1ai3VWSO3\nqhzo7OxgDLBAaz1bKXUN8CEV0c5fvPHGHdSv3xyAkJA6NGvWgYSEHoDt7NYbthMSelzU/uHhEbRs\nmUJc3B6KilqxaFEAU6fupHt3M6+80o+EBHX6zKdHj4qvl2333D7FVeoxavvMv3M3mzZtQil1esiv\nX7+eK664gs6dO1fr65OSkli4cCEAcXFxNarlghdeK4f2c1rrxMrtJwF95sVXpdR2oJ/W+mjl9l7g\naq111lnHkguvdmKxWNi7N52VK838+mt9WrSw8sADfowa5Yu/tFd4kOpceFXT7PPuPXqqfa59zZw5\nk2effZa8vDxSUlKYPn06jRo1on///gwbNuycX2fIm4YopXyA34HewDHgN2CM1jr1jH2+BhZrrf+l\nlGoL/Ki1blTFsWTIV9q2Len0GXtN5ebmsGbNCdauDSMzszYTJliZPDmAxo3tcniHS0pK+tPZmzeT\nXtic6oU73l1z8803k5KSws6dOwEYOXIk06dPv+DaNY4Y8heMa7TWFqXU/cAP2G6hTFVKTQM2aK2/\nAh4F3lVKPUTFRdjxl1KMuDTh4REMGRJB//6lpKSkkZSkaNeuAdddV85DDwXQp48JkzwRIYTTJCcn\nc+WVV57eTk1NpW3btobUImvXeCCtNcePZ5CUVMjatVEo5c+998KkSQFERBhdnRAXx93O5EtKSggN\nDWXWrFk8+OCDZGRkMHr06GpdX5C1a0S1KKVo0KAeY8Y0Y9YsEzfdlMZnn2UTG1vO8OFFrFhh9crl\nE4RwhtTUVKxW6+k1a9avX8+1115rWD0y5A1yvvvk7SkoKIRu3ZryzDP1efnlTHx80hk/voi4uDKm\nTy/j+HGnlHFe7noHhSNIL2zctRc7duzA19cXk8nEPffcw7p167jpppsMq0eGvJcwmUw0adKAceOa\n8tprcNNNR/j++wxatiwnMbGYr76yYrEYXaUQ7i81NZXLLruM4OBgGjVqREREhKHr10gm78WsVitp\naemsXVvMunURFBaGMH68lXvvDZB17oXLcLdMfsiQITRv3pzZs2df9NdKJi/symQy0ahRA26+uRkz\nZ/ozadIR1q9PJyHBzPXXF/PRR+V/estCIcSFpaSkcN111xldxmky5A3irEy+ugICguncuSkPP9yY\n2bNzad78OC++eIJ69coZN66YX37RDrtY667ZqyNIL2zcsRe5ubkcPnzYpYZ8dZY1EF5EKUV0dDTD\nh8PgwWXs2ZPGzz9bGDUqGj8/P8aPhzvu8Cc21uhKhXA9GzdupHXr1jRo0MDoUk6TTF5US2FhHlu2\nZPHrr75s3dqAyy6zcMcdvowe7UtIiNHVCU/mTpn8M888Q35+Pv/4xz8u6esNWdbAnmTIuz+r1Upm\nZgbr1hWyfn0tDhyIYNAgM3feGUj37kqerBV2505Dvm3btrz11luXvDSFXHj1IK6WyVeXyWSiXr36\nDB3anOnTw5k69RilpccZP76AmJgyHnighJQULiq/d8fs1VGkFzbu0oukpCRuueUWPvjgA0JDQ11u\n7SHJ5MUl8/X1o3XrxrRuDWPGFLJt21HWr1f07l2f8HDFrbcqxo/3p4YrpQrh0lq2bMmGDRtISUnh\n888/N7qcv5C4RthdXl4umzfnsGGDH8nJDWjVysJtt/lyyy2+1K1rdHXC3bhTXFNTkskLt2K1WsnJ\nyWLDhjw2bgxm585orr7azO23+zN8uI9csBXVIkNeMnm35K6Z/MUwmUzUrRtN//4tePrpesyalUFs\nbBqvvJJNvXrlDBlSzJIlFr77LsnoUl2Gu+TQziC9sA/J5IVT+Pj4EBPTkOHDYcgQM4cOHee338p4\n7rna7NtnITGxhLFj/RkwwERQkNHVCuE5JK4RhjKbSzl4MJ0NG8rYujWcw4fr0LevmbFj/enf30Rg\noNEVCqNJXCOZvPAQZWUllQPfzNat4Rw9Gka/fuWMHetPYqKJgACjKxRGkCEvmbxb8oZMvrpO9cLf\nP5CWLZtwyy0teP75YKZOPYqPz1GeeOIE0dHljB5dzNKlVoqLja3XkSSHtjnViyZNmqCU8oqPJg5Y\n/lUyeeGSAgKCaNWqCa1awciRxRw4cJTffrPw1FNhjB1bh549zdx8sz+DB5uoXdvoaoUjHThwwOgS\n3JrENcKtlJWVcOhQBps3l5GSUps9eyK5+mozN93kx403+hAdbXSFQtifZPLCK5nNZRw7ls7mzSWk\npISSmhpF+/bljBrly6hRvrJSpvAYksm7IcnkbS61F35+/sTGNmbYsJY8/XRdZs9Op2PHw3z++XES\nEswkJJQybZqZ1FT71utIksnbSC/sQzJ54RF8ff2oXz+G/v2hb18LubmZpKQU8NNPfrzxRn3CwmDo\nUBg1KoDOncHHx+iKhXAOiWuER9Nac/JkDjt2nCQ5WbFjRxSFhYEkJpYzalQAffooWV5BuDzJ5IWo\npqKiPPbsyWbrVgs7d9bhwIFwrr3WzPDhfgwd6kPDhkZXKMRfSSbvhiSTt3FmL4KDa3PZZU0ZN64F\n06eH8tJLR4mNPczChem0aVPO5ZeXMnWq+aLXxLcXyaFtpBf2IZm88Fr+/oE0bRpL06YweHA52dmZ\nbNtWwJo1/syfH42Pjy+DBlkZOdKf7t0V/vILqHBDEtcIcRatNfn5ufz+e25ljh9BRkYoPXuaGTEi\ngEGDTEREGF2l8CaSyQvhQCUlhRw8mMXWreVs3x7K7t11ad/ezLBhPgwf7kfr1kZXKDydDHk3tG1b\nEgkJPYwuwyW4Uy/Ky81kZmaSnFzMtm3+7NgRTa1aioEDrYwcGUDXrgo/v0s/flJSksu9R6hRPLEX\nFquFYwXHSMtPI6Mwg8zCTDKLMskszCS3JJcic9Hpj0JzIaXlpVi1lS13b7nkIS+ZvBAXwdfXjwYN\nGtKgAfTrVxHr7NhxkpQUmDAhgpycEHr1KmP4cH8GDfIlPNzoioWzlZaXsjtnN6mZqezK2sWe3D0c\nPHGQQycPcTT/KJFBkcTUjiE6JJqo4KiKj5AoWkW2IsQ/hGC/4NMfAT4B+Jh86HR3p0uuR87khbCT\n0tIi9u/PZutWM9u3h7JnTwQJCWaGDTMxYkQALVsaXaGwJ6u2sjt7N5uObSL5eDKpWamkZqVy+ORh\nmoY3pU3dNrSt25aWES1pUqcJsWGxNK7dmADfi18zW+IaIVyMxWIhMzOTzZuL2bbNjx07oggPh4ED\nLYwcGch115nwld+j3YZVW9mVtYuNaRvZfGwzm45tYuvxrUQFR9GpYSc61OtAu6h2tKnbhuYRzfH3\nse+MkyHvhtwph3Y0T++F1pqCgjy2bz9BcrJi+/Y6nDwZTO/epYwYEcDAgb6EhVXs64k59KUysheF\nZYX8dvQ3fj78M78c/oVfj/xKZFAkV8VcRcf6HenUsBNX1L+C8CDn5HE1GfJyLiGEgymlqFUrjC5d\nwujSpeItD/ftS2fTJjMvvhjMxIkRdOhQxtChJho3Nrpa75RXmsfqg6tZuX8law6tITUzlcvqXcZ1\nja9jUsdJLBi6gHqh9Ywu85JU60xeKZUIzKHiCdn3tdYzq9hnNDAVsALJWuuxVewjZ/JCnMFqtZKR\nkc3GjUVs2+bLzp11iYzUDBxoZfToQLp0Mcliag5QbC7ml8O/sHL/SlbsX8GOzB1cHXM1vZr2onuT\n7nRq2IlAX9d5g2GHxjVKKRPwB9AbSAM2ADdrrXedsU8L4BOgp9Y6TylVV2udVcWxZMgLcR4FBXmk\npJxg61bYvj2cgoJA+vQpY8SIAAYM8KVWLaMrdE9WbWXzsc18t+c7Vu5fyW9Hf+OyepfRu2lvejXt\nRZfGXVxqqJ/N0UP+GmCq1rp/5faTgD7zbF4pNRP4XWv9wQWOJUO+kqfn0BdDemFzZi/M5lL27Mli\n0yYz27eHcOBAOJ06lTFsmC8jRvjjgLcDdSk1zeRzi3P5Ye8PfLvnW77d8y0RQREkNk/khmY3cH2T\n66kV4D4/MR2dyccAh8/YPgJ0PmufVpWFrKUi0pmmtf7+UgoSQlTw8wugbdsY2ratuFsnPT2LTZuK\nWLLEn2nToqhfXzN4MIwZE0DHjqAuaQR4Dq01yenJfLP7G77Z/Q0p6Slc3+R6BrQcwNTuU2ka3tTo\nEg1RnTP5kUBfrfWdldtjgau01g+csc+XQBkwCogF1gDxWuu8s44lZ/JC1JDWmry8XJKTT5CSYiIl\nJRqTyZfhwy3cemsgV1+tvGbgmy1mVh1cxdJdS1n2+zICfQMZ0GIAA1oOoHtcd5eOYC6Go8/kj1Ax\nuE9pREU2f/Y+v2qtrcABpdTvQEtg09kHe+ONO6hfvzkAISF1aNasw+lfT08tOSvbsi3b598OC4sg\nPDyF7t3hb3+LYMeODH74YR0ffxyGr29vhg2zEB+/jvh4Ra9eFV9/auneUxGIu2536tKJ7/d+z/xP\n57PuyDriO8czrM0wZjSbQWxYrOH12WM7KSmJhQsXAhAXF0dNVOdM3gf4nYoLr8eA34AxWuvUM/bp\nV/l3tyul6lIx3DtorXPPOpacyVeSHNpGemFT016UlBSyY0cWGzdqkpMjKS0NZOjQcsaODaRrV+VW\nd+qcmcmnF6Tz5R9fsnTXUlYfXM21ja9lWJthDGk9hIa1PP+dXhx6Jq+1tiil7gd+wHYLZapSahqw\nQWv9ldb6e6VUX6XUDqAcePTsAS+EcLzAwBA6dQqhU6eKgZ+amsbGjVbGjYukuDiI0aPLmTAhkCuu\ncP1IJ6soizfWv8HinYvZlr6NxBaJjL1sLB8N/4iwwDCjy3Mb8sSrEF6gpKSQnTszWbdOsWlTxcqZ\nt94KEyYEUsM0wK6OFxzns52fsXjnYlLSUxjcajCj40fTp1mfS1rzxVPIsgZCiGorLMxn06YcfvvN\nh61bo2nRwsK4cT6MHetPZKTz68kozDg92Lcc28KgVoMYHT+avs37esyF05qSIe+GJIe2kV7YOLsX\nubk5rFt3kt9+C2TXrih69izj3nsD6NfPB5MD3wE6szCTL3Z9weIdi9mYtpEBLQcwOn40/Zr3I8gv\nCJB1fM4ka9cIIS5JeHgE/ftH0K+flaNH01m9upR77gmntDSIceMs3HtvkN0euiooK2DprqV8mPIh\n646sI7FFIvdedS/9W/Q/PdiF/cmZvBDiT8rKSti6NYPVq01s3lyPK64wc//9/gwf7nvRyyObLWZ+\n3PcjH237iK//+JqusV25NeFWhrQeQoh/iGP+AR5I4hohhENkZWWzenUea9bUobAwhDvvtDB5chB1\n6577a7TW/Hb0Nz5M+ZDFOxfTLLwZYxPGMjp+NFEhUc4r3oPIkHdDkkPbSC9sXLUXZnMZmzens2KF\nD9u2RTFoUBmPPx5Ex4624P6P7D/4KOUjPtr2ET4mH25NuJVbE26leUTzS3pNyeRtJJMXQjiUn58/\nV1/dmM6dNYcOZfDjj6X06RNNo9bH6DR+KdtNH3Po5CFubn8zH4/4mCsbXoly9RvxvYScyQshLorZ\nWsaGk1/xQ8Z7bM9fi9++/oTuH8v/je3FxAkh+Mt/b7uTuEYI4VBaa/YVb2VF9gJW535MbGA8vSNv\np0udEfha/Fm7Np3vvgsmJ6c2Dz1kZfLkQELkuqrd1GTIO/BOWHE+pxabEtKLM7laL06YM1iWPpsH\ndnXgpX03EuoTwazW63mxVRK9I28n2KcW/v4B9OoVy0sv1eHOOzNYvPgksbFlvPRSMcXFl/7apxbs\nEjUjmbwQ4k/M1jI25n3NiuyF7ChYReewoUxqNJf40OsxqXOfF/r4+NK5cyOuvNLKli3HWLTIj9de\n8+GJJ8qZMiWYAO9dlcBQEtcIIQDYV7SVFdkLWZ37XxoFtqV35O1cW2ckwT6X9g5KVquVjRuP88UX\nfmRm1mbaNCt33hnk0CdpPZVk8kKIS1JoOcmqnP/yQ9a75Fuy6R1xOz0jx9Eg4NJue6yKxWJh/fpj\nLF4cgp9fAHPn+tKvn8yAiyGZvBtytezVSNILG2f0QmvNzoK1zD4wnonbm7A9P4nxMTN5N34/tzSc\nZtcBD+Dj48O11zbi1VdD6dYth7Fjy+nVK5/UVMt5v04yefuQTF4IL3HCnMFPOf/mh6z3MCkTfSIn\nMiFmFmF+znkK1c/PjyFDGtGzZzFLlpykc+cAbr21iFmzahEa6pQSvJLENUJ4MIu2kJy/nB+y3iU5\nfznXhN1I37qTaBPSxfCHlY4cyWXhwjIOHKjDP/5h4aabgg2tx5VJJi+E+JPMssMsz/6A5dkfEOYb\nRd/ISXSLuJkQH9d6RyWtNatXp/Phh7WIj7fy/vvBxMW50XsUOolk8m5Icmgb6YVNTXpRrs38kvs5\n0/YM4MHUDpwsz+TpZkt5vc1GEqPucrkBDxXDq3v3+syd60tQUD6XXWbh5ZcLsFolk7cXyeSFcHOZ\nZYf4Pusdfsx6nwaBLekXOYknm31GgMl91mgPDg7grrsa0q1bFm+/7ccXX+Tz0EPnvzArqkfiGiHc\nkEVb2Jz3Hd9lvs2uwl/pETGWxLp30TiordGl1VhpqZl//zuLVasimTmzlLvuquXybzruaJLJC+El\ncs3H+THrfb7Pfodw3/okRt1Nt/CbCDB53kXLrVtzeOstf9q3t/LRR6FERXlvuiyZvBuSHNpGemFT\nVS+01iTnreDlfaO4d2dbMs2HeLrZF8xqs54bIv/mkQMewMcnhddf98diKaJduzKWL6/BQjheTDJ5\nIVxUXnk2K7P/xfdZ8/FV/vSPuofJTd5zyQuojhIS4s+DD9bnxx/TGTGiDpMnn2T69DBZGuEiSFwj\nhAvRWrOr8Fe+y3qb307+j85hQ0ise7dL3NdutEOH8pk1y0psrGLp0lDCw71n0ksmL4SbK7Lkk5Tz\nId9mvkWZLiax7t30jhxPbd/zvJmqFyouNjNvXi7794fyzTc+JCR4x9KWksm7Icmhbby5F0dKdjH/\n8GQmbm9Ccv5yeueN5612v3NjvUe8fsBX9X0RFOTHY49F0717Pl27KhYvznd+YW5GMnkhnMyiLWw8\n+TVfZ87jQHEKfetOZG7bZKL8G7OtMOm8a7aLCjfdVI+4uBzuvDOY1NQTTJ1ax+iSXJbENUI4SV55\nNj9mvc+3Wf+kjm99BkbdT9fwUfiZvCNycIQDBwqYMUORmFjKBx9EeOwFWcnkhXBhe4s281XmPNad\n+IKrw4YyMOo+WoZcZXRZHiMnp5Tnny8hLs7K//4XRnCw5016yeTdkDfn0GfzxF6YrWWsyvmYx3+/\njhl7hxET0Iq32/3Bg3ELzzvgPbEXl6q6vYiICODFF0PJzbVy7bUF5OaWO7YwNyOZvBB2lF2WxvdZ\n8/k+6x0aB7XjxnqP0jlsMD5K/qs5UlCQD3//ewSvvZZNly6KVasCqVfPz+iyXILENULUkNaa1MKf\n+SrjDbbm/8j14WMYEHUfsUHtjC7N62gNb7yRze7dAaxa5UeTJp5xvUMyeSEMYLaWsib3E77MmEux\nNZ9BUZPpGTnOq55IdUVaw/vvZ7FhQzBr15po1izQ6JJqTDJ5NyTZq4279eKEOYNFx6YzcXscSTkf\ncmvD5/lnu10Mip5c4wHvbr1wpEvthVIwcWJdrr66mOuvt3DwYIl9C3MzEhQKUU37i5L5X+Zc1p34\ngq7ho3m+5XJig+KNLkucw4QJkZSX53D99f78/HMJjRq5/xn9pZC4RojzsGgLG05+xZcZc0gr3c2A\nqPvoV3cdHQrtAAAYDUlEQVSS1z+N6i60hrffzmHXLj9+/TWAevXcc/Y4PJNXSiUCc6iId97XWs88\nx34jgcXAlVrrzVV8Xoa8cAtFljyWZy/gq4x/UMs3kiHRD3Fd+Eh8ldyx4W60hrlzszl+3Jd160Ko\nVcv9AgyHZvJKKRMwD+gHxANjlFJtqtgvFJgMrLuUQryNZK82rtSL46X7eO/IQ0za3pRdBb/wcNyH\nzGq9nu4RY5wy4F2pF0azVy+UgilTIgkJsZKYmEdZmdUux3UX1bnw2hnYrbU+qLU2A4uAoVXs9zww\nEyi1Y31COJzWmu35q3hx7408sqszfiqAOW238nizT2gTKkv8egKTCZ54og7Z2b6MGpWN1eq8mNpo\nF4xrlFIjgH5a6zsrt8cCnbXWU87YpwPwjNZ6lFLqJ+ARiWuEqyvXZtbkfMLSjNcosxYzJPpBekbc\nRqBPiNGlCQcpLLTy+ONFDBxYxFtvRRtdTrXVJK6pTjhV1YFP/2RQFac5s4HxF/gaIVxCoeUkP2S9\ny5cZc2kY2IrbGs6gY+1EWf3RC4SEmHjuuQAefdSXtm0zmDLFfQb9parOkD8CxJ6x3QhIO2O7FhVZ\nfVLlwK8PLFNKDanqbP6NN+6gfv3mAISE1KFZsw4kJPQAbBmcN2yfmTe6Qj1Gbp/6O0e/3prNn/DL\nic9IjlxOp9r9uank7zSkFQlhrtOPffu2MnTogy5Tj5Hby5bNcdh8eOqpcp54IpnS0mIee2wIAElJ\nFZ/v0aOH4dtJSUksXLgQgLi4OGqiOnGND/A70Bs4BvwGjNFap55j/5+Ah7XWW6r4nMQ1lbZtSzr9\nzeftHN2LvUWbWZr+GpvyvqV35N8YEv0AUf6xF/5CA8j3hY2je7FyZT7/+pcv69dbadnStSM6Z91C\nORfbLZQvK6WmARu01l+dte9K4FHJ5IWRtNZszvuOLzJmkVbyB4OjH6Bv3Umy5ID4kwULTrBtm4nN\nm4OoVct1b4+VtWuEqGS2lrIq978sTX8NH+XLsOhH6Ro+Gj+TfM+Jv9Iapk07SUxMCcuWRbvsnVSy\ndo0bkvuhbezRi/zyHBYff5FJO5qyJvcTJjaaw5w2W+gZOdatBrx8X9g4oxdKwaOP1mL9+trMnHnc\n4a9nBPd79EuIMxwv3c//MmbzU85/uDpsKFObf0fT4MuMLku4kdBQE08+6cP//V8k11yTTY8ekUaX\nZFcS1wi3tL8omc/SZ7Il73v61J3I4KgpRPrHGF2WcGNff53PsmUmkpN9iI52rcXMJJMXXkFrzY6C\nNXyW/jL7i7cyOOpBEqPukoupwm5eeeUkAQGl/PhjlEvl85LJuyHJXm0u1AurtrLuxDIe/+Na3jh0\nB9fUuZF34vcxov7jHjfg5fvCxoheTJ5cm23bavH668ec/tqOIpm8cFlmaxmrcz/m8/SZBJiCGVHv\nSa6pcyM+ysfo0oSHCgpSPPywieeeq0vfvidISKhjdEk1JnGNcDnFlgJ+yH6PZemvExPYihH1nuTy\nWr1d6tdn4dn+/e98UlPL2by5Fv7+xp8LS1wjPEJeeRb/TXuOO3c0I7XgZ55q9jnPt1xOh9o3yIAX\nTnXrrbUoLfXjkUfSjS6lxmTIG0SyV5s1mz/h3cMPcveOVuSYj/Jyq7U82WwJLUOuNLo0p5PvCxsj\ne+HjA489FsiCBVGsWZNpWB32YPzvIcJrpZXs5tP0l/n50BISo+/ijbbb5DZI4TIaNPBlxIhC7rjD\nj23bzAQEuO6yB+cjmbxwuoPF21ly/EW25v/IgLr3MTh6CrV8I4wuS4i/sFrhkUcKGTo0h9dea2xY\nHXKfvHALe4o2sfjYDHYV/sLQ6IfoH3UPwT61jS5LiPM6dKicJ56wsmpVPp06GfM0rFx4dUPelL3u\nLPiZ5/b0Z8beoSTU6sG77fcxov4Tpwe8N/XiQqQXNq7Si9hYXwYOLONvf4PycovR5Vw0GfLCIbTW\nbM1bzjN/9GT2gdsqH2Day+DoKQSYgo0uT4iLMmZMKLm5Qbz6qvs9JCVxjbArrTUb877mk2MvUGQ5\nyaj6T3N9xBh8lFzjF+5t+3YzM2da2L69jEaNnBszOvo9XoW4IIu28OuJz1lyfAYazej6z9KlznB5\nOlV4jPbt/WjfvpQHHsjn009ruc2zGxLXGMRV8saasmgLq3L+y+Sd7Vma/hpjG77A3DZb6Ro+qtoD\n3lN6YQ/SCxtX7MWdd4bwww9RLF/uPmvPy5m8uCQWbWFt7icsOjadWr6RTGo8lw61+rjN2Y0QlyI8\nXDFyZDlTpgSQnGzG39/1752XTF5clLOH+5gGz9Ghliw7ILyHxQJTphRxzz3ZPPWUc+6dl/vkhcOd\nPdxvaTBNFg0TXis52cycOeXs2mWhbt1Qh7+e3Cfvhlwxb6zKqcz9/p3xfJ35Jnc1nsfMVmvtumiY\nu/TCGaQXNq7ci8sv96NxYyvPPptjdCkXJJm8qNLZZ+53NZ4nZ+5CnGHSpCCeeKI+Dz2URevWdY0u\n55wkrhF/cuZwr+1blzENnpPhLsQ5zJlTSEhILsuWxTj0/4jcJy9q7OzhLmfuQlzY7beHcM89/qxd\ne5xu3RoYXU6VJJM3iKvkjVZt5Zfcz5iSetnpzP3lVmuc+kYdrtILVyC9sHGHXtSpAwMGmHnsMRNW\nq9XocqokZ/JeqmL5gW/4KO3vmJSJCTGv0bF2PzlzF+IijR4dzMSJvnz1VRpDhjQyupy/kEzeCyXn\nr+SjtGcpsuRxa8PnuSZsmAx3IWrgs89KSEnJ47ffIvD1tf+5s2Tyolp2FfzKh8eeJbPsELc0mEbX\n8JtkbRkh7GDw4ECWLTOxdGkaI0fGGl3On0gmbxBn5o17i7Ywfc8gXj1wM9eH38Kb7XbSPeIWlxnw\n7pC9Oov0wsadeuHvDyNHaqZPD8JsNhtdzp/IkPdgh4p38vK+UTy/dyAdayfydrs/6Fv3DnyV66+3\nIYS76d8/gKNHw/j0U9dac14yeQ90rHQvHx97ji1533NjvccYGHWfvFGHEE7w9delrFlTwNattfHz\ns9/JlCxrIADIKjvCvIN38uiuq2kY0Ir58XsYXu8xGfBCOEm/fgGkp4eyeLHrnM3LkDeIPfPG/PIc\nFhx5nAdSL6eWbwRvx//BzQ3+7jZvku1O2aujSS9s3LEXvr4wfDi8+mog5eXlRpcDyJB3a6XWIj49\n/jL37GxNkfUk/2i7jfExL1PLN8Lo0oTwWv36BbB/fzjffusaZ/OSybuhcm1mefYCPjk2nTYh13Jr\nw+dpFNja6LKEEJU+/riEgwdzWbOmHiZTzc+l5T55L6G15ucTn/JR2rNE+jfi6WZf0DLkKqPLEkKc\nZciQQCZOjOLnn4/RrVuMobVU60eMUipRKbVLKfWHUuqJKj7/kFJqh1Jqq1LqR6WUc94uxY1dbN6Y\nnLeCR37vzGfHX+bOxvN4vsVyjxnw7pi9Oor0wsadexESAr16mZkxQ+PMtKQqFzyTV0qZgHlAbyAN\n2KCUWqa13nXGbpuBTlrrEqXU3cCrwM2OKNjb7C3azL+OPsnxsn3c1mAG14WPwqTkUooQrm7kyCDu\nuace27dnkJBQz7A6LpjJK6WuAaZqrftXbj8JaK31zHPs3wF4Q2vdrYrPSSZfTWkle/jo2LNsz1/F\nTQ3+j751J8pDTEK4mdmzi4iOzmLRopotdeDo++RjgMNnbB+p/LtzuQP49lKKEZBXnsU7h6fw2O/X\nEBvYnrfjdzMg6h4Z8EK4oZEjg/j66/qkpRn3NoHVufBa1U+PKk//lVJjgU5A93Md7I037qB+/eYA\nhITUoVmzDiQk9ABsGZw3bJ+ZNyYk9KDMWsK7ax5gde4iena6jTfb7eTwrp3sydroEvU6cvvsnhhd\nj5Hb+/ZtZejQB12mHiO3ly2b4xHzIS6uM3PnnqR//xQAevSo+HxSUtI5t5OSkli4cGHl18dRE9WN\na57TWidWblcZ1yilbgDmAtdrrbPPcSyJaypt25ZEQkIPrNrK6tyP+U/a0zQP6si4mJe97nbIU70Q\n0oszeUovNmywsGBBIbt2+RIScmlPn9ckrqnOkPcBfqfiwusx4DdgjNY69Yx9rgCWAP201nvPcywZ\n8mdIyf+JBUcexaR8mNDoNeJD/3IZQwjh5qxWuOuuUl54IY2JE5te0jEcep+81tqilLof+IGKDP99\nrXWqUmoasEFr/RXwChACLFEV7z5xUGs97FIK8gaHinfyr6NPcLBkO+MavkTX8NFyx4wQHspkgkGD\nFG+9FcTtt5c75E1FzkeeeHWiXPNxPj72HL+c+IxrT4xkUtc5+JkCjC7LcJ7ya7k9SC9sPKkXxcUw\nYUI53313nK5dL/4tAmUVShdXYilk0bHnuX9nPAGmYN5q9zvdwm+SAS+ElwgKguuvN/Pqq85/s285\nk3cgq7byU85/+DDtGdqGXMdtMS/SIKC50WUJIQxw/Dg89JCZ1NSTNG5c96K+VtaucUE7C37mvSMP\nYsKHJ5ouoU1oF6NLEkIYqH59aNrUzPz5ebzwwsUN+ZqQuMbOMkoP8ur+m5m1/2aGRD/IK61/qXLA\nu/O6HPYmvbCRXth4Yi8GDgzgk0/qUFJS4rTXlCFvJ8WWAj5M+zsP7epITEAb/tluFz0ibpW7ZoQQ\np119tQ8ZGbVYsSLdaa8pmXwNWbWVpJwP+U/a07QP7c64mJeJ8pdFOIUQVfvggxJKSzP4+uvGVNxx\nfmGSyRskteAX3jvyIAolubsQoloGDAjkoYcacPhwJrGx0Q5/PckSLkFm2SFe3T+GV/aPZlDUFF5p\n/etFD3hPzBsvlfTCRnph46m9qF8f4uLKeffdAqe8ngz5i1BiKeSjtP/jwdQriAloxVvtfqdn5FjJ\n3YUQF2XgQH8WL65DaWmpw19LMvlq0FqzNncxC44+SrvQboyPeZko/5qtDy2E8F7l5XD77WaWLDlG\n374XniWSyTvQweLtvHN4MgWWXB5p+l9ZREwIUWO+vtCtWzlvv63p29exryU5wzkUlJ/g3cMP8uzu\nXlwXPorX22y064D31LzxUkgvbKQXNp7ei379gvjpp3pkZ+c69HVkyJ/Fqq0sz17AfTvbUmotYl7b\nHQyIuhcfJb/0CCHsJy4OatWCJUsc+65RksmfYXfhRuYfvh+AuxrPo2XIlQZXJITwZJ99VsaePZms\nWlUfHx+fc+4nmXwN5ZVn8e+jT7Ph5JeMi3mJnhHj5I4ZIYTD9erlz+LF9di3L52WLRs65DW8epJZ\ndDlfZ77JfTvbEWAK5s12qfSOvN0pA97T88aLIb2wkV7YeEMvwsOhZUszH3xQ5LDX8Noz+Z0Fa5l/\n+H5CfcJ5oeVKmgS1N7okIYQX6tPHn2XL6vDcc6UEBNj/PSa8LpM/Yc5gwdHH2Ja/kr/FzKJr+Ohq\nrx8hhBD2VloK48eXs3x5Gl26VH3PvLwzVDVYtZXvMuczObU9Yb5RzGu3k24RN8mAF0IYKiAArrrK\nzPz5Focc3yuG/N6iLTz++7WszPk301ssZ0KjWQT71DK0Jm/IG6tLemEjvbDxpl707RvIihWRFBYW\n2v3YHj3kiyx5vHv4QZ7bk0jfupN4udUamgZfZnRZQgjxJ+3aKfLzg1i9OsPux/bITP7UWjMfHH2E\nK2r34/aYmdT2dd7bbQkhxMWaN6+E8PBjLFrU9C+fk/vkz5BWspv5h+8nx5zGY00X0S60q9ElCSHE\nBfXqFcjcuVHk5eVRu3Ztux3XY+KaMmsJ/017jsd+70KH2n2Y3XazSw94b8obL0R6YSO9sPG2XrRp\nAyUlAaxcmWnX43rEkE/OW8Hk1AQOFKcwp+0Wbqz3KL7Kz+iyhBCi2kwm6NLFyocf2veOP7fO5PPK\ns/ngyCOk5K/k7sZv0rnOYLscVwghjPD77zBrViG7dpVTp07Y6b/3uvvktdYk5XzE/TvjCfGpw7x2\nO2TACyHcXqtWUF7ux4oV2XY7ptsN+eOl+3huTyKfp7/Cs83/x6TGcwy/5/1SeFveeD7SCxvphY03\n9kIpuPZabdfIxm2GvEWX83n6qzyyqzOX1erN62020iqks9FlCSGEXfXsGcDPP9fjxImTdjmeW2Ty\nuws3Mu/QJGr71uXe2LdpENDcQRUKIYSxtIYJE8qYP/8II0Y0Azw4ky+2FPD+kYd5fu8ghkU/zPQW\nP8iAF0J4NKWgc2crixfb53guO+Q3nfyWyantySvPZl677fSMvM2jFhPzxrzxXKQXNtILG2/uRbdu\ngfz6a13y8/NrfCyXe+I1vzyH9448yM6Ctdwf+x4dat9gdElCCOFUbdtCbm4wGzYcplevmt1Y4lKZ\n/C+5n/HOkclcV2c0tzWcQaBPiNNqE0IIVzJrVgmtWh3h3XdbuP/aNbnmdOYfvp+DxSk83nQJ7UKv\nM7okIYQwVNeuAXz6aTjFxcU1Ok61MnmlVKJSapdS6g+l1BNVfN5fKbVIKbVbKfWrUqrqtzc5y6mH\nmqakXkaDgObMabvVawa8N+eNZ5Ne2EgvbLy9F1dcoTh0qA6pqTVbfviCQ14pZQLmAf2AeGCMUqrN\nWbvdAeRorVsCc4BXLnTc7LKjvLB3CJ8dn8n/Nf+a8TEvE2AKuvh/gZvat2+r0SW4DOmFjfTCxtt7\nERAA8fFmFi0y1+g41TmT7wzs1lof1FqbgUXA0LP2GQr8q/LPnwK9z3UwrTU/ZL3HA7s60CLkSl5v\ns5GWIVdeSu1urbDwhNEluAzphY30wkZ6Addd58+PP4bW6BjVyeRjgMNnbB+hYvBXuY/W2qKUOqGU\nitBa55x9sGn7+1Os83i+xQp5lyYhhDiPzp1NvPNOzd7wqDpDvqorumffknP2PqqKfQCID+zGoMgp\n+GhfCgvt89iuO0pL+8Or//1nkl7YSC9spBfg4wOxsb788celH+OCt1Aqpa4BntNaJ1ZuPwlorfXM\nM/b5tnKf9UopH+CY1jq6imM5735NIYTwII68hXID0EIp1QQ4BtwMjDlrny+B8cB6YBSw0p5FCiGE\nuDQXHPKVGfv9wA9UXKh9X2udqpSaBmzQWn8FvA/8Rym1G8im4geBEEIIgzn1iVchhBDO5ZAFyhz1\n8JQ7qkYvHlJK7VBKbVVK/aiUamxEnc5woV6csd9IpZRVKdXRmfU5U3V6oZQaXfm9sU0p9aGza3SW\navwfaayUWqmU2lz5/6S/EXU6mlLqfaVUulIq5Tz7/KNybm5VSnWo1oG11nb9oOIHxx6gCeAHbAXa\nnLXPPcA/K/98E7DI3nW4wkc1e9EdCKz8893e3IvK/UKBVcAvQEej6zbw+6IFsAmoXbld1+i6DezF\nfOCuyj+3BfYbXbeDetEV6ACknOPz/YGvK/98NbCuOsd1xJm8XR+ecnMX7IXWepXWuqRycx0Vzxx4\noup8XwA8D8wESp1ZnJNVpxeTgDe11nkAWussJ9foLNXphRWoXfnnOsBRJ9bnNFrrtUDueXYZCvy7\nct/1QJhSqt6FjuuIIV/Vw1NnD64/PTwFnFBKRTigFqNVpxdnugP41qEVGeeCvaj89bOR1vobZxZm\ngOp8X7QCWiul1iqlflFK9XNadc5VnV5MA25TSh0GvgImO6k2V3N2r45SjZNCR6xCadeHp9xcdXpR\nsaNSY4FOVMQ3nui8vVAV7wgzm4pbcc/3NZ6gOt8XvlRENtcDscAapVT8qTN7D1KdXowBFmitZ1c+\nt/MhFetoeZtqz5MzOeJM/ggV35SnNALSztrnMNAYoPLhqdpa6/P9muKuqtMLlFI3AE8Bgyt/ZfVE\nF+pFLSr+4yYppfYD1wDLPPTia3W+L44Ay7TWVq31AeB3oKVzynOq6vTiDmAxgNZ6HRColKrZs/7u\n6QiVc7NSlfPkbI4Y8qcfnlJK+VNxz/z/ztrn1MNTcJ6HpzzABXuhlLoCeBsYorXONqBGZzlvL7TW\neVrraK11M611UyquTwzWWm82qF5Hqs7/kaVAL4DKgdYS2OfUKp2jOr04CNwAoJRqCwR48DUKxbl/\ng/0fMA5Or0RwQmudfqED2j2u0fLw1GnV7MUrQAiwpDKyOKi1HmZc1Y5RzV786Uvw0LimOr3QWn+v\nlOqrlNoBlAOPeuJvu9X8vngUeFcp9RAVF2HHn/uI7ksp9V+gBxCplDoETAX8qVhG5h2t9TdKqQFK\nqT1AIfC3ah238nYcIYQQHsghD0MJIYRwDTLkhRDCg8mQF0IIDyZDXgghPJgMeSGE8GAy5IUQwoPJ\nkBdCCA8mQ14IITzY/wMrDbByV/VblgAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7fe907b69c88>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "plot_x = np.linspace(0, 1, 200)\n",
        "\n",
        "pt.plot(plot_x, f(plot_x), label=\"f\")\n",
        "pt.fill_between(plot_x, 0*plot_x, f(plot_x),alpha=0.3)\n",
        "pt.plot(plot_x, int_f(plot_x), label=\"$\\int f$\")\n",
        "pt.grid()\n",
        "pt.legend(loc=\"best\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This here plots the function, the interpolant, and the area under the interpolant:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 133,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<matplotlib.legend.Legend at 0x7fe907a37d30>"
            ]
          },
          "execution_count": 133,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8TfcbwPHPN5GYCWIkQgixR+0ZI8Mm6LBXF6q12uqP\nUpuiRs0WXSiqRYdZo8Teo/YmhFixYkaS8/vjhotmXMm999x787xfr7yak3zvOc99Gk9OnvM936M0\nTUMIIYRjctI7ACGEEJYjRV4IIRyYFHkhhHBgUuSFEMKBSZEXQggHJkVeCCEcWLJFXin1g1LqqlLq\nYBJjpiilTimlDiilypk3RCGEECllypn8T0CDxL6plGoE+GmaVgToBswwU2xCCCFSKdkir2naFuBW\nEkOaA3Pjx+4EsiqlPM0TnhBCiNQwR08+L3Dxue1L8V8TQgihM3MUeZXA12StBCGEsAHpzLCPcMDn\nue18wOWEBiqlpPgLIUQKaJqW0Al1skw9k1ckfMYOsBToBKCUqgbc1jTtaqJ7GgoMADe/1/nw569o\nt7gdVWZVIddXucg+JjuN5jVi/NbxHLp6iLi4ODRNc8iPIUOG6B6DrXxILiQXkoukP1Ij2TN5pdQC\nIADIoZS6AAwBXAFN07RZmqatVEo1VkqdBu4D7yS5w2jw2u1FQKMu/NwzmNffPcP/3j7FY+0eF+5e\n4OSNk/x5/E/GbRtHTFwMdQrUoVHhRtQvXJ/8WfOn6s3akvPnz+sdgs2QXBhJLowkF+aRbJHXNK2d\nCWN6mHrAWuG16NivI175YghusYnpg8uy5W9veo06QKkSbpTKVYo4LY4HTx5w7vY5jl4/ysx9M+m7\nti/ZMmQjqGAQjYo0IrhgMB4ZPUw9rBBCpEkqtX8KvNLBlNKWHFmCi7PLs69pGqz/Ix+zx5ekfssL\ntOp+kvQZ4l54XZwWx93HdzkReYLjN45zIvIEZ2+dxS+7H02KNOH1Eq9TybsSTsp+buANDQ0lICBA\n7zBsguTCSHJhJLkwUkqhpbAnr3uRf+rW9fTMHFmas8ey0n3wIcrXvJ7ofmLiYrj58CaHrx/m8LXD\nHL9xnPvR96lbqC4tiregUeFGZM+Y3ZJvRQghrMYhivxTezbmZuaI0hR97Tbv9T+CR+7Hye73wZMH\nnL11loPXDnLk2hHO3DpDqVylaFq0Ka8Xf53SuUujVIryYzFylmIkuTAyVy58fX0JCwtLfUDCqgoU\nKJDgtYjUFHlzTKE0q0p1rlGmykZ+/bYIvZrXoV3PEzRoHYazc+KvyeSSidK5S1M6d2li42K58eAG\nB64eYP259UzfPR0n5URI0RDalGpDbd/apHOyubcthFmFhYWlelaGsD5LnIza3Jn88y6cysI3Q18j\n5okT3YcexK/k3Vc6nqZp3Iu+x/HI4xy4coCDVw9y5/EdGvo1pHXp1jQs3JAM6TKk5K0IYdPiz/z0\nDkO8osT+vzlUu+ZlcXHwzx8+zJ1YgoCm4bTteYJMWWJTdPzHMY85ffM0eyP2cvDaQS7dvUSgbyCt\nSrWiWbFmZM2QNUX7FcLWSJG3T2myyD9195YrP31Vgv1bc9G57zECQi6Rmr9snsQ+4eKdi+yO2M3B\nqwc5c+sMVfNWpX2Z9rxV8i2LF3zpQxtJLozMlQsp8vbJEkXebuYcumePpvfof+k/eS/L5hSif3t/\nzh5zT/H+XJxdKORRiNalWjM8YDjj643HN5svU3dNJe/EvDSc15D5B+dzP/q+Gd+FEALg5MmTVKhQ\ngaxZszJt2jS9w3FodnMm/7zYWFi3JD/zpxSjer0rtO91HPfsT8wSY2xcLJejLrPj0g72Ruwl7HYY\ndQvVpVPZTjQp0oT06dKb5ThCWJKtn8m///77ZM2alQkTJugdik1J0+2ahETddmH+lGJsXe1Nu54n\nqN8y6Vk4ryo2LpawO2HsCN/Bviv7iIiKoEmRJnQu15l6herh7GTGgwlhRrZe5OvVq0fbtm159913\n9Q7FpkiRT8S54+7MHFmaRw/S0fWLQ5SskNQzTlLmSewTzt46y85LO9kbsZeo6Cjalm5LlwpdKONZ\n5pX3J31oI8mFUVroyQcHB7Nx40ZcXFxwcXFh3759FC5cWO+wbIIlirxDTBgvWPwuo3/exqYV3oz/\npCLFy9+k06fH8Mr30GzHcHF2oVjOYhTLWYzWpVpz5PoRtl3cRsDsALyyePF2ubd5p/w75MyU02zH\nFMIR/fPPPwQGBtKxY0c5k7cChziTf97jh8788VMhls0pRL2WF2jZ7RSZ3WIsdryox1HsDN/J9kvb\nOXr9KDV8atClQheaF2+Oq7OrxY4rRFJMOZM31303KSkhUuQTlqZn15gqfcZY2nx4iilLQ7lz05Xu\njQJZtbAAsTGWWdbALb0bdf3q8kWtL5hQfwKemT35Yv0XeI735P2l73Mg4oBFjitEammaeT6EbXO4\nIv9UDs/H9P7yX4bM2smWVd70alGHvZtzWex4Sim83bxpVaoVo4NH81n1z4iIiqDuz3V57dvX+Gb3\nN9yLvvdsfGhoqMVisTeSCyPJhTA3h+jJJ8Wv5F1Gzt7Org2efDeyNMvyP6Bz36MULBZlsWO6OLtQ\nxrMMZTzL0L5Me7Ze3MqsvbPov64/LYq3oFfVXhY7thBCPM/hevJJeRKt+HuhL4tmFaa8/3Xa9TyB\npxkvziYlNi6Ws7fOEhoWyvbw7eTOnJsPKn7A2+Xexi29m1ViEGmHLc+uAQgKCqJDhw7Sk3+JTKE0\nkwf30vHnj4VYMb8gAc3DafXBKbJ6RFvt+Pej77Pt4jY2X9zM6ZunaVasGZ9W+5QK3hWsFoNwbLZe\n5EXC5MKrmWTKEkO7XieZtmIDcbGKDxsHsnB6ER7et87NTZldM+N1w4uhdYYyImAE9x7fo/68+lSY\nWYG5/84lOtZ6v3BsgfShjSQXwtzSZJF/KnvOaLoNOsyERZu5dC4LHzQMYvk8X55EW+cBI07KCT8P\nP96v8D5TGk6hsndlvtz8Jd4TvOm/rj8RURFWiUMI4bjSZLsmMWeOujN3YgkiwjLTvvdxajW+jJOV\nfw3GxMVw8NpB1p1dx/6I/QQXCuazGp/h7+Nvc0+3ErZL2jX2SXryVnJwRw7mTizB40fOtOt5gmp1\nr5jtxpFXcfnuZdaeW8uWC1vwyOhBn2p96FS2kzzoRCRLirx9kiJvRZoGuzd4Mn9qMZSCdj1OUDnw\nqtmK/aGdhyhT1bQ1bx48ecCWC1vYcH4DEfci6FaxG59U/8RhllCQtWuM0sLaNSJxcuHVipSCKkFX\nmfT7Jlp3P8nPk4rzacta7N2U2+p3+WVyyUR9v/qMDBxJ76q92Xh+I76TfOn8Z2dO3jhp3WCEEHZF\nzuRNFBcH29bk4ZdpxciU+Qntep2gXI0burRxNE3j/O3zrDq9im3h26jiXYWBtQdSp0Ad6dsLQM7k\n7ZW0a2xAbCxsWeXNwmnFyJrjMe16nuC1apG6xXPjwQ1Wn1nNhvMbyJ05N/38+9GmdBvSOTn8zcwi\nCVLk7ZO0a2yAszPUaXqZactDadDqAtMGl+XzDjXYvyXXK7VxDu08ZJZ4cmbKSfsy7ZnScAr+Pv4M\nDR1K/q/zM37beB4+sc7dvKklc8ONJBe2JTAwkB9//FHvMFJFinwKOafTCGwezrcrN1C/VRjfjy5F\n31Y12fGPJ3Fx1o8nk0smmhRpwrh64+hUthM/H/yZfF/nY0joEO4+vmv9gIRIQsGCBVm/fn2y4xyh\nyAJs3LgRHx8fXY4tRT6VnNNpBDa7xNRlobz5/hkWTi9G7xZ12LTCm9jYxF9n6syaV5XOKR3+Pv4M\nqzOM3lV6s+rUKny+9uGT1Z9w48ENixwztWRmjZHkwrzi9DjjSoCmabpdL5MibyZOTlCjQQRfL9lE\n50+PsXxeQT5qHMi6JT7EPLH+/1xnJ2fK5ynPF7W+YGDNgey6tIuCkwvSbVk3wu+EWz0eIRIyZ84c\natWqxWeffYaHhwd+fn6sXr0agC+++ILNmzfTo0cP3N3d6dXLsHrr8ePHqV+/Pjly5KBEiRIsWrTo\n2f7eeecdPvzwQ5o0aYKbmxuhoaG88847dO/enfr16+Pu7k5gYCAXLlx49ppt27ZRpUoVsmfPTtWq\nVdm+fXuCsZ49e5bg4GBy5sxJ7ty56dChA3fvGv9KLliwIBMmTKBs2bJkz56dNm3aEB0dzYMHD2jc\nuDGXL1/Gzc0Nd3d3rly5Yol0JkiKvJkpBZXqXGPsgq18NPwgG5fnpVuDIFYuKED0Y2O6zdWTTz4e\nRYlcJfhfjf8xImAEp2+epsQ3JejwewdOR562SgzJkT60UVrMxc6dOylRogSRkZF89tlnz1amHDly\nJLVq1WLatGncvXuXKVOm8ODBA+rXr0+HDh24ceMGv/zyCx9++CHHjh17tr9ffvmFQYMGERUVhb+/\nPwALFixgyJAhREZGUrZsWdq3bw/ArVu3aNq0KX369CEyMpKPP/6YJk2acOvWf58TrWkaAwYM4MqV\nKxw7dozw8HCGDh36wphFixaxZs0azp07x8GDB5k9ezaZMmVi1apVeHt7ExUVxd27d/Hy8rJQNv9L\npmBYiFJQpmokZapGcvxANhbNKMLCb4rStMN5GrU5r0M8Cj8PP/pU60P4nXCWnVpGhVkVqFuoLqOD\nR1MsZzGrxyT0pYaZ5y9MbUjqZvH4+vo+K+ydO3fmww8/5Nq1a+TOnfs/Y5cvX07BggXp1KkTAOXK\nlePNN99k8eLFDBo0CIDmzZtTrVo1ANKnTw9AkyZNnhX8UaNGkS1bNi5dusSGDRsoWrQo7dq1A6BN\nmzZMmTKFZcuWPTvGU35+fvj5+QGQI0cOPv74Y4YPH/7CmN69e+Pp6QlASEgIBw7o/2Q4KfJWULzc\nbQbN2M2FU1n440c/ujUIIrB5UTzznSV3XuvPgMmXNR/dK3Xn9eKvs+zkMqp8X4W6Besyuu5oiuYo\navV4pA9tZM1cpLY4m8vzZ7UZM2YE4N69ewkW+bCwMHbs2IGHhwdgOLuOjY19oSAndIHz+a9lzpyZ\n7Nmzc/nyZS5fvkyBAgVeGFugQAEuXbr0n31cv36dXr16sXnzZu7du0dsbOyzOJ56WuABMmXKRESE\n/osMSrvGivIXuUfv0f8ydelG0rnE8fEbtZnQtzxnj7nrEo9XFi+6VOjChHoTeBz7mEqzKtHyt5ac\nvmkbbRwhXr5Y6ePjQ0BAADdv3uTmzZvcunWLu3fvMm3atERfA3Dx4sVnn9+7d49bt27h7e2Nt7c3\n58+ff2HshQsXyJs373/28fnnn+Pk5MThw4e5ffs28+bNM/leBD1vUpQir4Mcno+oVPs3Zq37h4Il\n7jLigyoMerfaK8+1NxfPLJ50q9iNifUn8uDJA8rPLE/rRa05c/OMVY6fFvvQiZFcvMjT05OzZ88+\n227atCknT55k3rx5xMTE8OTJE/bs2cOJEyeS3M/KlSvZtm0b0dHRDBo0iKpVq5I3b14aN27MqVOn\nWLhwIbGxsfz6668cO3aMkJCQ/+wjKiqKLFmy4O7uzqVLlxg3btwrvY/IyMgXLtRaixR5HWV2i+GN\n984wa+0/BDQN54exJenzem02LM1rtTXtn+eZxZMPKn3AxPoTiYqOovzM8rRd0pZzt85ZPRbh2JI6\ns33+e71792bRokXkyJGDPn36kCVLFtasWcPChQufnYn379+fx48fJ3m8du3aMXToUHLkyMH+/fuZ\nP38+AB4eHixfvpzx48eTM2dOxo8fz4oVK8iePft/YhkyZAh79+4lW7ZshISE8Oabb5r8nooVK0bb\ntm0pVKgQHh4eVp1dY9KyBkqphsAkDL8UftA0bexL3/cB5gDZ4sd8rmnaqgT2Y/fLGliSpsHeTbn5\na3YhLp52o1G78zRsHWbVRxM+LyIqgiXHl7Dt4jZeL/46o4NHk9f9v3/GCtsjyxoYvfPOO/j4+Pzn\nIqkt0mVZA6WUEzANaACUAtoqpYq/NOwL4FdN0yoAbYFvUhJMWvd0+uWIn3Yw9PsdXL+ckQ8aBjJl\nYFnOnbD+w77zuOWhR+UejKs7jstRlyk+vTg9V/Xk1sP/Ti8TQtgmU9o1VYBTmqaFaZr2BFgINH9p\nTBzw9OphNuC/l6bFC5KbJ+9bLIoeIw4y4+8NePk8YFiXqnzxdjV2rbf+sgl53fPSu2pvRgSOYN/l\nffhO9mXwhsE8ePLALPuXPrSR5ML80vrKrMm2a5RSbwINNE3rGr/dAaiiaVqv58Z4AWuA7EAmoK6m\nafsT2Je0a+K9ykNDAJ5EK7at9mbp3IJE3XElpMM5gt+4QKYsSaydYAGapnHk2hEWHVvEhTsXGFBr\nAD2r9EzV/1N5aIiRPDQkbbNEu8aUefIJ7fjlKNoCP2ma9rVSqhowD0Nr5z+mDpiKVz7DvNjMbpkp\nVKLQs2L39Ow2LWyXqVrmlca7uGp45P6bzn3B1bU2S+cW5OfJEZStdocOvfNQoGiUVeMvmbskS1Yt\nYdLCSUzcPpGRQSPJfys/TsrpWZF6elYq26+2/VRq9yfsV2hoKLNnzwYMN4ulhiln8tWAoZqmNYzf\n7g9oz198VUodxnC2fyl++wxQVdO0Gy/tS87kzSjyagZW/5afNb8VIE+B+zRud55qdSNwcbXeGVxs\nXCybLmzij2N/4JrOlfH1xhNS7L/Tz4R1yZm8fdJrPfndQGGlVAGllCvQBlj60pgwoG58MCWA9C8X\nePEic6xdk8PzEe16nuT79eto2uEcf/9agPeD6zJvcjGuR1jnYd/OTs4E+gYyof4EAgoE0GVZF/x/\n9GdfxD6T9yF9aCPJhTC3ZNs1mqbFKqV6YOi5P51CeUwpNQzYrWnacqAv8J1S6mMMF2E7WzJo8aJ0\nLhr+DSPwbxjBhdNZ+HthAXq3qEPpypE0anuestVv4GThOyJcnF1oWrQpgb6B/H78dwLnBNLQryET\nG0yUaZc6KFCgQJq/4GiPXl5iwRzk8X8O6uF9ZzYuz8vKBb5EP3KmUdswglpcxC3bE6sc/9q9ayw8\nupAd4TvoXqk7Q+oMIbNrZqscWwhHI894FYnSNDi+PzsrF/iyZ6Mnlepcpd5bFyhTNdIqDyE/FXmK\nBYcXEHYnjGEBw+hWsRvOTs6WP7AQDkSKvB161SmU5nD3lguhy/KxZlF+nkQ7Ue/NiwS/fpHsuZK+\nJTy1NE1j56WdLDyykHQqHZMaTqJRkUbPvi9TKI0kF0aSCyNLT6EUDsI9+xOadTpHSMdznPg3G2sX\n5+fDJgGUqRJJ/ZYXKF/zGs4WOMlWSlEtXzUq5anEqjOr6PxnZ0rlLsW3jb+leK6Xb54WQpiTnMmn\ncQ/uObN5ZV7WLs7PzWsZCH7jIvXevGDRde7vPb7Hb8d+Y93ZdXR8rSNj647FLb31l20Qwl5Iu0aY\nxbkTbqxdlJ+Ny/NRuNRt6r55karBV3BNb5l1FMLuhDHv4DxO3zzNyKCRdK3YFSclC6MK8TIp8nZI\nj568qR4/cmL72jys/8OHM0eyUqNBBMGvX6RYuVtmv1iraRqLVi5is9Nm3NK7MaPJDPzz+5v3IHZE\n+tBGkgsj6ckLs0qfIY6AkEsEhFziekQGQpfmY/KAcmhxEPR6OIHNwsnlbZ52jlKKEjlL0KJSC5ae\nWErTBU2p61eXKQ2nkMctj1mOIURaJmfywiSaBif+zcb6P33YusqbQiXvENQinOr1IsiQyXyLpEU+\niOTnQz+z5/Ie+lbvS7+a/XB1djXb/oWwR9KuEVYV/diJXes9Wf+nD8f2Z6da8BWCWoRTqnKk2e6s\nPXLtCHMPzuXBkwfMaDqDhoUbmmfHQtghKfJ2yJZ78q/i5rX0bFyWl/V/+XA/yoXaTS5RJyQc36JR\nJvfvE8tFTGwMa8+t5dcjv+Lv48+MpjPwdvM28zuwLdKHNpJcGFl6gTIhEuWR+zGvv3eWqUs3MnjG\nTpTSGNm9Cj2b1WHRzMJcDc+Y4n2nc05Ho8KNmNRgEo9iHlFiegnGbBlDbJx119AXwp7Jmbwwu7g4\nOL7fg9Bledm2Og95C96nTtNwajaKwD17yp9Xe/jqYX7890eclTPfhXxHrQK1zBi1ELZL2jXCZj2J\nVhzYmouNy/Oyd5MnJSrcpE7TS1QNvpKiC7YxsTEsPbmU34//TtMiTZnaaCoemTwsELkQtkOKvB1y\nlJ78q3h435md673YuCwvx/Z5UKHWNWo2jCBDpnVUqJXgg8QSde3+NWYfmM3h64cZHTyabhW7OcTS\nutKHNpJcGMk8eWEXMmaOfTb//u4tV7av9WLVrwU4sb8SlQNLULPRZSrUukb6DMnfYZs7c27+5/8/\ndobv5MvNX/Ldvu+Y02IOpXOXtsI7EcJ+yJm80N2dm65sX5OHLX/n4cyRbFSqcxX/RpepUPO6SUsq\nPIp5xKKji1h1ehXdK3VnZNBImVsvHIq0a4TDuHXDle1r87B1lTdnj7tTOeAaNRtepnzN67i4Jl3w\nz946y3f7vuPBkwf82OxHAgoGWCdoISxMirwdSos9+cQklotb19OzbU0etqzyJuyUG5UDrlKjQQTl\nalxPtKUTGxfLytMr+e3IbzQr1oxpjabhnsHd0m/BbKQPbSS5MJJ58sIhZc/1mCbtzzN63jamLg2l\nSOnbLJ1diM416zOmd0U2Lvfmwb0XLys5OzkTUjSEifUnciryFIWnFua3I7/p9A6E0J+cyQu7c+em\nK7s2eLJ9TR6O7PGgZMWbVK93hSpBV8iWwzgPX9M0Nl/YzJx/51AxT0V+bP6jLHom7JK0a0Sa9eBe\nOvZuys32tV7s35Ib32J3qV4/gup1rzxbKTMqOorZB2az89JORgaOpEeVHg4x3VKkHVLk7ZD05I3M\nlYvox078uy0n29bmYfcGT3J7P6RavQiq17uCj9899l/Zzw/7fsDb3ZsFbyygYPaCZojevKQPbSS5\nMJJ58kIArunjqBx4jcqB14iNURzZ48H2tXkY/G41MmaOoWpwcbrU8WeXmkrZGWUZGjCUFl4tGPz1\nYC7dvURe97yM+GQEBX1tr/gLkVJyJi8cXlwcnD6cjZ3/eLFrvSe3b6ancL31nM3zIQ/2RPC41mNw\nBaLB718/1k5bK4Ve2BRp1wjxCq5czMSuDZ4sWjyYO2/9YSjwT0VD+6j2zJsyT7f4hHiZTKG0Q4d2\nHtI7BJth7Vx4+TygWadz5C97+sUCD+AK2w5e5uBBw9OwrC00NNT6B7VRkgvzkCIv0iyPTB7w8srH\n0RAVc4fmLeIoWBB69YJ16yA65SskC6EradeINOtK+BUGfzWYK5WuPOvJe+zwIKN/RvLky8OwMr9w\nYENhli6FEyegfn0ICYHGjcFDVjcWViQ9eSFS6Er4Feb/PJ+bD27ikcmD9h3bk90rOwsOL2Dd2XWM\nChpFjyo9uHpVsXw5LFsGGzZA+fKGgh8SAsWK6f0uhKOTIm+HZJ68ka3m4uDVg8zYM4PCOQqz4I0F\nz+6WffgQ1q83FPxlyyBzZmPB9/cHl1T8eMvccCPJhZFceBXCAl7zfI3x9cbj4uRCyW9KsuDQAgAy\nZoQmTWDGDAgPh4ULwd0d+vYFT09o1w5++QVu3dL5DQiBnMkLYZLtF7fz3f7vqJ2/Nj82/5GsGbIm\nOO7SJVixwnCGv3EjVKxoPMsvUsTKQQuHIe0aIazgzqM7fLv3W87cPMOcFnNoULhBkuMfPDDMzFm2\nDJYvh2zZjAW/enVIJ/ebCxNJu8YOyTx5I3vJRdYMWenv359WpVrRenFrui7ryqOYR4mOz5QJmjWD\n774znOHPnQsZMhimZXp5QceO8NtvcOeO8TUyN9xIcmEeUuSFeEX1CtVjfL3x7Ly0k5LTS7Ln0p5k\nX+PkBJUrw/DhsH+/4aNGDZg9G3x8oG5dmDwZLl+2fPwibTGpXaOUaghMwvBL4QdN08YmMKYVMASI\nA/7VNK1DAmOkXSMcRmxcLH8c/4M/jv9BP/9+DKg1IEVLGN+7Z2zrrFgBOXIYWjrNmkHVquDsbIHg\nhV2xaE9eKeUEnASCgcvAbqCNpmnHnxtTGPgVCNQ07a5SKqemaTcS2JcUeeFwTkWeYsquKeRzz8ei\nlotS9WCSuDjYvdtQ8JcuhStXDDdfhYQYbsZyczNj4MJuWLonXwU4pWlamKZpT4CFQPOXxnQBpmua\ndhcgoQIvXmQvfWhrsPdcFMlRhHH1xpHFNQulvinFkqNLUryvTZtCqVoVRo6EgwcNBb9SJZg1C/Lm\nhQYNYNo0CAsz4xuwUdKTNw9Tinxe4OJz2+HxX3teUaCYUmqLUmqbUirpaQdCOJgM6TLwUeWP6Fqh\nK12WdeGdP9/h0ZPEL8qaqkAB6NEDVq82XLzt2hX27DH09197DQYOhL179VlMTdgHU9o1bwH1NU3r\nGr/dAaisaVrv58Ysw7DUU0sgP7AZKPX0zP65cdKuEQ7vxoMbTN45mXvR91jcajHlvMqZ/RixsbBz\np6Gls2QJxMTAW28ZPqpUAXm6oWOx9JOhwjEU7qfyYejNvzxmu6ZpccB5pdQJoAiw9+WdTR0wFa98\nXgBkdstMoRKFnt3S/vTPdtmWbXvfHlZnGNMXT6fGoBoMf2c4n1b/lI0bNwI8u1X/aTsiJdvOzhAd\nHUrDhjB6dACHDsH48aG0bAkQwJtvgp9fKCVLQlBQ6o8n29bdDg0NZfbs2QD4+vqSGqacyTsDJzBc\neI0AdgFtNU079tyYBvFfe1splRNDcS+nadqtl/YlZ/LxbHW9Fj04ci6O3zjO1F1T8fPw47e3fiNX\n5lxJjjfHei1HjsDixYaPmzfhzTcNZ/j+/vY1U0fWrjGy6IVXTdNigR7AGuAIsFDTtGNKqWFKqabx\nY1YDkUqpI8A/QN+XC7wQaVHxnMUZX288TsqJkt+UZN3ZdRY/ZqlSMGQIHDoE//xjWE+nVy/In9+w\nvs7+/dLDT0tkWQMhrGTd2XXM/nc2H1T8gNF1R+OkrHsv4rFjMH++4SNTJmjf3rCYWiq7AcIKZO0a\nIezExTsXmbhjIrky5WJJ6yV4u3lbPQZNg23bDMV+0SLDevjt20OrVoYbsYTtkbVr7JC9zw03p7SU\nC5+sPowrlabtAAAbSElEQVQNHkvWDFkp820ZVp5a+cL3rTE3XClDf/6bbwzTMvv3N6yY6ednKPRr\n1hhuytKbzJM3DynyQliZazpXPqr8EW+Xe5t2S9rxyepPiI2L1ScWV2ja1LAm/vnzEBAAn38OBQvC\n0KFp46YrRyftGiF0dOnuJSbumIh7enf+bPMn+dzz6R0SYLg4+8MPhoefVKoE3boZ1tKR5ZH1Ie0a\nIexUXve8jA4ajWdmT1779jX+OvGX3iEBhmfYTptmePJVp04wYYKhnTNmDNyQRUvsihR5naSlPnRy\n0nouXNO58kGlD3i//Pu0n9Ce3qt669a+eVnGjIaLslu3wh9/wIkThidcvfceHDhg2WNLT948pMgL\nYSNqFahFt4rdWHV6FTV+qMGVqCt6h/SCChXgp5/g5EkoXNiwMmatWvDXX7ZxoVYkTHryQtiY6Jho\nZu2bxb4r+/j1rV8JKhikd0gJiokxnN2PHWt41OH//meYd+/qqndkjkd68kI4ENd0rvSo0oO2pdvS\nYmELRm4aiTVPxkyVLh20bGlYDnnqVMO8+8KFYdIkuH9f7+jEU1LkdZLW+9DPk1wYPZ+LeoXqMSJg\nBDP3zqTR/EbcfXw3iVfqRykIDoa1aw1n9lu2GIv9w4cp36/05M1DirwQNqyQRyHG1R3H3cd3Kf1N\naQ5csfDVzlSqWNGwMNrff0NoqKHYT58Ojx/rHVnaJT15IeyApmksObaEP47/wYT6E+hSsYveIZlk\n714YPNiwWNrw4YbpmE5yavnKZO0aIdKIg1cPMmnnJBoVbsT3zb7H1dk+rnJu3w6ffgqPHhnm3AcG\n6h2RfZELr3ZI+tBGkguj5HLxmudrjK83nv1X9lNhZgUu3L5gpchSp3p1w1z7/v0Nc+ybNzdMxUyK\n9OTNQ4q8EHbGI6MHwwOGUzBbQcrPKm+VNerNQSnDAmhHj0LNmlCjBnz2Gdy7p3dkjk3aNULYsXVn\n1/HTgZ8YVHsQn/l/pnc4r+TaNcNDTDZuNEzBbNZM74hsl/TkhUjDTkee5qvtX+Hv48+81+eRwSWD\n3iG9kvXroXt3KFkSpkwBHx+9I7I90pO3Q9KHNpJcGKUkF4VzFGZ83fGcvnmaCrMqcP72efMHZkFB\nQXDwoGFRtKcLo8XFSU/eXKTIC+EA3DO4M6zOMIrmKEqFmRVYc2aN3iG9kvTpDVMtt26FefOgQQO4\nfl3vqByDtGuEcDAbzm3g+/3fM6DWAPr590OpFP2Vr5uYGMN6OJMnw9dfG9bDsbO3YHbSkxdCvODM\nzTN8te0rquWrxvw35pPRJaPeIb2y/fuhY0coUQJmzEjbz5+Vnrwdkj60keTCyFy58PPwY1y9cZy/\nfZ6KsyrazXz65925E8qePZA/v6FXv3Wr3hHZJynyQjgo9/TuDK49mMIehan4XUU2hW3SO6RXliGD\n4Q7Zb76BN96AceNk7fpXJe0aIdKA1WdWM+ffOUyoP4GuFbvqHU6KXLgArVtDrlzw88+QNaveEVmP\ntGuEEElq4NeAgTUHMuCfAXy44kPiNPs7Hc6f33DjlI8PVK2a/LIIwkCKvE6kD20kuTCyZC5K5S7F\n2OCxrD6zmqA5QTa7Pv1TCc2Td3U1LF38ySeGpRFWrbJ+XPZGirwQaUjuLLkZEzSG6Nhoys4oy4kb\nJ/QOKUW6doXffzcsdjZlit7R2DbpyQuRBmmaxoLDC1h9ZjXz35hP4yKN9Q4pRc6dg8aNDR/jxjnu\nWvXSkxdCvBKlFO3LtOf98u/TZnEbxm8br3dIKVKwoGFq5e7d0KaNYb168SIp8jqRPrSR5MLI2rmo\nXaA2Q+sMZdy2cXT6vROxcbFWPX5STF27xsMD1qwx3BXbsCFERVk2LnsjRV6INK5IjiKMqzuOXZd3\nUXt2be48uqN3SK8sQwb45RcoXhzq1oWbN/WOyHZIT14IAcDjmMd8vfNrIqIiWN1hNUVyFNE7pFem\naYYHkaxeDWvXgpeX3hGZh/TkhRCplj5devrV6EfFPBWp+n1VQs+H6h3SK1PKcAG2ZUuoXRsuX9Y7\nIv1JkdeJ9KGNJBdGeudCKUWnsp3oUKYDIb+E8NP+n3SLJaXryStlWLb47bchOBiuXjVrWHYnnd4B\nCCFsTz2/enhm9uSTNZ9w/MZxxtQdY3dLFg8YANHRhh79hg2QM6feEelDevJCiERdvnuZEZtHUCVv\nFRa+uRDXdK56h/RKNA0GDjTcGbt+PWTPrndEKWPxnrxSqqFS6rhS6qRSql8S495SSsUppSqkJBgh\nhG3xdvfmq7pfcSryFNV/rM71+/b1uCalYNQoCAgwPCj84UO9I7K+ZIu8UsoJmAY0AEoBbZVSxRMY\nlwXoCewwd5COSO/eqy2RXBjZYi7c0rsxLGAY2dJno/zM8hy9ftQqxzXXM16VMixXnD8/tG1rePJU\nWmLKmXwV4JSmaWGapj0BFgLNExg3AhgLPDZjfEIIG+Di7EKfan2oU6AO/j/688/Zf/QO6ZU4OcFP\nP8H9+/DRR4Y2TlphSpHPC1x8bjs8/mvPKKXKAfk0TVtpxtgcWpmqZfQOwWZILoxsPRdtSrfh7bJv\n0+LXFsw+MNuixwoICDDr/lxdDYua7dkDI0aYddc2zZTZNQk1+5/9HlSGS+5fA52TeY0QwgEEFQwi\nZ8ac9Pm7D2G3wxgSMETvkEzm5gYrVhjWoy9eHFq10jsiyzOlyIcD+Z/bzgc8f4uBG4ZefWh8wfcC\n/lJKNdM0bd/LO5s6YCpe+Qy3oWV2y0yhEoWenb087Uemhe3ne6+2EI+e2y/nRO949Nw+e+wszd9u\nbjPxJLb9mtdrdM7amcm/Tubc7XP82PxHNm00PF7w6Rn40556SrcnTZpEuXLlzLa/57f/+gvq1Anl\n1i3o1s38+0/tdmhoKLNnzwbA19eX1Eh2CqVSyhk4AQQDEcAuoK2maccSGb8B+ETTtP0JfE+mUMY7\ntPOQzf9pbi2SCyN7y8XNhzcZsWkEvtl8Wdp2KZlcMplt36GhoWZv2Tzv99+hd2/YuRO8vS12GLNI\nzRRKk+bJK6UaApMx9PB/0DRtjFJqGLBb07TlL41dD/RN6CxeirwQjufhk4eM3TqW6Lho1nVch2cW\nT71DMtmoUfDXX7Bpk2GRM1tl8SJvLlLkhXBMsXGxfLPnG45eP8rajmspkauE3iGZRNMMDwf38IAZ\nM/SOJnGyQJkdssX50HqRXBjZay6cnZzpWaUndQrUocaPNcyyuJm55sknRSn4/nvDsgc//2zxw+lC\nirwQwmzalG5DxzIdCfklhAWHFugdjknc3WHxYsPDwQ/Z5+/YJEm7Rghhdvsi9jFxx0QG1R5E3xp9\n9Q7HJHPnGnr0u3cbCr8tkXaNEMKmVMhTgaF1hjJmyxg+Xf0p1jyZTKlOnQxr0PfooXck5iVFXif2\n2nu1BMmFkSPlorBHYUYHjWbhkYW0/739Kz8/1ho9+ZdNmmSYUvnrr1Y/tMVIkRdCWIyXmxdjg8ey\n+/JuGsxrwMMntr0MZObMMH8+9OwJFy8mP94eSE9eCGFxD588ZPSW0SilWNdpHR4ZPfQOKUlffml4\nRuw//xgWN9Ob9OSFEDYto0tGBtceTBbXLFT5rgpht8P0DilJ/fpBbKxhiWJ7J0VeJ47Ue00tyYWR\nI+cinXM6+lbvS8lcJan6fVUOXU36verRk3/K2dkw22bsWDh+XLcwzEKKvBDCapRSdKnQhfqF6lN7\ndm02hW3SO6RE+frC0KHw3nsQF6d3NCknPXkhhC7WnFnDnH/nMLvFbN4o8Ybe4SQoLg7q1DEsSdyz\np35xpKYnb8pSw0IIYXb1/eqTNX1W3v7zbSIfRNKlYhe9Q/oPJyfDsgf+/hASYji7tzfSrtGJI/de\nX5Xkwiit5aJqvqp8XvNzPlv7GeO2jnvhe3r25J9XrBj07QtdutjnYwOlyAshdFU6d2nD3bFbx9B/\nXX+9w0lQ375w4wb88ovekbw66ckLIWzC5buXGbJxCCFFQ5gVMgvDg+Zsx/bt8NZbcOyY9de2kXny\nQgi75+3uzZjgMaw+s5qWi1oSExujd0gvqF4dGjSAYcP0juTVSJHXSVrrvSZFcmGU1nORI1MOxgSP\n4d+r/1L1i6o8inmkd0gvGDPGsO784cN6R2I6KfJCCJvilt6NUYGjuPnoJgGzA4h6HKV3SM/kzg1D\nhhhWqrSXi7DSkxdC2KTo2GjGbh3Lo5hHbOi8gRyZcugdEmBY7qByZfjsM2jb1jrHlJ68EMLhuDq7\nMqDmADwyelDl+yqE3w3XOyTAsOTB5Mnw+efwyLa6SQmSIq+TtN57fZ7kwkhyYXRo5yGcnZz5pNon\nlMhZgqrfV+Vk5Em9wwKgVi0oVw6mTNE7kuRJkRdC2DSlFN0qdsPfx5+aP9bk8DXbuOo5diyMGweR\nkXpHkjTpyQsh7MaiI4tYcXoFqzusppJ3Jb3D4aOPwMXF8EQpS5KevBAiTWhZqiVvFH+DunPrsvXC\nVr3DYcgQmDcPTp/WO5LESZHXifRejSQXRpILo8Ry0axYM9qXaU/jBY1Zd3adlaN6Ue7c8Mknhouw\ntkqKvBDC7jQs3JD3yr/HG7++wfKTy3WNpU8f2LoV9u/XNYxESU9eCGG3tl3cxvTd0/mh2Q+0LNVS\ntzimTIF162DpUsvsX3ryQog0qYZPDT6u+jHvLX2POQfm6BZH166GM/ndu3ULIVFS5HUivVcjyYWR\n5MLI1FxUyluJfjX60XNVT77d862Fo0pYhgwwYAAMHqzL4ZMkRV4IYfde83qNL2p/wefrPmfi9om6\nxPDuu3D0KGzbpsvhEyU9eSGEwzgdeZoRm0fQt0ZfBtQaYPXjf/89LFxo6M+bk/TkhRACKJyjMMMC\nhjF+23iGhVp/4ffOneHsWdiyxeqHTpQUeZ1I79VIcmEkuTBKaS58s/kyInAEU3dN5Yv1X5g5qqS5\nuEC/foZ1522FFHkhhMPJnzU/IwJGMHPPTD5fZ907lTp3hr174eBBqx42UdKTF0I4rMtRlxm0YRCd\ny3ZmXP1xVjvu2LFw6JBhyQNzkJ68EEIkwNvNm5FBI5l7cC4f//2x1Y77wQewahWcO2e1QybKpCKv\nlGqolDqulDqplOqXwPc/VkodUUodUEqtVUr5mD9UxyK9VyPJhZHkwshcuciTJQ9fBn7JgsML6Lmy\np1n2mZysWQ03SE2YYJXDJSnZIq+UcgKmAQ2AUkBbpVTxl4btAypqmlYOWAJY7+8iIYRIRu4sufky\n6EsWH11M9+XdsUabundvWLAArl2z+KGSlGxPXilVDRiiaVqj+O3+gKZp2thExpcDpmqaViuB70lP\nXgihmxv3bzBww0CaFm3KzKYzUSpFbW6Tde8OOXLAyJGp24+le/J5gYvPbYfHfy0x7wGrUhKMEEJY\nUs7MORkVNIoVp1bw3tL3LH5G36cPzJql77Ng05kwJqHfHglmRinVAagI1ElsZ1MHTMUrnxcAmd0y\nU6hEIcpULQMYe3BpYfv5fqMtxKPn9ss50TsePbfPHjtL87eb20w8em7/Nfsvi9WH0UGj+fS7T2l0\nqBGrBq5CKUVoaCgAAQEBAGbbrlgxgAULoFAh018fGhrK7NmzAfD19SU1TG3XDNU0rWH8doLtGqVU\nXWAyUFvTtASfeijtGqNDOw89++FL6yQXRpILI0vn4tbDWwxYP4AGhRvwfcj3Fmvd/P039O9vWKUy\npYdITbvGlCLvDJwAgoEIYBfQVtO0Y8+NKQ8sAhpomnYmiX1JkRdC2IybD28yYP0AmhRpYrEefVwc\nlCwJM2dCnUR7HEmzaE9e07RYoAewBjgCLNQ07ZhSaphSqmn8sK+AzMAipdR+pdSfKQlGCCGsySOj\nB6MCR7H85HK6r+hukWM4OUHPnjB5skV2nyy541Un8me5keTCSHJhZM1c3Lh/gwEbBvBG8TeY3mS6\n2fd/7x4UKGBY7iAlLXa541UIIVIhZ+acjAocxZJjS+i1qpfZ958lC7z9Nkw3/++PZMmZvBBCxLt2\n7xoDNgygXel2fN3wa7Pu++xZqFIFLl6EjBlf7bVyJi+EEGaQO0tuRgWOYv6h+fRd09es+y5UCCpX\nhsWLzbrbZEmR14msUWIkuTCSXBjplQvPLJ6MChrFnH/n0H9tf7Puu2tXw81R1iRFXgghXuKVxYuR\ngSP54cAPDPxnoNn227QpnD5teBastUhPXgghEnEp6hJfrP+Cjyp/xLBA8zxOcOBAePAAvn6Flr/0\n5IUQwgLyuuVleMBwpu6aypgt5nmm3/vvw88/W289GynyOpHeq5HkwkhyYWQrufDJ6sOwOsMYs2UM\nU3ZMSfX+ChaEihVhyRIzBGcCKfJCCJEM3+y+DK49mEGhg5i1N/VXTrt2he++M0NgJpCevBBCmOjo\n9aOM2jyK6Y2n07FsxxTv58kT8PGBjRuhWLHkx0tPXgghrKBkrpL8r8b/+GjlRyw6sijF+3FxgY4d\nYc4cMwaXCCnyOrGVfqMtkFwYSS6MbDUXZb3K8nG1j3l36bssO7Esxfvp1MlwATYuzozBJUCKvBBC\nvKJK3pXoUbkH7X9vz9oza1O0jzJlIGdOiH9WiMVIT14IIVIo9Hwo3+/7nuXtllOrwH8ea52sr7+G\nf/+F+IdAJUp68kIIoYMA3wA6le1EyC8h7Arf9cqvb9sW/vwT7t+3QHDxpMjrxFb7jXqQXBhJLozs\nJRf1/erTqlQrGs5vyIErB17ptV5e4O8Pf/xhoeCQIi+EEKkWUjSEZsWaUe/nehy/fvyVXtupE8yd\na6HAkJ68EEKYzYJDC9hwfgPb3ttGoeyFTHrNw4eQNy8cOmT4b0KkJy+EEDagXZl2VM9XncA5gURE\nRZj0mowZ4c03Yd48y8QkRV4n9tJvtAbJhZHkwshec/FOuXcolqMYAXMCuPXwlkmv6dABfvnFMvFI\nkRdCCDNSStGjcg9yZcpF0Nwg7kcnP3WmZk24dg1OnLBAPNKTF0II84uJjWHE5hFkds3MP53+wdXZ\nNcnxvXoZbo4aPPi/35OevBBC2Jh0zukYWHMgNx7coMXCFsTGxSY5vnVr+O0388chRV4n9tpvtATJ\nhZHkwsgRcuGazpWhdYZy/MZxOvzegbgkFqqpXh3u3IEjR8wbgxR5IYSwoEwumRgeMJxNFzbx0cqP\nEh3n5AQtW5r/bF568kIIYQU3Htyg37p+vF/+fb6s+2WCY3buhM6d4dgxUM914KUnL4QQNi5nppwM\nDxjOt3u/ZfzW8QmOqVIFHj+GgwfNd1wp8jpxhH6juUgujCQXRo6Yi7zueRlSewjDNw3nh30//Of7\nSkGrVuZt2UiRF0IIK/Lz8KN/zf70Wd2HJUf/+zTvVq3g11/BXJ106ckLIYQOdl3axeSdk1nSagn1\n/Oo9+7qmgZ+fYQni114zfE168kIIYWeq5K3C++Xf561Fb7H70u5nX1cKWrQwFHlzkCKvE0fsN6aU\n5MJIcmGUFnIRWDCQ1iVb02h+I07cMK5pIEVeCCEcREixEIIKBhE8N5hLdy8BUKMGXLwIYWGp37/0\n5IUQwgZM2zWNsDth7HhvB9kyZuPdd6FcOcOaNtKTF0IIO/dR5Y/InjE79efV59GTR2Zr2ZhU5JVS\nDZVSx5VSJ5VS/RL4vqtSaqFS6pRSartSKn/qQ3NsaaHfaCrJhZHkwiit5UIpxWfVP+Phk4e0+LUF\nQcGx7N0LkZGp22+yRV4p5QRMAxoApYC2SqniLw17D7ipaVoRYBLwVerCcnxnj53VOwSbIbkwklwY\npcVcuDi7MKj2IE5GnuSD1Z0JDoYVK1K3T1PO5KsApzRNC9M07QmwEGj+0pjmwJz4zxcDwakLy/Hd\nj0r+QQJpheTCSHJhlFZzkdElI8PqDGP9ufVEFuvBwIkdUrU/U4p8XuDic9vh8V9LcIymabHAbaWU\nR6oiE0KINCpbxmz0KNKDzZu+JbzJ/FTtK50JYxK6ovvylJyXx6gExgDwMOYh0bHRJhzWsV2+eNmk\nx4KlBZILI8mFUVrPxcrFK9EC4iDpB0olK9kplEqpasBQTdMaxm/3BzRN08Y+N2ZV/JidSilnIELT\ntNwJ7Mt68zWFEMKBpHQKpSln8ruBwkqpAkAE0AZo+9KYZUBnYCfQElhvziCFEEKkTLJFXtO0WKVU\nD2ANhh7+D5qmHVNKDQN2a5q2HPgB+FkpdQqIxPCLQAghhM6seserEEII67LIHa9y85SRCbn4WCl1\nRCl1QCm1Vinlo0ec1pBcLp4b95ZSKk4pVcGa8VmTKblQSrWK/9k4pJSaZ+0YrcWEfyM+Sqn1Sql9\n8f9OGukRp6UppX5QSl1VSiX6XCil1JT4unlAKVXOpB1rmmbWDwy/OE4DBQAX4ABQ/KUx3YFv4j9v\nDSw0dxy28GFiLuoAGeI//yAt5yJ+XBZgI7ANqKB33Dr+XBQG9gLu8ds59Y5bx1zMBLrFf14COKd3\n3BbKRU2gHHAwke83AlbEf14V2GHKfi1xJi83TxklmwtN0zZqmvYofnMH/70HwVGY8nMBMAIYCzy2\nZnBWZkouugDTNU27C6Bp2g0rx2gtpuQiDnCP/zwbcMmK8VmNpmlbgFtJDGkOzI0fuxPIqpTyTG6/\nlijycvOUkSm5eN57wCqLRqSfZHMR/+dnPk3TVlozMB2Y8nNRFCimlNqilNqmlGpgteisy5RcDAM6\nKqUuAsuBnlaKzda8nKtLmHBSaMoUyldl1pun7JwpuTAMVKoDUBFD+8YRJZkLpZQCvsYwFTep1zgC\nU34u0mFo2dQG8gOblVKlnp7ZOxBTctEW+EnTtK/j79uZh2EdrbTG5HryPEucyYdj+KF8Kh9w+aUx\nFwEfgPibp9w1TUvqzxR7ZUouUErVBT4HQuL/ZHVEyeXCDcM/3FCl1DmgGvCXg158NeXnIhz4S9O0\nOE3TzgMngCLWCc+qTMnFe8BvAJqm7QAyKKVyWic8mxJOfN2Ml2A9eZklivyzm6eUUq4Y5swvfWnM\n05unIImbpxxAsrlQSpUHZgDNNE1L5aKiNi3JXGiadlfTtNyaphXSNK0ghusTIZqm7dMpXksy5d/I\nn0AQQHxBKwI44rKMpuQiDKgLoJQqAaR34GsUisT/gl0KdIJnKxHc1jTtanI7NHu7RpObp54xMRdf\nAZmBRfEtizBN01roF7VlmJiLF16Cg7ZrTMmFpmmrlVL1lVJHgBigryP+tWviz0Vf4Dul1McYLsJ2\nTnyP9ksptQAIAHIopS4AQzCsXKNpmjZL07SVSqnGSqnTwH3gHZP2Gz8dRwghhAOSx/8JIYQDkyIv\nhBAOTIq8EEI4MCnyQgjhwKTICyGEA5MiL4QQDkyKvBBCODAp8kII4cD+D25pxSSRQ0yZAAAAAElF\nTkSuQmCC\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7fe907ad0198>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# fix nodes\n",
        "h = 1\n",
        "x = nodes * h\n",
        "\n",
        "# find interpolant\n",
        "coeffs = la.solve(vander(x), f(x))\n",
        "\n",
        "# evaluate interpolant\n",
        "plot_x = np.linspace(0, h, 200)\n",
        "interpolant = vander(plot_x, len(coeffs)) @ coeffs\n",
        "\n",
        "# plot\n",
        "pt.plot(plot_x, f(plot_x), label=\"f\")\n",
        "pt.plot(plot_x, interpolant, label=\"Interpolant\")\n",
        "pt.fill_between(plot_x, 0*plot_x, interpolant, alpha=0.3, color=\"green\")\n",
        "pt.plot(x, f(x), \"og\")\n",
        "pt.grid()\n",
        "pt.legend(loc=\"best\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Compute the following:\n",
        "\n",
        "*   The true integral as `true_val` (from `int_f`)\n",
        "*   The quadrature result as `quad` (using `x` and `weights` and `h`)\n",
        "*   The error as `err` (the difference of the two)\n",
        "\n",
        "    (Do not be tempted to compute a relative error--that has one order lower.)\n",
        "\n",
        "Compare the error for $h=1,0.5,0.25$. What order of accuracy do you observe?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 134,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "1 0.6456349206349207 0.610677083333 0.0349578373016\n",
            "0.5 0.4054346478174603 0.40550104777 6.6399952722e-05\n",
            "0.25 0.22314353367638967 0.223148116221 4.58254464089e-06\n",
            "0.125 0.11778303564688754 0.117783224384 1.88736665879e-07\n",
            "0.0625 0.06062462181642996 0.0606246286393 6.8228956096e-09\n"
          ]
        }
      ],
      "source": [
        "\n",
        "errors = []\n",
        "\n",
        "for h in [1, 0.5, 0.25, 0.125, 0.125*0.5]:\n",
        "    true_val = int_f(h)\n",
        "    quad = h * weights @ f(h * nodes)\n",
        "    error = abs(quad - true_val)\n",
        "    print(h, true_val, quad, error)\n",
        "    \n",
        "    errors.append(error)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Estimate the order of accuracy:\n",
        "\n",
        "We assume that the error depends on the mesh spacings $h$ as\n",
        "$E(h)\\approx C h^p$ for some unknown power $p$.  Taking the $\\log$\n",
        "of this approximate equality reveals a linear function in $p$:\n",
        "$$\n",
        "E(h) \\approx C h^p \\quad \\iff \\quad \\log E(h) \\approx \\log(C) +\n",
        "p\\log(h).\n",
        "$$\n",
        "You can now either do a least-squares fit for $\\log C$ and $p$ from\n",
        "a few data points $(h,E(h))$ (more accurate, more robust), or you\n",
        "can use just two grid sizes $h_1$ and $h_2$, and estimate the slope:\n",
        "(less accurate, less robust)\n",
        "$$\n",
        "   p \\approx \\frac{ \\log(\\frac{E(h_2)}{E(h_1)}) } {\\log(\\frac{h_2}{h_1})}.\n",
        "$$\n",
        "This is called the *empirical order of convergence* or EOC.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 135,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "9.04021800427\n",
            "3.85696137438\n",
            "4.60170230716\n",
            "4.78984676915\n"
          ]
        }
      ],
      "source": [
        "for i in range(len(errors)-1):\n",
        "    print(np.log(errors[i+1]/errors[i])/np.log(1/2))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": []
    }
  ],
  "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.5.1+"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}