{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "DEVOIR MAISON 3\n",
    "\n",
    "---\n",
    "# Formule de quadrature de Fejér\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.special import roots_legendre\n",
    "from scipy.integrate import quad\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le but de ce DM est de programmer la méthode de quadrature de Fejér pour calculer l'intégrale \n",
    "\\begin{equation*}\n",
    "I(f) = \\int\\limits_{-1}^{1} f( x ) dx,\n",
    "\\end{equation*}\n",
    "où $ f $ est une fonction continue sur $ [ -1, 1 ] $, et de la comparer avec la méthode de Gauss et la formule composée de Simpson."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous allons tester les différentes méthodes de quadrature sur les fonctions suivantes:\n",
    "\\begin{align*}\n",
    "& f_0( x ) = e^x, \\quad f_1( x ) = \\frac{ 1 }{ 1 + 16 x^2 }, \\\\\n",
    "& f_2( x ) = | x^2 - 0.25 |^3, \\quad f_3( x ) = | x + 0.5 |^{1/2}. \n",
    "\\end{align*}\n",
    "\n",
    "On remarque que $ f_0 $ est une fonction entière, $ f_1 $ est analytique au voisinage de $ [ -1, 1 ] $, $ f_2 \\in C^3 $ et $ f_3 \\in C $.\n",
    "\n",
    "**Question 1.**\n",
    "> Définir les fonctions $ f_i $, $ i = 0,  1, 2, 3 $.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def f0(x): \n",
    "    return np.exp(x)\n",
    "\n",
    "def f1(x):\n",
    "    return 1 / (1 + 16*np.power(x, 2))\n",
    "\n",
    "def f2(x):\n",
    "    return np.power(np.abs(x**2 - 1/4), 3)\n",
    "\n",
    "def f3(x):\n",
    "    return np.power(np.abs(x+1/2), 1/2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Formule de quadrature de Gauss à $ N $ points\n",
    "\n",
    "On rappelle que la formule de quadrature de Gauss à $ N $ points est la formule qui s'écrit sous la forme \n",
    "\\begin{equation*}\n",
    "J_G( f ) = \\sum\\limits_{ k = 1 }^N \\lambda_k f( x_k ),\n",
    "\\tag{1}\n",
    "\\end{equation*}\n",
    "où les points $ x_k $, $ k = 1, \\ldots, N $, sont les racines du polynome de Legendre $ P_N $ défini par recurrence via la formule:\n",
    "\\begin{align*}\n",
    "& P_0( x ) = 1, \\\\\n",
    "& P_1( x ) = x, \\\\\n",
    "& P_{ N }( x ) = \\frac{ ( 2 N - 1 ) x P_{ N - 1 }( x ) - ( N - 1 ) P_{ N - 2 }( x ) }{ N }\n",
    "\\end{align*}\n",
    "et $ \\lambda_k $ sont donnés par la formule\n",
    "\\begin{equation*}\n",
    "\\lambda_k = \\frac{ 1 - x_k^2 }{ N P^2_{ N-1 }( x_k ) }.\n",
    "\\end{equation*}\n",
    "\n",
    "\n",
    "**Question 2.** \n",
    "> 1. Programmez une fonction `gauss` qui prend en argument une fonction `f` et un entier `N` et qui retourne une valeur approchée de $ I(f) $ obtenue par la formule de quadrature de Gauss (1). On pourra utiliser la fonction `roots_legendre` du module `scipy.special`.\n",
    "> 2. Tester votre fonction en comparant la valeur de $ \\int\\limits_{-1}^1 f_i(x)dx $ calculée à l'aide de la fonction `quad` de `scipy.integrate` et à l'aide de votre fonction `gauss` pour les fonction $ f_i $, $ i = 0, 1, 2, 3 $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gauss(f, N):\n",
    "    roots, weights = roots_legendre(N)\n",
    "    return np.sum(weights * f(roots))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calcule de I(f) par la méthode de gauss et par la formule quadratique pour la fonction f0\n",
      "Pour n = 1, gauss = 2.0 et quad = 2.3504023872876028\n",
      "Pour n = 2, gauss = 2.3426960879097307 et quad = 2.3504023872876028\n",
      "Pour n = 3, gauss = 2.3503369286800115 et quad = 2.3504023872876028\n",
      "Pour n = 4, gauss = 2.3504020921563766 et quad = 2.3504023872876028\n",
      "Pour n = 5, gauss = 2.350402386462826 et quad = 2.3504023872876028\n",
      "Pour n = 6, gauss = 2.3504023872860342 et quad = 2.3504023872876028\n",
      "Pour n = 7, gauss = 2.3504023872876014 et quad = 2.3504023872876028\n",
      "Pour n = 8, gauss = 2.3504023872876028 et quad = 2.3504023872876028\n",
      "Pour n = 9, gauss = 2.3504023872876023 et quad = 2.3504023872876028\n",
      "Pour n = 10, gauss = 2.350402387287602 et quad = 2.3504023872876028\n",
      "\n",
      "Calcule de I(f) par la méthode de gauss et par la formule quadratique pour la fonction f1\n",
      "Pour n = 1, gauss = 2.0 et quad = 0.6629088318340162\n",
      "Pour n = 2, gauss = 0.31578947368421056 et quad = 0.6629088318340162\n",
      "Pour n = 3, gauss = 0.9937106918238988 et quad = 0.6629088318340162\n",
      "Pour n = 4, gauss = 0.5118212522733177 et quad = 0.6629088318340162\n",
      "Pour n = 5, gauss = 0.7721547547946116 et quad = 0.6629088318340162\n",
      "Pour n = 6, gauss = 0.6029222321445539 et quad = 0.6629088318340162\n",
      "Pour n = 7, gauss = 0.7019205771618546 et quad = 0.6629088318340162\n",
      "Pour n = 8, gauss = 0.6400200424327264 et quad = 0.6629088318340162\n",
      "Pour n = 9, gauss = 0.6772078574245298 et quad = 0.6629088318340162\n",
      "Pour n = 10, gauss = 0.6543127409862085 et quad = 0.6629088318340162\n",
      "\n",
      "Calcule de I(f) par la méthode de gauss et par la formule quadratique pour la fonction f2\n",
      "Pour n = 1, gauss = 0.03125 et quad = 0.09375000000000001\n",
      "Pour n = 2, gauss = 0.0011574074074074067 et quad = 0.09375000000000001\n",
      "Pour n = 3, gauss = 0.06152777777777782 et quad = 0.09375000000000001\n",
      "Pour n = 4, gauss = 0.08579899983141438 et quad = 0.09375000000000001\n",
      "Pour n = 5, gauss = 0.09724206349206342 et quad = 0.09375000000000001\n",
      "Pour n = 6, gauss = 0.09293246353076502 et quad = 0.09375000000000001\n",
      "Pour n = 7, gauss = 0.09347309930500289 et quad = 0.09375000000000001\n",
      "Pour n = 8, gauss = 0.09415589118481105 et quad = 0.09375000000000001\n",
      "Pour n = 9, gauss = 0.0935821376879856 et quad = 0.09375000000000001\n",
      "Pour n = 10, gauss = 0.09370357485096373 et quad = 0.09375000000000001\n",
      "\n",
      "Calcule de I(f) par la méthode de gauss et par la formule quadratique pour la fonction f3\n",
      "Pour n = 1, gauss = 1.4142135623730951 et quad = 1.4604471317871044\n",
      "Pour n = 2, gauss = 1.3160740129524924 et quad = 1.4604471317871044\n",
      "Pour n = 3, gauss = 1.5468727439784977 et quad = 1.4604471317871044\n",
      "Pour n = 4, gauss = 1.4734441274369583 et quad = 1.4604471317871044\n",
      "Pour n = 5, gauss = 1.4158420574869128 et quad = 1.4604471317871044\n",
      "Pour n = 6, gauss = 1.4926838373564866 et quad = 1.4604471317871044\n",
      "Pour n = 7, gauss = 1.467885304909335 et quad = 1.4604471317871044\n",
      "Pour n = 8, gauss = 1.4370627889357577 et quad = 1.4604471317871044\n",
      "Pour n = 9, gauss = 1.4783803880783275 et quad = 1.4604471317871044\n",
      "Pour n = 10, gauss = 1.4652835957768406 et quad = 1.4604471317871044\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for f in [f0, f1, f2, f3]:\n",
    "    print(f'Calcule de I(f) par la méthode de gauss et par la formule quadratique pour la fonction {f.__name__}')\n",
    "    for n in range(1, 11):\n",
    "        print(f\"Pour n = {n}, gauss = {gauss(f, n)} et quad = {quad(f, -1, 1)[0]}\")\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Formule composée de Simpson à $ N $ points  \n",
    "\n",
    "On remarque que la formule de quadrature de Gauss à $ N $ points (1) utilise exactement $ N $ *points de quadrature*, i.e. $ N $ points où $ f $ est évaluée. \n",
    "\n",
    "La formule composée de Simpson utilise toujours un nombre *impair* de points de quadrature. Plus particulièrment, la formule composée de Simpson à $ N $ points de quadrature, pour $ N $ impair, s'écrit\n",
    "\\begin{equation*}\n",
    "J_S( f ) = \\frac{ h }{ 3 } \\sum\\limits_{ k = 1 }^{ m }[ f( x_{ 2k - 2 } ) + 4 f( x_{2k - 1} ) + f( x_{ 2k } ) ]\n",
    "\\tag{2}\n",
    "\\end{equation*}\n",
    "où $ m = N // 2 $, $ h = \\frac{ b - a }{ 2m } $, $ x_k = a + k h $.\n",
    "\n",
    "**Question 3.** \n",
    "> 1. Programmez une fonction `simpson` qui prend en argument une fonction `f` et un entier impair `N` et qui retourne une valeur approchée de $ I(f) $ obtenue par la formule composée de Simpson à $ N $ points de quadrature (2). *On privilégiera une implémentation qui n'utilise pas de boucle for.* \n",
    "> 2. Tester votre fonction en comparant la valeur de $ \\int\\limits_{-1}^1 f_i(x)dx $ calculée à l'aide de la fonction `quad` de `scipy.integrate` et à l'aide de votre fonction `simpson` pour les fonction $ f_i $, $ i = 0, 1, 2, 3 $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simpson(f, N):\n",
    "    if N % 2 == 0:\n",
    "        raise ValueError(\"N doit est impair.\")\n",
    "    \n",
    "    h = 2 / (2 * (N - 1) // 2)\n",
    "    fx = f(np.linspace(-1, 1, N))\n",
    "    \n",
    "    return (h / 3) * (fx[0] + 4 * fx[1:-1:2].sum() + 2 * fx[2:-1:2].sum() + fx[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calcule de I(f) par la méthode de simpson et par la formule quadratique pour la fonction f0\n",
      "Pour n = 3, simpson = 2.362053756543496 et quad = 2.3504023872876028\n",
      "Pour n = 5, simpson = 2.3511948318802554 et quad = 2.3504023872876028\n",
      "Pour n = 7, simpson = 2.350561486811035 et quad = 2.3504023872876028\n",
      "Pour n = 9, simpson = 2.3504530172422795 et quad = 2.3504023872876028\n",
      "Pour n = 11, simpson = 2.3504231806814837 et quad = 2.3504023872876028\n",
      "Pour n = 13, simpson = 2.350412429522249 et quad = 2.3504023872876028\n",
      "Pour n = 15, simpson = 2.3504078125830676 et quad = 2.3504023872876028\n",
      "\n",
      "Calcule de I(f) par la méthode de simpson et par la formule quadratique pour la fonction f1\n",
      "Pour n = 3, simpson = 1.372549019607843 et quad = 0.6629088318340162\n",
      "Pour n = 5, simpson = 0.6196078431372549 et quad = 0.6629088318340162\n",
      "Pour n = 7, simpson = 0.7271053809651714 et quad = 0.6629088318340162\n",
      "Pour n = 9, simpson = 0.6431372549019607 et quad = 0.6629088318340162\n",
      "Pour n = 11, simpson = 0.6738214805917826 et quad = 0.6629088318340162\n",
      "Pour n = 13, simpson = 0.6583227633851206 et quad = 0.6629088318340162\n",
      "Pour n = 15, simpson = 0.6650749418306187 et quad = 0.6629088318340162\n",
      "\n",
      "Calcule de I(f) par la méthode de simpson et par la formule quadratique pour la fonction f2\n",
      "Pour n = 3, simpson = 0.3020833333333333 et quad = 0.09375000000000001\n",
      "Pour n = 5, simpson = 0.14583333333333331 et quad = 0.09375000000000001\n",
      "Pour n = 7, simpson = 0.10842001981405272 et quad = 0.09375000000000001\n",
      "Pour n = 9, simpson = 0.09765625 et quad = 0.09375000000000001\n",
      "Pour n = 11, simpson = 0.09526680000000004 et quad = 0.09375000000000001\n",
      "Pour n = 13, simpson = 0.09473593964334705 et quad = 0.09375000000000001\n",
      "Pour n = 15, simpson = 0.09423859015254817 et quad = 0.09375000000000001\n",
      "\n",
      "Calcule de I(f) par la méthode de simpson et par la formule quadratique pour la fonction f3\n",
      "Pour n = 3, simpson = 1.5867595924414424 et quad = 1.4604471317871044\n",
      "Pour n = 5, simpson = 1.2243442024918718 et quad = 1.4604471317871044\n",
      "Pour n = 7, simpson = 1.4840004641551812 et quad = 1.4604471317871044\n",
      "Pour n = 9, simpson = 1.4401918987573803 et quad = 1.4604471317871044\n",
      "Pour n = 11, simpson = 1.4713568906975611 et quad = 1.4604471317871044\n",
      "Pour n = 13, simpson = 1.4149243202984132 et quad = 1.4604471317871044\n",
      "Pour n = 15, simpson = 1.4670273366115245 et quad = 1.4604471317871044\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for f in [f0, f1, f2, f3]:\n",
    "    print(f'Calcule de I(f) par la méthode de simpson et par la formule quadratique pour la fonction {f.__name__}')\n",
    "    for n in range(3, 16, 2):\n",
    "        print(f\"Pour n = {n}, simpson = {simpson(f, n)} et quad = {quad(f, -1, 1)[0]}\")\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Formule de quadrature de Fejér à $N$ points\n",
    "La formule de quadrature de Fejér fait intervenir les points et les polynomes de Tchebychev. On rappelle que les polynomes de Tchebychev $ T_N $ sont des polynomes définis sur $ [ -1, 1 ] $ via la formule \n",
    "$$\n",
    "T_N( \\cos \\theta ) = \\cos( N \\theta ).\n",
    "$$\n",
    "On peut également les définir via la formule recurrente suivante:\n",
    "\\begin{align*}\n",
    "& T_0( x ) = 1, \\\\\n",
    "& T_1( x ) = x, \\\\\n",
    "& T_N( x ) = 2 x T_{ N-1 }( x ) - T_{ N-2 }( x ).\n",
    "\\end{align*}\n",
    "\n",
    "Les points de Tchebychev $ x_{ k, N } $, $ k = 1, \\ldots, N $, sur $ [ -1, 1 ] $ sont définis par la formule\n",
    "\\begin{equation*}\n",
    "x_{ k, N } = \\cos\\left( \\frac{ 2k - 1 }{ 2 N } \\pi \\right), \\quad k = 1, 2, \\ldots, N.\n",
    "\\end{equation*}\n",
    "\n",
    "**Question 4.** \n",
    "> Rappeler pourquoi les points de Tchebychev $ x_{k, N} $, $ k = 1, \\ldots, N $, sont racines du polynome de Tchebychev $ T_N $.\n",
    "\n",
    "**Question 5.**\n",
    "> 1. Définir une fonction `poly_tchebychev` qui prend en argument un `ndarray` `x` et un entier positif `N` et qui retourne l'évaluation de $ T_N $ en $ x $. \n",
    "> 2. Définir une fonction `points_tchebychev` qui prend en argument un entier positif $ N $ et qui retourne un `ndarray` des points de Tchebychev $ x_{ k, N } $, $ k = 1, \\ldots, N $. *On privilégiera une implémentation qui n'utilise pas de boucle for.*\n",
    "> 3. Tester vos fonctions en évaluant un polynome de Tchebychev $ T_N $ en les points de Tchebychev $ x_{ k, N } $, $ k = 1, \\ldots, N $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Montrons que les points de Tchebychev T_N( x_{k, N}  ) sont racines du polynomes de Tchebychev $T_N$\n",
    "\n",
    "On sait que $\\cos(A-B) = \\cos(A)\\cos(B) - \\sin(A)\\sin(B)$\n",
    "\n",
    "Pour commencer, on a $$\n",
    "T_N( \\cos \\theta ) = \\cos( N \\theta ).\n",
    "$$\n",
    "Donc $T_N( x_{k, N}  ) = T_N( \\cos\\left( \\frac{ 2k - 1 }{ 2 N } \\pi \\right) ) = \\cos\\left( \\frac{ N(2k - 1) }{ 2 N } \\pi \\right) = \\cos\\left( \\frac{ 2k - 1 }{ 2 } \\pi \\right) = \\cos\\left( k\\pi - \\frac{\\pi}{2} \\right) = \\cos(k\\pi)cos(\\frac{\\pi}{2}) - \\sin(k\\pi)\\sin(\\frac{\\pi}{2})$ pour $k=1,{\\dots},N$\n",
    "\n",
    "Or $cos(\\frac{\\pi}{2}) = 0$ et $\\sin(k\\pi)$ pour $k=1,{\\dots},N$\n",
    "\n",
    "Donc $T_N( x_{k, N}  ) = 0$ donc  $x_{k, N}$ sont racines du polynomes de Tchebychev $T_N$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "def poly_tchebychev(x, N):\n",
    "    if N == 0:\n",
    "        return np.ones_like(x)\n",
    "    elif N == 1:\n",
    "        return x\n",
    "    else:\n",
    "        return 2 * x * poly_tchebychev(x, N-1) - poly_tchebychev(x, N-2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def points_tchebychev(N):\n",
    "    k = np.arange(1, N+1)\n",
    "    return np.cos((2 * k - 1) * np.pi / (2 * N))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pour N = 2\n",
      "Les points de Tchebychev sont [ 0.70710678 -0.70710678]\n",
      "L'evaluation du polynome de Tchebychev Tn en ces points est [ 2.22044605e-16 -2.22044605e-16]\n",
      "\n",
      "Pour N = 3\n",
      "Les points de Tchebychev sont [ 8.66025404e-01  6.12323400e-17 -8.66025404e-01]\n",
      "L'evaluation du polynome de Tchebychev Tn en ces points est [ 3.33066907e-16 -1.83697020e-16 -3.33066907e-16]\n",
      "\n",
      "Pour N = 4\n",
      "Les points de Tchebychev sont [ 0.92387953  0.38268343 -0.38268343 -0.92387953]\n",
      "L'evaluation du polynome de Tchebychev Tn en ces points est [-2.22044605e-16 -2.22044605e-16  1.11022302e-16 -2.22044605e-16]\n",
      "\n",
      "Pour N = 5\n",
      "Les points de Tchebychev sont [ 9.51056516e-01  5.87785252e-01  6.12323400e-17 -5.87785252e-01\n",
      " -9.51056516e-01]\n",
      "L'evaluation du polynome de Tchebychev Tn en ces points est [-4.44089210e-16  0.00000000e+00  3.06161700e-16 -7.77156117e-16\n",
      "  4.44089210e-16]\n",
      "\n",
      "Pour N = 6\n",
      "Les points de Tchebychev sont [ 0.96592583  0.70710678  0.25881905 -0.25881905 -0.70710678 -0.96592583]\n",
      "L'evaluation du polynome de Tchebychev Tn en ces points est [ 7.77156117e-16 -3.33066907e-16 -2.22044605e-16 -7.77156117e-16\n",
      "  6.66133815e-16 -1.55431223e-15]\n",
      "\n",
      "Pour N = 7\n",
      "Les points de Tchebychev sont [ 9.74927912e-01  7.81831482e-01  4.33883739e-01  6.12323400e-17\n",
      " -4.33883739e-01 -7.81831482e-01 -9.74927912e-01]\n",
      "L'evaluation du polynome de Tchebychev Tn en ces points est [ 2.22044605e-16  0.00000000e+00  4.44089210e-16 -4.28626380e-16\n",
      "  4.44089210e-16 -3.55271368e-15 -4.88498131e-15]\n",
      "\n",
      "Pour N = 8\n",
      "Les points de Tchebychev sont [ 0.98078528  0.83146961  0.55557023  0.19509032 -0.19509032 -0.55557023\n",
      " -0.83146961 -0.98078528]\n",
      "L'evaluation du polynome de Tchebychev Tn en ces points est [-8.88178420e-16  1.11022302e-16  5.55111512e-16 -5.55111512e-16\n",
      "  6.10622664e-16 -2.44249065e-15 -1.44328993e-15 -8.88178420e-16]\n",
      "\n",
      "Pour N = 9\n",
      "Les points de Tchebychev sont [ 9.84807753e-01  8.66025404e-01  6.42787610e-01  3.42020143e-01\n",
      "  6.12323400e-17 -3.42020143e-01 -6.42787610e-01 -8.66025404e-01\n",
      " -9.84807753e-01]\n",
      "L'evaluation du polynome de Tchebychev Tn en ces points est [-1.99840144e-15 -8.88178420e-16  4.44089210e-16 -8.88178420e-16\n",
      "  5.51091060e-16 -2.44249065e-15 -4.44089210e-16 -3.21964677e-15\n",
      "  1.99840144e-15]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for n in range(2, 10):\n",
    "    xk = points_tchebychev(n)\n",
    "    Tn = poly_tchebychev(xk, n)\n",
    "    print(f\"Pour N = {n}\")\n",
    "    print(f\"Les points de Tchebychev sont {xk}\")\n",
    "    print(f\"L'evaluation du polynome de Tchebychev Tn en ces points est {Tn}\")\n",
    "    print(\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La formule de quadrature de Fejér à $ N $ points s'écrit sous la forme \n",
    "\\begin{equation*}\n",
    "J_F( f ) = \\sum\\limits_{ k = 1 }^N \\lambda_k f( x_k ),\n",
    "\\tag{1}\n",
    "\\end{equation*}\n",
    "où $ x_k = x_{ k, N } $, $ k = 1, \\ldots, N $, sont les points de Tchebychev\n",
    "et $ \\lambda_k $ sont donnés par la formule\n",
    "\\begin{equation*}\n",
    "\\lambda_k = \\frac{ 2 }{ N }\\left( 1 - \\sum_{ m = 1 }^{ \\lfloor N/2 \\rfloor }  \\frac{ 2T_{ 2m }( x_k ) }{ 4 m^2 - 1 }  \\right).\n",
    "\\end{equation*}\n",
    "\n",
    "**Question 6.** \n",
    "> 1. Programmez une fonction `fejer_points_weights` qui prend en argument un entier `N` et qui retourne un ndarray `x` des points de Tchebychev $ x_{k, N} $ et un ndarray `lam` des poids $ \\lambda_k $ de la formule de Fejér.\n",
    "> 2. Programmez une fonction `fejer` qui prend en argument une fonction `f` et un entier `N`et qui retourne une valeur approchée de $ I $ obtenue par la formule de quadrature de Fejér. \n",
    "> 3. Tester votre fonction en comparant la valeur de $ \\int\\limits_{-1}^1 f_i(x)dx $ calculée à l'aide de la fonction `quad` de `scipy.integrate` et à l'aide de votre fonction `fejer` pour les fonction $ f_i $, $ i = 0, 1, 2, 3 $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fejer_points_weights(N):\n",
    "    xk = points_tchebychev(N)\n",
    "    lamk = np.zeros(N)\n",
    "    for k in range(N):\n",
    "        s = 0\n",
    "        for m in range(1, N//2+1):\n",
    "            T = poly_tchebychev(xk[k], 2 * m)\n",
    "            s += 2 * T / (4 * np.power(m, 2) - 1)\n",
    "        lamk[k] = 2 / N * (1 - s)\n",
    "    return xk, lamk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fejer(f, N):\n",
    "    xk, lamk = fejer_points_weights(N)\n",
    "    return np.sum(lamk * f(xk))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calcule de I(f) par la méthode de fejer et par la formule quadratique pour la fonction f0\n",
      "Pour n = 1, fejer = 2.0 et quad = 2.3504023872876028\n",
      "Pour n = 2, fejer = 2.521183673042712 et quad = 2.3504023872876028\n",
      "Pour n = 3, fejer = 2.35469453390679 et quad = 2.3504023872876028\n",
      "Pour n = 4, fejer = 2.351164449305366 et quad = 2.3504023872876028\n",
      "Pour n = 5, fejer = 2.3504110924399146 et quad = 2.3504023872876028\n",
      "Pour n = 6, fejer = 2.3504049906823328 et quad = 2.3504023872876028\n",
      "Pour n = 7, fejer = 2.350402405080618 et quad = 2.3504023872876028\n",
      "Pour n = 8, fejer = 2.350402393654548 et quad = 2.3504023872876028\n",
      "Pour n = 9, fejer = 2.3504023873162794 et quad = 2.3504023872876028\n",
      "Pour n = 10, fejer = 2.350402387298774 et quad = 2.3504023872876028\n",
      "\n",
      "Calcule de I(f) par la méthode de fejer et par la formule quadratique pour la fonction f1\n",
      "Pour n = 1, fejer = 2.0 et quad = 0.6629088318340162\n",
      "Pour n = 2, fejer = 0.22222222222222224 et quad = 0.6629088318340162\n",
      "Pour n = 3, fejer = 1.179487179487179 et quad = 0.6629088318340162\n",
      "Pour n = 4, fejer = 0.4761904761904763 et quad = 0.6629088318340162\n",
      "Pour n = 5, fejer = 0.7960396039603961 et quad = 0.6629088318340162\n",
      "Pour n = 6, fejer = 0.5849607182940518 et quad = 0.6629088318340162\n",
      "Pour n = 7, fejer = 0.7170436790978874 et quad = 0.6629088318340162\n",
      "Pour n = 8, fejer = 0.6339761526632889 et quad = 0.6629088318340162\n",
      "Pour n = 9, fejer = 0.6801443554906159 et quad = 0.6629088318340162\n",
      "Pour n = 10, fejer = 0.6516647293527041 et quad = 0.6629088318340162\n",
      "\n",
      "Calcule de I(f) par la méthode de fejer et par la formule quadratique pour la fonction f2\n",
      "Pour n = 1, fejer = 0.03125 et quad = 0.09375000000000001\n",
      "Pour n = 2, fejer = 0.03125 et quad = 0.09375000000000001\n",
      "Pour n = 3, fejer = 0.12847222222222227 et quad = 0.09375000000000001\n",
      "Pour n = 4, fejer = 0.11785113019775793 et quad = 0.09375000000000001\n",
      "Pour n = 5, fejer = 0.10458333333333332 et quad = 0.09375000000000001\n",
      "Pour n = 6, fejer = 0.09359684369075255 et quad = 0.09375000000000001\n",
      "Pour n = 7, fejer = 0.09403994748859118 et quad = 0.09375000000000001\n",
      "Pour n = 8, fejer = 0.09415848041941287 et quad = 0.09375000000000001\n",
      "Pour n = 9, fejer = 0.09340379593991482 et quad = 0.09375000000000001\n",
      "Pour n = 10, fejer = 0.0937828936137298 et quad = 0.09375000000000001\n",
      "\n",
      "Calcule de I(f) par la méthode de fejer et par la formule quadratique pour la fonction f3\n",
      "Pour n = 1, fejer = 1.4142135623730951 et quad = 1.4604471317871044\n",
      "Pour n = 2, fejer = 1.5537739740300371 et quad = 1.4604471317871044\n",
      "Pour n = 3, fejer = 1.5740169694012405 et quad = 1.4604471317871044\n",
      "Pour n = 4, fejer = 1.4306412648377405 et quad = 1.4604471317871044\n",
      "Pour n = 5, fejer = 1.4523326100552618 et quad = 1.4604471317871044\n",
      "Pour n = 6, fejer = 1.4976002258596133 et quad = 1.4604471317871044\n",
      "Pour n = 7, fejer = 1.4523308297718713 et quad = 1.4604471317871044\n",
      "Pour n = 8, fejer = 1.4578810195719196 et quad = 1.4604471317871044\n",
      "Pour n = 9, fejer = 1.4808716293756385 et quad = 1.4604471317871044\n",
      "Pour n = 10, fejer = 1.4553892699618816 et quad = 1.4604471317871044\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for f in [f0, f1, f2, f3]:\n",
    "    print(f'Calcule de I(f) par la méthode de fejer et par la formule quadratique pour la fonction {f.__name__}')\n",
    "    for n in range(1, 11):\n",
    "        print(f\"Pour n = {n}, fejer = {fejer(f, n)} et quad = {quad(f, -1, 1)[0]}\")\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparaison des méthodes\n",
    "\n",
    "**Question 7.**\n",
    "> 1. Pour chacune des fonction $ f_i $, $ i = 0, 1, 2, 3 $, tracez sur un graphique le logarithme d'erreur entre la valeur \"exacte\" de $ I $ (calculée via `quad`) et l'intégrale approchée calculée via les méthodes de Simpson composée, la méthode de Gauss et la méthode de Fejér à $ N $ points en fonction de $N$ (on prendra que des valeurs impairs de $N$). Ajoutez à chacun des quatre graphiques un titre et une légende. \n",
    "> 2. Commentez les résultats obtenus. Quelle méthode choisiriez-vous pour l'intégration de fonctions entières? Pour l'intégration de fonctions peu regulières?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 78752 characters]\n",
      "text/plain": [
       "<Figure size 1500x1000 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 20\n",
    "\n",
    "figure = plt.figure(figsize=(15, 10))\n",
    "for fi, f in enumerate([f0, f1, f2, f3]):\n",
    "    error_gauss = np.zeros((N,))\n",
    "    error_simp = np.zeros(((N-1)//2,))\n",
    "    error_fejer = np.zeros((N,))\n",
    "    I_quad, _ = quad(f, -1, 1)\n",
    "    \n",
    "    for n in range(1, N+1):\n",
    "        I_gauss = gauss(f, n)\n",
    "        error_gauss[n-1]  = np.abs(I_gauss - I_quad)\n",
    "        I_fejer = fejer(f, n)\n",
    "        error_fejer[n-1]  = np.abs(I_fejer - I_quad)\n",
    "    \n",
    "    for n in range( 3, N+1, 2 ):\n",
    "        I_simp = simpson(f, n)\n",
    "        error_simp[(n-2)//2] = np.abs(I_simp - I_quad)\n",
    "    \n",
    "    ax = figure.add_subplot(2, 2, fi + 1)\n",
    "    ax.scatter(np.arange(1, N+1), np.log10(error_gauss, out = -16. * np.ones(error_gauss.shape), \n",
    "                                          where = (error_gauss > 1e-16)), label = 'Gauss', marker=\"+\")\n",
    "    ax.scatter(np.arange(3, N+1, 2), np.log10( error_simp ), label = 'Simpson', marker=\"+\")\n",
    "    ax.scatter(np.arange(1, N+1), np.log10(error_fejer, out = -16. * np.ones(error_fejer.shape), \n",
    "                                          where = (error_fejer > 1e-16)), label = 'Fejer', marker=\"+\")\n",
    "    ax.legend()\n",
    "    ax.set_title(f\"Erreur de différentes méthodes de quadrature pour {f.__name__}\")\n",
    "    ax.set_xlabel(\"n\")\n",
    "    ax.set_ylabel(\"log10(e)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Commentaires des résultats obtenus\n",
    "\n",
    "Globalement, pour nos quatres fonctions, on observe que plus N, le nombre de points augmente, plus nos trois méthodes sont précises car les erreurs diminuent. Cependant, en comparant les trois méthodes entres-elles, on remarque que, pour une fonction entière telle que f0 ou f1, la méthode de simpson est la moins précise des trois quelque surtout quand n augmente. Pour une fonction peu régulière comme f2 ou f3, les trois méthodes sont a peu près équivalentes mais restent moins précises que pour des fonctions régulières.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Degré de précision de la formule de Fejér\n",
    "\n",
    "Nous allons maintenant trouver numériquement le dégré de précision de la formule de Fejer.\n",
    "\n",
    "On rappelle qu'une formule de quadrature $ J(f) $ est dite précise à l'ordre $ d $ si $ J( p ) = I( p ) $ pour tout $ p \\in \\mathbb{R}_d[ X ] $ et $ J( x^{ d + 1 } ) \\neq I( x^{ d + 1 } ) $, où $ I( p ) $ est la valeur exacte de l'intégrale de $ p $. \n",
    "\n",
    "Nous allons tout d'abord vérifier numériquement que la formule $ J_G $ sur $ [-1,1] $ est précise de dégré $ 2N - 1 $.\n",
    "\n",
    "*Les questions théoriques Question 9 et Question 11 sont facultatives (questions bonus). Vous pouvez utiliser les résultats de ces questions pour jusifier vos réponses aux Questions 10 et 12.*\n",
    "\n",
    "**Question 9 (bonus).**\n",
    "> 1. Expliquez pourquoi il suffit de montrer que $ J_G( x^k ) = I( x^k ) $, $ k = 0, 1, \\ldots, 2N-1 $ et $  J_G( x^{ 2N } ) \\neq I( x^{ 2N } ) $ pour montrer que $ J_G $ est précise de dégré $ 2N - 1 $.\n",
    "> 2. Montrer que $ P_{ 2k } $ est une fonction paire et $ P_{ 2k + 1 } $ est une fonction impaire, $ k = 0, 1, 2, \\ldots $, où $ P_N $ est le $ N $-ième polynome de Legendre. En déduire que $ \\forall N \\in \\mathbb{N} $ si $ x $ est racine de $ P_N $, alors $ -x $ est aussi racine de $ P_N $. En déduire que pour $ k $ impair on a toujours $ J_G( x^k ) = I( x^k ) = 0 $.\n",
    "\n",
    "**Question 10.**\n",
    "> 1. Pour les valeurs $k = \\lbrace 0, 2, 4, 6, 8, 10 \\rbrace$ calculez l'erreur d'approximation de l'intégrale $ I(x^k) $ par la formule de Gauss à $ N $ points avec $N =  1,2,3, \\ldots, 10$. Affichez l'erreur sous le format ci-dessous.\n",
    "> 2. Expliquez comment les résultats obtenus montrent numériquement que la formule de quadrature de Gauss à $ N $ points est précise de degré $ 2N - 1 $.\n",
    "\n",
    "\n",
    "```\n",
    "------------------------------------------------------------------------------\n",
    "     N        x^0      x^2      x^4      x^6      x^8      x^10      \n",
    "------------------------------------------------------------------------------\n",
    "     1        0.000    0.667    0.400    0.286    0.222    0.182  \n",
    "     2        0.000    0.000    0.178    0.212    0.198    0.174  \n",
    "    ...        ...      ...      ...      ...      ...      ...\n",
    "```\n",
    "\n",
    "**Question 11 (bonus).**\n",
    "> Montrez que pour les points de Tchebychev $ x_{k, N} $ on a $ x_{ N+1-k, N } = -x_{k, N} $. En déduire que pour $ k $ impair on a toujours $ J_F( x^k ) = I( x^k ) = 0 $.\n",
    "\n",
    "**Question 12.**\n",
    "> 1. Reprendre la question 10.1 pour la formule de Fejér. \n",
    "> 2. En utilisant les résultats de la question 12.1 trouver numériquement le degré de précision de la formule de Fejér. Comparer avec le degré de précision de la formule de Newton-Cotes à $N$ points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----------------------------------------------------------------------\n",
      "    N     |     x^0       x^2       x^4       x^6       x^8      x^10\n",
      "-----------------------------------------------------------------------\n",
      "    1     |    0.000     0.667     0.400     0.286     0.222     0.182    \n",
      "    2     |    0.000     0.000     0.178     0.212     0.198     0.174    \n",
      "    3     |    0.000     0.000     0.000     0.046     0.078     0.095    \n",
      "    4     |    0.000     0.000     0.000     0.000     0.012     0.026    \n",
      "    5     |    0.000     0.000     0.000     0.000     0.000     0.003    \n",
      "    6     |    0.000     0.000     0.000     0.000     0.000     0.000    \n",
      "    7     |    0.000     0.000     0.000     0.000     0.000     0.000    \n",
      "    8     |    0.000     0.000     0.000     0.000     0.000     0.000    \n",
      "    9     |    0.000     0.000     0.000     0.000     0.000     0.000    \n",
      "   10     |    0.000     0.000     0.000     0.000     0.000     0.000    \n"
     ]
    }
   ],
   "source": [
    "def f(x, k):\n",
    "    return x**k\n",
    "\n",
    "print(\"-----------------------------------------------------------------------\")\n",
    "print(\"{:>5s}     | {:>7s} {:>9s} {:>9s} {:>9s} {:>9s} {:>9s}\".format(\"N\", \"x^0\", \"x^2\", \"x^4\", \"x^6\", \"x^8\", \"x^10\"))\n",
    "print(\"-----------------------------------------------------------------------\")\n",
    "\n",
    "for N in range(1, 11):\n",
    "    approx_errors = []\n",
    "    for k in [x for x in range(0, 11, 2)]:\n",
    "        I_approx = gauss(lambda x: f(x, k), N)\n",
    "        I_exact = 2 / (k + 1) if k % 2 == 0 else 0\n",
    "        approx_error = np.abs(I_approx - I_exact)\n",
    "        approx_errors.append(approx_error)\n",
    "    print(\"{:5d}     |    \".format(N) + \" \".join(\"{:.3f}    \".format(e) for e in approx_errors))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Les résultats obtenus dans le tableau d'erreur montrent numériquement que la formule de quadrature de Gauss à $N$ points est précise de degré $2N-1$ car les erreurs d'approximation pour les polynômes de degré inférieur à $2N-1$ sont toutes nulles. En effet, les erreurs d'approximation pour les polynômes de degré $0, 2, \\ldots, 2N-2$ sont toutes nulles, comme on peut le voir dans le tableau. En revanche, les erreurs d'approximation pour les polynômes de degré $2N, 2N+2, \\ldots$ sont non nulles. Cela correspond exactement à la propriété de précision de la formule de quadrature de Gauss à $N$ points qui est de degré $2N-1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----------------------------------------------------------------------\n",
      "    N     |     x^0       x^2       x^4       x^6       x^8      x^10\n",
      "-----------------------------------------------------------------------\n",
      "    1     |    0.000     0.667     0.400     0.286     0.222     0.182    \n",
      "    2     |    0.000     0.333     0.100     0.036     0.097     0.119    \n",
      "    3     |    0.000     0.000     0.100     0.089     0.059     0.029    \n",
      "    4     |    0.000     0.000     0.017     0.048     0.059     0.058    \n",
      "    5     |    0.000     0.000     0.000     0.006     0.017     0.027    \n",
      "    6     |    0.000     0.000     0.000     0.002     0.005     0.010    \n",
      "    7     |    0.000     0.000     0.000     0.000     0.001     0.002    \n",
      "    8     |    0.000     0.000     0.000     0.000     0.000     0.001    \n",
      "    9     |    0.000     0.000     0.000     0.000     0.000     0.000    \n",
      "   10     |    0.000     0.000     0.000     0.000     0.000     0.000    \n"
     ]
    }
   ],
   "source": [
    "def f(x, k):\n",
    "    return x**k\n",
    "\n",
    "print(\"-----------------------------------------------------------------------\")\n",
    "print(\"{:>5s}     | {:>7s} {:>9s} {:>9s} {:>9s} {:>9s} {:>9s}\".format(\"N\", \"x^0\", \"x^2\", \"x^4\", \"x^6\", \"x^8\", \"x^10\"))\n",
    "print(\"-----------------------------------------------------------------------\")\n",
    "\n",
    "for N in range(1, 11):\n",
    "    approx_errors = []\n",
    "    for k in [x for x in range(0, 11, 2)]:\n",
    "        I_approx = fejer(lambda x: f(x, k), N)\n",
    "        I_exact = 2 / (k + 1) if k % 2 == 0 else 0\n",
    "        approx_error = np.abs(I_approx - I_exact)\n",
    "        approx_errors.append(approx_error)\n",
    "    print(\"{:5d}     |    \".format(N) + \" \".join(\"{:.3f}    \".format(e) for e in approx_errors))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dans les formules de Newton-Côtes fermé de pas, le degré de précision est $n + 3$ si $n$ est pair $n + 2$ si $n$ est impair.\n",
    "La formule de Fejér est une formule de quadrature numérique qui utilise une moyenne des polynômes de degré $n$ de la formule de quadrature de Gauss à $n$ points. Elle est précise pour des polynômes de degré $2n-1$ ou moins.\n",
    "\n",
    "En observant le tableau d'erreur donné, on peut voir que pour la formule de Fejér, l'erreur pour $x^0$ est nulle à chaque itération, tandis que l'erreur pour $x^{2n}$ diminue à chaque itération. Cela suggère que la formule de Fejér est précise pour les polynômes de degré $2n-1$ ou moins.\n",
    "\n",
    "D'autre part, la formule de Newton-Cotes à $n$ points est précise pour les polynômes de degré $n$ ou moins. Ainsi, pour la formule de Newton-Cotes à $n$ points, on peut voir que l'erreur diminue à chaque itération jusqu'à ce qu'elle atteigne zéro pour $x^0$ à la quatrième itération, pour $x^2$ à la cinquième itération, et ainsi de suite. Cela suggère que la formule de Newton-Cotes à $n$ points est précise pour les polynômes de degré $n$ ou moins.\n",
    "\n",
    "En comparant les deux, on peut voir que la formule de Fejér est plus précise que la formule de Newton-Cotes à $n$ points, car elle est précise pour les polynômes de degré $2n-1$, tandis que la formule de Newton-Cotes à $n$ points n'est précise que pour les polynômes de degré $n$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
