248
7 Equations différentielles ordinaires
fournit l’expression du second membre de (7.69). On suppose que les
conditions initiales sont données dans le vecteur y0=[0,1,0,.8,0,1.2]
et que l’intervalle d’intégration est tspan=[0,25]. On exécute la méthode d’Euler explicite de la manière suivante :
[t , y ]= feuler ( @fvinc , tspan , y0 , nt );
(on procède de même pour les méthodes d’Euler implicite beuler et de
Crank-Nicolson cranknic), où nt est le nombre d’intervalles (de longueur
constante) utilisés pour discrétiser l’intervalle [tspan(1),tspan(2)].
Les graphiques de la Figure 7.17 montrent les trajectoires obtenues avec
10000 et 100000 noeuds de discrétisation. La solution ne semble raisonnablement précise que dans le second cas. En effet, bien qu’on ne connaisse
pas la solution exacte du problème, on peut avoir une idée de la précision
en remarquant que la solution vérifie r(y) ≡ |y
2
1 + y
2
2 + y
2
3 − 1| = 0. On
peut donc mesurer la valeur maximale du résidu r(y n ) quand n varie,
y n étant l’approximation de la solution exacte construite au temps t n .
En utilisant 10000 noeuds de discrétisation, on trouve r = 1.0578, tandis
qu’avec 100000 noeuds on a r = 0.1111, ce qui est en accord avec le résultat théorique prédisant une convergence d’ordre un pour la méthode
d’Euler explicite.
En utilisant la méthode d’Euler implicite avec 20000 pas on obtient
la solution tracée sur la Figure 7.18, tandis que la méthode de CrankNicolson (d’ordre 2) donne, avec seulement 1000 pas, la solution tracée
sur la même figure (à droite) qui est visiblement plus précise. On trouve
en effet r = 0.5816 pour la méthode d’Euler implicite et r = 0.0928 pour
la méthode de Crank-Nicolson.
A titre de comparaison, résolvons le même problème avec les méthodes adaptatives explicites de Runge-Kutta ode23 et ode45 de MATLAB. Celles-ci adaptent le pas d’intégration afin d’assurer que l’erreur
−1
−0.5
0
0.5
1
−1
−0.5
0
0.5
1
−1
−0.5
0
y
1
y 2
y
3
−1
−0.5
0
0.5
1
−1
−0.5
0
0.5
1
−1
−0.5
0
y 1
y 2
y
3
Figure 7.17. Trajectoires obtenues avec la méthode d’Euler explicite pour
h = 0.0025 (à gauche), et pour h = 0.00025 (à droite). Le point noir désigne
la donnée initiale
7 Equations différentielles ordinaires
fournit l’expression du second membre de (7.69). On suppose que les
conditions initiales sont données dans le vecteur y0=[0,1,0,.8,0,1.2]
et que l’intervalle d’intégration est tspan=[0,25]. On exécute la méthode d’Euler explicite de la manière suivante :
[t , y ]= feuler ( @fvinc , tspan , y0 , nt );
(on procède de même pour les méthodes d’Euler implicite beuler et de
Crank-Nicolson cranknic), où nt est le nombre d’intervalles (de longueur
constante) utilisés pour discrétiser l’intervalle [tspan(1),tspan(2)].
Les graphiques de la Figure 7.17 montrent les trajectoires obtenues avec
10000 et 100000 noeuds de discrétisation. La solution ne semble raisonnablement précise que dans le second cas. En effet, bien qu’on ne connaisse
pas la solution exacte du problème, on peut avoir une idée de la précision
en remarquant que la solution vérifie r(y) ≡ |y
2
1 + y
2
2 + y
2
3 − 1| = 0. On
peut donc mesurer la valeur maximale du résidu r(y n ) quand n varie,
y n étant l’approximation de la solution exacte construite au temps t n .
En utilisant 10000 noeuds de discrétisation, on trouve r = 1.0578, tandis
qu’avec 100000 noeuds on a r = 0.1111, ce qui est en accord avec le résultat théorique prédisant une convergence d’ordre un pour la méthode
d’Euler explicite.
En utilisant la méthode d’Euler implicite avec 20000 pas on obtient
la solution tracée sur la Figure 7.18, tandis que la méthode de CrankNicolson (d’ordre 2) donne, avec seulement 1000 pas, la solution tracée
sur la même figure (à droite) qui est visiblement plus précise. On trouve
en effet r = 0.5816 pour la méthode d’Euler implicite et r = 0.0928 pour
la méthode de Crank-Nicolson.
A titre de comparaison, résolvons le même problème avec les méthodes adaptatives explicites de Runge-Kutta ode23 et ode45 de MATLAB. Celles-ci adaptent le pas d’intégration afin d’assurer que l’erreur
−1
−0.5
0
0.5
1
−1
−0.5
0
0.5
1
−1
−0.5
0
y
1
y 2
y
3
−1
−0.5
0
0.5
1
−1
−0.5
0
0.5
1
−1
−0.5
0
y 1
y 2
y
3
Figure 7.17. Trajectoires obtenues avec la méthode d’Euler explicite pour
h = 0.0025 (à gauche), et pour h = 0.00025 (à droite). Le point noir désigne
la donnée initiale