M62 Étude de schémas numériques pour l’approximation de problèmes de Cauchy - 2021 - Contrôle Terminal session 1

2021 - Contrôle Terminal session 1

Auteur·rice

Gloria FACCANONI

Date de publication

25 mars 2024

1 Exercice 1 : étude d’un schéma RK

Soit le schéma de Runge-Kutta dont la matrice de Butcher est \( \begin{array}{c|cc} 0 & 0 & 0\\ \dfrac{\eta}{4} & \dfrac{\eta}{4} & 0\\ \hline & \dfrac{3-\mu}{3} & \dfrac{\mu}{3} \end{array} \)

  1. Écrire le schéma sous la forme d’une suite définie par récurrence.
    Le schéma est explicite, semi-implicite ou implicite ?
  2. Sans faire de calculs, indiquer quel est l’ordre maximale du schéma.
  3. Étudier théoriquement l’ordre du schéma en fonction de \(\eta\) et \(\mu\).
  4. Étudier théoriquement la A-stabilité du schéma pour les valeurs de \(\eta\) et \(\mu\) qui donnent l’ordre maximal.
  5. Implémenter le schéma avec les valeurs de \(\eta\) et \(\mu\) trouvées avec vos nom et prénom et approcher la solution du problème de Cauchy suivant avec \(N+1\) points (quelles valers de \(N\) peut-on choisir? justifier la réponse) \( \begin{cases} y'(t)=-50y(t), &t\in[0;1]\\ y(0)=1 \end{cases} \) Vérifier ensuite empiriquement l’ordre de convergence.
Code 
%reset -f
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

import sympy 
sympy.init_printing()

from IPython.display import display, Math, Markdown
prlat= lambda *args: display(Math(''.join( sympy.latex(arg) for arg in args )))

Q1
Écrire le schéma sous la forme d’une suite définie par récurrence.
Le schéma est explicite, semi-implicite ou implicite ?

Le schéma associé à cette matrice est explicite (car la matrice \(\mathbb{A}\) est triangulaire inférieure à diagonale nulle).
On peut calculer \(u_{n+1}\) à partir de \(u_n\) par la formule de récurrence

\( \begin{cases} u_0 = y_0 \\ K_1 = \varphi\left(t_n,u_n\right)\\ K_2 = \varphi\left(t_{n}+\dfrac{\eta}{4}h,u_n+\dfrac{\eta}{4}hK_1\right)\\ u_{n+1} = u_n + \dfrac{h}{3}\left((3-\mu) K_1+\mu K_2\right) & n=0,1,\dots N-1 \end{cases} \)

Q2
Sans faire de calculs, indiquer quel est l’ordre maximale du schéma.

Soit \(\omega\) l’ordre de la méthode et \(s=2\) le nombre d’étages. Étant un schéma explicite, \(\omega\le s=2\).

Q3
Étudier théoriquement l’ordre du schéma en fonction de \(\eta\) et \(\mu\).

  • Si \( \begin{cases} \displaystyle\sum_{j=1}^s b_{i}=1& \\ \displaystyle c_i=\sum_{j=1}^s a_{ij}&i=1,\dots,s \end{cases} \) alors \(\omega\ge1\)

  • Si \(\omega\ge1\) et \(\displaystyle\sum_{j=1}^s b_j c_j=\frac{1}{2}\) alors \(\omega\ge2\)

D’après les calculs ci-dessous il est toujours au moins d’ordre 1 et il est d’ordre 2 ssi \(\eta\mu=6\).

Code 
def ordre_RK(s, A=None, b=None, c=None):

    j = sympy.symbols('j')
    
    if A is None:
        A = sympy.MatrixSymbol('a', s, s)
    else:
        A = sympy.Matrix(A)
    if c is None:
        c = sympy.symbols(f'c_0:{s}')    
    if b is None:
        b = sympy.symbols(f'b_0:{s}')
    
    
    display(Markdown("**Matrice de Butcher**"))
    But = matrice_Butcher(s, A, b, c)
    
    Eqs = []

    display(Markdown(f"**On a {s+1} conditions pour avoir consistance = pour être d'ordre 1**"))
    ordre_1(s, A, b, c, j)
    
    display(Markdown("**On doit ajouter 1 condition pour être d'ordre 2**"))
    Eq_2 = ordre_2(s, A, b, c, j)
    Eqs.append(Eq_2)

    return Eqs


def matrice_Butcher(s, A, b, c):
    But = sympy.Matrix(A)
    But = But.col_insert(0, sympy.Matrix(c))   
    last = [sympy.Symbol(" ")]
    last.extend(b)
    But = But.row_insert(s,sympy.Matrix(last).T)
    display(Math(sympy.latex(sympy.Matrix(But))))
    return But

def ordre_1(s, A, b, c, j):
    texte = "\sum_{j=1}^s b_j =" + f"{sum(b).simplify()}"
    texte += r"\text{ doit être égale à }1"
    display(Math(texte))
    for i in range(s):
        somma = sympy.summation(A[i,j],(j,0,s-1)).simplify()
        texte = r'\sum_{j=1}^s a_{'+str(i)+r'j}=' + sympy.latex( somma )
        texte += r"\text{ doit être égale à }"+sympy.latex(c[i])
        display(Math( texte ))
    return None    

def ordre_2(s, A, b, c, j):
    texte = '\sum_{j=1}^s b_j c_j=' +sympy.latex(sum([b[i]*c[i] for i in range(s)]).simplify())
    texte += r"\text{ doit être égale à }\frac{1}{2}"
    display(Math(texte))
    Eq_2 = sympy.Eq( sum([b[i]*c[i] for i in range(s)]).simplify(), sympy.Rational(1,2) )
    return Eq_2




##################################################
# Conditions avec s étages
##################################################

eta = sympy.Symbol(r'\eta')
mu = sympy.Symbol(r'\mu')

c = [0,sympy.S(eta)/4]
b = [sympy.S(3-mu)/3,sympy.S(mu)/3]
A = [[0,0],[sympy.S(eta)/4,0]]
s = len(c)


Eqs = ordre_RK(s, A, b, c)

display(Markdown("**Solutions**"))

sol = sympy.solve(Eqs)[0]
sol
# for eq in Eqs:
#     display(eq)

Matrice de Butcher

\(\displaystyle \left[\begin{matrix}0 & 0 & 0\\\frac{\eta}{4} & \frac{\eta}{4} & 0\\ & 1 - \frac{\mu}{3} & \frac{\mu}{3}\end{matrix}\right]\)

On a 3 conditions pour avoir consistance = pour être d’ordre 1

\(\displaystyle \sum_{j=1}^s b_j =1\text{ doit être égale à }1\)

\(\displaystyle \sum_{j=1}^s a_{0j}=0\text{ doit être égale à }0\)

\(\displaystyle \sum_{j=1}^s a_{1j}=\frac{\eta}{4}\text{ doit être égale à }\frac{\eta}{4}\)

On doit ajouter 1 condition pour être d’ordre 2

\(\displaystyle \sum_{j=1}^s b_j c_j=\frac{\eta \mu}{12}\text{ doit être égale à }\frac{1}{2}\)

Solutions

\(\displaystyle \left\{ \eta : \frac{6}{\mu}\right\}\)

Q4
Étudier théoriquement la A-stabilité du schéma pour les valeurs de \(\eta\) et \(\mu\) qui donnent l’ordre maximal.

Soit \(\beta>0\) un nombre réel positif et considérons le problème de Cauchy \(y'(t)=-\beta y(t)\) pour \(t>0\) et \(y(0)=1\). Sa solution est \(y(t)=e^{-\beta t}\) donc \(\lim_{t\to+\infty}|y(t)|=0.\)

Le schéma appliqué à ce problème de Cauchy s’écrit
\( \begin{cases} u_0 = y_0 \\ K_1 = -\beta u_n \\ K_2 = -\beta\left(u_n+\dfrac{\eta}{4}hK_1\right)\\ u_{n+1} = u_n + \dfrac{h}{3}\left((3-\mu)K_1+\mu K_2\right) & n=0,1,\dots N-1 \end{cases} \)

Dans le cas optimal on a \(\eta\mu=6\), ainsi
\(u_{n+1} = \left(\dfrac{1}{2}(\beta h)^2 -(\beta h)+1\right)u_n\)

Code 
u_n = sympy.Symbol('u_n')
dt = sympy.Symbol('h') # je change le nom de la variable "h" en "dt" pour ne pas effacer l'affectation de "h"
beta = sympy.Symbol(r'\beta')

K1 = -beta*u_n
display(Math('K_1='+sympy.latex(K1)))
K2 = -beta*(u_n+sympy.S(eta)/4*dt*K1).factor()
display(Math('K_2='+sympy.latex(K2)))
u_np1 = (u_n+dt*(sympy.S(3-mu)/3*K1+sympy.S(mu)/3*K2)).factor()
display(Math('u_{n+1}='+sympy.latex(u_np1)))

display(Markdown("Sachant que"))
display(sol)
display(Markdown("on a"))
u_np1 = u_np1.subs(sol)
display(Math('u_{n+1}='+sympy.latex(u_np1)))

\(\displaystyle K_1=- \beta u_{n}\)

\(\displaystyle K_2=\frac{\beta u_{n} \left(\beta \eta h - 4\right)}{4}\)

\(\displaystyle u_{n+1}=\frac{u_{n} \left(\beta^{2} \eta \mu h^{2} - 12 \beta h + 12\right)}{12}\)

Sachant que

\(\displaystyle \left\{ \eta : \frac{6}{\mu}\right\}\)

on a

\(\displaystyle u_{n+1}=\frac{u_{n} \left(6 \beta^{2} h^{2} - 12 \beta h + 12\right)}{12}\)

On note \(x=\beta h\) et on étudie la fonction \(q\colon \mathbb{R}^+\to\mathbb{R}\) définie par \(q(x)=\frac{1}{2}x^2-x+1\) car le schéma est A-stable ssi \(|q(x)|<1\).
Il s’agit d’une parabole convexe de sommet en \(x_s=1\) et \(q(x_s)=\frac{1}{2}>-1\) donc \(q(x)<1\) pour \(0<x<2\) et \(q(x)>0\) pour tout \(x\). On conclut que le schéma est A-stable pour tout \(h<\dfrac{2}{\beta}\).

On peut vérifier les calculs avec sympy:

Code 
x = sympy.Symbol('x',nonnegative=True)

q = (u_np1/u_n).subs(beta*dt,x)
display(Math('q(x)='+sympy.latex(q)))

display(Markdown("$|q(x)|<1$ pour"))
display(sympy.solve(abs(q)<1))

sympy.plot(q,xlim=[0,3],ylim=[-1.5,1.5]);


# display(Math('q(0)='+sympy.latex(q.subs(x,0))))
# lim=sympy.Limit(q,x,sympy.oo)
# display(Math(sympy.latex(lim)+'='+sympy.latex(lim.doit())))

# print('On cherche les points stationnaires dans R^+')
# dq=(q.diff(x)).factor()
# display(Math("q'(x)="+sympy.latex(dq)))
# print("q est une parabole dont le sommet est")
# sol=sympy.solve(dq)
# display(Math("q'(x)=0 \iff x\in"+sympy.latex(sol)))
# minimum=sol[0]
# display(Math("x="+sympy.latex(minimum)+"\\text{ est un minimum et}"))
# qmin=q.subs(x,sol[0])
# display(Math("q("+sympy.latex(minimum)+")="+sympy.latex(qmin)))
# sympy.plot(q,1,-1,xlim=[-1,15],ylim=[-1.5,1.5]);
# print("Conclusion: -1<q(x)<1 pour tout x<2. Vérifions ce calcul:")
# print("q(x)<1")
# display(sympy.solve(q<1))
# print("q(x)>-1")
# display(sympy.solve(q>-1))
# print("Notons que q(x)>0 pour tout 0<x<2")
# #display(sympy.solve(q>0))

\(\displaystyle q(x)=\frac{x^{2}}{2} - x + 1\)

\(|q(x)|<1\) pour

\(\displaystyle 0 < x \wedge x < 2\)

Pour la question 5 on va fixer \(\eta\) et \(\mu\). Pour cela, écrire votre nom et prénom dans la cellule ci-dessous et vous obtiendrez un choix pour les paramètres.

Code 
nom = "Faccanoni"
prenom = "Gloria"

L = [1,2,3,6]
idx = sum([ord(c) for c in nom+prenom])%len(L)
display(Math(r'\eta='+sympy.latex(L[idx])))
display(Math(r'\mu='+sympy.latex(sympy.S(6)/L[idx])))

eta, mu = L[idx], 6/L[idx]

\(\displaystyle \eta=1\)

\(\displaystyle \mu=6\)

Q5
Implémenter le schéma avec les valeurs de \(\eta\) et \(\mu\) trouvées avec vos nom et prénom et approcher la solution du problème de Cauchy suivant avec \(N+1\) points (quelles valers de \(N\) peut-on choisir? justifier la réponse) \( \begin{cases} y'(t)=-50y(t), &t\in[0;1]\\ y(0)=1 \end{cases} \)
Vérifier ensuite empiriquement l’ordre de convergence.

Code 
# ==================================================
# SCHEMA DE RUNGE-KUTTA
# ==================================================

# def RK(tt,phi,y0): # comme suite geometrique
#     x = beta*h
#     q = x**2/2-x+1
#     uu = y0*q**np.arange(len(tt))
#     return uu

def RK(tt,phi,y0):
    uu = np.zeros_like(tt)
    uu[0] = y0
    h = tt[1]-tt[0]
    for n in range(len(tt)-1):
        K1 = phi(tt[n],uu[n])
        K2 = phi(tt[n]+eta/4*h,uu[n]+eta/4*h*K1)
        uu[n+1] = uu[n]+h/3*((3-mu)*K1+mu*K2) 
    return uu


# ==================================================
# EDO
# ==================================================
t0, y0, tfinal = 0, 1, 1
beta = 50

phi = lambda t,y : -beta*y
sol_exacte = lambda t : np.exp(-beta*t)


# ==================================================
# COMPARAISON POUR N DONNE
# ==================================================
plt.figure(figsize=(20,5))

plt.subplot(1,2,1)

N = 50
tt = np.linspace(t0,tfinal,N+1)
h = tt[1]-tt[0]

if h>=2/beta:
    print('Attention: h > limite de stabilité')
    
yy = sol_exacte(tt)   
# uu = RK(tt)
uu = RK(tt,phi,y0)
plt.plot(tt,yy,'b-',label=("Exacte"))
plt.plot(tt,uu,'ro',label=("RK"))
plt.title(f'${N=}$, ${h=} < {2/beta}$')
plt.xlabel('$t$')
plt.ylabel('$u$')
plt.legend() 
plt.grid(True)



# ==================================================
# ORDRE DE CONVERGENCE
# ==================================================

plt.subplot(1,2,2)

H = []
err = []
N = 50
for k in range(7):
    N += 100
    tt = np.linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = sol_exacte(tt)
    uu = RK(tt,phi,y0)
    H.append(h)
    err.append(np.linalg.norm(yy-uu, np.inf))
    
omega = np.polyfit(np.log(H),np.log(err), 1)[0]
plt.plot(np.log(H), np.log(err), 'c-o')
plt.xlabel('$\ln(h)$')
plt.ylabel('$\ln(e)$')
plt.title(rf'Ordre $\omega$={omega:1.2f}')
plt.grid(True);

2 Exercice 2 : choix d’un schéma

Soit le problème de Cauchy \( \begin{cases} y'(t)=-10^2 \Big(y(t)-\cos(t)\Big), &t\in[0;2]\\ y(0)=0 \end{cases} \)

  1. Calculer la solution exacte.

  2. Choisir le schéma le plus adapté parmi les suivants (en justifiant la réponse)

    • EE (Euler Explicite),
    • EI (Euler Implicite),
    • EM (Euler Modifié),
    • H (Heun).

  3. Approcher ensuite la solution du problème donné avec \(N=40\) avec le schéma choisi (on affichera la solution exacte et la solution approchée sur le même repère).

Code 
%reset -f
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

import sympy 
sympy.init_printing()

from IPython.display import display, Math, Markdown
prlat= lambda *args: display(Math(''.join( sympy.latex(arg) for arg in args )))
  • On définit la solution exacte (calculée en utilisant le module sympy):
Code 
t  = sympy.Symbol('t')
y  = sympy.Function('y')

t0     = 0
tfinal = 2
y0     = 0

# NB il faut utiliser le cos du module sympy et non de numpy
edo = sympy.Eq( sympy.diff(y(t),t) , -50*(y(t)-sympy.cos(t)) )
display(edo)
solpar = sympy.dsolve(edo, y(t), ics={y(t0):y0})
display(solpar)

sympy.plot(solpar.rhs, (t, t0, tfinal), ylabel='$y(t)$', xlabel='$t$');

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

\(\displaystyle \frac{d}{d t} y{\left(t \right)} = - 50 y{\left(t \right)} + 50 \cos{\left(t \right)}\)

\(\displaystyle y{\left(t \right)} = \frac{50 \sin{\left(t \right)}}{2501} + \frac{2500 \cos{\left(t \right)}}{2501} - \frac{2500 e^{- 50 t}}{2501}\)

  • Il s’agit d’une EDO stiff (le terme en exponentielle décroît très rapidement): on ne peut utiliser ni les schémas EE et EM ni le schéma H car la condition de stabilité est trop contraignante. Le schéma EI, en revanche, est inconditionnellement A-stable: on peut donc choisir le pas \(h\) sans contrainte de stabilité.

  • On calcule la solution approchée

Code 
def EE(phi, tt, y0):
    h = tt[1] - tt[0]
    uu =np.zeros_like(tt)
    uu[0] = y0
    for n in range(len(tt) - 1):
        uu[n+1] = uu[n] + h * phi(tt[n], uu[n])
    return uu
    
def EI(phi, tt, y0):
    h = tt[1] - tt[0]
    uu =np.zeros_like(tt)
    uu[0] = y0
    for n in range(len(tt) - 1):
        temp = fsolve(lambda x: -x + uu[n] + h * phi(tt[n + 1], x), uu[n])
        uu[n+1] = temp[0]
    return uu


phi = lambda t,y : -50*(y-np.cos(t))

N = 40
tt = np.linspace(t0, tfinal, N + 1)
h = tt[1] - tt[0]
yy = sol_exacte(tt)

# uu = EE(phi, tt, y0) # pour montrer les pb de stabilité avec le même nombre de points
uu = EI(phi, tt, y0)
err = np.linalg.norm(yy-uu, np.inf)

plt.figure(figsize=(8,5))
plt.plot(tt, uu, 'o', label='EI')
plt.plot(tt, yy, 'r-', label='Exacte')
plt.xlabel('$t$')
plt.ylabel('$y$')
plt.grid(True)
plt.legend()
plt.title(f'{h = }')
plt.show();

3 Exercice 3 : choix d’un schéma

Soit le problème de Cauchy \( \begin{cases} y'(t)=\alpha (y(t)-t), &t\in[0;8]\\ y(0)=\dfrac{1}{\alpha}+\varepsilon \end{cases} \)

Pour fixer \(\alpha\), écrivez votre nom et prénom dans la cellule ci-dessous et vous obtiendrez un choix pour ce paramètre.

  1. Calculer la solution exacte.

  2. Soit \(\varepsilon=0\). Parmi les schémas suivants, choisir le plus adapté (plusieurs réponses possibles, choix à justifier):

    • EE (Euler Explicite),
    • EI (Euler Implicite),
    • EM (Euler Modifié),
    • H (Heun).

  3. Soit \(\varepsilon=0\). Approcher la solution du problème donné avec un des schémas choisis (on affichera la solution exacte et la solution approchée sur le même repère).

Code 
%reset -f
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

import sympy 
sympy.init_printing()

from IPython.display import display, Math, Markdown
prlat= lambda *args: display(Math(''.join( sympy.latex(arg) for arg in args )))

nom = "Faccanoni"
prenom = "Gloria"

L = [2,3,4,5]
idx = sum([ord(c) for c in nom+prenom])%len(L)
alpha = L[idx]
display(Math(r'\alpha='+sympy.latex(alpha)))

# alpha = 5

\(\displaystyle \alpha=2\)

  • On calcule la solution exacte (calculée en utilisant le module sympy) :
Code 
t       = sympy.Symbol('t')
epsilon = sympy.Symbol(r'\varepsilon')
y       = sympy.Function('y')

t0 = 0
y0 = sympy.Rational(1,alpha)+epsilon
display(Math('y(0)='+sympy.latex(y0)))

edo = sympy.Eq( sympy.diff(y(t),t) , alpha*(y(t)-t) )
display(edo)
solpar = sympy.dsolve(edo, y(t), ics={y(t0):y0})
display(solpar)

tt = np.linspace(0,8,101)
plt.grid()
for val in [1.e-7,0,-1.e-7]:
    sol_exacte = sympy.lambdify(t,solpar.rhs.subs(epsilon,val),'numpy')
    plt.plot(tt,sol_exacte(tt), label=f'$\epsilon={val}$')
plt.legend();

\(\displaystyle y(0)=\varepsilon + \frac{1}{2}\)

\(\displaystyle \frac{d}{d t} y{\left(t \right)} = - 2 t + 2 y{\left(t \right)}\)

\(\displaystyle y{\left(t \right)} = \varepsilon e^{2 t} + t + \frac{1}{2}\)

  • On remarque que c’est un problème de Cauchy numériquement mal posé (si \(\varepsilon=0\) la fonction est linéaire, avec \(\varepsilon>0\) elle croît comme une exponentielle et avec \(\varepsilon<0\) elle decroît comme une exponentielle): on doit donc choisir le schéma le plus précis parmi les schémas proposés, à savoir le schéma de Heun ou le schéma d’Euler Modifié qui sont d’ordre 2.
Code 
phi = lambda t,y : alpha*(y-t)

t0     = 0
tfinal = 8
y0     = 1/alpha

N = 5
tt = np.linspace(t0, tfinal, N + 1)

sol_exacte = sympy.lambdify(t,solpar.rhs.subs(epsilon,0),'numpy')
yy = sol_exacte(tt)
Code 
def EM(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for n in range(len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        uu[n+1] = uu[n]+h*phi(tt[n]+h/2,uu[n]+k1*h/2)
    return uu

def heun(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for n in range(len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        k2 = phi( tt[n+1], uu[n] + h*k1  )
        uu[n+1] = uu[n] + (k1+k2)*h/2 
    return uu


fig,ax = plt.subplots(nrows=1,ncols=2,figsize=(16,5))
ax[0].plot(tt, yy, 'r-', label='Exacte')
ax[1].plot(tt, yy, 'r-', label='Exacte')
for i,schema in enumerate([EM, heun]):
    uu = schema(phi, tt, y0)
    ax[i].plot(tt, uu, 'o', label=schema.__name__)
    ax[i].set_xlabel('$t$')
    ax[i].set_ylabel('$y$')
    ax[i].grid(True)
    ax[i].legend();

En effet, si on utilise les schémas d’ordre 1 proposés avec le même nombre de points de discrétisation on obtient

Code 
def EE(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for n in range(len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        uu[n+1] = uu[n]+h*k1
    return uu

def EI(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for n in range(len(tt)-1):
         uu[n+1] = fsolve(lambda x: -x+uu[n]+h*phi(tt[n+1],x) , uu[n] )[0]
    return uu


fig,ax = plt.subplots(nrows=1,ncols=2,figsize=(16,5))
ax[0].plot(tt, yy, 'r-', label='Exacte')
ax[1].plot(tt, yy, 'r-', label='Exacte')
for i,schema in enumerate([EE, EI]):
    uu = schema(phi, tt, y0)
    ax[i].plot(tt, uu, 'o', label=schema.__name__)
    ax[i].set_xlabel('$t$')
    ax[i].set_ylabel('$y$')
    ax[i].grid(True)
    ax[i].legend();

Remarque : dans les anciennes versions de python, les calculs avec des flottants donnaient des résultats différents.

Retour au sommet