70
M´ ethodes directes pour la r´ esolution des syst` emes lin´ eaires
Programme 3 - backwardcol : Substitution r´ etrograde : version orient´ ee
colonne
function [b]=backwardcol(U,b)
% BACKWARDCOL substitution r´ etrograde: version orient´ ee colonne.
% X=BACKWARDCOL(U,B) r´ esout le syst` eme triangulaire sup´ erieur
% U*X=B avec la m´ ethode de substitution r´ etrograde dans sa
% version orient´ ee colonne.
[n,m]=size(U);
if n ˜= m, error(’Seulement des syst` emes carr´ es’); end
if min(abs(diag(U))) == 0, error(’Le syst` eme est singulier’); end
for j = n:-1:2,
b(j)=b(j)/U(j,j); b(1:j-1)=b(1:j-1)-b(j)*U(1:j-1,j);
end
b(1) = b(1)/U(1,1);
return
Quand on r´ esout de grands syst` emes triangulaires, seule la partie triangulaire de la matrice doit ˆ etre stock´ ee, ce qui permet une ´ economie de m´ emoire
consid´ erable.
3.2.2 Analyse des erreurs d’arrondi
Dans l’analyse effectu´ ee jusqu’` a pr´ esent, nous n’avons pas consid´ er´ e la pr´ esence des erreurs d’arrondi. Quand on prend celles-ci en compte, les algorithmes de substitution directe et r´ etrograde ne conduisent plus aux solutions
exactes des syst` emes Lx=b et Uy=b, mais fournissent des solutions approch´ ees
x qu’on peut voir comme des solutions exactes des syst` emes perturb´ es
(L + δL) x = b, (U + δU) x = b,
o` u δL = (δl ij ) et δU = (δu ij ) sont des matrices de perturbation. En vue
d’appliquer l’estimation (3.9) ´ etablie ` a la Section 3.1.2, on doit estimer les
matrices de perturbation δL et δU en fonction des coefficients des matrices L
et U, de leur taille, et des caract´ eristiques de l’arithm´ etique `
a virgule flottante.
On peut montrer que
|δT| ≤
nu
1 − nu
|T|,
(3.21)
o` u T est ´ egal ` a L ou U et o` u u est l’unit´ e d’arrondi d´ efinie en (2.34). Clairement,
si nu < 1, en utilisant un d´ eveloppement de Taylor, il d´ ecoule de (3.21) que
|δT| ≤ nu|T| + O(u
2 ). De plus, d’apr` es (3.21) et (3.9), si nuK(T) < 1 alors
−
x
x
≤
nuK(T)
1 − nuK(T)
= nuK(T) + O(u
2 ),
(3.22)
pour les normes · · 1 , · · ∞ et la norme de Frobenius. Si la valeur de u est
assez petite (comme c’est typiquement le cas), les perturbations introduites
M´ ethodes directes pour la r´ esolution des syst` emes lin´ eaires
Programme 3 - backwardcol : Substitution r´ etrograde : version orient´ ee
colonne
function [b]=backwardcol(U,b)
% BACKWARDCOL substitution r´ etrograde: version orient´ ee colonne.
% X=BACKWARDCOL(U,B) r´ esout le syst` eme triangulaire sup´ erieur
% U*X=B avec la m´ ethode de substitution r´ etrograde dans sa
% version orient´ ee colonne.
[n,m]=size(U);
if n ˜= m, error(’Seulement des syst` emes carr´ es’); end
if min(abs(diag(U))) == 0, error(’Le syst` eme est singulier’); end
for j = n:-1:2,
b(j)=b(j)/U(j,j); b(1:j-1)=b(1:j-1)-b(j)*U(1:j-1,j);
end
b(1) = b(1)/U(1,1);
return
Quand on r´ esout de grands syst` emes triangulaires, seule la partie triangulaire de la matrice doit ˆ etre stock´ ee, ce qui permet une ´ economie de m´ emoire
consid´ erable.
3.2.2 Analyse des erreurs d’arrondi
Dans l’analyse effectu´ ee jusqu’` a pr´ esent, nous n’avons pas consid´ er´ e la pr´ esence des erreurs d’arrondi. Quand on prend celles-ci en compte, les algorithmes de substitution directe et r´ etrograde ne conduisent plus aux solutions
exactes des syst` emes Lx=b et Uy=b, mais fournissent des solutions approch´ ees
x qu’on peut voir comme des solutions exactes des syst` emes perturb´ es
(L + δL) x = b, (U + δU) x = b,
o` u δL = (δl ij ) et δU = (δu ij ) sont des matrices de perturbation. En vue
d’appliquer l’estimation (3.9) ´ etablie ` a la Section 3.1.2, on doit estimer les
matrices de perturbation δL et δU en fonction des coefficients des matrices L
et U, de leur taille, et des caract´ eristiques de l’arithm´ etique `
a virgule flottante.
On peut montrer que
|δT| ≤
nu
1 − nu
|T|,
(3.21)
o` u T est ´ egal ` a L ou U et o` u u est l’unit´ e d’arrondi d´ efinie en (2.34). Clairement,
si nu < 1, en utilisant un d´ eveloppement de Taylor, il d´ ecoule de (3.21) que
|δT| ≤ nu|T| + O(u
2 ). De plus, d’apr` es (3.21) et (3.9), si nuK(T) < 1 alors
−
x
x
≤
nuK(T)
1 − nuK(T)
= nuK(T) + O(u
2 ),
(3.22)
pour les normes · · 1 , · · ∞ et la norme de Frobenius. Si la valeur de u est
assez petite (comme c’est typiquement le cas), les perturbations introduites