Appunti di Calcolo Numerico

J −(L

N = + U )

J 21

La generica equazione del sistema (3.11) diventa

n

(k+1) (k)

X

− ∀i

a x = a x + b = 1, 2, . . . , n

ii ij i

i j

j=1

j6 = i

dove a è il generico elemento sulla diagonale, mentre gli elementi a all’in-

ii ij

terno della sommatoria sono i termini sulla matrice triangolare superiore ed

inferiore. (k)

Sommando e sottraendo il termine a x

ii i

n

(k+1) (k) (k) (k)

X

− −

a x = a x + b + a x a x

ii ij i ii ii

i j i i

j=1

j6 = i

Portando il termine sottratto all’interno della sommatoria questa risulta

n n

(k) (k) (k)

X X

−

a x a x = a x

ij ii ij

j i j

j=1 j=1

j6 =

i

Allora, risulta n

(k+1) (k) (k)

X

−

a x = a x a x + b

ii ii ij i

i i j

j=1

(k+1)

e isolando x dividendo per a ii

i (k)

n

P −

a x b

ij j

j=1 j

(k+1) (k) − (3.15)

x = x

i i a

ii

Metodo di Gauss - Seidel

Perché la formula generale (3.14) soddisfi il sistema (3.12) deve essere

M = L + D (matrice triangolare inferiore)

G −U

N =

G

Metodo S.O.R.

Nel metodo S.O.R. si deve risultare

M = D + ω L

ω − −

N = (1 ω) D ω U

ω

Condizione necessaria per la convergenza è che il raggio spettrale (l’autova-

lore di massimo modulo) della matrice di iterazione sia strettamente minore

di 1 −1

ρ(M N ) < 1

22

Per il metodo S.O.R. si deve scegliere il valore di ω in modo tale che il raggio

spettrale sia il più piccolo possibile.

Inoltre, per semplificare i calcoli, è conveniente che la matrice M sia

facilmente invertibile. Tutti e tre i metodi soddisfano questo requisito: in

Jacobi M è una matrice diagonale mentre in Gauss - Seidel e nel S.O.R. M

è triangolare.

4 Sistemi debolmente non lineari

I sistemi debolmente non lineari sono sistemi composti da una parte

lineare e una parte che può essere non lineare.

 · · ·

a x + a x + + a x + f (x ) = 0

11 1 12 2 1n n 1 1







 · · ·

a x + a x + + a x + f (x ) = 0

 21 1 22 2 2n n 2 2





 · · ·

a x + a x + + a x + f (x ) = 0 (4.1)

31 1 32 2 3n n 3 3

..



 .









 · · ·

a x + a x + + a x + f (x ) = 0

 n1 1 n2 2 nn n n n

Teorema 1.4. ”Se la matrice A della parte lineare del sistema (4.1) è

tridiagonale a diagonale dominante con

∀i

a < 0 = 1, 2, . . . , n

ii ∀j

a > 0, a > 0 = 1, 2, . . . , n

j,j+1 j+1,j 0

∀i ≤ ∀x ∈

e se inoltre f (x) = 1, 2, . . . , n sono differenziabili e f (x) 0 R

i i

allora la soluzione del sistema (4.1) esiste ed è unica.

Dimostrazione: unicità della soluzione (per assurdo). Si supponga che esi-

n

∈

stano due soluzioni x, y .

R

Si assuma che la generica equazione sia soddisfatta dal vettore x. Si

supponga poi la medesima cosa per il vettore y. Questo implica che

n

X ∀i

a x + f (x ) = 0 = 1, 2, . . . , n

ij j i i

j=1

n

X ∀i

a y + f (y ) = 0 = 1, 2, . . . , n

ij j i i

j=1

Sottraendole tra loro si ottiene

n

X − −

a (y x ) + f (y ) f (x ) = 0

ij j j i i i i

j=1 23

Si definisca il valore −

e = y x

j j j

Per il teorema del valore medio si ha

0 0

− −

f (y ) f (x ) = f (ξ ) (y x ) = f (ξ ) e

i i i i i i i i i

i i

dove ξ è un valore a metà tra y e x .

i i i

Allora diventa n

X 0 ∀i

a e + f (ξ ) e = 0 = 1, 2, . . . , n (4.2)

ij j i i

i

j=1 ∀j

dove le incognite sono e = 1, 2, . . . , n. Il sistema (4.2) risulta lineare e

j

omogeneo. 0 ≤

Dato che f (ξ ) 0 il sistema soddisfa il teorema 1.2: sulla diagonale

i

i 0 ≤

a + f (ξ ) a < 0

ii i ii

i

allora la soluzione esiste − ∀j

e = y x = 0 = 1, 2, . . . , n

j j j

Questo implica che le soluzioni x e y coincidono



x = y

1 1









x = y

 2 2

. . .









x = y

 n n

Si è dimostrato che la soluzione esiste ed è unica.

5 Metodi del gradiente

Si consideri il sistema lineare Ax = b T

e si supponga che la matrice A sia simmetrica (A = A ) e che sia definita

T n T

≥ ∀x ∈

positiva (x A x 0 e inoltre x A x = 0 se e solo se x = 0).

R

Si consideri ora il sistema lineare omogeneo

Ax = 0

24

6

soddisfatto per un vettore x = 0. Allora

T T

x A x = x 0 = 0 (Assurdo!)

6

Non esiste alcun vettore x = 0 che soddisfi il sistema omogeneo. Allora, la

soluzione esiste ed è unica! n

∈

Si supponga ora che, dato un vettore y R

T T

x A y = y A x (possibile perché A è simmetrica)

infatti T T T T T

x A y = (A y) x = y A x = y A x

Si consideri ora il seguente sistema

Ax̄ = b (5.1)

dove x̄ è la soluzione esatta del sistema.

−

Posto che x x̄ = e si consideri la funzione

T T

− −

Φ(x) = e A e = (x x̄) A(x x̄)

Notazione:

il prodotto scalare può essere riscritto come

T

x A y = (x, A y)

T

x y = (x, y)

n

∈

per ogni vettore x, y .

R

Detto questo, la funzione Φ(x) si riscrive come

− − − − − − ≥

Φ(x) = ((x x̄), A(x x̄)) = (x x̄, A x A x̄) = (x x̄, A x b) 0

≥

Poiché A è definita positiva, risulta Φ(x) 0 con uguaglianza se e solo se

e = 0 e cioè x = x̄.

A questo punto, se x è soluzione risulta che A x−b = 0; nel caso contrario,

si può definire il vettore dei residui come

−

A x b = r(x)

25

x̄ (0)

x

Figura 3: Esempio di sistema a due variabili

5.1 Metodo del gradiente semplice

Si supponga di essere nel piano (n = 2). In due dimensioni, la funzione

Φ(x) risulta Φ(x) = c

che rappresenta i punti sull’ellisse di centro x̄ (vedi figura 3).

La risoluzione del sistema A x = b equivale a minimizzare la funzione

Φ(x), cioè calcolare il minimo della funzione (x̄). (0)

Per trovare x̄ si parte da un generico punto x sull’ellisse e si indivi-

dua a direzione del gradiente (o di massima pendenza) di Φ(x) (che risulta

ortogonale alla curva) lungo cui muoversi.

∇Φ(x) − − − − −

= A x b + A(x x̄) = A x b + A x b = 2(A x b) = 2 r(x)

(0) n

∈

Scelto x si ricava

R (k+1) (k) (k)

x = x + α r

k

(k) (k) (k) −

dove r = r(x ) = A x b, mentre il parametro α indica la posizione

k

in cui fermarsi lungo la direzione del gradiente. (k+1)

Si prosegue scegliendo α in modo da minimizzare la funzione Φ(x )

k (k+1)

dΦ(x ) =0

dα

k

(k+1) (k+1) (k+1)

− −

Φ(x ) =((x x̄), A x b) =

(k) (k) (k) (k)

− −

=((x + α r x̄), A(x + α r ) b) =

k k

(k) (k) (k) (k)

− −

=((x + α r x̄), A x b + α A r ) =

k k

(k) (k) (k) (k)

−

=((x + α r x̄), r + α A r )

k k

26

Allora, il minimo si calcola come

(k+1)

dΦ(x ) (k) (k) (k) (k) (k) (k)

−

0= =(r , r + α A r ) + ((x + α r x̄), A r ) =

k k

dα k (k) (k) (k) (k) (k) (k) (k)

−

=(r , r ) + α (r , A r ) + (A(x + α r x̄), r ) =

k k

(k) (k) (k) (k) (k) (k) (k)

−

=(r , r ) + α (r , A r ) + ((A x A x̄ + α A r ), r ) =

k k

(k) (k) (k) (k) (k) (k) (k)

=(r , r ) + α (r , A r ) + (r + α A r , r ) =

k k

(k) (k) (k) (k) (k) (k) (k) (k)

=(r , r ) + α (r , A r ) + (r , r ) + α (r , A r ) =

k k

(k) (k) (k) k

=2[(r , r ) + α (r , A r )]

k

Risolvendo, risulta (k) (k)

(r , r )

− (5.2)

α =

k (k) (k)

(r , A r )

Data la convergenza assicurata se A è simmetrica e definita positiva, la

soluzione è (k) (k)

(r , r )

(k+1) (k) (k)

− ∀k

x = x r = 0, 1, . . . , K (5.3)

(k) (k)

(r , A r )

(k) (k) −

con r = A x b

5.2 Metodo del gradiente coniugato

Quanto vale K?

(0)

Se x è sulla stessa direzione del gradiente di Φ(x) passante per x̄, il

sistema si risolve in un’unica iterazione.

Se cosı̀ non fosse, si utilizzano delle direzioni coniugate per velocizzare

la convergenza.

(0) n

∈

Scelto x si pone

R (0) (0) (0) −

p = r = A x b

Allora, la soluzione generica è

(k+1) (k) (k)

x = x + α p

k

Il vettore p deve soddisfare la seguente condizione

(i) (j) ∀i 6

(p , A p ) = 0 = j

cioè il prodotto scalare tra i due elementi deve risultare identicamente nullo;

(i)

questo implica che i due vettori siano ortogonali! Allora, si dice che p e

(j)

p sono A-coniugati. 27

Anche in questo caso si deve calcolare il minimo della funzione

(k+1) (k) (k) (k) (k)

− −

Φ(x ) = ((x + α p x̄), A(x + α p ) b)

k k

(k) (k)

−

dove A x b = r , allora

(k+1) (k) (k) (k) (k)

−

Φ(x ) = ((x + α p x̄), r + α A p )

k k

La derivata prima rispetto a α risulta

k

(k+1)

dΦ(x ) (k) (k) (k) (k) (k) (k)

−

=(p , r + α A p ) + ((x + α p x̄), A p ) =

0= k k

dα k (k) (k) (k) (k) (k) (k)

−

=(p , r + α A p ) + (A(x + α p x̄), p ) =

k k

(k) (k) (k) (k) (k) (k)

−

=(p , r + α A p ) + ((A x A x̄ + α A p ), p )

k k

(k) (k)

−

considerando che A x̄ = b e che A x b = r , si ottiene

(k+1)

dΦ(x ) (k) (k) (k) (k) (k) (k)

· · · −

= = (p , r + α A p ) + ((A x b + α A p ), p ) =

0= k k

dα k (k) (k) (k) (k) (k) (k)

=(p , r + α A p ) + ((r + α A p ), p ) =

k k

(k) (k) (k) (k) (k) (k) (k) (k)

=(p , r ) + α (p , A p ) + (p , r ) + α (p , A p ) =

k k

(k) (k) (k) (k)

=2[(p , r ) + α (p , A p )]

k

Il valore di α per cui si annulla è

k (k) (k)

(p , r )

−

α =

k (k) (k)

(p , A p )

(k+1)

Mentre, il valore di r si calcola come

(k+1) (k+1) (k) (k) (k) (k)

− − −

r = A