import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from numpy.linalg import inv
$\newcommand{\R}{\mathbb{R}}$ $\newcommand{\N}{\mathbb{N}}$
L'objectif de ce TP est de mettre en oeuvre différents schémas numériques pour approcher la solution d'une équation différentielle ordinaire. Nous allons considérer différents schémas numériques, dont certains nous n'avons pas encore vus en cours. L'objectif est de comprendre ce que c'est une solution approchée donnée par une méthode numérique, mais aussi de comprendre des notions comme l' erreur entre la solution exacte et la solution approchée, la convergence d'une méthode numérique et l' ordre de précision d'une méthode numérique
Dans ce TP on considère des équations différentielles ordinaires, ou des systèmes d'équations différentielles ordinaires, de la forme \begin{equation} \label{EDO} \left\{\begin{aligned} y'(t) & = f(t,y(t)), \\ y(t_0) & = y_0, \end{aligned} \right. \tag{EDO} \end{equation}
où $f$ : $I \times \R^n \longrightarrow \R^n$ est une fonction dans les conditions du théorème de Cauchy-Lipschitz, avec $I\subseteq\R$ un intervalle ouvert de $\R$ et où $(t_0,y_0)\in I\times\R^n$ est donné. Le problème de Cauchy (EDO) admet alors une unique solution maximale définie dans un intervalle ouvert $J\subseteq I$.
On souhaite calculer une solution approchée de ce problème dans un intervalle de la forme $[t_0,t_f]=[t_0,t_0+T]\subseteq J$, avec $T>0$. Pour ce faire, on se donne $N\in\N$ et on considère une discrétisation de $[t_0,t_0+T]$ de pas $h=\frac TN,$ donnée par les points $$ t_0< t_1< \cdots <t_N=t_0+T,\ \ \ t_n=nh,\ n=0,\dots,N, $$ et l'on utilise un schéma numérique afin de calculer une approximation $y^n$ de $y(t_n)$ pour $n=0,\dots,N$.
Les schémas explicites qui suivent sont fréquemment utilisés pour discrétiser le problème (EDO). On pose $y^0=y_0,\,$ et, pour $n=0,\dots,N-1,$
Euler explicite : $y^{n+1}=y^{n}+h f(t_n,y^{n})$ où $f^{n}=f(t_n,y^{n})$ ;
Heun : $y^{n+1}=y^{n}+\frac{h}{2}\big(f(t_n,y^{n})+f\big(t_n+h,y^{n}+h f(t_n,y^{n})\big)\big)$ où $p_1=f(t_n,y^{n})$ et $p_2=f(t_n+h,y^{n}+h p_1)$ ;
Q1) Écrire une fonction python de la forme euler_exp(fct, t0, tf, y0, h)
prenant en argument fct
la fonction $f(t,y)$ de (EDO), les extrémités t0
et tf
de l'intervalle de temps $[t_0,t_f]=[t_0,t_0+T],$ la donnée initiale y0
et le pas de temps h
. Cette fonction devra retourner deux tableaux :
[t_0,\, t_1,\, ...,\, t_N]
, tableau unidimensionnel de taille $(N+1)\times 1$ représentant la subdivision de l'intervalle $[t_0,t_f]$ de pas $h$ considérée,
[y^0,\, y^1,\, ...,\, y^N]
tableau de taille $(N+1)\times n$ représentant la solution approchée aux instants $t_n,\ n=0,\dots,N.$
#Il est plus simple de définir d'abord une version scalaire des schémas (c'est à dire un schéma pour approcher une EDO scalaire) :
def euler_exp_s(f,t0,tf,y0,h):
T=np.arange(t0,tf+h,h) # T est le vecteurs (t_0,...,t_NN) avec t_n=t_0+n*h
N=T.size # Pour connaître la taille de T. !!! ATTENTION !!! N=taille de T correspond à NN+1 ci-dessus, le vecteur (t_0,...,t_NN) a NN+1 points. Ce qu'on appelle ici N n'est pas ce qu'on appelle N dans les notes de cours : N(ici) est égal à N(cours)+1
y=np.zeros(N) # y est le vecteur (y_0,...,y_N) qui a donc la même taille de T, i. e. N
y[0]=y0
for n in range(N-1):
y[n+1]=y[n]+h*f(y[n],T[n])
return T,y
# La version qui suit est la version vectorielle ; elle peut être aussi utilisée dans le cas scalaire si on définit y0 comme un numpy array ( par exemple, si y0=1 il faut alors définir y0=np.array([1.]) ). On va par la suite toujours utiliser cette version vectorielle
def euler_exp(f,t0,tf,y0,h):
T=np.arange(t0,tf+h,h)
N=T.size
y=np.zeros((N,y0.size))
y[0,:]=y0
for n in range(N-1):
y[n+1,:]=y[n,:]+h*f(y[n,:],T[n])
return T,y
Q2) Testez d'abord votre fonction euler_exp
sur un exemple simple d'une EDO scalaire. Vous pouvez par exemple considérer le problème
\begin{equation}\label{edoexp}
\begin{cases}
y'(t)=y(t),\\
y(0)=1,
\end{cases}
\end{equation}
dont la solution est $y(t)=e^t.$
Q2.a) Tracer sur une même fenêtre graphique :
La solution exacte sur l'intervalle $[0,1]$ discrétisé avec un pas de temps de $10^{-4}$ (autrement dit, en remplaçant l'intervalle $[0,1]$ par l'ensemble discret correspondant aux points d'une sub-division de $[0,1]$ de pas $10^{-4}$) ;
Les solutions approchées sur le même intervalle obtenues avec la méthode d'Euler avec des pas de temps $h$ égaux à $1/4,\ 1/10,\ 1/100.$ Observer que lorsque le pas diminue, la solution approchée est plus proche de la solution exacte.
def f(y,t):
return y
def u(t):
return np.exp(t)
t0=0
tf=1
# représentation de la solution exacte avec un pas de 10^(-4). REMARQUE IMPORTANTE : pour représenter la solution exacte, qui est une fonction définie dans un continuum de points, il faut utiliser beaucoup de points, sinon la courbe que l'on obtient ne sera pas proche de son graphique.
pas=10**(-4)
vect=np.arange(t0,tf+pas,pas)
plt.figure(1)
plt.plot(vect,u(vect),label='sol. exacte')
plt.xlabel('t')
plt.ylabel('y')
# représentation des solutions approchées obtenues avec le schéma d'EE avec des pas 1/4, 1/10 et 1/100
y0=np.array([1.])
# pour stocker les valeurs approchées de e^1 (Question 2b) et l'erreur |e-y_app(1)|
valfin=[np.exp(1)]
errvalfin=[]
H=[1./4, 1./10, 1./100, 1./1000]
plt.figure(1)
for h in H:
T,y=euler_exp(f,t0,tf,y0,h)
plt.plot(T,y,label=('sol app h=%s' %h))
plt.legend()
plt.title('y\'=y : solution exacte et solutions approchées par Euler exp')
Text(0.5,1,"y'=y : solution exacte et solutions approchées par Euler exp")
Q2.b) Donner les valeurs approchées de $y(1)=e$ obtenues avec $h=1/10,\ 1/100,\ 1/1000$, et donner les erreurs entre la valeur exacte de $e$ et les valeurs approchées obtenues avec ces pas de temps. Observer la décroissance de l'erreur lorsque le pas de temps diminue.
Q2.c) Tracer, sur une autre figure, la différence en valeur absolue entre la solution exacte et la solution approchée, en fonction du temps (i.e. $|y(t_n)-y^n|$ en fonction de $t_n$), en utilisant des discrétisations de $[0,1]$ avec des pas de temps $h$ égaux à $1/10,\ 1/100,\ 1/1000.$ Vous devez alors tracer trois courbes et observer que lorsque le pas diminue, cette différence s'approche de plus en plus de 0.
plt.figure(2)
for h in H:
T,y=euler_exp(f,t0,tf,y0,h)
plt.plot(T,np.abs(u(T)-y[:,0]),label=('h=%s' %h))
valfin=np.append(valfin,y[-1]) # la valeur approchée de e=y(1) est la dernière position du vecteur y
errvalfin=np.append(errvalfin,np.exp(1)-y[-1])
plt.legend()
plt.title('|sol exacte - sol approchée| pour différentes valeurs du pas h')
# Affichage des valeurs approchées de exp(1) (solution exacte en t=1)
print('\n\n')
print('e^1 et valeurs approchées avec h=1/10, 1/100 et 1/1000')
print(valfin[0], valfin[1], valfin[2], valfin[3])
print('\n\n')
# Affichage de l'erreur entre val approché et valeur exacte en t=1
print('erreur en t=1 : |y(tN) -y^N| avec h=1/10, 1/100 et 1/1000')
print(abs(errvalfin[0]), abs(errvalfin[1]), abs(errvalfin[2]))
print('\n\n')
e^1 et valeurs approchées avec h=1/10, 1/100 et 1/1000 2.718281828459045 2.44140625 2.5937424601 2.704813829421526 erreur en t=1 : |y(tN) -y^N| avec h=1/10, 1/100 et 1/1000 0.2768755784590451 0.124539368359045 0.013467999037519274
Q3) Testez votre fonction euler_exp
sur le modèle logistique
$$
(P_1)\ \ \ \
\begin{cases}
y'(t)=c y (1 - \frac{y}{b}),\\
y(0)=a,
\end{cases}
$$
dont la solution exacte est
$$
y(t) = \frac{b}{1 + \frac{b-a}{a} e^{-ct}},
$$
avec les données $c=1$, $b=2$, $a=0.1$, dans l'intervalle $[0,15]$, avec un pas $h=0.2$. Tracer sur la même fenêtre la solution exacte et la solution approchée, obtenue avec le pas $h$ et avec un pas égal à $\frac h2$.
def f1(y,t,b=2.,c=1.):
return c*y*(1.-y/b)
# Solution exacte du modèle logistique (P1) (exo 2, Q3)
def u1(t,a=0.1,b=2.,c=1.):
return b/(1+(np.exp(-c*t))*(b-a)/a)
t0=0
tf=15
y0=np.array([0.1])
h=0.2
# solutions approchées par E. exp, avec pas h et avec pas h/2
T1,y1=euler_exp(f1,t0,tf,y0,h)
T2,y2=euler_exp(f1,t0,tf,y0,h/2)
# pour représenter la solution exacte
pas=10**(-4)
vect=np.arange(t0,tf+pas,pas)
plt.figure(3)
plt.plot(vect,u1(vect),label='sol exacte')
plt.plot(T1, y1,label='sol app h=0.2')
plt.plot(T2, y2,label='sol app h/2')
plt.legend()
plt.xlabel('t')
plt.ylabel('y')
plt.title('Eq. logistique : solution exacte et solutions approchées par E. explicite')
Text(0.5,1,'Eq. logistique : solution exacte et solutions approchées par E. explicite')
Q4) Testez ensuite votre fonction dans le cas vectoriel $n>1$ ($n=2$) sur le problème $$ (P_2)\ \ \ \ \begin{cases} y''(t) = -y(t) + \cos(t) \\ y(0) = 5,\ \ y'(0) = 1, \end{cases} $$ dont la solution exacte est
$$ y(t) = \frac{1}{2} \sin(t) t + 5 \cos(t) + \sin(t), $$dans l'intervalle $[0,15]$, avec un pas $h=0.2$.\ Pour ce faire, il faut écrire l'équation d'ordre 2 de $(P_2)$ sous la forme d'un système de deux équations d'ordre 1 dans les nouvelles variables $u(t)=y(t)$ et $v(t)=y'(t)$. On se ramènera alors à la résolution d'une équation de la forme
$$ X'=F(t,X),\ \ \ \textrm{avec}\ X=(u,v)=(y,y')^T. $$Représenter à nouveau la solution exacte et la solution approchée dans une même fenêtre graphique.
Remarque : La solution $y$ de $(P_2)$ correspond à la première composante du vecteur $X$ ci-dessus. Votre fonction euler_ex
! retournera dans ce cas, si vous avez respecté la structure conseillée, un tableau de taille $(N+1)\times2$, $N$ étant le nombre de points de la discrétisation. Ce tableau donne les valeurs approchées de $X$ au points de la discrétisation considérée, la première colonne correspondant à la première composante de $X$, la seconde à la seconde composante de $X$. La solution approchée de $(P_2)$ que l'on cherche correspond alors à la première colonne de ce tableau.
# EDO d'ordre 2 (P2) y''=-y+cos(t) sous forme d'un système de 2 éq. d'ordre 1 (exo 2, Q4)
def f2(y,t): # Écriture vectorielle u=y, v=y'
return np.array([y[1],-y[0]+np.cos(t)])
def u2(t):
return np.sin(t)*t/2+5*np.cos(t)+np.sin(t)
# Le système s'écrit U'=F(t,U), avec U=(y,y') et F(t,U)=(y',-y+cos(t)) : VOIR AU DÉBUT COMMENT DÉFINIR LA FONCTION F.
t0=0
tf=15
h=0.2
y0=np.array([5.,1.])
# pour représenter la solution exacte
vect=np.linspace(t0,tf+pas,1000)
T,y=euler_exp(f2,t0,tf,y0,h)
plt.figure(4)
plt.plot(vect,u2(vect),label='sol exacte')
plt.plot(T,y[:,0],label='sol app h=0.2') # on représente y, donc uniquement la première colonne de la solution (y[:,0]=y,y[:,1]=y')
plt.title('y\'\'=-y+cos(t) : solution exacte et solutions approchées par E. explicite')
plt.legend()
plt.xlabel('t')
plt.ylabel('y')
# COMMENTAIRE :
# L'approximation n'est pas très bonne (autrement dit l'erreur est grande), on voit dans la figure, cela veut dire qu'on n'a pas utilisé un pas suffisamment petit. Pour voir si le schéma est bien convergent, on divise le pas par 10, puis par 100
h=0.02
T,y=euler_exp(f2,t0,tf,y0,h)
plt.figure(5)
plt.plot(vect,u2(vect),label='sol exacte')
plt.plot(T,y[:,0],label='sol app h=0.02')
h=0.002
T,y=euler_exp(f2,t0,tf,y0,h)
plt.plot(T,y[:,0],label='sol app h=0.002')
plt.legend()
plt.xlabel('t')
plt.ylabel('y')
plt.title('y\'\'=-y+cos(t) : solution exacte et solutions approchées par E. explicite')
# COMMENTAIRE :
# Les résultats sont rassurants, le schéma semble bien converger.
Text(0.5,1,"y''=-y+cos(t) : solution exacte et solutions approchées par E. explicite")
Dans cet exercice on essaie de comprendre ce que c'est l'erreur associé à un schéma numérique et ce que c'est un schéma numérique convergent.
On se donne un problème de Cauchy de la forme (EDO) et une méthode numérique pour approcher la solution de (EDO) dans un intervalle de la forme $[t_0,t_0+T]$. Pour un certain pas de temps $h$ fixé (ou pour un certain nombre de points de la discrétisation $N$ fixé), l'erreur globale entre la solution approchée associée à une discrétisation de pas $h$ de l'intervalle $[t_0,t_0+T]$ et la solution exacte est donnée par :
\begin{equation*} E_h = \max_{k=0,\cdots,N}( |y(t_k)-y^{k}|). \end{equation*}Ci dessus, $y(t_k)$ est la solution exacte à l'instant $t_k$ et $y^{k}$ la valeur approchée de $y(t^k),$ donnée par le schéma numérique.
Remarque 1 : l'erreur globale $E$ dépend du pas $h$, ou, de manière équivalente, du nombre de points $N$ de la discrétisation. Pour des valeurs de $N$ différentes, ou de manière équivalente pour des valeurs de $h$ différentes, les discrétisations de l'intervalle $[t_0,t_0+T]$ sont différentes (elles ont un nombre de points différent), et les solutions approchées sont différentes. On s'attend à que, lorsque $N$ augmente, ou de manière équivalente lorsque $h$ diminue, l'erreur $E_h$ diminue, puisque l'on considère dans ce cas de plus en plus de points dans la discrétisation de l'intervalle $[t_0,t_0+T]$ et les approximations que l'on a faites pour construire le schéma numérique deviennent alors de plus en plus précises.
Remarque 2 : Parfois dans la littérature on ne spécifie pas la dépendance de $E$ par rapport à $h$, mais on doit toujours garder à l'esprit cette dépendance et que l'erreur est donc une fonction de $h$
Considérons le problème $$ (P_3)\ \ \ \ \begin{cases} y'(t)=\frac{\cos(t)-y(t)}{1+t},\\ y(0)=-\frac14, \end{cases} $$ dont la solution exacte est $$ y(t) = \frac{\sin(t)-1/4}{1 + t}. $$
Q1) Calculez les solutions approchées de $(P_3)$ obtenues avec le schéma d'Euler explicite, avec $h=1/2^s$ pour $s = 1,2,...,8$ ; représentez dans la même figure la différence en valeur absolue entre la solution exacte et la solution approchée, en fonction du temps, pour chaque valeur de $h$.
# EDO (P3) y'=(cos(t)-y)/(1+t) (exo 3)
def f3(y,t):
return (np.cos(t)-y)/(1+t)
# Solution exacte du problème (P3) y'=(cos(t)-y)/(1+t), y(0)=-1/4 (exo 3)
def u3(t):
return (np.sin(t)-1./4)/(1+t)
t0=0
tf=20
y0=np.array([-1./4])
h=0.5
vect=np.linspace(t0,tf+pas,1000)
# listes ou on va garder les pas de temps et les erreurs globales correspondantes a chaque pas de temps
vec_h=[]
vec_err=[]
# figure où l'on va afficher la différence en valeur absolue entre sol exacte et sol approchée
plt.figure(6)
for k in range(8):
T,y=euler_exp(f3,t0,tf,y0,h)
plt.plot(T,np.abs(u3(T)-y[:,0]),label=h) # !!!! ATTENTION !!!! il faut faire np.abs(u3(T)-y) si on utilise la version scalaire d'Euler explicite
E=max(abs(y[:,0]-u3(T))) # !!!! ATTENTION !!!! E=max(abs(y-u2(T))) si on utilise la version scalaire d'Euler explicite
vec_h.append(h)
vec_err.append(E)
print('erreur global pour h=',h,' est égale à ',E)
print('\n')
h=h/2
plt.legend()
plt.title('y\'=(cos(t)-y)/(1+t) : |sol. exacte - sol. approchée| pour différentes valeurs du pas h')
erreur global pour h= 0.5 est égale à 0.24679202822617596 erreur global pour h= 0.25 est égale à 0.11040093896816727 erreur global pour h= 0.125 est égale à 0.05280190074448804 erreur global pour h= 0.0625 est égale à 0.02585727613212707 erreur global pour h= 0.03125 est égale à 0.012797869725736155 erreur global pour h= 0.015625 est égale à 0.006366858911482087 erreur global pour h= 0.0078125 est égale à 0.003175485838083414 erreur global pour h= 0.00390625 est égale à 0.0015857718084631434
Text(0.5,1,"y'=(cos(t)-y)/(1+t) : |sol. exacte - sol. approchée| pour différentes valeurs du pas h")
Q2) Calculer, pour chaque valeur de $h=1/2^s,\ s = 1,2,...,8,$ l'erreur globale $E_h$ correspondante. Vérifiez que $E_h$ diminue et s'approche de $0$ lorsque $h$ diminue.
On dira qu'une méthode numérique converge si $\displaystyle{\lim_{h\to0}E_h}=0$.
Représentez ensuite, en échelle logarithmique, l'erreur en fonction du pas de temps $h$, autrement dit, représentez $\log(E_h)$ en fonction de $\log(h)$. Vous devez obtenir des points qui sont à peu près alignés sur une droite de pente 1. Vérifiez graphiquement que c'est le cas, en estimant la pente de la droite passant au plus prêt des points (ou en représentant une droite de pente $1$ qui passe par un des points et en vérifiant que tous les points sont à peu près sur cette droite).
Ceci signifie que $\log(E_h) \sim C+\log(h)$ et donc que $E_h\sim \widetilde{C}h$, pour certaines constantes $C$ et $\widetilde{C}$. On dit que la méthode d'Euler explicite est d' ordre 1 : c'est l'ordre de la puissance de $h$ dans cette relation. On a donc que l'erreur globale $E_h$ tend vers 0 comme $h$. L'ordre d'une méthode donne une indication sur sa vitesse de convergence. Une méthode d'ordre $p$ est une méthode dont l'erreur globale tend vers $0$ comme $h^p$. Donc plus l'ordre est élevé, plus la méthode converge plus vite
Remarque : pour étudier numériquement l'ordre d'une méthode, on utilise souvent l'échelle logarithmique pour tracer l'erreur en fonction du pas de discrétisation $h$. La pente de la droite obtenue donne l'ordre $p$ de la méthode : si $E_h \sim Ch^p$ alors $\log(E_h)\sim \log(C) + p\log(h)$.
plt.figure(7)
plt.plot(vec_h,vec_err,'x-',label='erreur')
plt.plot(vec_err,vec_err,'x-',label='droite de pente 1') #Ceci affiche une droite de pente 1.
plt.xscale('log')
plt.yscale('log')
plt.legend()
plt.title('Erreur global en fonction de h, en échelle logarithmique')
plt.xlabel('log(h)')
plt.ylabel('log($E_h$)')
Text(0,0.5,'log($E_h$)')
Q1) À l'image de ce que vous avez fait pour le schéma d'Euler explicite, écrire une fonction python de la forme Heun(fct, t0, tf, y0, h)
prenant les mêmes arguments que la fonction euler_exp
et retournant les mêmes tableaux.
def Heun(f,t0,tf,y0,h):
T=np.arange(t0,tf+h,h)
N=T.size
y=np.zeros((N,y0.size))
y[0,:]=y0
for n in range(N-1):
p1=f(y[n,:],T[n])
p2=f(y[n,:]+h*p1,T[n+1])
y[n+1,:]=y[n,:]+h*(p1+p2)/2
return T,y
Q2) Pour le problème $(P_1),$ comparez les solutions approchées obtenues avec la méthode d'Euler explicite et avec la méthode de Heun avec le même pas de temps. Pour ce faire, représenter dans la même fenêtre graphique la solution exacte et les solutions approchées obtenues avec chacune des méthodes. Quelle méthode semble être plus précise ? Pour le confirmer, représenter dans une autre figure la différence, en valeur absolue, entre la solution exacte et la solution approchée, pour les deux méthodes.
t0=0
tf=15
y0=np.array([0.1])
h=0.5
Teul,yeul=euler_exp(f1,t0,tf,y0,h)
Theun,yheun=Heun(f1,t0,tf,y0,h)
vect=np.linspace(t0,tf,1000)
plt.figure(8)
plt.plot(vect,u1(vect),label='sol. exacte')
plt.plot(Teul, yeul,label='sol. Euler exp.')
plt.plot(Theun, yheun,label='sol. Heun')
plt.legend()
plt.xlabel('t')
plt.ylabel('y')
plt.title('(P1) : sol. exacte et sol. approchées par E. exp et Heun')
plt.figure(9)
plt.plot(Teul, abs(yeul[:,0]-u1(Teul)),label='erreur Eul. exp.')
plt.plot(Theun, abs(yheun[:,0]-u1(Theun)),label='erreur Heun')
plt.legend()
plt.title('|sol. exacte - sol. approchee| pour Eul. exp et Heun')
# COMMENTAIRE :
# Pour la même valeur du pas, la solution approchée obtenue avec le schéma d'Heun est plus précise (plus proche de la solution exacte) que celle obtenue par le schéma d'Euler explicite.
Text(0.5,1,'|sol. exacte - sol. approchee| pour Eul. exp et Heun')
Q3) Reprenez l'exercice 2 pour la méthode de Heun. Vérifiez graphiquement que cette méthode est d'ordre $2,$ en représentant l'erreur globale en fonction du pas de discrétisation en échelle logarithmique et en estimant la pente de la droite passant au plus prêt des points correspondant à cette représentation ou en représentant une droite de pente $2$ qui passe par un des points.
t0=0
tf=20
y0=np.array([-1./4])
h=0.5
# listes ou on va garder les pas de temps et les erreurs globales correspondantes a chaque pas de temps
vec_h=[]
vec_err=[]
for k in range(8):
T,y=Heun(f3,t0,tf,y0,h)
E=max(abs(y[:,0]-u3(T)))
vec_h.append(h)
vec_err.append(E)
h=h/2
plt.figure(10)
plt.plot(vec_h,vec_err,'x-',label='erreur')
plt.plot(vec_h,vec_h,'x-',label='droite de pente 1') #Ceci affiche une droite de pente 1.
plt.plot(vec_h,np.array(vec_h)**2,'x-',label='droite de pente 2') #Ceci affiche une droite de pente 2.
plt.xscale('log')
plt.yscale('log')
plt.legend()
plt.title('sch. Heun : erreur global en fonction de h, en échelle logarithmique')
plt.xlabel('log(h)')
plt.ylabel('log($E_h$)')
# COMMENTAIRE :
#En échelle logarithmique, les points correspondant à l'erreur globale du schéma d'Heun en fonction de h semblent bien à
#peu prés alignés sur une droite de pente 2 puisque l'on a représenté une droite de pente 2 à coté qui semble parallèle à
#la ligne qui contient les points. Cela montre que le schéma d'Heun est d'ordre 2. Il est donc plus précis que le schéma d'Euler
#(car l'erreur tend vers 0 comme h^2, alors que pour Euler l'erreur tend vers 0 comme h).
Text(0,0.5,'log($E_h$)')