Corso di Calcolo Numerico - old · maggior parte dei suoi elementi `e non nullo. Deﬁnizione Una...

Facolta di Ingegneria - Laurea Triennale in Ingegneria Meccanica

Corso di Calcolo Numerico

Dott.ssa M.C. De Bonis

Universita degli Studi della Basilicata, Potenza

Facolta di Ingegneria

Corso di Laurea in Ingegneria Meccanica

Corso di Calcolo Numerico - Dott.ssa M.C. De Bonis


Risoluzione di Sistemi LineariVogliamo trovare un’approssimazione delle soluzioni disistemi di n equazioni lineari in n incognite:

a1,1x1 +a1,2x2 + · · · +a1,nxn = b1

a2,1x1 +a2,2x2 + · · · +a2,nxn = b2

......

......

...

an−1,1x1 +an−1,2x2 + · · · +an−1,nxn = bn−1

an,1x1 +an,2x2 + · · · +an,nxn = bn

(1)

Ponendo

A =

a1,1 a1,2 · · · a1,n−1 a1,n

a2,1 a2,2 · · · an,n−1 a2,n

......

. . ....

...

an−1,1 an−1,2 · · · an−1,n−1 an−1,n

an,1 an,2 · · · an,n−1 an,n

b =

b1

b2

...

bn−1

bn

e x = [x1x2 · · ·xn−1xn]T , il sistema (1) diventa:

Ax = b, A ∈ Rn×n, b, x ∈ Rn



E noto che il sistema (1) ammette un’unica soluzione se esolo se la matrice A e non singolare, cioe

det(A) 6= 0.

In questo caso possiamo scrivere

x = A−1b.

Il metodo di risoluzione che certamente conoscete e laRegola di Cramer:

xi =det(Ai)det(A)

, i = 1, . . . , n,

dove Ai denota la matrice ottenuta da A sostituendo lacolonna i−esima con il vettore b.

I determinanti coinvolti potrebbero essere calcolatiutilizzando la Regola di Laplace scritta nella forma

det(A) =n∑

j=1

(−1)j+1a1,j det(A1,j),

dove A1,j rappresenta la matrice di ordine n− 1 ottenutada A eliminando la prima e la j−esima colonna.



La Regola di Cramer non e ammissibile come metodonumerico in quanto il suo costo computazionale edell’ordine di [(n + 1)(n− 1)n]!

Supponendo di poter effettuare ogni operazionemoltiplicativa in 7 · 10−8 secondi, per risolvere un sistema20× 20 occorrono circa 2 · 106 anni!!

Si potrebbe pensare di calcolare A−1 e risolvere il sistemautilizzando la formula

x = A−1b.

Tra tutti i metodi che si conoscono per il calcolo dellamatrice inversa, quello meno “costoso” ha un costocomputazionale dell’ordine di n3. Inoltre e necessarioeseguire il prodotto A−1b che ha un costo computazionaledi n2.

In seguito vedremo dei metodi numerici per la risoluzionedi un sistema lineare che sono molto piu efficienti dal puntodi vista computazionale.



I metodi numerici per la risoluzione di un sistema lineare sisuddividono in due classi:

• Metodi diretti: in assenza di errori diarrotondamento la soluzione viene calcolataesattamente in un numero finito di passi;

• Metodi iterativi: attraverso un processo iterativoviene generata una successione infinita di vettoriconvergente, sotto opportune condizioni, alla soluzionecercata. Anche in assenza di errori di arrotondamentonon si ottiene mai la soluzione esatta perche ad uncerto punto e necessario arrestare il processo.

I criteri di scelta di un algoritmo sono:

• stabilita

• occupazione di memoria

• costo computazionale

Non esiste un algoritmo “migliore” degli altri in assoluto.

La scelta viene fatta in funzione delle caratteristiche dellamatrice A del sistema (“sparsita”, dimensione, struttura).



Definizione Una matrice A ∈ Rn×n si dice densa se lamaggior parte dei suoi elementi e non nullo.

Definizione Una matrice A ∈ Rn×n si dice sparsa se ilnumero dei suoi elementi non nulli e dell’ordine di n.

Osservazione Per memorizzare in un calcolatore unamatrice sparsa (soprattutto se n e molto grande) e possibileutilizzare solo tre vettori di lunghezza n: uno permemorizzare gli elementi non nulli, uno per memorizzare icorrispondenti indici di riga ed uno per memorizzare icorrispondenti indici di colonna.



I metodi diretti risolvono il problema inizialetrasformandolo in problemi equivalenti aventi la stessasoluzione ma matrici dei coefficienti differenti.

In generale, se la matrice A del sistema e densa i metodidiretti sono piu efficienti.

Se, invece, la matrice A del sistema e sparsa, i metodidiretti non preservano la sparsita delle matrici coinvolte nelprocedimento, anzi con il procedere del metodo il numerodegli elementi non nulli presenti in queste matricigeneralmente cresce e (se n e molto grande) puo ben prestosaturare la memoria centrale del calcolatore.

In questi casi e utile e, spesso indispensabile, utilizzare imetodi iterativi perche, invece, quest’ultimi lascianoinalterata la matrice A del sistema.



Condizionamento del problema

Supponiamo di introdurre una perturbazione ∆A ∈ Rn×n

sulla matrice A e una perturbazione ∆b ∈ Rn sul vettoredei termini noti b.

La soluzione che calcoleremo, supponendo che i numerisiano rappresentabili con infinite cifre e che tutte leoperazioni aritmetiche siano fatte esattamente, non sara x

ma y = x + ∆x, ∆x ∈ Rn, cioe la soluzione sara a sua voltaperturbata. Risolveremo allora il sistema

(A + ∆A)y = (b + ∆b).

Anche se la matrice A e non singolare, A + ∆A potrebberisultare singolare.

Tuttavia, supponendo che A + ∆A e non singolare e che‖A−1‖‖∆A‖ < 1

2 otteniamo la seguente maggiorazione perl’errore relativo indotto sulla soluzione x:

‖δx‖ ≤ 2 cond(A) (‖δA‖+ ‖δb‖) ,

dove cond(A) = ‖A‖‖A−1‖ e ‖ · ‖ denota una fissata normamatriciale indotta.



La quantita cond(A) viene chiamata indice o numero dicondizionamento della matrice A.

cond(A) rappresenta il coefficiente di amplificazione delleperturbazioni indotte sulla matrice A e sul termine noto b.

Se cond(A) e “piccolo” il problema e ben condizionato,altrimenti e mal condizionato.

Osservazione Per qualunque norma matriciale indotta‖ · ‖ si ha

cond(A) = ‖A‖‖A−1‖ ≥ ‖AA−1‖ = ‖I‖ = 1.



Il seguente teorema precisa meglio il significato di cond(A)

Teorema Per ogni matrice non singolare A ∈ Rn×n e perogni norma matriciale indotta, la quantita (cond(A))−1

rappresenta la distanza relativa minima di A dall’insieme ditutte le matrici singolari di ordine n, cioe

1cond(A)

= min{‖A−B‖

‖A‖ : B ∈ Rn×n singolare}



Esempi di matrici mal condizionate sono:

Matrici di Hilbert

Hn =

1 12

13 · · · 1

n

12

13

14 · · · 1

n+1

13

14

15 · · · 1

n+2...

......

......

1n

1n+1

1n+2 · · · 1

2n−1



Matrici di Vandermonde

Vn =

1 1 1 · · · 1

x1 x2 x3 · · · xn

x21 x2

2 x23 · · · x2

n

......

......

...

xn−11 xn−1

2 xn−13 · · · xn−1

n

xi 6= xj per i 6= j



Richiami sulle Norme di Vettori e Matrici

Una norma vettoriale e una funzione che associa ad unvettore x ∈ Rn un numero reale ‖x‖ con le seguentiproprieta:

1. ‖x‖ > 0, ∀x 6= 0 e ‖x‖ = 0 ⇔ x = 0;

2. ‖cx‖ = |c|‖x‖, ∀c ∈ R;

3. ‖x + y‖ ≤ ‖x‖+ ‖y‖, ∀y ∈ Rn.

Le norme vettoriali piu frequentemente utilizzate sono:

• norma infinito:

‖x‖∞ = max1≤i≤n

|xi|

• norma 1:

‖x‖1 =n∑

i=1

|xi|

• norma euclidea:

‖x‖2 =

(n∑

i=1

|xi|2) 1

2

=√

xT x



Analogamente, una norma matriciale e una funzione cheassocia ad una matrice A ∈ Rn×n il numero reale ‖A‖ taleche

1. ‖A‖ > 0, ∀A 6= 0 e ‖A‖ = 0 ⇔ A = 0;

2. ‖cA‖ = |c|‖A‖, ∀c ∈ R;

3. ‖A + B‖ ≤ ‖A‖+ ‖B‖, ∀B ∈ Rn×n;

4. ‖AB‖ ≤ ‖A‖‖B‖, ∀B ∈ Rn×n.

Le norme matriciali piu note sono:

• norma infinito: ‖A‖∞ = max1≤i≤n

n∑

j=1

|ai,j |

• norma 1: ‖A‖1 = max1≤j≤n

n∑

i=1

|ai,j |

• norma spettrale: ‖A‖2 =√

ρ(AT A),dove ρ(A) e il raggio spettrale della matrice A, cioe ilmodulo massimo degli autovalori di A;

• norma di Frobenius: ‖A‖F =

n∑

i=1

n∑

j=1

|ai,j |2

12

.



Data una norma vettoriale e una norma matriciale, si diceche le due norme sono compatibili se

‖Ax‖ ≤ ‖A‖‖x‖, ∀A ∈ Rn×n, x ∈ Rn.

Una norma vettoriale puo risultare compatibile con piunorme matriciali.

Ad ogni norma vettoriale possiamo associare una normamatriciale nel seguente modo:

‖A‖ = supx 6=0

‖Ax‖‖x‖ = max

‖x‖=1‖Ax‖.

Una norma cosı definita viene chiamata naturale oindotta da quella del vettore.

Per le norme matriciali naturali abbiamo inoltre

‖I‖ = max‖x‖=1

‖Ix‖ = 1.

Le norme matriciali ‖ · ‖∞, ‖ · ‖1 e ‖ · ‖2 sono naturali. Lanorma ‖ · ‖F non e naturale in quanto

‖I‖F = n12 .



Si ha

‖x‖∞ induce ‖A‖∞

‖x‖1 induce ‖A‖1

‖x‖2 induce ‖A‖2



Metodi Diretti

Sistemi Triangolari

Siano

x = (x1x2 · · ·xn)T e b = (b1b2 · · · bn)T

Denotiamo con

L =

l1,1 0 0 · · · 0

l2,1 l2,2 0 · · · 0...

.... . . . . .

...

ln−1,1 ln−1,2 · · · ln−1,n−1 0

ln,1 ln,2 · · · ln,n−1 ln,n

matrice triangolare inferiore (Low)



e con

U =

u1,1 u1,2 · · · u1,n−1 u1,n

0 u2,2 · · · u2,n−1 u2,n

0. . . . . .

......

.... . . 0 un−1,n−1 un−1,n

0 · · · 0 0 un,n

matrice triangolare superiore (Up)

Vogliamo risolvere i sistemi

Lx = b e Ux = b

Osservazione

det(L) = l1,1l2,2 · · · ln,n 6= 0 ⇔ li,i 6= 0 ∀i = 1, . . . , n

det(U) = u1,1u2,2 · · ·un,n 6= 0 ⇔ ui,i 6= 0 ∀i = 1, . . . , n



Algoritmo di sostituzione in avanti (forwardsubstitution)

Lx = b

l1,1x1 = b1

l2,1x1 + l2,2x2 = b2

l3,1x1 + l3,2x2 + l3,3x3 = b3

......

ln,1x1 + ln,2x2 + ln,3x3 + · · ·+ ln,nxn = bn

x1 =b1

l1,1

x2 =b2 − l2,1x1

l2,2

x3 =b3 − l3,1x1 − l3,2x2

l3,3

...

xi =

bi −i−1∑

k=1

li,kxk

li,i, i = 4, . . . , n



x1 = b1/l1,1; 1 operazione

for i=2:n

xi = bi;

for k=1:i-1

xi = xi − li,k ∗ xk; i-1 operazioni

end

xi = xi/li,i; 1 operazione

end

costo computazionale

n∑

i=1

i =n(n + 1)

2' n2

2



Algoritmo di sostituzione all’indietro (backwardsubstitution)

Ux = b

u1,1x1 + u1,2x2 + u1,3x3 + · · ·+ u1,nxn = b1

u2,2x2 + u2,3x3 + · · ·+ u2,nxn = b2

......

un−1,n−1xn−1 + un−1,nxn = bn−1

un,nxn = bn

xn =bn

un,n

xn−1 =bn−1 − un−1,nxn

un−1,n−1

xn−2 =bn−2 − un−2,n−1xn−1 − un−2,nxn

un−2,n−2

...

xi =

bi −n∑

k=i+1

ui,kxk

ui,i, i = n− 3, n− 4, . . . , 2, 1



xn = bn/un,n; 1 operazione

for i=n-1:-1:1

xi = bi;

for k=i+1:n

xi = xi − ui,k ∗ xk; n-i operazioni

end

xi = xi/ui,i; 1 operazione

end

costo computazionale

1∑

i=n

n− i + 1 =n∑

i=1

i =n(n + 1)

2' n2

2



Osservazione

Gli Algoritmi di sostituzione in avanti e all’indietro sonoben definiti se le matrici L e U sono non singolari, cioe

li,i 6= 0, ui,i 6= 0, ∀i = 1, . . . , n



Metodo di eliminazione di GaussSia A ∈ Rn×n una generica matrice di ordine n e x, b ∈ Rn

due vettori colonna di lunghezza n. Scriviamoesplicitamente il sistema Ax = b.

a1,1x1 + a1,2x2 + a1,3x3 · · ·+ a1,nxn = b1

a2,1x1 + a2,2x2 + a2,3x3 · · ·+ a2,nxn = b2

......

......

an,1x1 + an,2x2 + an,3x3 · · ·+ an,nxn = bn

Supponendo che a1,1 6= 0. Possiamo eliminare l’incognita x1

dalla 2a, 3a, . . . , n−esima equazione sommando all’i−esimaequazione i = 2, . . . , n la prima equazione moltiplicata per

mi,1 = − ai,1

a1,1, i = 2, . . . , n.

Eseguiamo, cioe, le seguenti combinazioni lineari

a(2)i,j = ai,j + mi,1a1,j

i, j = 2, . . . , n

b(2)i = bi + mi,1b1



Il nuovo sistema, equivalente a quello precedente, diventa

a1,1x1 + a1,2x2 + a1,3x3 + · · ·+ a1,nxn = b1

a(2)2,2x2 + a

(2)2,3x3 + · · ·+ a

(2)2,nxn = b

(2)2

a(2)3,2x2 + a

(2)3,3x3 + · · ·+ a

(2)3,nxn = b

(2)3

......

......

a(2)n,2x2 + a

(2)n,3x3 + · · ·+ a

(2)n,nxn = b

(2)n

Supponendo che a2,2 6= 0. Possiamo eliminare l’incognita x2

dalla 3a, 4a, . . . , n−esima equazione sommando all’i−esimaequazione i = 3, . . . , n la seconda equazione moltiplicata per

mi,2 = −a(2)i,2

a(2)2,2

, i = 3, . . . , n.

Eseguiamo, cioe, le seguenti combinazioni lineari

a(3)i,j = a

(2)i,j + mi,2a

(2)2,j

i, j = 3, . . . , n

b(3)i = b

(2)i + mi,2b

(2)2



Il nuovo sistema, equivalente a quello precedente, diventa

a1,1x1 + a1,2x2 + a1,3x3 + · · ·+ a1,nxn = b1

a(2)2,2x2 + a

(2)2,3x3 + · · ·+ a

(2)2,nxn = b

(2)2

a(3)3,3x3 + · · ·+ a

(3)3,nxn = b

(2)3

......

......

a(3)n,3x3 + · · ·+ a

(3)n,nxn = b

(2)n

Gli elementi a1,1, a(2)2,2, a

(3)3,3, . . . , che compaiono durante le

successive eliminazioni vengono detti elementi pivot.

Dopo (n− 1) passi, supponendo che tutti gli elementi pivotsiano non nulli, si ottiene il seguente sistema triangolaresuperiore equivalente ai precedenti

a1,1x1 + a1,2x2 + a1,3x3 + · · ·+ a1,nxn = b1

a(2)2,2x2 + a

(2)2,3x3 + · · ·+ a

(2)2,nxn = b

(2)2

a(3)3,3x3 + · · ·+ a

(3)3,nxn = b

(3)3

......

a(n)n,nxn = b

(n)n

che puo essere risolto con l’Algoritmo di sostituzioneall’indietro.



Esempio Consideriamo il seguente sistema

2x1 − x2 + x3 − 2x4 = 0

2x2 − x4 = 1

x1 − 2x3 + x4 = 0

2x2 + x3 + x4 = 4

applicando il metodo di Gauss con

m2,1 = 0, m3,1 = −12

m4,1 = 0

si ottiene il sistema

2x1 − x2 + x3 − 2x4 = 0

2x2 − x4 = 112x2 − 5

2x3 + 2x4 = 0

2x2 + x3 + x4 = 4



Conm3,2 = −1

4m4,2 = −1

si ha

2x1 − x2 + x3 − 2x4 = 0

2x2 − x4 = 1

− 52x3 + 9

4x4 = − 14

x3 + 2x4 = 3

Infine, con

m4,3 =25

si ottiene

2x1 − x2 + x3 − 2x4 = 0

2x2 − x4 = 1

− 52x3 + 9

4x4 = − 14

2910x4 = 29

10

Da cui, sostituendo all’indietro, si ha

x4 = 1, x3 =(−1

4− 9

4x4

)(−2

5

)= 1, x2 =

1 + x4

2= 1,

x1 =x2 − x3 + 2x4

2= 1



Ponendo

ai,j = a(1)i,j e bi = b

(1)i , i, j = 1, . . . , n,

il seguente schema di calcolo riassume il metodo dieliminazione di Gauss:

Per k = 1, 2, . . . , n− 1

mi,k = −a(k)i,k

a(k)k,k

a(k+1)i,j = a

(k)i,j + mi,ka

(k)k,j , i, j = k + 1, . . . , n

b(k+1)i = b

(k)i + mi,kbk

OsservazioneIl metodo di Gauss si puo applicare se e soltanto se ad ognipasso k = 1, . . . , n− 1 l’elemento pivot a

(k)k,k 6= 0.



Consideriamo le sottomatrici di A

Ak = (ai,j)i,j=1,...,k, k = 1, . . . , n,

dette matrici principali di testa di A. Vale il seguenteteorema

Teorema Sia A ∈ Rn×n. Gli elementi pivota(k)k,k, k = 1, . . . , n, sono tutti diversi da zero se e soltanto se

tutte le matrici principali di testa di A sono non singolari,cioe

det(Ak) 6= 0, ∀k = 1, . . . , n.



Esempi di matrici per le quali la condizione

det(Ak) 6= 0, ∀k = 1, . . . , n,

e verificata sono:

• matrici a diagonale dominante per righe, cioe

|ai,i| >n∑

j=1j 6=i

|ai,j |, i = 1, . . . , n;

• matrici a diagonale dominante per colonne, cioe

|aj,j | >n∑

i=1j 6=i

|ai,j |, j = 1, . . . , n, ;

• matrici simmetriche definite positive, cioe

AT = A e xT Ax > 0, ∀x ∈ Rn, x 6= 0.



· · ·for k=1:n-1

for i=k+1:n

ai,k = −ai,k/ak,k; 1 operazione

for j=k+1:n

ai,j = ai,j + ai,k ∗ ak,j; n-k operazioni

end

bi = bi + ai,k ∗ bk; 1 operazione

end

end

U=triu(A);

Si risolve il sistema triangolare superiore

Ux = b con l’algoritmo di sostituzione

all’indietro

costo computazionalen−1∑

k=1

[2(n− k) + (n− k)2] +n2

2= 2

n−1∑

i=1

i +n−1∑

i=1

i2 +n2

2

= 2n(n− 1)

2+

n(n− 1)(2n− 1)6

+n2

2

=n3

3+ n2 − 5

6n ≤ n3

3+ n2 ' n3

3



Pivoting e stabilita

Se al k−esimo passo si trova a(k)k,k = 0, il metodo di Gauss,

cosı come e stato descritto, non puo proseguire.

Se a(k)k,k = 0 necessariamente qualche altro elemento

a(k)i,k , i = k + 1, . . . , n, della colonna k−esima della matrice

dei coefficienti deve essere non nullo, altrimenti la matricedei coefficienti sarebbe singolare.

Se, ad esempio, a(k)r,k 6= 0, basta scambiare l’equazione

k−esima con la r−esima e poi procedere con leeliminazioni.

Dunque, ogni sistema non singolare puo essere semprericondotto alla forma triangolare superiore con il metodo diGauss piu eventuali scambi di equazioni.



Per assicurare una migliore stabilita numerica al metodo dieliminazione di Gauss e spesso necessario permutarel’ordine delle equazioni anche quando l’elemento pivot none esattamente zero, ma e molto piccolo (in valore assoluto)rispetto agli altri elementi.

Per cercare di evitare catastrofiche propagazioni di errori edi solito necessario scegliere, al generico passo k−esimo,l’elemento pivot seguendo una delle seguenti strategie:

• Pivoting parziale: si sceglie r uguale al piu piccolointero ≥ k tale che

|a(k)r,k | = max

k≤i≤n|a(k)

i,k |

e, se r 6= k, si scambia l’equazione k−esima conl’r−esima;

• Pivoting totale: si sceglie la coppia (r, s), con r, s ≥ k

tale che|a(k)

r,s | = maxk≤i,j≤n

|a(k)i,j |

e si scambiamo l’equazione k−esima con l’r−esima el’incognita k−esima (con il suo coefficiente) conl’s−esima.



La strategia di pivoting parziale e la meno “costosa” e,poiche, in generale, risulta soddisfacente, essa e la strategiapiu utilizzata.

Osservazione Il metodo di eliminazione di Gauss senzapivoting e comunque numericamente stabile quando:

• la matrice A del sistema e simmetrica e a diagonaledominante (la strategia di pivot non produce scambi);

• la matrice A del sistema e simmetrica e definitapositiva.



Esempio Consideriamo il sistema lineare Ax = b di ordinen = 18, dove

ai,j = cos(

(j − 1)2i− 1

2nπ

), i, j = 1, . . . , n,

e

bi =n∑

j=1

ai,j , i = 1, . . . , n,

la cui soluzione esatta e x = (1, 1, . . . , 1)T .

cond(A) = ‖A‖∞‖A−1‖∞ = 16.90251471518910

Risolvendo il sistema con il metodo di eliminazione diGauss senza e con la strategia di pivoting parziale, siottengono i seguenti risultati:



Gauss Gauss + pivoting parziale

9.999999968425205e-001 9.999999999999993e-001

1.000000005848771e+000 1.000000000000000e+000

9.999999953310822e-001 1.000000000000001e+000

1.000000003291504e+000 1.000000000000001e+000

9.999999977773912e-001 9.999999999999997e-001

1.000000001703849e+000 1.000000000000000e+000

9.999999983501708e-001 1.000000000000000e+000

1.000000001767299e+000 1.000000000000001e+000

9.999999982508074e-001 9.999999999999998e-001

1.000000001433063e+000 9.999999999999999e-001

9.999999991425489e-001 9.999999999999993e-001

1.000000000200525e+000 1.000000000000000e+000

1.000000000347843e+000 9.999999999999996e-001

9.999999993163441e-001 9.999999999999990e-001

1.000000000800448e+000 9.999999999999996e-001

9.999999992564082e-001 1.000000000000000e+000

1.000000000564552e+000 1.000000000000000e+000

9.999999996959291e-001 9.999999999999999e-001



Metodo di Gauss e Fattorizzazione LU

Il metodo di eliminazione di Gauss consiste nel costruireuna successione di matrici

[A|b] = [A(1)|b(1)] . . . [A(k)|b(k)] . . . [A(n)|b(n)]

in modo tale che, per k = 1, . . . , n− 1,

[A(k+1)|b(k+1)] = M (k)[A(k)|b(k)]

dove

M (k) =

1 0 . . . . . . 0

0 1 0...

.... . . . . . . . .

0 0 1...

0 0 mk+1,k 1. . .

0 0 mk+2,k 0. . .

......

......

. . . . . . 0

0 . . . 0 mn,k 0 . . . 0 1

e la matrice elementare di Gauss.



Si ha

A = [M (1)]−1A(2) = [M (1)]−1[M (2)]−1A(3) = · · ·= [M (1)]−1 · · · [M (n−1)]−1A(n)

eb = [M (1)]−1 · · · [M (n−1)]−1b(n).

Ponendo

L = [M (1)]−1 · · · [M (n−1)]−1, U = A(n) e y = b(n),

otteniamoA = LU e b = Ly,

dove L e una matrice triangolare inferiore e U e unamatrice triangolare superiore.

Dunque, il sistemaAx = b

diventaLUx = Ly

e, dopo aver applicato il metodo di eliminazione di Gauss,rimane da risolvere il sistema triangolare superiore

Ux = y.



Osservazione Il metodo di eliminazione di Gauss calcolaesplicitamente la matrice U . La matrice L e

L =

1 0 . . . . . . 0

−m2,1 1 0...

.... . . . . . . . .

. . . 1...

... −mk+1,k 1. . .

−mk+2,k. . . . . .

......

. . . 1 0

−mn,1 . . . . . . −mn,k . . . . . . −mn,n−1 1

Calcolate le matrici L e U , per risolvere il sistema Ax = b,basta risolvere i sistemi

Ly = b e Ux = y,

il primo triangolare inferiore e il secondo triangolaresuperiore.

Il costo computazionale e n3

3 + n2 di cui n3

3 − n3 ≤ n3

3 per lasola fattorizzazione A = LU .



Eseguire la fattorizzazione LU di A puo essere utile se:

• e necessario risolvere p sistemi

Ax1 = b1, Ax2 = b2, . . . , Axp = bp,

cioe il sistema

AX = B, con A ∈ Rn×n, X,B ∈ Rn×p

Il costo computazionale p(

n3

3 + n2)

si riduce an3

3 + pn2.

Se p = n e B = I, risolvere il sistema AX = B, eequivalente a calcolare A−1.

• si vuole calcolare il det A. Infatti

det(A) = det(LU) = det(L) det(U) = det(U) =n∏

i=1

ui,i.



Il metodo di Gauss con la variante del pivot esegue

ΠA = LU e Πb = Ly

dove Π ∈ Rn×n, detta matrice di permutazione,contiene le informazioni relative agli scambi di righe.

Vale il seguente

Teorema Per ogni matrice A ∈ Rn×n esiste una matricedi permutazione Π ∈ Rn×n tale che

ΠA = LU.



Metodo di Cholesky

Sia A ∈ Rn×n una matrice simmetrica definita positiva.

Vale il seguente

Teorema Se A ∈ Rn×n e una matrice simmetrica definitapositiva esiste ed e unica la fattorizzazione

A = LLT .

Posto A = (ai,j)i,j=1,...,n e L = (li,j)i,j=1,...,n, si ha

ai,j =n∑

k=1

li,klTk,j =n∑

k=1

li,klj,k.

Poiche la matrice e simmetrica, possiamo considerare sologli elementi di A con j ≤ i. Otteniamo

ai,j =j−1∑

k=1

li,klj,k + li,j lj,j , i = 1, . . . , n, j = 1, . . . , i− 1

e

ai,i =i−1∑

k=1

l2i,k + l2i,i, i = 1, . . . , n.



Da cui

li,i =

√√√√ai,i −i−1∑

k=1

l2i,k

i = 1, . . . , n

j = 1, . . . , i− 1

li,j =1

lj,j

[ai,j −

j−1∑

k=1

li,klj,k

]



· · ·for i=1:n

li,i = 0;for k=1:i-1

li,i = li,i + l2i,k; i− 1 operazioni

end

li,i =√

ai,i − li,i; 1 operazione

for j=1:i-1

li,j = 0;for k=1:j-1

li,j = li,j + li,k ∗ lj,k; j − 1 operazioni

end

li,j = (ai,j − li,j)/lj,j; 1 operazione

end

end

Si risolve il sistema triangolare inferiore

Ly = b con l’algoritmo di sostituzione

in avanti e il sistema triangolare superiore

LT x = y con l’algoritmo di sostituzione

all’indietro



Costo computazionale

n∑

i=1

i +

i−1∑

j=1

j

+ n2 =

n∑

i=1

i +n∑

i=1

i(i− 1)2

+ n2

=n∑

i=1

i +12

n∑

i=1

i2 − 12

n∑

i=1

i + n2

=12

n∑

i=1

i +12

n∑

i=1

i2 + n2

=n(n + 1)

4+

(n + 1)n(2n + 1)12

+ n2

=n3

6+

32n2 +

n

3

' n3

6


Corso di Calcolo Numerico - old · maggior parte dei suoi elementi `e non nullo. Deﬁnizione Una...

Documents

Transcript of Corso di Calcolo Numerico - old · maggior parte dei suoi elementi `e non nullo. Deﬁnizione Una...