{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1769d872",
   "metadata": {},
   "source": [
    "Quiz 8, MCS 471, Friday 21 October 2022, due 10:50am"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "828b9eda",
   "metadata": {},
   "source": [
    "Consider the quadrature rule\n",
    "\n",
    "$$\n",
    "   \\int_{-1}^{+1} f(x) dx \\approx w_0 f(x_0) + w_1 f(x_1) + w_{2} f(x_2).\n",
    "$$\n",
    "\n",
    "1. Write the conditions on the weights and the points of his rule\n",
    "   so the degree of precision is as high as possible.\n",
    "\n",
    "2. What is the degree of precision of this rule?  Justify your answer."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7935042b",
   "metadata": {},
   "source": [
    "# Answer to Question 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4974286",
   "metadata": {},
   "source": [
    "We have three unknown weights and three unknown interpolation points, so we may enforce six conditions, requiring the first six basis functions to be exact."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "deaccc23",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(w0, w1, w2, x0, x1, x2)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using SymPy\n",
    "w0, w1, w2, x0, x1, x2 = Sym(\"w0, w1, w2, x0, x1, x2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "35ea849c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6-element Vector{var\"#b0#3\"{Int64}}:\n",
       " (::var\"#b0#3\"{Int64}) (generic function with 1 method)\n",
       " (::var\"#b0#3\"{Int64}) (generic function with 1 method)\n",
       " (::var\"#b0#3\"{Int64}) (generic function with 1 method)\n",
       " (::var\"#b0#3\"{Int64}) (generic function with 1 method)\n",
       " (::var\"#b0#3\"{Int64}) (generic function with 1 method)\n",
       " (::var\"#b0#3\"{Int64}) (generic function with 1 method)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B = [b0(x) = x^k for k=0:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "21d4756b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\left[ \\begin{array}{r}1\\\\\\frac{1}{2}\\\\\\frac{1}{3}\\\\\\frac{1}{4}\\\\\\frac{1}{5}\\\\\\frac{1}{6}\\end{array} \\right]$\n"
      ],
      "text/plain": [
       "6-element Vector{Sym}:\n",
       "   1\n",
       " 1/2\n",
       " 1/3\n",
       " 1/4\n",
       " 1/5\n",
       " 1/6"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z = Sym(\"z\")\n",
    "r = [integrate(B[k](z), (z, 0, 1)) for k=1:6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "422cbaae",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\begin{align*}f\\end{align*}$\n"
      ],
      "text/plain": [
       "f"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f = SymFunction(\"f\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "dcbefd8c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "rule (generic function with 1 method)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rule(f) = w0*f(x0) + w1*f(x1) + w2*f(x2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6517b73f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\left[ \\begin{array}{r}w_{0} + w_{1} + w_{2}\\\\w_{0} x_{0} + w_{1} x_{1} + w_{2} x_{2}\\\\w_{0} x_{0}^{2} + w_{1} x_{1}^{2} + w_{2} x_{2}^{2}\\\\w_{0} x_{0}^{3} + w_{1} x_{1}^{3} + w_{2} x_{2}^{3}\\\\w_{0} x_{0}^{4} + w_{1} x_{1}^{4} + w_{2} x_{2}^{4}\\\\w_{0} x_{0}^{5} + w_{1} x_{1}^{5} + w_{2} x_{2}^{5}\\end{array} \\right]$\n"
      ],
      "text/plain": [
       "6-element Vector{Sym}:\n",
       "                w0 + w1 + w2\n",
       "       w0*x0 + w1*x1 + w2*x2\n",
       " w0*x0^2 + w1*x1^2 + w2*x2^2\n",
       " w0*x0^3 + w1*x1^3 + w2*x2^3\n",
       " w0*x0^4 + w1*x1^4 + w2*x2^4\n",
       " w0*x0^5 + w1*x1^5 + w2*x2^5"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R = [rule(B[k]) for k=1:6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "7394623d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\left[ \\begin{array}{r}w_{0} + w_{1} + w_{2} - 1\\\\w_{0} x_{0} + w_{1} x_{1} + w_{2} x_{2} - \\frac{1}{2}\\\\w_{0} x_{0}^{2} + w_{1} x_{1}^{2} + w_{2} x_{2}^{2} - \\frac{1}{3}\\\\w_{0} x_{0}^{3} + w_{1} x_{1}^{3} + w_{2} x_{2}^{3} - \\frac{1}{4}\\\\w_{0} x_{0}^{4} + w_{1} x_{1}^{4} + w_{2} x_{2}^{4} - \\frac{1}{5}\\\\w_{0} x_{0}^{5} + w_{1} x_{1}^{5} + w_{2} x_{2}^{5} - \\frac{1}{6}\\end{array} \\right]$\n"
      ],
      "text/plain": [
       "6-element Vector{Sym}:\n",
       "                  w0 + w1 + w2 - 1\n",
       "       w0*x0 + w1*x1 + w2*x2 - 1/2\n",
       " w0*x0^2 + w1*x1^2 + w2*x2^2 - 1/3\n",
       " w0*x0^3 + w1*x1^3 + w2*x2^3 - 1/4\n",
       " w0*x0^4 + w1*x1^4 + w2*x2^4 - 1/5\n",
       " w0*x0^5 + w1*x1^5 + w2*x2^5 - 1/6"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eqs = [R[k] - r[k] for k=1:6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e1c88046",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6-element Vector{NTuple{6, Sym}}:\n",
       " (5/18, 5/18, 4/9, 1/2 - sqrt(15)/10, sqrt(15)/10 + 1/2, 1/2)\n",
       " (5/18, 5/18, 4/9, sqrt(15)/10 + 1/2, 1/2 - sqrt(15)/10, 1/2)\n",
       " (5/18, 4/9, 5/18, 1/2 - sqrt(15)/10, 1/2, sqrt(15)/10 + 1/2)\n",
       " (5/18, 4/9, 5/18, sqrt(15)/10 + 1/2, 1/2, 1/2 - sqrt(15)/10)\n",
       " (4/9, 5/18, 5/18, 1/2, 1/2 - sqrt(15)/10, sqrt(15)/10 + 1/2)\n",
       " (4/9, 5/18, 5/18, 1/2, sqrt(15)/10 + 1/2, 1/2 - sqrt(15)/10)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve(eqs, (w0, w1, w2, x0, x1, x2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1b4eb7b",
   "metadata": {},
   "source": [
    "The above solution applied ``SymPy`` and is computational.  Alternatively, a solution that does not involve computations, but that is also correct is below."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb72410d",
   "metadata": {},
   "source": [
    "Because we have six unknowns, we can impose six conditions, requiring that $f = 1, x, x^2, x^3, x^4, x^5$ are integrated correctly.\n",
    "\n",
    "$$\n",
    "   \\int_{-1}^{+1} 1 dx = w_0 + w_1 + w_2\n",
    "$$\n",
    "\n",
    "$$\n",
    "   \\int_{-1}^{+1} x dx = w_0 x_0 + w_1 x_1 + w_2 x_2\n",
    "$$\n",
    "\n",
    "$$\n",
    "   \\int_{-1}^{+1} x^2 dx = w_0 x_0^2 + w_1 x_1^2 + w_2 x_2^2\n",
    "$$\n",
    "\n",
    "$$\n",
    "   \\int_{-1}^{+1} x^3 dx = w_0 x_0^3 + w_1 x_1^3 + w_2 x_2^3\n",
    "$$\n",
    "\n",
    "$$\n",
    "   \\int_{-1}^{+1} x^4 dx = w_0 x_0^4 + w_1 x_1^4 + w_2 x_2^4\n",
    "$$\n",
    "\n",
    "$$\n",
    "   \\int_{-1}^{+1} x^5 dx = w_0 x_0^5 + w_1 x_1^5 + w_2 x_2^5\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2885b6a9",
   "metadata": {},
   "source": [
    "# Answer to Question 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3815b287",
   "metadata": {},
   "source": [
    "The degree of precision is 5 because the first 5 basis functions are integrated correctly and the integral is a linear operator."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.8.0",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}