Gauss-Seidel
Scopo
Risolve un sistema lineare del tipo Ax = b attraverso l’algoritmo iterativo di
Gauss-Seidel, con A matrice sparsa e b vettore dei termini noti.
Specifiche e Comportamento
Per la realizzazione di tale algoritmo abbiamo utilizzato un’unica funzione,
gausiedel.m che utilizza le seguenti sintassi possibili.
Sintassi
[x] = gausiedel(A, b)
[x, NITER] = gausiedel(A, b)
[x] = gausiedel(A, b, TOL)
[x] = gausiedel(A, b, TOL, MAXITER)
[x, NITER] = gausiedel(A, b, TOL)
[x, NITER]= gausiedel(A, b, TOL, MAXITER)
Parametri di input
1. A (obbligatorio) : matrice di numeri con m righe ed n colonne.
2. b (obbligatorio) : vettore dei termini noti.
3. k (facoltativo) : numero di cifre significative utilizzate per il calcolo della
soluzione (10^-k), se omesso la tolleranza è pari all’epsilon macchina.
Bisogna indicare il numero di cifre significative che si vogliono ottenere nel
calcolo (k positivo).
4. MAXITER (facoltativo) : numero massimo di iterazioni
all’algoritmo, se omesso MAXITER si imposta al valore 1000.
consentito
Parametri di output
1. x (obbligatorio) : vettore soluzione del sistema Ax = b.
2. NITER (facoltativo) : numero di iterazioni effettuate dall’algoritmo per il
calcolo della soluzione ottenuta.
Descrizione
Analizziamo ora una descrizione formale dei vari modi in cui può essere utilizzata
tale funzione dall’utente :
[x] = gausiedel(A, b)
risolve il sistema lineare Ax = b, con A matrice sparsa, b vettore dei termini
noti, restituendo il vettore soluzione. Il metodo iterativo terminerà con una
tolleranza pari all'epsilon macchina o ad un numero massimo
di iterazioni pari
a 1000.
[x, NITER] = gausiedel(A, b)
risolve il sistema lineare Ax = b, restituendo il vettore x, soluzione calcolata del
sistema, ed il numero di iterazioni effettuate per il calcolo di tale soluzione. Il
metodo iterativo terminerà con una tolleranza pari all'epsilon macchina o ad un
numero massimo di iterazioni pari a 1000.
[x, NITER] = gausiedel(S, b,k)
risolve il sistema con una tolleranza TOL=10^-k (dove k è il numero di cifre
significative esatte). Il metodo iterativo terminerà se si giunge alla tolleranza
imposta oppure se si giunge ad un numero massimo di iterazioni pari a 1000.
[x, NITER] = gausiedel(S, b, k, MAXITER)
risolve il sistema limitando a MAXITER il numero di iterazioni da effettuare. l
metodo iterativo terminerà se si giunge alla tolleranza imposta oppure se si
giunge ad un numero massimo di iterazioni pari a MAXITER.
ESEMPIO D'USO
Analizziamo ora il risultato di un’esempio d’uso basato su tali parametri di ingresso : matrice
1500x1500 con elementi pari a 9 sulla diagonale principale e 1 sulle due diagonali adiacenti.
Richiamando il seguente script file:
format long
n=1500;
A=9*eye(n)+diag(ones(n-1,1),1)+diag(ones(n-1,1),-1);
A = sparse(A);
spy(A)
indice_condizionamento=condest(A)
xsol = rand(length(A),1);
b = A*xsol;
[x, NITER] = gausiedel(A, b);
errore=norm(x-xsol)/norm(x)
numero_iterazioni=NITER
Otteniamo:
indice_condizionamento =
1.571428571428571
errore =
1.094258536091669e-016
numero_iterazioni =
18
Richiamando sulla stessa matrice A,
[x, NITER] = gausiedel(A, b,10);
Otteniamo :
indice_condizionamento =
1.571428571428571
errore =
3.563365286007379e-012
numero_iterazioni =
12
Richiamando la function gausiedel nella seguente maniera:
[x, NITER] = gausiedel(A, b,eps,10);
Otteniamo:
indice_condizionamento =
1.571428571428571
errore =
2.595270241660143e-010
numero_iterazioni =
10
Richiamando la function nella seguente maniera:
[x, NITER] = gausiedel(A, b,50,10);
Otteniamo:
indice_condizionamento =
1.571428571428571
Warning: la tolleranza indicata è troppo elevata, è
posta pari a eps
> In gausiedel at 47
In gauss_errore at 20
errore =
2.507364375013947e-010
numero_iterazioni =
10
Matrice 100x100 risultante dalla discretizzazione attraverso l'equazione di
Poisson.
format long
A = gallery('poisson', 10); %di ordine n^2
spy(A)
indice_condizionamento=condest(A)
xsol = rand(length(A),1);
b = A*xsol;
[x, NITER] = gausiedel(A,b);
errore=norm(x-xsol)/norm(x)
numero_iterazioni=NITER
Otteniamo:
indice_condizionamento =
69.863370896510318
errore =
2.177630186169446e-015
numero_iterazioni =
405
Situazioni di errore
Le varie situazioni di errore che possono verificarsi sono le seguenti :
Omissione del vettore dei termini noti:
gausiedel(A);
Otteniamo:
termini noti
??? Error using ==> gausiedel at 21
Bisogna inserire la matrice A, e
Error in ==> gauss_errore at 21
gausiedel(A);
Omissione parametri di output:
il
vettore
dei
gausiedel(A, b)
Otteniamo:
??? Error using ==> gausiedel at 63
Indicare almeno il termine che risulterà essere il
vettore soluzione
Error in ==> gauss_errore at 21
gausiedel(A,b);
Inserimento di una matrice non sparsa:
format long
n=1500;
A=rand(n);
spy(A)
indice_condizionamento=condest(A)
xsol = rand(length(A),1);
b = A*xsol;
[x, NITER] = gausiedel(A, b,16);
errore=norm(x-xsol)/norm(x)
numero_iterazioni=NITER
Otteniamo:
??? Error using ==> gausiedel at 24
matrice non sparsa
Error in ==> gauss_errore at 20
[x, NITER] = gausiedel(A, b,50,10);
MAXITER è pari a NaN
[x, NITER] = gausiedel(A, b,15,NaN);
Otteniamo:
indice_condizionamento =
1.571428571428571
Warning: il numero di iterazioni deve
numero positivo,sarà
impostato a 1000
> In gausiedel at 50
In gauss_errore at 23
errore =
1.450996268306523e-016
numero_iterazioni =
17
essere
un
La diagonale contiene elementi nulli:
A=9*eye(n)+diag(ones(n-1,1),1)+diag(ones(n-1,1),-1);
A = sparse(A);
A(1)=0;
Otteniamo:
??? Error using ==> gausiedel at 27
sulla diagonale ci sono elementi nulli, per applicare
Gauss-Siedel
trasformare la matrice
Error in ==> gauss_errore at 20
[x, NITER] = gausiedel(A, b,12,10);
La matrice sparsa non è quadrata:
A=rand(n,10); %per una matrice non quadrata
A=sparse(A);
[x, NITER] = gausiedel(A, b,16);
Otteniamo:
??? Error using ==> gausiedel at 30
La matrice A deve essere quadrata
Error in ==> gauss_errore at 20
[x, NITER] = gausiedel(A, b,12,10);
La variabile b non è un vettore:
b=rand(3,2);
[x, NITER] = gausiedel(A, b,16);
Otteniamo:
??? Error using ==> gausiedel at 33
b deve essere il vettore dei termini noti
Error in ==> gauss_errore at 20
[x, NITER] = gausiedel(A, b,12,10);
La lunghezza di S non corrisponde con la lunghezza di b:
b=rand(3,1);
[x, NITER] = gausiedel(A, b,16);
Otteniamo:
??? Error using ==> gausiedel at 36
La matrice A e il vettore b hanno diversa lunghezza
Error in ==> gauss_errore at 20
[x, NITER] = gausiedel(A, b,12,10);
Codice e documentazione interna
Gauss-Seidel
%Metodo di gauss-siedel per la risoluzione di un sistema lineare
Ax=b
%con A sparsa
%[x] = gausiedel(A, b)
%parametri di input della funzione
%A (obbligatorio) – matrice sparsa di numeri con m righe ed n
colonne.
%b (obbligatorio) – vettore dei termini noti.
%function [x,Niter]=gausiedel(A,b,k,MaxIter)
%k (facoltativo) – numero di cifre significative utilizzate per
il calcolo
%della soluzione
%MAXITER (facoltativo) – numero massimo di iterazioni consentito
all’algoritmo.
%Parametri di output:
%x (obbligatorio) – vettore soluzione del sistema Ax = b.
%NITER
(facoltativo)
–
numero
di
iterazioni
effettuate
dall’algoritmo per il calcolo della soluzione ottenuta.
function [x,Niter]=gausiedel(A,b,k,MaxIter)
%controllo sul numero di parametri in ingresso
if nargin < 2
error('Bisogna inserire la matrice A, e il vettore dei
termini noti');
end
if(issparse(A)==0)
error('matrice non sparsa')
end
if(find(diag(A)==0))
error('sulla diagonale ci sono elementi nulli, per applicare
Gauss-Siedel trasformare la matrice')
end
if size(A,1) ~= size(A,2)
%numero di righe diverso numero di
colonne
error('La matrice A deve essere quadrata')
end
if(isvector(b)==0)
error('b deve essere il vettore dei termini noti')
end
if length(A) ~= length(b)
error('La matrice A e il vettore b hanno diversa lunghezza')
end
if nargin<4
MaxIter=1000;
if nargin<3
TOL=eps;
end
end
if nargin>2
TOL=max(10^(-k),eps);
if(TOL==eps)
warning('la tolleranza indicata è troppo elevata, è
posta pari a eps')
end
if (MaxIter<1 |isnan(MaxIter))
warning('il numero di iterazioni deve essere un numero
positivo,sarà impostato a 1000')
MaxIter=1000;
end
end
x=zeros(length(A), 1);
xk=x;
Niter=0;
%controllo sul numero di parametri in uscita
if nargout < 1
error('Indicare almeno il termine che risulterà essere il
vettore soluzione')
end
while(((norm(x-xk)>TOL*norm(xk))| Niter==0) && (Niter<MaxIter))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%
% Niter==0 è fondamentale affinchè non si esca dal ciclo alla
prima
% iterazione, infatti inizialmente norm(xk) è = 0, norm(x-xk)
= 0 si
% esce.
%if(Niter==0)
%
norm(x-xk)
% norm(xk)
%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%
xk=x;
for i=1:length(x)
x(i)=(-(A(i,1:i-1)*x(1:i1)+A(i,i+1:length(x))*x(i+1:length(x)))+b(i))/A(i,i);
end
Niter=Niter+1;
end