{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "#  Computing the Weights in Newton-Cotes Rules"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import numpy.linalg as la"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We start by choosing our *quadrature nodes*, the maximum degree which will be exact, as well as the interval $(a,b)$ on which we integrate:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nodes = [0, 1]\n",
        "#nodes = [0, 0.5, 1]\n",
        "#nodes = [3, 3.5, 4]\n",
        "#nodes = [0, 1, 2]\n",
        "#nodes = np.linspace(0, 1, 12)\n",
        "\n",
        "max_degree = len(nodes)-1\n",
        "\n",
        "a = nodes[0]\n",
        "b = nodes[-1]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we compute the transpose of the Vandermonde matrix $V^T$ and the integrals $\\int_a^b x^i$ as `rhs`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[[1 1]\n",
            " [0 1]]\n"
          ]
        }
      ],
      "source": [
        "nodes = np.array(nodes)\n",
        "powers = np.arange(max_degree+1)\n",
        "\n",
        "Vt = nodes ** powers.reshape(-1, 1)\n",
        "\n",
        "rhs = 1/(powers+1) * (b**(powers+1) - a**(powers+1))\n",
        "\n",
        "if len(nodes) <= 4:\n",
        "    print(Vt)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Set up the linear system for the weights:\n",
        "\n",
        "$$\n",
        "\\begin{align*}\n",
        "\\alpha_0 x_0^0 + \\cdots + \\alpha_{n-1} x_{n-1}^{0} &= \\int_a^b x^0\\\\\n",
        "\\vdots &= \\vdots \\\\\n",
        "\\alpha_0 x_0^{n-1} + \\cdots + \\alpha_{n-1} x_{n-1}^{n-1} &= \\int_a^b x^{n-1}\n",
        "\\end{align*}\n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[ 0.5  0.5]\n"
          ]
        }
      ],
      "source": [
        "weights = la.solve(Vt, rhs)\n",
        "\n",
        "print(weights)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we test our quadrature rule by integrating the monomials $\\int_a^b x^i dx$ and comparing quadrature results to the true answers:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Error at degree 0: 0\n",
            "Error at degree 1: 0\n",
            "Error at degree 2: 0.166667\n"
          ]
        }
      ],
      "source": [
        "for i in range(len(nodes) + 1):\n",
        "    approx = weights @ nodes**i\n",
        "    true = 1/(i+1)*(b**(i+1) - a**(i+1))\n",
        "    \n",
        "    print(\"Error at degree %d: %g\" % (i, approx-true))\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {
        "collapsed": false
      },
      "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
}