ISTITUTO DI ISTRUZIONE ISTITUTO DI ISTRUZIONE SUPERIORESUPERIORE

““A. EINSTEINA. EINSTEIN””

CORSO AM08CORSO AM08Approfondimenti di Approfondimenti di

matematicamatematica

Sistemi lineari e Sistemi lineari e Metodi iterativiMetodi iterativi

Prof. Fernando Prof. Fernando DD’’AngeloAngelo

CosCos’è’è un metodo iterativo?un metodo iterativo?

I metodi iterativi consentono di I metodi iterativi consentono di ottenere una soluzione approssimata ottenere una soluzione approssimata del sistema lineare del sistema lineare AxAx = b mediante = b mediante una successione di vettori {una successione di vettori {xxkk,k ,k ≥≥ 0}, 0}, a partire da un vettore iniziale xa partire da un vettore iniziale x00

AffinchAffinchèè il metodo iterativo sia efficace il metodo iterativo sia efficace èè necessario che la successione necessario che la successione converga verso la soluzione x del converga verso la soluzione x del sistema. sistema.

La struttura generale di un La struttura generale di un metodo iterativometodo iterativo

Ax=bAx=bSi decompone la matrice A dei Si decompone la matrice A dei

coefficienticoefficientiA=MA=M--NNIn modo cheIn modo che

e M sia facilmente invertibilee M sia facilmente invertibile

0)det( ≠M

Il sistema si riscrive:Il sistema si riscrive:

a)a) Ax=bAx=bb)b) (M(M--N)N)x=bx=bc)c) Mx=Nx+bMx=Nx+bd)d) x=Mx=M--11Nx+MNx+M--11bbe)e) x=Px+qx=Px+q

Iterazione di x=Px+qIterazione di x=Px+q

xx(o)(o) vettore inizialevettore inizialexx(1)(1)=P=P xx(o)(o) +q+qxx(2)(2)=P=P xx(1)(1) +q+q……………………xx(i+1)(i+1)=P=P xx(i)(i) +q+q

Convergenza di un metodo iterativoConvergenza di un metodo iterativo

Un metodo iterativo si dice Un metodo iterativo si dice convegenteconvegente se la se la successione successione {x{x(i)(i)} } èè convergente per ogni convergente per ogni possibile scelta del vettore iniziale possibile scelta del vettore iniziale xx(o)(o)

Condizione necessaria e sufficiente per la Condizione necessaria e sufficiente per la convergenza di un metodo iterativo convergenza di un metodo iterativo èè che il che il raggio spettrale della matrice di iterazioneraggio spettrale della matrice di iterazioneP=MP=M--11N sia minore di 1N sia minore di 1

Il raggio spettrale di una matrice si calcola in Il raggio spettrale di una matrice si calcola in OctaveOctave nel modo seguente:nel modo seguente:

ρρ(P)=(P)=max(abs(eig(Pmax(abs(eig(P))))))

SoluzioneSoluzione

xx**=x=x¶¶=lim=limii-->>¶¶xx(i)(i)

I Metodi di I Metodi di JacobiJacobi e e GaussGauss--SeidelSeidel

Il metodo di Il metodo di JacobiJacobi

⎩⎨⎧

≤>−

=jijia

b ijij 0 ⎩

⎨⎧

<−≥

=jiaji

cij

ij

0

Data la matrice quadrata A definiamo Data la matrice quadrata A definiamo tre nuove matrici:tre nuove matrici:

D:matrice dei valori diagonali di AD:matrice dei valori diagonali di AB,C:matrici triangolari inferiori e B,C:matrici triangolari inferiori e superiori definite dasuperiori definite da

A=DA=D--((B+CB+C))M=DM=D

N=B+CN=B+CPPJJ=M=M--11N=DN=D--11((B+CB+C))

qqJJ =M=M--11b=Db=D--11bb

Iterazione di Iterazione di x=x= PPJJ x+x+ qqJJ

xx(o)(o) vettore inizialevettore inizialexx(1)(1)==PPJJ xx(o)(o) + + qqJJ

xx(2)(2)==PPJJ xx(1)(1) + + qqJJ

……………………xx(i+1)(i+1)==PPJJ xx(i)(i) + + qqJJ

Il metodo di Il metodo di GaussGauss--SeidelSeidel

⎩⎨⎧

≤>−

=jijia

b ijij 0 ⎩

⎨⎧

<−≥

=jiaji

cij

ij

0

Data la matrice quadrata A definiamo Data la matrice quadrata A definiamo tre nuove matrici:tre nuove matrici:

D:matrice dei valori diagonali di AD:matrice dei valori diagonali di AB,C:matrici triangolari inferiori e B,C:matrici triangolari inferiori e superiori definite dasuperiori definite da

A=DA=D--((B+CB+C)=(D)=(D--B)B)--CCM=DM=D--BBN=CN=C

PPGSGS=M=M--11N=(DN=(D--B)B)--11CCqqGSGS =M=M--11b=(Db=(D--B)B)--11bb

Iterazione di Iterazione di x=x=PPGSGSx+ x+ qqGSGS

xx(o)(o) vettore inizialevettore inizialexx(1)(1)==PPGSGS xx(o)(o) + + qqGSGS

xx(2)(2)==PPGSGS xx(1)(1) + + qqGSGS

……………………xx(i+1)(i+1)==PPGSGS xx(i)(i) + + qqGSGS

ApplicazioneApplicazione deidei metodimetodi in un in un casocasoconcretoconcreto utilizzandoutilizzando OctaveOctave

Il Il programmaprogramma principaleprincipale sisi trovatrova nelnelfile file jac_gs.mjac_gs.m;;

DalDal programmaprogramma principaleprincipale vengonovengonochiamatechiamate due functions :due functions :

(a)(a) matrice.mmatrice.m cheche èè semplicementesemplicementeunauna proceduraprocedura per per inserireinserire dadatastieratastiera matricimatrici o o vettorivettori con con comoditcomoditàà;;

(b)(b) met_jac_gs.mmet_jac_gs.m cheche èè la la proceduraproceduracheche svolgesvolge i i calcolicalcoli..

## Listato del file ## Listato del file jac_gs.mjac_gs.m## ## AuthorAuthor: Prof. Fernando : Prof. Fernando DD’’AngeloAngelo## ## CreatedCreated: 2009: 2009--0101--1111

dispdisp("Questo script consente di risolvere iterativamente un sistema ("Questo script consente di risolvere iterativamente un sistema lineare lineare Ax=bAx=b ");");dispdisp("quale metodo iterativo vuoi usare ?")("quale metodo iterativo vuoi usare ?")n_met=input("scrivi 1 per n_met=input("scrivi 1 per jacobijacobi oppure 2 per oppure 2 per gaussgauss--seidelseidel ");");dispdisp("definisci la matrice A dei coefficienti mediante procedura gui("definisci la matrice A dei coefficienti mediante procedura guidata: ");data: ");nr=inputnr=input("inserisci il numero di righe ")("inserisci il numero di righe ")nc=inputnc=input("inserisci il numero di colonne ")("inserisci il numero di colonne ")nome_elemento="a";nome_elemento="a";A=matriceA=matrice((nrnr,,ncnc,,nome_elementonome_elemento))dispdisp("definisci il vettore colonna b dei termini noti mediante proce("definisci il vettore colonna b dei termini noti mediante procedura dura guidata: ");guidata: ");b=matriceb=matrice((nrnr,1,"b"),1,"b")x0=rand(x0=rand(nrnr,1),1)tol=inputtol=input("quanto deve valere la tolleranza ? ");("quanto deve valere la tolleranza ? ");[sol,iter,errore]=met_jac_gs(A,b,x0,[sol,iter,errore]=met_jac_gs(A,b,x0,toltol,,n_metn_met););dispdisp("la soluzione e':");("la soluzione e':");solsoldispdisp("il numero di iterazioni impiegato e' :");("il numero di iterazioni impiegato e' :");iteriterdispdisp("guarda il grafico dell'errore in funzione del numero di iteraz("guarda il grafico dell'errore in funzione del numero di iterazioni");ioni");plot(1:plot(1:lengthlength(errore),(errore),erroreerrore,"1,"1--*")*")

## Listato del file ## Listato del file met_jac_gs.mmet_jac_gs.m## ## AuthorAuthor: Prof. Fernando : Prof. Fernando DD’’AngeloAngelo## ## 20092009--0101--1111

functionfunction [sol,iter,errore]=met_jac_gs(A,b,x0,[sol,iter,errore]=met_jac_gs(A,b,x0,toltol,,n_metn_met))itermax=100;itermax=100;D=diagD=diag((diagdiag(A));(A));ifif (n_met==1)(n_met==1)

M1=inv(D);M1=inv(D);N=N=--(A(A--D);D);Pj=M1*N;Pj=M1*N;dispdisp("raggio spettrale matrice iterazione di ("raggio spettrale matrice iterazione di JacobiJacobi");");maxmax((absabs((eigeig((PjPj))))))qj=M1*b;qj=M1*b;x1=Pj*x0+qj;x1=Pj*x0+qj;iter=1;iter=1;errore(iter)errore(iter)=norm=norm(x1(x1--x0);x0);whilewhile((normnorm(x1(x1--x0,2)>x0,2)>=tol=tol & iter <& iter <=itermax=itermax))

x0=x1;x0=x1;x1=Pj*x0+qj;x1=Pj*x0+qj;iter=iter+1;iter=iter+1;errore(iter)errore(iter)=norm=norm(x1(x1--x0);x0);

endend((segue nella prossima pagina)segue nella prossima pagina)

elseelseL=trilL=tril(A)(A)--D; D; U=triuU=triu(A)(A)--D; D; DI=invDI=inv((D+LD+L); ); GS=GS=--DI*U;DI*U;dispdisp("raggio spettrale matrice iterazione di ("raggio spettrale matrice iterazione di GaussSeidelGaussSeidel");");maxmax((absabs((eigeig(GS)))(GS)))b1=DI*b;b1=DI*b;x1=GS*x0+b1;x1=GS*x0+b1;iter=1;iter=1;errore(iter)errore(iter)=norm=norm(x1(x1--x0);x0);whilewhile((normnorm(x1(x1--x0,2)>x0,2)>=tol=tol & iter<& iter<=itermax=itermax))

x0=x1;x0=x1;x1=GS*x0+b1;x1=GS*x0+b1;iter=iter+1;iter=iter+1;errore(iter)errore(iter)=norm=norm(x1(x1--x0);x0);

endendendendifif iter>iter>itermaxitermax

errorerror('non converge nel numero di iterazioni prefissate')('non converge nel numero di iterazioni prefissate')elseelse

sol=x1;sol=x1;iter=iteriter=iter--1;1;

endendendend

## Listato del file ## Listato del file matrice.mmatrice.m## ## AuthorAuthor: Prof. Fernando : Prof. Fernando DD’’AngeloAngelo## 2009## 2009--0101--1111

functionfunction [y][y]=matrice=matrice((nrnr,,ncnc,,labellabel))%Funzione%Funzione per inserire una matrice in per inserire una matrice in OctaveOctave%da%da usare con la seguente sintassi:usare con la seguente sintassi:% [y]% [y]=matrice=matrice((nrnr,,ncnc))%in%in cui si cui si avra'avra' che:che:% % nrnr = numero di righe della matrice da inserire= numero di righe della matrice da inserire% % ncnc = numero di colonne della matrice da inserire= numero di colonne della matrice da inserire% % labellabel èè l'etichetta della matrice o vettore passata dal l'etichetta della matrice o vettore passata dal programma chiamanteprogramma chiamanteforfor i=1:i=1:nrnrforfor j=1:j=1:ncncy(i,j)y(i,j)=input=input((sprintfsprintf([([labellabel,"[,"[%o%o,,%o%o]="],i,j));]="],i,j));end;end;end;end;endend

Applicazione dei metodi di Applicazione dei metodi di JacobiJacobi e e GaussGauss--SeidelSeidel in un caso in un caso

concreto.concreto.

bXA =⋅

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

−−−

−−−

=

4011141111401114

A

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

=

4

3

2

1

xxxx

X

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

−=

01

21

b

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛−

=

4000040000400004

D A=DA=D--((B+CB+C))

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

−

=

0011001100000000

B

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

−−−−

=

0000100011001110

C

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

−−−

−−−

=

4011141111401114

A

Risoluzione con il metodo della matrice inversaRisoluzione con il metodo della matrice inversa

>> >> A=A=[[4 [[4 --1 1 1];[0 1 1 1];[0 --4 4 --1 1];[1 1];[--1 1 --11 4 1];[1 4 1];[1 --1 0 4]]1 0 4]]A =A =

4 4 --1 1 11 1 10 0 --4 4 --1 11 1--1 1 --11 4 14 11 1 --1 0 41 0 4

>>>>b=b=[1 2 [1 2 --1 0]'1 0]'b =b =

1122--1100

>> >> xstar=invxstar=inv(A)(A)*b*bxstarxstar ==

0.241380.24138--0.478930.47893--0.264370.26437--0.180080.18008

Adesso invece utilizziamo il metodo di Adesso invece utilizziamo il metodo di JacobiJacobi..octaveoctave--3.0.2.exe:1:C:3.0.2.exe:1:C:\\ProgrammiProgrammi\\OctaveOctave\\3.0.2_gcc3.0.2_gcc--4.3.04.3.0\\binbin> > jac_gsjac_gsQuesto script consente di risolvere iterativamente un sistema liQuesto script consente di risolvere iterativamente un sistema lineare neare Ax=bAx=bquale metodo iterativo vuoi usare ?quale metodo iterativo vuoi usare ?scrivi 1 per scrivi 1 per jacobijacobi oppure 2 per oppure 2 per gaussgauss--seidelseidel 11definisci la matrice A dei coefficienti mediante procedura guidadefinisci la matrice A dei coefficienti mediante procedura guidata:ta:inserisci il numero di righe 4inserisci il numero di righe 4nrnr = 4= 4inserisci il numero di colonne 4inserisci il numero di colonne 4ncnc = 4= 4a[1,1]=4a[1,1]=4a[1,2]=a[1,2]=--11a[1,3]=1a[1,3]=1a[1,4]=1a[1,4]=1a[2,1]=0a[2,1]=0a[2,2]=a[2,2]=--44a[2,3]=a[2,3]=--11a[2,4]=1a[2,4]=1a[3,1]=a[3,1]=--11a[3,2]=a[3,2]=--11a[3,3]=4a[3,3]=4a[3,4]=1a[3,4]=1a[4,1]=1a[4,1]=1a[4,2]=a[4,2]=--11a[4,3]=0a[4,3]=0a[4,4]=4a[4,4]=4A =A =

4 4 --1 1 11 1 10 0 --4 4 --1 11 1

--1 1 --11 4 14 11 1 --1 0 41 0 4

definisci il vettore colonna b dei termini noti mediante procedudefinisci il vettore colonna b dei termini noti mediante procedura guidata:ra guidata:b[1,1]=1b[1,1]=1b[2,1]=2b[2,1]=2b[3,1]=b[3,1]=--11b[4,1]=0b[4,1]=0b =b =

1122

--1100

x0 =x0 =

0.805020.805020.263730.263730.730390.730390.928850.92885

quanto deve valere la tolleranza ? quanto deve valere la tolleranza ? 1.e1.e--88raggio spettrale matrice iterazione di raggio spettrale matrice iterazione di JacobiJacobiansans = 0.39369= 0.39369la soluzione e':la soluzione e':sol =sol =

0.241380.24138--0.478930.47893--0.264370.26437--0.180080.18008

il numero di iterazioni impiegato e' :il numero di iterazioni impiegato e' :iter = 20iter = 20guarda il grafico dell'errore in funzione del numero di iterazioguarda il grafico dell'errore in funzione del numero di iterazionini

Andamento dellAndamento dell’’errore in funzione errore in funzione del numero di iterazioni per il del numero di iterazioni per il

metodo di metodo di JacobiJacobi

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

5 10 15 20

Adesso invece utilizziamo il metodo di Adesso invece utilizziamo il metodo di GaussGauss--SeidelSeidel..ctavectave--3.0.2.exe:1:C:3.0.2.exe:1:C:\\ProgrammiProgrammi\\OctaveOctave\\3.0.2_gcc3.0.2_gcc--4.3.04.3.0\\binbin> > jac_gsjac_gsQuesto script consente di risolvere iterativamente un sistema liQuesto script consente di risolvere iterativamente un sistema lineare neare Ax=bAx=bquale metodo iterativo vuoi usare ?quale metodo iterativo vuoi usare ?scrivi 1 per scrivi 1 per jacobijacobi oppure 2 per oppure 2 per gaussgauss--seidelseidel 22definisci la matrice A dei coefficienti mediante procedura guidadefinisci la matrice A dei coefficienti mediante procedura guidata:ta:inserisci il numero di righe 4inserisci il numero di righe 4nrnr = 4= 4inserisci il numero di colonne 4inserisci il numero di colonne 4ncnc = 4= 4a[1,1]=4a[1,1]=4a[1,2]=a[1,2]=--11a[1,3]=1a[1,3]=1a[1,4]=1a[1,4]=1a[2,1]=0a[2,1]=0a[2,2]=a[2,2]=--44a[2,3]=a[2,3]=--11a[2,4]=1a[2,4]=1a[3,1]=a[3,1]=--11a[3,2]=a[3,2]=--11a[3,3]=4a[3,3]=4a[3,4]=1a[3,4]=1a[4,1]=1a[4,1]=1a[4,2]=a[4,2]=--11a[4,3]=0a[4,3]=0a[4,4]=4a[4,4]=4A =A =

4 4 --1 1 11 1 10 0 --4 4 --1 11 1

--1 1 --11 4 14 11 1 --1 0 41 0 4

definisci il vettore colonna b dei termini noti mediante procedudefinisci il vettore colonna b dei termini noti mediante procedura guidata:ra guidata:b[1,1]=1b[1,1]=1b[2,1]=2b[2,1]=2b[3,1]=b[3,1]=--11b[4,1]=0b[4,1]=0b =b =

1122

--1100

x0 =x0 =0.3330530.3330530.5439480.5439480.0743530.0743530.0809330.080933

quanto deve valere la tolleranza ? 1.equanto deve valere la tolleranza ? 1.e--88raggio spettrale matrice iterazione di raggio spettrale matrice iterazione di GaussSeidelGaussSeidelansans = 0.17678= 0.17678la soluzione e':la soluzione e':sol =sol =

0.241380.24138--0.478930.47893--0.264370.26437--0.180080.18008

il numero di iterazioni impiegato e' :il numero di iterazioni impiegato e' :iter = iter = 1111guarda il grafico dell'errore in funzione del numero di iterazioguarda il grafico dell'errore in funzione del numero di iterazioninioctaveoctave--3.0.2.exe:2:C:3.0.2.exe:2:C:\\ProgrammiProgrammi\\OctaveOctave\\3.0.2_gcc3.0.2_gcc--4.3.04.3.0\\binbin>>

Andamento dellAndamento dell’’errore in funzione errore in funzione del numero di iterazioni per il del numero di iterazioni per il

metodo di metodo di GaussGauss--SeidelSeidel

0

0.2

0.4

0.6

0.8

1

1.2

2 4 6 8 10 12

Cosa succede applicando i metodi Cosa succede applicando i metodi iterativi al seguente sistema?iterativi al seguente sistema?(lieve modifica alla matrice A)(lieve modifica alla matrice A)

bXA =⋅'

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

−−−

−−−

=

1411141111401114

'A

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

=

4

3

2

1

xxxx

X

⎟⎟⎟⎟⎟

⎠

⎞

⎜⎜⎜⎜⎜

⎝

⎛

−=

01

21

b