Add documentation on Gaussian quadrature

doc/notebooks/Gaussian_Quadrature.ipynb

@@ -0,0 +1,170 @@
+{
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "#### Gaussian Quadrature\n",
+        "\n",
+        "The general idea of Gaussian quadrature is to find, in an interval, a set of points and weights such that the integral of a function in the interval will be approximate to the sum of the weighted function values at the points. That is:\n",
+        "<a name=\"eqn:num_int\"></a>\n",
+        "\\begin{equation*}\n",
+        "        \\int_{a}^{b}f(x) \\simeq \\sum_{i=1}^{n}w_{i}f(x_{i})\n",
+        "        \\tag{1}\n",
+        "\\end{equation*} \n",
+        "        \n",
+        "\n",
+        "This method is designed to give highly accurate results for $n$ points given that the function to integrate can be approximated accurately as a multiple of a polynomial of grade $m \\leq 2n-1$ by a weight function $w(x)$.\n",
+        "\n",
+        "<a name=\"eqn:poly_approx\"></a>\n",
+        "\\begin{equation}\n",
+        "        f(x)\\simeq L^{m}w(x)\n",
+        "        \\tag{2}\n",
+        "\\end{equation}\n",
+        "\n",
+        "The way to select the points and weights can be understood as follows:\n",
+        "\n",
+        "Given a set of $n$ polynomials $\\{p^{i}(x)\\}$ orthogonal in the interval $[a,b]$ with respect to a weight function $w(x)$:\n",
+        "<a name=\"eqn:poly_ort\"></a>\n",
+        "\\begin{equation}\n",
+        "        \\int_{a}^{b} p^{i}(x) p^{j}(x) w(x) = C_{i} \\delta_{ij}\n",
+        "        \\tag{3}\n",
+        "\\end{equation}\n",
+        "\n",
+        "Where $\\delta_{ij}$ is the Kronecker delta:\n",
+        "\n",
+        "\\begin{equation}\n",
+        "        \\delta_{ij} = \\begin{cases}\n",
+        "                0 & \\text{if } i \\neq j, \\\\\n",
+        "                1 & \\text{if } i = j.\n",
+        "        \\end{cases}\n",
+        "        \\tag{4}\n",
+        "\\end{equation}\n",
+        "\n",
+        "Any polynomial of order $m\\leq2n-1$ $(L^{m}(x))$ divided by $p^{n}(x)$ will give a quotient $Q(x)$ and remainder $R(x)$ of order at most $n-1$. Then $L^{m}(x)$ can be expressed as:\n",
+        "\\begin{equation}\n",
+        "        L^{m}(x) = Q(x)p^{n}(x) + R(x) \n",
+        "        \\tag{5}\n",
+        "\\end{equation}\n",
+        "\n",
+        "Also, as the remainder is of order lower or equal to $n-1$, it can be expressed as a linear combination of the orthogonal polynomials $\\{p^{j}(x)\\}$ for $j=0\\ldots n-1$\n",
+        "<a name=\"eqn:remaind_expanx\"></a>\n",
+        "\\begin{equation}\n",
+        "        R(x) = \\sum_{j=0}^{n-1}c_{j}p^{j}(x)\n",
+        "        \\tag{6}\n",
+        "\\end{equation}\n",
+        "\n",
+        "Substituting [5](#eqn:poly_div) in [2](#eqn:poly_approx) gives:\n",
+        "<a name=\"eqn:fcn_aprox_poly_div\"></a>\n",
+        "\\begin{equation}\n",
+        "        f(x)\\simeq [Q(x)p^{n}(x) + R(x)]w(x)\n",
+        "        \\tag{7}\n",
+        "\\end{equation}\n",
+        "\n",
+        "Substituting [6](#eqn:remaind_expanx):\n",
+        "<a name=\"eqn:fcn_aprox_poly_div_exp\"></a>\n",
+        "\\begin{equation}\n",
+        "        f(x)\\simeq [Q(x)p^{n}(x) + \\sum_{j=0}^{n-1}c_{j}p^{j}(x)]w(x)\n",
+        "        \\tag{8}\n",
+        "\\end{equation}\n",
+        "\n",
+        "\n",
+        "Then the approximate function integral can be expressed as:\n",
+        "\n",
+        "<a name=\"eqn:aprox_int\"></a>\n",
+        "<a name=\"eqn:aprox_int_exp\"></a>\n",
+        "<a name=\"eqn:aprox_int_exp_reord\"></a>\n",
+        "\\begin{align}\n",
+        "        \\int_{a}^{b}f(x) &\\simeq \\int_{a}^{b}[Q(x)p^{n}(x) + \\sum_{j=0}^{n-1}c_{j}p^{j}(x)]w(x)dx \\tag{9} \\\\\n",
+        "        &\\simeq \\int_{a}^{b}Q(x)p^{n}(x)w(x) dx + \\int_{a}^{b}\\sum_{j=0}^{n-1}c_{j}p^{j}(x)w(x)dx \\tag{10} \\\\ \n",
+        "        &\\simeq \\int_{a}^{b}Q(x)p^{n}(x)w(x) dx + \\sum_{j=0}^{n-1}c_{j}\\int_{a}^{b}p^{j}(x)w(x)dx \\tag{11} \t\t\n",
+        "\\end{align} \n",
+        "\n",
+        "where the first integral on the right side is 0 because $p^{n}(x)$ is orthogonal to any polynomial of lower order in the $[a,b]$ interval with respect to the weight function $w(x)$:\n",
+        "<a name=\"eqn:orth_quot\"></a>\n",
+        "\\begin{equation}\n",
+        "        \\int_{a}^{b}Q(x)p^{n}(x)w(x) dx = 0\n",
+        "        \\tag{12}\n",
+        "\\end{equation}\n",
+        "\n",
+        "Also substituting [8](#eqn:fcn_aprox_poly_div_exp) in [1](#eqn:num_int):\n",
+        "\n",
+        "\\begin{align}\n",
+        "        \\int_{a}^{b}f(x) \\simeq \\sum_{i=1}^{n}w_{i}f(x_{i}) &= \\sum_{i=1}^{n}w_{i} \\left[Q(x_{i})p^{n}(x_{i}) + \\sum_{j=0}^{n-1}c_{j}p^{j}(x_{i})\\right]w(x_{i}) \\tag{13}\\\\\n",
+        "        &= \\sum_{i=1}^{n}w_{i} Q(x_{i})p^{n}(x_{i}) + \\sum_{i=1}^{n}w_{i}\\sum_{j=0}^{n-1}c_{j}p^{j}(x_{i})w(x_{i}) \\tag{14}\\\\\n",
+        "        &= \\sum_{i=1}^{n}w_{i} Q(x_{i})p^{n}(x_{i}) + \\sum_{j=0}^{n-1}c_{j}\\sum_{i=1}^{n}w_{i}p^{j}(x_{i})w(x_{i}) \\tag{15}\n",
+        "\\end{align} \n",
+        "\n",
+        "To make numeric the integral of the function as close to the exact value as possible the $\\{x_{i},w_{i}\\}$ values are selected such that the sums in [15](#eq:num_int_dec) are equal to the integrals in [11](#eqn:aprox_int_exp_reord). The $\\{x_{i}\\}$ vales are then selected to be the roots of $p^{n}(x)$. In this way the leftmost sum is zero (the exact value of the integral [12](#eqn:orth_quot)) for any possible set of weights $\\{w_{i}\\}$.\n",
+        "\\begin{equation}\n",
+        "        \\sum_{i=1}^{n}w_{i} Q(x_{i})p^{n}(x_{i}) = 0 \\tag{16}\n",
+        "\\end{equation}\n",
+        "\n",
+        "The corresponding weights are then selected such that for all $n-1$ polynomials $p^{j}(x)$, the sum give their exact integral value:\n",
+        "\\begin{equation}\n",
+        "        \\sum_{i=1}^{n}w_{i}p^{j}(x_{i})w(x_{i}) = \\int_{a}^{b}p^{j}(x)w(x)dx \\tag{17}\n",
+        "\\end{equation}\n",
+        "\n",
+        "This condition guarantees that the rightmost sum in [15](#eq:num_int_dec) gives the exact value of the rightmost integral in [11](#eqn:aprox_int_exp_reord). \n",
+        "\n",
+        "The weights that fulfill [17](#eqn:int_sum_weights) can be then found by solving the following linear system.\n",
+        "\n",
+        "\\begin{equation}\n",
+        "        \\begin{bmatrix}\n",
+        "                p^{0}(x_{1})w(x_{1})&p^{0}(x_{2})w(x_{2})&p^{0}(x_{3})&\\dots &p^{0}(x_{n})w(x_{n})\\\\\n",
+        "                p^{1}(x_{1})w(x_{1})&p^{1}(x_{2})w(x_{2})&p^{1}(x_{3})&\\dots &p^{1}(x_{n})w(x_{n})\\\\\n",
+        "                p^{2}(x_{1})w(x_{1})&p^{2}(x_{2})w(x_{2})&p^{2}(x_{3})&\\dots &p^{2}(x_{n})w(x_{n})\\\\\n",
+        "                \\vdots &\\vdots &\\vdots &\\ddots &\\vdots \\\\\n",
+        "                p^{(n-1)}(x_{1})w(x_{1})&p^{(n-1)}(x_{2})w(x_{2})&p^{(n-1)}(x_{3})&\\dots &p^{(n-1)}(x_{n})w(x_{n})\\\\\n",
+        "        \\end{bmatrix}\n",
+        "        \\cdot \n",
+        "        \\begin{bmatrix}\n",
+        "                w_{0}\\\\w_{1}\\\\w_{2}\\\\\\vdots \\\\w_{n}\n",
+        "        \\end{bmatrix}\n",
+        "        =\n",
+        "        \\begin{bmatrix}\n",
+        "                \\int_{a}^{b} p^{0}(x)w(x)\\\\\n",
+        "                \\int_{a}^{b} p^{1}(x)w(x)\\\\\n",
+        "                \\int_{a}^{b} p^{2}(x)w(x)\\\\\n",
+        "                \\vdots \\\\\n",
+        "                \\int_{a}^{b} p^{(n-1)}(x)w(x)\\\\\n",
+        "        \\end{bmatrix}\n",
+        "        \\tag{18}\n",
+        "\\end{equation}\n",
+        "\n",
+        "Given a family of orthogonal polynomials with respect to a weight function in a determined interval, this procedure allows obtaining a set of $n$ optimal points $\\{x_{i}\\}$ and their corresponding weights $\\{w_{i}\\}$ for integrating a function in this interval. The function needs to be approximateable as a polynomial of grade at most $2n-1$. Also, if the function can be represented exactly as a polynomial of grade at most $2n-1$, the numerical integration is exact.\n",
+        "\n",
+        "There are several Gauss quadrature rules for different choices of orthogonal polynomial sets, intervals and weight functions. Grid offers several one-dimensional quadratures for one-dimensional integration such as:\n",
+        "\n",
+        "\t\n",
+        "-   Gauss-Chebyshev\n",
+        "\n",
+        "-   Gauss-Legendre\n",
+        "\n",
+        "-   Gauss-Laguerre"
+      ]
+    }
+  ],
+  "metadata": {
+    "kernelspec": {
+      "display_name": "grid_env",
+      "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.10.12"
+    },
+    "orig_nbformat": 4
+  },
+  "nbformat": 4,
+  "nbformat_minor": 2
+}