{"metadata": {"kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"}, "celltoolbar": "None", "language_info": {"pygments_lexer": "ipython3", "nbconvert_exporter": "python", "version": "3.5.5", "mimetype": "text/x-python", "codemirror_mode": {"version": 3, "name": "ipython"}, "file_extension": ".py", "name": "python"}}, "nbformat_minor": 2, "nbformat": 4, "cells": [{"metadata": {}, "cell_type": "code", "source": ["import numpy as np\n", "import scipy as sp\n", "\n", "import matplotlib.pyplot as plt\n", "from scipy.integrate import odeint\n", "from numpy.linalg import inv"], "outputs": [], "execution_count": 1}, {"metadata": {"deletable": false}, "cell_type": "markdown", "source": ["$\\newcommand{\\R}{\\mathbb{R}}$\n", "$\\newcommand{\\N}{\\mathbb{N}}$\n", "\n", "#  TP 3 : Sch\u00e9mas explicites pour les EDO\n", "\n", "L'objectif de ce TP est de mettre en oeuvre diff\u00e9rents sch\u00e9mas num\u00e9riques pour approcher la solution d'une \u00e9quation diff\u00e9rentielle ordinaire. Nous allons consid\u00e9rer diff\u00e9rents sch\u00e9mas num\u00e9riques, dont certains  nous n'avons pas encore vus en cours. L'objectif est de comprendre ce que c'est une solution approch\u00e9e donn\u00e9e par une m\u00e9thode num\u00e9rique, mais aussi de comprendre des notions comme l' erreur entre la solution exacte et la solution approch\u00e9e, la  convergence d'une m\u00e9thode num\u00e9rique et l' ordre de pr\u00e9cision d'une m\u00e9thode num\u00e9rique\n", "\n", "\n", "\n", "Dans ce TP on consid\u00e8re des \u00e9quations diff\u00e9rentielles ordinaires, ou des syst\u00e8mes d'\u00e9quations diff\u00e9rentielles ordinaires, de la forme\n", "\\begin{equation}\n", "\\label{EDO}\n", "\\left\\{\\begin{aligned}\n", "  y'(t) & =  f(t,y(t)), \\\\\n", "  y(t_0) & =  y_0,\n", "\\end{aligned} \\right.\n", "\\tag{EDO}\n", "\\end{equation}\n", "\n", "o\u00f9 $f$ : $I \\times \\R^n \\longrightarrow \\R^n$ est une\n", "fonction dans les conditions du th\u00e9or\u00e8me de Cauchy-Lipschitz, 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 (EDO) admet alors une unique solution maximale d\u00e9finie dans un intervalle ouvert $J\\subseteq I$.\n", "\n", "\n", "On souhaite calculer une solution approch\u00e9e de ce probl\u00e8me dans un intervalle de la forme $[t_0,t_f]=[t_0,t_0+T]\\subseteq J$, avec $T>0$. Pour ce faire, on se donne $N\\in\\N$ et on consid\u00e8re une discr\u00e9tisation de $[t_0,t_0+T]$ de pas $h=\\frac TN,$ donn\u00e9e par les points \n", "$$\n", "t_0< t_1< \\cdots <t_N=t_0+T,\\ \\ \\ t_n=nh,\\ n=0,\\dots,N,\n", "$$\n", "et l'on utilise un sch\u00e9ma num\u00e9rique afin de calculer une approximation $y^n$ de $y(t_n)$ pour $n=0,\\dots,N$. \n", "\n", "Les sch\u00e9mas explicites qui suivent sont fr\u00e9quemment utilis\u00e9s pour discr\u00e9tiser le probl\u00e8me (EDO). On pose $y^0=y_0,\\,$ et, pour $n=0,\\dots,N-1,$ \n", "\n", "\n", "- **Euler explicite** :   $y^{n+1}=y^{n}+h f(t_n,y^{n})$ o\u00f9   $f^{n}=f(t_n,y^{n})$ ;\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)$ o\u00f9 $p_1=f(t_n,y^{n})$ et $p_2=f(t_n+h,y^{n}+h p_1)$ ;\n", "\n", "\n"]}, {"metadata": {"deletable": false}, "cell_type": "markdown", "source": ["### Exercice 1\n", "**Q1)** \u00c9crire une fonction python de la forme `euler_exp(fct, t0, tf, y0, h)`\n", "\n", "prenant en argument  `fct` la fonction $f(t,y)$ de (EDO), 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 tableaux :  \n", "\n", "- `[t_0,\\, t_1,\\, ...,\\, t_N]`, tableau  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", "\n", "- `[y^0,\\, y^1,\\, ...,\\, y^N]` tableau de taille $(N+1)\\times n$ repr\u00e9sentant la solution approch\u00e9e aux instants $t_n,\\ n=0,\\dots,N.$"]}, {"metadata": {}, "cell_type": "code", "source": ["# <completer>\n", "\n", "\n"], "outputs": [], "execution_count": 2}, {"metadata": {"deletable": false}, "cell_type": "markdown", "source": ["**Q2)** Testez d'abord votre fonction `euler_exp` sur un exemple simple d'une EDO scalaire. Vous pouvez par exemple consid\u00e9rer le probl\u00e8me\n", "\\begin{equation}\\label{edoexp}\n", "\\begin{cases}\n", "y'(t)=y(t),\\\\\n", "y(0)=1,\n", "\\end{cases}\n", "\\end{equation}\n", "dont la solution est $y(t)=e^t.$ \n", "\n", "**Q2.a)** Tracer sur une m\u00eame fen\u00eatre graphique :\n", "\n", "-  La solution exacte sur l'intervalle $[0,1]$ discr\u00e9tis\u00e9 avec un pas de temps de $10^{-4}$ (autrement dit, en rempla\u00e7ant l'intervalle $[0,1]$ par l'ensemble discret correspondant aux points d'une sub-division de $[0,1]$ de pas $10^{-4}$) ;\n", "\n", "- Les solutions approch\u00e9es sur le m\u00eame intervalle obtenues avec la m\u00e9thode d'Euler avec des pas de temps $h$ \u00e9gaux \u00e0 $1/4,\\ 1/10,\\ 1/100.$ Observer que lorsque le pas diminue, la solution approch\u00e9e est  *plus proche* de la solution exacte."]}, {"metadata": {}, "cell_type": "code", "source": ["# <completer>\n", "\n"], "outputs": [], "execution_count": 3}, {"metadata": {"deletable": false}, "cell_type": "markdown", "source": ["**Q2.b)** Donner les valeurs approch\u00e9es de $y(1)=e$ obtenues avec $h=1/10,\\ 1/100,\\ 1/1000$, et donner les erreurs entre la valeur exacte de $e$ et les valeurs approch\u00e9es obtenues avec ces pas de temps. Observer la d\u00e9croissance de l'erreur lorsque le pas de temps diminue.\n", "\n", "**Q2.c)** Tracer, sur une autre figure, la diff\u00e9rence en valeur absolue entre la solution exacte et la solution approch\u00e9e, en fonction du temps **(i.e. $|y(t_n)-y^n|$ en fonction de $t_n$)**, en utilisant des discr\u00e9tisations de $[0,1]$ avec des pas de temps $h$ \u00e9gaux \u00e0 $1/10,\\ 1/100,\\ 1/1000.$ **Vous devez alors tracer trois courbes et observer que lorsque le pas diminue, cette diff\u00e9rence s'approche de plus en plus de 0.**"]}, {"metadata": {}, "cell_type": "code", "source": ["# <completer>\n"], "outputs": [], "execution_count": 47}, {"metadata": {"deletable": false}, "cell_type": "markdown", "source": ["**Q3)** Testez votre fonction `euler_exp` 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 la solution approch\u00e9e, obtenue avec le pas $h$ et avec un pas \u00e9gal \u00e0 $\\frac h2$."]}, {"metadata": {}, "cell_type": "code", "source": ["# <completer>\n"], "outputs": [], "execution_count": 48}, {"metadata": {"deletable": false}, "cell_type": "markdown", "source": ["**Q4)** Testez ensuite votre fonction 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", "$$\n", "y(t) = \\frac{1}{2} \\sin(t) t  + 5 \\cos(t) + \\sin(t),\n", "$$\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", "$$\n", "X'=F(t,X),\\ \\ \\ \\textrm{avec}\\ X=(u,v)=(y,y')^T.\n", "$$\n", "\n", "Repr\u00e9senter \u00e0 nouveau la solution exacte et la solution approch\u00e9e dans une m\u00eame fen\u00eatre graphique.\n", "\n", "**Remarque :** La solution $y$ de $(P_2)$ correspond \u00e0 la premi\u00e8re composante du vecteur $X$ ci-dessus. Votre fonction `euler_ex`! 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."]}, {"metadata": {}, "cell_type": "code", "source": ["# <completer>\n"], "outputs": [], "execution_count": 49}, {"metadata": {"deletable": false}, "cell_type": "markdown", "source": ["### Exercice 2 : \u00c9tude de l'erreur et convergence du sch\u00e9ma d'Euler explicite\n", "\n", " Dans cet exercice on essaie de comprendre ce que c'est l'erreur associ\u00e9 \u00e0 un sch\u00e9ma num\u00e9rique et ce que c'est un sch\u00e9ma num\u00e9rique convergent.\n", " \n", "On se donne un probl\u00e8me de Cauchy de la forme (EDO) et une m\u00e9thode num\u00e9rique pour approcher la solution de (EDO) 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", "\n", "\\begin{equation*}\n", "E_h =  \\max_{k=0,\\cdots,N}( |y(t_k)-y^{k}|). \n", "\\end{equation*}\n", "\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 1 :**  l'erreur globale $E$ d\u00e9pend du pas $h$, ou, de mani\u00e8re \u00e9quivalente, du nombre de points $N$ de la discr\u00e9tisation. Pour des valeurs de $N$ diff\u00e9rentes, ou de mani\u00e8re \u00e9quivalente pour des valeurs de $h$ diff\u00e9rentes, les  discr\u00e9tisations de l'intervalle $[t_0,t_0+T]$ sont diff\u00e9rentes (elles ont un nombre de points diff\u00e9rent), et les solutions approch\u00e9es sont diff\u00e9rentes. On s'attend \u00e0 que, lorsque $N$ augmente, ou de mani\u00e8re \u00e9quivalente lorsque $h$ diminue, l'erreur $E_h$ diminue, puisque l'on consid\u00e8re dans ce cas de plus en plus de points dans la discr\u00e9tisation de l'intervalle $[t_0,t_0+T]$ et les approximations que l'on a faites pour construire le sch\u00e9ma num\u00e9rique deviennent alors de plus en plus pr\u00e9cises.\n", "\n", "**Remarque 2 :** 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$"]}, {"metadata": {"deletable": false}, "cell_type": "markdown", "source": ["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", "**Q1)** 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$ ; repr\u00e9sentez dans la m\u00eame figure la diff\u00e9rence en valeur absolue entre la solution exacte et la solution approch\u00e9e, en fonction du temps, pour chaque valeur de $h$."]}, {"metadata": {}, "cell_type": "code", "source": ["# <completer>\n"], "outputs": [], "execution_count": 50}, {"metadata": {"deletable": false}, "cell_type": "markdown", "source": ["**Q2)** Calculer, pour chaque valeur de $h=1/2^s,\\ s = 1,2,...,8,$ l'erreur globale $E_h$ correspondante. V\u00e9rifiez que $E_h$ diminue et s'approche de $0$ lorsque $h$ diminue.\n", "\n", "**On dira qu'une m\u00e9thode num\u00e9rique converge si $\\displaystyle{\\lim_{h\\to0}E_h}=0$.**\n", "\n", "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", "\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' 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)$."]}, {"metadata": {}, "cell_type": "code", "source": ["# <completer>\n"], "outputs": [], "execution_count": 51}, {"metadata": {"deletable": false}, "cell_type": "markdown", "source": ["### Exercice 3 : Un sch\u00e9ma d'ordre 2\n", "\n", "**Q1)** \u00c0 l'image de ce que vous avez fait pour le sch\u00e9ma d'Euler explicite, \u00e9crire une fonction python de la forme `Heun(fct, t0, tf, y0, h)` prenant les m\u00eames arguments que la fonction `euler_exp` et retournant les m\u00eames tableaux. "]}, {"metadata": {}, "cell_type": "code", "source": ["# <completer>\n"], "outputs": [], "execution_count": 52}, {"metadata": {"deletable": false}, "cell_type": "markdown", "source": ["**Q2)** Pour le probl\u00e8me $(P_1),$ comparez les solutions approch\u00e9es obtenues avec la m\u00e9thode d'Euler explicite et avec la m\u00e9thode de Heun avec le m\u00eame pas de temps. Pour ce faire, repr\u00e9senter dans la m\u00eame fen\u00eatre graphique la solution exacte et les solutions approch\u00e9es obtenues avec chacune des m\u00e9thodes. Quelle m\u00e9thode semble \u00eatre plus pr\u00e9cise ? Pour le confirmer, repr\u00e9senter dans une autre figure la diff\u00e9rence, en valeur absolue, entre la solution exacte et la solution approch\u00e9e, pour les deux m\u00e9thodes."]}, {"metadata": {}, "cell_type": "code", "source": ["# <completer>\n"], "outputs": [], "execution_count": 53}, {"metadata": {"deletable": false}, "cell_type": "markdown", "source": ["**Q3)** Reprenez l'exercice **2** pour la m\u00e9thode de Heun. V\u00e9rifiez graphiquement que cette m\u00e9thode est d'ordre $2,$ en repr\u00e9sentant l'erreur globale en fonction du pas de discr\u00e9tisation en \u00e9chelle logarithmique et en estimant la pente de la droite passant au plus pr\u00eat des points correspondant \u00e0 cette repr\u00e9sentation ou en repr\u00e9sentant une droite de pente $2$ qui passe par un des points."]}, {"metadata": {}, "cell_type": "code", "source": ["# <completer>\n"], "outputs": [], "execution_count": 54}, {"metadata": {}, "cell_type": "code", "source": [], "outputs": [], "execution_count": null}, {"metadata": {}, "cell_type": "code", "source": [], "outputs": [], "execution_count": null}, {"metadata": {}, "cell_type": "code", "source": [], "outputs": [], "execution_count": null}]}