"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from IPython.core.display import HTML\n",
"css_file = './custom.css'\n",
"HTML(open(css_file, \"r\").read())"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Python version 3.8.10 (default, Nov 26 2021, 20:14:08) \n",
"[GCC 9.3.0]\n"
]
}
],
"source": [
"import sys #only needed to determine Python version number\n",
"print('Python version ' + sys.version)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"%reset -f\n",
"%matplotlib inline\n",
"\n",
"from matplotlib.pylab import *\n",
"from scipy.optimize import fsolve"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# M62_CM2 : schémas \"classiques\" à un pas"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Considérons le problème de Cauchy\n",
"\n",
"\n",
"trouver une fonction $y \\colon I\\subset \\mathbb{R} \\to \\mathbb{R}$ définie sur un intervalle $I=[t_0,T]$ telle que\n",
"$$\\begin{cases}\n",
"y'(t) = \\varphi(t,y(t)), &\\forall t \\in I=[t_0,T],\\\\\n",
"y(t_0) = y_0,\n",
"\\end{cases}$$\n",
"avec $y_0$ une valeur donnée et supposons que l'on ait montré l'existence et l'unicité d'une solution $y$ pour $t\\in I$.\n",
"
\n",
"\n",
"Pour $h>0$ soit $t_n\\equiv t_0+nh$ avec $n=0,1,2,\\dots,N$ une suite de $N+1$ nœuds de $I$ induisant une discrétisation de $I$ en $N$ sous-intervalles $I_n=[t_n;t_{n+1}]$ chacun de longueur $h=\\frac{T-t_0}{N}>0$ (appelé le *pas de discrétisation*).\n",
"\n",
"Pour chaque nœud $t_n$, on cherche la valeur inconnue $u_n$ qui approche la valeur exacte $y_n\\equiv y(t_n)$. \n",
"- L'ensemble de $N+1$ valeurs $\\{t_0, t_1=t_0+h,\\dots , t_{N}=T \\}$ représente les points de la *discrétisation*. \n",
"- L'ensemble de $N+1$ valeurs $\\{y_0, y_1,\\dots , y_{N} \\}$ représente la *solution exacte discrète*. \n",
"- L'ensemble de $N+1$ valeurs $\\{u_0 = y_0, u_1,\\dots , u_{N} \\}$ représente la *solution numérique* obtenue **en construisant une suite récurrente**."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Les schémas que l'on va construire permettent de calculer (explicitement ou implicitement) $u_{n+1}$ à partir de $u_n, u_{n-1}, ..., u_{n-k}$ et il est donc possible de calculer successivement $u_1$, $u_2$,..., en partant de $u_0$ par une formule de récurrence de la forme\n",
"$$\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"\\vdots\\\\\n",
"u_\\kappa=y_\\kappa,\\\\\n",
"u_{n+1}=\\Phi(u_{n+1},u_n, u_{n-1}, \\dots, u_{n-k}),&\\forall n=\\kappa,\\kappa+1,\\dots,N-1.\n",
"\\end{cases}$$\n",
"
\n",
"\n",
"**Méthodes explicites et méthodes implicites** \n",
"Une méthode est dite *explicite* si la valeur $u_{n+1}$ peut être calculée directement à l'aide des valeurs précédentes $u_k$, $k\\le n$ (ou d'une partie d'entre elles). \n",
"Une méthode est dite *implicite* si $u_{n+1}$ n'est défini que par une relation implicite faisant intervenir la fonction $\\varphi$.\n",
"\n",
"\n",
"**Méthodes à un pas et méthodes multi-pas** \n",
"Une méthode numérique pour l'approximation du problème de Cauchy est dite *à un pas* si pour tout $n\\in\\mathbb{N}$, $u_{n+1}$ ne dépend que de $u_n$ et éventuellement de lui-même. \n",
"Autrement, on dit que le schéma est une méthode *multi-pas* (ou à pas multiples)."
]
},
{
"cell_type": "markdown",
"metadata": {
"toc": true
},
"source": [
"
\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Construction de schémas à un pas\n",
"Si nous intégrons l'EDO $y'(t)=\\varphi(t,y(t))$ entre $t_n$ et $t_{n+1}$ nous obtenons\n",
"$$\n",
"\\int_{t_n}^{t_{n+1}}y'(t)\\mathrm{d}t=\\int_{t_n}^{t_{n+1}} \\varphi(t,y(t))\\mathrm{d}t\n",
"$$\n",
"c'est-à-dire\n",
"$$\n",
"y_{n+1}-y_n=\\int_{t_n}^{t_{n+1}} \\varphi(t,y(t))\\mathrm{d}t.\n",
"$$\n",
"On peut construire différents schémas selon la formule d'approximation utilisée pour approcher le membre de droite. \n",
"Cette solution approchée sera obtenue en construisant une suite récurrente comme suit:\n",
"\n",
"$$\\begin{cases}\n",
"u_0=y_0,\\\\\n",
"u_{n+1}=u_n+\\displaystyle\\int_{t_n}^{t_{n+1}} \\tilde f(t) \\mathrm{d}t \n",
"\\quad \\text{où $\\tilde f(t)$ est un polynôme interpolant }\\varphi(t,y(t)) \n",
"\\end{cases}\n",
"$$\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Schéma d'Euler explicite\n",
"Si on remplace une fonction $f$ par une constante égale à la valeur de $f$ à borne gauche de l'intervalle $[a;b]$ \n",
"(polynôme qui interpole $f$ au point $(a,f(a))$ et donc de degré $0$), \n",
"on a\n",
"\\begin{align*}\n",
"\\tilde f(x)&=f(a)\\\\\n",
"\\int_a^b f(x)\\mathrm{d}x &\\approx \\int_a^b \\tilde f(x)\\mathrm{d}x=(b-a)f(a).\n",
"\\end{align*}\n",
"Cette formule est dite *formule de quadrature du rectangle à gauche*. \n",
"\n",
"En utilisant cette formule pour approcher la fonction $t\\mapsto \\varphi(t,y(t))$ on a\n",
"$$\n",
"\\int_{t_n}^{t_{n+1}} \\varphi(t,y(t)) \\mathrm{d}t \\approx h \\varphi(t_n,y(t_n))\n",
"$$\n",
"et on reconnait le **schéma d'Euler progressif**\n",
"\n",
"$$\n",
"\\text{(EE)}\\quad\n",
"\\begin{cases}\n",
"u_0=y(t_0)=y_0,\\\\\n",
"u_{n+1}=u_n+h \\varphi(t_n,u_n)& n=0,1,2,\\dots N-1\n",
"\\end{cases}$$\n",
"
\n",
"Il s'agit d'un schéma à 1 pas explicite car il permet d'expliciter $u_{n+1}$ en fonction de $u_n$."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def EE(phi,tt,y0):\n",
" h = tt[1]-tt[0]\n",
" uu = [y0]\n",
" for i in range(len(tt)-1):\n",
" uu.append( uu[i]+h*phi(tt[i],uu[i]) )\n",
" return uu\n",
"\n",
"# def euler_progressif(phi,tt,y0):\n",
"# h = tt[1]-tt[0]\n",
"# N = len(tt)-1 # tt contient N+1 points numerotés de 0 à N\n",
"# uu = [y0 for i in range(N+1)]\n",
"# for i in range(N): # i = 0,1,...,N-1 \n",
"# uu[i+1] = uu[i]+h*phi(tt[i],uu[i])\n",
"# return uu"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Rappels sur la boucle.
\n",
" \n",
"- `len(tt)` = nombre d'éléments de la liste `tt` = $N+1$ \n",
"- les indices des éléments de `tt` vont de $0$ à $N$\n",
"- `range(M)` produit les nombres entiers de $0$ à $M-1$\n",
"\n",
"Conclusion : `range(len(tt)-1)` donne $0,1,2,\\dots,N-1$ comme souhaité\n",
" \n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Schéma d'Euler implicite\n",
"Si on remplace une fonction $f$ par une constante égale à la valeur de $f$ à borne droite de l'intervalle $[a;b]$ \n",
"(polynôme qui interpole $f$ au point $(b,f(b))$ et donc de degré $0$), on a\n",
"\\begin{align*}\n",
"\\tilde f(x)&=f(b)\\\\\n",
"\\int_a^b f(x)\\mathrm{d}x &\\approx \\int_a^b \\tilde f(x)\\mathrm{d}x=(b-a)f(b).\n",
"\\end{align*}\n",
"Cette formule est dite *formule de quadrature du rectangle à droite*. \n",
"\n",
"En utilisant cette formule pour approcher la fonction $t\\mapsto \\varphi(t,y(t))$ on a\n",
"$$\n",
"\\int_{t_n}^{t_{n+1}} \\varphi(t,y(t))\\mathrm{d}t\\approx (t_{n+1}-t_n) \\varphi(t_{n+1},y(t_{n+1})) = h \\varphi(t_{n+1},y(t_{n+1}))\n",
"$$\n",
"et on obtient le **schéma d'Euler rétrograde**\n",
"\n",
"$$\n",
"\\text{(EI)}\\quad\n",
"\\begin{cases}\n",
"u_0=y(t_0)=y_0,\\\\\n",
"u_{n+1}=u_n+h \\varphi(t_{n+1},u_{n+1})& n=0,1,2,\\dots N-1\n",
"\\end{cases}\n",
"$$\n",
"
\n",
"Il s'agit d'un schéma à 1 pas implicite car il ne permet pas de calculer directement $u_{n+1}$ en fonction de $u_n$ lorsque la fonction $\\varphi$ n'est pas triviale.\n",
"\n",
"Pour calculer $u_{n+1}$ il faudra utiliser un schéma pour le calcul du zéro d'une fonction quelconque. En effet, $u_{n+1}$ est solution de l'équation $x=u_n+h\\varphi(t_{n+1},x)$, c'est-à-dire un zéro de la fonction (en général non linéaire) $$x\\mapsto -x+u_n+h\\varphi(t_{n+1},x).$$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"def EI(phi,tt,y0):\n",
" h = tt[1]-tt[0]\n",
" uu = [y0]\n",
" for i in range(len(tt)-1):\n",
" uu.append( fsolve(lambda x: -x+uu[i]+h*phi(tt[i+1],x), uu[i]) )\n",
" return uu"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Schéma d'Euler modifié et RK1_M\n",
"Si on remplace une fonction $f$ par une constante égale à la valeur de $f$ au milieu de l'intervalle $[a;b]$ \n",
"(polynôme qui interpole $f$ au point $\\left(\\tfrac{a+b}{2},f\\left(\\tfrac{a+b}{2}\\right)\\right)$ et donc de degré $0$), on a\n",
"\\begin{align*}\n",
"\\tilde f(x)&=f\\left(\\tfrac{a+b}{2}\\right)\\\\\n",
"\\int_a^b f(x)\\mathrm{d}x &\\approx \\int_a^b \\tilde f(x)\\mathrm{d}x=(b-a)f\\left(\\tfrac{a+b}{2}\\right).\n",
"\\end{align*}\n",
"Cette formule est dite *formule de quadrature du rectangle* ou *du point milieu*. \n",
"\n",
"En utilisant cette formule pour approcher la fonction $t\\mapsto \\varphi(t,y(t))$ on a\n",
"$$\n",
"\\int_{t_n}^{t_{n+1}} \\varphi(t,y(t))\\mathrm{d}t\\approx h \\varphi\\left(t_n+\\frac{h}{2},y\\left(t_n+\\frac{h}{2}\\right)\\right)\n",
"$$\n",
"et on obtient \n",
"$$\\begin{cases}\n",
"u_0=y(t_0)=y_0,\\\\\n",
"u_{n+1}=u_n+h \\varphi\\left(t_n+\\frac{h}{2},u_{n+1/2}\\right)& n=0,1,2,\\dots N-1\n",
"\\end{cases}$$\n",
"où $u_{n+1/2}$ est une approximation de $y(t_n+h/2)$. \n",
"\n",
"Cependant, $u_{n+1/2}$ est une valeur inconnue, on doit alors en calculer une approximation $\\tilde u_{n+1/2}$.\n",
"Pour cela nous pouvons utiliser deux stratégies:\n",
"\n",
"1. Utiliser une prédiction d'Euler progressive sur un demi-pas: $\\tilde u_{n+1/2}=u_n+(h/2) \\varphi(t_{n},u_{n})$.\n",
" \n",
" On peut ainsi remplacer $\\varphi(t_{n}+h/2,u_{n+1/2})$ par $\\varphi(t_{n}+h/2,\\tilde u_{n+1/2})$. \n",
"\n",
" Nous avons construit ainsi un nouveau schéma appelé **schéma d'Euler modifié** qui s'écrit comme suit et qui est à 1 pas explicite car il permet d'expliciter $u_{n+1}$ en fonction de $u_n$.\n",
"\n",
"$$\n",
"\\text{(EM)}\\quad\n",
"\\begin{cases}\n",
"u_0=y(t_0)=y_0,\\\\\n",
"\\tilde u_{n+1/2}=u_n+\\frac{h}{2}\\varphi(t_{n},u_{n}),\\\\\n",
"u_{n+1}=u_n+h \\varphi\\left(t_n+\\frac{h}{2},\\tilde u_{n+1/2}\\right)& n=0,1,2,\\dots N-1\n",
"\\end{cases}$$\n",
"
\n",
"\n",
"1. Utiliser une moyenne: $\\tilde u_{n+1/2}=\\dfrac{u_n+u_{n+1}}{2}$.\n",
" \n",
" On peut ainsi remplacer $\\varphi(t_{n}+h/2,u_{n+1/2})$ par $\\varphi(t_{n}+h/2,\\tilde u_{n+1/2})$. \n",
"\n",
" Nous avons construit ainsi un nouveau schéma qu'on appelera **schéma RK1_M** qui s'écrit comme suit et qui est à 1 pas implicite car il ne permet pas d'expliciter $u_{n+1}$ en fonction de $u_n$. Pour calculer $u_{n+1}$ il faudra utiliser un schéma pour le calcul du zéro d'une fonction quelconque. En effet, $u_{n+1}$ est solution de l'équation $x=u_n+h\\varphi(t_{n}+\\frac{h}{2},\\frac{u_n+x}{2})$, c'est-à-dire un zéro de la fonction (en général non linéaire) $$x\\mapsto -x+u_n+h\\varphi(t_{n}+\\frac{h}{2},\\frac{u_n+x}{2}).$$\n",
"\n",
"$$\n",
"\\text{(RK1_M)}\\quad\n",
"\\begin{cases}\n",
"u_0=y(t_0)=y_0,\\\\\n",
"u_{n+1}=u_n+h \\varphi\\left(\\frac{t_n+t_{n+1}}{2},\\dfrac{u_n+u_{n+1}}{2}\\right)& n=0,1,2,\\dots N-1\n",
"\\end{cases}$$\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"def EM(phi,tt,y0):\n",
" h = tt[1]-tt[0]\n",
" uu = [y0]\n",
" for i in range(len(tt)-1):\n",
" u_tilde = uu[i]+0.5*h*phi(tt[i],uu[i])\n",
" uu.append( uu[i]+h*phi( tt[i]+h/2,u_tilde ) )\n",
" return uu\n",
"\n",
"def RK1_M(phi,tt,y0):\n",
" h = tt[1]-tt[0]\n",
" uu = [y0]\n",
" for i in range(len(tt)-1):\n",
" uu.append( fsolve(lambda x: -x+uu[i]+h*phi( (tt[i]+tt[i+1])/2,(uu[i]+x)/2 ), uu[i]) )\n",
" return uu"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Schéma du trapèze ou de Crank-Nicolson \n",
"Si on remplace une fonction $f$ par le segment qui relie $(a,f(a))$ à $(b,f(b))$ \n",
"(polynôme qui interpole $f$ en les points $(a,f(a))$ et $(b,f(b))$ et donc de degré $1$), on a\n",
"\\begin{align*}\n",
"\\tilde f(x)&=\\dfrac{f(b)-f(a)}{b-a}(x-a)+f(a)\\\\\n",
"\\int_a^b f(x)\\mathrm{d}x &\\approx \\int_a^b \\tilde f(x)\\mathrm{d}x=\\frac{b-a}{2}\\left(f(a)+f(b)\\right).\n",
"\\end{align*}\n",
"Cette formule est dite *formule de quadrature du trapèze*. \n",
"\n",
"En utilisant cette formule pour approcher la fonction $t\\mapsto \\varphi(t,y(t))$ on a\n",
"$$\n",
"\\int_{t_n}^{t_{n+1}} \\varphi(t,y(t))\\mathrm{d}t\\approx \\frac{h}{2}\\Big(\\varphi(t_{n},y(t_{n}))+\\varphi(t_{n+1},y(t_{n+1}))\\Big)\n",
"$$\n",
"et on obtient le **schéma du trapèze ou de Crank-Nicolson**\n",
"\n",
"$$\n",
"\\text{(CN)}\\quad\n",
"\\begin{cases}\n",
"u_0=y(t_0)=y_0,\\\\\n",
"u_{n+1}=u_n+\\dfrac{h}{2} \\Big(\\varphi(t_{n},u_{n}) + \\varphi(t_{n+1},u_{n+1})\\Big)& n=0,1,2,\\dots N-1\n",
"\\end{cases}\n",
"$$\n",
"
\n",
"En fait, ce schéma fait la moyenne des schémas d'Euler progressif et rétrograde.\n",
"\n",
"Il s'agit à nouveau d'un schéma à 1 pas implicite car il ne permet pas d'expliciter directement $u_{n+1}$ en fonction de $u_n$ lorsque la fonction $\\varphi$ n'est pas triviale. \n",
"Pour calculer $u_{n+1}$ il faudra utiliser un schéma pour le calcul du zéro d'une fonction quelconque. En effet, $u_{n+1}$ est solution de l'équation $x=u_n+\\frac{h}{2}\\left(\\varphi(t_n,u_n)+\\varphi(t_{n+1},x)\\right)$, c'est-à-dire un zéro de la fonction (en général non linéaire) $$x\\mapsto -x+u_n+\\frac{h}{2}\\Big(\\varphi(t_n,u_n)+\\varphi(t_{n+1},x)\\Big).$$"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"def CN(phi,tt,y0):\n",
" h = tt[1]-tt[0]\n",
" uu = [y0]\n",
" for i in range(len(tt)-1):\n",
" uu.append( fsolve(lambda x: -x+uu[i]+0.5*h*( phi(tt[i],uu[i])+phi(tt[i+1],x) ), uu[i]) )\n",
" return uu"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Schéma de Heun \n",
"Pour éviter le calcul implicite de $u_{n+1}$ dans le schéma du trapèze, nous pouvons utiliser une prédiction d'Euler progressive et remplacer le $u_{n+1}$ dans le terme $\\varphi(t_{n+1},u_{n+1})$ par $\\tilde u_{n+1}=u_n+h \\varphi(t_{n},u_{n})$. \n",
"Nous avons construit ainsi un nouveau schéma appelé **schéma de Heun**. \n",
"Plus précisément, la méthode de Heun s'écrit\n",
"\n",
"$$\n",
"\\text{(H)}\\quad\n",
"\\begin{cases}\n",
"u_0=y(t_0)=y_0,\\\\\n",
"\\tilde u_{n+1}=u_n+h \\varphi(t_{n},u_{n}),\\\\\n",
"u_{n+1}=u_n+\\dfrac{h}{2} \\Big(\\varphi(t_{n},u_{n}) + \\varphi(t_{n+1},\\tilde u_{n+1})\\Big) & n=0,1,2,\\dots N-1\n",
"\\end{cases}\n",
"$$\n",
"
\n",
"Il s'agit à nouveau d'un schéma à 1 pas explicite."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"def heun(phi,tt,y0):\n",
" h = tt[1]-tt[0]\n",
" uu = [y0]\n",
" for i in range(len(tt)-1):\n",
" k1 = phi( tt[i], uu[i] )\n",
" k2 = phi( tt[i+1], uu[i] + h*k1 )\n",
" uu.append( uu[i] + 0.5*h*(k1+k2) )\n",
" return uu"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Remarques\n",
"\n",
"1. À première vue, il semble que les schémas d'**Euler progressif** et de **Heun** soient préférable aux schémas d'**Euler rétrograde** et de **Crank-Nicolson** puisque ces derniers ne sont pas explicites. Cependant, on verra que les méthodes d'Euler implicite et de Crank-Nicolson sont inconditionnellement A-stables. C'est aussi le cas de nombreuses autres méthodes implicites. Cette propriété rend les méthodes implicites attractives, bien qu'elles soient plus coûteuses que les méthodes explicites.\n",
"\n",
"2. Pour la mise en application d'un schéma il faut aussi prendre en compte l'influence des erreurs d'arrondi. En effet, afin de minimiser l'erreur globale théorique, on pourrait être tenté d'appliquer une méthode avec un pas très petit, par exemple de l'ordre de $10^{-16}$, mais ce faisant, outre que le temps de calcul deviendrait irréaliste, très rapidement les erreurs d'arrondi feraient diverger la solution approchée. En pratique il faut prendre $h$ assez petit pour que la méthode converge, mais pas trop petit non plus pour que les erreurs d'arrondi ne donnent pas lieu à des résultats incohérents et pour que les calculs puissent être effectués en un temps raisonnable."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Annexes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Schémas basés sur la formule de quadrature de Simpson \n",
"\n",
"#### Schéma de Simpson implicite \n",
"Si on remplace une fonction $f$ par la parabole segment qui passe par $(a,f(a))$, $(b,f(b))$ et $\\left(\\tfrac{a+b}{2},f\\left(\\tfrac{a+b}{2}\\right)\\right)$ \n",
"(polynôme qui interpole $f$ en les points $(a,f(a))$, $(b,f(b))$ et $\\left(\\tfrac{a+b}{2},f\\left(\\tfrac{a+b}{2}\\right)\\right)$ et donc de degré $2$), on a\n",
"\\begin{align*}\n",
"\\tilde f(x)&=f(a)\\frac{(x-b)\\left(x-\\tfrac{a+b}{2}\\right)}{(a-b)\\left(a-\\tfrac{a+b}{2}\\right)}\n",
" +f\\left(\\tfrac{a+b}{2}\\right)\\frac{(x-a)(x-b)}{\\left(\\tfrac{a+b}{2}-a\\right)\\left(\\tfrac{a+b}{2}-b\\right)}\n",
" +f(b)\\frac{(x-a)\\left(x-\\tfrac{a+b}{2}\\right)}{(b-a)\\left(b-\\tfrac{a+b}{2}\\right)}\n",
"\\\\\n",
"\\int_a^b f(x)\\mathrm{d}x &\\approx \\int_a^b \\tilde f(x)\\mathrm{d}x=\\frac{h}{6} \\left( f(a)+4f\\left(\\tfrac{a+b}{2}\\right)+ f(b)\\right).\n",
"\\end{align*}\n",
"Cette formule est dite *formule de Simpson*. \n",
"\n",
"En utilisant cette formule pour approcher la fonction $t\\mapsto \\varphi(t,y(t))$ on a\n",
"$$\n",
"\\int_{t_n}^{t_{n+1}} \\varphi(t,y(t))\\mathrm{d}t\\approx \\frac{h}{6} \\left( \\varphi(t_n,y(t_n))+4\\varphi\\left(t_n+\\frac{h}{2},y\\left(t_n+\\frac{h}{2}\\right)\\right)+ \\varphi(t_{n+1},y(t_{n+1})) \\right)\n",
"$$\n",
"et on obtient \n",
"$$\\begin{cases}\n",
"u_0=y(t_0)=y_0,\\\\\n",
"u_{n+1}=u_n+\\frac{h}{6} \\left(\\varphi(t_n,u_n)+4\\varphi\\left(t_n+\\frac{h}{2},u_{n+1/2}\\right)+\\varphi(t_{n+1},u_{n+1})\\right)& n=0,1,2,\\dots N-1\n",
"\\end{cases}$$\n",
"où $u_{n+1/2}$ est une approximation de $y(t_n+h/2)$. \n",
"\n",
"Cependant, $u_{n+1/2}$ est une valeur inconnue, on doit alors en calculer une approximation $\\tilde u_{n+1/2}$.\n",
"Pour cela nous pouvons utiliser une prédiction d'Euler progressive sur un demi-pas: $\\tilde u_{n+1/2}=u_n+(h/2) \\varphi(t_{n},u_{n})$.\n",
"On peut ainsi remplacer $\\varphi(t_{n}+h/2,u_{n+1/2})$ par $\\varphi(t_{n}+h/2,\\tilde u_{n+1/2})$. \n",
"\n",
"Nous avons construit ainsi un nouveau schéma appelé **schéma de Simpson implicite** qui s'écrit\n",
"\n",
"$$\n",
"\\text{(SI)}\\quad\n",
"\\begin{cases}\n",
"u_0=y(t_0)=y_0,\\\\\n",
"\\tilde u_{n+1/2}=u_n+\\frac{h}{2}\\varphi(t_{n},u_{n}),\\\\\n",
"u_{n+1}=u_n+\\frac{h}{6} \\left(\\varphi(t_n,u_n)+4\\varphi\\left(t_n+\\frac{h}{2},\\tilde u_{n+1/2}\\right)+\\varphi(t_{n+1},u_{n+1})\\right)& n=0,1,2,\\dots N-1\n",
"\\end{cases} \n",
"$$\n",
"
\n",
"Il s'agit d'un schéma à 1 pas implicite car il ne permet pas d'expliciter $u_{n+1}$ en fonction de $u_n$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Schéma de Simpson explicite \n",
"Pour éviter le calcul implicite de $u_{n+1}$ dans le schéma de Simpson implicite, nous pouvons utiliser une prédiction d'Euler progressive et remplacer le $u_{n+1}$ dans le terme $\\varphi(t_{n+1},u_{n+1})$ par $\\check u_{n+1}=u_n+h \\varphi(t_{n},u_{n})$. \n",
"Nous avons construit ainsi un nouveau schéma qu'on appellera **schéma de Simpson explicite** et qui s'écrit\n",
"\n",
"$$\n",
"\\text{(SE)}\\quad\n",
"\\begin{cases}\n",
"u_0=y(t_0)=y_0,\\\\\n",
"\\tilde u_{n+1/2}=u_n+\\frac{h}{2}\\varphi(t_{n},u_{n}),\\\\\n",
"\\check u_{n+1}=u_n+h \\varphi(t_{n},u_{n}),\\\\\n",
"u_{n+1}=u_n+\\frac{h}{6} \\left(\\varphi(t_n,u_n)+4\\varphi\\left(t_n+\\frac{h}{2},\\tilde u_{n+1/2}\\right)+\\varphi(t_{n+1},\\check u_{n+1})\\right)& n=0,1,2,\\dots N-1\n",
"\\end{cases}\n",
"$$\n",
"
\n",
"Il s'agit à nouveau d'un schéma à 1 pas explicite."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Schéma de Simpson explicite-variante\n",
"Notons qu'on aurait pu remplacer le $u_{n+1}$ dans le terme $\\varphi(t_{n+1},u_{n+1})$ par une approximation utilisant $\\tilde u_{n+1/2}$ comme par exemple une prédiction d'Euler progressive à partir de $t_n+h/2$, ce qui donne $\\hat u_{n+1}=\\tilde u_{n+1/2}+\\frac{h}{2}\\varphi(t_{n}+h/2,\\tilde u_{n+1/2})$: \t\n",
"\n",
"$$\n",
"\\text{(SEV)}\\quad\n",
"\\begin{cases}\n",
"u_0=y(t_0)=y_0,\\\\\n",
"\\tilde u_{n+1/2}=u_n+\\frac{h}{2}\\varphi(t_{n},u_{n}),\\\\\n",
"\\hat u_{n+1}=\\tilde u_{n+1/2}+\\frac{h}{2}\\varphi(t_{n}+\\frac{h}{2},\\tilde u_{n+1/2}),\\\\\n",
"u_{n+1}=u_n+\\frac{h}{6} \\left(\\varphi(t_n,u_n)+4\\varphi\\left(t_n+\\frac{h}{2},\\tilde u_{n+1/2}\\right)+\\varphi(t_{n+1},\\hat u_{n+1})\\right)& n=0,1,2,\\dots N-1\n",
"\\end{cases}\n",
"$$\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Schémas de Simpson multipas\n",
"Au lieu de travailler sur $[t_n,t_{n+1}]$ on peut répéter la même construction sur $[t_{n-1},t_{n+1}]$ (et donc intégrer entre $t_{n-1}$ et $t_{n+1}$).\n",
"\n",
"Le SI devient alors\n",
"\n",
"$$\n",
"\\text{(SMI)}\\quad\n",
"\\begin{cases}\n",
"u_0=y(t_0)=y_0,\\\\\n",
"u_{1}=u_0+h\\varphi(t_{0},u_{0}),\\\\\n",
"u_{n+1}=u_{n-1}+\\frac{h}{3} \\left(\\varphi(t_{n-1},u_{n-1})+4\\varphi(t_n,u_{n})+\\varphi(t_{n+1},u_{n+1})\\right)& n=1,2,3,\\dots N-1\n",
"\\end{cases} \n",
"$$\n",
"
\n",
"\n",
"On peut remplacer le $u_{n+1}$ dans le terme $\\varphi(t_{n+1},u_{n+1})$ par une approximation utilisant une prédiction d'Euler progressive à partir de $t_n$, ce qui donne $\\tilde u_{n+1}=u_{n}+h\\varphi(t_{n},u_{n})$: \t\n",
"\n",
"$$\n",
"\\text{(SME)}\\quad\n",
"\\begin{cases}\n",
"u_0=y(t_0)=y_0,\\\\\n",
"u_{1}=u_0+h\\varphi(t_{0},u_{0}),\\\\\n",
"\\tilde u_{n+1}=u_{n}+h\\varphi(t_{n},u_{n}),\\\\\n",
"u_{n+1}=u_{n-1}+\\frac{h}{3} \\left(\\varphi(t_{n-1},u_{n-1})+4\\varphi(t_n,u_{n})+\\varphi(t_{n+1},\\tilde u_{n+1})\\right)& n=1,2,3,\\dots N-1\n",
"\\end{cases} \n",
"$$\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Calcul exact des polynomes interpolants et leur intégration\n",
"Pour vérifier nos calculs d'interpolation puis intégration, nous pouvons utiliser le package de calcul formel `SymPy`.\n",
"\n",
"On notera\n",
"- `t_n` : $t_n$\n",
"- `t_nm1` : $t_{n-1}=t_n-h$\n",
"- `t_np1` : $t_{n+1}=t_n+h$\n",
"- `t_np12` : $t_{n+1/2}=t_n+\\frac{h}{2}$\n",
"- `phi_n` : $\\varphi(t_n,y(t_n))$\n",
"- `phi_nm1` : $\\varphi(t_{n-1},y(t_{n-1}))$\n",
"- `phi_np1` : $\\varphi(t_{n+1},y(t_{n+1}))$\n",
"- `phi_np12` : $\\varphi(t_{n+1/2},y(t_{n+1/2}))$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"EE:\n"
]
},
{
"data": {
"text/latex": [
"$\\displaystyle t\\mapsto \\tilde f(t)=\\phi_{n}$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle \\int_{t_n}^{t_{n+1}}\\tilde f(t) \\mathrm{d}t=\\phi_{n} h$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"\n",
"EI:\n"
]
},
{
"data": {
"text/latex": [
"$\\displaystyle t\\mapsto \\tilde f(t)=\\phi_{n+1}$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle \\int_{t_n}^{t_{n+1}}\\tilde f(t) \\mathrm{d}t=\\phi_{n+1} h$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"\n",
"EM-init:\n"
]
},
{
"data": {
"text/latex": [
"$\\displaystyle t\\mapsto \\tilde f(t)=\\phi_{n+\\frac{1}{2}}$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle \\int_{t_n}^{t_{n+1}}\\tilde f(t) \\mathrm{d}t=\\phi_{n+\\frac{1}{2}} h$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"\n",
"CN:\n"
]
},
{
"data": {
"text/latex": [
"$\\displaystyle t\\mapsto \\tilde f(t)=\\frac{\\phi_{n+1} t}{h} - \\frac{\\phi_{n+1} t_{n}}{h} + \\phi_{n} - \\frac{\\phi_{n} t}{h} + \\frac{\\phi_{n} t_{n}}{h}$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle \\int_{t_n}^{t_{n+1}}\\tilde f(t) \\mathrm{d}t=\\frac{h \\left(\\phi_{n+1} + \\phi_{n}\\right)}{2}$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"\n",
"SI:\n"
]
},
{
"data": {
"text/latex": [
"$\\displaystyle t\\mapsto \\tilde f(t)=- \\frac{\\phi_{n+1} t}{h} + \\frac{\\phi_{n+1} t_{n}}{h} + \\frac{2 \\phi_{n+1} t^{2}}{h^{2}} - \\frac{4 \\phi_{n+1} t t_{n}}{h^{2}} + \\frac{2 \\phi_{n+1} t_{n}^{2}}{h^{2}} + \\frac{4 \\phi_{n+\\frac{1}{2}} t}{h} - \\frac{4 \\phi_{n+\\frac{1}{2}} t_{n}}{h} - \\frac{4 \\phi_{n+\\frac{1}{2}} t^{2}}{h^{2}} + \\frac{8 \\phi_{n+\\frac{1}{2}} t t_{n}}{h^{2}} - \\frac{4 \\phi_{n+\\frac{1}{2}} t_{n}^{2}}{h^{2}} + \\phi_{n} - \\frac{3 \\phi_{n} t}{h} + \\frac{3 \\phi_{n} t_{n}}{h} + \\frac{2 \\phi_{n} t^{2}}{h^{2}} - \\frac{4 \\phi_{n} t t_{n}}{h^{2}} + \\frac{2 \\phi_{n} t_{n}^{2}}{h^{2}}$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle \\int_{t_n}^{t_{n+1}}\\tilde f(t) \\mathrm{d}t=\\frac{h \\left(\\phi_{n+1} + 4 \\phi_{n+\\frac{1}{2}} + \\phi_{n}\\right)}{6}$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"\n",
"SMI:\n"
]
},
{
"data": {
"text/latex": [
"$\\displaystyle t\\mapsto \\tilde f(t)=\\frac{\\phi_{n+1} t}{2 h} - \\frac{\\phi_{n+1} t_{n}}{2 h} + \\frac{\\phi_{n+1} t^{2}}{2 h^{2}} - \\frac{\\phi_{n+1} t t_{n}}{h^{2}} + \\frac{\\phi_{n+1} t_{n}^{2}}{2 h^{2}} - \\frac{\\phi_{n-1} t}{2 h} + \\frac{\\phi_{n-1} t_{n}}{2 h} + \\frac{\\phi_{n-1} t^{2}}{2 h^{2}} - \\frac{\\phi_{n-1} t t_{n}}{h^{2}} + \\frac{\\phi_{n-1} t_{n}^{2}}{2 h^{2}} + \\phi_{n} - \\frac{\\phi_{n} t^{2}}{h^{2}} + \\frac{2 \\phi_{n} t t_{n}}{h^{2}} - \\frac{\\phi_{n} t_{n}^{2}}{h^{2}}$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle \\int_{t_{n-1}}^{t_{n+1}}\\tilde f(t) \\mathrm{d}t=\\frac{h \\left(\\phi_{n+1} + \\phi_{n-1} + 4 \\phi_{n}\\right)}{3}$"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"\n"
]
}
],
"source": [
"import sympy as sym\n",
"from IPython.display import display, Math\n",
"sym.init_printing()\n",
"\n",
"phi_n=sym.Symbol('\\phi_{n}')\n",
"phi_np12=sym.Symbol(r'\\phi_{n+\\frac{1}{2}}')\n",
"phi_np1=sym.Symbol('\\phi_{n+1}')\n",
"phi_nm1=sym.Symbol('\\phi_{n-1}')\n",
"sym.var('h,t,t_n')\n",
"\n",
"t_np1 = t_n+h\n",
"t_np12 = t_n+h/2\n",
"t_nm1 = t_n-h\n",
"\n",
"print(\"EE:\")\n",
"p = sym.interpolate([(t_n,phi_n)], t)\n",
"display(Math(r't\\mapsto \\tilde f(t)='+sym.latex(p)))\n",
"EE = sym.integrate(p,(t,t_n,t_np1)).simplify()\n",
"display(Math(r'\\int_{t_n}^{t_{n+1}}\\tilde f(t) \\mathrm{d}t='+sym.latex(EE)))\n",
"print(\"\\n\")\n",
"\n",
"print(\"EI:\")\n",
"p = sym.interpolate([(t_np1,phi_np1)], t)\n",
"display(Math(r't\\mapsto \\tilde f(t)='+sym.latex(p)))\n",
"EI = sym.integrate(p,(t,t_n,t_np1)).simplify()\n",
"display(Math(r'\\int_{t_n}^{t_{n+1}}\\tilde f(t) \\mathrm{d}t='+sym.latex(EI)))\n",
"print(\"\\n\")\n",
"\n",
"print(\"EM-init:\")\n",
"p = sym.interpolate([(t_np12,phi_np12)], t)\n",
"display(Math(r't\\mapsto \\tilde f(t)='+sym.latex(p)))\n",
"PM = sym.integrate(p,(t,t_n,t_np1)).simplify()\n",
"display(Math(r'\\int_{t_n}^{t_{n+1}}\\tilde f(t) \\mathrm{d}t='+sym.latex(PM)))\n",
"print(\"\\n\")\n",
"\n",
"print(\"CN:\")\n",
"p = sym.interpolate([(t_np1,phi_np1),(t_n,phi_n)], t)\n",
"display(Math(r't\\mapsto \\tilde f(t)='+sym.latex(p)))\n",
"CN = sym.integrate(p,(t,t_n,t_np1)).simplify()\n",
"display(Math(r'\\int_{t_n}^{t_{n+1}}\\tilde f(t) \\mathrm{d}t='+sym.latex(CN)))\n",
"print(\"\\n\")\n",
"\n",
"print(\"SI:\")\n",
"p = sym.interpolate([(t_np1,phi_np1),(t_np12,phi_np12),(t_n,phi_n)], t)\n",
"display(Math(r't\\mapsto \\tilde f(t)='+sym.latex(p)))\n",
"SI = sym.integrate(p,(t,t_n,t_np1)).simplify()\n",
"display(Math(r'\\int_{t_n}^{t_{n+1}}\\tilde f(t) \\mathrm{d}t='+sym.latex(SI)))\n",
"print(\"\\n\")\n",
"\n",
"print(\"SMI:\")\n",
"p = sym.interpolate([(t_np1,phi_np1),(t_n,phi_n),(t_nm1,phi_nm1)], t)\n",
"display(Math(r't\\mapsto \\tilde f(t)='+sym.latex(p)))\n",
"SI = sym.integrate(p,(t,t_nm1,t_np1)).simplify()\n",
"display(Math(r'\\int_{t_{n-1}}^{t_{n+1}}\\tilde f(t) \\mathrm{d}t='+sym.latex(SI)))\n",
"print(\"\\n\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"hide_input": false,
"kernelspec": {
"display_name": "Python 3",
"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.6"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": true,
"title_cell": "",
"title_sidebar": "Contents",
"toc_cell": true,
"toc_position": {
"height": "calc(100% - 180px)",
"left": "10px",
"top": "150px",
"width": "373.594px"
},
"toc_section_display": true,
"toc_window_display": true
}
},
"nbformat": 4,
"nbformat_minor": 4
}