{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from __future__ import division\n",
      "\n",
      "import numpy as np\n",
      "import numpy.linalg as la\n",
      "import scipy.special as sps\n",
      "import matplotlib.pyplot as pt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Set parameters\n",
      "\n",
      "n = 5\n",
      "left_end_nodes = 0\n",
      "right_end_nodes = 0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Ignore (don't even evaluate) this cell at the outset.\n",
      "\n",
      "# If enabled, it computes weights for variants of Gauss \n",
      "# quadratures that include interval ends.\n",
      "\n",
      "# Enable (only) this for Gauss-Radau:\n",
      "if 1:\n",
      "    left_end_nodes = 1\n",
      "    \n",
      "# Enable (only) this for Gauss-Lobatto:\n",
      "if 0:\n",
      "    left_end_nodes = 1\n",
      "    right_end_nodes = 1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "deg_reduction = left_end_nodes+right_end_nodes\n",
      "\n",
      "nodes = np.empty(n)\n",
      "nodes[left_end_nodes:n-right_end_nodes] = \\\n",
      "    sps.legendre(n - deg_reduction).weights[:, 0] \n",
      "\n",
      "nodes[:left_end_nodes] = -1\n",
      "nodes[n-right_end_nodes:] = 1\n",
      "\n",
      "mesh = np.linspace(-1, 1, 300)\n",
      "\n",
      "pt.plot(mesh, sps.eval_legendre(n - deg_reduction, mesh))\n",
      "pt.plot(nodes, 0*nodes, \"o\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "[<matplotlib.lines.Line2D at 0x38fac50>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD7CAYAAABt0P8jAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlclWXaB/DfUVDcl1JUoFBAcSgRl2xzIhVNTMZ2rSaz\nZaxUfFudZqYGP2OljuOUWi69jqOVS/U64YJrdbRUJAV1XFjdAJVUtNBStuf94xJXQDjPOed+lt/3\n8+mj6MN5rh6fc3Gf67nv63ZomqaBiIhso47qAIiIyLuY+ImIbIaJn4jIZpj4iYhshomfiMhmmPiJ\niGzGR3UAFbp27YqdO3eqDoOIyDQiIyOxY8eOWn+fwyjz+B0OBwwSiuklJCQgISFBdRiWwevpXrye\n7uNq3mSph4jIZpj4iYhshonfgqKjo1WHYCm8nu7F66kea/xERCbFGj8REdUIEz8Rkc0w8RMR2QwT\nPxGRzTDxExHZDBM/EZHNMPETEdkMEz8Rkc0w8RMR2QwTPxGRzTDxExHZDBM/EZHNMPETEdkMEz8R\nkc3oSvzPPPMM/P39ceutt1Z5THx8PMLCwhAZGYm0tDQ9pyMiIjfQlfhHjBiB1atXV/n3SUlJyM7O\nRlZWFubMmYMXX3xRz+mIiMgNdCX+3r17o0WLFlX+/bJlyzB8+HAAQK9evXD69GkUFBToOSUREenk\n0Rp/fn4+goKCLn4dGBiIvLw8T56SiMgWystd/14f94VRuau3BXM4HFUem5CQcPH30dHR3JuTiOgy\nTqcTTqcTAHDqlOuv49HEHxAQgNzc3Itf5+XlISAgoMrjL0/8RER0pcsHxGvXAtOmjXfpdTxa6omL\ni8OCBQsAAMnJyWjevDn8/f09eUoiIlvIznb9e3WN+IcNG4YNGzbgxIkTCAoKwvjx41FSUgIAGDly\nJGJjY5GUlITQ0FA0atQI8+bN03M6IiK6QE/id2hXF+EVcTgc1zwPICKiysXFAcuXu5Y3uXKXiMiE\nLDPiLy/XUM2kHyIiAlBWBjRuDJw7Z4ER/7FjqiMgIjK+vDzghhtc/35DJX49H12IiOwiOxsIDXX9\n+5n4iYhMhomfiMhmLJX4s7JUR0BEZHxZWUBYmOvfb6jEzxE/EdH16R3xG2o6Z5MmGn76CZzSSURU\nhfJymcr5449AkyYWmM5Zrx5w/LjqKIiIjOvIEaBZM0n+rjJU4g8NZbmHiKg6ess8ABM/EZGpWDLx\nc2YPEVHVsrIsmPg54iciqlp2tr6pnAATPxGRqbij1GOo6ZwnTmgICZG9JDmlk4joSpoms3kqZva4\nuo+JoUb8LVtKwi8sVB0JEZHxHDsGNGokSV8PQyV+h4MPeImIquKOB7uAwRI/wDo/EVFV3FHfBwyY\n+MPCgMxM1VEQERmP3uZsFQyX+MPDgYwM1VEQERlPerrkSL0MmfjT01VHQURkPO5K/IaazqlpGs6c\nAfz9gaIioI7hfiwREalRUgI0aQKcPg34+cmfWWI6JyBzVFu2BA4fVh0JEZFxHDgAtGt3KenrYbjE\nD7DcQ0R0NXeVeQAmfiIiU7B84u/UiTN7iIgul5Fh8cTPET8R0ZUsP+Jn4iciukTTgH37pBriDoZM\n/AEBwJkzwE8/qY6EiEi9Eyck+bdu7Z7XM2TidziAjh1Z5yciAi7V993Vrt6QiR9guYeIqII76/sA\nEz8RkeHZKvGz1ENEJLnQXQ92AQMn/k6dOOInIgLcP+I3XJO2Cr/+Kj17iooAHx+FgRERKXT+vGy1\nWFQE+Ppe+XeWadJWoUEDoG1baUxERGRX2dlAcPC1SV8PwyZ+gHV+IiJ31/cBgyd+1vmJyO7cXd8H\nDJ74OaWTiOyOiZ+IyGbS091f6jHsrB4AOH5c/odPnnTfUmUiIrPQNKBpUyA3F2je/Nq/t9ysHgBo\n1UqeZB89qjoSIiLvO3xYEn9lSV8PQyd+ALjlFmD3btVREBF53549kgPdjYmfiMigdu8GIiLc/7pM\n/EREBrVnDxM/EZGteCrxG3pWDyC7cAUEAD//DNQx/I8pIiL3KC+XB7tHjsivlbHkrB5AmhO1bAkc\nPKg6EiIi7zl4UHJfVUlfD8MnfoDlHiKyH0+VeQAmfiIiQ/LUjB7ARIl/zx7VURAReY+n5vADJkr8\nHPETkZ0YutSzevVqhIeHIywsDJMmTbrm751OJ5o1a4aoqChERUVhwoQJtT5H585AVhZQUqI3WiIi\n4ysrkz78nTt75vV1bWpYVlaG0aNHY/369QgICEDPnj0RFxeHzldFe88992DZsmUun6dBAyAwUHai\n8dSFICIyipwcwN8faNzYM6+va8SfkpKC0NBQBAcHw9fXF0OHDkViYuI1x7ljqUBEBMs9RGQPO3cC\nkZGee31diT8/Px9BQUEXvw4MDER+fv4VxzgcDmzevBmRkZGIjY3F3r17XToX6/xEZBe7dnk28esq\n9Thq0CS/W7duyM3NRcOGDbFq1SoMGTIEmZmZlR6bkJBw8ffR0dGIjo6++PUttwCff64nWiIic9i1\nC3jqqWv/3Ol0wul06n59XS0bkpOTkZCQgNWrVwMA3nvvPdSpUwfjxo2r8nvat2+P7du3o2XLllcG\ncp2lx3v2AA8+yM3Xicj62rcH1q4FwsKqP05Jy4YePXogKysLBw8eRHFxMZYsWYK4uLgrjikoKLgY\nWEpKCjRNuybp10RYmGxK8OuveiImIjK2n36S3Qc7dPDcOXSVenx8fDBjxgwMGDAAZWVlePbZZ9G5\nc2fMnj0bADBy5Eh8+eWXmDlzJnx8fNCwYUMsXrzYpXPVqweEhgL79gHduumJmojIuP77Xylt163r\nuXMYvjvn5X7/e6BPH2DECC8FRUTkZR99BOzYAcyZc/1jLdud83Jdu8oFISKyql27gC5dPHsOJn4i\nIgPxRuI3VannxAmp8586BdRgJikRkamUl8seJLm5QPPm1z/eFqWeG2+UJczclIWIrOjAAdl8pSZJ\nXw9TJX6A5R4isq6dOz1f5gFMmvh37lQdBRGR+6WlAVFRnj+PKRM/R/xEZEWpqd5Zp8TET0RkEBzx\nV6FDB6CwUGb2EBFZxbFjwPnzwE03ef5cpkv8derIww/W+YnIStLSpMzjjanqunr1qNK1q1yky7o2\nE+miaUBBAbB3r0ypO3oUOHlStvssK5NpxE2bAu3ayafOzp2BNm1UR01WkprqnTIPYNLE37078PXX\nqqMgszt0CFixAnA6ge++kyQfEQGEhABt28p2n/XqSbOss2ela+J33wELFsimQH5+wO23A/37A/fd\nB9x8s+r/IzKz1FTgkUe8cy5TrdytsGsX8OijQHq6h4Miyzl1CvjkE+Df/5bVkYMGAX37Ar17S+Ku\n6cdsTZOFhJs2AWvWyH/t2wPDhgFPPimLDYlqo0MHYNUqoFOnmn+Pqyt3TZn4S0tlZduRI/Lxm+h6\n0tKADz4AvvoKiI0FnnsO+O1vAR83feYtLQW++Qb49FNg+XIZmIwdC/zmN+55fbK2U6fkoe5PP8lz\nzJqyRcuGCj4+wK23ypuZqDq7dsnObYMGSRLOygIWLpT23u5K+oC8Vv/+UgZKT5dnAX36AHFxsnsc\nUXV27JA9dmuT9PUwZeIHpM6/fbvqKMio8vKk7NK/P3D33UBODvDGG0CrVp4/t78/8Ne/SimoTx/g\n3nuBZ5+VmIgq462FWxVMm/h79AC2bVMdBRlNaSkwdarM/AoLA7KzgVdeARo08H4sfn7A//wPkJkJ\ntG4tMU2dKrOEiC63bZt3E78pa/yAbE/28MPcfJ0uSU2V3dlatwY+/BDo2FF1RFfKzgaef15mCM2d\nK+VKIkDazS9bVvtnQraq8QMyjzovD/j5Z9WRkGrl5cDkyTKl8vXXgbVrjZf0AXlzf/21JP8+fYCJ\nEyV2sreTJ4Eff6zdbB69TJv4fXzkYUhqqupISKX8fKBfP5mP/8MPMpXSyJv01KkjiT81VWKOjZU3\nPdnXtm3yzNKTm6tfzbSJH+ADXrv77jugZ08ZPX/7rbkWUAUFSczdusl/GzaojohU+eEHuY+9ydSJ\nv0cPuWhkL5oGzJwpz3jmzQP+8hfvjpbcxdcXePddqfc/+igwZ47qiEiFlBTgttu8e07TPtwFZL70\nwIHSW4XsoaQEGD0a+P57IDFR6uZWkJUF3H+/3M//+Ic5f5BR7WmatAfZutW1T6y2e7gLyAO806dZ\nI7WLs2eB3/1OHuonJ1sn6QMy9TQ5WRZ7DR4MnDmjOiLyhrw8Sf7eaMV8OVMn/jp15CPS1q2qIyFP\nO35cFkK1aSNtF5o0UR2R+7VoASQlyarffv1ktgdZW0V939sTEkyd+AGgVy8ZKZF1HTwI3HWXrMKd\nO1dq41bl6wt8/DFwzz3SSyg/X3VE5EkqHuwCFkj8t9/OEb+V5eRIEhw9GpgwwdhTNd3F4QAmTQKe\neupSuwmyppQUNYnf1A93AeDECemfXljIB2JWk5UlLZP//Gdg5EjV0agxezbwzjsy9TMkRHU05E5l\nZVLeO3AAuOEG117D1bxpyo1YLnfjjbJEPz1dNtEga8jIkKQ/frw0OLOrih94FWsVOnRQGw+5z+7d\nQECA60lfD9MnfuBSnZ+J3xoOHJCHmxMmAE8/rToa9UaOlNFhnz6yW1hwsOqIyB02bwbuuEPNuU1f\n4wdY57eSY8eAmBjgj39k0r/cSy8Br70myf/IEdXRkDts2QLceaeac1si8XNmjzWcOgUMGCAJf9Qo\n1dEYz+jRsnPYgAFyrcjcNm9Wl/hN/3AXAIqLgZYtZepbs2ZuDoy84uxZGenffrusXLXD7B1XaBrw\n6qvyCXfdOqBhQ9URkSsqunGePKlv1y1brtytUK+eTInaskV1JOSK4mLZHrFTJ2DKFCb96jgcco1C\nQqRXUUmJ6ojIFVu2yCDHW1stXs0SiR+Q+c7ff686CqotTZOHl35+snBJ1RvBTOrUkYVsDoeUxIzx\nmZ1qQ2WZB2DiJ8Xee092U1u40L2bn1udry+weLEsAJoyRXU0VFsqZ/QAFqnxA7ITV7t2spCrXj03\nBkYes3gxMG6cfOxt1051NOaUlyclgw8+AB56SHU0VBMVzySPHAGaNtX3Wrau8QNyAcPCuCOXWWze\nDMTHA8uXM+nrERgoe7W+8IKM/sn4tm2T51l6k74elkn8AMs9ZpGTI6PT+fOBLl1UR2N+3bpJzX/I\nEGloR8a2YYM04FOJiZ+8qqgIiIuTXbMGDlQdjXXExQFvvCG/nj2rOhqqzsaN6hO/ZWr8gNQ7u3aV\n3u2cEmg85eUy0m/VitsMeoKmASNGAL/+Ks9P+B4wntJS6c2TkyN9xvSyfY0fkHpnkybAvn2qI6HK\nvPMOUFAATJ+uOhJrcjiAWbOk19Hkyaqjocrs2AEEBbkn6ethqcQPANHR0siKjGXFCmkx/H//B9Sv\nrzoa6/LzA5YulVk+a9aojoautmGD7C+hmuUS/733SvtaMo6MDOCZZ4AvvpCNpcmzAgOBJUtkIxdu\n4mIsRqjvAxar8QNAbq7Mcigo4CpQIygqkn2RX31VGoyR93z0ETBzpjQwbNRIdTRUXi4lnj173DcA\nYo3/gqAgadS2Z4/qSEjTgOefB3r3ZtJX4cUXge7d5VdjDO/sbfduSfxG+NRrucQPsNxjFDNnSpln\n2jTVkdiTwyGj/rQ0medPan3zjTyDNAImfvKIbduAhASp6/v5qY7Gvho2lH+DN9+UGSWkzvr10nrc\nCCxX4weAo0dlG8YTJ1jnV+HUKSkxTJ4srYNJvUWLgLfflh/I3LPC+0pKpMzjrvn7FVjjv0zbtrIB\n+86dqiOxH02THbQGD2bSN5Jhw2S0+dxzrPerkJIChIaqn79fwZKJHwD69pWPVuRdU6fKvrl//7vq\nSOhqU6cC+/cDM2aojsR+1q+XnGQUlk38/fsDa9eqjsJeNm2S8s7nn7M1thH5+Um9/29/YydPb1u/\nHujXT3UUl1iyxg9If/6AAJnPz31JPa+wEIiKAj78ELj/ftXRUHW++gp4+WV52Mt6v+edOQO0aSP7\n7Lo7Fymr8a9evRrh4eEICwvDpEmTKj0mPj4eYWFhiIyMRFpamt5T1kjTprKQa8MGr5zO1irm6z/w\nAJO+GQwZAgwaJFteGmPYZ20bN8oiRiMNQHWN+MvKytCpUyesX78eAQEB6NmzJxYtWoTOnTtfPCYp\nKQkzZsxAUlIStm7dirFjxyI5OfnaQBwO9H+6P+Ifj8egmEGuhnSFp55biQ17p6F9x/Oo76jv1te2\nu5XrVmLawmk4r53H8fz6+CU3Huk7B7EPj0n8+qsko76xK7HvR/l35HvEvSreI3uyzqNB3fp4/y/u\nv7aujvh17XKakpKC0NBQBAcHAwCGDh2KxMTEKxL/smXLMHz4cABAr169cPr0aRQUFMDf3/+a11sb\nvBY5H0pzEb0XaOW6lfg2byzyBuTg8IU/c9dr293KdSsx9sOxyIm60AimPRC0NQfrN/LamkWDBsCL\nY1dizPSxKH/wUkMfvkfc44r3SLD82VgDXVtdpZ78/HwEBQVd/DowMBD5+fnXPSYvL6/K18yJysH0\nRfr79k5bOA15d1zZocpdr2130xZOu5T0L8jtxWtrNombpl2R9AG+R9ylsveIka6trhG/o4Y7PVz9\nUaTK77uw2jb9dDqcTieidaxvPq+dr/TPz5Wfc/k1SfDaWgP/HT3HU9fW6XTC6Ya+87oSf0BAAHJz\ncy9+nZubi8DAwGqPycvLQ0BAQOUveK/8En4oXFfSB4D6jsqLzX512D9AL15ba+C/o+d46tpGR0df\nkRvHjx/v0uvoKvX06NEDWVlZOHjwIIqLi7FkyRLExcVdcUxcXBwWLFgAAEhOTkbz5s0rre9XCEkN\nwZhhY/SEBQCIfzweIWkhHnltu3vyvnjUXcpra3aVvUduTuG/ozvEPx6PDtuN+x7RNeL38fHBjBkz\nMGDAAJSVleHZZ59F586dMXv2bADAyJEjERsbi6SkJISGhqJRo0aYN29ela834NAAjBk9xi0PPype\nY/qi6diy/RzCbvbD+DHueW07Ky8HPpk7CA93BU4fmo5z5efgV8fPbf9u5D2Xv0fOlZ9D3n4/ND43\nBgP78t9Rr0Exg5CSAsxYPB23djfee8SyC7guN3kycPCgtKglfaZMAf7zH1kf4aNr2EBGU1oqbYPj\n4oA33lAdjfm9+CIQEgK89prnzuFq3rRF4t+3T1o4HD4sPcrJNdu3AwMHynL/CzN4yWIOHQJ69pQ9\nkm+7TXU05lVeDtx0E7BuHXDZ7Ha3Y3fOaoSHS58S9iN33Zkz0uFx+nQmfSu7+Wb5ZDxsmLQ9Idds\n2wY0aeLZpK+HLRK/wyFtgpcvVx2JeY0ZA9x9N/DYY6ojIU97+GGgTx9g9GjVkZjXf/4jLUyMyhaJ\nH5C65bJlqqMwp8WLpfMmt1C0j/ffl5LeZ5+pjsScjJ74bVHjB+TBVbt2wNatQPv2HjuN5Rw4ILXe\n1atlVy2yj7Q0eTa2dSvQoYPqaMzDm88UWeO/Dh8f+Qn85ZeqIzGP0lLgiSeAceOY9O0oKgr405+k\n3l9Sojoa81i6VDqgGnkiiW0SPwA88ohsREE1M3480Lgx8MorqiMhVcaOBVq2BP76V9WRmIfRyzyA\njUo9gIxg27aV2iXLPdXbsAEYOlQ+7rdpozoaUqmgQEb/n34qD32pavv3A7ffDhw54p11Liz11ICP\nD/Dggyz3XE9hIfD73wNz5zLpE+DvD8ybBzz1FHDihOpojG3xYpkVZfTFjbZK/ADLPdejacBzzwEP\nPQTExqqOhoxiwACZyvvcc9y1qzoLFwKPP646iuuzXeKPjpb2DTk51zvSnubMkZk8EyeqjoSM5t13\nZabKrFmqIzGm//4XKCoC7rxTdSTXZ7vE7+MjtetPP1UdifHs3Qv85S/AokXgFop0jfr15d54+21g\n927V0RjPwoWSW+qYIKuaIET3e+opYMECfmS93LlzctO+9560uCCqTKdOwKRJMsXz119VR2McmiY/\nFM1Q5gFsmvi7d5fePZs3q47EON54Q97Uzz6rOhIyuhEjgN/8Bnj9ddWRGMemTUDDhkCXLqojqRlb\nJn6H49Kon6QTY2Ki1PeNvOiEjMHhAGbPlvuGbVDEv/4lPxDN8v6x1Tz+y+XlyU/n/HygQQOvndZw\njh6VOdpffilN2IhqatMmmf21fTtQ1W6qdvDzz9LVND1dpr56E+fx11JgoJR8EhNVR6JOWZnM13/h\nBSZ9qr277gJeekk+PZeXq45GnSVLZLagt5O+HrZN/ADwhz/Ye2ra3/4mb9i33lIdCZnVn/8sfXwm\nT1YdiTpz58r6BjOxbakHkBv25puB9evlYZWdrFsHDB8OpKZydS7pc/iwdHC1Y7lwzx7pxHnokJrV\nuiz1uMDXV2ax2G3Uf+SIfDz/7DMmfdLvppukpcPQodLXx05mz5aHukZv0XA1W4/4ARmtREXJr40a\nef30XldaKo22+veXxVpE7vLWWzJFeu1aoG5d1dF43k8/SbPHXbvkmaEKHPG76Kab5CHVokWqI/GO\nt96S+cZ/+pPqSMhqEhKu/NXq5s2TAZSqpK+H7Uf8gNS7X35Zem2YZR6uK1aulBk8qalAq1aqoyEr\n+vFHmS03Zw4wcKDqaDynrAzo2FFav9xxh7o4OOLXoV8/qfevWqU6Es85dAh45hn5ZMOkT57SurXc\nY08/LfecVSUlyQY1t9+uOhLXMPFDRvmvvWbdKWm//CI7Ar3xhv1mXZD33X233GuPPAKcP686Gs+Y\nOlV2JzNrhYClngtKSoDQUJmS1rOnsjDcTtOAJ5+UG/STT8x7o5K5aJqs6r3xRpn5YqX7btMmWfiY\nkSGVApVY6tHJ11fq/FYb9f/jH7KU/OOPrfXmI2NzOID582WWz0cfqY7GvSZMAP74R/VJXw+O+C9z\n5oyM+teuNU+XveqsWSO11q1bZfYSkbft3y8bk3z2GdC3r+po9PvhB9m+NTvbGHtWcMTvBo0by09y\nK7QwyM6WRVpLljDpkzodOlzqU5+drToa/SZMkOcXRkj6enDEf5Vz54CwMKn19+qlOhrXnD4taxNG\njZImWkSqzZwJTJ8OJCcDTZuqjsY1W7YAjz4KZGYap6Ovq3mTib8Sc+bIhuzr1qmOpPaKi2X+dEQE\n8MEHrOuTcbz0kkzxTEw0X4sDTZPZSn/4g/S4MgqWetxoxAi5QZOSVEdSO5oGPP+8lKz++U8mfTKW\nDz6QliGjRplv29OlS4GzZ2WGnBUw8VfC1xeYNk3m6ZppHvL48bJh+sKF9uiVQubi6ysl1G3bpCW4\nWRQXA+PGAVOmWOd9xcRfhfvuk1bNU6eqjqRm5syRrSRXrLBHszkypyZNpHXIv/8N/O//qo6mZv7+\nd6BzZ1nhbxWs8Vdj/35ZzJWWZuyZMZ99JiOSDRuAkBDV0RBdX2YmcM89Msf/gQdUR1O1jAyZKJGa\naswcwBq/B3ToALzyivTsN+rWcomJwKuvypx9Jn0yi44dLzUNXLFCdTSVKy+Xh7lvv23MpK8HE/91\njBsHFBXJdDSjWbtWHuauWCGzeIjMpFs3YPlyaR64Zo3qaK41a5ZM7x41SnUk7sdSTw1kZsrHvU2b\nZKRiBMuWyT6fS5ey8RqZ2+bNwO9+Jwu9jFJH37lTYjHSe74yLPV4UMeOMgvhkUekrYNqixbJR9Ck\nJCZ9Mr8775QBzOOPy6wf1YqKZKHW++8bO+nrwRF/DWma1PqLioDPP1c3R/7jj2WHozVrgFtuURMD\nkSekpQH33y8tU154QU0M5eXAsGGyFmbuXDUx1AZX7nrB+fMyEyE2Vh74eFN5ueyRu2SJJP3QUO+e\nn8gbcnKAAQOAJ56QAY63B1ivvy6lp/XrjdOWoTqu5k2TLZxWq359+UjauzfQvDkQH++d8545I102\nCwqk0+aNN3rnvETeFhIidfW4OGDfPuBf/5LRtzdMmyYPmzdtMkfS14M1/lpq1w74+mtZ2DVrlufP\nt3u3rCVo3lxGIUz6ZHX+/rImpXFj2c82Pd3z55w6VVbmrl4N3HCD58+nGhO/C4KDJflPnCglH0/M\n8S8vBz78ELj3XuDNN2WVo9lbwRLVlJ+f1NjHjJEJDDNneqa/j6ZJCXXOHOD77+W9bQes8etQUCCb\nMrRpA8yb5752sxkZMj+/tFQ+6oaHu+d1icwoI0OaozVtCsyYIe0T3OHECVlD8OOPUuJp1co9r+tN\nnM6pgL8/8M038mtEBPDVV/per6BAFovcdZfsV/rdd0z6RJ06SS/8uDh5vvbyy8CxY66/nqbJiveo\nKHl/bdxozqSvBxO/TvXrS7+RTz+VVb4xMdLHvzY/hHftknn54eHyehkZ0hnUKp0AifTy8ZH3xO7d\nQFmZNFAcPRrYsaPmr1FeLgO13/5Wyjvz5ske2/XqeS5uo2Kpx42Ki6Ul8pQpQEkJMHiw/CDo1AkI\nDJREXlIC5ObKDbxpk3zELCoCRo6U5O/vr/r/gsj4jh2TZ2ALFgDNmkk33bvvlrUtbdvKM4LiYuD4\ncWmwtmULsHixPDB+5RXZltQKAyvO4zcQTZOe46tWAd9+K10+jxyRkYqPj8wMiogAuneXHw7duwN1\n+NmLqNbKy2UA5XRKaTQzEzh6VN5rdesCLVoAXbvKzLiHHgIiI621QRETv8GVlckNxwRP5FmaJhMj\nfH1VR+J5TPxERDbDWT1ERFQjTPxERDbjcq+ewsJCPPbYYzh06BCCg4Px+eefo3nz5tccFxwcjKZN\nm6Ju3brw9fVFSkqKroCJiEgfl0f8EydORExMDDIzM9G3b19MnDix0uMcDgecTifS0tKY9ImIDMDl\nxL9s2TIMHz4cADB8+HB8Vc2yVT60JSIyDpcTf0FBAfwvrDby9/dHQUFBpcc5HA7069cPPXr0wMcf\nf+zq6YiIyE2qrfHHxMTgWCVNMd55550rvnY4HHBUsSpi06ZNaNu2LY4fP46YmBiEh4ejd+/elR6b\nkJBw8ffR0dGIjo6+TvhERPbhdDrhdDp1v47L8/jDw8PhdDrRpk0bHD16FPfeey/Sr9M4e/z48Wjc\nuDFeffUhaTnmAAAEc0lEQVTVawPhPH4iolrx+jz+uLg4zJ8/HwAwf/58DBky5JpjfvnlFxQVFQEA\nzp49i7Vr1+LWW2919ZREROQGLo/4CwsL8eijj+Lw4cNXTOc8cuQInn/+eaxcuRL79+/Hgw8+CAAo\nLS3FE088gTfffLPyQDjiJyKqFbZsICKyGbZsICKiGmHiJyKyGSZ+IiKbYeInIrIZJn4LcscCD7qE\n19O9eD3VY+K3IL6x3IvX0714PdVj4icishkmfiIimzHMAq6uXbti586dqsMgIjKNyMhI7Nixo9bf\nZ5jET0RE3sFSDxGRzTDxExHZjJLE/8UXXyAiIgJ169ZFampqlcetXr0a4eHhCAsLw6RJk7wYobkU\nFhYiJiYGHTt2RP/+/XH69OlKjwsODkaXLl0QFRWF2267zctRGl9N7rf4+HiEhYUhMjISaWlpXo7Q\nPK53LZ1OJ5o1a4aoqChERUVhwoQJCqI0h2eeeQb+/v7VtrSv9X2pKbBv3z4tIyNDi46O1rZv317p\nMaWlpVpISIh24MABrbi4WIuMjNT27t3r5UjN4fXXX9cmTZqkaZqmTZw4URs3blylxwUHB2snT570\nZmimUZP7beXKldrAgQM1TdO05ORkrVevXipCNbyaXMtvv/1WGzx4sKIIzWXjxo1aamqqdsstt1T6\n967cl0pG/OHh4ejYsWO1x6SkpCA0NBTBwcHw9fXF0KFDkZiY6KUIzYUb3+tXk/vt8uvcq1cvnD59\nusq9pu2spu9d3os107t3b7Ro0aLKv3flvjRsjT8/Px9BQUEXvw4MDER+fr7CiIyLG9/rV5P7rbJj\n8vLyvBajWdTkWjocDmzevBmRkZGIjY3F3r17vR2mZbhyX1a72boeVW3U/u6772Lw4MHX/f6qNm+3\nK29vfG83Nb3frh6l8j69Vk2uSbdu3ZCbm4uGDRti1apVGDJkCDIzM70QnTXV9r70WOJft26dru8P\nCAhAbm7uxa9zc3MRGBioNyzTqu56+vv749ixYxc3vm/dunWlx7Vt2xYA0KpVKzzwwANISUlh4r+g\nJvfb1cfk5eUhICDAazGaRU2uZZMmTS7+fuDAgXjppZdQWFiIli1bei1Oq3DlvlRe6qmqztejRw9k\nZWXh4MGDKC4uxpIlSxAXF+fl6MyBG9/rV5P7LS4uDgsWLAAAJCcno3nz5hdLbHRJTa5lQUHBxfd+\nSkoKNE1j0neRS/ele547187SpUu1wMBAzc/PT/P399fuu+8+TdM0LT8/X4uNjb14XFJSktaxY0ct\nJCREe/fdd1WEagonT57U+vbtq4WFhWkxMTHaqVOnNE278nrm5ORokZGRWmRkpBYREcHrWYnK7rdZ\ns2Zps2bNunjMqFGjtJCQEK1Lly5Vzkij61/LGTNmaBEREVpkZKR2xx13aFu2bFEZrqENHTpUa9u2\nrebr66sFBgZqc+fO1X1fsmUDEZHNKC/1EBGRdzHxExHZDBM/EZHNMPETEdkMEz8Rkc0w8RMR2QwT\nPxGRzTDxExHZzP8DJ60Em2UYStAAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x38fae10>"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Use method of undetermined coefficients to find\n",
      "# Newton-Cotes rule for Gauss nodes.\n",
      "\n",
      "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)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for i in range(2*len(nodes) + 1):\n",
      "    approx = np.dot(weights, nodes**i)\n",
      "    true = 1/(i+1)*(1. - (-1)**(i+1))\n",
      "    \n",
      "    print \"Error at degree %d: %g\" % (i, approx-true)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Error at degree 0: 0\n",
        "Error at degree 1: -1.11022e-16\n",
        "Error at degree 2: -1.11022e-16\n",
        "Error at degree 3: -2.498e-16\n",
        "Error at degree 4: 0\n",
        "Error at degree 5: -4.71845e-16\n",
        "Error at degree 6: 2.22045e-16\n",
        "Error at degree 7: -5.68989e-16\n",
        "Error at degree 8: -0.01161\n",
        "Error at degree 9: -6.66134e-16\n",
        "Error at degree 10: -0.0257832\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}