1 Exercice 1 : implémentation d’un schéma pour un système
Le modèle S.Z.R., proposée par Munz, Hudea, Imad et Smith en 2009, est une variation du modèle S.I.R pour modéliser une apocalypse zombie.
À chaque instant on décide de diviser la population en trois catégories (qu’on appelle « compartiments » dans le langage de l’épidémiologie): - les individus « Susceptibles » ou « Sains » (S) : ceux qui n’ont jamais eu la maladie et peuvent la contracter ; - les individus « Zombies » (Z) ; - les individus « Rétablis » (R, comme « Recovered » en anglais) : les personnes décédées qui ont la faculté de revenir à la vie (terme \(\zeta R\)) mais aussi les individus débarassés de leur infection (suite à un coup de hache bien placé).
On travaillera avec des variables normalisées: leur valeur est comprise entre \(0\) et \(1\) et représente une fraction de la population totale.
Soit \(P(t)=S(t)+Z(t)+R(t)\). Montrer que \(P\) est un invariant.
Calculer la solution approchée obtenue par la méthode d’Euler Progressif avec \(61\) points. Afficher ensuite sur le même repère \(t\mapsto S(t)\), \(t\mapsto Z(t)\), \(t\mapsto R(t)\) et \(t\mapsto P(t)\).
Même exercice avec la méthode d’Euler Régressif.
Même exercice avec la fonction odeint du module scipy.
Testons ce modèle avec les coefficients
\(\alpha=1\) [humains sains moyennement violents]
\(\beta=0.5\) [zombies peu agressifs]
\(\zeta=0.25\) [résurrection des zombies non négligeable]
Q2
Calculer la solution approchée obtenue par la méthode d’Euler Progressif avec \(61\) points.
Afficher ensuite sur le même repère \(t\mapsto S(t)\), \(t\mapsto Z(t)\), \(t\mapsto R(t)\) et \(t\mapsto P(t)\).
On notera \(S_n\approx S(t_n)\), \(Z_n\approx Z(t_n)\) et \(R_n\approx R(t_n)\).
# S0, Z0, R0, h variables globalesdef euler_progressif(phi0,phi1,phi2,tt): SS = [S0] ZZ = [Z0] RR = [R0]for i inrange(len(tt)-1): SS.append(SS[i]+h*phi0(SS[i],ZZ[i],RR[i],tt[i])) ZZ.append(ZZ[i]+h*phi1(SS[i],ZZ[i],RR[i],tt[i])) RR.append(RR[i]+h*phi2(SS[i],ZZ[i],RR[i],tt[i]))return [SS,ZZ,RR][SS_ep, ZZ_ep, RR_ep] = euler_progressif(phi0,phi1,phi2,tt)figure(figsize=(10,5))PP_ep=[SS_ep[i]+ZZ_ep[i]+RR_ep[i] for i inrange(len(tt))]plot(tt,SS_ep,'--',tt,ZZ_ep,'.-',tt,RR_ep,':',tt,PP_ep)xlabel('t')legend(['S(t)','Z(t)','R(t)','P(t)'])title('Euler progressif') grid()axis([0,60,0,1.1]);
Q3
Calculer la solution approchée obtenue par la méthode d’Euler Régressif avec \(61\) points.
Afficher ensuite sur le même repère \(t\mapsto S(t)\), \(t\mapsto Z(t)\), \(t\mapsto R(t)\) et \(t\mapsto P(t)\).
from scipy.optimize import fsolve# S0, Z0, R0, h variables globalesdef euler_regressif(phi0,phi1,phi2,tt): SS = [S0] ZZ = [Z0] RR = [R0]for i inrange(len(tt)-1): sys =lambda z : [-z[0]+SS[i]+h*phi0(z[0],z[1],z[2],tt[i+1]) , -z[1]+ZZ[i]+h*phi1(z[0],z[1],z[2],tt[i+1]) , -z[2]+RR[i]+h*phi2(z[0],z[1],z[2],tt[i+1]) ] Stemp,Ztemp,Rtemp = fsolve( sys , (SS[i],ZZ[i],RR[i]) ) SS.append(Stemp) ZZ.append(Ztemp) RR.append(Rtemp)return [SS,ZZ,RR][SS_er, ZZ_er, RR_er] = euler_regressif(phi0,phi1,phi2,tt)figure(figsize=(10,5))PP_er=[SS_er[i]+ZZ_er[i]+RR_er[i] for i inrange(len(tt))]plot(tt,SS_er,'--',tt,ZZ_er,'.-',tt,RR_er,':',tt,PP_er)xlabel('t')legend(['S(t)','Z(t)','R(t)','P(t)'])title('Euler regressif') grid()axis([tt[0],tt[-1],0,1.1]);
Q4
Calculer la solution approchée obtenue par la fonction odeint du module scipy.
Afficher sur le même repère \(t\mapsto S(t)\), \(t\mapsto Z(t)\), \(t\mapsto R(t)\) et \(t\mapsto P(t)\).
Code
from scipy.integrate import odeintpphi =lambda yy,t : [ phi0(yy[0],yy[1],yy[2],t), phi1(yy[0],yy[1],yy[2],t), phi2(yy[0],yy[1],yy[2],t) ]CI = [S0,Z0,R0]sol = odeint(pphi,CI,tt)SS,ZZ,RR = sol[:,0],sol[:,1],sol[:,2]figure(figsize=(10,5))PP = [SS[i]+ZZ[i]+RR[i] for i inrange(len(tt))]plot(tt,SS,'--',tt,ZZ,'.-',tt,RR,':',tt,PP)xlabel('t')legend(['S(t)','Z(t)','R(t)','P(t)'])title('Odeint') grid()axis([0,60,0,1.1]);
2 Exercice 2 : étude d’un schéma multipas
On notera \(\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)\). Soit la méthode multipas \(
u_{n+1}
=2(1-\gamma) u_{n}
+\left(2\gamma-1 \right)u_{n-1}
+h(1-\gamma) \varphi_{n+1}
+h\left(3\gamma-1\right)\varphi_{n}.
\)
Q5
Pour quelles valeurs de $$ est-elle consistante?
\(p=1\): c’est une méthode à \(q=p+1=2\) : - $a_0=2(1-) $ et \(a_1=2\gamma-1\) - $b_0=3 $ et \(b_1=0\) - $b_{-1}=1-$
La méthode est implicite si \(\gamma\neq1\).
La méthode est consistante pour tout \(\gamma\) car \(
\begin{cases}
\displaystyle\sum_{j=0}^p a_j=a_0+a_1=2(1-\gamma) +(2\gamma-1)=1,
\\
\displaystyle-\sum_{j=0}^p ja_j+\sum_{j=-1}^p b_j=-\big(0a_0+1a_1\big)+\big(b_{-1}+b_0+b_1\big)=-(2\gamma-1)+(1-\gamma+3\gamma-1)=1
\end{cases}
\) Vérifions ces calculs ci-dessous:
Code
import sympy as symsym.init_printing()sym.var('gamma')p=1a0=2*(1-gamma)a1=2*gamma-1b0=3*gamma-1b1=0bm1=1-gammacond11=a0+a1cond12=-(0*a0+1*a1)+(bm1+b0+b1)display(cond11)display(cond12.factor())display(sym.solve([sym.Eq(cond11,1),sym.Eq(cond12,1)],gamma))
\(\displaystyle 1\)
\(\displaystyle 1\)
\(\displaystyle \left[ \right]\)
Q6
Pour quelles valeurs de \(\gamma\) est-elle convergente?
Le premier polynôme caractéristique est \(
\varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}=r^2-a_0r-a_1
{}=
r^2-2(1-\gamma) r-(2\gamma-1)
\) dont les racines sont \(
r_0=1,\quad r_1=1-2\gamma.
\)
\(\displaystyle \left(r - 1\right) \left(2 \gamma + r - 1\right)\)
\(\displaystyle \left[ 1, \ 1 - 2 \gamma\right]\)
La méthode est zéro-stable ssi \(
\begin{cases}
|r_j|\le1 \quad\text{pour tout }j=0,\dots,p
\\
\varrho'(r_j)\neq0 \text{ si } |r_j|=1,
\end{cases}
\) donc ssi \(0<\gamma\le1\).
Code
cond=sym.solve([abs(c)<=1for c in sol],gamma)display(cond)
La méthode est convergente ssi elle est consistante et zéro-stable donc ssi \(0<\gamma\le1\).
Q7
Quel ordre de convergente a-t-elle?
La méthode est d’ordre au moins 2 car \(
\sum_{j=0}^p (-j)^{2}a_j+2\sum_{j=-1}^p (-j)^{1}b_j
{}=\big(0^2a_0+1^2a_1\big)+2\big(-(-1)b_{-1}+0b_0+1b_1\big)
{}=\big(2\gamma-1\big)+2\big(1-\gamma\big)=1
\)
La méthode est d’ordre 3 ssi \(\gamma=\frac{2}{5}\) car \(
\sum_{j=0}^p (-j)^{3}a_j+3\sum_{j=-1}^p (-j)^{2}b_j=-a_1+3\big(b_{-1}+b_1\big)=-(2\gamma-1)+3\big(1-\gamma\big)=4-5\gamma=1
\iff \gamma=\frac{3}{5}
\)
On conclut que - le schéma est convergent au moins d’ordre 2 pour tout \(\gamma\in]0,1]\) - le schéma est convergent à l’ordre 3 pour \(\gamma=\frac{3}{5}\) (et on peut montrer qu’il n’est pas d’ordre 4)
Q8
Pour deux valeurs différentes de \(\gamma\) pour lesquelles le schéma est convergent (dont une qui donne l’ordre de convergence maximale), vérifier empiriquement l’ordre de convergence sur le problème de Cauchy \(
\begin{cases}
y'(t) = y(t)(1-y(t)), &\forall t \in I=[0,1],\\
y(0) = 2,
\end{cases}
\) après avoir calculé la solution exacte.
On définit la solution exacte (calculée en utilisant le module sympy) pour estimer les erreurs: