Calcolo numerico - metodo di Gauss Seidel e formula di Simpson

METODI ITERATIVI PER LA RISOLUZIONE DI UN SISTEMA LINEARE

Molti problemi conducono alla risoluzione di un sistema lineare Ax=b di dimensioni molto grandi con

matrice A sparsa, ovvero con pochi elementi non nulli. Se ad un tale sistema si applica un metodo diretto,

le matrici dei sistemi intermedi o di arrivo possono diventare matrici dense, cioè con un elevato numero di

elementi non nulli.

Sorgono così seri problemi di costo computazionale e di ingombro di memoria. In questi casi può giovare il

ricorso ai metodi iterativi, in cui ogni iterazione richiede il prodotto di una matrice P per un vettore. Poiché

la densità di P è paragonabile a quella di A, se questa è una matrice sparsa, ogni iterazione comporta una

mole di calcoli relativamente modesta ed un ingombro di memoria limitato.

Dato il sistema lineare Ax=b

con detA 0, scomposta la matrice A in A=M-N,

con M non degenere, il problema dato può equivalentemente porsi nella forma

-1 1

x = M Nx M b

la cui soluzione è un vettore, detto punto unito, che soddisfa ovviamente l’uguaglianza.

Il problema posto nella forma precedente suggerisce il procedimento iterativo:

( i 1

) 1 ( i ) 1 ;

x M Nx M b

Gael MASULLO – Università di Pisa – Mat. 450569 – Ing. Delle Telecomunicazioni Pagina 1

Metodi Iterativi per la risoluzione di un sistema lineare: il Metodo di Gauss - Seidel

Integrazione Numerica: Le formule di quadratura: La formula di Simpson

(0)

cioè, assegnato un vettore iniziale x , si determina

( 1

) 1 ( 0 ) 1

x M Nx M b

e quindi ( 2 ) 1 ( 1

) 1

x M Nx M b

e così via.

Si osserva che se la successione così costruita è convergente ad un vettore x allora x è sicuramente

soluzione. (0)

( i )

È interessante però sapere a priori quando la successione , generata a partire da un x arbitrario,

{

x }

risulta convergente. -1 -1

A tale scopo posto P=M N e q= M b il problema

x=Px+q

ha soluzione unica se e solo se (P)<1

ossia se e solo se P è una matrice convergente.

In pratica il processo iterativo, supposto convergente, si arresta quando

Gael MASULLO – Università di Pisa – Mat. 450569 – Ing. Delle Telecomunicazioni Pagina 2

Metodi Iterativi per la risoluzione di un sistema lineare: il Metodo di Gauss - Seidel

Integrazione Numerica: Le formule di quadratura: La formula di Simpson

( i 1

) ( i ) .

x x

Tuttavia, non è detto che la soluzione sia approssimata con la precisione .

Si conclude che, l’errore relativo sulla soluzione, può essere grande se, il numero di condizionamento di A è

grande nonostante l’errore relativo sul residuo, sia piccolo.

Metodo iterativo di Gauss-Seidel

La matrice A si scompone nella forma specifica A=D-L-U dove D=diag{ }, , L={-a se i>j, 0 se i j

a , a ,..., a ij

11 22 nn

} e U={-a se i<j, 0 se i j}.

ij

In questo caso M=D-L, N=U e il metodo in forma matriciale si scrive

( i 1

) 1 ( i ) 1

x ( D L

) Ux ( D L

) b

mentre esplicitamente in termini di componenti diviene

( i 1

) 1 ( i 1

) ( i )

x D ( Lx Ux b

)

ossia j 1 n

1

( i 1

) ( i 1

) (

ki ) , j=1,...,n .

x ( b a x a x )

j j jk k jk

a k 1 k j 1

jj

Diremo poi che una matrice è a predominanza diagonale forte se

n

a a , i 1

,..., n ,

ii ij

j 1

j i

Gael MASULLO – Università di Pisa – Mat. 450569 – Ing. Delle Telecomunicazioni Pagina 3

Metodi Iterativi per la risoluzione di un sistema lineare: il Metodo di Gauss - Seidel

Integrazione Numerica: Le formule di quadratura: La formula di Simpson

All’atto pratico può accadere che il metodo di Jacobi sia convergente e quello di Gauss-Seidel no e

viceversa. Entrambi i metodi convergono se la matrice A è a predominanza diagonale forte per righe o

per colonne.

Se alcuni elementi della diagonale di A sono nulli è sempre possibile determinare un sistema equivalente

permutando le righe e/o le colonne di A in modo che la nuova matrice abbia gli elementi diagonali non

nulli.

Nel caso la matrice A sia hermitiana il metodo di Gauss-Seidel è convergente se e solo se A è definita

positiva.

Nelle pagine successive sono riportati degli esempi di applicazione.

Gael MASULLO – Università di Pisa – Mat. 450569 – Ing. Delle Telecomunicazioni Pagina 4

Metodi Iterativi per la risoluzione di un sistema lineare: il Metodo di Gauss - Seidel

Integrazione Numerica: Le formule di quadratura: La formula di Simpson

MATRICE HERMITIANA/SIMMETRICA DEFINITA POSITIVA:

(Raggio spettrale matrice di iterazione=0,837<1)

Gael MASULLO – Università di Pisa – Mat. 450569 – Ing. Delle Telecomunicazioni Pagina 5

Metodi Iterativi per la risoluzione di un sistema lineare: il Metodo di Gauss - Seidel

Integrazione Numerica: Le formule di quadratura: La formula di Simpson

Gael MASULLO – Università di Pisa – Mat. 450569 – Ing. Delle Telecomunicazioni Pagina 6

Metodi Iterativi per la risoluzione di un sistema lineare: il Metodo di Gauss - Seidel

Integrazione Numerica: Le formule di quadratura: La formula di Simpson

MATRICE A PREDOMINANZA DIAGONALE FORTE:

(Raggio spettrale matrice di iterazione=0,50<1)

Gael MASULLO – Università di Pisa – Mat. 450569 – Ing. Delle Telecomunicazioni Pagina 7

Metodi Iterativi per la risoluzione di un sistema lineare: il Metodo di Gauss - Seidel

Integrazione Numerica: Le formule di quadratura: La formula di Simpson

Gael MASULLO – Università di Pisa – Mat. 450569 – Ing. Delle Telecomunicazioni Pagina 8

Metodi Iterativi per la risoluzione di un sistema lineare: il Metodo di Gauss - Seidel

Integrazione Numerica: Le formule di quadratura: La formula di Simpson

DIVERGENZA METODO DI GAUSS SEIDEL:

(Raggio spettrale matrice di iterazione=2>1)

Gael MASULLO – Università di Pisa – Mat. 450569 – Ing. Delle Telecomunicazioni Pagina 9

Metodi Iterativi per la risoluzione di un sistema lineare: il Metodo di Gauss - Seidel

Integrazione Numerica: Le formule di quadratura: La formula di Simpson

CODICE SORGENTE

% Progetto Risoluzione Sistemi Lineari con metodo di Gauss-Seidel

%

% Programma elaborato da

%

% Gael MASULLO

% Mat 450569

% Università di Pisa

% Corso di Calcolo Numerico - Prof.ssa L. ACETO

%

% Informazioni sul programma

%

% Scopo: Risolve un sistema lineare tramite il metodo di Gauss-Seidel, con

% la fattorizzazione A=D-L-U. La funzione prosegue nel processo iterativo

% fino a quando:

% 1) Il margine di errore raggiunge la tolleranza richiesta;

% 2) Si sono effettuate nmax valutazioni;

%

% Parametri:

%

% Input: A = Matrice del sistema Ax=b

% b = Termine noto del sistema Ax=b

% x0 = vettore di partenza per la risoluzione del sistema

% toll = Tolleranza richiesta

% nmax = Massimo numero consentito di iterazioni

%

% Output: x= soluzione approssimata del sistema

% k= numero di iterazioni utilizzate

% D,L,U= Matrici della fattorizzazione di A

%

function [x,k,D,L,U]=gs(A,b,x0,toll,nmax)

%%%%%%%%%%%%%%%%%%% Definizione D L U %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Parte diagonale

D=diag(diag(A));

% Parte triangolare inferiore (esclusa la diagonale principale)

L=-tril(A,-1);

% Parte triangolare superiore (esclusa la diagonale principale)

U=-triu(A,1);

%%%%%%%%%%%%%%%%%%%%%%%% inizializzazione %%%%%%%%%%%%%%%%%%%%%%%%%%

P=D-L;

x=x0;

rs=max(eig(P\U))

Gael MASULLO – Università di Pisa – Mat. 450569 – Ing. Delle Telecomunicazioni Pagina 10

Metodi Iterativi per la risoluzione di un sistema lineare: il Metodo di Gauss - Seidel

Integrazione Numerica: Le formule di quadratura: La formula di Simpson

%%%%%%%%%%%%%%%%%%% iterazioni di Gauss-Seidel %%%%%%%%%%%%%%%%%%%%%%

for k=1:nmax

z=(P\U)*x+P\b;

if (norm(z-x)<=toll)

break

end

x=z;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (k==nmax)

warning('Gauss-Seidel non converge')

else

fprintf('\nGauss-Seidel e'' converge in %d iterazioni\n',k);

end

return

**********************************************************

Gael MASULLO – Università di Pisa – Mat. 450569 – Ing. Delle Telecomunicazioni Pagina 11

Metodi Iterativi per la risoluzione di un sistema lineare: il Metodo di Gauss - Seidel

Integrazione Numerica: Le formule di quadratura: La formula di Simpson

INTEGRAZIONE NUMERICA E FORMULE DI QUADRATURA

Il problema è di approssimare l’integrale b f ( x ) dx

a

in cui f(x) è una funzione continua in [a,b].

Una formula d’integrazione numerica è della forma

n

I w f ( x )

n 1 i i

i 0

 

dove, sono detti nodi e pesi.

x , x , , x w , w , , w

0 1 n 0 1 n i 

Questa formula permette di calcolare il grado di precisione p se per f(x) = , il resto è nullo,

x i 0

, 1

, , p

p 1

mentre per f(x) = è diverso da zero.

x

La differenza fra l’approssimazione dell’integrale ed il risultato di una qualsiasi formula di integrazione

r I

n 1 n 1

è detta resto e rappresenta l’errore analitico.

Gael MASULLO – Università di Pisa – Mat. 450569 – Ing. Delle Telecomunicazioni Pagina 12

Metodi Iterativi per la risoluzione di un sistema lineare: il Metodo di Gauss - Seidel

Integrazione Numerica: Le formule di quadratura: La formula di Simpson

Formula di Simpson

La regola di Simpson prevede la suddivisione dell'intervallo di integrazione in sottointervalli e la

sostituzione in questi sottointervalli della funzione integranda mediante archi di parabola, cioè mediante

polinomi quadratici. La formula di Simpson è un caso particolare delle formule di Newton-Cotes,

nell’ipotesi che n=2: x x

N 1 N 1

b a j j 1

( N )

J f ( x ) 2 f ( x ) 4 f f ( x )

3 0 j N

6 N 2

j 1 j 0

con l’errore analitico o resto espresso dalla formula 5

( b a )

( N ) ( 4 ) .

R f ( )

3 4

2880 N

Per calcolare l’errore, nel sorgente in Matlab, è stato utilizzato un altro metodo, ovvero quello di

Richardson. ( 2 N ) ( N )

J J

( 2 N )

J s

2 1

L’algoritmo presentato esegue i seguenti passi:

1. Inizializzazione delle variabili.

2. Valutazione della funzione negli estremi a et b, e salvataggio di tali valori che verranno poi

riutilizzati in seguito.

3. Suddivisione dell’intervallo di integrazione (a,b) in due sottointervalli e calcolo della prima stima

dell’integrale usando la formula di Simpson semplice su 3 punti.

4. Inizio del ciclo while che da