9.7 Chapitre 7
343
Solution 7.17 On utilise les instructions suivantes :
t0 =0; y0 =0; f = inline ( ’ cos (2* y) ’, ’t ’, ’y ’);
y = inline ( ’ 0.5* asin (( exp (4* t ) -1)./( exp (4*t )+1)) ’ ,’t ’ );
T =1; N =2; for k =1:10;
[ tt ,u ]= rk2(f ,[t0 ,T ] ,y0 ,N );
e (k )= abs( u( end) - feval (y , tt ( end ))); N =2* N; end
p = log ( abs (e (1: end -1)./ e (2: end )))/ log (2); p (1:2: end )
2.4733
2.1223
2.0298
2.0074
2.0018
Comme prévu, on trouve que l’ordre de convergence de la méthode est 2. Le
coût du calcul est pourtant comparable à celui de la méthode d’Euler explicite,
qui n’est que d’ordre 1.
Solution 7.18 L’équation différentielle du second ordre de cet exercice est
équivalente au système du premier ordre suivant
x
(t) = z(t), z
(t) = −5z(t) − 6x(t),
avec x(0) = 1, z(0) = 0. On utilise la méthode de Heun :
t0 =0; y0 =[1 0]; T =5;
[t , u ]= rk2 ( @fspring ,[ t0 ,T ] ,y0 ,N );
où N est le nombre de noeuds et fspring.m est la fonction suivante :
function fn = fspring (t , y)
b =5;
k =6;
[n , m ]= size ( y );
fn= zeros (n , m );
fn (1)= y (2);
fn (2)= -b* y (2) -k * y (1);
Sur la Figure 9.10, on représente les deux composantes de la solution, calculées
avec N=20,40 et on les compare à la solution exacte x(t) = 3e
−2t − 2e
−3t et à
sa dérivée première.
0
1
2
3
4
5
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
0
1
2
3
4
5
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
Figure 9.10. Approximations de x(t) (trait plein) et x
(t) (trait discontinu)
calculées avec N=20 (à gauche) et N=40 (à droite). Les petits cercles (resp.
carrés) représentent la solution exacte x(t) (resp. x
(t))