5.9 Méthodes itératives
163
d’itérations nmax et la tolérance tol pour le test d’arrêt. On stoppe les
itérations si le rapport entre la norme euclidienne du résidu courant et
celle du résidu initial est inférieur ou égal à la tolérance tol (pour une
justification de ce critère d’arrêt, voir la Section 5.12).
Programme 5.2. itermeth : méthode itérative générale
function [x , iter ]= itermeth (A ,b , x0 , nmax , tol , P)
% ITERMETH
Méthode i t é rat ive générale
% X = ITERMETH (A ,B ,X0 , NMAX , TOL ,P ) tente de résoudre le
% système d ’ é q u a tions l i n é aires A *X =B d ’ inconnue X .
% La matrice A , de taille NxN , doit etre i n v e rsi ble et
% le second membre B doit être de longueur N.
% P = ’J ’ s é l e ctio nne la methode de Jacobi , P = ’G ’ celle
% de Gauss - Seidel . Autrement , P est une matrice N x N
% qui joue le rôle de p r é c o nd itio nne ur dans la methode
% de R i c h ard son d y n a mique.
% Les i t é r ations s ’ arrêtent quand le rapport entre la
% norme du k - ème residu et celle du résidu initial est
% i n f ér ieure ou égale à TOL , le nombre d ’ i t é r ati ons
% e f f ec tuées est alors renvoyé dans ITER .
% NMAX est le nombre maximum d ’ i t é r ation s. Si P
% n ’ est pas défini , c ’ est la méthode du Gradient à
% pas optimal qui est utilisée
[n , n ]= size ( A );
if nargin == 6
if ischar ( P )==1
if P == ’J ’
L= diag ( diag( A )); U= eye( n );
beta =1; alpha =1;
elseif P == ’G ’
L= tril (A ); U = eye( n );
beta =1; alpha =1;
end
else
[L ,U ]= lu(P );
beta = 0;
end
else
L = eye (n ); U = L;
beta = 0;
end
iter = 0;
x = x0;
r = b - A * x0;
r0 = norm (r );
err = norm (r );
while err > tol & iter < nmax
iter = iter + 1;
z = L\r; z = U\z;
if beta == 0
alpha = z ’* r /(z ’* A* z );
end
x = x + alpha * z;
r = b - A * x;
err = norm (r ) / r0;
end
return