Calcolo Numerico II – Elaborato

Elaborato 1 calcolo numerico 2, A.A. 2009/2010

Cozzolino Davide 885000466

De Simone Luigi 885000538

jacvet: Risoluzione di sistemi lineari sparsi

Il file jacvet jacvet risolve un sistema lineare sparso con l'algoritmo di Jacobi.

Sintassi

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

Descrizione

La funzione jacvet(A, b) risolve il sistema lineare A*x=b dove A è una matrice quadrata sparsa dei coefficienti e b è un vettore 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 A (vedi la funzione condest), oltre a restituire il vettore soluzione, restituisce anche il numero di iterazioni effettuate per il calcolo della soluzione.

jacvet(A, b, TOL) risolve il sistema con una tolleranza assegnata. Se k è il numero di cifre significative esatte decimali desiderate nel risultato (a meno delle cifre perse per il condizionamento della matrice A), allora TOL deve essere pari a 10^-k. Se la tolleranza viene omessa, si utilizza la massima precisione possibile per il sistema aritmetico, pari all'epsilon macchina (vedi la funzione eps).

jacvet(A, b, TOL, MAXITER) risolve il sistema limitando a MAXITER il numero di iterazioni da effettuare. Se MAXITER è omesso, sarà impostato al valore 1000.

Esempio

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

Indicatori d'errore

La routine genera warning:

• Se il parametro TOL è minore dell'epsilon macchina. In tal caso viene impostato pari all'epsilon macchina.
• Se il metodo non converge dopo MAXITER iterazioni.

La routine genera errore:

• Se il numero di parametri è insufficiente.
• Se il parametro A non è una matrice quadrata sparsa di numeri.
• Se il parametro b non è un vettore di numeri.
• Se le dimensioni dei parametri A e b non sono compatibili. La lunghezza di b deve essere pari al numero di colonne di A.
• Se il parametro TOL non è uno scalare reale positivo.
• Se il parametro MAXITER non è uno scalare reale maggiore o uguale a 1.
• Se la diagonale della matrice A contiene elementi nulli.

Algoritmo

jacvet usa il metodo iterativo di Jacobi per risolvere sistemi lineari di grandi dimensioni con matrice dei coefficienti sparsa. Condizione sufficiente per la convergenza del metodo di Jacobi è che la matrice dei coefficienti del sistema sia a diagonale strettamente dominante.

Vedi anche: gseidel, condest, eps

gseidel: Risoluzione di sistemi lineari sparsi

Il file gseidel risolve un sistema lineare sparso 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(A, b) risolve il sistema lineare A*x=b dove A è una matrice quadrata sparsa dei coefficienti e b è un vettore 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 A (vedi la funzione condest), oltre a restituire il vettore soluzione, restituisce anche il numero di iterazioni effettuate per il calcolo della soluzione.

gseidel(A, b, TOL) risolve il sistema con una tolleranza assegnata. Se k è il numero di cifre significative esatte decimali desiderate nel risultato (a meno delle cifre perse per il condizionamento della matrice A), allora TOL deve essere pari a 10^-k. Se la tolleranza viene omessa, si utilizza la massima precisione possibile per il sistema aritmetico, pari all'epsilon macchina (vedi la funzione eps).

gseidel(A, b, TOL, MAXITER) risolve il sistema limitando a MAXITER il numero di iterazioni da effettuare. Se MAXITER è omesso, sarà impostato al valore 1000.

Esempio

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

Indicatori d'errore

La routine genera warning:

• Se il parametro TOL è minore dell'epsilon macchina. In tal caso viene impostato pari all'epsilon macchina.
• Se il metodo non converge dopo MAXITER iterazioni.

La routine genera errore:

• Se il numero di parametri è insufficiente.
• Se il parametro A non è una matrice quadrata sparsa di numeri.
• Se il parametro b non è un vettore di numeri.
• Se le dimensioni dei parametri A e b non sono compatibili. La lunghezza di b deve essere pari al numero di colonne di A.
• Se il parametro TOL non è uno scalare reale positivo.
• Se il parametro MAXITER non è uno scalare reale maggiore o uguale a 1.
• Se la diagonale della matrice A contiene elementi nulli.

Algoritmo

gseidel 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: Con matrice 100x100 ben-condizionata

n = 100;
A = gallery('poisson',sqrt(n));
spy(A);
condest(A)
xsol = ones(n,1);
b = A*xsol;
[x, NITER] = jacvet(A, b, 10^(-7), 1000);
norm(x-xsol)/norm(x)
NITER
[x, NITER] = gseidel(A, b, 10^(-7), 1000);
norm(x-xsol)/norm(x)
NITER
t=cputime; [x, NITER] = jacvet(A,b,10^(-7),1000); temp=cputime-t;
temp
t=cputime; [x, NITER] = gseidel(A,b,10^(-7),1000); temp=cputime-t;
temp

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);