{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# L3 Analyse Numérique – TP d'évaluation\n",
    "\n",
    "### Fin de séance :\n",
    "A l'issue de la séance, cet énoncé-questionnaire nous sera remis sous le format `NomPrenom.ipynb` via un mail adressé à [Pierre Le Barbenchon](mailto:pierre.le-barbenchon@ens-rennes.fr) et à [Antoine Dequay](mailto:antoine.dequay@ens-rennes.fr), et ayant pour objet [TP Mag].\n",
    "\n",
    "\n",
    "### Instructions générales :\n",
    "Vous prendrez garde à rester concis dans vos programmes, à éviter les boucles inutiles et à incorporer les commentaires que vous jugerez utiles (à l'aide de `#` en Python).\n",
    "\n",
    "Les programmes seront rendus utilisables tels quels par le correcteur qui n'aura donc pas à définir lui-même de variables supplémentaires pour les tester. Les résultats numériques seront facilement retrouvés dans vos programmes. De même, les présentations graphiques seront idéalemet titrées et légendées (Rappel : `plt.title(\"Titre de la figure\")` et `plt.legend()` en ayant mis des `plt.plot(X,Y,label = \"fonction f\")`).\\\n",
    "Les programmes seront exempts de messages d'erreurs à l'éxecution. Dans le cas de résultats abérrants à l'exécution n'ayant pas pu être débuggués, nous vous invitons à indiquer en commentaire votre reflexion sur la difficulté survenue et les moyens qui pourraient permettre de la résoudre.\\\n",
    "Les recherches internet et échanges d'information ne sont pas autorisés, vous pouvez bien entendu utiliser le mémo Python qui vous aurez préalablement imprimé.\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercice 0 : Choix personnel de représentation du produit matriciel\n",
    "\n",
    "Pour faire des produits matriciels entre une matrice et un vecteur, avez-vous choisi une représentation avec un vecteur 1D ou une matrice colonne 2D pour représenter le vecteur ? Autrement dit, vous écrivez \n",
    "\n",
    "```np.array([[1,2],[5,4]]).dot(np.array([[1],[1]]))```ou ```np.array([[1,2],[5,4]]).dot(np.array([1,1]))```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<i>Réponse :</i>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Une fois ce choix fait, tenez-vous y pendant tout le sujet !"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercice 1\n",
    "\n",
    "Dans tout cet exercice, on résolvera les systèmes avec la commande `np.linalg.solve`. On considère le problème de Dirichlet homogène pour le laplacien : étant donné $f:\\R\\rightarrow\\R$, une fonction régulière, il s'agit de déterminer une fonction régulière $u:[0,1]\\rightarrow\\R$ telle que :$$\\begin{cases}-u''=f\\\\ u(0)=u(1)=0\\end{cases}$$\n",
    "\n",
    "On veut approcher cette solution $u$ à l'aide d'un *schéma aux différences finies*. Etant donné un entier $N\\geq4$, on considère une grille de discrétisation régulière $x_j^N=hj$, où $j\\in\\Z$ et $h=\\frac1{N+1}$ est le pas de la grille. On propose deux schémas permettant de déterminer une approximation $u_j^N\\simeq u\\left(x_j^N\\right)$, $j=1,\\dots,N$, de la solution du problème.\n",
    "\n",
    "**Schéma 1** : Pour $2\\leq j \\leq N-1$ :\n",
    "$$\\left\\{ \\begin{array}{lllllllll}& && 2 u^N_1 &- &u^N_2 &=& h^2 f(x_1^N) \\\\ -&u^N_{j-1} &+& 2 u^N_j &- &u^N_{j+1} &=& h^2 f(x^N_j) \\\\ -&u^N_{N-1} &+& 2 u^N_N  && &=& h^2 f(x_N^N) \\end{array} \\right.$$\n",
    "\n",
    "**Schéma 2** : Pour $3\\leq j \\leq N-2$ :\n",
    "$$\\left\\{ \\begin{array}{llllllllllllllll}& & &  & &120 u^N_1 &- &60 u^N_2 & & &=& 5h^2 f(x_{0}^N) &+& 50 h^2 f(x_1^N) &+& 5 h^2 f(x_2^N) \\\\ & &-& 48u^N_{1} &+&102 u^N_2 & -& 48u^N_{3} & - & 3u^N_{4} &=& 8h^2 f(x_{1}^N) &+& 44 h^2 f(x_2^N) &+& 8 h^2 f(x_{3}^N) \\\\ -& 3u^N_{j-2} &-& 48u^N_{j-1} &+&102 u^N_j & -& 48u^N_{j+1} & - & 3u^N_{j+2} &=& 8h^2 f(x_{j-1}^N) &+& 44 h^2 f(x_j^N) &+& 8 h^2 f(x_{j+1}^N) \\\\ -& 3u^N_{N-3} &-& 48u^N_{N-2} &+&102 u^{N}_{N-1} & -& 48u^N_{N} &  &  &=& 8h^2 f(x_{N-2}^N) &+& 44 h^2 f(x_{N-1}^N) &+& 8 h^2 f(x_{N}^N) \\\\ &  &-&60 u^N_{N-1} &+&120 u^N_N & & & & &=& 5h^2 f(x_{N-1}^N) &+& 50 h^2 f(x_N^N) &+& 5 h^2 f(x_{N+1}^N) \\end{array} \\right.\n",
    "$$\n",
    "\n",
    "\n",
    "### Question 1\n",
    "\n",
    "Pour chacun de ces deux schémas, déterminer des matrices $D^N$, $S^N$ et un vecteur $f^N$ tels qu'ils s'écrivent de manière équivalente sous la forme $$D^N u^N = h^2 S^N f^N.$$\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vous pouvez écrire dans des cellules de texte ou nous rendre une copie avec vos noms et prénoms sur laquelle vous aurez déterminé les matrices voulues"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On note $E_N$ l'erreur numérique obtenue pour $N$ points :$$E_N := \\max_{1\\leq j \\leq N} |u^N_j - u(x_j^N)|.$$\n",
    "\n",
    "\n",
    "### Question 2\n",
    "\n",
    "La fonction $u(x) = x(1-x)$ est solution du problème pour la donnée $f(x)= - u''(x)$. On peut démontrer (on ne le demande pas) que l'erreur numérique est, pour cet exemple, rigoureusement nulle. Vérifier-le numériquement (calculer $E_4$ pour chaque schéma)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Schéma 1 : $E_4 = ...$ <i>à remplir</i>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Schéma 2 : $E_4 = ...$ <i>à remplir</i>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 3\n",
    "\n",
    "On considère à présent $u(x) = x(1-x)e^{-4 \\cos (41 x)}$ et la donnée correspondante $f(x) = -u''(x)$. Avec ces nouvelles données, comparer numériquement l'erreur de convergence de ces deux schémas : calculer $E_{100}$ pour chaque schéma et commentez les résultats.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Schéma 1 : $E_{100} = ...$ <i>à remplir</i>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Schéma 2 : $E_{100} = ...$ <i>à remplir</i>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 4\n",
    "\n",
    "Toujours pour ce même exemple, déterminer numériquement l'ordre de convergence de chacun des deux schémas, et expliquer en quelques mots la méthodologie qui vous a permis de l'obtenir. On précise que dans ce contexte :\n",
    "$$\\textrm{ ordre } = \\max \\{ \\alpha\\in \\mathbb{R} \\ | \\ \\exists C>0, \\forall N\\geq 4, \\ E_N \\leq C h^{\\alpha}. \\}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercice 2\n",
    "\n",
    "Soit $f$ une fonction à  valeurs réelles définie sur $[0,1]$. Afin d'approcher l'intégrale $I_f = \\int_0^1 f(x)\\,dx$, on utilise la méthode de *quadrature des trapèzes* décrite ci-après.\n",
    "\n",
    "Pour tout entier naturel non-nul $N$, on note $h=1/N$ ainsi que $x_i= ih$ pour $i=0,\\ldots,N$, qui constitue une subdivision uniforme de $[0,1]$. On considère alors :\n",
    "\\begin{equation*}\n",
    " T_f(h) = \\sum_{i=0}^{N-1} h\\dfrac{f(x_i)+f(x_{i+1})}{2} = h\\left(\\dfrac{f(x_0)}{2}+f(x_1)+\\ldots+f(x_{N-1})+\\dfrac{f(x_N)}{2}\\right).\n",
    "\\end{equation*}\n",
    "\n",
    "\n",
    "### Question 1\n",
    "\n",
    "Déterminer par cette méthode la valeur approchée pour $N=16$ de $I=\\int_0^1 \\mathrm{e}^x \\,dx$ et $J=\\int_0^1 \\mathrm{e}^{-10x(1-x)}(1+\\cos( 100x)) \\,dx$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La méthode précédente converge au moins à l'ordre 2 au sens où pour toute fonction $f\\in\\mathcal{C}^{2}([0,1])$, il existe une constante $C>0$ telle que pour tout $N$ on ait l'estimation d'erreur $\\left| I_f - T_f(h) \\right| \\leq C h^2$. Plus précisément, on peut démontrer (on ne le demande pas !) que si $f$ est assez régulière, alors il existe une constante $\\alpha$ telle que pour $h>0$ suffisamment petit : $$T_f(h) = I_f + \\alpha h^2 + o(h^2).$$\n",
    "\n",
    "\n",
    "### Question 2\n",
    "\n",
    "Soit deux réels $a$ et $b$. On pose la fonction $T_f^{(1)}$ définie pour tout $h>0$ par\n",
    "$$T_f^{(1)}(h) := aT_f(2h) + bT_f(h).$$ \n",
    "\n",
    "Déterminer par le calcul exact les deux réels $a$ et $b$, tels que pour $h>0$ suffisamment petit : $$T_f^{(1)}(h) = I_f + o(h^2).$$\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vous pouvez écrire dans des cellules de texte ou écrire sur la copie que vous avez utilisée précédemment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 3\n",
    "\n",
    "Puisque le calcul effectif de $T_f(2h)$ utilise uniquement des abscisses $x_i$ intervenant déjà dans le calcul de $T_f(h)$, déterminer la nouvelle approximation $T_f^{(1)}(h)$ de l'intégrale $J$ obtenue pour $N=16$, ceci en limitant le nombre d'évaluations de la fonction $f$ dans votre programme au strict minimum. Vous afficherez le nombre d'évaluations de f nécessaires à votre approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 4\n",
    "\n",
    "En utilisant en guise de comparaison la valeur approchée de $J$ obtenue via la fonction `scipy.integrate.quad(f,0,1)[0]` de Python (en mettant `import scipy` au début du fichier), identifier numériquement l'ordre de convergence de la méthode $T_f^{(1)}(h)$ ainsi définie. On expliquera en quelques mots la méthodologie employée."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dans la suite on envisage d'obtenir des méthodes plus efficaces. Pour ce faire, on s'appuie sur le théorème suivant, déduit directement de la formule d'Euler Mac-Laurin, et qui sera admis. Les coefficients réels non-nuls $b_i$ qui interviennent sont les nombres de Bernoulli dont aucune connaissance supplémentaire n'est nécessaire pour traiter la suite du sujet.\n",
    "\n",
    "**Théorème**\n",
    "Soient $m$ un entier naturel non-nul et $f\\in\\mathcal{C}^{2m+2}([0,1])$, alors pour tout entier naturel non-nul $N$ (avec $h=1/N$ le pas associé), il existe $\\xi\\in(0,1)$ tel que :\n",
    "\\begin{equation*}\n",
    " T_f(h) = I_f +\\sum_{i=1}^{m}\\dfrac{b_{i}}{(2i)!}h^{2i}\\left(f^{(2i-1)}(1)-f^{(2i-1)}(0)\\right) + \\dfrac{b_{m+1}}{(2m+2)!}h^{2m+2}f^{(2m+2)}(\\xi) .\n",
    "\\end{equation*}\n",
    "\n",
    "\n",
    "### Question 5\n",
    "\n",
    "Considérant $N=32$, utiliser la méthode $T_f(h)$ pour approcher l'intégrale $K:=\\displaystyle\\int_0^1 \\cos(20\\sin(\\pi x)) \\,dx$ suivante. Déterminer l'approximation et l'erreur numérique ($|K-T_f(h)|$) alors obtenue, en calculant la valeur approchée de $K$ via la fonction `scipy.integrate.quad(f,0,1)[0]` de Python. Commenter au vu du théorème précédent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $n$ un entier naturel. On recherche une nouvelle méthode de quadrature prennant la forme suivante :$$T_f^{(n)}(h) := \\sum_{j=0}^n c_{j} T_f(2^{n-j} h),$$\n",
    "où les coefficients $\\bold{c}=(c_j)_{0\\leq j\\leq n}$ sont à déterminer, et seront bien entendu indépendants de $f$ et de $h$.\n",
    "\n",
    "\n",
    "### Question 6\n",
    "\n",
    "Déterminer une matrice $\\bold{A}\\in\\mathrm{GL}_{n+1}(\\mathbb R)$ et un vecteur $\\bold{v}\\in\\mathbb R^{n+1}$ tels que sous la condition $\\bold{A}\\bold{c}=\\bold{v}$, pour toute fonction $f$ définie sur $[0,1]$ et suffisamment régulière et pour $h$ suffisamment petit : $$T_f^{(n)}(h)-I_f = O(h^{2n+2}).$$\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vous pouvez écrire dans des cellules de texte ou écrire sur la copie que vous avez utilisée précédemment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 7\n",
    "\n",
    "Pour $N=2048$, utiliser $T_f^{(4)}(h)$ pour approcher de nouveau l'intégrale $J$, tout en limitant comme à la question 3, le nombre d'évaluations de $f$ au strict minimum. Quelle approximation obtenez-vous ? Comparez-la au résultat obtenu à l'aide de `scipy.integrate` et affichez le nombre d'évaluations de $f$ nécessaires à votre approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.9.13"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}