{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "547b4273",
   "metadata": {},
   "source": [
    "# Gauss quadrature rules"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ae3f9e1",
   "metadata": {},
   "source": [
    "Evaluating the integrand at the $n$ roots of an orthogonal polynomial of degree $n$, the algebraic degree of precision of the Gauss quadrature rule is $2n-1$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc44d2c5",
   "metadata": {},
   "source": [
    "## 1. Gauss-Legendre quadrature"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce8750a1",
   "metadata": {},
   "source": [
    "We use SymPy to compute the expanded form of the Legendre polynomials.\n",
    "The Legendre polynomial are defined by the 3-terms recursion:\n",
    "\n",
    "$$\n",
    "   L_0(x) = 1, ~~ L_1(x) = x, ~~\n",
    "   (n+1) L_{n+1}(x) - (2n+1) x L_n(x) + n L_{n-1}(x) = 0.\n",
    "$$\n",
    "\n",
    "The Legendre polynomials are defined over $[-1,+1]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ae4a63db",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$x$"
      ],
      "text/plain": [
       "x"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using SymPy\n",
    "x = Sym(\"x\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5041f835",
   "metadata": {},
   "source": [
    "The variable ``x`` is a global variable to represent the symbol ``x``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d2873741",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "legendre"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    legendre(d::Int)\n",
    "\n",
    "returns the Legendre polynomial of degree d,\n",
    "as a SymPy expression.\n",
    "\"\"\"\n",
    "function legendre(d::Int)\n",
    "    if d == 0\n",
    "        return 1\n",
    "    elseif d == 1\n",
    "        return x\n",
    "    end\n",
    "    L0 = 1\n",
    "    L1 = Sym(\"x\")\n",
    "    L2 = Sym(\"0\")\n",
    "    for n = 2:d\n",
    "        L2 = expand(((2*n - 1)*x*L1 - (n-1)*L0)/n)\n",
    "        (L0, L1) = (L1, L2)\n",
    "    end\n",
    "    return L2\n",
    "end"
   ]
  },
  {
   "cell_type": "raw",
   "id": "ab7b1539",
   "metadata": {},
   "source": [
    "Let us check the first six Legendre polynomials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "98f097dd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6-element Vector{Number}:\n",
       "                            1\n",
       "                            x\n",
       "                3*x^2/2 - 1/2\n",
       "              5*x^3/2 - 3*x/2\n",
       "    35*x^4/8 - 15*x^2/4 + 3/8\n",
       " 63*x^5/8 - 35*x^3/4 + 15*x/8"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L6 = [legendre(i-1) for i=1:6]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a8b3e1a",
   "metadata": {},
   "source": [
    "To compute the roots of the polynomial, we extract the coefficient vectors with array comprehensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f5adc1b5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\left[ \\begin{array}{r}0\\\\\\frac{15}{8}\\\\0\\\\- \\frac{35}{4}\\\\0\\\\\\frac{63}{8}\\end{array} \\right]$\n"
      ],
      "text/plain": [
       "6-element Vector{Sym}:\n",
       "     0\n",
       "  15/8\n",
       "     0\n",
       " -35/4\n",
       "     0\n",
       "  63/8"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = legendre(5)\n",
    "Cp = [p.coeff(x,k) for k=0:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e55474e6",
   "metadata": {},
   "source": [
    "We need floating-point numbers, we convert to ``Float64`` type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fe6ef972",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6-element Vector{Float64}:\n",
       "  0.0\n",
       "  1.875\n",
       "  0.0\n",
       " -8.75\n",
       "  0.0\n",
       "  7.875"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C5 = [Float64(c) for c in Cp]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22d80eaf",
   "metadata": {},
   "source": [
    "The roots of the polynomial are computed via the eigenvalues of the companion matrix of the polynomial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "5278cf6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "0fe3490c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "rootsCompanion"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    rootsCompanion(cff::Array{Float64,1})\n",
    "\n",
    "returns the roots of the polynomial with coefficients cff,\n",
    "by computing the eigenvalues of the companion matrix.\n",
    "The last coefficient should not be zero.\n",
    "\"\"\"\n",
    "function rootsCompanion(cff::Array{Float64,1})\n",
    "    lead = cff[end] # leading coefficient\n",
    "    dim = length(cff) - 1\n",
    "    companion = zeros(dim, dim)\n",
    "    for k = 1:dim-1\n",
    "        companion[k+1, k] = 1\n",
    "    end\n",
    "    for k = 1:dim\n",
    "        companion[k, dim] = -cff[k]/lead\n",
    "    end\n",
    "    return eigvals(companion)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3af3dc30",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Vector{Float64}:\n",
       " -0.9061798459386623\n",
       " -0.5384693101056837\n",
       "  0.0\n",
       "  0.5384693101056831\n",
       "  0.9061798459386653"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rootsL5 = rootsCompanion(C5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "784ea81f",
   "metadata": {},
   "source": [
    "Let us check the residuals of the roots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e12f587a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-9.0617984593866230e-01 : 1.14e-14\n",
      "-5.3846931010568366e-01 : 1.42e-15\n",
      " 0.0000000000000000e+00 : 0.00e+00\n",
      " 5.3846931010568311e-01 : 9.83e-18\n",
      " 9.0617984593866530e-01 : 8.53e-15\n"
     ]
    }
   ],
   "source": [
    "using Printf\n",
    "for r in rootsL5\n",
    "    sr = @sprintf(\"%23.16e\", r)\n",
    "    err = evalpoly(r, C5)\n",
    "    serr = @sprintf(\"%.2e\", err)\n",
    "    println(\"$sr : $serr\")\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b614cc6d",
   "metadata": {},
   "source": [
    "The weights are the integrals of the Lagrange polynomials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "dfc39316",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "weights"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    weights(pts::Array{Float64,1})\n",
    "\n",
    "returns the weights of the Gauss-Legendre rule\n",
    "at the given points in pts,\n",
    "using the Lagrange polynomials.\n",
    "\"\"\"\n",
    "function weights(pts::Array{Float64,1})\n",
    "    dim = length(pts)\n",
    "    result = zeros(dim)\n",
    "    for i=1:dim\n",
    "        L = 1.0\n",
    "        for j=1:dim\n",
    "            if j != i\n",
    "                L = L*(x - pts[j])/(pts[i] - pts[j])\n",
    "            end\n",
    "        end\n",
    "        result[i] = integrate(L, (x, -1.0, 1.0))\n",
    "    end\n",
    "    return result\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "0e701b77",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Vector{Float64}:\n",
       " 0.2369268850561907\n",
       " 0.47862867049936353\n",
       " 0.5688888888888899\n",
       " 0.4786286704993681\n",
       " 0.23692688505618775"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "W5 = weights(rootsL5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2d56b50",
   "metadata": {},
   "source": [
    "The Legendre quadrature rule with 5 points can be defined in one line:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e8f38d19",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Legendre5rule (generic function with 1 method)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Legendre5rule(f) = sum(W5[i]*f(rootsL5[i]) for i=1:5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25a023ab",
   "metadata": {},
   "source": [
    "Let us test this on some polynomial of a degree 9."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "95876f99",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "some_p (generic function with 1 method)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "some_p(z) = z^9 + 12.3*z^3 - z + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "be5731d2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0000000000000018"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "approximation = Legendre5rule(some_p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "455426df",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$2.0$"
      ],
      "text/plain": [
       "2.00000000000000"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "exact = integrate(some_p(x), (x,-1,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c54df29",
   "metadata": {},
   "source": [
    "Let us check the error:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "7ca03493",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$2.22044604925031 \\cdot 10^{-15}$"
      ],
      "text/plain": [
       "2.22044604925031e-15"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "approximation - exact"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48f178e9",
   "metadata": {},
   "source": [
    "The error is close to the machine precision."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f86f0489",
   "metadata": {},
   "source": [
    "## 2. Reduction to an Eigenvalue Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84fb7b12",
   "metadata": {},
   "source": [
    "Following a 1969 paper of Golub and Welsch,\n",
    "published in Mathematics of Computation (volume 23, pages 221-230),\n",
    "the construction of a Gauss quadrature is reduced to an eigenvalue problem\n",
    "of a tridiagonal matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98ca98b2",
   "metadata": {},
   "source": [
    "If $p_n$ is an orthogonal polynomial of degree $n$,\n",
    "with the three terms recursion denoted as\n",
    "$$ p_{-1}(x) = 0, \\quad p_0(x) = 1, \\quad \\mbox{for } j>1:\n",
    "   p_j(x) = (a_j x + b_j) p_{j-1}(x) - c_j p_{j-2}(x),$$\n",
    "then the roots of $p_n$ are the eigenvalues of"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a71f1a93",
   "metadata": {},
   "source": [
    "$$ J =\n",
    "   \\left[\n",
    "      \\begin{array}{cccccc}\n",
    "         \\alpha_1 & \\beta_1 & & & & \\\\\n",
    "         \\beta_1 &\\alpha_2 & \\beta_2 & & & \\\\\n",
    "          & \\ddots &\\ddots & \\ddots & & \\\\\n",
    "          & & \\beta_{n-2} & \\alpha_{n-1} & \\beta_{n-1} \\\\\n",
    "          & & & \\beta_{n-1} & \\alpha_{n} \\\\\n",
    "      \\end{array}\n",
    "   \\right],\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05f09eda",
   "metadata": {},
   "source": [
    "with $\\displaystyle \\alpha_i = -\\frac{b_i}{a_i}$, $i=1,2,\\ldots,n$,\n",
    "     $\\displaystyle \\beta_i = \\sqrt{\\frac{c_{i+1}}{a_i a_{i+1}}}$,\n",
    "     $i = 1,2,\\ldots,n-1$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49abec80",
   "metadata": {},
   "source": [
    "If $\\bf q$ is the first row of $Q$, of the orthogonal matrix\n",
    "with the eigenvectors of $J$ in its columns, then\n",
    "$$\n",
    "   w_i = q_i^2 \\times \\int_{a}^{b} w(x) dx\n",
    "$$\n",
    "is the weight of the $i$-th point in the Gauss quadrature rule\n",
    "with weight function $w(x)$ over the interval $[a,b]$, as in\n",
    "$$\n",
    "   \\int_{a}^b w(x) f(x) dx \\approx \\sum_{i=1}^n w_i f(x_i),\n",
    "   \\quad \\mbox{with } p_n(x_i) = 0, i=1,2,\\ldots,n.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3a9ebc5",
   "metadata": {},
   "source": [
    "## 3. Application to Gauss-Legendre Quadrature"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bef9226",
   "metadata": {},
   "source": [
    "For Gauss-Legendre quadrature\n",
    "the weight function $w(x)$ is $1$ and $[a,b]$ is $[-1,+1]$.\n",
    "Following the notation of the previous section,\n",
    "the coefficients of the three terms recursion are\n",
    "$a_j = (2j-1)/j$, $b_j = 0$, and $c_j = (j-1)/j$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6dc3777",
   "metadata": {},
   "source": [
    "The ``eigvals`` and ``eigvecs`` of the ``LinearAlgebra`` module\n",
    "are used in the code below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "47461603",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "alphabetaLegendre"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Printf\n",
    "using LinearAlgebra\n",
    "\n",
    "\"\"\"\n",
    "    alphabetaLegendre(degree) computes the alpha\n",
    "and the beta coefficients for the tridiagonal matrix\n",
    "of the Legendre polynomials.\n",
    "\"\"\"\n",
    "function alphabetaLegendre(degree)\n",
    "    A = [(2*n-1)/Float64(n) for n=1:degree]\n",
    "    B = [0 for n=1:degree]\n",
    "    alphas = [-B[k]/A[k] for k=1:degree]\n",
    "    C = [(n-1)/Float64(n) for n=1:degree]\n",
    "    betas = [sqrt(C[k+1]/(A[k]*A[k+1])) for k=1:degree-1]\n",
    "    return alphas, betas\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "eefe9ddb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "makeGaussLegendre"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    makeGaussLegendre(degree::Int)\n",
    "\n",
    "computes the points and weights of the\n",
    "Gauss-Legendre quadrature rule with degree number of points\n",
    "via the eigenvalues and eigenvectors\n",
    "of a tridiagonal Jacobi matrix.\n",
    "\"\"\"\n",
    "function makeGaussLegendre(degree::Int)\n",
    "    a, b = alphabetaLegendre(degree)\n",
    "    J = diagm(a) + diagm(-1 => b) + diagm(+1 => b)\n",
    "    println(\"The triadiagonal matrix J :\")\n",
    "    for k=1:degree\n",
    "        println(J[k,:])\n",
    "    end\n",
    "    pts = eigvals(J)\n",
    "    println(\"The points :\")\n",
    "    for i=1:degree\n",
    "        spts = @sprintf(\"%23.16e\", pts[i])\n",
    "        println(\"x[\", i, \"] : $spts\")\n",
    "    end\n",
    "    Q = eigvecs(J)\n",
    "    weights = 2*(Q[1,:].^2)\n",
    "    println(\"The weights :\")\n",
    "    for i=1:degree\n",
    "        sweight = @sprintf(\"%23.16e\", weights[i])\n",
    "        println(\"w[\", i, \"] : $sweight\")\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "55d17247",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The triadiagonal matrix J :\n",
      "[0.0, 0.5773502691896257, 0.0, 0.0, 0.0]\n",
      "[0.5773502691896257, 0.0, 0.5163977794943222, 0.0, 0.0]\n",
      "[0.0, 0.5163977794943222, 0.0, 0.50709255283711, 0.0]\n",
      "[0.0, 0.0, 0.50709255283711, 0.0, 0.5039526306789697]\n",
      "[0.0, 0.0, 0.0, 0.5039526306789697, 0.0]\n",
      "The points :\n",
      "x[1] : -9.0617984593866407e-01\n",
      "x[2] : -5.3846931010568311e-01\n",
      "x[3] :  2.6689322308431607e-17\n",
      "x[4] :  5.3846931010568344e-01\n",
      "x[5] :  9.0617984593866419e-01\n",
      "The weights :\n",
      "w[1] :  2.3692688505618881e-01\n",
      "w[2] :  4.7862867049936775e-01\n",
      "w[3] :  5.6888888888888889e-01\n",
      "w[4] :  4.7862867049936575e-01\n",
      "w[5] :  2.3692688505618845e-01\n"
     ]
    }
   ],
   "source": [
    "makeGaussLegendre(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f45c11a6",
   "metadata": {},
   "source": [
    "We recognize the points and the weights of the\n",
    "Gauss-Legendre quadrature rule with 5 points\n",
    "that we computed above."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f795be8e",
   "metadata": {},
   "source": [
    "## 4. Gauss-Chebyshev quadrature"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d0ca68a",
   "metadata": {},
   "source": [
    "In Gauss-Chebyshev quadrature, the interval $[a,b]$ is also $[-1,+1]$,\n",
    "but the weight function $w(x)$ is $1/\\sqrt{1-x^2}$.\n",
    "This weight is useful in the calculation of improper integrals,\n",
    "for functions with vertical asymptotes at the end points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "290e2772",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "alphabetaChebyshev"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    alphabetaChebyshev(degree) computes the alpha\n",
    "and the beta coefficients for the tridiagonal matrix\n",
    "of the Cheybyshev polynomials.\n",
    "\"\"\"\n",
    "function alphabetaChebyshev(degree)\n",
    "    A = [2.0 for n=1:degree]\n",
    "    A[1] = 1.0;\n",
    "    B = [0.0 for n=1:degree]\n",
    "    alphas = [-B[k]/A[k] for k=1:degree]\n",
    "    C = [1.0 for n=1:degree]\n",
    "    betas = [sqrt(C[k+1]/(A[k]*A[k+1])) for k=1:degree-1]\n",
    "    return alphas, betas\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "90fe3219",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "makeGaussChebyshev"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    makeGaussChebyshev(degree::Int)\n",
    "\n",
    "computes the points and weights of the\n",
    "Gauss-Chebyshev quadrature with degree\n",
    "number of points via the eigenvalues and \n",
    "eigenvectors of a tridiagonal Jacobi matrix.\n",
    "\"\"\"\n",
    "function makeGaussChebyshev(degree::Int)\n",
    "    a, b = alphabetaChebyshev(degree)\n",
    "    J = diagm(a) + diagm(-1 => b) + diagm(+1 => b)\n",
    "    println(\"The triadiagonal matrix J :\")\n",
    "    for k=1:degree\n",
    "        println(J[k,:])\n",
    "    end\n",
    "    pts = eigvals(J)\n",
    "    println(\"The points :\")\n",
    "    for i=1:degree\n",
    "        spts = @sprintf(\"%23.16e\", pts[i])\n",
    "        root = -cos(((2.0*i-1)*pi)/(2.0*degree))\n",
    "        sroot = @sprintf(\"%23.16e\", root)\n",
    "        err = abs(pts[i] - root)\n",
    "        serr = @sprintf(\"%.2e\", err)\n",
    "        println(\"x[\", i, \"] : $spts : $sroot : $serr\")\n",
    "    end\n",
    "    Q = eigvecs(J)\n",
    "    # integrate(1/sqrt(1-x**2),(x,-1,1)) equals pi\n",
    "    println(\"The first row of the matrix of eigenvectors :\")\n",
    "    for i in 1:degree\n",
    "        println(@sprintf(\"%23.16e\", Q[1,i]))\n",
    "    end\n",
    "    weights = pi*(Q[1,:].^2)\n",
    "    println(\"The weights :\")\n",
    "    for i=1:degree\n",
    "        sweight = @sprintf(\"%23.16e\", weights[i])\n",
    "        println(\"w[\", i, \"] : $sweight\")\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "91f2bda2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The triadiagonal matrix J :\n",
      "[0.0, 0.7071067811865476, 0.0, 0.0, 0.0]\n",
      "[0.7071067811865476, 0.0, 0.5, 0.0, 0.0]\n",
      "[0.0, 0.5, 0.0, 0.5, 0.0]\n",
      "[0.0, 0.0, 0.5, 0.0, 0.5]\n",
      "[0.0, 0.0, 0.0, 0.5, 0.0]\n",
      "The points :\n",
      "x[1] : -9.5105651629515353e-01 : -9.5105651629515353e-01 : 0.00e+00\n",
      "x[2] : -5.8778525229247314e-01 : -5.8778525229247314e-01 : 0.00e+00\n",
      "x[3] : -5.8019460467133651e-17 : -6.1232339957367660e-17 : 3.21e-18\n",
      "x[4] :  5.8778525229247347e-01 :  5.8778525229247303e-01 : 4.44e-16\n",
      "x[5] :  9.5105651629515320e-01 :  9.5105651629515353e-01 : 3.33e-16\n",
      "The first row of the matrix of eigenvectors :\n",
      "-4.4721359549995793e-01\n",
      " 4.4721359549995809e-01\n",
      " 4.4721359549995804e-01\n",
      "-4.4721359549995793e-01\n",
      " 4.4721359549995754e-01\n",
      "The weights :\n",
      "w[1] :  6.2831853071795862e-01\n",
      "w[2] :  6.2831853071795907e-01\n",
      "w[3] :  6.2831853071795896e-01\n",
      "w[4] :  6.2831853071795862e-01\n",
      "w[5] :  6.2831853071795751e-01\n"
     ]
    }
   ],
   "source": [
    "makeGaussChebyshev(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d9138a0",
   "metadata": {},
   "source": [
    "We see that the eigenvalues agree with the explicit formulas for the roots of the Chebyshev polynomial."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "919f3c93",
   "metadata": {},
   "source": [
    "Let us apply the rule to $\\displaystyle \\int_{-1}^1 \\frac{1}{\\sqrt{1-x^2}} dx$."
   ]
  }
 ],
 "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"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}