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(f'Python version: {sys.version}')
import numpy as np #only needed to determine Numpy version number
print(f"Numpy version: {np.__version__}")
from datetime import datetime
print('Last Updated On: ', datetime.now())
M62_CM1 Introduction à l'approximation numérique d'EDO¶
La résolution formelle d’équations différentielles s’avère très compliquée et limitée: la plupart des problèmes ne peuvent être qu’approchés.
On ne peut expliciter des solutions analytiques que pour des équations différentielles ordinaires très particulières. Par exemple :
- dans certains cas, on ne peut exprimer la solution que sous forme implicite. C'est le cas par exemple de l'EDO $y'(t)=\dfrac{y(t)-t}{y(t)+t}$ dont les solutions vérifient la relation implicite $$ \frac{1}{2}\ln(t^2+y^2(t))+\arctan\left( \frac{y(t)}{t} \right)=C, $$ où $C$ est une constante arbitraire.
- dans d'autres cas, on ne parvient même pas à représenter la solution sous forme implicite.
C'est le cas par exemple de l'EDO $y'(t)=e^{-t^2}$ dont les solutions ne peuvent pas s'écrire comme composition de fonctions élémentaires.
Lorsqu’on ne connait pas de solution exacte à un problème de Cauchy qui admet une et une seule solution, on essaye d’en avoir une bonne approximation par des méthodes numériques.
- 1 Position du problème
- 2 Construction élémentaire des méthodes d'Euler explicite et implicite
- 3 Implémentation des schémas d'Euler explicite et implicite
- 4 Convergence théorique des schémas d'Euler
- 5 Convergence empirique des schémas d'Euler
- 6 Annexe: implémentation d'un schéma "moyenne" des deux schémas précédents
Position du problème¶
L'idée la plus simple pour résoudre de manière approchée un problème de Cauchy est de discrétiser l'intervalle de temps avec un pas $h$ et d'approcher la dérivée temporelle sur chaque sous-intervalle de longieur $h$.
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)$.
- 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.
- L'ensemble de $N+1$ valeurs $\{y_0, y_1,\dots , y_{N} \}$ représente la solution exacte discrète.
- L'ensemble de $N+1$ valeurs $\{u_0 = y_0, u_1,\dots , u_{N} \}$ représente la solution numérique.
Construction élémentaire des méthodes d'Euler explicite et implicite¶
Une méthode classique, la méthode d'Euler explicite (ou progressive, de l'anglais forward), est obtenue en considérant l'équation différentielle en chaque nœud $t_n$ et en remplaçant la dérivée exacte $y'(t_n)$ par le taux d'accroissement $$ \varphi(t_n,y(t_n))=y'(t_n)=\lim_{h\to0}\frac{y(t_n+h)-y(t_n)}{h}\simeq\frac{y(t_{n+1})-y(t_n)}{h}. $$ Cela permet de construire une solution numérique par une suite récurrente:
De même, en utilisant le taux d'accroissement $$ \varphi(t_{n+1},y(t_{n+1}))=y'(t_{n+1})\simeq\frac{y(t_{n+1})-y(t_n)}{h} $$ pour approcher $y'(t_{n+1})$, on obtient la méthode d'Euler implicite (ou rétrograde, de l'anglais backward)
Ces deux méthodes sont dites à un pas: pour calculer la solution numérique $u_{n+1}$ au nœud $t_{n+1}$, on a seulement besoin des informations disponibles au nœud précédent $t_n$. Plus précisément,
- pour la méthode d'Euler progressive, $u_{n+1}$ ne dépend que de la valeur $u_n$ calculée précédemment,
- pour la méthode d'Euler rétrograde, $u_{n+1}$ dépend aussi "de lui-même" à travers la valeur de $\varphi(t_{n+1},u_{n+1})$.
C'est pour cette raison que la méthode d'Euler progressive est dite explicite tandis que la méthode d'Euler rétrograde est dite implicite.
Les méthodes implicites sont plus coûteuses que les méthodes explicites car, si la fonction $\varphi$ est non linéaire, un problème non linéaire doit être résolu à chaque temps $t_{n+1}$ pour calculer $u_{n+1}$.
Néanmoins, nous verrons que les méthodes implicites jouissent de meilleures propriétés de stabilité que les méthodes explicites.
Implémentation des schémas d'Euler explicite et implicite¶
Voyons un exemple complet.
On commence par importer
- la fonction
fsolvedu modulescipy.optimizepour résoudre les équations implicites présentes dans le schéma implicite. - le module
matplotlibet le modulenumpy:- soit en important séparemment
matplotlib.pyplotetnumpy(avec les alias classiques); - soit en important juste
matplotlib.pylab, ce qui importe à la foismatplotlibetnumpy. Dans ce cas on utilise la méthode du "paresseau" (i.e. on utilise*).
- soit en important séparemment
Rappels: numpy rédéfinit toutes les fonctions mathématiques du module math et ces fonctions sont vectorisées (e.g. on pourra écrire directement yy=sin(xx) avec xx une liste/tuple/array au lieu d'écrire yy=[sin(x) for x in xx]).
%reset -f
%matplotlib inline
from scipy.optimize import fsolve
#from matplotlib.pylab import *
# rcdefaults()
# rc('font', size=16)
import numpy as np
import matplotlib.pyplot as plt
# plt.rcdefaults()
plt.rcParams.update({'font.size': 16})
On initialise le problème de Cauchy
t0 = 0
tfinal = 1
y0 = 1
On définit l'équation différentielle : phi est une fonction python qui contient la fonction mathématique $\varphi\colon [t_0,T]\times\mathbb{R}\to\mathbb{R}$ définie par $\varphi(t, y)=2ty$. Bien noter que $\varphi$ dépend de deux variables: $t$ et $y$.
Dans l'EDO on a $\varphi(t,y(t))$: il s'agit d'une fonction d'une seule variable obtenue par composition de la fonction $\varphi$ avec la fonction $t\mapsto y(t)$.
phi = lambda t,y : 2*y*t
On introduit la discrétisation: les nœuds d'intégration $[t_0,t_1,\dots,t_{N}]$ sont contenus dans le vecteur tt.
On a $N+1$ points espacés de $h=\frac{t_N-t_0}{N}$.
N = 8
tt = np.linspace(t0,tfinal,N+1)
On écrit les schémas numériques : les valeurs $[u_0,u_1,\dots,u_{N}]$ pour chaque méthode sont contenues dans le vecteur uu.
Schéma d'Euler progressif : $$\begin{cases} u_0=y_0,\\ u_{n+1}=u_n+h\varphi(t_n,u_n)& n=0,1,2,\dots N-1 \end{cases}$$
def euler_progressif(phi,tt,y0):
h = tt[1]-tt[0]
N = len(tt)-1 # len(tt) = N+1 car tt contient N+1 points
uu = np.array([y0]) # array à partir d'une liste qui contient un seul element
for n in range(N): # range(N) = 0,1,...,N-1
uu = np.append( uu, uu[n]+h*phi(tt[n],uu[n]) )
return uu
# même fonction sans utiliser les array de numpy mais juste des listes
# def euler_progressif(phi,tt,y0):
# h = tt[1]-tt[0]
# uu = [y0]
# for n in range(len(tt)-1):
# uu.append( uu[-1]+h*phi(tt[-1],uu[-1]) )
# return uu # on renvoie une liste au lieu d'un np.array
Schéma d'Euler régressif : $$\begin{cases} u_0=y_0,\\ u_{n+1}=u_n+h\varphi(t_{n+1},u_{n+1})& n=0,1,2,\dots N-1 \end{cases}$$
Attention :
- $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 general non linéaire) $$x\mapsto -x+u_n+h\varphi(t_{n+1},x)$$
- la fonction
fsolvedu modulescipy.optimizerequiert deux paramètres : une fonction et un point de départ. Elle renvoie une liste avec la valeur du zéro approché (c'est une liste carfsolveresout un système, ici avec une seule équation).
def euler_regressif(phi,tt,y0):
h = tt[1]-tt[0]
uu = np.array([y0])
for n in range(len(tt)-1):
temp = fsolve( lambda x: -x+uu[n]+h*phi(tt[n+1],x) , uu[n] )
uu = np.append( uu, temp[0] )
return uu
On calcule les solutions approchées:
uu_ep = euler_progressif(phi,tt,y0)
uu_er = euler_regressif(phi,tt,y0)
Comme on la connait, on définit la solution exacte pour calculer les erreurs:
sol_exacte = lambda t : y0*np.exp(t**2)
yy = sol_exacte(tt) # = [sol_exacte(t) for t in tt]
On compare les graphes des solutions exacte (en bleu) et approchées (en rouge) et on affiche le maximum de l'erreur:
plt.figure(figsize=(18, 7))
plt.subplot(1,2,1)
plt.plot(tt,yy ,'b-o',label='Exacte')
plt.plot(tt,uu_ep,'r-D',label='Approchée')
erreur = np.abs(uu_ep-yy) # = [abs(uu_ep[i]-yy[i]) for i in range(N)] = [abs(u-y) for u,y in zip(uu_ep,yy)]
plt.title(f'Euler explicite\nmax(|erreur|) = {max(erreur):g}') # idem que np.norm(uu_ep-yy,inf)
plt.grid()
plt.legend();
plt.subplot(1,2,2)
plt.plot(tt,yy ,'b-o',label='Exacte')
plt.plot(tt,uu_er,'r-D',label='Approchée')
erreur = np.abs(uu_er-yy)
plt.title(f'Euler implicite\nmax(|erreur|) = {max(erreur):g}')
plt.grid()
plt.legend();
Convergence théorique des schémas d'Euler¶
Remarque: $N\to+\infty$ lorsque $h\to0$.
Soit $u_{n+1}^*$ la solution numérique au temps $t_{n+1}$ qu'on obtiendrait en insérant de la solution exacte dans le schéma (par exemple, pour la méthode d'Euler explicite on a $u_{n+1}^*\equiv y_{n} + h\varphi(t_{n}, y_{n})$). Pour vérifier qu'une méthode converge, on écrit l'erreur ainsi \begin{equation}\label{quart7.9} e_n \equiv y_n - u_n = (y_n - u_n^*) + (u_n^* - u_n). \end{equation} Si les deux termes $(y_n - u_n^*)$ et $(u_n^* - u_n)$ tendent vers zéro quand $h \to 0$ alors la méthode converge.
Remarque: la propriété de consistance est nécessaire pour avoir la convergence. En effet, si elle n'était pas consistante, la méthode engendrerait à chaque itération une erreur qui ne tendrait pas vers zéro avec $h$. L'accumulation de ces erreurs empêcherait l'erreur globale de tendre vers zéro quand $h \to 0$.
Preuve
On étudie séparément l'erreur de consistance et l'accumulation de ces erreurs. On en déduit ensuite la convergence.
- Terme $y_n - u_n^*$ (erreur de consistance).
Il représente l'erreur engendrée par une seule itération de la méthode d'Euler explicite.
En supposant que la dérivée seconde de $y$ existe et est continue, on écrit le développement de Taylor de $y$ au voisinage de $t_n$: $$ y(t_{n+1})=y(t_n)+hy'(t_n)+\frac{h^2}{2}y''(\eta_n) $$ où $\eta_n$ est un point de l'intervalle $]t_n;t_{n+1}[$. Donc il existe $\eta_n \in ]t_n, t_{n+1}[$ tel que $$ y_{n+1}-u_{n+1}^* =y_{n+1}-\Big(y_{n} + h\varphi(t_{n}, y_{n})\Big)=y_{n+1}-y_{n} - hy'(t_n)=\frac{h^2}{2}y''(\eta_n). $$ L'erreur de troncature de la méthode d'Euler explicite est donc de la forme $$ \tau(h)=M\frac{h}{2}, \qquad M\equiv\max_{t\in [t_0,T]}|y''(t)|. $$ On en déduit que $\lim_{h\to0} \tau (h) = 0$: la méthode est consistante.
Terme $u_{n+1}^* - u_{n+1}$ (accumulation des erreurs).
Il représente la propagation de $t_{n}$ à $t_{n+1}$ de l'erreur accumulée au temps précédent $t_{n}$. On a$$u_{n+1}^* - u_{n+1} = \Big(y_n + h \varphi(t_n, y_n)\Big)-\Big(u_n + h \varphi(t_n, y_n)\Big)=e_n+h\Big( \varphi(t_n,y_n)-\varphi(t_n,u_n) \Big). $$
Comme $\varphi$ est lipschitzienne par rapport à sa deuxième variable, on a $$ |u_{n+1}^{*} - u_{n+1}|\le (1 + hL)|e_{n}|. $$
- Convergence.
Comme $e_0 = 0$, les relations précédentes donnent \begin{align*} |e_n| &\le |y_n - u_n^*| + |u_n^* - u_n| \\ &\le h|\tau_n(h)| + (1+hL)|e_{n-1}| \\ &\le h|\tau_n(h)| + (1+hL)\left( h|\tau_{n-1}(h)| + (1+hL)|e_{n-2}| \right)| \\ &\le \left( 1+(1+hL)+\dots+(1+hL)^{n-1} \right)h\tau(h) \\ &=\left(\sum_{i=0}^{n-1}(1+hL)^i\right)h\tau(h) \\ &=\frac{(1+hL)^n-1}{hL}h\tau(h) \\ &\le\frac{(e^{hL})^n-1}{hL}h\tau(h) \qquad\text{car }(1+x)\le e^x \\ &=\frac{(e^{hL})^{(t_n-t_0)/h}-1}{L}\tau(h)\qquad\text{car } t_n-t_0=nh \\ &=\frac{e^{L(t_n-t_0)}-1}{L}\tau(h) \\ &=\frac{e^{L(t_n-t_0)}-1}{L}\frac{M}{2}h=\mathcal{O}(h) \end{align*} On peut conclure que la méthode d'Euler explicite est convergente d'ordre 1.
On remarque que l'ordre de cette méthode coïncide avec l'ordre de son erreur de troncature. On retrouve cette propriété dans de nombreuses méthodes de résolution numérique d'équations différentielles ordinaires.
Remarque: l'estimation de convergence est obtenue en supposant seulement $\varphi$ lipschitzienne. On peut établir une meilleure estimation si $\partial_y\varphi$ existe et est non négative pour tout $t \in [t_0;T]$ et tout $y\in\mathbb{R}$. En effet dans ce cas $$\begin{align*} u_n^* - u_n &= \left(y_{n-1} + h\varphi(t_{n-1}, y_{n-1})\right) - \left( u_{n-1} + h\varphi(t_{n-1},u_{n-1})\right) \\ &= e_{n-1}+h\left( \varphi(t_{n-1},y_{n-1})-\varphi(t_{n-1},u_{n-1}) \right) \\ &= e_{n-1}+h\left( e_{n-1}\partial_y\varphi(t_{n-1},\eta_n) \right) \\ &= \left( 1+h\partial_y\varphi(t_{n-1},\eta_n) \right)e_{n-1} \end{align*}$$ où $\eta_n$ appartient à l'intervalle dont les extrémités sont $y_{n-1}$ et $u_{n-1}$. Ainsi, si $$ 0<h<\frac{2}{\max\limits_{t\in[t_0,T]}\partial_y\varphi(t,y(t))} $$ alors $$ |u_n^* - u_n |\le |e_{n-1}|. $$ On en déduit $|e_n| \le |y_n - u_n^*| + |e_{n-1}| \le nh\tau(h) + |e_0|$ et donc $$ |e_n|\le M \frac{h}{2}(t_n-t_0). $$ La restriction sur le pas de discrétisation $h$ est une condition de stabilité, comme on le verra dans la suite.
Convergence empirique des schémas d'Euler¶
Considérons le même problème de Cauchy.
On se propose d'estimer l'ordre de convergence des méthodes d'Euler.
Pour chaque schéma, on calcule la solution approchée avec différentes valeurs de $h$ (ce qui correspond à différentes valeurs de $N$). On sauvegarde les valeurs de $h$ dans le vecteur H.
Pour chaque valeur de $h$, on calcule le maximum de la valeur absolue de l'erreur et on sauvegarde toutes ces erreurs dans le vecteur err_schema de sort que err_schema[k] contient $e_k=\max_{i=0,\dots,N_k}|y(t_i)-u_{i}|$ avec $N_k$.
H = []
err_ep = []
err_er = []
for N in range(8,7*5+9,5) :
tt = np.linspace(t0,tfinal,N+1)
H.append(tt[1]-tt[0])
yy = sol_exacte(tt)
uu_ep = euler_progressif(phi,tt,y0)
uu_er = euler_regressif(phi,tt,y0)
err_ep.append( max( abs(uu_ep-yy) ) )
err_er.append( max( abs(uu_er-yy) ) )
Pour afficher l'ordre de convergence on affiche les points (h[k],err_ep[k]) en échelle logarithmique: on représente $\ln(h)$ sur l'axe des abscisses et $\ln(\text{err})$ sur l'axe des ordonnées. En effet, si $\text{err}=Ch^p$ alors $\ln(\text{err})=\ln(C)+p\ln(h)$, autrement dit, en échelle logarithmique, $p$ représente la pente de la ligne droite $\ln(\text{err})$.
plt.figure(figsize=(10, 7))
plt.loglog(H,err_ep, 'r-o',label='Euler Explicite')
plt.loglog(H,err_er, 'g-+',label='Euler Implicite')
plt.xlabel('$h$')
plt.ylabel('erreur')
plt.legend(bbox_to_anchor=(1.04,1),loc='upper left')
plt.grid(True);
Pour estimer l'ordre de convergence on doit estimer la pente de la droite qui relie l'erreur au pas $k$ à l'erreur au pas $k+1$ en échelle logarithmique.
Pour estimer la pente globale de cette droite (par des moindres carrés) on peut utiliser la fonction polyfit (du module numpy que nous avons déjà importé avec matplotlib.pylab).
# ln(e) = a ln(h) + b
a_ep, b_ep = np.polyfit(np.log(H),np.log(err_ep), 1) # polyfit ( [liste des abscisses], [liste des ordonnées], degré du polynome)
print (f'Euler progressif: ordre {a_ep :1.2f}')
a_er, b_er = np.polyfit(np.log(H),np.log(err_er), 1)
print (f'Euler regressif: ordre {a_er:1.2f}')
On peut bien-sûr afficher la droite obtenue par régression linéaire en même temps que les points.
Soit on affiche les logarithmes des données et l'équation de la droite avec la commande plot:
plt.figure(figsize=(18, 7))
plt.subplot(1,2,1)
plt.plot(np.log(H),np.log(err_ep), 'ro',label='Euler Explicite')
plt.plot(np.log(H),[a_ep*np.log(h)+b_ep for h in H])
plt.xlabel('$\ln(h)$')
plt.ylabel('$\ln(e)$')
plt.legend(loc='upper left')
plt.grid(True);
plt.subplot(1,2,2)
plt.plot(np.log(H),np.log(err_er), 'ro',label='Euler Implicite')
plt.plot(np.log(H),[a_er*np.log(h)+b_er for h in H])
plt.xlabel('$\ln(h)$')
plt.ylabel('$\ln(e)$')
plt.legend(loc='upper left')
plt.grid(True);
Soit on affiche les données en échelle logarithmique et l'équation de l'erreur $Ch^p$ avec $C=\ln(b)$ et $p=a$ avec l'instruction loglog:
plt.figure(figsize=(18, 7))
plt.subplot(1,2,1)
plt.loglog(H,err_ep, 'ro',label='Euler Explicite')
plt.loglog(H,[np.exp(b_ep)*(h**a_ep) for h in H])
plt.xlabel('$h$')
plt.ylabel('$e$')
plt.legend(loc='upper left')
plt.grid(True);
plt.subplot(1,2,2)
plt.loglog(H,err_er, 'ro',label='Euler Implicite')
plt.loglog(H,[np.exp(b_er)*(h**a_er) for h in H])
plt.xlabel('$h$')
plt.ylabel('$e$')
plt.legend(loc='upper left')
plt.grid(True);
Annexe: implémentation d'un schéma "moyenne" des deux schémas précédents¶
On obtient la méthode de Crank-Nicolson (ou du trapèze)
%reset -f
%matplotlib inline
from scipy.optimize import fsolve
from matplotlib.pylab import *
#rcdefaults()
rc('font', size=14)
t0 = 0
tfinal = 1
y0 = 1
phi = lambda t,y : 2*y*t
N = 8
tt = linspace(t0,tfinal,N+1)
def CN(phi,tt,y0):
h = tt[1]-tt[0]
uu = [y0]
for i in range(len(tt)-1):
temp = fsolve( lambda x: -x+uu[i]+h*(phi(tt[i],uu[i])+phi(tt[i+1],x))/2 , uu[i] )
uu.append(temp[0])
return uu
uu_cn = CN(phi,tt,y0)
sol_exacte = lambda t : y0*exp(t**2)
yy = [sol_exacte(t) for t in tt]
figure(1, figsize=(18, 7))
subplot(1,2,1)
plot(tt,yy,'b-o',label='Exacte')
plot(tt,uu_cn,'r-D',label='Approchée')
erreur = [abs(uu_cn[i]-yy[i]) for i in range(N)]
title(f'CN\nmax(|erreur|) = {max(erreur):g}')
grid()
legend();
H = []
err_cn = []
for k in range(7):
N += 5
tt = linspace(t0,tfinal,N+1)
h = tt[1]-tt[0]
yy = [sol_exacte(t) for t in tt]
uu_cn = CN(phi,tt,y0)
H.append(h)
err_cn.append( max([abs(u-y) for u,y in zip(uu_cn,yy) ]) )
# ln(e) = a ln(h) + b
a_cn, b_cn = polyfit(log(H),log(err_cn), 1)
print (f'Crank Nicolson: ordre {a_cn :1.2f}')