Séance 1 : résolution numérique des équations différentielles ordinaires¶
%pylab inline
Exercice 1 (hypothèses du théorème de Cauchy-Lipchitz et résolution approchée).¶
On considère l'équation différentielle $$ x'(t)=\sqrt{x(t)},\ t>0, $$ associée à la condition initiale $x(0)=0$.
Question : Expliquer en quoi ce problème ne relève pas du théorème de Cauchy-Lipchitz.
Réponse : La fonction racine carrée n'est pas lipschitzienne sur tout intervalle contenant l'origine et ne satisfait par conséquent pas les hypothèses du théorème de Cauchy-Lipshitz.
Question : Écrire un code de résolution numérique de ce problème de Cauchy sur l'intervalle $[0,10]$ utilisant la méthode d'Euler explicite sur une grille uniforme.
Pour cela, écrire une fonction euler_explicite dont les arguments d'entrée sont la fonction $f$ définissant le membre de droite de l'équation différentielle, les valeurs initiale et finale de la variable de temps $t$, la valeur de la donnée initiale du problème, le nombre $N$ de pas de la grille de discrétisation uniforme de l'intervalle de temps considéré, et dont les arguments de sorties sont le tableau contenant les valeurs de la variable de temps aux instants auxquels la solution approchée a été calculée et le tableau des valeurs de la solution approchée.
Réponse : On rappelle que la méthode d'Euler explicite pour la résolution numérique de l'équation différentielle ordinaire $x'(t)=f(t,x(t))$ sur l'intervalle $[t_0,t_0+T]$ est définie par la relation de récurrence, ou schéma, $$ x_{n+1}=x_n+h\,f(t_n,x_n),\ n=0,\dots,N-1, $$ avec $t_n=t_0+n\,h$ et $h=\frac{T}{N}$, où $N$ est le nombre de pas de la grille de discrétisation uniforme de $[t_0,t_0+T]$.
def euler_explicite(f,ti,tf,x0,N):
h=(tf-ti)/N
t,x=linspace(ti,tf,N+1),[]
x.append(x0)
for n in range(N):
x.append(x[n]+h*f(t[n],x[n]))
return t,x
Effectuer alors une résolution approchée du problème avec $N=100$ pas de discrétisation en exécutant le code suivant.
def f(t,x): return sqrt(x)
x0=array([0.])
t,x=euler_explicite(f,0.,10.,x0,100)
plot(t,x)
La solution calculée est identiquement nulle.
Question : Remplacer la valeur de la donnée initiale par l'évaluation numérique de $\lvert\left(\sqrt{2}\right)^2-2\rvert$ et effectuer la résolution approchée. Qu'observe-t-on ? Comment interpréter ce comportement de la méthode ?
x0=array([abs(sqrt(2.)**2-2.)])
t,x=euler_explicite(f,0.,10.,x0,100)
plot(t,x)
Réponse : En raison des erreurs d'arrondi ayant lieu lors des calculs en arithmétique en précision finie, la valeur initiale utilisée dans cette question peut être vue comme une petite perturbation de la valeur initiale précédemment employée. Bien que ce soit une perturbation minime (on pourra afficher la valeur numérique de la variable x0 pour s'en assurer), la solution approchée fournie par la méthode est très différente de celle précédemment calculée. La fonction $f$ considérée n'étant pas lipschitzienne, la stabilité de la méthode à l'égard des perturbations n'est en effet pas garantie.
Exercice 2 (modèle SIR).¶
On considère le système du modèle SIR de Kermack-McKendrick, décrivant l'évolution d'une population touchée par une maladie infectieuse, $$ \left\{ \begin{align*} S'(t)&=-r\,S(t)I(t),\\ I'(t)&=r\,S(t)I(t)-a\,I(t),\\ R'(t)&=a\,I(t), \end{align*} \right.\ t>0, $$ associé à la donnée initiale $$ S(0)=S_0,\ I(0)=I_0,\ R(0)=R_0. $$ On rappelle que les fonctions $S$, $I$ et $R$ représentent le nombre au cours du temps de personnes respectivement susceptibles de contracter la maladie, infectées et immunisées après avoir été malades, le réel $r$ est le taux d'infection et le réel $r$ est le taux de guérison. On utilisera les valeurs suivantes des paramètres du problème :
$$ S_0=762,\ I_0=1,\ R_0=0,\ r=0.00218,\ a=0.44036,\ T=14, $$
obtenues par un étalonnage du modèle effectué à partir des données d'une épidémie de grippe dans une école de garçons publiées dans la revue médicale britannique The Lancet le 4 mars 1978.
Question : Calculer une solution approchée du système sur l'intervalle de temps $[0,T]$ en utilisant la méthode d'Euler explicite (on se servira de la fonction écrite dans l'exercice précédent) et $N=10^5$ pas de discrétisation. Représenter l'évolution des valeurs des approximations obtenues des nombres $S(t)$, $I(t)$ et $R(t)$ au cours du temps.
r,a=0.00218,0.44036
def f(t,x) : return array([-r*x[0]*x[1],r*x[0]*x[1]-a*x[1],a*x[1]])
x0=array([762.,1.,0])
t,x=euler_explicite(f,0.,14.,x0,10**5)
[p1,p2,p3]=plot(t,x)
legend([p1,p2,p3],["S(t)","I(t)","R(t)"])
Le pic de l'épidémie est observé au septième jour.
Question : On note $N$ le nombre de pas de discrétisation employés. Vérifier théoriquement et numériquement que l'approximation obtenue vérifie, pour tout $0\leq n\leq N$, $$ S_n+I_n+R_n=S_0+I_0+R_0. $$
Réponse : En sommant les équations du système du modèle SIR, il vient $$ S'(t)+I'(t)+R'(t)=0, t>0, $$ dont on déduit que la somme $S(t)+I(t)+R(t)$ est constante au cours du temps et donc égale à $S_0+I_0+R_0$. On vérifie que c'est aussi le cas dans la simulation numérique du problème en superposant aux précédents graphes celui de l'évolution de la quantité $S(t)+I(t)+R(t)$ au cours du temps.
[p1,p2,p3,p4]=plot(t,x,t,sum(x,1))
legend([p1,p2,p3,p4],["S(t)","I(t)","R(t)","S(t)+I(t)+R(t)"])
Il y a bien conservation de la quantité par le schéma numérique.
Question : En prenant comme solution de référence celle obtenue avec $N=10^5$ pas de discrétisation, représenter (avec une échelle logarithmique) les erreurs au temps $T$ trouvées en utilisant successivement $N=10,10^2,\dots,10^4$ pas de discrétisation.
t,x=euler_explicite(f,0.,14.,x0,10**5)
xref,err,h=x[-1],[],[]
for i in range(4) :
N=10**(i+1)
h.append(14/N)
t,x=euler_explicite(f,0.,14.,x0,N)
err.append(norm(xref-x[-1]))
loglog(h,err,'+')
Retrouver alors l'ordre de convergence effectif de la méthode en utilisant la fonction polyfit de numpy.
Réponse : Pour déterminer l'ordre de la méthode, il faut revenir à la définition de cette notion. L'erreur globale d'une méthode d'ordre $p$, avec $p$ un entier, décroît avec le pas de discrétisation $h$ comme $O(h^p)$ et l'on a donc, pour $h$ suffisamment petit,
$$
\text{erreur}\simeq C\,h^p,
$$
avec $C$ une constante. En prenant le logarithme de cette relation, il vient
$$
\ln(\text{erreur})\simeq p\,\ln(h)+ln(C).
$$
Il suffit donc, à partir des valeurs de l'erreur approchée obtenue avec les différents choix de pas, déterminer la pente de la droite obtenue par une régression linéaire au moyen de la fonction polyfit.
slope=polyfit(log(h),log(err),1)
loglog(h,err,'+',h,exp(slope[1])*(h**(slope[0])))
print("L'ordre de convergence effectif vaut",slope[0])
La valeur de la pente calculée est proche de l'ordre théorique de la méthode.
Question : Écrire, sur le modèle de la fonction euler_explicite, une fonction rk4 pour la méthode de Runge--Kutta d'ordre quatre «classique». On rappelle que cette méthode a pour tableau de Butcher
$$
\begin{array}{c|cccc}
0\\
\frac{1}{2}&\frac{1}{2}\\
\frac{1}{2}&0&\frac{1}{2}\\
1&0&0&1\\
\hline
&\frac{1}{6}&\frac{1}{3}&\frac{1}{3}&\frac{1}{6}
\end{array}.
$$
def rk4(f,ti,tf,x0,N):
h=(tf-ti)/N
t,x=linspace(ti,tf,N+1),[]
x.append(x0)
for n in range(N):
k1=f(t[n],x[n])
k2=f(t[n]+0.5*h,x[n]+0.5*h*k1)
k3=f(t[n]+0.5*h,x[n]+0.5*h*k2)
k4=f(t[n]+h,x[n]+h*k3)
x.append(x[n]+h*(k1+2.*(k2+k3)+k4)/6.)
return t,x
Répondre à la question précédente en utilisant cette méthode.
t,x=rk4(f,0.,14.,x0,10**5)
xref,err,h=x[-1],[],[]
for i in range(4):
N=10**(i+1)
h.append(14/N)
t,x=rk4(f,0.,14.,x0,N)
err.append(norm(xref-x[-1]))
slope=polyfit(log(h),log(err),1)
loglog(h,err,'+',h,exp(slope[1])*(h**(slope[0])))
print("L'ordre de convergence effectif vaut",slope[0])
Là encore, la valeur de la pente calculée est proche de l'ordre théorique de la méthode.
Exercice 3 (équation raide et stabilité absolue).¶
On cherche à approcher numériquement la solution du problème de Cauchy $$ x'(t)=-50\,\left(x(t)-\cos(t)\right),\ t>0,\text{ et }x(0)=0, $$ sur l'intervalle $[0,2]$.
Question : Utiliser la méthode d'Euler explicite pour résoudre le problème avec un pas de discrétisation de longueur constante, successivement choisie
- strictement supérieure à 0,04,
- comprise entre 0,02 et 0,04,
- strictement inférieure à 0,02.
Qu'observe-t-on et quelle explication peut-on fournir ?
def f(t,x): return -50*(x-cos(t))
x0=array([0.])
figure()
t,x=euler_explicite(f,0,2,x0,int(2/0.05))
plot(t,x)
title("h=0.05")
figure()
t,x=euler_explicite(f,0,2,x0,int(2/0.03))
plot(t,x)
title("h=0.03")
figure()
t,x=euler_explicite(f,0,2,x0,int(2/0.01))
plot(t,x)
title("h=0.01")
Réponse : Pour un pas de longueur $h=0.005$, on observe un important phénomène d'instabilité numérique dans l'approximation calculée. Avec un pas de longueur $h=0.03$, le phénomène est moins marqué, mais clairement visible aux premiers instants de la simulation. Enfin, avec un pas de longueur $h=0.01$, l'approximation obtenue ne présente plus d'instabilité numérique.
On qualifie de raide tout système d'équations différentielles ordinaires pour lequel une méthode de résolution numérique explicite présente des problèmes d'instabilité numérique, à moins que la longueur du pas de discrétisation ne soit extrêmement petite (ou, pour être plus précis, beaucoup plus petite que ne l'exigerait en principe la précision demandée lorsqu'on utilise un mécanisme d'adaptation du pas). Les simulations effectuées indiquent ici qu'on est en présence d'une équation de ce type.
Question : Effectuer les calculs précédents avec la méthode d'Euler implicite. Pour cela, écrire une fonction euler_implicite dans laquelle la méthode de Newton-Raphson est employée pour résoudre le système d'équations non linéaires (on ajoutera dans les arguments d'entrée la fonction donnant la dérivée par rapport à $x$ de la fonction $f$ définissant le membre de droite de l'équation différentielle, ainsi que des paramètres d'arrêt pour la méthode de Newton-Raphson). Relancer les simulations et conclure.
def euler_implicite(f,df,ti,tf,x0,N,tol=1e-6,maxiter=50):
h=(tf-ti)/N
t,x=linspace(ti,tf,N+1),[]
x.append(x0)
neq=len(x0)
for n in range(N):
xnew=x[n]
for j in range(maxiter) :
rhs=(xnew-x[n])-h*f(t[n]+h,xnew)
jac=eye(neq)-h*df(t[n]+h,xnew)
inc=solve(jac,rhs)
xnew=xnew-inc
if norm(inc)<tol:
break
x.append(xnew)
return t,x
def df(t,x): return -50
x0=array([0.])
figure()
t,x=euler_implicite(f,df,0,2,x0,int(2/0.05))
plot(t,x)
title("h=0.05")
figure()
t,x=euler_implicite(f,df,0,2,x0,int(2/0.03))
plot(t,x)
title("h=0.03")
figure()
t,x=euler_implicite(f,df,0,2,x0,int(2/0.01))
plot(t,x)
title("h=0.01")
Réponse : La méthode d'Euler implicite possède un domaine de stabilité absolue non borné, qui contient l'ensemble des nombres complexes dont la partie réelle est strictement négative. On dit qu'elle est A-stable. Une méthode A-stable n'est pas sujette aux problèmes d'instabilité lorsqu'on l'utilise pour la résolution d'une équation raide, comme on peut le constater ici.