from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())
import sys #only needed to determine Python version number
print('Python version ' + sys.version)
%reset -f
%matplotlib inline
from matplotlib.pylab import *
from scipy.optimize import fsolve
Considérons le problème de Cauchy
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).
Pour chaque nœud $t_n$, on cherche la valeur inconnue $u_n$ qui approche la valeur exacte $y_n\equiv y(t_n)$.
Méthodes explicites et méthodes implicites
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).
Une méthode est dite implicite si $u_{n+1}$ n'est défini que par une relation implicite faisant intervenir la fonction $\varphi$.
Méthodes à un pas et méthodes multi-pas
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.
Autrement, on dit que le schéma est une méthode multi-pas (ou à pas multiples).
Si nous intégrons l'EDO $y'(t)=\varphi(t,y(t))$ entre $t_n$ et $t_{n+1}$ nous obtenons $$ \int_{t_n}^{t_{n+1}}y'(t)\mathrm{d}t=\int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t $$ c'est-à-dire $$ y_{n+1}-y_n=\int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t. $$ On peut construire différents schémas selon la formule d'approximation utilisée pour approcher le membre de droite. Cette solution approchée sera obtenue en construisant une suite récurrente comme suit:
Si on remplace une fonction $f$ par une constante égale à la valeur de $f$ à borne gauche de l'intervalle $[a;b]$
(polynôme qui interpole $f$ au point $(a,f(a))$ et donc de degré $0$),
on a
\begin{align*}
\tilde f(x)&=f(a)\\
\int_a^b f(x)\mathrm{d}x &\approx \int_a^b \tilde f(x)\mathrm{d}x=(b-a)f(a).
\end{align*}
Cette formule est dite formule de quadrature du rectangle à gauche.
En utilisant cette formule pour approcher la fonction $t\mapsto \varphi(t,y(t))$ on a $$ \int_{t_n}^{t_{n+1}} \varphi(t,y(t)) \mathrm{d}t \approx h \varphi(t_n,y(t_n)) $$ et on reconnait le schéma d'Euler progressif
def EE(phi,tt,y0):
h = tt[1]-tt[0]
uu = [y0]
for i in range(len(tt)-1):
uu.append( uu[i]+h*phi(tt[i],uu[i]) )
return uu
# def euler_progressif(phi,tt,y0):
# h = tt[1]-tt[0]
# N = len(tt)-1 # tt contient N+1 points numerotés de 0 à N
# uu = [y0 for i in range(N+1)]
# for i in range(N): # i = 0,1,...,N-1
# uu[i+1] = uu[i]+h*phi(tt[i],uu[i])
# return uu
Si on remplace une fonction $f$ par une constante égale à la valeur de $f$ à borne droite de l'intervalle $[a;b]$
(polynôme qui interpole $f$ au point $(b,f(b))$ et donc de degré $0$), on a
\begin{align*}
\tilde f(x)&=f(b)\\
\int_a^b f(x)\mathrm{d}x &\approx \int_a^b \tilde f(x)\mathrm{d}x=(b-a)f(b).
\end{align*}
Cette formule est dite formule de quadrature du rectangle à droite.
En utilisant cette formule pour approcher la fonction $t\mapsto \varphi(t,y(t))$ on a $$ \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})) $$ et on obtient le schéma d'Euler rétrograde
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).$$
def EI(phi,tt,y0):
h = tt[1]-tt[0]
uu = [y0]
for i in range(len(tt)-1):
uu.append( fsolve(lambda x: -x+uu[i]+h*phi(tt[i+1],x), uu[i]) )
return uu
Si on remplace une fonction $f$ par une constante égale à la valeur de $f$ au milieu de l'intervalle $[a;b]$
(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
\begin{align*}
\tilde f(x)&=f\left(\tfrac{a+b}{2}\right)\\
\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).
\end{align*}
Cette formule est dite formule de quadrature du rectangle ou du point milieu.
En utilisant cette formule pour approcher la fonction $t\mapsto \varphi(t,y(t))$ on a $$ \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) $$ et on obtient $$\begin{cases} u_0=y(t_0)=y_0,\\ 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 \end{cases}$$ où $u_{n+1/2}$ est une approximation de $y(t_n+h/2)$.
Cependant, $u_{n+1/2}$ est une valeur inconnue, on doit alors en calculer une approximation $\tilde u_{n+1/2}$. Pour cela nous pouvons utiliser deux stratégies:
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})$.
On peut ainsi remplacer $\varphi(t_{n}+h/2,u_{n+1/2})$ par $\varphi(t_{n}+h/2,\tilde u_{n+1/2})$.
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$.
Utiliser une moyenne: $\tilde u_{n+1/2}=\dfrac{u_n+u_{n+1}}{2}$.
On peut ainsi remplacer $\varphi(t_{n}+h/2,u_{n+1/2})$ par $\varphi(t_{n}+h/2,\tilde u_{n+1/2})$.
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}).$$
def EM(phi,tt,y0):
h = tt[1]-tt[0]
uu = [y0]
for i in range(len(tt)-1):
u_tilde = uu[i]+0.5*h*phi(tt[i],uu[i])
uu.append( uu[i]+h*phi( tt[i]+h/2,u_tilde ) )
return uu
def RK1_M(phi,tt,y0):
h = tt[1]-tt[0]
uu = [y0]
for i in range(len(tt)-1):
uu.append( fsolve(lambda x: -x+uu[i]+h*phi( (tt[i]+tt[i+1])/2,(uu[i]+x)/2 ), uu[i]) )
return uu
Si on remplace une fonction $f$ par le segment qui relie $(a,f(a))$ à $(b,f(b))$
(polynôme qui interpole $f$ en les points $(a,f(a))$ et $(b,f(b))$ et donc de degré $1$), on a
\begin{align*}
\tilde f(x)&=\dfrac{f(b)-f(a)}{b-a}(x-a)+f(a)\\
\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).
\end{align*}
Cette formule est dite formule de quadrature du trapèze.
En utilisant cette formule pour approcher la fonction $t\mapsto \varphi(t,y(t))$ on a $$ \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) $$ et on obtient le schéma du trapèze ou de Crank-Nicolson
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. 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).$$
def CN(phi,tt,y0):
h = tt[1]-tt[0]
uu = [y0]
for i in range(len(tt)-1):
uu.append( fsolve(lambda x: -x+uu[i]+0.5*h*( phi(tt[i],uu[i])+phi(tt[i+1],x) ), uu[i]) )
return uu
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})$. Nous avons construit ainsi un nouveau schéma appelé schéma de Heun. Plus précisément, la méthode de Heun s'écrit
def heun(phi,tt,y0):
h = tt[1]-tt[0]
uu = [y0]
for i in range(len(tt)-1):
k1 = phi( tt[i], uu[i] )
k2 = phi( tt[i+1], uu[i] + h*k1 )
uu.append( uu[i] + 0.5*h*(k1+k2) )
return uu
À 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.
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.
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)$
(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
\begin{align*}
\tilde f(x)&=f(a)\frac{(x-b)\left(x-\tfrac{a+b}{2}\right)}{(a-b)\left(a-\tfrac{a+b}{2}\right)}
+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)}
+f(b)\frac{(x-a)\left(x-\tfrac{a+b}{2}\right)}{(b-a)\left(b-\tfrac{a+b}{2}\right)}
\\
\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).
\end{align*}
Cette formule est dite formule de Simpson.
En utilisant cette formule pour approcher la fonction $t\mapsto \varphi(t,y(t))$ on a
$$
\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)
$$
et on obtient
$$\begin{cases}
u_0=y(t_0)=y_0,\\
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
\end{cases}$$
où $u_{n+1/2}$ est une approximation de $y(t_n+h/2)$.
Cependant, $u_{n+1/2}$ est une valeur inconnue, on doit alors en calculer une approximation $\tilde u_{n+1/2}$. 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})$. On peut ainsi remplacer $\varphi(t_{n}+h/2,u_{n+1/2})$ par $\varphi(t_{n}+h/2,\tilde u_{n+1/2})$.
Nous avons construit ainsi un nouveau schéma appelé schéma de Simpson implicite qui s'écrit
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})$. Nous avons construit ainsi un nouveau schéma qu'on appellera schéma de Simpson explicite et qui s'écrit
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})$:
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}$).
Le SI devient alors
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})$:
Pour vérifier nos calculs d'interpolation puis intégration, nous pouvons utiliser le package de calcul formel SymPy.
On notera
t_n : $t_n$t_nm1 : $t_{n-1}=t_n-h$t_np1 : $t_{n+1}=t_n+h$t_np12 : $t_{n+1/2}=t_n+\frac{h}{2}$phi_n : $\varphi(t_n,y(t_n))$phi_nm1 : $\varphi(t_{n-1},y(t_{n-1}))$phi_np1 : $\varphi(t_{n+1},y(t_{n+1}))$phi_np12 : $\varphi(t_{n+1/2},y(t_{n+1/2}))$import sympy as sym
from IPython.display import display, Math
sym.init_printing()
phi_n=sym.Symbol('\phi_{n}')
phi_np12=sym.Symbol(r'\phi_{n+\frac{1}{2}}')
phi_np1=sym.Symbol('\phi_{n+1}')
phi_nm1=sym.Symbol('\phi_{n-1}')
sym.var('h,t,t_n')
t_np1 = t_n+h
t_np12 = t_n+h/2
t_nm1 = t_n-h
print("EE:")
p = sym.interpolate([(t_n,phi_n)], t)
display(Math(r't\mapsto \tilde f(t)='+sym.latex(p)))
EE = sym.integrate(p,(t,t_n,t_np1)).simplify()
display(Math(r'\int_{t_n}^{t_{n+1}}\tilde f(t) \mathrm{d}t='+sym.latex(EE)))
print("\n")
print("EI:")
p = sym.interpolate([(t_np1,phi_np1)], t)
display(Math(r't\mapsto \tilde f(t)='+sym.latex(p)))
EI = sym.integrate(p,(t,t_n,t_np1)).simplify()
display(Math(r'\int_{t_n}^{t_{n+1}}\tilde f(t) \mathrm{d}t='+sym.latex(EI)))
print("\n")
print("EM-init:")
p = sym.interpolate([(t_np12,phi_np12)], t)
display(Math(r't\mapsto \tilde f(t)='+sym.latex(p)))
PM = sym.integrate(p,(t,t_n,t_np1)).simplify()
display(Math(r'\int_{t_n}^{t_{n+1}}\tilde f(t) \mathrm{d}t='+sym.latex(PM)))
print("\n")
print("CN:")
p = sym.interpolate([(t_np1,phi_np1),(t_n,phi_n)], t)
display(Math(r't\mapsto \tilde f(t)='+sym.latex(p)))
CN = sym.integrate(p,(t,t_n,t_np1)).simplify()
display(Math(r'\int_{t_n}^{t_{n+1}}\tilde f(t) \mathrm{d}t='+sym.latex(CN)))
print("\n")
print("SI:")
p = sym.interpolate([(t_np1,phi_np1),(t_np12,phi_np12),(t_n,phi_n)], t)
display(Math(r't\mapsto \tilde f(t)='+sym.latex(p)))
SI = sym.integrate(p,(t,t_n,t_np1)).simplify()
display(Math(r'\int_{t_n}^{t_{n+1}}\tilde f(t) \mathrm{d}t='+sym.latex(SI)))
print("\n")
print("SMI:")
p = sym.interpolate([(t_np1,phi_np1),(t_n,phi_n),(t_nm1,phi_nm1)], t)
display(Math(r't\mapsto \tilde f(t)='+sym.latex(p)))
SI = sym.integrate(p,(t,t_nm1,t_np1)).simplify()
display(Math(r'\int_{t_{n-1}}^{t_{n+1}}\tilde f(t) \mathrm{d}t='+sym.latex(SI)))
print("\n")