{"cells": [{"cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": ["import numpy as np\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt"]}, {"cell_type": "markdown", "metadata": {"deletable": false}, "source": ["# Cours acc\u00e9l\u00e9r\u00e9 analyse num\u00e9rique - M2 AMS - 2023/2024\n", "\n", "## Tout d'abord : Notebooks \n", "\n", "Ce fichier est un notebook. Il permet d'\u00e9crire du texte math\u00e9matique et d'\u00e9crire du code python. Il s'ouvre avec jupyter.\n", "\n", "Pour lancer Jupyter, il faut, dans un terminal, lancer la commande\n", "\n", "`jupyter-notebook`\n", "\n", "Cela ouvrira automatiquement un navigateur Web dans lequel on peut travailler. L'onglet principal repr\u00e9sente l\u2019arborescence des fichiers depuis le fichier o\u00f9 la commande \u00e0 \u00e9t\u00e9 lanc\u00e9.\n", "\n", "Les notebooks sont compos\u00e9s de cellules contenant du code (en python) ou du texte (simple ou format\u00e9 avec les balisages Markdown). Ces notebooks permettent de faire des calculs interactifs en Python, et constituent un outil de choix pour l'enseignement.\n", "\n", "On peut \u00e9diter une cellule en double-cliquant dessus, et l\u2019\u00e9valuer en tapant **Ctrl+Entr\u00e9e** (on utilisera aussi souvent **Maj.+Entr\u00e9e** pour \u00e9valuer et passer \u00e0 la cellule suivante). Les boutons dans la barre d\u2019outil vous seront tr\u00e8s utiles, survolez-les pour faire appara\u00eetre une infobulle si leur pictogramme n\u2019est pas assez clair. N\u2019oubliez pas de sauvegarder de temps en temps votre travail, m\u00eame si Jupyter fait des sauvegardes automatiques r\u00e9guli\u00e8res.\n", "\n"]}, {"cell_type": "markdown", "metadata": {"deletable": false}, "source": ["### Cellule Markdown\n", "Le *Markdown* est un format de texte qui accepte un minimum de formattage. Il permet de rapidement :\n", "- faire des liste.\n", "- faire des *italiques*, du **gras**.\n", "- faire des [liens](https://fr.wikipedia.org/wiki/Markdown)\n", "- \u00e9crire des formules math\u00e9matiques avec  $\\LaTeX$ ."]}, {"cell_type": "markdown", "metadata": {"deletable": false}, "source": ["### Cellule Python\n", "La cellule suivante est une cellule python. On peut l'\u00e9valuer avec **Maj.+Entr\u00e9e**."]}, {"cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["a = 8 et 2*a = 16\n"]}], "source": ["# Ceci est un commentaire Python\n", "# On prendra l'habitude de commenter son code :\n", "# Un code doit \u00eatre lisible et compr\u00e9hensible par d'autres utilisateurs\n", "#\n", "a = 8\n", "print(\"a =\", a, \"et 2*a =\", 2*a)"]}, {"cell_type": "markdown", "metadata": {"deletable": false}, "source": ["$$\\newcommand{\\tn}{\\mathrm}\n", "\\newcommand{\\eref}[1]{(\\ref{#1})}\n", "\\newcommand{\\dep}[2]{\\partial_{#2}{#1}}\n", "\\newcommand{\\dis}{\\displaystyle}\n", "\\newcommand{\\R}{\\mathbb{R}}\n", "\\newcommand{\\C}{\\mathbb{C}}\n", "\\newcommand{\\Z}{\\mathbb{Z}}\n", "\\newcommand{\\N}{\\mathbb{N}}\n", "$$\n", "# Sch\u00e9mas de diff\u00e9rences finies pour les \u00e9quations diff\u00e9rentielles.\n", "\n", "## L'approximation num\u00e9rique des \u00e9quations diff\u00e9rentielles.\n", "\n", "On s'int\u00e9resse \u00e0 l'approximation num\u00e9rique des \u00e9quations diff\u00e9rentielles, ordinaires et aux d\u00e9riv\u00e9es partielles. C'est surtout les deuxi\u00e8mes qui nous int\u00e9resseront mais pour introduire la m\u00e9thode des diff\u00e9rences finies, on commencera par les premi\u00e8res.\n", "\n", "Les \u00e9quations diff\u00e9rentielles correspondent souvent \u00e0 des mod\u00e8les math\u00e9matiques d\u00e9crivant des ph\u00e9nom\u00e8nes physiques. Dans beaucoup de cas, on ne sait pas calculer une solution explicite d'une \u00e9quation donn\u00e9e, et on doit utiliser des techniques de r\u00e9solution approch\u00e9e. On s'int\u00e9resse alors \u00e0 la **discr\u00e9tisation** et \u00e0 la **r\u00e9solution num\u00e9rique** de l'\u00e9quation, autrement dit on remplace l'\u00e9quation, qui est pos\u00e9e sur un domaine continu, par un probl\u00e8me discret, en discr\u00e9tisant le domaine. Si l'\u00e9quation originale est lin\u00e9aire on obtient, en discr\u00e9tisant, un syst\u00e8me lin\u00e9aire, si l'\u00e9quation originale est non lin\u00e9aire, on peut obtenir une \u00e9quation non lin\u00e9aire \u00e0 r\u00e9soudre, par exemple par une m\u00e9thode de r\u00e9solution approch\u00e9e d'\u00e9quations non lin\u00e9aires comme la m\u00e9thode de Newton.   \n", "Les principales m\u00e9thodes de discr\u00e9tisation pour une \u00e9quation diff\u00e9rentielle sont les m\u00e9thodes des diff\u00e9rences finies, la m\u00e9thode des \u00e9l\u00e9ments finis et la m\u00e9thode des volumes finis. Il y a aussi d'autres classes de m\u00e9thodes, comme les m\u00e9thodes spectrales. Au d\u00e9but de ce cours on s'int\u00e9resse \u00e0 la m\u00e9thode des diff\u00e9rences finies.\n", "\n", "### La m\u00e9thode des diff\u00e9rences finies - approximation des d\u00e9riv\u00e9es.\n", "\n", "On consid\u00e8re une \u00e9quation diff\u00e9rentielle, dont l'inconnue est une fonction $u$, pos\u00e9e sur un domaine physique $\\Omega\\subseteq\\R^d,\\ d\\geq1$.\n", "\n", "Le principe de la m\u00e9thode des diff\u00e9rences finies est d'approcher la solution $u$ de cette \u00e9quation en un ensemble discret (mais **grand**) $\\{x_1,\\dots,x_N\\}$ de points du domaine $\\Omega,$ en rempla\u00e7ant les d\u00e9riv\u00e9s de $u$ aux points $x_i$ par des quotients de diff\u00e9rences faisant intervenir les points voisins de $x_i$. Cette approche sera valable si la solution du probl\u00e8me est r\u00e9guli\u00e8re.\n", "\n", "Par exemple, on a pour une fonction $u$ d'une seule variable que sa d\u00e9riv\u00e9 en un point $x$ de son domaine v\u00e9rifie\n", "$$\n", "u'(x)=\\lim_{h\\to0}\\frac{u(x+h)-u(x)}{h}=\\lim_{h\\to0}\\frac{u(x)-u(x-h)}{h}=\\lim_{h\\to0}\\frac{u(x+h)-u(x-h)}{2h}=\\dots,\n", "$$\n", "et on peut donc approcher $u'(x)$ par les quotients\n", "$$\n", "\\frac{u(x+h)-u(x)}{h},\\ \\frac{u(x)-u(x-h)}{h},\\ \\frac{u(x+h)-u(x-h)}{2h},\\dots,\n", "$$\n", "avec $h\\,$ **petit**.\n", "\n", "Ces quotients de diff\u00e9rences finies font intervenir des valeurs de $u$ en un nombre fini (dans chacun de ces exemples, deux) de points du domaine proche de $x.$ On va v\u00e9rifier un peu plus en bas que l'on peut aussi approcher $u'(x)$ par exemple par la quantit\u00e9\n", "$$\n", "\\frac{2u(x+h)+3u(x)-6u(x-h)+u(x-2h)}{6h},\n", "$$\n", "pour $h$ petit.\n", "\n", "On a aussi que si $u$ est deux fois d\u00e9rivable au point $x,$ alors\n", "$$\n", "u''(x)=\\lim_{h\\to0}\\frac{u(x+h)-2u(x)+u(x-h)}{h^2}\n", "$$\n", "et on peut donc aussi approcher $u''(x)$ par le quotient de diff\u00e9rences finies\n", "$$\n", "\\frac{u(x+h)-2u(x)+u(x-h)}{h^2},\n", "$$\n", "avec $h$ **petit**.\n", "\n", "Une question que l'on se posera souvent dans le cadre de la m\u00e9thode des diff\u00e9rences finies est celle de l'**ordre** de l'approximation. On va expliciter cette notion plus loin, mais essayons de la comprendre d\u00e8s maintenant. Elle mesure l'erreur commise dans ces approximations.\n", "\n"]}, {"cell_type": "markdown", "metadata": {"deletable": false}, "source": ["Supposons $a,\\ b\\in\\R$ et $u:[a,b]\\longrightarrow\\R$ une fonction de classe $\\mathcal C^4$.\n", "Soit $x\\in]a,b[$.\n", "\n", "**Exercice 1**\n", "\n", "1. Montrer, en utilisant des d\u00e9veloppements de Taylor de $u$ autour de $x$, que, lorsque $h\\to0$,\n", "  $$\n", "\\left|u'(x)-\\frac{u(x+h)-u(x)}{h}\\right|=\\mathcal{O}(h),\\ \\ \\ \\left|u'(x)-\\frac{u(x)-u(x-h)}{h}\\right|=\\mathcal{O}(h), \n", "  $$\n", "  $$\n", "  \\left|u'(x)-\\frac{u(x+h)-u(x-h)}{2h}\\right|=\\mathcal{O}(h^2).\n", "  $$\n", "  et que \n", "  $$\n", "  \\left|u'(x)-\\frac{2u(x+h)+3u(x)-6u(x-h)+u(x-2h)}{6h}\\right|=\\mathcal{O}(h^3).\n", "  $$\n", "\n", "2. *Construction d'une approximation de diff\u00e9rences finies \u00e0 un ordre pr\u00e9cis.* Chercher une approximation par diff\u00e9rences finies de $u'(x)$ bas\u00e9e sur les points $x,\\ x-h$ et $x-2h$ de la forme $au(x)+bu(x-h)+c u(x-2h)$ tel que\n", "  $$\n", "  \\left|u'(x)-\\big(au(x)+bu(x-h)+c u(x-2h)\\big)\\right|=\\mathcal{O}(h^2).\n", "  $$\n", "\n", "Les diff\u00e9rences ci-dessus s'appellent **erreurs de troncature** ou de **consistance** au point $x$. On dira que cette erreur est d'ordre 1 si elle se comporte en $\\mathcal{O}(h)$, d'ordre $p>0$ si elle se comporte en $\\mathcal{O}(h^p).$ Plus $p$ est grand, plus l'approximation choisie de la d\u00e9riv\u00e9e de $u$ est pr\u00e9cise.\n", "\n", "L'exercice 1.2 nous fait penser que l'on peut approcher la d\u00e9riv\u00e9e d'une fonction $u$ en un point $x$ par une formule de diff\u00e9rences finies \u00e0 un ordre souhait\u00e9, en utilisant des valeurs de $u$ en points autour de $x$, en cherchant les bons poids pour chaque valeur. Mais on voit facilement que le plus \u00e9lev\u00e9 est l'ordre souhait\u00e9, le plus de points il est n\u00e9c\u00e9ssaire d'utiliser.\n", "\n", "**Exercice 2. Interpolation polynomial**\n", "\n", "On peut approcher la d\u00e9riv\u00e9e d'une fonction $u$ en un point $x$ en approchant la fonction $u$ par un polyn\u00f4me $X\\mapsto P(X)$ et en calculant $P'(x)$. V\u00e9rifier que si l'on approche $u$ par ses polyn\u00f4mes d'interpolation de d\u00e9gr\u00e9 1 aux points $x$ et $x-h$, ou $x$ et $x+h$, ou $x-h$ et $x+h$, et que l'on calcule les d\u00e9riv\u00e9es de ces p\u00f4lyn\u00f4mes au point $x$, on obtient les formules de diff\u00e9rences finies correspondantes que l'on a vu pr\u00e9c\u00e9damment. On peut faire le m\u00eame exercice en approchant $u$ par le polyn\u00f4me d'interpolation (de degr\u00e9 2) aux points $x,\\ x-h$ et $x-2h$ et par le polyn\u00f4me (de degr\u00e9 3) qui interpole $u$ aux points $x,\\ x+h,\\ x-h$ et $x-2h$.\n", "\n", "\n", "#### Approximation des d\u00e9riv\u00e9es dans un ensemble de points d'un intervale.\n", "\n", "Soit $[a,b]\\subseteq\\R$ et $h>0$ tel qu'il existe $N\\in\\N$ tel que $\\frac{b-a}{N}=h$. Consid\u00e9rons alors les $N+1$ points $x_n=a+nh,\\ n=0,\\dots,N$. L'ensemble discret form\u00e9 par les points ${x_0,\\dots,x_N}$ s'appelle une discr\u00e9tisation r\u00e9guli\u00e8re ou uniforme de $[a,b]$ de pas $h$.\n", "\n", "Soit $u$ une fonction r\u00e9guli\u00e8re d\u00e9finie sur $[a,b]$ et soit $U\\in\\R^{N+1}$ le vecteur $(U_0,\\dots,U_N)$ avec $U_n=u(x_n),\\ n=0,\\dots,N$.\n", "\n", "On souhaite ici approcher les d\u00e9riv\u00e9es de $u$ en chaque point $x_n$ de la discr\u00e9tisation de l'intervalle $[a,b]$ au m\u00eame ordre, en utilisant le m\u00eame type de diff\u00e9rences finies en chaque point. La difficult\u00e9 en faire ceci consiste dans l'approximation des d\u00e9riv\u00e9es aux points $\\{a,\\ b\\}$ du bord de l'intervalle et en g\u00e9n\u00e9ral dans ces points on doit utiliser une approximation diff\u00e9rente de celle utilis\u00e9e aux points de l'int\u00e9rieur du domaine. Ce probl\u00e8me se pose lors de la discr\u00e9tisation des conditions aux limites pour un probl\u00e8me aux limites pour une \u00e9quation aux d\u00e9riv\u00e9es partielles.\n", "\n", "**Exercice 3**\n", "1. D\u00e9terminer en fonction de $h$ et de $U$ un vecteur $V\\in\\R^{N+1}$ v\u00e9rifiant\n", "    $$\n", "\\max_{n=0,\\dots,N}|V_n-u'(x_n)|=\\mathcal{O}(h).\n", "    $$\n", "2. D\u00e9terminer en fonction de $h$ et de $U$ un vecteur $W\\in\\R^{N+1}$ v\u00e9rifiant\n", "    $$\n", "\\max_{n=0,\\dots,N}|W_n-u'(x_n)|=\\mathcal{O}(h^2).\n", "    $$\n", "3. **A faire dans la cellule code ci-dessous**.\\\\\n", "  Illustrer num\u00e9riquement les deux approximations de la d\u00e9riv\u00e9e de $u$ pour la fonction $u(x)=\\sin(x)$ dans l'intervalle $[a,b]=[0,2\\pi]$. Pour ce faire repr\u00e9senter dans la m\u00eame figure la fonction $u'$, le vecteur $V$ et le vecteur $W$, en fonction de $x$, avec $h=\\frac{2\\pi}{10}$ et puis avec $h=\\frac{2\\pi}{100}$. Repr\u00e9senter dans une autre figure les diff\u00e9rences $|u'-V|$ et $|u'-W|$. \n", "\n", "On peut construire une discr\u00e9tisation non uniforme de l'intervale $[a,b]$, c'est-\u00e0-dire telle que la diff\u00e9rence entre chaque deux points cons\u00e9cutifs de la discr\u00e9tisation ne soit pas toujours \u00e9gale. Dans ce cas on appelera pas de la discr\u00e9tisation la plus grande distance entre chaque deux points cons\u00e9cutifs de la discr\u00e9tisation."]}, {"cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": ["#CODER ICI\n", "# <completer>\n"]}, {"cell_type": "markdown", "metadata": {"deletable": false}, "source": ["#### Analyse de l'ordre de pr\u00e9cision.\n", "\n", "La m\u00e9thode des diff\u00e9rences finies pour approcher la solution $u$ d'une \u00e9quation diff\u00e9rentielle pos\u00e9e dans un domaine $\\Omega$ consiste \u00e0 remplacer l'\u00e9quation diff\u00e9rentielle par une \u00e9quation aux diff\u00e9rences finies, en rempla\u00e7ant les d\u00e9riv\u00e9es de $u$ par des quocients de diff\u00e9rences finies en des points obtenues en discr\u00e9tisant le domaine $\\Omega$ avec un **pas** $h>0$.\n", "\n", "Nous avons vu qu'en utilisant des d\u00e9veloppements de Taylor relativement simples, on peut obtenir l'ordre de l'erreur commmise entre la d\u00e9riv\u00e9e et son approximation par diff\u00e9rences finies. Dans les exemples que l'on a vu, cette erreur se comporte en $\\mathcal{O}(h^p)$, avec $p\\in\\N$. En g\u00e9neral on ne sait pas exactement comment se comportent ses erreurs. Une mani\u00e8re d'avoir une id\u00e9e de ce comportement est de repr\u00e9senter graphiquement la diff\u00e9rence entre la fonction et ses d\u00e9riv\u00e9es, pour des valeurs de $h$ diff\u00e9rentes, et de chercher \u00e0 comprendre comment se comporte cette erreur en fonction de $h$.\n", "\n", "  Plus pr\u00e9cis\u00e9ment, consid\u00e9rons toujours le cas d'une fonction $u$ d\u00e9finie sur un intervalle $[a,b]$ et, pour $h>0$, une discr\u00e9tisation de l'intervalle $[a,b]$ de pas $h$, donn\u00e9e par $N+1$ points $x_0,\\dots,x_N$ ($N$ d\u00e9pend de $h$).\n", "\n", "  Supposons de l'on approche la d\u00e9riv\u00e9e $u'(x_n)$ en chaque point $x_n$ de la discr\u00e9tisation par un quocient de diff\u00e9rences finies donn\u00e9e par une formule $Du(x_n).$ On note $E(h)$ l'erreur globale commise (par exemple en norme infinie) avec le pas $h$ : $E(h)=\\max_{n=0,\\dots,N}|u'(x_n)-Du(x_n)|$. Supposons que l'on s'attend \u00e0 ce que $E(h)$ se comporte comme $h^p$ pour un certain $p\\in\\N$, c'est-\u00e0-dire qu'il existe $C>0$ ind\u00e9pendant de $h$ tel que $E(h)\\sim C h^p$. Alors on s'attend \u00e0 ce que $\\log\\big(E(h)\\big)$ se comporte comme\n", "  $$\n", "\\log\\big(E(h)\\big)\\sim \\log(C)+p\\log(h).\n", "  $$\n", "On peut alors repr\u00e9senter dans un graphique $\\log(E(h)$ en fonction de $\\log(h)$, pour des valeurs de $h$ de plus en plus petits. Si le comportement est l'attendu, on obtiendra que, dans cette \u00e9chelle log-log, le graphique de $E(h)$ en fonction de $h$ est une droite de pente $p$, et cette pente de la droite nous donne l'ordre de la m\u00e9thode.\n", "\n", "**Exercice 4**\n", "1. Choisir un point quelconque $x$ et repr\u00e9senter, pour plusieurs valeurs de $h$ de plus en plus petites, la diff\u00e9rence en valeur absolue entre $u'(x)$ et, respectivement $\\frac{u(x+h)-u(x)}{h}$, $\\frac{u(x+h)-u(x-h)}{2h}$ et $\\frac{2u(x+h)+3u(x)-6u(x-h)+u(x-2h)}{6h}$, en fonction de $h$, en \u00e9ch\u00e8lle logarithmique. Repr\u00e9senter des droites de pente 1, 2 et 3 sur le m\u00eame graphique pour visualiser comment se comporte cette diff\u00e9rence en fonction de $h$.\n", "\n", "2. Faire le m\u00eame exercice mais en consid\u00e9rant maintenant des discr\u00e9tisations uniformes de pas $h$ de par exemple l'intervale $[0,\\pi]$, et calculer l'erreur globale des m\u00eames approximations sur cet intervale. Pour approcher les points de bord dans la derni\u00e8re formule on peut utiliser les valeurs exactes ici."]}, {"cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": ["# <completer>\n"]}, {"cell_type": "markdown", "metadata": {"deletable": false}, "source": ["## Approximation des EDO.\n", "\n", "On s'int\u00e9resse \u00e0 des m\u00e9thodes num\u00e9riques pour approcher les \u00e9quations diff\u00e9rentielles ordinaires.\\\\\n", "\n", "On consid\u00e8re le probl\u00e8me de Cauchy pour une \u00e9quation diff\\'erentielle ordinaire, ou pour un syst\\`eme d'\u00e9quations diff\\'erentielles ordinaires, de la forme\n", "\\begin{equation}\n", "\\label{sys}\n", "\\tag{sys}\n", "\\left\\{\\begin{aligned}\n", "  y'(t) & =  f(t,y(t)), \\\\\n", "  y(t_0) & =  y_0,\n", "\\end{aligned} \\right.\n", "\\end{equation}\n", "o\\`u $f$ : $I \\times \\R^n \\longrightarrow \\R^n$ est une\n", "fonction continue et localement Lipschitzienne par rapport \u00e0 sa deuxi\u00e8me variable, avec $I\\subseteq\\R$ un intervalle ouvert de $\\R$ et o\u00f9 $(t_0,y_0)\\in I\\times\\R^n$ est donn\u00e9. Le probl\u00e8me de Cauchy \\eqref{sys} admet alors une unique solution maximale d\u00e9finie dans un intervalle ouvert $J\\subseteq I$ : ce r\u00e9sultat est une cons\u00e9quence du **th\u00e9or\u00e8me de Cauchy-Lipschitz**.\n", "\n", "On souhaite calculer une solution approch\u00e9e (**discr\u00e8te**) de ce probl\u00e8me dans un intervalle de la forme $[t_0,t_f]=[t_0,t_0+T]\\subseteq J$, avec $T>0$.\n", "\n", "Pour cela on proc\u00e8de comme suit.\n", "\n", "1. On discr\u00e9tise (c'est-\u00e0-dire on remplace du continu par du discret) l'intervalle $[t_0,t_0+T]$. Pour ce faire on se donne $N\\in\\N$ et on d\u00e9finit une sub-division de $[t_0,t_0+T]$ en $N$ sous-intervalles d\u00e9finis par les $N+1$ points, appel\u00e9es points de la discr\u00e9tisation,\n", "$$\n", "t_0<t_1<\\dots<t_N=t_0+T.\n", "$$\n", "On appelle aussi les points $t_n$ les **instants temporels**, car pour de nombreuses EDO mod\u00e9lisant des ph\u00e9nom\u00e8nes physiques, la variable $t$ repr\u00e9sente le temps}.\n", "Le point $t_0$ est appel\u00e9 l'**instant initial** et le point $t_N$ l'**instant final**. On appelle aussi :\n", " - $h_n:=t_{n+1}-t_n$ le **pas de temps** ou **le pas de discr\u00e9tisation** **entre $t_n$ et $t_{n+1}$**, pour $n=0,\\dots,N-1$ ;\n", " - $h:=\\displaystyle{\\max_{n=0,\\dots,N-1} h_n}$ le **pas de temps** ou **pas de la discr\u00e9tisation**.\n", " \n", " La discr\u00e9tisation est dite **uniforme** si $h_n=h,\\ $ pour tout $n=0,\\dots,N-1.$ Dans ce cas la pas $h$ de la discr\u00e9tisation est d\u00e9fini par\n", "$$\n", "h=\\frac TN\n", "$$\n", "et les $N+1$ points de discr\u00e9tisation sont d\u00e9finis par \n", "$$\n", "t_n=t_0+nh,\\ \\ \\ n=0,\\dots,N.\n", "$$\n", "2. On construit $N+1$ valeurs $y^0,\\dots,y^N$ qui approchent la solution exacte $y$ de \\eqref{sys}, aux points $t_0,\\dots,t_N$, autrement dit on construit des valeurs $y^0,\\dots,y^N$ tels que\n", "  $$\n", "  y^n\\approx y(t_n),\\ \\textrm{pour }  n=0,\\dots,N.\n", "  $$\n", "  \n", "  \n", "Une **m\u00e9thode num\u00e9rique** ou un **sch\u00e9ma num\u00e9rique** est donc la donn\u00e9e d'une telle construction. Plus pr\u00e9cis\u00e9ment, c'est la donn\u00e9e d'une suite $(y^0,\\dots,y^N)_{N\\in\\N}$, index\u00e9e par $N\\in\\N$, de valeurs approchant les valeurs exactes $(y(t_0),\\dots,y(t_N))$ de la solution $y$ de \\eqref{sys} aux points de la discr\u00e9tisation $(t_0,\\dots,t_N)$. On appelle souvent la suite $(y^0,\\dots,y^N)_{N\\in\\N}$ de **solution approch\u00e9e** ou **solution num\u00e9rique** et, pour $N\\in\\N$ donn\u00e9e, le terme d'ordre $N$ de cette suite solution approch\u00e9e (ou solution num\u00e9rique) associ\u00e9e \u00e0 une discr\u00e9tisation de l'intervalle $[t_0,t_0+T]$ de $N+1$ points ou de pas $h=\\displaystyle{\\max_{0,\\dots,N}}h_n$.\n", "\n", "Comme la valeur de la solution exacte \u00e0 l'instant initial $y(t_0)=y_0$ est connue, on prendra $y^0=y_0$.\n", "\n", "La d\u00e9pendance de la suite $(y^0,\\dots,y^N)$ par rapport \u00e0 $N$ n'est pas toujours sp\u00e9cifi\u00e9e mais il faut la garder \u00e0 l'esprit, ainsi que la d\u00e9pendance de la suite $(t_0,\\dots,t_N)$ d\u00e9finissant les points de la discr\u00e9tisation de l'intervalle $[t_0,t_0+T]$.\n", "\n", "**Lorsque $N$ augment, le nombre de points de la discr\u00e9tisation augmente et on s'attend \u00e0 ce que les valeurs approch\u00e9es $y^0,\\dots,y^N$  approchent de mani\u00e8re de plus en plus pr\u00e9cise les valeurs exactes $y(t_0),\\dots,y(t_N)$.**\n", "\n", "Nous allons nous restreindre \u00e0 des discr\u00e9tisations uniformes de l'intervalle $[t_0,t_0+T]$ ({\\it le cas d'une discr\u00e9tisation \u00e0 pas $h_n$ variable est souvent utilis\u00e9 lorsque l'on veut adapter le pas de la discr\u00e9tisation pour am\u00e9liorer la pr\u00e9cision du sch\u00e9ma}). Dans ce cas, comme $h=\\frac TN$, on peut de mani\u00e8re \u00e9quivalente voir la suite approch\u00e9e $(y^0,\\dots,y^N)$ comme une suite index\u00e9e par $N\\in\\N$ ou comme une suite index\u00e9e par $h\\in]0,1]$ ; on a choisi l'intervalle $]0,1]$ pour $h$ mais on aurait pu choisir n'importe quel autre intervalle de la forme $]0,h_{max}]$. Ce qui nous int\u00e9ressera sera le comportement de la solution approch\u00e9e $(y^0,\\dots,y^N)$ lorsque $N$ tend vers l'infini ({\\it i.e., lorsque le nombre de points de la discr\u00e9tisation tend vers l'infini, ce qui correspond \u00e0 ce que l'ensemble discret qui repr\u00e9sente l'intervalle $[t_0,t_0+T]$ tend vers cet intervalle}), ou de mani\u00e8re \u00e9quivalente, lorsque $h$ tend vers 0 (**i.e. lorsque la distance entre deux points cons\u00e9cutifs de la discr\u00e9tisation tend vers 0**).\n", "\n", "Dans cette feuille on va consid\u00e9rer les sch\\'emas num\\'eriques suivants afin de calculer une approximation $y^n$ de $y(t_n)$ pour $n=0,\\dots,N$. \n", "\n", "On pose $y^0=y_0,\\,$ et, pour $n=0,\\dots,N-1,$ \n", "\n", "- **Euler explicite** :   $y^{n+1}=y^{n}+h f(t_n,y^{n})$|\n", "-  **Heun** :  $y^{n+1}=y^{n}+\\frac{h}{2}\\big(f(t_n,y^{n})+f\\big(t_n+h,y^{n}+h f(t_n,y^{n})\\big)\\big)$ \n", "- **Runge Kutta 4** :  $y^{n+1}=y^{n}+\\frac{h}{6}(p_1+2p_2+2p_3+p_4)$  o\u00f9\n", "\n", " $p_1=f(t_n,y^{n})$,\n", " \n", " $p_2=f(t_n+\\frac{h}{2},y^{n}+\\frac{h}{2} p_1)$,\n", " \n", " $p_3=f(t_n+\\frac{h}{2},y^{n}+\\frac{h}{2} p_2)$,\n", " \n", " $p_4=f(t_n+h,y^{n}+h p_3)$ \n", "\n", "Ces trois m\u00e9thodes sont des m\u00e9thodes explicites (la valeur $y^{n+1}$ se calcule explicitement \u00e0 partir de valeurs ant\u00e9rieures), \u00e0 un pas (la valeur $y^{n+1}$ se calcule uniquement \u00e0 partir de la valeur pr\u00e9c\u00e9dente $y^n$). Ces m\u00e9thodes peuvent s'\u00e9crire sous une forme g\u00e9n\u00e9rique\n", "$$\n", "(S)\\ \\ \\ \\ y^{n+1}=y^n+h\\Phi(t_n,y^n,h),\n", "$$\n", "pour $n=0,\\dots,N-1$, avec $y^0$ donn\u00e9 (nous consid\u00e9rerons toujours $y^0=y_0$),\n", "o\u00f9 $\\Phi:[t_0,t_0+T]\\times \\R^n\\times[0,1]\\longrightarrow\\R^n$ est une fonction continue. Par exemple, dans le cas de la m\u00e9thode d'Euler, $\\Phi(t,y,h)=f(t,y).$\n", "\n", "\n", "### Questions th\u00e9oriques - convergence d'une m\u00e9thode num\u00e9rique.\n", "\n", "La convergence d'une m\u00e9thode num\u00e9rique pour une EDO est li\u00e9e \u00e0 deux notions : la **consistance** et la **stabilit\u00e9**. Nous allons retrouver ces notions dans le cadre des m\u00e9thodes de diff\u00e9rences finies pour les EDP.\n", "\n", "Consid\u00e9rons une m\u00e9thode num\u00e9rique d\u00e9finie par $(S)$. Soit $y$ la solution exacte du probl\u00e8me de Cauchy \\eqref{sys}. \n", "On appelle **erreur locale \u00e0 l'instant $t_n$** de la m\u00e9thode, la quantit\u00e9\n", "$$\n", "e^n=y(t_n)-y^n.\n", "$$\n", "La m\u00e9thode d\u00e9finie par $(S)$ est dite **convergente** sur l'intervalle $[t_0,t_0+T]$ si pour toute donn\u00e9e initiale $y_0$, \n", "$$\n", "\\lim_{h\\to0}\\Big(\\max_{n=0,\\dots,N}\\|e^n\\|\\Big)=0,\n", "$$\n", "autrement dit si pour toute donn\u00e9e initiale $y_0$,\n", "$$\n", "\\lim_{h\\to0}\\Big(\\max_{n=0,\\dots,N}\\|y(t_n)-y^n\\|\\Big)=0.\n", "$$\n", "\n", "On remarque que la solution exacte $y$ de \\eqref{sys} v\u00e9rifie\n", "\\begin{equation*}\n", "y(t_{n+1})%&=y(t_n)+h\\Phi(t_n,y(t_n),h)+y(t_{n+1})-y(t_n)-h\\Phi(t_n,y(t_n),h)\\\\\n", "           %&\n", "            =y(t_n)+h\\Phi(t_n,y(t_n),h)+\\varepsilon^n,\\ \\ \\ n=0,\\dots,N-1,%\\\\ y(t_0)&=y_0,\n", "\\end{equation*}\n", "o\u00f9 la quantit\u00e9\n", "$$\n", "\\varepsilon^n=y(t_{n+1})-y(t_n)-h\\Phi(t_n,y(t_n),h)\n", "$$\n", "s'appelle **erreur de consistance** de la m\u00e9thode $(S)$ \u00e0 l'\u00e9tape $n$, relative \u00e0 la solution $y$. C'est l'erreur commise par une it\u00e9ration du sch\u00e9ma sur la solution exacte. La m\u00e9thode est dite **consistante** si la somme des erreurs de consistance est petite :\n", "$$\n", "\\lim_{h\\to0}\\sum_{n=0}^{N-1}\\|\\varepsilon^n\\|=0 ;\n", "$$\n", "la m\u00e9thode est dite **consistante d'ordre au moins $p$** s'il existe une constante $C>0$, ind\u00e9pendante de $h$, tel que \n", "$$\n", "\\forall\\ N\\in\\N,\\ \\ \\sum_{n=0}^{N-1}\\|\\varepsilon^n\\|\\leq Ch^p.\n", "$$\n", "\n", "**La m\u00e9thode est consistante si l'erreur commise par une \u00e9tape du sch\u00e9ma est petite car la discr\u00e9tisation est coh\u00e9rente avec l'EDO**.\n", "\n", "La m\u00e9thode est dite **stable** s'il existe une constante $M>0$, ind\u00e9pendante de $h$, tel que pour toutes suites $y^n$ et $\\widetilde{y}^n$ d\u00e9finies par\n", "\\begin{align*}\n", "  y^{n+1}&=y^n+h\\Phi(t_n,y^n,h),\\ \\ n=0,\\dots,N-1,\\ \\ \\ y^0\\ \\textrm{donn\\'e},\\\\\n", "  \\widetilde{y}^{n+1}&=\\widetilde{y}^n+h\\Phi(t_n,\\widetilde{y}^n,h)+\\rho^n,\\ \\ n=0,\\dots,N-1,\\ \\ \\ \\widetilde{y}^0\\ \\textrm{donn\\'e},\n", "\\end{align*}\n", "avec $\\rho^n\\in\\R,\\ n=0,\\dots,N-1$, on ait\n", "$$\n", "\\max_{n=0,\\dots,N}\\|y^n-\\widetilde{y}^n\\|\\leq M\\Big(\\|y^0-\\widetilde{y}^0\\|+\\sum_{n=0}^{N-1}\\|\\rho^n\\|\\Big).\n", "$$\n", "\n", "**Une m\u00e9thode stable est alors une m\u00e9thode telle que, lorsque l'on introduit une perturbation $\\rho^n$ \u00e0 chaque \u00e9tape du sch\u00e9ma, l'erreur totale de la m\u00e9thode est contr\u00f4l\u00e9e par le cumul des perturbations**.\n", "\n", "**Exercice 5.**\n", "  Montrer que si la m\u00e9thode est stable et consistante, alors elle est convergente. Montrer que si la m\u00e9thode est stable et consistante d'ordre au moins $p$ alors\n", "  il existe une constante $\\tilde{C}$ ind\u00e9pendante de $h$ tel que\n", "$$\n", "\\max_{n=0,\\dots,N}\\|y(t_n)-y^n\\|\\leq \\tilde{C}h^p.\n", "$$\n", "\n", "\n", "**On utilise souvent la formulation  \"m\u00e9thode d'ordre $p$\"  au lieu de \"m\u00e9thode consistante d'ordre $p$\" pour une m\u00e9thode dont l'erreur de consistance est d'ordre $p$. Ce dernier exercice montre que si la m\u00e9thode  est consistante d'ordre $p$, l'erreur globale du sch\u00e9ma se comporte aussi comme $\\mathcal{O}(h^p)$.**\n", "\n", "**Exercice 6**\n", "Supposons la solution $y$ du probl\u00e8me de Cauchy \\eqref{sys} r\u00e9guli\u00e8re. On peut montrer, en utilisant des d\u00e9veloppements de Taylor appropri\u00e9s, que la m\u00e9thode d'Euler explicite est consistante d'ordre 1 et que la m\u00e9thode de Heun est consistante d'ordre 2 (celle de RK4 est consistante d'ordre 4 mais les calculs sont plus compliqu\u00e9s).\n", "\n"]}, {"cell_type": "markdown", "metadata": {"deletable": false}, "source": ["**Exercice 7.  M\u00e9thodes d'Euler, Heun et RK4.**\n", "1. Pour chacune des trois m\u00e9thodes num\u00e9riques donn\u00e9es, \u00e9crire une fonction python de la forme\n", "   `meth(fct, t0, tf, y0, h)`\n", "\n", "  avec `meth = euler_exp, Heun` ou `RK4`,\n", " prenant en argument  `fct` la fonction $f(t,y)$, les extr\u00e9mit\u00e9s `t0` et `tf` de l'intervalle de temps  $[t_0,t_f]=[t_0,t_0+T],$ la donn\u00e9e initiale `y0` et le pas de temps `h`. Cette fonction devra retourner deux \n", " tableaux :  \n", " - $[t_0,\\, t_1,\\, ...,\\, t_N],$ tableau numpy unidimensionnel de taille $(N+1)\\times 1$ repr\u00e9sentant la subdivision de l'intervalle $[t_0,t_f]$ de pas $h$ consid\u00e9r\u00e9e, \n", " - $[y^0,\\, y^1,\\, ...,\\, y^N],$ tableau numpy de taille $(N+1)\\times n$ repr\u00e9sentant la solution approch\u00e9e aux instants $t_n,\\ n=0,\\dots,N$, donn\u00e9e par la m\u00e9thode choisie.\n", "2. Testez vos trois fonctions sur le mod\u00e8le logistique\n", "$$\n", "(P_1)\\ \\ \\ \\ \n", "\\begin{cases}\n", "y'(t)=c y (1 - \\frac{y}{b}),\\\\\n", "y(0)=a,\n", "\\end{cases}\n", "$$\n", "dont la solution exacte est\n", "$$\n", "y(t) = \\frac{b}{1 + \\frac{b-a}{a} e^{-ct}},\n", "$$\n", "avec les donn\u00e9es $c=1$, $b=2$, $a=0.1$, dans l'intervalle $[0,15]$, avec un pas $h=0.2$. Tracer sur la m\u00eame fen\u00eatre la solution exacte et les solutions approch\u00e9es, obtenue avec le pas $h$ et dans une autre fen\u00eatre avec un pas \u00e9gal \u00e0 $\\frac h2$.\n", "\n", "3. Testez ensuite vos fonctions dans le cas vectoriel $n>1$ ($n=2$) sur le probl\u00e8me\n", "$$\n", "(P_2)\\ \\ \\ \\\n", "\\begin{cases}\n", "y''(t)  =    -y(t) + \\cos(t)  \\\\\n", "y(0)  =   5,\\ \\ y'(0)  =   1,  \n", "\\end{cases}\n", "$$\n", "dont la solution exacte est \n", "$$\n", "y(t) = \\frac{1}{2} \\sin(t) t  + 5 \\cos(t) + \\sin(t),\n", "$$\n", "dans l'intervalle $[0,15]$, avec un pas  $h=0.2$.\\\\\n", "Pour ce faire, il faut \u00e9crire l'\u00e9quation d'ordre 2 de $(P_2)$ sous la forme d'un syst\u00e8me de deux \u00e9quations d'ordre 1 dans les nouvelles variables $u(t)=y(t)$ et $v(t)=y'(t)$. On se ram\u00e8nera alors \u00e0 la r\u00e9solution d'une \u00e9quation de la forme\n", "$$\n", "X'=F(t,X),\\ \\ \\ \\textrm{avec}\\ X=(u,v)=(y,y')^T.\n", "$$\n", " Repr\u00e9senter \u00e0 nouveau la solution exacte et les solutions approch\u00e9es dans une m\u00eame fen\u00eatre graphique.\n", " **La solution $y$ de $(P_2)$ correspond \u00e0 la premi\u00e8re composante du vecteur $X$ ci-dessus. Votre fonction `euler_exp` retournera dans ce cas, si vous avez respect\u00e9 la structure conseill\u00e9e, un tableau de taille $(N+1)\\times2$, $N$ \u00e9tant le nombre de points de la discr\u00e9tisation. Ce tableau donne les valeurs approch\u00e9es de $X$ au points de la discr\u00e9tisation consid\u00e9r\u00e9e, la premi\u00e8re colonne correspondant \u00e0 la premi\u00e8re composante de $X$, la seconde \u00e0 la seconde composante de $X$. La solution approch\u00e9e de $(P_2)$ que l'on cherche correspond alors \u00e0 la premi\u00e8re colonne de ce tableau.**\n"]}, {"cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": ["# <completer>\n"]}, {"cell_type": "markdown", "metadata": {"deletable": false}, "source": ["**Exercice 8**\n", "On se donne un probl\u00e8me de Cauchy de la forme **(sys)** et une m\u00e9thode num\u00e9rique pour approcher la solution de **(sys)** dans un intervalle de la forme $[t_0,t_0+T]$. Pour un certain pas de temps $h$ fix\u00e9 (ou pour un certain nombre de points de la discr\u00e9tisation $N$ fix\u00e9), l'erreur globale entre la solution approch\u00e9e associ\u00e9e \u00e0 une discr\u00e9tisation de pas $h$ de l'intervalle $[t_0,t_0+T]$ et la solution exacte est donn\u00e9e par :\n", "\\begin{equation*}\n", "E_h =  \\max_{k=0,\\cdots,N}( |y(t_k)-y^{k}|). \n", "\\end{equation*}\n", "Ci dessus, $y(t_k)$ est la solution exacte \u00e0 l'instant $t_k$ et $y^{k}$ la valeur approch\u00e9e de $y(t^k),$ donn\u00e9e par le sch\u00e9ma num\u00e9rique.\n", "\n", "**Remarque** : l'erreur globale $E$ d\u00e9pend du pas $h$, ou, de mani\u00e8re \u00e9quivalente, du nombre de points $N$ de la discr\u00e9tisation.  Parfois dans la litt\u00e9rature on ne sp\u00e9cifie pas la d\u00e9pendance de $E$ par rapport \u00e0 $h$, mais on doit toujours garder \u00e0 l'esprit cette d\u00e9pendance et que l'erreur est donc une fonction de $h$.\n", "\n", "Consid\u00e9rons le probl\u00e8me\n", "$$\n", "(P_3)\\ \\ \\ \\\n", "\\begin{cases}\n", "y'(t)=\\frac{\\cos(t)-y(t)}{1+t},\\\\\n", "y(0)=-\\frac14, \n", "\\end{cases}\n", "$$\n", "dont la solution exacte est \n", "$$\n", "y(t) = \\frac{\\sin(t)-1/4}{1 + t}.\n", "$$\n", "\n", "1. Calculez les solutions approch\u00e9es de $(P_3)$ obtenues avec le sch\u00e9ma d'Euler explicite, avec $h=1/2^s$ pour $s = 1,2,...,8$ ; \n", "  Calculer, pour chaque valeur de $h=1/2^s,\\ s = 1,2,...,8,$ l'erreur globale $E_h$ correspondante. Repr\u00e9sentez ensuite, en \u00e9chelle logarithmique, l'erreur en fonction du pas de temps $h$, autrement dit, repr\u00e9sentez $\\log(E_h)$ en fonction de $\\log(h)$. Vous devez obtenir des points qui sont \u00e0 peu pr\u00e8s align\u00e9s sur une droite de pente 1. V\u00e9rifiez graphiquement que c'est le cas, en estimant la pente de la droite passant au plus pr\u00eat des points (ou en repr\u00e9sentant une droite de pente $1$ qui passe par un des points et en v\u00e9rifiant que tous les points sont \u00e0 peu pr\u00e8s sur cette droite).\n", "  Ceci signifie que $\\log(E_h) \\sim C+\\log(h)$ et donc que $E_h\\sim \\widetilde{C}h$, pour certaines constantes $C$ et $\\widetilde{C}$. On dit que la m\u00e9thode d'Euler explicite est d'{\\color{blue} ordre} 1 : c'est l'ordre de la puissance de $h$ dans cette relation. On a donc que l'erreur globale $E_h$ tend vers 0 comme $h$. L'ordre d'une m\u00e9thode donne une indication sur sa vitesse de convergence. Une m\u00e9thode d'ordre $p$ est une m\u00e9thode dont l'erreur globale tend vers $0$ comme $h^p$. Donc plus l'ordre est \u00e9lev\u00e9, plus la m\u00e9thode converge plus vite.\n", "  \n", " **Remarque** : pour \u00e9tudier num\u00e9riquement l'ordre d'une m\u00e9thode, on utilise souvent l'\u00e9chelle logarithmique pour tracer l'erreur en fonction du pas de discr\u00e9tisation $h$. La pente de la droite obtenue donne l'ordre $p$ de la m\u00e9thode : si $E_h \\sim Ch^p$ alors $\\log(E_h)\\sim \\log(C) + p\\log(h)$."]}, {"cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": ["#Euler\n", "# <completer>\n"]}, {"cell_type": "markdown", "metadata": {"deletable": false}, "source": ["2. Vous pouvez refaire l'exercice pour la m\u00e9thode de Heun et de RK4. Vous allez conclure que la m\u00e9thode de Heun est d'ordre 2 et la m\u00e9thode de RK4 d'ordre 4.\n"]}, {"cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": ["# Heun\n", "# <completer>\n"]}, {"cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [], "source": ["# <completer>\n"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"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.8.16"}, "celltoolbar": "None"}, "nbformat": 4, "nbformat_minor": 2}