Calcolo Numerico II – Elaborato

Descrizione del metodo di Jacobi

Il metodo di Jacobi è un metodo iterativo utilizzato per risolvere sistemi lineari di grandi dimensioni con una matrice dei coefficienti sparsa. La condizione sufficiente per la convergenza di questo metodo è che la matrice dei coefficienti sia diagonale strettamente dominante.

La routine genera un errore se si verificano le seguenti condizioni:

• Il metodo non converge dopo un numero massimo di iterazioni.
• Il numero di parametri è insufficiente.
• Il parametro non è una matrice quadrata sparsa di numeri.
• Il parametro non è un vettore di numeri.
• Le dimensioni dei parametri non sono compatibili. La lunghezza di A deve essere pari al numero di colonne di b.
• Il parametro non è uno scalare reale positivo.
• Il parametro non è uno scalare reale maggiore o uguale a TOL.
• La diagonale della matrice contiene elementi nulli.

La funzione gseidel può essere utilizzata come alternativa al metodo di Jacobi.

Per ulteriori informazioni, consulta la documentazione relativa alla funzione gseidel.

con l'algoritmo di Gauss Seidel.

Sintassi

x = gseidel(A, b)
x = gseidel(A, b, TOL)
x = gseidel(A, b, TOL, MAXITER)
[x, NITER] = gseidel(A, b)
[x, NITER] = gseidel(A, b, TOL)
[x, NITER] = gseidel(A, b, TOL, MAXITER)

Descrizione

La funzione gseidel risolve il sistema lineare con matrice quadrata A di dimensione m-by-m sparsa dei coefficienti e vettore b dei termini noti di lunghezza m. Il risultato, se il metodo converge, è il vettore soluzione calcolato con un'accuratezza dipendente dal parametro TOL utilizzato e dal condizionamento della matrice (vedi la funzione condest).
La funzione gseidel, oltre a restituire il vettore soluzione calcolata del sistema, restituisce anche il numero di iterazioni effettuate per il calcolo della soluzione.

NITER, risolve il sistema con una tolleranza assegnata, se TOL è il numero di cifre significative esatte k desiderate nel risultato (a meno delle cifre perse per il condizionamento della matrice A).

A = gallery('poisson',10);
b = rand(size(A,1),1);
[x,NITER] = gseidel(A,b,10^-7,1000);

parametro non è uno scalare reale positivo.

TOLSe il parametro non è uno scalare reale maggiore o uguale ad .MAXITER 1

Se la diagonale della matrice contiene elementi nulli.

Algoritmo usa il metodo iterativo di Gauss Seidel per risolvere sistemi lineari di grandi dimensioni con matrice dei coefficienti sparsa.

Condizione sufficiente per la convergenza del metodo di Gauss Seidel è che la matrice dei coefficienti del sistema sia a diagonale strettamente dominante o sia simmetrica e definita positiva.

Vedi anche| |jacvet condest eps

Esempi di accuratezza, funzioni JACVET e GSEIDEL

Esempio 1 Esempio 2 Esempio 3

Matrice 100x100 Matrice 100x100 Matrice 1024x1024

Tabella Riassuntiva

Ben-condizionata Mal-condizionata Ben-condizionata

100*100 100*100 1024*1024

Dimensione matrice Poisson A Bande Poisson

Descrizione matrice 69.8634 5.1882e+009 640.3620

Condizionamento matrice 1.0000e-007 eps=2.2204e-016 1.0000e-007

Tolleranza richiesta 6.9863e-006 6.4036e-0051.1520e-006

Tol. * Cond.

311 2324101NITER per JACVET 165 2 1240NITER per GSEIDEL 2.2869e-006 2.1925e-0053.2839e-008ERR-REL per JACVET 1.1024e-006 1.0874e-0053.2839e-008ERR-REL per GSEIDEL 0.0312 sec 0.0312 sec 0.1560 secTempo per JACVET 0.2184 sec 0.0312 sec 60.1228 secTempo per GSEIDELEsempio 1: Con matrice 100x100 ben-condizionata>> n = 100;>> A = gallery('poisson',sqrt(n));>> spy(A); 0102030405060708090100 0 20 40 60 80 100nz = 460>> condest(A)ans =69.8634>> xsol = ones(n,1);>> b = A*xsol;>> [x, NITER] = jacvet(A, b, 10^(-7), 1000);>> norm(x-xsol)/norm(x)ans =2.2869e-006>> NITERNITER =311 1>> [x, NITER] = gseidel(A, b, 10^(-7), 1000);>> norm(x-xsol)/norm(x)ans =1.1024e-006>> NITERNITER =165>> t=cputime; [x, NITER] = jacvet(A,b,10^(-7),1000); temp=cputime-t;>> temptemp =0.0312>> t=cputime; [x, NITER] = gseidel(A,b,10^(-7),1000); temp=cputime-t;>> temptemp =0.2184>>Esempio 2: Con matrice 100x100

mal-condizionata>> n = 100; d = ones(100,1);>> A = spdiags([-d -0.9*d 2.1*d],[-1 0 -1],n,n);>> spy(A); 
0102030405060708090100 
0 20 40 60 80 100
nz = 199>> condest(A)
ans = 5.1882e+009>> det(A)
ans = 2.6561e-005>> xsol = rand(n,1);>> b = A*xsol;>> [x, NITER] = jacvet(A, b, eps, 1000); 
2>> norm(x-xsol)/norm(x)
ans = 3.2839e-008>> NITER
NITER = 101>> [x, NITER] = gseidel(A, b, eps, 1000);>> norm(x-xsol)/norm(x)
ans = 3.2839e-008>> NITER
NITER = 2>> t=cputime; [x, NITER] = jacvet(A,b,eps,1000); temp=cputime-t;>> temptemp = 0.0312>> t=cputime; [x, NITER] = gseidel(A,b,eps,1000); temp=cputime-t;>> temptemp = 0.0312>>Esempio 3: Con matrice 1024x1024 ben-condizionata>> n = 1024;>> A = gallery('poisson',sqrt(n));>> spy(A); 
010020030040050060070080090010000 
200 400 600 800 1000
nz = 4992>> condest(A)
ans = 640.3620 3>> xsol = ones(n,1);>> b = A*xsol;>> [x,

NITER] = jacvet(A, b, 10^(-7), 5000);
>> norm(x-xsol)/norm(x)
ans = 2.1925e-005
>> NITER
NITER = 2324
>> [x, NITER] = gseidel(A, b, 10^(-7), 2000);
>> norm(x-xsol)/norm(x)
ans = 1.0874e-005
>> NITER
NITER = 1240
>> t=cputime; [x, NITER] = jacvet(A,b,10^(-7),5000); temp=cputime-t;
>> temp
temp = 0.1560
>> t=cputime; [x, NITER] = gseidel(A,b,10^(-7),2000); temp=cputime-t;
>> temp
temp = 60.1228
>> 4Diagrammi efficienza temporale, funzioni JACVET e GSEIDEL-14NITER ERROR RELx 10
5000 JACVET JACVET
4000 GSEIDEL GSEIDEL
13000
2000 0.5
1000 0
0 200 400 600 0 200 400 600
TIME(sec) COND
50 500
JACVET
40 400
GSEIDEL
30 300
20 200
10 100
0 0 200 400 600 0 200 400 600

I diagrammi riportati sono generati utilizzando come matrice dei coefficienti la matrice di
riportate sull'asse dell'ascisse.
POISSON, a vari ordini di grandezza Processore Intel® Core™ 2 Duo T7300 2GHz,
Il calcolatore utilizzato per i test è un laptop HP con
Memoria centrale 2048 MB,

Sistema Operativo Windows 7 Professional, Matlab 7.9.0(R2009b)

Di seguito è riportato lo script per generare i grafici:

data = zeros(2,7);
dataj = zeros(3,7);
datag = zeros(3,7);
ord = round(logspace(0.5,1.4,7));
for index = 1:7,
    A = gallery('poisson',ord(index));
    dim = ord(index).^2;
    xsol = rand(dim,1);
    b = A*xsol;
    data(:,index) = [dim; condest(A)];
    start = cputime;
    [jx jNITER] = jacvet(A,b,eps,5000);
    jtime = cputime - start;
    dataj(:,index) = [jNITER; norm(jx-xsol)/norm(jx); jtime];
    start = cputime;
    [gx gNITER] = gseidel(A,b,eps,5000);
    gtime = cputime - start;
    datag(:,index) = [gNITER; norm(gx-xsol)/norm(gx); gtime];
    disp(['passo: ', num2str(index),' di 7']);
end;
figure;
subplot(2,2,1); plot(data(1,:), dataj(1,:), 'b-', data(1,:), datag(1,:), 'r-');
axis([0,650,0,5000]); legend('JACVET','GSEIDEL','Location','NorthWest');
title('NITER');
subplot(2,2,2);plot(data(1,:), dataj(2,:), 'b-', data(1:),

datag(2,:), 'r-');
axis([0,650,0,1.5*10^(-14)]);
legend('JACVET','GSEIDEL','Location','NorthWest');
title('ERROR REL');
subplot(2,2,3);
plot(data(1,:), dataj(3,:), 'b-', data(1,:), datag(3,:), 'r-');
axis([0,650,0,50]);
legend('JACVET','GSEIDEL','Location','NorthWest');
title('TIME(sec)');
subplot(2,2,4);
plot(data(1,:), data(2,:), 'g-');
axis([0,650,0,500]);
title('COND');
5Test di Robustezza, funzione JACVETVerifica della funzione jacvet nei confronti dei parametri di ingresso:Script del test Output ottenuto EsitoOmissione parametri di input>> x = jacvet; ??? Error usingI parametri A e b sono obbligatori.TJ1 Positivo>> x = jacvet(sprand(2,2,0.2)); ??? Error usingI parametri A e b sono obbligatori.TJ2 PositivoCon matrice A non numerica, sparsa, quadrata>> x = jacvet(‘ciao’,[1;2]); ??? Error usingParametro A

invalido: deve contenere solo numeri.

>> x = jacvet(rand(3,3),[ 1 2 3]); Parametro A invalido: deve essere una matrice sparsa.

>> x=jacvet(sprand(5,3,0.08),[ 1 2 3]); Parametro A invalido: deve essere una matrice quadrata.

Con matrice A avente elementi NaN o Inf

>> A = sprand(30,30,0.08); ??? Error using

>> A(3) = nan; A(6) = nan; Parametro A invalido: non deve contenere valori NaN.

>> x = jacvet(A,1:30); Parametro A invalido: non deve contenere valori NaN.

>> A = sprand(30,30,0.08); ??? Error using

>> A(6) = inf; Parametro A invalido: non deve contenere valori Inf.

>> x = jacvet(A,1: