Mise en évidence des ordres de convergence
Pour les schémas d'Euler, point au milieu, Heun et Runge-Kutta d'ordre 4.
- Numpy : vectorisation, fonction sur les vecteurs
- Matplotlib : légende, échelle logarithmique
Rappel
Considérons l'équation différentielle ordinaire
\[y'(t)=f(t,y(t)),\quad (t,y)\in U\]\[y(t_0)=y_0\]
avec \(f\), \(U\), \(y_0\), \(t\) et \(T\) comme il faut (existence, unicité, etc.). Le plus souvent \(U\) de la forme \([t_0,t_0+T]\times V\) avec \(V\) ouvert de \(\mathbb{R}^m\).
Toujours sous certaines hypothèses, nous savons que les schémas d'Euler, point au milieu, Heun et Runge-Kutta convergent avec une estimation de l'erreur vers la solution exacte.
Rappel
Les schémas (à pas fixe) consistent en une discrétisation du problème continue. Le pas est noté \(h=T/N\), la discrétisation en temps \(t_i=t_0+hi\). La suite \(y_n\) sensée approcher \(y(t_n)\) est initialisée par la donnée initiale.
Euler (explicite)
\[y_{n+1}=y_n + h f(t_n,y_{n}).\]
Point au milieu
\[y_{n+1}=y_n + h f(t_n+h/2,y_n+h/2f(t_n,y_n))\]
Méthode Heun
\[y_{n+1}=y_n+\frac{h}{2}\big(f(t_n,y_n)+f(t_{n+1},y_n+hf(t_n,y_n))\big)\]
Rappel
Méthode de Runge-Kutta d'ordre 4
\[p_{n,1}=f(t_n,y_n) \quad t_{n+1/2}=t_n+h/2\]\[y_{n,2}=y_n+hp_{n,1}/2, \quad p_{n,2}=f(t_{n+1/2},y_{n,1})\]\[y_{n,3}=y_n + hp_{n,2}/2,\quad p_{n,3}=f(t_{n+1/2},y_{n,3})\]\[y_{n,4}=y_n + hp_{n,3}, \quad p_{n,4}=f(t_{n+1},y_{n,4})\]\[y_{n+1}=y_{n}+\frac{h}{6}\big(p_{n,1}+2p_{n,2}+2p_{n,3}+p_{n,4}\big)\]
Python
Il n'y aucune difficulté à programmer les quatre méthodes. La vectorisation de Numpy permet de gérer les systèmes différentiels. Par exemple Runge-Kutta d'ordre 4 s'écrit à peu près
def myrk4_all(f,t0,T,y0,n):
# retourne toutes les valeurs
dt = T/n
Y=np.zeros((np.size(y0),n+1))
Y[:,0]=y0
tpt=np.linspace(t0,t0+T,n+1)
for i in range(n):
p1=f(tpt[i],Y[:,i])
Y2=Y[:,i]+dt/2*p1
p2=f(tpt[i]+dt/2,Y2)
Y3=Y[:,i]+dt/2*p2
p3=f(tpt[i]+dt/2,Y3)
Y4=Y[:,i]+dt*p3
p4=f(tpt[i+1],Y4)
Y[:,i+1]=Y[:,i]+dt*(p1+2*p2+2*p3+p4)/6
return Y
Erreur
Euler : \(C/n\)
Point au milieu et Heun : \(C/n^2\)
Runge-Kutta d'ordre 4 : \(C/n^4\)
Si ypt contient les approximations en la subdivision et si y est la solution exacte
np.max(np.abs(ypt-y(tpt)))
retourne le maximum de l'erreur discrète
\[\max_{0\leq n \leq N} |y_n-y(t_n)|\]
Pour la visualisation un graphique à double échelle logarithmique est l'idéal
Graphique 1
plt.xscale('log')
plt.yscale('log')
Graphique 2 (avec légende et marqueurs)
plt.plot(VN,err1,'+',label='Euler')
plt.plot(VN,err2,'*',label='Pt au milieu')
plt.legend()
Graphique 3 (avec deux légendes)
Graphique 3 (avec deux légendes)
p1, = plt.plot(VN,err1,'+')
p2, = plt.plot(VN,err2,'*')
p3, = plt.plot(VN,err3,'o')
p4, = plt.plot(VN,err4,'x')
q1, = plt.plot(VN,1/VN)
q2, = plt.plot(VN,1/(VN**2))
q3, = plt.plot(VN,1/(VN**4))
l1 = plt.legend([p1 , p2, p3, p4], ["Euler", "Pt au milieu", "Heun", "Runge Kutta 4"],loc=3)
l2 = plt.legend([q1, q2, q3],["1/N", "1/N^2", "1/N^4"],loc=1)
plt.gca().add_artist(l1)
plt.xscale('log')
plt.yscale('log')
plt.show()
Graphique 4
On choisit les entiers \(n\) de tests d'erreur avec une échelle exponentielle
In [2]: VN=(1.3)**np.arange(10,32)
In [3]: VN=VN.astype(int)
In [4]: print(VN)
[ 13 17 23 30 39 51 66 86 112 146 190 247 321 417 542
705 917 1192 1550 2015 2619 3405]
La conversion permet de faire une boucle selon VN (ce n'est pas seule solution).
Graphique 4
