{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Finding the Weights in Gaussian Quadrature"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import numpy.linalg as la\n",
        "import scipy.special as sps\n",
        "import matplotlib.pyplot as pt"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here's a utility routine to do biseciton, to use below:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "def bisection(f, a, b, tol=1e-14):\n",
        "    assert np.sign(f(a)) != np.sign(f(b))\n",
        "    while b-a > tol:\n",
        "        m = a + (b-a)/2\n",
        "        fm = f(m)\n",
        "        if np.sign(f(a)) != np.sign(fm):\n",
        "            b = m\n",
        "        else:\n",
        "            a = m\n",
        "            \n",
        "    return m"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Set the number of nodes:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n = 5"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Node selection: Plain Gauss"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Gauss nodes are the roots of the $n$th Legendre polynomial $P_n$:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.5/dist-packages/IPython/core/formatters.py:92: DeprecationWarning: DisplayFormatter._ipython_display_formatter_default is deprecated: use @default decorator instead.\n",
            "  def _ipython_display_formatter_default(self):\n",
            "/usr/local/lib/python3.5/dist-packages/IPython/core/formatters.py:669: DeprecationWarning: PlainTextFormatter._singleton_printers_default is deprecated: use @default decorator instead.\n",
            "  def _singleton_printers_default(self):\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7f4a1c40d438>]"
            ]
          },
          "execution_count": 5,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEACAYAAABVtcpZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XeYFeXZx/HvDQjYADvSEY0tNux9jcFusMUSUbFE7Mb2\n2oKiRCOoRGONDcFG7FETI1jQIIoKK6ggLLLAAgIqAgooyD7vH/eiiLuwu2fOmZkzv8917cXZ3dmZ\n+wxz5p6nWwgBERHJpgZxByAiIvFREhARyTAlARGRDFMSEBHJMCUBEZEMUxIQEcmwSJKAmT1oZrPM\nbMxKtvm7mZWZ2Ydmtn0UxxURkdxEVRLoDxxY0y/N7GCgUwhhM6AHcG9ExxURkRxEkgRCCMOAr1ey\nSVdgYNW2I4DmZrZRFMcWEZH6K1SbQGugYrnvp1f9TEREYlSoJGDV/EzzVYiIxKxRgY4zDWi73Pdt\ngBkrbmRmSgwiIvUQQqjuYXuVoiwJGNU/8QO8AJwMYGa7AXNDCLOq2zCEoK+Ivq699trYYyimL51P\nnc8kft18c27PzpGUBMzscaAEWM/MpgLXAo2BEEK4L4TwHzM7xMwmAguAU6M4rohI1k2fntvfR5IE\nQgh/qMU250VxLBER+UmuSUAjhotYSUlJ3CEUFZ3PaOl8RiPXJGAhJKct1sxCkuIREUm6Dh1gyhQj\n1LNhWElARCSlKith9dVh8eL6JwFVB4mIpNSXX8Jaa+W2DyUBEZGUmjEDWuc494KSgIhISk2friQg\nIpJZSgIiIhmmJCAikmFKAiIiGTZ9OrRqlds+lARERFJKJQERkQxTEhARyajvvoNvv4X1189tP0oC\nIiIpNGMGbLwxNMjxLq4kICKSQlE0CoOSgIhIKkXRHgBKAiIiqaQkICKSYdOmQZs2ue9HSUBEJIWm\nToV27XLfj5KAiEgKVVQoCYiIZNbUqdC2be770fKSIiIp8/330KwZLFwIDRuCmZaXFBHJjGnTfIxA\nw4a570tJQEQkZaKqCgIlARGR1ImqURiUBEREUkclARGRDFNJQEQkw6IaKAZKAiIiqaPqIBGRjApB\nJQERkcyaNw/MoHnzaPanJCAikiIVFV4VZPUaH/xLSgIiIikSZVUQKAmIiKTKspJAVJQERERSRCUB\nEZEMUxIQEcmwKMcIgJKAiEiqlJdDx47R7U+LyoiIpMSyxWQWLIBGjX76uRaVERHJgIoKaN365wkg\nV0oCIiIpUV4OHTpEu08lARGRlIi6PQCUBEREUkNJQEQkwyZPVnWQiEhmqSQgIpJhSgIiIhm1YAHM\nnw8tW0a730iSgJkdZGafmtkEM7u8mt+fYmazzWxU1ddpURxXRCQrpkzxOYMaRPzonvOQAzNrANwJ\n7A/MAN43s3+FED5dYdNBIYQLcj2eiEgW5aMqCKIpCewClIUQpoQQlgCDgK7VbBfROjgiItmT5CTQ\nGqhY7vtpVT9b0VFm9qGZPWlmbSI4rohIZuQrCUQxA0V1T/grzgL3AvB4CGGJmfUABuDVR7/Qq1ev\nH1+XlJRQUlISQYgiIuk2eTLsuqu/Hjp0KEOHDo1kvznPImpmuwG9QggHVX1/BRBCCH1q2L4BMCeE\n0KKa32kWURGRanTuDP/4B+y88y9/F/csou8Dm5pZezNrDByPP/kvH+DynZq6AmMjOK6ISGYktjoo\nhLDUzM4DBuNJ5cEQwjgzuw54P4TwEnCBmf0OWALMAbrnelwRkayYMwcqK2G99aLftxaVERFJuBEj\n4JxzYOTI6n8fd3WQiIjk0cSJsNlm+dm3koCISMKVlcGmm+Zn30oCIiIJp5JAzBYu9Hk7Pv/cG2dE\nRApJJYEYfPUV3HADbLcdrL8+7LWXv27eHA46CAYMgEWL4o5SRLJAJYECCgHuuQc23xwmTYK774Z5\n86CiAmbPhqlT4fTTYdAg+NWv4Ikn/G9ERPJhzhxYsgQ22CA/+1cX0eV89x2ccgp89hkMHAhbbbXy\n7YcNgwsugGbN4LHHoHV1MyaJiOTgvffg7LNr7h4K6iIaiUWL4NBD/fWwYatOAOBVRO+/D126+FDu\niKbyEBH50cSJ+WsPACUBwBt7TzgBNtwQHn8cmjat/d82bAhXXw0PPwzHHQePPJK3MEUkg8rK8tce\nANHMIpp6vXp5Q/Brr/lNvT4OOADeeMP/XbgQevSINEQRyaiJE2H/audcjkbmk8Cbb8IDD0BpKTRu\nnNu+ttrK9/fb33pj8VlnRROjiGRXWVl+7yWZTgILFkD37nD//bDRRtHss1MnL1Hsvbe35h99dDT7\nFZFsynebQKaTQO/e3ri7rEE4KptsAi+9BAce6GMM9t032v2LVCcEmDDBqyVLS+HTT32Q47x5XkXZ\ntCmstRa0b+91zDvu6NfmNttEv3i5ROPrr2HxYm+vzJfMdhGdMAH23BM+/ji6UsCKXn0VunWDd97J\nzzzgIuBdB594Ap58Esxgv/28t9oWW/h116IFrLmm94D75htfoWr8eO96+OabMH++d2o46SRPDJIc\n773nVUGjRq18u1y6iGY2CRxzjH9QLr88v8e5/XZ46CEYPtw/iCJR+OEHePppuO02mDXLHzaOP97b\npayOt4Lx471X3IAB0K4dXHYZHHZY3fcj0Rs4EP77X///WRklgToaORK6dvUGl9VXz++xQoBTT/WB\naE88oQ+W5KayEp56Cq65xkuwl1ziN+z69mpb3rLEctNNsMYacOutsPvuue9X6u/KK/3/omfPlW+n\nwWJ11Levf3jynQDAb/r33uv1sw88kP/jSfEaOdJvyrfcAnfd5VU5XbtGkwAAGjXy0sSoUd7F+fe/\n9weYr7+OZv9Sd+PGwZZb5vcYmUsCkyZ5750zzijcMZs29bmGrroKxmp1ZamjhQvhwgvhkEO8fnjE\nCO+GnK9SZYMGPn3KuHH+FLrNNl4lIYWnJJAH/frBH/8Ia69d2ONusYUXs48/3quGRGqjtBR22skn\nLxw71p/MC9WTZ+21vcTxyCP+mbn6aq8yksJYvNh7d+WzeyhkrE3gyy+9a9zYsbDxxnk7TI1C8CSw\n4YZwxx2FP76kR2Wl18n37euNvyeeGG88s2fDH/7g1/DTT8M668QbTxaMHQtHHOE9GVdFbQK1dP/9\ncOSR8SQA8OL7P/7hYwj+8594YpDk+/Zb77323HM+QWHcCQD8weWVV2D77b1dYtKkuCMqfoWoCoIM\nJYEQoH9/OPPMeONo0QIefNAb3ubOjTcWSZ5Jk/wmu+66PuirQ4e4I/pJw4ZeOrngAh9kWVoad0TF\nTUkgYsOHe13qrrvGHQn85jdw+OFw8cVxRyJJ8tZbsMce/qBy//3QpEncEVXvnHPgzjt9hb133407\nmuI1bpy3JeZbZpJA//7eqJaUfvp9+viTnqqFBOD5570K6NFH4fzzk3Od1uSoo/wzdfjhvv6GRK9Q\nJYFMNAwvWABt2sAnn0CrVpHvvt5ef9274n30kVcTSTY9+KAPBnrxxfRN2zB4sI9WHjzY2wskGpWV\n3jtrxgxf13xV1DC8Cs8+6/WsSUoA4NVChx0Gl14adyQSlz594C9/8YFfaUsA4Otn3H23j2EoK4s7\nmuJRUeEPhrVJALnKxCyijz3mVUFJdNNNsPXW8L//+fTTkh1/+Ytfm2+/nbwHlLo45hjv5HDAAX4d\nt2kTd0TpV6j2AMhASWDuXG8Ujnq66Kg0b+6TzPXo4YNDJBuWJYDXX093AljmjDO8wfjgg31WUslN\nodoDIANJ4N//hpISn0c9qY46yhejufnmuCORQlg+AcQ1ZiUfLr3UezedeCIsXRp3NOn20Uew7baF\nOVbRJ4HnnvMBYklm5l3u/vY3X0VIile/ft4DqNgSAPx0HX/7rc+TJfU3ZkzhkkBR9w5atAhatoTP\nPvMVvpLu1lt9oq7Bg5PfRVDqbsAAnwJ62DBo2zbuaPLnq69gl13g2mvh5JPjjiZ9fvgBmjXzqTpq\nW4Oh3kE1GDIEOndORwIAnynyiy9WvYCEpM8LL8AVV/jUC8WcAADWW8+7u156KXzwQdzRpM/Eid5O\nVKgq7KJOAmmoClpeo0Zw333+4ZkzJ+5oJCpvvgmnn+6JoFA9PuK21VbedfS44zQ9Sl0VsioIijgJ\nVFb6RG1du8YdSd3ssos3FKtOtTh8+KEvzvLEE76caZYcc4z3yjvtNJ+7S2pHSSAipaVeDdS+fdyR\n1N0NN/hT44gRcUciuZg61QcD3nWXLwKTRTffDNOmeTdoqR0lgYgMHgxdusQdRf20aOEfnrPO0iIe\naTV/vieAiy7ykkBWNWkC//wn3HijHmpqS0kgIkOGpDcJgC/gse66/hQp6fLDD3DssT7dsmaKhY4d\nfZ3tbt28+6jUbN48X/xqk00Kd8yi7CK6cKEvgvH554VfRjJKn37qU0mMHl0co0qzIAQ4+2xfFvDF\nF72xX9ypp8Jqq3nnB6nesGHeMaSuU3Sri+gK3nrLu4amOQGA9yTp0UNPk2nSrx+8845XgSgB/Nzt\nt8Orr3p7l1Sv0FVBUKRJIO1VQcu7+mp47z1/T5JsL77oo75feskH+8jPNWsGAwf6g82sWXFHk0xK\nAhFJc6PwilZf3Yfin3MOfPdd3NFITcaN87EAzzxT/IPBcrHXXl4tdOaZ6jZandJS2G67wh6z6JLA\nzJneJW2nneKOJDqHHALbbONzz0vyzJ3r41H69k3G8qVJ16uXr6U8aFDckSTL4sXw8cewww6FPW7R\nJYFhw/xpo9jqY2+/He64QxPMJc3SpXDCCT6FcvfucUeTDo0bw0MPeffZ2bPjjiY5Pv7Ye1IVesbj\nokwCe+4ZdxTRa9vW554591wVo5Pk6qv9Ce6WW+KOJF123tknl7vggrgjSY4PPoinBqMok8Bee8Ud\nRX5ceKGvOfrUU3FHIuBTQTz5pPcEWm21uKNJn+uug1Gj4Pnn444kGeJKAkU1TuDbb2GjjXwq26ZN\nIwwsQd5+2yflGjtWPVDiVFrqyym+9lrhe3MUk7fe8uq0jz+GddaJO5p4de7sk+7ttlvd/1bjBKqM\nGOGNKsWaAMCrug480Oell3jMnu2z0959txJArvbZx8/lJZfEHUm8vvvOB4cWumcQRJQEzOwgM/vU\nzCaY2eXV/L6xmQ0yszIze8fM2kVx3BUVc1XQ8vr08aqIDz+MO5LsWbLE5wLq1i3bcwJF6a9/9ZXW\nXn017kjiM2YMbL65dwkvtJyTgJk1AO4EDgS2Bk4wsxVnTT8dmBNC2Ay4Deib63GXVz65nG4XdOO2\n5/fj3QndKJ9cHuXuE2f99X1CrlNPK+fE87uxX/f96HZB8b/vuCy7vvbrvh9b79ONho3Kuf76uKMq\nHmuv7XNknfHHck44N3vXc/nkcs69phuzmsTzvnNuEzCz3YBrQwgHV31/BRBCCH2W2+a/VduMMLOG\nwMwQwgbV7KvObQLlk8vpcl4XPtvuM2gMLIZOozsx5M4hdOzQMZe3lmifTSpn66O68P2h2XrfhVbd\n9dWxtBOv3a3zHKXyyeX8+pguLDwwW9dzVPevuNsEWgMVy30/repn1W4TQlgKzDWzdSM4Nj379fzp\nBAI0hs+2+4ye/XpGsfvEuva2nj8lAMjM+y606q6v8h10nqPWs1/PnxIAZOZ6TsL9K4ohVdVlnxUf\n51fcxqrZBoBevXr9+LqkpISSkpKVHnz6/Omw3go/bAwz5s9Y6d+lXVbfd6HpPBdGVs9zfd/30KFD\nGTp0aCQxRJEEpgHLN/S2AVZ8BxVAW2BGVXVQsxDC19XtbPkkUButm7WGxfyUSQEWQ6tmxT33clbf\nd6G1XFPnuRCyej3X932v+IB83XXX1TuGKKqD3gc2NbP2ZtYYOB5YcbLYF4FTql7/Hng9guMC0Pvi\n3nQa3clPJPxYp9b74t5RHSKRqnvfm3xY/O+7kEKA72f2Zs1Xsnd9FVp113PH0uI/z70v7s06b8R7\nfUUyWMzMDgJux5PKgyGEm8zsOuD9EMJLZtYEeATYAfgKOD6EMLma/dRrsNjHn5Sz/W96stdBM2jT\nvBW9L+5d1I1Jy5RPLqdnv55Mnz+D8SNb0e2Q3vTtU/zvu1Buuw0efhgee7ycv97bkxnzZ9CqWXau\nr0Jbdj3PmD+DinGt2H3z3gwcWPznefc9ymnaqie2Vv2vr1wahotixPD//uer8WR5DdNJk2CXXWDk\nSGjfPu5o0m/wYDjlFF8gpkOHuKPJnjlzYOutfQGanXeOO5r8WbLEl5GtqPC1xesr7t5BsYtrzo0k\n2WQT+NOfNCFXFMrK4KSTfE4gJYB4rLuuT8191lk+U2uxKi31z24uCSBXSgJF5LLLYPx4Ld+Xi/nz\nfW2A66/3KQ0kPt26+UCye++NO5L8ScKsx0oCRaRJE7jnHp9ueu7cuKNJn6VL4cQToaTEl0CUeJn5\nSOJevYp3OcokTHWT+jaBefOgdWu/6RXbQjL1de65sGCBN2pK7V11FQwf7us5a2ro5Pi///MVAwcO\njDuSaIUALVvC++9DuxxnU8t0m8CoUbD99koAy+vTx58wVC1Ue4MG+aR8Tz2lBJA011wDQ4fCm2/G\nHUm0ysp8lbVcE0CuUp8E3n9fVUErWmst6N/fG9W+/DLuaJJv1Cg4/3xf3GSDX8xoJXFbay3vrnvO\nOd6bpli8+irsv3/cURRBEhg1yhdjkJ/be29frOPcc+OOJNlmzvT57O+9N5653KV2jjzSn5hvuy3u\nSKIzeDB06RJ3FEXQJrDVVl6M1wf4lxYt8gTZq5evRiY/t3Ah7Lsv/O530LO45ykrChMn+qpbpaW+\n5naa/fCDTwk/fryvhpirzA4W++47X5Ju7lzvGSO/9N57cPjh/sFpVdzTsNTJ0qW+KMzaa3sDutXr\n4yOF1qsXfPQRPPNM3JHkZvhwOPtsGD06mv1ltmF47FjYdFMlgJXZZRevSz3ppOIedFNXl18OX38N\n99+vBJAml1/uN86XX447ktwMGeJrVCdBqpPAmDFa47U2/vxnTwA33RR3JMlwzz3w0kv+NNm48aq3\nl+RYfXW44w5vyP/uu7ijqb8hQ5LRHgBKApnQsCE8+qh/eIYNizuaeL38so8G/ve/fWoCSZ+DD/Y2\nwD59Vr1tEs2f76WZvfeOOxKnJJARbdrAAw/4iNg5c+KOJh7vvgsnn+wlgE6d4o5GcnHbbf5QM3Fi\n3JHU3RtvwK67xrOofHVSnQQ++khJoC4OOwyOOspvhJWVcUdTWJ98AkccAQMGwB57xB2N5KptWx9J\nfP75PvI2TV58EQ45JO4ofpLaJDBrlnezUo+XuunbF775xntZZMXkyXDQQXDrrcn68Elu/vQnmDoV\nnnsu7khqb+lSH8l/5JFxR/KT1CaBZVVB6tlRN6ut5lMjPPwwPPts3NHk3+zZ3gvjssu8KkyKR+PG\ncPfdngy+/TbuaGrn7bf9wbVjgtbKSX0SkLrbcENPAD16eDVJsZo3zxsRjz9e6ywUq3339a/rr487\nktp57rlklQJASSCzdtrJq0e6doUvvog7mujNm+clgD32gBzW4JYUuPlmnysr6Q80ISgJRGrMGNhm\nm7ijSLeTT/an5EMP9amni8WyBLDzzvD3v6vKsNi1bAnXXuuDIpPcSPzhh95dO2n3rVQmgaVLfc6N\nrbaKO5L0693b13I99lhvaE+7efPgwAM9AdxxhxJAVpx9ts8F9eCDcUdSs2WlgKRdk6lMAlOm+ORL\na60VdyTpZwb33edPUD16JPtJalVmz/apeZUAsqdhQ08AV14J06fHHc0vhQCPP+4PW0mTyiQwbhxs\nuWXcURSP1VaDJ5/0cRcXXZTORDB5si/Td+ihqgLKqm239SqhJFYLDRvmc5ztvHPckfySkoAAXqoa\nPBjeeSd9A3A++siH4F9wgTcCKwFk11VX+Sjif/4z7kh+rn9/OPXUZF6bSgLyoxYtPBGMHOlPU2kY\nVfzyy14F1LcvnHde3NFI3Jo0gYce8rEDSen19u233h7QrVvckVRPSUB+pnlzeOUVf7o+7TRYvDju\niKoXgt/4zzjDl4U84YS4I5Kk2HVXnzr97LOTUaJ95hmvqmzZMu5Iqpe6JBACfPqpkkA+NWvmiWDO\nHB9sNXdu3BH93MKF/lT15JM+KZzmApIV9e4NEybAwIFxR+JVQd27xx1FzVKXBGbP9nq19dePO5Li\ntuaaXoT99a9h99299JUEI0f6kpmNGsFbb6V/mUHJj6ZN4bHH4NJLYdKk+OIoLYWyMl/dL6lSlwSW\nVQUlsYGl2DRsCLff7vPu7LOPd3GLy7JFcQ4+2Ce/GzAA1lgjvngk+bbZxhuKTzopvjEwt9zi7RNJ\nXrwodWsM33MPjBrlywJK4Ywe7Wvy7rSTd8EsZEls9GhvqG7UCB55BNq1K9yxJd0qK33w4D77QM+e\nhT32lCleap00ydva8ilTawyrUTge223nw95btvT+2I88kv/eQ198ARde6Mvwde/ui3EoAUhdNGjg\nM+befbdfP4X0t7/B6afnPwHkSklAam2NNaBfP5+B9K67fBH7IUOi74HxxRdwzTWwxRZeDfTJJ/DH\nP/oHWqSuWrf2h5YTTyzcaOKvvvJG6QsvLMzxcpG6j5V6BsVvt918UNmll/pF3rmzf8hymYQuBO/p\nc8YZ8KtfwYwZ8MEHcOedsMEG0cUu2fTb33qV4rHHwvff5/94N9zgkzO2bp3/Y+UqVW0C33zj1RHf\nfKOnwqSorIT//MeL28OHe8PtAQdASQl06LDyBvyvvvJFNt54w/v6N23qXT/PPFM3foleZSUcfbRX\nz/Tvn7/OJeXl3nb2ySeFGxuQS5tAqpLAqFE+9Hr06AIGJbX2xRferfSNN+DNN71ksNlmvpLSmmt6\nw+78+T7+YMIEWLTIB/bst58v+7jddur1Jfm1YIE3Eh99tPcciloI8LvfeWn56quj339NckkCjaIO\nJp8mToRNN407CqnJBhv4U/yZZ/r3c+Z4H+mZM/3D98MPPhBtnXX8/7FVK930pbDWXNMXet9tN2jT\nxtfUiNKzz/o1//TT0e43n1KVBMrK/MlS0mHddf1JXyRJWrXyEfH77+9zDR13XDT7nTkTzj3XE0GT\nJtHssxBSlQQmToQ994w7ChFJuy23hP/+19uvGjTwMTC5WLLE568688z0TWOSquZVlQREJCrbbuuJ\n4KKLvE9/fZtHQ/AZbFdf3Ze5TJtUJQG1CYhIlLbf3nu1PfQQnHWWT05YF5WVcMkl3mll0CCfaiVt\nUpME5s/3rqGtWsUdiYgUk3btfOWvBQs8KbzzTu3+bv58rwJ6910vUTRrlt848yU1SeCzz6BTJ/Um\nEZHoNW8Ojz4KN94IRx0FRxwBI0ZUv+2SJT5D6bbb+kJMr70G661X2HijlJqG4bIyVQWJSH4dc4yP\nWXnoIe811LSpLwjTrp1PYTJhArz+uk9pMmAA7Ltv3BHnLjWDxW68EebNgz59ChyUiGRSZaWvB/De\ne/D5596LqH17v/Fvsknc0f1cJgaLlZWlr+uViKRXgwaw447+VcxS0yYwcaK6h4qIRC01SUBtAiIi\n0UtFEvjmG++Ope6hIiLRyikJmNk6ZjbYzMab2StmVu0aOma21MxGmVmpmT1f1+NMnOjdQzV9tIhI\ntHK9rV4BvBpC2Bx4Hbiyhu0WhBA6hxB2CCEcUdeDaKSwiEh+5JoEugIDql4PAGq6wec0xKu8PHld\nskREikGuSWDDEMIsgBDCTKCm9aCamNl7ZjbczLrW9SDl5dCxYy5hiohIdVY5TsDMhgAbLf8jIAB/\nrsNx2oUQZppZR+B1MxsTQiivbsNevXr9+LqkpISSkhImT4ZDD63D0UREitjQoUMZOnRoJPvKacSw\nmY0DSkIIs8ysJfBGCGGly8CbWX/gxRDCs9X8rtoRw1tsAc88A1tvXe9QRUSKVi4jhnOtDnoB6F71\n+hTgXytuYGYtzKxx1ev1gT2AsbU9QGUlTJnii5aLiEi0ck0CfYAuZjYe+C1wE4CZ7Whm91VtsyXw\ngZmVAq8Bfw0hfFrbA8yc6VO0rrlmjpGKiMgvJH4CueHDfeWfmqZ1FRHJujirg/JOPYNERPJHSUBE\nJMMSnwQmT1YSEBHJl8QngfJy9QwSEcmXVCQBlQRERPIj0b2DfvjBu4bOnw9NmsQYmIhIghVt76Dp\n02HDDZUARETyJdFJQO0BIiL5lfgkoPYAEZH8URIQEcmwRCeBKVOgffu4oxARKV6JTgIVFdCuXdxR\niIgUr0QngalTlQRERPIpseMEKithjTVgzhz/V0REqleU4wS++ALWXlsJQEQknxKbBKZOhbZt445C\nRKS4JTYJqFFYRCT/EpsE1CgsIpJ/iU4Cqg4SEcmvxCYBVQeJiORfYpOASgIiIvmX2CSgkoCISP4l\ncrDY99/7GIFFi6Bhw7ijEhFJtqIbLDZ9Omy8sRKAiEi+JTIJqCpIRKQwEpkE1CgsIlIYiUwCKgmI\niBRGIpOASgIiIoWR2CSgkoCISP4lMglUVKgkICJSCIlMAtOmQZs2cUchIlL8EpcEFi70QWLrrRd3\nJCIixS9xSWD6dGjdGqxeY99ERKQuEpcEVBUkIlI4iUsCy0oCIiKSf4lLAioJiIgUTuKSgEoCIiKF\nk8gkoJKAiEhhJC4JTJumkoCISKEkLgmoOkhEpHASt7LYaqsFFi6ERo3ijkZEJB2KamWx9dZTAhAR\nKZTEJYFWreKOQEQkO5QEREQyTElARCTDEpcE1DNIRKRwckoCZnaMmX1sZkvNrPNKtjvIzD41swlm\ndvnK9qmSgIhI4eRaEvgIOBJ4s6YNzKwBcCdwILA1cIKZbVHT9koC0Rk6dGjcIRQVnc9o6XwmQ05J\nIIQwPoRQBqysf+ouQFkIYUoIYQkwCOha08ZKAtHRhyxaOp/R0vlMhkK0CbQGKpb7flrVz6qlJCAi\nUjirHJZlZkOAjZb/ERCAq0MIL9biGNWVEmocprz++rXYo4iIRCKSaSPM7A3gkhDCqGp+txvQK4Rw\nUNX3VwAhhNCnmm2TM4eFiEiK1HfaiCgnaKgpgPeBTc2sPfA5cDxwQnUb1vdNiIhI/eTaRfQIM6sA\ndgNeMrPRGkiQAAADFklEQVSXq36+sZm9BBBCWAqcBwwGPgEGhRDG5Ra2iIhEIVGziIqISGHFOmI4\nH4PNsszM1jGzwWY23sxeMbPmNWy31MxGmVmpmT1f6DiTblXXm5k1NrNBZlZmZu+YWbs44kyDWpzL\nU8xsdtX1OMrMTosjzjQwswfNbJaZjVnJNn+vui4/NLPta7PfuKeNiHywWcZdAbwaQtgceB24sobt\nFoQQOocQdgghHFG48JKvltfb6cCcEMJmwG1A38JGmQ51+OwOqroeO4cQHipokOnSHz+X1TKzg4FO\nVddlD+De2uw01iSQj8FmGdcVGFD1egBQ0w1eDfA1q831tvx5fhrYv4DxpUltP7u6HmshhDAM+Hol\nm3QFBlZtOwJobmYbrWR7IP6SQG3UabBZxm0YQpgFEEKYCWxQw3ZNzOw9MxtuZkqoP1eb6+3Hbao6\nPsw1s3ULE16q1Paze1RV9cWTZtamMKEVpRXP93Rqca/M+xpehR5sVuxWcj7/XIfdtAshzDSzjsDr\nZjYmhFAeZZwpVpvrbcVtrJptpHbn8gXg8RDCEjPrgZewVLKqn3rdK/OeBEIIXXLcxTRg+Ya3NsCM\nHPeZWis7n1WNRhuFEGaZWUtgdg37mFn1b7mZDQV2AJQEXG2utwqgLTDDzBoCzUIIKyumZ9Uqz+UK\n5+1+4BeDSKXWpuHX5TK1ulcmqTpolYPNzKwxPtjshcKFlSovAN2rXp8C/GvFDcysRdV5xMzWB/YA\nxhYqwBSozfX2In5+AX6PN8LLL63yXFY9rCzTFV2Lq2LUfK98ATgZfpypYe6y6uGVCiHE9oU3XFYA\ni/DRxC9X/Xxj4KXltjsIGA+UAVfEGXOSv4B1gVerztUQoEXVz3cE7qt6vTswBigFRgPd4447aV/V\nXW/AdcBhVa+bAE9W/f5doEPcMSf1qxbn8kbg46rr8TXgV3HHnNQv4HH8yf57YCpwKt4L6MzltrkT\nmFj12e5cm/1qsJiISIYlqTpIREQKTElARCTDlARERDJMSUBEJMOUBEREMkxJQEQkw5QEREQyTElA\nRCTD/h/3DoSYWzrvZwAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7f4a1c159b00>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "nodes = sps.legendre(n).weights[:, 0]\n",
        "\n",
        "mesh = np.linspace(-1, 1, 300)\n",
        "\n",
        "pt.plot(mesh, sps.eval_legendre(n, mesh))\n",
        "pt.plot(nodes, 0*nodes, \"o\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Node Selection: Gauss-Lobatto"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Gauss-Lobatto nodes are (except for the endpoints) the roots of $P_{n-1}'$:\n",
        "\n",
        "(See [here](https://en.wikipedia.org/wiki/Legendre_polynomials#Additional_properties_of_Legendre_polynomials) or [here](http://dlmf.nist.gov/18.9.17) for a formula for $P_n'$.)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 27,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.5/dist-packages/ipykernel/__main__.py:5: RuntimeWarning: invalid value encountered in true_divide\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7f4a13d9ad30>]"
            ]
          },
          "execution_count": 27,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHP9JREFUeJzt3Xl8VOW9x/HPIwWqoICtgJQtLAIvyiIuF4WroUpxoRf3\nBREQb6ULorW3RauCFuteKBRtRQX1sqlQLahUQIgIgstFBEUQ2ZeCtYggSELIc//4BQmQwCSznGW+\n79drXkkmZ875zZkzv3nmWZ33HhERyR7HBB2AiIhklhK/iEiWUeIXEckySvwiIllGiV9EJMso8YuI\nZJmUJH7n3NPOua3OuSUl7qvlnJvhnFvhnHvdOVcjFccSEZHkpKrEPxbodsh9twOzvPctgNnAHSk6\nloiIJMGlagCXc64RMM1737b47+XAud77rc65ukCe975lSg4mIiIVls46/tre+60A3vstwElpPJaI\niCRIjbsiIlnmO2nc91bnXJ0SVT2fl7aRc06TBYmIVID33lXkcaks8bvi235Tgb7Fv/cB/l7WA733\nuqXoNmTIkMBjiNNN51PnMqy3ZKSkxO+cmwDkAt9zzq0HhgAPAi865/oB64ErU3EsEZFsNm4cNGuW\n3D5Skvi99z3L+Nf5qdi/iIiA9zB0KPzv/ya3HzXuxkxubm7QIcSKzmfq6Fwmb/Fi2LsXzjgjuf2k\nrB9/hQNwzgcdg4hIFAwaBJUqwf33g3MOX8HGXSV+EZEI8B4aN4Zp06Bt2+QSv6p6REQiYOFCqFYN\n2rRJfl9K/CIiETBpElxzDbgKlfEPpqoeEZGQKyyEBg0gLw9atLD7VNUjIhJjb7xhiX9/0k+WEr+I\nSMiNGwe9eqVuf6rqEREJsa+/hvr14dNPoXbtA/erqkdEJKZefhk6dTo46SdLiV9EJMTGjYPrr0/t\nPlXVIyISUlu2QKtWsGkTHHfcwf9TVY+ISAxNnAg9ehye9JOlxC8iElLpqOYBJX4RkVBatgy2boV0\nTGqqxC8iEkLjxkHPnjYbZ6qpcVdEJGSKiiAn58BMnKVR466ISIy89RbUrFl20k+WEr+ISMiMHQt9\n+qRv/6rqEREJkR07oFEjm6LhpJPK3k5VPSIiMTFpEpx33pGTfrKU+EVEQuTpp6Ffv/QeQ4lfRCQk\nPvrIpmfo1i29x1HiFxEJiTFjoG/f9PTdL0mNuyIiIVBQYPPuL1gATZsefXs17oqIRNzUqdC6dWJJ\nP1lK/CIiIfD003DjjZk5lqp6REQCtmEDtG8PGzfCsccm9hhV9YiIRNizz8LVVyee9JOlEr+ISID2\n7YNmzeDFF+H00xN/nEr8IiIR9frr8L3vlS/pJ0uJX0QkQH/5C/z855k9pqp6REQCsm4ddOgA69dD\ntWrle6yqekREImj0aOjVq/xJP1kq8YuIBKCgABo2hDlzoFWr8j9eJX4RkYh56SVL+BVJ+slS4hcR\nCUAQjbr7qapHRCTDli2zxVbWrYMqVSq2D1X1iIhEyJ//DP37VzzpJ0slfhGRDPryS2jSBD75BOrW\nrfh+VOIXEYmIp56C7t2TS/rJUolfRCRDCgttvv0pU5KfokElfhGRCJg61VbZyuS8PKVR4hcRyZAR\nI2DgwKCjUOIXEcmIxYth9Wq47LKgI1HiFxHJiD/9CX7xC6hcOehI1LgrIpJ2mzZBmzawahXUqpWa\nfSbTuPud1IRQNufcWuAroAjY670/M93HFBEJkxEjoHfv1CX9ZKW9xO+cWw2c5r3/soz/q8QvIrG1\nYwfk5MCiRdCoUer2G/bunC5DxxERCZ3Ro6Fbt9Qm/WRlqsS/DfDAaO/9k4f8XyV+EYmlggKbnmHa\nNDj11NTuO9R1/MDZ3vstzrmTgJnOuU+89/NKbnDPPfd8+3tubi65ubkZCEtEJL0mToSWLVOT9PPy\n8sjLy0t+R2S4V49zbgiw03s/rMR9KvGLSOwUFUG7dvDoo1bVk2qhreN3zh3nnKte/Hs14MfAR+k8\npohIGEybZn32f/zjoCM5XLqreuoALznnfPGxxnvvZ6T5mCIigfIe7rsP7roLXIXK5OmV1sTvvV8D\ntE/nMUREwmbGDNi9Gy65JOhISqduliIiKeQ9DB0Kd94Jx4Q0w4Y0LBGRaJo7F7ZsgauuCjqSsinx\ni4ik0H33wR13wHcy0Vm+gpT4RURSZMECWLECrr8+6EiOTIlfRCRFBg+2uv0qVYKO5MiU+EVEUmDu\nXJt2+YYbgo7k6JT4RUSS5L312R8yJPylfVDiFxFJ2qxZsHUrXHdd0JEkRolfRCQJ3sPdd8O994a7\nJ09JSvwiIkl49VXYtSvc/fYPpcQvIlJB+/ZZL57f/z68o3RLE6FQRUTCZdw4qF49vHPylCWj8/GX\nGoDm4xeRCPrmG2jRAiZNgrPPzvzxQzsfv4hIXP35z3D66cEk/WSpxC8iUk7//rctqThvnpX6g5BM\niV+JX0SknG67DfbsgccfDy4GJX4RkQxZuRLOOgs+/hjq1AkuDtXxi4hkyG23wW9+E2zST1ZExpmJ\niATvH/+A5cth8uSgI0mOSvwiIgkoKIBbb4Vhw6Bq1aCjSY4Sv4hIAkaNgsaNoXv3oCNJnhp3RUSO\nYutWaN3aum+2bBl0NEa9ekRE0qh3b6hdGx59NOhIDkgm8atxV0TkCN54A/LyYNmyoCNJHdXxi4iU\nYc8e+PnPrX6/evWgo0kdJX4RkTI8+KDV7f/XfwUdSWqpjl9EpBQrVkCnTvDBB9CgQdDRHE4jd0VE\nUqioCH72M1tkJYxJP1lK/CIih/jLX6x+f+DAoCNJD1X1CPv2wb/+BZs3w5df2vqhX38N+fm2nNwx\nx0ClSlCjBpx4ItSqZT9PPDE6i0uLJGr1ajjzzHD12S+NunNKQgoKrL5y0SL45BPrnrZihQ1OqVUL\nTj7Zknn16lCtmg1L996+9u7dCzt2wLZt9uGwbZv93aAB5ORAkybQtCm0aQOnngp16wb9bEXKr6gI\nbrwRbr893Ek/WSrxx9jevbBgAUyfDm+9ZUm/eXNbNah1a2jVyi7uevWgSpXy7z8/H9atgzVrrJT0\n2Wfw4Yd2nMqVoX17Kzmdc45NY1utWuqfo0gqPfaYraM7b559yw0zjdyVb+3ZA6+9ZuuAzpxpJfGL\nLoLcXEvCxx+f/hi8h40b7QNgwQKYOxcWL4a2beHcc+GCC6y3ROXK6Y9FJFH759kPexXPfkr8Wc57\neOcdePpp+NvfoF07uPZam0zq5JODjs7s3m0xzpljH0yrVsH558PFF9vtpJOCjlCyWUGBFUb69IEB\nA4KOJjFK/FkqPx+ef94Wfd62Dfr3h549oX79oCM7ui1bbG7zV1+1byZnnglXXw2XXmrtDCKZdMcd\nsHQpTJsGrkKpNPOU+LPMN9/AE0/Aww9bY+rAgXDhhdb7Jop277YPgOeftw+Bzp3tQ+Cyy+I1TF7C\nac4c6NXLqiZr1w46msQp8WeJ/HxL+A89ZCXkIUOsATVOdu60UtfEiVbXevnl1suiY8folMQkOrZt\ns/fQk09Ct25BR1M+Svwx5z288oqt9dmiBQwdal0m4+6f/4TnnoMxY6yHRb9+B6bHFUlWURH06AGn\nnAJ//GPQ0ZSfEn+MLV8Ot9wCGzbA8OHRK5Wkgvcwf741Xr/8sjVaDxwIZ5wRdGQSZfffbx0N5syJ\nZg8zzdUTQ4WF8MADVt994YXWPz4bkz5YFU/nzjB2rPUGatcOrrzSqn/Gj7ceGSLlMWuWTbX8/PPR\nTPrJUok/hJYuhRtusN4to0fbOp9ysH37rPpr5Egbgdy/v82bXqdO0JFJ2G3caN8WJ0yALl2Cjqbi\nVOKPCe8tkf3oR5bEXn9dSb8slSpZ/ewbb1jpbcsWG3Tzs5/ZCGKR0uzZA1dcYdWnUU76yVKJPyS2\nbbPeKxs22KjbZs2Cjih6Pv/cxjT89a/2ph40CE47LeioJCy8h+uus58TJkS/l5hK/BH3/vvQoYOV\n7ufPV9KvqNq1rcfT6tU29P6SS+C882DGDHuzS3a77z5rIxozJvpJP1kq8Qds0iS4+WYrpV5+edDR\nxEtBgZ3fhx+2Seh+9zsbFBbVgW5ScS++CL/+tU0bEpZpTJKl7pwRVFQEgwfbTIB//7v1VJH0KCqy\nkcFDh9oo4bvusl5BYZ99UVJjwQJbM3fGjHiNf1Hij5g9e+D6661BcsoUDUjKFO/tzX/vvdamctdd\ncM01WkwmzpYts/aesWNtlto4CXUdv3PuAufccufcp865Qek+Xtjt2HHgApw1S0k/k5yzsRDz59u8\n66NHW0+gsWNt7QKJlw0bbArwRx+NX9JPVloTv3PuGGAU0A1oDVzrnDtsputeA3uxZu2adIYSmDVr\n19BrYC+69O3C5Tf14qyz19CihdU9V60adHTZyTlr9J0710YDjxtnw/affPLAYLCSr1ucr884Kfma\nXdG/F11+tIZbb7Vv13Gy/3kmI61VPc65jsAQ7/2FxX/fDnjv/UMltvH8Dpp+2JSZo2aS0zgnbfFk\n2pq1a+g6oCur2q2CKkAB1JrTlPfHz6RJTnyeZxzMn29tAMuWwX//dA3PLuzK6vYHXrc4Xp9xUtp7\nrebspiyaGK/X7KDneT+hrer5AbChxN8bi+87WBVY1W4Vdw+7O83hZNbdw+4+cCECVIEvu6xi8PB4\nPc846NTJ1gd48UV44sW7DyR9iO31GSelvde2/yh+r9lhz7OC0p34S/s0OvwrxhxgPsyfO5+8vLw0\nh5Q5m3ZsOvwFqgKbd2wOJB45uv/4Dzilg163qMmG91peXh7z5s6D+VjOTEK6+zNsBBqW+Ls+cPgr\n0QUogE47O5Gbm5vmkDKnVuUfQAEHX5AFUO+EekGFJAn4wQmlv27bN9Zj1y4tGh9G36sa//dabm4u\nnc/pzLrj19nzfLPi+0p3if89oJlzrpFzrgpwDTD1sK2K61CH3jY0zeFkzpdfwkdzh3JiXlO7ICGW\nzzOOht42lKYfHvy61X+vKSd/dyhNmsCDD9qCMRIO69bBollDqTk7/u+1w67NCkpr4vfe7wMGADOA\nj4FJ3vtPDt3uup3Xxarh7JtvbMDIxRfl8N64mVy38zq6rOkSu+cZVzmNc5g56uDXbe7ombz6Sg6z\nZ8OSJdC0KfzhD/DVV0FHm90+/hj+8z/h1ltyWDQx/u+1ktdmMjSAK8UKC21agBNOsNWjND1APC1f\nbgt5TJ8OAwbYbI81awYdVXaZOxeuugqGDYOePYOOJvNCPYArm3hv0ykXFNhEUEr68dWypX2wv/02\nrF1rE+sNHmwjgiX9Ro+2aTfGjcvOpJ8spaYUGjYM3nsPJk+2ScEk/po3t5G/775rawQ3b26TwX3+\nedCRxVNhoU1qOHw4zJsH558fdETRpMSfItOmWeKfOhWqVw86Gsm0Jk1s5O+iRbB9+4FFYVauDDqy\n+Ni0yRL9Z5/BwoX2ISsVo8SfAkuWQL9+8Le/QcOGR99e4qtRI3j8cWsDqF0bzj7bptt+552gI4u2\n6dPh9NMt8b/yCtSoEXRE0abG3SR98YVdkA88ANdeG3Q0Eja7dll7z7Bh0KAB/Pa3NmGY2n8Ss2cP\n3H23zW01fjycc07QEYWHpmUOyL59cOGF0L69LfYhUpbCQmv7efhhS2YDBtjkYccfH3Rk4bVwoX2T\nbtnSGnO///2gIwoXJf6ADB5sXcpmzdKc7pIY7yEvz6aFnj3b1oD95S8tuYnZvfvAIkUjR1rvnWxf\nKrE06s4ZgFdfta/wkyYp6UvinLOFQSZPtrahmjUhNxe6doWXX7ZvBtnKe5g40T4EN2+GpUutn76S\nfuqpxF8Ba9faZF5TpkDnzkFHI1GXn28fBKNG2eIhvXtD3762RkC2eP99uPVWK+2PGGGjceXIVOLP\noMJCGzDym98o6UtqVK1qVT4LFljvlfx8S3ydOsFTT9mqbXG1aJFNb9KjB/TpY+NglPTTTyX+choy\nxBqdpk9XzwxJn7177Rp75hlrC+ja1eq6L7oo+uNEvIe33oI//tFK+oMGwU03wXe/G3Rk0aLG3Qx5\n6y17833wAZx8ctDRSLb44gur/5882aaIOP98uw4vvtjmhIqKXbusTWzkSPtWc/PN1mvn2GODjiya\nlPgzYPt267Y5ahR07x50NJKttm2z0eGTJ8Obb8IZZ9iC4t26Qdu24WsILSyEN96wHjrTpln16MCB\n9uGlb8zJUeLPgJ494cQTLfGLhMHXX8OcOfD667Zs5O7dcO651jbQuTO0aQOVKmU+rm3bLKbXXrO4\ncnKsDePqq6Fu3czHE1dK/Gk2ZYpNvLV4sb6WSnitWmXjSubPt9vmzTaqvG1b+xBo0wZat4bjjkvd\nMfPz4dNPrZH27betgXrtWuuietFFNsCxUaPUHU8OUOJPo3/9y944U6bYvCsiUfHFF9ZLZunSA7cV\nK+yba8OGdmvQAOrVs/EENWtam0HlyvZN4ZhjrJF51y77drF9u81AunmzTZi2YoV1P83JgXbt7P1x\n1ln2u2anTT8l/jS6+mp7czz6aNCRiCSvsNCS9/r1lrTXr7dE/tVXlth37LBkX1RkU5JUrmy9iKpV\nsw+FevUO3Jo3t3UIlOSDocSfJpMnw513qopHRMJHiT8N/v1v+OEPVcUjIuGkxJ8G/frZzIkjRgQd\niYjI4ZJJ/JperBRvvgkzZ8LHHwcdiYhI6mkIxSHy823JvBEjojUqUkQkUUr8h3jkEeupcOmlQUci\nIpIequMv4bPPoGNH+L//06ATEQk3TcucIrfcYmuiKumLSJypcbfYa6/BypXw0ktBRyIikl4q8QMF\nBfCrX8Hw4RqFKCLxp8SPzQ/etKnNby4iEndZ37i7ZYuN0J0/H1q0CCwMEZFy0cjdJNx4I9SqpUnY\nRCRaNHK3gpYsgVdesfnERUSyRVbX8Q8aBHfdBTVqBB2JiEjmZG3inzXLum/27x90JCIimZWVib+o\nyAZqPfCAum+KSPbJysQ/caKtLHTFFUFHIiKSeVnXqyc/37ptPvccnHNOxg4rIpJSmqunHJ54Atq0\nUdIXkeyVVSX+3bttyuXXXoP27TNySBGRtFCJP0GPPQadOinpi0h2y5oS/86dVtqfPRtat0774URE\n0kol/gSMGAFduyrpi4hkRYl/+3Zo3hzeftt+iohEnUr8RzF8OPzkJ0r6IiKQBSX+r76yufbffRea\nNEnbYUREMkol/iN47DG46CIlfRGR/WJd4t+1yxJ+Xh60apWWQ4iIBEIl/jI88QSce66SvohISWkr\n8TvnhgA/BT4vvut33vt/lLJdWkr8e/ZYaV+jdEUkjsK8Atcw7/2wNB+jVGPGwGmnKemLiBwq3Ym/\nQp9Gydq7Fx56CJ5/Poiji4iEW7rr+H/pnFvsnHvKOZexBQ5feAFycqBjx0wdUUQkOpIq8TvnZgJ1\nSt4FeOBO4HHg995775y7DxgG3Fjafu65555vf8/NzSU3N7fCMXkPjzwC999f4V2IiIROXl4eeXl5\nKdlXRrpzOucaAdO8921L+V9KG3dnzoRf/QqWLgUXSEWTiEj6hbI7p3Oubok/LwM+StexSnrkEfif\n/1HSFxEpSzq7cz4HtAeKgLVAf+/91lK2S1mJf/Fi6N4dVq/WIuoiEm/JlPhjNXK3Vy9o2xZ++9uU\n7E5EJLSU+IH16+HUU620XyNj/YdERIIRyjr+THvsMejdW0lfRORoYlHi370bGjWChQttCmYRkbjL\n+hL/hAk2WEtJX0Tk6CKf+L2HkSNh4MCgIxERiYbIJ/65c6GwEM4/P+hIRESiIfKJf+RIGDBAA7ZE\nRBIV6cbddeugQwf7Wb16igMTEQmxrG3cffxx6NNHSV9EpDwiW+JXF04RyWZZWeJXF04RkYqJZOL3\nHkaNgptvDjoSEZHoiWTif+892LlTXThFRCoikol/9Gj46U/hmEhGLyISrMg17u7YYY26n3wCdese\nfXsRkTjKqsbd8eOtikdJX0SkYiKV+L2HJ56Am24KOhIRkeiKVOJ//31r1D3vvKAjERGJrkglfjXq\niogkLzKNu2rUFRE5ICsadydOtCoeJX0RkeREIvGrUVdEJHUikfgXLYLt2zVSV0QkFSKR+J95Bvr2\nVaOuiEgqhL5xNz8f6te3+XkaN85cXCIiYRbrxt1XX4Uf/lBJX0QkVUKf+PdX84iISGqEuqpn61Zo\n0QI2btTyiiIiJcW2qmfCBLjkEiV9EZFUCnXif+YZW0xdRERSJ7SJf/Fi+OorOPfcoCMREYmX0Cb+\nZ56B3r3Vd19EJNVC2bhbUAANGsDbb0PTpgEFJiISYrFr3J0+HU45RUlfRCQdQpn4n31WjboiIukS\nuqqe7dtt3v1166BmzQADExEJsVhV9bz0ks27r6QvIpIeoUv8EybAtdcGHYWISHyFqqpnyxZo1Qo2\nb4Zjjw00LBGRUItNVc8LL8BPfqKkLyKSTqFK/BMnQs+eQUchIhJvoanqWb0aOnaETZugcuVAQxIR\nCb1YVPVMmgRXXqmkLyKSbqFJ/OrNIyKSGaFI/EuXws6dcPbZQUciIhJ/oUj8EyfCNddoJk4RkUxI\nKtU6565wzn3knNvnnOtwyP/ucM6tdM594pz78ZH2o948IiKZk2wZeylwKfBmyTudc62Aq4BWwIXA\n4865MlufTzgB2rZNMhIBIC8vL+gQYkXnM3V0LsMjqcTvvV/hvV8JHJrUewCTvPeF3vu1wErgzLL2\n8/77UPbHgpSH3lyppfOZOjqX4ZGuWvUfABtK/L2p+L5SqQuniEjmfOdoGzjnZgJ1St4FeOBO7/20\nsh5Wyn3BjhQTEREgRSN3nXNzgF977xcV/3074L33DxX//Q9giPf+nVIeqw8EEZEKqOjI3aOW+Muh\nZABTgfHOueFYFU8z4N3SHlTRwEVEpGKS7c55iXNuA9AReMU5Nx3Ae78MeAFYBrwG/OKwFdVFRCQQ\ngU/SJiIimZXxsbJHGvR1yHYXOOeWO+c+dc4NymSMUeKcq+Wcm+GcW+Gce905V6OM7fY55xY55z5w\nzr2c6TjD7GjXmnOuinNuUvGAxAXOuYZBxBkVCZzPPs65z4uvx0XOuX5BxBkFzrmnnXNbnXNLjrDN\nyOJrc7Fzrn0i+w1ikoRSB32V5Jw7BhgFdANaA9c651pmJrzIuR2Y5b1vAcwG7ihju13e+w7e+1O9\n95dkLrxwS/BauxHY5r1vDvwJeDizUUZHOd67k4qvxw7e+zEZDTJaxmLnslTOuQuBpsXXZn/gr4ns\nNOOJ/wiDvko6E1jpvV/nvd8LTMIGhcnhegDPFv/+LFBWUlcjeukSudZKnuPJwHkZjC9qEn3v6npM\ngPd+HvDlETbpATxXvO07QA3nXJ0jbA+EZJK2Uhw6AGwjRxgAluVqe++3AnjvtwAnlbFdVefcu865\nt51z+hA9IJFr7dttvPf7gO3OuRMzE17kJPrevay4auIF51z9zIQWS+UaLLtfKrtzfquCg74O2kUp\n92VtK/QRzudd5dhNQ+/9FudcDjDbObfEe78mlXFGVCLX2qHbuFK2EZPI+ZwKTPDe73XO9ce+Telb\nVMVUKFemJfF777smuYuNQMkGtPrA5iT3GVlHOp/FDT91vPdbnXN1gc/L2MeW4p9rnHN5wKmAEn9i\n19oGoAGw2TlXCTjBe3+kr9/Z7Kjn85Bz9yTwUAbiiquN2LW5X0K5MuiqnrLq+d4DmjnnGjnnqgDX\nYKUEOdxUoG/x732Avx+6gXOuZvF5xDn3feBsbIyFJHatTcPOLcCVWCO6lO6o57O4gLJfD3QtHo2j\n7Fw5FegN4JzrCGzfX/V7RN77jN6wxscNwDfAP4HpxfefDLxSYrsLgBXYzJ63ZzrOqNyAE4FZxedq\nJlCz+P7TgNHFv58FLAE+AD4E+gYdd5hupV1rwL1A9+Lfq2IDElcCC4HGQccc5lsC5/N+4KPi6/EN\n4JSgYw7rDZiAleDzgfXADVjvnZtKbDMK+Kz4vd0hkf1qAJeISJYJuqpHREQyTIlfRCTLKPGLiGQZ\nJX4RkSyjxC8ikmWU+EVEsowSv4hIllHiFxHJMv8PTNjedI8P5SQAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7f4a13d9a358>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "def eval_legendre_deriv(n, x):\n",
        "    return (\n",
        "        (x*sps.eval_legendre(n, x) - sps.eval_legendre(n-1, x))\n",
        "        /\n",
        "        ((x**2-1)/n))\n",
        "\n",
        "brackets = sps.legendre(n-1).weights[:, 0]\n",
        "\n",
        "nodes = np.zeros(n)\n",
        "nodes[0] = -1\n",
        "nodes[-1] = 1\n",
        "\n",
        "from functools import partial\n",
        "\n",
        "# Use the fact that the roots of P_{n-1} bracket the roots of P_{n-1}':\n",
        "for i in range(n-2):\n",
        "    nodes[i+1] = bisection(\n",
        "        partial(eval_legendre_deriv, n-1),\n",
        "        brackets[i], brackets[i+1])\n",
        "\n",
        "mesh = np.linspace(-1, 1, 300)\n",
        "\n",
        "pt.plot(mesh, eval_legendre_deriv(n-1, mesh))\n",
        "pt.plot(nodes, 0*nodes, \"o\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Node Selection: Gauss-Radau"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For Gauss-Radau (with the left endpoint included), the nodes are the roots of the following function:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 31,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.5/dist-packages/ipykernel/__main__.py:5: RuntimeWarning: invalid value encountered in true_divide\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7f4a13b73588>]"
            ]
          },
          "execution_count": 31,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHSFJREFUeJzt3Xe0VOW5x/HvQxVCUUQpYogNUBRBEEFBjxIVNBdjYjfG\nRE3MjV5NzI3Gm0RxJWZdzVJjTZPYXaDxqlgRgmPUIIKA9GIQBWmGroKU894/nnOUHE+fPbPL/D5r\nzeKIm70fNnueeectz2shBEREJFuaxB2AiIhET8ldRCSDlNxFRDJIyV1EJIOU3EVEMkjJXUQkg5pF\ncRIzWwpsBMqB7SGEgVGcV0REGieS5I4n9bIQwvqIziciInmIqlvGIjyXiIjkKaqEHIDxZjbVzL4X\n0TlFRKSRouqWOTqEsMrM9gImmNn8EMJrEZ1bREQaKJLkHkJYVfHrh2b2JDAQ+LfkbmYqYiMi0ggh\nBGvon8m7W8bMWptZm4qfvwScBMyp7tgQgl4Rva6//vrYY8jKS/dS9zPJr8aKouXeCXiyomXeDHgk\nhPBSBOcVEZFGyju5hxDeBfpGEIuIiERE0xdTqqysLO4QMkP3Mlq6n8lg+fTpNOhCZqFY1xIRyQoz\nI8QxoNoQzzxTzKuJiJSuoib3Rx4p5tVEREpXUZP7/PnFvJqISOkqap97q1aBzZuhadOiXFJEJPVS\n0ee+116wdGkxrygiUpqKmtx79YIFC4p5RRGR0lTU5H7wwep3FxEpBiV3EZEMUnIXEcmgoif3BQtA\nC1VFRAqrqMm9Y0cwgzVrinlVEZHSU9TkbqauGRGRYih6VUgldxGRwosluWuuu4hIYRU9uffqpZa7\niEihqVtGRCSDir5ZR3k5tG0Lq1b5ryIiUrNUFA4DaNIEevSAhQuLfWURkdIRyx6q6poRESmsWJK7\nBlVFRApLLXcRkQyKLblrrruISOEUfbYMwKefwu67w4YN0LJlUS4vIpJKqZktA57Q99sPFi2K4+oi\nItkXS3IHOPRQmDMnrquLiGRbbMm9d28ldxGRQoksuZtZEzObbmbj6nO8Wu4iIoUTZcv9SmBefQ9W\nchcRKZxIkruZdQNOAe6t75854ABYuRI+/jiKCEREZFdRtdxvA34K1HteZbNm0LOnFjOJiBRC3snd\nzE4FVocQZgJW8aoXdc2IiBRGswjOcQww0sxOAVoBbc3swRDCt6seOGrUqM9+Lisro3fvMiV3EZFd\n5HI5crlc3ueJdIWqmR0H/CSEMLKa/xeqXuvZZ+Guu+DFFyMLQUQkU1K1QrWSumVERAojltoylcrL\noX17WLbMa82IiMi/S2XLvUkTX6k6d26cUYiIZE+syR1UhkBEpBBiT+7qdxcRiZ6Su4hIBiUiuc+e\nDUUa1xURKQmxJ/fOnT2xr1kTdyQiItkRe3I3U9eMiEjUYk/uAH36wNtvxx2FiEh2JCK5H364kruI\nSJSU3EVEMijW8gOVtmyBDh1g40Zo0aIo4YiIpEIqyw9UatUK9ttPG3eIiEQlEckd1DUjIhIlJXcR\nkQxSchcRyaDEJPe+fWHmTJUhEBGJQmKSe+fOXt99xYq4IxERSb/EJHczdc2IiEQlMckdlNxFRKKS\nqOTet6+Su4hIFBKV3A8/3AdVRUQkP4koP1Bp+3Zo3x7+9S9o3booYYmIJFqqyw9Uat4cevVSbXcR\nkXwlKrmDBlVFRKKQyOSufncRkfwkLrn36wfTp8cdhYhIuiVqQBVg0ybo0sVruzdrVoTAREQSLBMD\nqgDt2kG3bqrtLiKSj8Qld4ABA+Ctt+KOQkQkvfJO7mbW0symmNkMM5ttZtfne87+/ZXcRUTykXdy\nDyF8ChwfQugH9AVGmNnAfM6p5C4ikp9IumVCCJ9U/NgSaAbkNUrbr5/Pdd+xI+/QRERKUiTJ3cya\nmNkMYBUwIYQwNZ/zaVBVRCQ/kUw2DCGUA/3MrB3wlJkdEkKYV/W4UaNGffZzWVkZZWVlNZ6zsmvm\nsMOiiFBEJB1yuRy5XC7v80Q+z93MrgM+CiHcWuX36zXPvdItt8DSpXDnnZGGJyKSKrHNczezjmbW\nvuLnVsBXgQX5nleDqiIijRdFt0wX4AEza4J/WIwNITyf70n79YNZs3xQVStVRUQaJnHlB3bVowc8\n8YT63UWkdGWm/MCutFJVRKRxEp3c1e8uItI4Su4iIhmU6D73jRthn31gwwYNqopIacpkn3v79rDv\nvtpTVUSkoRKd3AEGDYIpU+KOQkQkXRKf3I86Ct54I+4oRETSJfHJXS13EZGGS/SAKvgK1d13h+XL\n/VcRkVKSyQFV8Fky/fvD1LyKCIuIlJbEJ3fwfnd1zYiI1F8qkvugQRpUFRFpiMT3uQN88AH07Qtr\n1oA1uOdJRCS9MtvnDr5KtWVLWLIk7khERNIhFckdNCVSRKQhUpPcNagqIlJ/qUnuGlQVEam/VAyo\nAnzyCey1F6xb5/3vIiKlINMDqgCtW0PPnjBjRtyRiIgkX2qSO8DgwTB5ctxRiIgkX6qS+5Ah8Npr\ncUchIpJ8qelzB1i2zOvMrF6txUwiUhoy3+cOvitTq1aweHHckYiIJFuqkjuoa0ZEpD5Sl9yHDoVX\nX407ChGRZEtdclfLXUSkbqlL7occAmvXwqpVcUciIpJcqUvuTZrA0Uer9S4iUpu8k7uZdTOzSWY2\nz8xmm9kVUQRWm6FDldxFRGoTRct9B3BVCOEQYDBwmZn1iuC8NVK/u4hI7fJO7iGEVSGEmRU/fwTM\nB/bJ97y1GTAAFiyAzZsLeRURkfSKtM/dzL4C9AUKWnm9ZUvo108lgEVEatIsqhOZWRvgr8CVFS34\nLxg1atRnP5eVlVFWVtbo61XOdz/xxEafQkQkcXK5HLlcLu/zRFJbxsyaAc8CL4QQbq/hmLxry+zq\n+efh5pshgnsgIpJYja0tE1VyfxD4VwjhqlqOiTS5b9oEXbvChx96vRkRkSyKrXCYmR0DnA+cYGYz\nzGy6mQ3P97x1adcODjtM9d1FRKqTd597COF1oGkEsTTYCSfApEn+q4iIfC51K1R3VZncRUTk36Vq\ns46qtmzxTbNXroS2bSM9tYhIIpTEZh1VtWoFRx6p1aoiIlWlOrmDumZERKqT+uR+/PFK7iIiVaW6\nzx1g2zbo2BGWLoUOHSI/vYhIrEqyzx2gRQuv7/7KK3FHIiKSHKlP7uD97i+/HHcUIiLJkZnkrn53\nEZHPpb7PHWDnTp/vPm8edO5ckEuIiMSiZPvcAZo2hWHDYPz4uCMREUmGTCR3gOHD4cUX445CRCQZ\nMtEtA7B8OfTtC6tXe0teRCQLSrpbBqBbN+jSBaZNizsSEZH4ZSa5g7pmREQqZSq5jxih5C4iAhnq\ncwf49FPYe29YsgT23LOglxIRKYqS73MHaNkSjjsOJkyIOxIRkXhlKrmD+t1FRCBj3TLgXTJHHw0r\nVkCTzH10iUipUbdMhf33h3btYNasuCMREYlP5pI7+KyZ556LOwoRkfhkMrmPHAnjxsUdhYhIfDLX\n5w6wfTt06gRz5kDXrkW5pIhIQajPfRfNm8Mpp6j1LiKlK5PJHeC00+Dpp+OOQkQkHpnslgHYvBn2\n2cerRbZrV7TLiohESt0yVbRtC0OGaEGTiJSmSJK7mY02s9VmlqjZ5eqaEZFSFUm3jJkNAT4CHgwh\n9KnhmKJ2y4CvUj30UN/Ao3nzol5aRCQSsXbLhBBeA9ZHca4ode0KPXrAK6/EHYmISHFlts+9krpm\nRCStVq9u/J9tFl0YdRs1atRnP5eVlVFWVlbwa55+OgwbBrffrkJiIpJ8uVyOXC5HCPD73zf+PJFN\nhTSz7sAzSepzr3T44XDXXTB0aCyXFxFpsClT4FvfgnfeiX8qpFW8Euecc2DMmLijEBGpv3vvhYsv\nbvyfj2q2zKNAGbAnsBq4PoRwX5VjYmu5L1kCgwfDBx9As6J2RImINNxHH8G++8K8edC1a+Na7pGk\nuhDCeVGcp1D23x+6d4eXX4YTT4w7GhGR2j32GBx7LHTp0vhzlMwQ4znnwNixcUchIlK30aPhkkvy\nO0dma8tUtWwZ9O0LK1dCixaxhSEiUqvZs30v6Pfe825k1Zapw777wiGHwEsvxR2JiEjN7r4bLr00\n//HBkmm5g9+0yZPh4YdjDUNEpFobN8JXvuIDqZX97Wq518MZZ8Czz8Inn8QdiYjIFz3wAJx8cn4D\nqZVKKrl36gRHHaVyBCKSPCHAPffAD38YzflKKrkDfPe7cN99dR8nIlJMEyd69dqoVtKXVJ87wNat\nvkPTjBnw5S/HHY2IiBs+HM4+2xugu1Kfez3ttpvPeX/ggbgjERFxs2fDrFlwXoTLQUsuucPnXTPl\n5XFHIiICt94Kl18OLVtGd86SrLTSvz+0aQN//zsUoeqwiEiNVqyAp56Cf/4z2vOWZMvdzFvvf/lL\n3JGISKm74w44/3zo0CHa85bcgGqlDz+Egw6C99+Hdu3ijkZEStHatb4V6PTpXtywOhpQbaC99oIT\nTlAxMRGJzx13+G5xNSX2fJRsyx3ghRfgF7+AadO8q0akvsrLva90xQovRrdypS8d/+gj2LwZtm3z\n48ygaVNo2xbat/dviXvvDd26+ZTczp39/0vp2bgRDjjAd1w64ICaj2tsy72kk3t5uX8levhhGDQo\n7mgkqdatg6lT4c03fcrawoWweLEn6332ga5dPUnvsYcP1Ldp47MeQvDXjh2e9Ddt8jf0qlW+ccwH\nH8CGDf4MHnywF7br3RsGDvRCd5JtN94ICxbAQw/VfpySeyPdcgvMnFn3DZbSsWKFVw+dONFbVatX\n+wyrgQN9P95evTwht2mT/7U+/tg/LObN89fs2X7Nli3hmGPg6KPhuOOgTx99u8yS9ev9GXr1VX+e\naqPk3kjr1vlXooUL/euylJ7ycnj9dRg3DsaP9xb1sGFw0km+PWOvXsXtOgnBp8W9/rq/Jk2CLVvg\nlFP8NWyYJgGk3f/8jzcaRo+u+1gl9zxcfDEceCBce23ckUixhOAt5LFjfUuzjh19YGv4cDjyyOT1\ngy9eDM8/76/Jk701f845MHKk9+dLeqxcCYceWv8SKErueZgxw98kS5Z44R7JrhUrfH3D6NGfl6I4\n++y6vxonyaZN/i1jzBj/Wn/SSfCtb8Gpp2oD+DS47DLvdrv11vodr+SepxNO8Bb8+efHHYlErbzc\n+9D/+EfI5eCss+D734cjjkh/P/a6dfDEE14rackSX5x3ySWw335xRybVWbjQx1Lmz/fp2PWh5J6n\n557zaZHTp6f/DS9uyxZ48EEfNG/TBn7wAzj33Ox2Y8ydC/fe67O/jjgCrrgCRoyAJiW7miV5Tj3V\nG5I/+Un9/4ySe57Ky70f7M47fcBK0mvtWt/04O67vf/8pz/1Gtml8qG9dSs8/jjcdpt/wP34x3DB\nBdCqVdyRlbYXXoArr4Q5c6BFi/r/Oa1QzVOTJv5pesstcUcijbV2rQ+K9+gBS5f6LJNnnoFjjy2d\nxA4+lnDBBfDWW/D73/vWkt27w3XX+QwNKb7t2+Gqq7yfvSGJPR9K7rs4/3yf8z5zZtyRSEOsXw+/\n/KUn9Q0b/N9v9GhfFFTKzLzq6bhxPvD64Ye+WOrKK326pxTPrbf6xtennlq8ayq572K33eC//xt+\n/eu4I5H6+OQTX+V30EG+6rOyparVnV/Us6ffm7lzfZrnYYf5rI333487suz75z/ht7/1rsJifoNU\ncq/i0ku9lTN3btyRSE3Ky31Fca9e3kp/4w3485+9ZSS169LFW5ELFvggc9++PnNISb4wQoD//E+4\n+uriz2BScq/iS1/yAagbb4w7EqlOLueDpPfc4/O8H3/cF6BJw+y9N9x0Eyxa5Au4+vXz5/7DD+OO\nLFsefhjWrPF7W2yRJHczG25mC8xskZldE8U543TZZTBhgrduJBmWLIGvf93ncV99NfzjH153RfLT\nsSP85jf+TXXHDv82dP31vlBK8rNsmU/SuO++eBZH5p3czawJcBdwMtAbONfMUrTe74vatvV/lOuu\nizsS2boVfvUrL9p11FG++OPss0tr9ksxdO7s04CnTYN33/VxjFtv9fsvDVdeDt/5DvzoR/6tKA5R\ntNwHAotDCO+FELYDY4DTIjhvrK64wos2TZsWdySla/x4H/ibPt0HS6+91ge9pXD2288Xfv3tb77H\ncI8e3vLcuTPuyNLljjv8g/Hqq+OLIYrkvg+wbJf/Xl7xe6nWurVPr1MxseJbvhzOPNMHon73O3jy\nycLsVCM1O/RQ37R57FhP7ocf7msGErwOMTGmTPGuroceirfWTxTJvbovyJl4BC6+2BfDTJwYdySl\nYedO7wro29fnqM+dW9x5wfJFgwfDK6/A//6vl6k99lgf75DqrV3rtYv+9CfYf/94Y4nic2U5sGvh\nym7AiuoOHDVq1Gc/l5WVUVZWFsHlC6d5c3+or7rKuwZUca9w5syBiy7y6XmTJ3ufrySDGXzta16n\n5qGHvJJm//7eOj344LijS46dO70651ln+eB/Y+VyOXK5XN7x5F1bxsyaAguBYcBK4E3g3BDC/CrH\nJbq2TE1C8EI/Z54JP/xh3NFkz7Zt/gF6550+/fR739NgadJt2eJ1e266yZPYqFG+3WCp+/GPYdYs\nePHFaGfHxFZbJoSwE7gceAmYC4ypmtjTzAxuv90f4HXr4o4mW6ZNgwEDfG/SGTN8MY0Se/K1auUr\nuRctgj339C0Ar73WSz+Uqnvu8aT+178mZ08IVYWsp8su81/vvjveOLJgyxb/sLz/fu9jP+88JfU0\nW77c/z3HjYNrrvH3SinNanrqKR/8f+0137IzaqoKWWC/+hX83//5UndpvFdf9ZkX777rX2HPP1+J\nPe26dfM68rmc//v27OnTKUth+uRLL/k3zmefLUxiz4da7g0wZoz3C7/1VvHKdmbF5s3+1f3JJ+Gu\nu3y/Usmm11/3FvymTT6eMmJENj/AK3f1evJJ312pUNRyL4Kzz/YNbX/727gjSZeXXvLFSB9/7LNi\nlNiz7ZhjvAX/61/7RinHH+9zv7Nk3DhP7GPHFjax50Mt9wZ6/32fBjZpkicsqdn69V7GYdIk37/0\n5JPjjkiKbccO76K5/novH3Hjjd5tk2b33w8/+5kv6jryyMJfTy33IqlsuZ93nupu1Oapp3yVY+vW\nMHu2EnupatbM1y8sWuT1gYYM8bngc+bEHVnD7djx+X4PL79cnMSeD7XcGyEE/0q2774+20M+t3q1\n1+WZMcMH2Y49Nu6IJEk2boQ//MHLSgwY4OMwaajuuWoVfPvb/t4fOxY6dCjetdVyLyIzf0Affxye\nfz7uaJIhBF+92KePF596+20ldvmi9u19sHXJEjjlFG/FH3ccPPdccmfXjBvnlR0HDvRNrouZ2POh\nlnseXnsNvvlNXy4fdx2JOL33HvzgB7Bype9d2r9/3BFJWuzYAY895t+A163z5+iii7zOfNyWL/cx\no6lTfdxgyJB44lDLPQZDhsDPf+4JfsuWuKMpvvJyX9Q1YIC30qdOVWKXhmnWzMevpk3z7o75872u\n0AUXeNnhOFrzGzbADTd4AbuePX18IK7Eng+13PMUgn+13LkTHn0UmpTIx+XChXDJJf73v/de38FH\nJArr1sEDD/gWdatWeaGys87yAcxCvr+WLfNZXX/8oxdK++Uvk/GNXC33mJh5V8Ty5V4SNeu2bvWl\n5kOG+Jvu739XYpdodejgRbjeesvLbbdu7eW3O3f23Y3GjoUV1dadbbj1673L5T/+w1dOb9zoJY3v\nuy8ZiT0farlHZO1ar339ox9lt3rkiy/C5Zf74NJtt/myc5FiWbrUB17Hj/cE3KaNv+f69PEdo3r2\n9KnKbdt+cUVsCJ64lyzxrp8ZM7xhMn8+DBsG3/iGv9q0ieWvVqvGttyV3CO0ZImP/N9wgw8KZcWy\nZd6SmjnTSwcMHx53RFLqQvC585Mnw7x53k24aJE/q1u3+qyctm19XGjbNm+h77ab7+h18MH+gTB0\nqHf1tGoV99+mdo1N7tp+IkL77++DQCecAE2bwoUXxh1Rfj791Msd33yzt9gffri0qv1Jcpl5S726\n1a7bt3tdm02b/H3YvDnssUfpPbtK7hHr0cP7Cb/6VR91v/LKuCNquBC8LvU113iJBe2MJGnSvLnX\nmd9zz7gjiZeSewH06uWV8YYP94HWm25KzyyaKVN8W8FPPvGB4uOPjzsiEWmMlKSc9One3RP8lCk+\nEp/0XZzefturNX7zm77V3bRpSuwiaabkXkAdOngffI8evtAniWVPZ83yhD5ihA8GL17s082aNo07\nMhHJh5J7gTVv7tMGb74ZTjvN+7HjriYZgtfbPv10r9Z4zDHwzjs+jTPpMwdEpH6U3IvkjDO8lfzu\nu9C7tw9YFntm6Nat8MgjPv3r4ot90Pedd7yPvXXr4sYiIoWlee4xmDjR60LvtpsX/R85snADriF4\nzZf77/eVff37w3/9F5x6anoGeUVKmRYxpczOnfDEE77xx8aNvujpvPN8hV2+duzwipVPP+3lSsH7\n0S+8MJrzi0jxKLmnVAi+lPqhh7yrpnt37y4ZOtRrXXTrVvvmwtu3+6q8hQt9PvrkyfDmmz4vfeRI\n7+fv0yebGxSLlAIl9wzYvt1n1EyYAG+84X30H38MXbtCp04+2Nm8ua8c/egj3/Vo+XLo0gUOPND3\nqBw8GAYNSkY9bBHJn5J7Rm3Y4JtgrFrlA6LbtkHLll43Y++9vaXfokXcUYpIoSi5i4hkkOq5i4jI\nZ5TcRUQyKK/kbmZnmNkcM9tpZkdEFZSIiOQn35b7bOB04JUIYpEGyOVycYeQGbqX0dL9TIa8knsI\nYWEIYTGgWdRFpjdQdHQvo6X7mQzqcxcRyaA6N+swswlAp11/CwjAz0MIzxQqMBERabxI5rmb2cvA\nT0II02s5RpPcRUQaIe4Nsmu9eGOCExGRxsl3KuTXzWwZMAh41sxeiCYsERHJR9HKD4iISPEUbLZM\nfRc4mdlwM1tgZovM7JpCxZN2ZraHmb1kZgvNbLyZta/huJ1mNt3MZpjZU8WOM8nqetbMrIWZjTGz\nxWY22cxU/b4W9bifF5rZmorncbqZXRRHnGlgZqPNbLWZzarlmDsqns2ZZta3rnMWcipknQuczKwJ\ncBdwMtAbONfMehUwpjT7GTAxhNATmARcW8NxH4cQjggh9AshfL144SVbPZ+1i4F1IYSDgN8BNxc3\nyvRowHt3TMXzeEQI4S9FDTJd7sPvZbXMbARwQMWzeSnwh7pOWLDkXs8FTgOBxSGE90II24ExwGmF\niinlTgMeqPj5AaCmxK2B6+rV51nb9R7/FRhWxPjSpr7vXT2P9RBCeA1YX8shpwEPVhw7BWhvZp1q\nOT72RUz7AMt2+e/lFb8nX7R3CGE1QAhhFbBXDce1NLM3zewfZqYPys/V51n77JgQwk5gg5l1KE54\nqVPf9+43KroRHjOzbsUJLZOq3u8PqCNX5jUVMoIFTtV9qpfsCG8t9/MXDTjNl0MIq8xsP2CSmc0K\nIbwbZZwpVZ9nreoxVs0x4upzP8cBj4YQtpvZpfi3In0bapwG58q8knsI4cR8/jz+ab/roFU3YEWe\n50yt2u5nxWBLpxDCajPrDKyp4RyrKn5918xyQD9Ayb1+z9oyYF9ghZk1BdqFEGr7qlzK6ryfVe7d\nn4GbihBXVi3Hn81KdebKYnXL1NTvNhU40My6m1kL4Bz8016+aBzwnYqfLwSernqAme1ecR8xs47A\n0cC8YgWYcPV51p7B7y3AmfjAtVSvzvtZ0QipdBp6Futi1JwrxwHfBjCzQcCGym7aGoUQCvLCB/yW\nAVuAlcALFb/fBXh2l+OGAwuBxcDPChVP2l9AB2Bixb2aAOxe8fv9gT9V/DwYmAXMAN4GvhN33El6\nVfesATcAX6v4uSXwWMX/fwP4StwxJ/lVj/v5G2BOxfP4N6BH3DEn9QU8irfEPwXeB76Lz4r5/i7H\n3AW8U/HePqKuc2oRk4hIBsU9W0ZERApAyV1EJIOU3EVEMkjJXUQkg5TcRUQySMldRCSDlNxFRDJI\nyV1EJIP+H1tMyOTGOngBAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7f4a13c11908>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "def radau_func(n, x):\n",
        "    return (\n",
        "        (sps.eval_legendre(n-1, x) + sps.eval_legendre(n, x))\n",
        "        /\n",
        "        (1+x))\n",
        "\n",
        "nodes = None\n",
        "# Root finding left as an exercise for the reader. :)\n",
        "\n",
        "mesh = np.linspace(-1, 1, 300)\n",
        "pt.plot(mesh, radau_func(n, mesh))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Finding the weights"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Use method of undetermined coefficients to find the interpolatory quadrature rule for `nodes`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "max_degree = len(nodes) - 1\n",
        "powers = np.arange(max_degree+1)\n",
        "\n",
        "Vt = nodes ** powers.reshape(-1, 1)\n",
        "\n",
        "a, b = -1, 1\n",
        "rhs = 1/(powers+1) * (b**(powers+1) - a**(powers+1))\n",
        "\n",
        "weights = la.solve(Vt, rhs)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now compare the approximate integrals of monomials from our rule with the true answers:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Error at degree 0: 0\n",
            "Error at degree 1: -1.66533e-16\n",
            "Error at degree 2: 1.11022e-16\n",
            "Error at degree 3: -2.498e-16\n",
            "Error at degree 4: 2.22045e-16\n",
            "Error at degree 5: -3.21965e-15\n",
            "Error at degree 6: -8.88178e-16\n",
            "Error at degree 7: -4.4964e-15\n",
            "Error at degree 8: 0.0145125\n",
            "Error at degree 9: -5.02376e-15\n",
            "Error at degree 10: 0.0339253\n"
          ]
        }
      ],
      "source": [
        "for i in range(2*n + 1):\n",
        "    approx = weights @ nodes**i\n",
        "    true = 1/(i+1)*(1. - (-1)**(i+1))\n",
        "    \n",
        "    print(\"Error at degree %d: %g\" % (i, approx-true))"
      ]
    },
    {
      "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
}