2021 Contrôle Continu – M64 – Data Science & Scientific Computing 6 : Schémas numériques pour les problèmes de Cauchy

2021 Contrôle Continu

Auteur·rice

Gloria Faccanoni

Date de publication

06 mars 2026

Code 
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
from scipy.optimize import fsolve
from scipy.integrate import solve_ivp, odeint
from IPython.display import display, Markdown, Math

EXERCICE 1 - sujet A : 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.

Le modèle s’écrit \( \begin{cases} S'(t)=-\beta S(t)Z(t)\\ Z'(t)=\beta S(t)Z(t)+\zeta R(t)-\alpha S(t) Z(t)\\ R'(t)=\alpha S(t)Z(t)-\zeta R(t) \end{cases} \)

  1. Soit \(P(t)=S(t)+Z(t)+R(t)\). Montrer que \(P\) est un invariant.
  2. 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)\).
  3. Même exercice avec la méthode d’Euler Régressif.
  4. 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]

Dans la suite, considérons les données initiales :

  • \(S(0)=0.95\)
  • \(Z(0)=0.05\)
  • \(R(0)=0\)

On simulera pour \(t\in[0;60]\) jours.

Code 
alpha = 1
beta  = 0.5
zeta  = 0.25  

# yy=[S,I,R]
phi0 = lambda S,Z,R,t :  -beta*S*Z
phi1 = lambda S,Z,R,t :   (beta-alpha)*S*Z+zeta*R
phi2 = lambda S,Z,R,t :   alpha*S*Z-zeta*R

S0 = 0.95
Z0 = 1-S0
R0 = 0 

tt = np.linspace(0,60,61) # en jours
h = tt[1]-tt[0]

\(P'(t)=S'(t)+Z'(t)+R'(t)=-\beta S(t)Z(t)+\beta S(t)Z(t)+\zeta R(t)-\alpha S(t) Z(t)+\alpha S(t)Z(t)-\zeta R(t)=0\).

On notera \(S_n\approx S(t_n)\), \(Z_n\approx Z(t_n)\) et \(R_n\approx R(t_n)\).

Euler explicite \( \begin{cases} S_{n+1}=S_n+h\varphi_0(S_n,Z_n,R_n,t_n),\\ Z_{n+1}=Z_n+h\varphi_1(S_n,Z_n,R_n,t_n),\\ R_{n+1}=R_n+h\varphi_2(S_n,Z_n,R_n,t_n). \end{cases} \)

Code 
# S0, Z0, R0, h variables globales
def euler_progressif(phi0,phi1,phi2,tt):
    SS = [S0]
    ZZ = [Z0]
    RR = [R0]
    for i in range(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)

plt.figure(figsize=(10,5))
PP_ep = [SS_ep[i]+ZZ_ep[i]+RR_ep[i] for i in range(len(tt))]
plt.plot(tt,SS_ep,'--',tt,ZZ_ep,'.-',tt,RR_ep,':',tt,PP_ep)
plt.xlabel('t')
plt.legend(['S(t)','Z(t)','R(t)','P(t)'])
plt.title('Euler progressif') 
plt.grid()
plt.axis([0,60,0,1.1]);

Euler implicite \( \begin{cases} S_{n+1}=S_n+h\varphi_0(S_{n+1},Z_{n+1},R_{n+1},t_{n+1}),\\ Z_{n+1}=Z_n+h\varphi_2(S_{n+1},Z_{n+1},R_{n+1},t_{n+1}),\\ R_{n+1}=R_n+h\varphi_3(S_{n+1},Z_{n+1},R_{n+1},t_{n+1}). \end{cases} \)

Code 
# S0, Z0, R0, h variables globales
def euler_regressif(phi0,phi1,phi2,tt):
    SS = [S0]
    ZZ = [Z0]
    RR = [R0]
    for i in range(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)

plt.figure(figsize=(10,5))
PP_er=[SS_er[i]+ZZ_er[i]+RR_er[i] for i in range(len(tt))]
plt.plot(tt,SS_er,'--',tt,ZZ_er,'.-',tt,RR_er,':',tt,PP_er)
plt.xlabel('t')
plt.legend(['S(t)','Z(t)','R(t)','P(t)'])
plt.title('Euler regressif') 
plt.grid()
plt.axis([tt[0],tt[-1],0,1.1]);

odeint

Code 
pphi = 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]

plt.figure(figsize=(10,5))
PP = [SS[i]+ZZ[i]+RR[i] for i in range(len(tt))]
plt.plot(tt,SS,'--',tt,ZZ,'.-',tt,RR,':',tt,PP)
plt.xlabel('t')
plt.legend(['S(t)','Z(t)','R(t)','P(t)'])
plt.title('Odeint') 
plt.grid()
plt.axis([0,60,0,1.1]);

EXERCICE 1 - sujet B : implémentation d’un schéma pour un système

Le modèle S.I.R. a été présenté pour la première fois par KERMACK & McKENDRICK à Londres et Cambridge en 1927 pour expliquer a posteriori l’évolution de l’épidémie de peste à Bombay en 1905-1906.

À 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 « Infectés » (I) : les malades, ce sont aussi les contagieux (c’est une hypothèse de ce modèle) ; - les individus « Rétablis » (R, comme « Recovered » en anglais) : ceux qui ont déjà eu la maladie et sont désormais immunisés contre cette maladie. Dans ce modèle on inclut dans ce groupe les personnes décédées (puisqu’elles ne peuvent plus contracter la maladie).

On travaillera avec des variables normalisées : leur valeur est comprise entre \(0\) et \(1\) et représente une fraction de la population totale (P).

Le modèle s’écrit \( \begin{cases} S'(t)=-\beta S(t)I(t)\\ I'(t)=\beta S(t)I(t)-\nu I(t)\\ R'(t)=\nu I(t) \end{cases} \)

  1. Soit \(P(t)=S(t)+I(t)+R(t)\). Montrer que \(P\) est un invariant.
  2. Calculer la solution approchée obtenue par la méthode d’Euler Progressif avec \(91\) points. Afficher ensuite sur le même repère \(t\mapsto S(t)\), \(t\mapsto I(t)\), \(t\mapsto R(t)\) et \(t\mapsto P(t)\).
  3. Même exercice avec la méthode d’Euler Régressif.
  4. Même exercice avec la fonction odeint du module scipy.

Dans la suite, considérons les données initiales :

  • \(S(0)=0.99\)
  • \(I(0)=0.01\)
  • \(R(0)=0\)

Pour la rougeole on prendra

  • \(\beta=0.8\) = taux de transmission (plus il est grand, plus la maladie est infectieuse)
  • \(\nu=0.05\) = taux de rémission ou de décès du virus ( = 1 / (duréé moyenne de l’infection); plus \(\nu\) est petit, plus le temps durant lequel le malade est contagieux est grand )

On simulera pour \(t\in[0;90]\) jours.

Code 
beta = 0.8
nu   = 0.05  

# yy=[S,I,R]
phi0 = lambda S,I,R,t :  -beta*S*I
phi1 = lambda S,I,R,t :   beta*S*I-nu*I
phi2 = lambda S,I,R,t :   nu*I

S0=0.99
I0=1-S0
R0=0 

tt = np.linspace(0,90,91) # en jours
h = tt[1]-tt[0]

\(P'(t)=S'(t)+I'(t)+R'(t)=-\beta S(t)I(t)+\beta S(t)I(t)-\nu I(t)+\nu I(t)=0\).

On notera \(S_n\approx S(t_n)\), \(I_n\approx I(t_n)\) et \(R_n\approx R(t_n)\).

Euler explicite \( \begin{cases} S_{n+1}=S_n+h\varphi_0(S_n,I_n,R_n,t_n),\\ I_{n+1}=I_n+h\varphi_1(S_n,I_n,R_n,t_n),\\ R_{n+1}=R_n+h\varphi_2(S_n,I_n,R_n,t_n). \end{cases} \)

Code 
# S0, I0, R0, h variables globales
def euler_progressif(phi0,phi1,phi2,tt):
    SS = [S0]
    II = [I0]
    RR = [R0]
    for i in range(len(tt)-1):
        SS.append(SS[i]+h*phi0(SS[i],II[i],RR[i],tt[i]))
        II.append(II[i]+h*phi1(SS[i],II[i],RR[i],tt[i]))
        RR.append(RR[i]+h*phi2(SS[i],II[i],RR[i],tt[i]))
    return [SS,II,RR]

[SS_ep, II_ep, RR_ep] = euler_progressif(phi0,phi1,phi2,tt)

plt.figure(figsize=(10,5))
PP_ep=[SS_ep[i]+II_ep[i]+RR_ep[i] for i in range(len(tt))]
plt.plot(tt,SS_ep,'--',tt,II_ep,'.-',tt,RR_ep,':',tt,PP_ep)
plt.xlabel('t')
plt.legend(['S(t)','I(t)','R(t)','P(t)'])
plt.title('Euler progressif') 
plt.grid()
plt.axis([tt[0],tt[-1],0,1.1]);

Euler implicite \( \begin{cases} S_{n+1}=S_n+h\varphi_0(S_{n+1},I_{n+1},R_{n+1},t_{n+1}),\\ I_{n+1}=I_n+h\varphi_1(S_{n+1},I_{n+1},R_{n+1},t_{n+1}),\\ R_{n+1}=R_n+h\varphi_2(S_{n+1},I_{n+1},R_{n+1},t_{n+1}). \end{cases} \)

Code 
from scipy.optimize import fsolve
# S0, I0, R0, h variables globales
def euler_regressif(phi0,phi1,phi2,tt):
    SS = [S0]
    II = [I0]
    RR = [R0]
    for i in range(len(tt)-1):
        sys = lambda z : [ -z[0]+SS[i]+h*phi0(z[0],z[1],z[2],tt[i+1]) , -z[1]+II[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,Itemp, Rtemp = fsolve( sys , (SS[i],II[i],RR[i]) ) 
        SS.append(Stemp)
        II.append(Itemp)
        RR.append(Rtemp)
    return [SS,II,RR]

[SS_er, II_er, RR_er] = euler_regressif(phi0,phi1,phi2,tt)

plt.figure(figsize=(10,5))
PP_er=[SS_er[i]+II_er[i]+RR_er[i] for i in range(len(tt))]
plt.plot(tt,SS_er,'--',tt,II_er,'.-',tt,RR_er,':',tt,PP_er)
plt.xlabel('t')
plt.legend(['S(t)','I(t)','R(t)','P(t)'])
plt.title('Euler regressif') 
plt.grid()
plt.axis([tt[0],tt[-1],0,1.1]);

odeint

Code 
pphi = 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,I0,R0]
sol  = odeint(pphi,CI,tt)

SS,II,RR = sol[:,0],sol[:,1],sol[:,2]

plt.figure(figsize=(10,5))
PP = [SS[i]+II[i]+RR[i] for i in range(len(tt))]
plt.plot(tt,SS,'--',tt,II,'.-',tt,RR,':',tt,PP)
plt.xlabel('t')
plt.legend(['S(t)','I(t)','R(t)','P(t)'])
plt.title('Odeint') 
plt.grid()
plt.axis([tt[0],tt[-1],0,1.1]);

EXERCICE 2 - sujet A : é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\gamma u_{n} +\left(1-2\gamma \right)u_{n-1} +h\left(\gamma \varphi_{n+1} +\left(2-3\gamma \right)\varphi_{n}\right). \)

  1. Pour quelles valeurs de \( est-elle consistante ?
  2. Pour quelles valeurs de \(\gamma\) est-elle convergente ?
  3. Quel ordre de convergente a-t-elle ?
  4. 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.

\(p=1\): c’est une méthode à \(q=p+1=2\) : - $a_0=2$ et $a_1=1-2$ - $b_0=2-3$ et \(b_1=0\) - $b_{-1}=$

La méthode est implicite si \(\gamma\neq0\).

La méthode est consistante pour tout \(\gamma\) car \( \begin{cases} \displaystyle\sum_{j=0}^p a_j=a_0+a_1=2\gamma+(1-2\gamma)=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)=-(1-2\gamma)+(\gamma+2-3\gamma)=1 \end{cases} \)

Vérifions ces calculs ci-dessous :

Code 
sp.var('gamma')
p=1
a0=2*gamma
a1=1-2*gamma
b0=2-3*gamma
b1=0
bm1=gamma

cond11=a0+a1
cond12=-(0*a0+1*a1)+(bm1+b0+b1)
display(cond11)
display(cond12.factor())
display(sp.solve([sp.Eq(cond11,1),sp.Eq(cond12,1)],gamma))

\(\displaystyle 1\)

\(\displaystyle 1\)

[]

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\gamma r+(2\gamma-1) \) dont les racines sont \( r_0=1,\quad r_1=2\gamma-1. \)

Code 
sp.var('r')
rho=r**2-a0*r-a1
display(sp.factor(rho))
sol=sp.solve(rho,r)
display(sol)

\(\displaystyle \left(r - 1\right) \left(- 2 \gamma + r + 1\right)\)

[1, 2*gamma - 1]

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\le\gamma<1\).

Code 
cond=sp.solve([abs(c)<=1 for c in sol],gamma) 
print('La première condition impose')
display(cond) 
La première condition impose

\(\displaystyle 0 \leq \gamma \wedge \gamma \leq 1\)

Conclusion : la méthode est convergente ssi elle est consistante et zéro-stable donc ssi \(0\le\gamma<1\).

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(1-2\gamma\big)+2\big(\gamma\big) {}= 1 \)

Code 
cond2=(0*a0+1*a1)+2*(-(-1)*bm1+0*b0+1*b1)
display(cond2.factor())
display(sp.solve(sp.Eq(cond2,1),gamma))

\(\displaystyle 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) {}=-(1-2\gamma)+3\big(\gamma\big) {}=5\gamma-1=1 \iff \gamma=\frac{2}{5} \)

Code 
cond3=(0*a0-1*a1)+3*(1*bm1+0*b0+1*b1)
display(cond3.factor())
display(sp.solve(sp.Eq(cond3,1),gamma))

\(\displaystyle 5 \gamma - 1\)

[2/5]

Conclusion

  • 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{2}{5}\) (et on peut montrer qu’il n’est pas d’ordre 4)

On définit la solution exacte (calculée en utilisant le module sympy) pour estimer les erreurs :

Code 
# variables globales
t0     = 0
tfinal = 1
y0     = 2

phi = lambda t,y : y*(1-y)


t  = sp.Symbol('t')
y  = sp.Function('y')
edo = sp.Eq( sp.diff(y(t),t) , phi(t,y(t)) )
solpar = sp.dsolve(edo,y(t),ics={y(t0):y0})
display(Markdown(f'La solution exacte est : \){sp.latex(solpar)}\('))

sol_exacte = sp.lambdify(t,solpar.rhs,'numpy')

La solution exacte est : \(y{\left(t \right)} = \frac{1}{1 - \frac{e^{- t}}{2}}\)

On calcule la solution approchée pour différentes valeurs de \(N\):

Code 
def multipas(phi, tt, gamma):
    h  = tt[1] - tt[0]
    uu = [y0]
    uu.append(sol_exacte(tt[1]))
    for i in range(1,len(tt) - 1):
        temp = fsolve( lambda x : -x+2*gamma*uu[i]+(1-2*gamma)*uu[i-1]+h*(gamma*phi(tt[i+1],x)+(2-3*gamma)*phi(tt[i],uu[i]))   ,uu[i])
        uu.append( temp[0] )
    return uu


H = np.array([])
err_mp_2 = np.array([])
err_mp_3 = np.array([])
N = 10
for k in range(7):
    N += 20
    tt = np.linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = sol_exacte(tt)
    uu_mp_2 = multipas(phi, tt, 1/5)
    uu_mp_3 = multipas(phi, tt, 2/5)
    H = np.append(H,h)
    err_mp_2 = np.append(err_mp_2,np.linalg.norm(uu_mp_2-yy,np.inf))
    err_mp_3 = np.append(err_mp_3,np.linalg.norm(uu_mp_3-yy,np.inf))
    
ordre_mp_2 = np.polyfit(np.log(H),np.log(err_mp_2), 1)[0]
ordre_mp_3 = np.polyfit(np.log(H),np.log(err_mp_3), 1)[0]
print ('Multipas gamma=1/5 ordre=%1.2f' %(ordre_mp_2))
print ('Multipas gamma=2/5 ordre=%1.2f' %(ordre_mp_3))


plt.figure(figsize=(8,5))
plt.plot(np.log(H), np.log(err_mp_2), 'r-o', label=r'MPL avec $\gamma=\frac{1}{5}$')
plt.plot(np.log(H), np.log(err_mp_3), 'b-o', label=r'MPL avec $\gamma=\frac{2}{5}$')
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
plt.legend(bbox_to_anchor=(1.04, 1), loc='upper left')
plt.grid(True)
plt.show()

Multipas gamma=1/5 ordre=1.97
Multipas gamma=2/5 ordre=2.94

EXERCICE 2 - sujet B : é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}. \)

  1. Pour quelles valeurs de \) est-elle consistante ?
  2. Pour quelles valeurs de \(\gamma\) est-elle convergente ?
  3. Quel ordre de convergente a-t-elle ?
  4. 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.

\(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 
sp.var('gamma')
p=1
a0=2*(1-gamma)
a1=2*gamma-1
b0=3*gamma-1
b1=0
bm1=1-gamma

cond11=a0+a1
cond12=-(0*a0+1*a1)+(bm1+b0+b1)
display(cond11)
display(cond12.factor())
display(sp.solve([sp.Eq(cond11,1),sp.Eq(cond12,1)],gamma))

\(\displaystyle 1\)

\(\displaystyle 1\)

[]

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. \)

Code 
sp.var('r')
rho=r**2-a0*r-a1
display(sp.factor(rho))
sol=sp.solve(rho,r)
display(sol)

\(\displaystyle \left(r - 1\right) \left(2 \gamma + r - 1\right)\)

[1, 1 - 2*gamma]

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=sp.solve([abs(c)<=1 for c in sol],gamma)
display(cond)

\(\displaystyle 0 \leq \gamma \wedge \gamma \leq 1\)

Conclusion : la méthode est convergente ssi elle est consistante et zéro-stable donc ssi \(0<\gamma\le1\).

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 \)

Code 
cond2=(0*a0+1*a1)+2*(-(-1)*bm1+0*b0+1*b1)
display(cond2.factor())
display(sp.solve(sp.Eq(cond2,1),gamma))

\(\displaystyle 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} \)

Code 
cond3=(0*a0-1*a1)+3*(1*bm1+0*b0+1*b1)
display(cond3.factor())
display(sp.solve(sp.Eq(cond3,1),gamma))

\(\displaystyle 4 - 5 \gamma\)

[3/5]

Conclusion

  • 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)

On définit la solution exacte (calculée en utilisant le module sympy) pour estimer les erreurs:

Code 
# variables globales
t0     = 0
tfinal = 1
y0     = 2

phi = lambda t,y : y*(1-y)


t  = sp.Symbol('t')
y  = sp.Function('y')
edo= sp.Eq( sp.diff(y(t),t) , phi(t,y(t)) )
solpar = sp.dsolve(edo, y(t), ics={y(t0):y0})
display(Markdown(f'La solution exacte est : \({sp.latex(solpar)}\)'))

sol_exacte = sp.lambdify(t,solpar.rhs,'numpy')

La solution exacte est : \(y{\left(t \right)} = \frac{1}{1 - \frac{e^{- t}}{2}}\)

On calcule la solution approchée pour différentes valeurs de \(N\):

Code 
def multipas(phi, tt, gamma):
    h  = tt[1] - tt[0]
    uu = np.empty_like(tt)
    uu[0] = y0
    uu[1] = sol_exacte(tt[1])
    for i in range(1,len(tt) - 1):
        temp = fsolve( lambda x : -x+2*(1-gamma)*uu[i]+(2*gamma-1)*uu[i-1]+h*(1-gamma)*phi(tt[i+1],x)+h*(3*gamma-1)*phi(tt[i],uu[i])   ,uu[i])
        uu[i+1] = temp[0]
    return uu

H = np.array([])
err_mp_2 = np.array([])
err_mp_3 = np.array([])
N = 10
for k in range(7):
    N += 20
    tt = np.linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = sol_exacte(tt)
    uu_mp_2 = multipas(phi, tt, 1/5)
    uu_mp_3 = multipas(phi, tt, 3/5)
    H = np.append(H,h)
    err_mp_2 = np.append(err_mp_2,np.linalg.norm(uu_mp_2-yy,np.inf))
    err_mp_3 = np.append(err_mp_3,np.linalg.norm(uu_mp_3-yy,np.inf))
    
ordre_mp_2 = np.polyfit(np.log(H),np.log(err_mp_2), 1)[0]
ordre_mp_3 = np.polyfit(np.log(H),np.log(err_mp_3), 1)[0]

print (f'MPL avec gamma=1/5 ordre={ordre_mp_2:.2f}')
print (f'MPL avec gamma=3/5 ordre={ordre_mp_3:.2f}')

plt.figure(figsize=(8,5))
plt.plot(np.log(H), np.log(err_mp_2), 'r-o', label=r'Multipas $\gamma=\frac{1}{5}$')
plt.plot(np.log(H), np.log(err_mp_3), 'b-o', label=r'Multipas $\gamma=\frac{3}{5}$')
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
plt.legend(bbox_to_anchor=(1.04, 1), loc='upper left')
plt.grid(True)
plt.show()

MPL avec gamma=1/5 ordre=1.95
MPL avec gamma=3/5 ordre=2.94