Che materia stai cercando?

Anteprima

ESTRATTO DOCUMENTO

Proof :

h1i1. Dim 1.

Osserviamo che kbk

⇒ kbk kAxk ≤ kAkkxk ⇒ ≤ kAk

b = Ax = (1)

kxk

|{z}

[3.8,pg.10] −1

⇔ ⇔ ⇔

A(x + ∆x) = b + ∆b Ax + A∆x = b + ∆b A∆x = ∆b ∆x = A ∆b kbk

k∆bk k∆bk

k∆xk

−1 −1 −1

−1

k∆xk kA ≤ ≤ kA k kA k ≤

kA kk∆bk ⇔

= ∆bk =

kxk kxk kbk kxk

|{z}

[3.8,pg.10]

k∆bk k∆bk

−1

≤ kA kkAk = µ(A)

kbk kbk

|{z}

(1)

h1i2. Dim 2. −1

kAk kA k

Poiche’ µ (A) = , basta usare le prop [3.10,pg.11], [3.11,pg.12].

2 2 2

4.1 Fattorizzazione LU

n×n

Theorem 4.5. Data A , si ha

R

6 ∀k ⇒

D = 0 = 1, . . . , n A = LU

k

dove L, U sono le seguenti matrici triangolari :

   

1 a a a ... a

11 12 13 1n

2 2

m 1 0 a a ... a

21 23

 

 

22 2n

 

 3 3

m 1

m 0 0 a ... a

L = , U =

31 32 

 

 33 3n

 

... . . .

 

 n

m m . . . m 1 0 0 0 ... a

n1 n2 n,n−1 nn

k

Dove m , a sono ottenuti tramite l’algoritmo [4.2,pg.13].

ij ij 6 ∀k

Proof : Per la prop [4.3,pg.16], la condizione D = 0 = 1, . . . , n assicura che e’ possibile

k

utilizzare l’algoritmo [4.2,pg.13], ovvero

k 6 ∀k

a = 0 = 1, . . . , n

kk

Sia

0 0 ... 0 m m . . . m

i+1,i i+2,i n,i

m = |{z}

i posto i-esimo

( ≤

0 j i

= (m ) = (1)

m j

ij i j>i

m ji

t

M = I m e

i i

i

n

{e }

con base canonica di .

R

i

h1i1. Si ha 

1 i = j

−m

(M ) = i > h, j = h

h ij ji

0 altrimenti

Ovvero, M e’ la matrice triangolare che ha tutti 1 nella diagonale principale, come colonna

h 18

t

h ha m , tutti gli altri elementi sono 0:

i 1

 

0 1

 

..

 

 

. 

 

0 0 ... 0 1

 

 

−m

0 0 ... 0

M = k+1,k

h  

 −m

0 0 ... 0 k+2,k

 

..

 

 

.

 

 −m 1

0 0 ... 0 n,k

 

{z }

|

colonna k

n×1 1×n t n×n

t ∈ ∈ ∈

, e , si ha m e . Ed e’

Osserviamo che poiche’ m R R R

h h

h h

( (

0 i h m i > h, j = h

ih

t t t

( m e ) = m e = m δ = =

h ij hj hj

h hi hi m δ i > h 0 altrimenti

|{z} ih hj

(1) 

1 i = j

( 

m i > h, j = h 

ih

t

− − −m

(M ) = (I m e ) = δ = i > h, j = h

h ij h ij ij ih

h 0 altrimenti 

0 altrimenti

h

h1i2. Si ha M . . . M A = D

h−1 1 n h

h e’ il minore principale di ordine k di A .

Dove D

k A

Consideriamo M 1 

A i = r

rj

n

n 

X X −m

(M ) A =

(M A) = A i > 1, r = 1

1 ir rj

1 ij ri rj

|{z} 

r=1 r=1 0 altrimenti

h1i1 

passo

   

1 a ... a

11 1n

 

a . . . a

11 1n

−m − −

1 a m a . . . a m a

21 21 21 11 2n 21 1n

   

a . . . a

21 2n

  

  

−m − −

0 1 a m a . . . a m a

M A = =

31 31 31 11 3n 31 1n

   

 

..

1    

 

.. ..

.

   

 

. .

   

. . . a

a

n1 nn

−m − −

1 a m a . . . a m a

n1 n1 n1 11 nn n1 1n

E ricordando che m e’ il moltiplicatore di posto i, j si ha

ij 2 2 2

2 

 a a . . . a

a

a a ... a 11 12 13 1n

11 12 1n 2 2 2

0 a a . . . a

− −

0 a m a . . . a m a 

 22 23 2n

22 21 12 2n 21 1n

  2

 2 2 2

0 a a . . . a

M A = = = D

 

.

1  

32 33 3n n

.

 

. 

 . . .

  

− −

0 a m a . . . a m a 2 2 2

0 a a . . . a

n2 n1 12 nn n1 1n nn

n2 n3

Cosi’ procedendo si arriva a h

M . . . M A = D

h−1 1 n

. . . M , Sia U = N A. Si ha

Sia N = M n−1 1 n

U = D

n

|{z}

h1i2

passo

n Quindi, U coincide

D e’ il minore di ordine n dell’ultima matrice data dall’algoritmo [4.2,pg.13].

n

con l’U della tesi.

−1

h1i3. ∃N

Abbiamo gia’ visto che M e’ una matrice triangolare. Poiche’ gli elementi della sua diagonale

i 19

sono tutti 1, si ha −1

|M | ⇒ ∃M ∀i

= 1 = 1, . . . , n

i i

−1 −1 −1

−1

N = M M . . . M

1 2 n−1

−1

Sia L = N , si ha −1

U = N A A = N U = LU

h1i4. Resta da controllare che L sia proprio quello della tesi

−1 t

h2i1. = I + m

M e

i i

i

Infatti, t t t t t

− −

M (I + m e ) = (I m e )(I + m e ) = I m e e =

m

i i i i i i

i i i i i

t ·m

e = e = m = 0

m

i i

i i ii |{z}

(1)

= I t e )M = I

(I + m i i

i

n−1 t

P

h2i2. L = I + e

m i

i

i=1

−1 −1 −1

−1 t t t

L = N = M M . . . M = (I + m e )(I + m e ) . . . (I + m e )

1 2 n−1 1 2 n−1

1 2 n−1

t t t t t t

(I + m e )(I + m e ) = I + m e + m e + m e e =

m

1 2 1 2 1 2

1 2 1 2 1 2

t ·m

e m = e = m = 0

1 1

2 2 21 |{z}

(1)

t t

= I + m e + m e

1 2

1 2 n−1 n−2

X X

t t t t t t

· · ·

(I + m e + + (I + m e )(I + m e ) = I + m e + m e m e =

1 n−2 n−1 i i n−1

1 n−2 n−1 i i n−1

i=1 i=1

t ·m

e m = e = m = =0

i i

n−1 n−1 n−1,i |{z}

n−1≥i, (1)

n−1

X t

m

= I + e

i

i

i=1

h2i3. Q.E.D. t

E considerando che le matrici m e sono del tipo:

i

i

 

0

0 0

 

.. 

 

.

 

 

0 0 . . . 0 0

 

t 

0 0 . . . 0 m

m e = i+1,i

i  

i  

0 0 . . . 0 m i+2,i

 

 

.. 

 .

 

 

0 0 . . . 0 m 0

ni 

 |{z}

colonna i

n−1 t

P e coincide con l’L della tesi.

si ha che L = I + m i

i

i=1 k 6

Proposition 4.6. Poniamoci nel contesto della prop [4.2,pg.13], senza supporre che a =

kk

∀k −

0 = 1, . . . , n 1.

Sia s := (1, 2, . . . , n), inizializzato all’inizio dell’algoritmo.

20

Definiamo le seguenti procedure: 0

Sia scambio righe i, i :

T := a

ij

 ∀j := 1, . . . , n

:= a

a 0

ij i j

 := T

a

 0 j

i 0

Sia scambio colonne j, j :

T := a

ij

 ∀j := 1, . . . , n

:= a

a 0

ij ij

a := T

 0

ij

t := s

j

 := s

s 0

j j

s := t

 0

j

Sia pivot riga k : | |a |

|a := max

Sia i tale che ik

k i:=k,...,n

scambio righe i, k

colonna k :

Sia pivot |a | |a |

Sia j tale che := max kj

kj j:=k,...,n

colonne j, k

scambio

Sia pivot rigcol k :

pivot riga k

colonna k

pivot

Sia pivot scalato k :

|a |

ik ∀i = 1, . . . , n

d :=

i |a |

max j=k,...,n ij

Sia = max d

i tale che d i

i i=1,...,n

scambio riga i, k

Sia ripristina colonne :

( 6

Se i = s :

i ∀i = 1, . . . , n

scambio colonna i, s

i

21

L’algoritmo di eliminazione di Gauss generale, e’ il seguente:

s := (1, 2, . . . , n)

∀k −

= 1, . . . , n 1

pivot k

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

k

a

ik

m := k

a

kk

∀j = k, . . . , n + 1

k+1 k k

a := a ma

ij kj

ij ∨

Se pivot = pivot colonna pivot = pivot rigcol :

ripristina colonne

Nota: pivot puo’ essere una qualsiasi delle procedure di pivot sopra definite.

Infine, l’algoritmo che permette di usare una sola matrice in memoria, e nessuna variabile ausil-

iaria per i moltiplicatori, durante tutto il corso dell’esecuzione e’:

s := (1, 2, . . . , n)

∀k −

= 1, . . . , n 1

pivot k

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

a

ik

a :=

ik a

kk

∀j = k, . . . , n + 1

a := a a a

ij ij ik kj

colonna pivot = pivot rigcol :

Se pivot = pivot

ripristina colonne k = 0 con

j), permette di scambiare gli a

Proof : Per costruzione. la scelta dei vari pivot (i, kk

altri coefficienti non nulli. Questo avviene tramite lo scambio di righe e/o colonne. Ovviamente,

lo scambio di righe o colonne effettuato da pivot non altera il sistema.

Queste tecniche di pivoting, permettono anche di evitare di incorrere in nella situazione in cui il

moltiplicatore m diventa molto grande. In questo caso nella differenza

a := a ma (1)

ij ij kj

si avrebbe una perdita sostanziale di cifre significative.

rigcol, pivot scalato permettono anche di evitare la situazione in cui nella (1), a

I pivot pivot ij

e’ molto piu’ piccolo rispetto a a .

kj

Infine, l’algoritmo ottimizzato per la memoria usa la cella della matrice che sara’ sovrascritta da

k+1 n

k · · ·

uno zero per memorizzare il moltiplicatore, e inoltre, sfrutta il fatto che a = a = = a .

ij ij

ij

n×n

∈ |A| 6

Theorem 4.7. Data A , = 0, si ha

R

∃P matrice di permutazione: P A = LU che si ottengono pero’

dove L, U sono le matrici triangolari descritte nel thm [4.5,pg.18],

usando l’algoritmo [4.6,pg.20] con un pivot che non scambia colonne

22

Una matrice P , si dice matrice di permutazione se e’ ottenuta dalla matrice identica I, mediante

scambi di righe (o, equivalentemente, di colonne).

Proof : Ricordiamo che P A e’ quella matrice ottenuta da A effettuando lo stesso scambio di

righe che sono state utilizzate per ottenere P da I (vedi [AL,1.4.1,pg.1]).

Rileggendo l’algoritmo [4.6,pg.20], possiamo esprimere lo scambio di righe (o di colonne) di A

. I passi compiuti dall’algoritmo si possono

mediante una opportuna matrice di permutazione P i

allora esprimere cosi’:

M P Ax = M P b

1 1 1 1

M P M P Ax = M P M P b

2 2 1 1 2 2 1 1

...

M P . . . M P M P Ax = M P M P M P b

n−1 n−1 2 2 1 1 n−1 n−1 2 2 1 1

Considerando la matrice M P . . . M P M P A, e ponendo P = P P . . . P si di-

n−1 n−1 2 2 1 1 n−1 n−2 1

mostra che P A = LU

4.1.1 Applicazioni n×n

∈ |A| 6

Proposition 4.8. Data la matrice A , = 0, e la matrice U ottenuta dalla fattoriz-

R

zazione [4.7,pg.22], si ha n

Y

s

det A = (−1) u ii

i=1

dove s e’ il numero di scambi di righe effettuati durante la fattorizzazione [4.7,pg.22].

Il costo di questo calcolo e’ quindi 3

n 2

− ≈

costo di [4.7,pg.22] + n 1 + O(n )

3

Proof : Si ha det P A = det P det A = det P det A n

Y

det P A = det LU = det L det U = det U = u ii

|{z} i=1

U e’ triangolare

n

Y

det P det A = u ii

i=1

Ricordando che uno scambio di riga cambia di segno il determinante, si ha

s

det P = (−1) n

Y

s

det A = (−1) u ii

i=1

dove s e’ il numero di righe che sono state scambiate per ottenere P da I. Ma dato che

P = P . . . P , s si puo’ derivare direttamente considerando il numero di scambi di righe

n−1 1

effettuati durante l’uso dell’algoritmo [4.7,pg.22].

23

n×n

∈ |A| 6

Proposition 4.9. Sia A , = 0. Siano dati m sistemi lineari del tipo:

R Ax = b k = 1, . . . , m

k k

Per risolverli si puo’ procedere cosi’:

Si fattorizza A in P A = LU

∀k = 1, . . . , m

y := U x

k

b = P b

k k

Si risolve con [4.1,pg.12] il sistema triangolare Ly = b k

Si risolve con [4.1,pg.12] il sistema triangolare U x = y

k

Il costo di questo algoritmo e’ 3

n 2

+ mn + O(n)

3

Proof : Fissiamo k Ax = b

k k

P Ax = P b

k k b

LU x = P b = k

k k

Sia y = U x

k

Ly = b k 2

n

Il costo dell’algoritmo [4.1,pg.12] usato per risolvere il sistema triangolare e’ + O(n). Poiche’

2

2

stiamo risolvendo due sistemi triangolari, abbiamo n + O(n).

−1

n×n

∈ |A| 6

Proposition 4.10. Sia A , = 0. Sia C = A .

R

Per calcolare C si puo’ procedere cosi’:

∀i = 1, . . . , n i

Si risolve con [4.9,pg.23] il sistema AC = e

i

i

dove C e’ la colonna i-esima di C.

Il costo e’ 4 3

n + O(n)

3

−1 i i

Proof : AC = AA = I, e poiche’ AC = I , si deve avere

i ∀i

AC = e = 1, . . . , n

i

con m = n, si ricava il costo dell’algoritmo.

Applicando la [4.9,pg.23],

4.2 Fattorizzazione di Cholesky

n×n

Theorem 4.11. Sia A una matrice simmetrica e definita positiva, allora

R n×n t

∃R ∈ triangolare inferiore: A = R R

R 24

L’algoritmo per ricavare R e’: ∀j = 1, . . . , n

v j−1

u

u X 2

r := a r

t

jj jj jk

k=1

∀i = j + 1, . . . , n j−1

P

a r r

ij ik jk

k=1

r :=

ij r jj

t

Proof : La condizione A = R R si traduce nel sistema

n

X t ∀j ∀i ⇔

a = R ( R) = 1, . . . , n, = 1, . . . , n

ij ik kj

k=1 n

X

⇔ ∀j ∀i

a = R R = 1, . . . , n, = 1, . . . , n

ij ik jk

k=1

dove r sono gli elementi di R e le incognite del sistema.

hk

Tenendo conto del fatto che R e’ una matrice triangolare inferiore, possiamo scrivere

min{i,j}

X t ∀j ∀i

a = R ( R) = 1, . . . , n, = 1, . . . , n

ij ik kj

k=1

Poiche’ A e’ simmetrica, a = a , quindi le uniche equazioni che ci interessano sono:

ij ji

min{i,j}

X ∀j ∀i

a = r r = 1, . . . , i, = 1, . . . , n

ij ik jk

k=1

Adesso basta scrivere le equazioni nel seguente ordine:

2

1. a = r

11 11

 ∀i

2. a = r r = 2, . . . , n

 i1 i1 11

 2

2 + r

3. a = r

22 21 22 ∀i

4. a = r r + r r = 3, . . . , n

 i2 i1 21 i2 22

..

.

e ricavare le equazioni dall’alto verso il basso. Per ricordare l’ordine delle equazioni si puo’ usare

questo disegno:  

1.

2. 3. 

 

2. 4. 5.

 

2. 4. 6.

Le incognite sono: √

=

r a

11 11

a

i1 ∀i = 2, . . . , n

r =

i1 r 11

q 2

r = a r

22 22 21

a r r

i2 i1 21 ∀i

r = = 3, . . . , n

i2 r 22

... 25

Le radici quadrate sono lecite perche’ A e’ definita positiva.

In generale, v

j j−1 j−1

u

u

X X X

2 2 2 2

⇔ −

a = r = r + r r = a r

t

jj jj jj

jk jj jk jk

k=1 k=1 k=1

j j−1 j−1

P

a r r

ij ik jk

X X k=1 ∀i

⇔ = j + 1, . . . , n

a = r r = r r + r r = r =

ij ik jk ij jj ik jk ij r jj

k=1 k=1

L’algoritmo completo e’ quindi

∀j = 1, . . . , n

v j−1

u

u X 2

r := a r

t

jj jj jk

k=1

∀i = j + 1, . . . , n j−1

P

a r r

ij ik jk

k=1

r :=

ij r jj

4.3 Autovalori m×n

∈ ≥

Theorem 4.12. Data A , con m n,

R

m×m m×n n×n t

∃U ∈ ∈ ∈

, Σ V ortogonali, tali che A = U Σ V

R R R

 σ 1 σ 2 

 

σ 3

 

. . .

 

 σ

Σ= n

 

 

0 0 0 ... 0 

 

0 0 0 ... 0 

 

..

 

. 

 0 0 0 ... 0

≥ ≥ · · · ≥ ≥

σ σ σ 0

1 2 n

I σ si chiamano valori singolari.

i

Proposition 4.13. Nel contesto del thm [4.12,pg.26], se r e’ il rango di A, allora

r

X it i

A = σ U V

i

i=1

i i

dove U e’ la colonna i-esima di U , V la colonna i-esima di V . Inoltre,

kAk

σ =

1 2

Proof : Per il thm [4.12,pg.26] n

X

t t i

A = U Σ V = σ U V

i i

i=1

26

Si ha 6 ⇒ 6

σ = 0 σ , . . . , σ = 0

h 1 h

|{z}

≥σ ≥···≥σ ≥0

σ

1 2 n n×n

Theorem 4.14. di Gershgorin-Hadamard. Sia A C , e siano

{z ∈ | |z − | ≤ }

D a R [e’ un cerchio in di centro a , raggio R ]

C C

i ii i ii i

X |a |

R =

i ij

j6 =

i

n

[ D

D = i

i=1

Si ha: ⇒ ∈

λ autovalore di A λ D

n×1

Proof : Sia u l’autovettore di A, associato a λ. Si ha

C ⇒ 6

u autovettore u = 0

k k

X X X

⇔ ⇔ ⇔ −

Au = λu a u = λu a u + a u = λu a u = (λ a )u

ij j i ii i ij j i ij j ii i

j6 =

i j6 =

i

X

| | |(λ − | |λ − ||u |

a u = a )u = a

ij j ii i ii i

j6 =

i X X

|λ − ||u | | | ≤ |a ||u |

a = a u (1)

ii i ij j ij j

j6 =

i j6 =

i

|u | {|u |}

Sia k : = max . Ponendo i = k nella (1), si ha

k i i X X

X |a ||u | ≤ |a ||u | |u | |a |

|λ − ||u | ≤ =

a kj j kj k k kj

kk k |{z} j6 =

k j6 =

k

j6 =

k |u |=max{|u |}

i

k

X X

|λ − ||u | ≤ |u | |a | ⇒ |a | ⇒ ∈ ⊆D

|λ − | ≤

a = R λ D

a

kk k k kj kj k k

kk

|{z}

j6 =

k j6 =

k

6

u =0

k

Proposition 4.15. Metodo delle potenze per il calcolo dell’autovalore massimo.

n×n

, . . . , λ tutti gli autovalori di A , matrice diagonalizzabile. Sia

Siano λ C

1 n n,1

∈ 6

x , x = 0

C

0 0

( x

y = k

k kx k ∀k = 0, 1, 2, . . .

k

x = Ay

k+1 k

t

x Ax

k k

σ =

k t

x x

k k

t

xAx

σ(x) = e’ chiamato rapporto di Rayleigh. Supponiamo che

t

xx |λ | |λ | ≥ |λ | ≥ · · · ≥ |λ |

>

1 2 3 n

27

Allora si ha: t

y x

k k+1

σ =

k t

y y

k k

lim σ = λ

k 1

k→∞

Supponendo che |λ | |λ | |λ | · · · |λ | |λ | ≥ |λ | ≥ · · · ≥ |λ |

> > > > >

1 2 3 h h+1 h+2 n

per ricavare λ poniamo

h h−1

X

0 −

x = x α u

0 i i

0 i=1

,dove α , u sono descritti nella dimostrazione, e riapplichiamo lo stesso algoritmo. Il nuovo σ

i i k

convergera’ a λ .

h

Infine, scelto un errore ε > 0, possiamo arrestare la successione quando

k(A − k

σ )y

k k <ε

ky k

σ k k

Proof : Sia z = x , z = Az , si ha

0 0 k+1 k n×n

⇒ ∃B {u }

A diagonalizzabile = , . . . , u base di : Au = λ u (1)

C

1 n i i i

n

X

B base x = α u

0 i i

i=1

n n

X X

z = Az = α Au = α λ u

1 0 i i i i i

|{z}

i=1 i=1

(1)

n n

X X 2

z = Az = α λ Au = α λ u

2 1 i i i i i

i

|{z}

i=1 i=1

(1)

... n

X ki

z = Az = α λ u (1.1)

k k−1 i i

i=1 6

Con [v] indichiamo la componente i-esima del vettore v secondo la base canonica. Se α u = 0,

i 1 1l

si ha k+1

λ

n

P

α u + α u

i

n k+1 k+1

P 1 1l i il

α λ u

[z ] λ k+1

i=2 λ

i il

k+1 l i 1

i=1 1 =

= = ki

ki k

P

[z ] λ

α λ u λ n

P

α u + α u

k l i il 1 1 1l i il

k

i=2 λ

1

k+1

n λ

P

α u + α u

i

i il

1 1l i=2 λ

1

= λ (2)

1 k

n λ

P

α u + α u

i

1 1l i il

i=2 λ

1

Poiche’ λ > λ si ha che

1 i k

λ i

lim =0

λ

k→∞ 1

28

e quindi dalla (2) lim z = λ

k 1

k→∞

Per evitare di dover scegliere ogni volta una componente opportuna, si puo’ usare σ ,

k

!

n n

n 2k

|λ |

k i

X

X

X 2 2

τ ki 2 2k 2k |α | |α |

|α | |λ | |λ |

α +

z z = λ = = (3)

α λ

i i i 1 i

k k i i 1 |λ |

1

|{z} i=1 i=2

i=1

(1.1) n

n k X

X k+1

τ τ 2 2 2k

⇒ |α | |α | |λ |

Az = z z Az = z z = λ λ = λ =

k k+1 k k k k+1 i i i 1

i i

|{z} i=1

i=1

(1.1) !

n 2k

|λ | 1

i

X

2k 2 2

|λ | |α | |α |

= λ + (3.1)

1 1 1 i |λ | λ

1 1

i=2 2k

|λ |

n 1

2 2

P

|α | |α | i

+

τ 1 i

z Az i=2 |λ | λ

k k 1 1 0

0 ⇒ lim σ = λ

σ = = λ 1

1 k

k 2k

τ

z z k→∞

|λ |

n

k k |{z} P

2 2

|α | |α | i

+

(3),(3.1) 1 i

i=2 |λ |

1

dove con α abbiamo indicato il complesso coniugato di α .

i i

Poiche’ !

n k

λ i

X

k ⇒

α u + α

[z ] = λ u

1 1l i

k l il

1 λ 1

i=2

 |λ |

+∞ > 1

1

⇒ |[z | |

lim ] = u λ = 1

k l i il 1

k→∞  ≤ |λ |

0 0 < 1

 1

nell’implementazione macchina rischiamo di incorrere in errori di overflow o di underflow. Per

questo motivo si usano gli x , y della tesi, in modo tale da avere gli y della successione con

k k k

ky k = 1. I calcoli svolti fin adesso non variano, infatti, si vede facilmente che

k z k

y =

k kz k

k z k+1

x = Ay =

k+1 k k

kz k t

z z

k+1 k+1

A

t t

x Ax z Az

kz k kz k

k+1 k+1 k+1 k+1 0

k k

= = σ

σ = =

k+1 k+1

t

t t

z z

x x z z

k+1 k+1

k+1 k+1 k+1 k+1

kz k kz k

k k

τ

y x

h1i1. k k+1

Dim. = σ k

τ

y y

k k

Infine abbiamo τ

x x

k

τ τ

x x

y x k+1

kxk k k+1

k k+1 k

= = =

τ x

x x

τ τ

y y x k

k k

k k k kxk

kxk kxk k

k k

τ

x

τ

x A k

τ τ

x Ax

x Ay k kxk k k

k k k

= = = σ

= k

x x

τ τ τ

x x x x

k k

k k k k

kxk kxk

k k

h1i2. Calcolo di λ h

Abbiamo n

h−1

X X

0 − α u

x = x α u =

i i i i

0

0 i=1 i=h

In pratica, da x abbiamo eliminato tutti i termini precedenti a α u . Adesso, il λ dominante

0 h h i

e’ λ . Ripetendo in modo analogo i ragionamenti di prima, si ha la tesi.

h

h1i3. Criterio di arresto 29

Il criterio di arresto, si basa sulle seguenti osservazioni:

⇔ −

Au = λ u (A λ I)u = 0

1 1 1 1 1

lim y = u

k 1

k→∞

lim σ = λ

k 1

k→∞ −

lim (A σ I)y = 0

k k

k→∞ |{z}

vettore nullo

k(A − k

lim σ I)y = 0

k k

k→∞

Proposition 4.16. Metodo delle potenze inverse per il calcolo dell’autovalore minimo.

n×n

, . . . , λ tutti gli autovalori di A , matrice diagonalizzabile. Sia

Siano λ C

1 n n,1

∈ 6

x , x = 0

C

0 0

( x

y = k

k kx k ∀k = 0, 1, 2, . . .

k −1

x = A y

k+1 k

t

x Ax

k k

σ =

k t

x x

k k

Supponiamo che −1

∃A

|λ | ≥ |λ | ≥ · · · ≥ |λ | |λ |

>

1 2 n−1 n

Allora si ha: 1

lim σ =

k λ

k→∞ n

−1

Al posto di calcolare A , per poi calcolare i termini della successione:

−1

x = A y

k+1 k

conviene calcolare i termini risolvendo di volta in volta il sistema

Ax = y

k+1 k

nell’incognita x .

k+1

Poiche’ i singoli sistemi hanno la stessa matrice dei coefficienti, si puo’ usare la prop. [4.9,pg.23].

Proof : Basta osservare che 1 −1

λ autovalore di A autovalore di A

i λ i

1 1 1

1 ≤ ≤ ··· ≤ <

|λ | |λ | |λ | |λ |

1 2 n−1 n

1 1

−1 −1

λ (A ) = λ (A ) =

max max

(A)

λ λ

min n

(vedi la dimostrazione di [3.10,pg.11]) e applicare la prop. [4.15,pg.27].

30 ∈

Proposition 4.17. Nelle ipotesi della prop [4.16,pg.30], preso un µ e ponendo

C

0

A = A µ I

µ k

k

(k)

σ approssimazione di σ di A calcolato con la prop. [4.16,pg.30], con errore ε

i µ

k

1

µ = µ +

k+1 k (k)

σ

si ha lim µ = λ

k µ

k→∞

e’ l’autovalore di A piu’ vicino a µ, ovvero

dove λ µ − ≤ |λ − ∀i

|λ µ| µ| = 1, . . . , n

µ i

Proof : Si ha ⇔ ⇔ − − ⇔ − − ⇔ − ⇔ −

λ autovalore di A Au = λu Au µu = λu µu (A µI)u = (λ µ)u A u = (λ µ)u λ µ autovalor

µ

⇒ {λ − |

Gli autovalori di A sono µ i = 1, . . . , n} (1)

µ i

Quindi applicando la prop [4.16,pg.30] a A , abbiamo che

µ 1

1

(0) ≈ =

σ −

λ (A ) λ µ

min µ t

|{z}

(1)

|λ − ≤ |λ − ∀i

dove µ| µ| = 1, . . . , n

t i

1 1 1

(0) ≈ ⇔ ≈ − ⇔ ≈

σ λ µ + µ λ

t t

(0) (0)

λ µ σ σ

t

1 + µ, e reiterando, otteniamo

Ponendo µ =

1 (0)

σ 1 ≈

+ µ λ

1 t

(1)

σ

che e’ una approssimazione migliore di λ .

t

5 Interpolazione e approssimazione 2

∈ ∀i

Definition 5.1. Supponiamo di avere gli n+1 punti (x , y ) = 0, 1, . . . , n. L’interpolazione

R

i i

−→

consiste nel trovare una funzione f : tale che

R R ∀i

f (x ) = y = 0, 1, . . . , n

i i

Interpolazione polinomiale Si richiede che f (x) sia un polinomio in di grado minimo.

R[x]

Interpolazione razionale Si richiede che f (x) sia una funzione razionale, ovvero che si del tipo

g(x) ∈

, g(x), h(x) R[x]

h(x)

Interpolazione trigonometrica Si richiede che f (x) sia una funzione di periodo L. Si ricerca

una soluzione tra le funzioni del tipo

n

a 2π 2π

0 X

+ b sin k x + a cos x

k

k k

2 L L

k=1 31

5.1 Interpolazione polinomiale n

· · ·

Proposition 5.2. Sia f (x) = a + a x + + a x . Un algoritmo per computare f (x) in

0 1 n

maniera efficace e’ il seguente: t := a

n

∀i = 1, . . . , n

t := a + tx

n−i

Il costo dell’algoritmo e’ quindi n moltiplicazioni, n addizioni

Proposition 5.3.

∃!p(x) ∈ ∀i ⇔ 6 ∀i

: p(x ) = y = 0, 1, . . . , n x = x > j

R[x]

n i i i j

dove e’ lo spazio vettoriale dei polinomi di grado n a coefficienti reali.

R[x]

n ≤

Proof : Sia p(x) il generico polinomio di grado n incognita:

n

X j

p(x) = a x

j

j=0

n

X ji

∀i ⇔ ∀i

p(x ) = y = 0, . . . , n a x = y = 0, . . . , n

i i j i

j=0

Quest’ultimo e’ un sistema lineare nelle incognite a , . . . , a , in n + 1 equazioni. La matrice dei

0 n

coefficienti e’ 10 n

0 

 x . . . x

x

0 0

01 11 n

x x . . . x 

 1

V = 

 .

.. 

 

 0 1 nn

x x . . . x

n n

Il suo determinante e’ il noto determinante di Vandermonde:

−x

x

i j

Y

|V | = i>j

⇔ |V | 6

Il sistema ammette un unica soluzione = 0, ovvero

|V | 6 ⇔ − 6 ∀i

= 0 x x = 0 > j

i j

Theorem 5.4. Forma di Lagrange.

Sia −

x x

i

Y

λ (x) =

j −

x x

j i

0≤i≤n, i6 =

j

allora la soluzione dell’interpolazione e’ n

X

p(x) = y λ (x)

j j

j=0

32

Proof : Nella dimostrazione [5.3,pg.32], abbiamo visto che trovare p(x) equivale a risolvere il

sistema Va = y  

 

a y

0 0

a y

1 1

   

con a = , y =

   

. .

.. ..

 

 

   

a y n

n

Poiche’ p(x) e e’ uno spazio vettoriale, possiamo considerare una

Proof sketch: R[x] R[x]

n 0

{λ (x), λ (x), . . . , λ (x)} in modo tale che la nuova matrice dei coefficienti V sia I e che

base 0 1 n

deg λ (x) n Ovvero:

j n

X 0

p(x) = a λ (x) (1)

j

j

j=0

n

X 0 ∀i

p(x ) = a λ (x ) = 0, 1, . . . , n (1.1)

i j i

j

j=0

0

(V ) = λ (x )

ij j i

⇔ ⇔

V = I (V ) = δ λ (x ) = δ

ij ij j i ij

Ricaviamoci i λ (x):

j ∀i ⇒ ∀i 6

λ (x ) = δ λ (x ) = 0 = j

j i ij j i

( ∀i −

λ (x ) = 0 = 0, 1, . . . , j 1, j + 1, . . . , n

j i Y −

⇒ x x

λ (x) = c i

j

deg λ (x) n |{z}

j i6 =

j

Ruffini

1

Y

⇔ − ⇔

λ (x ) = 1 c x x = 1 c =

j j j i Q −

x x

j i

i6 =

j

i6 =

j

Q −

x x −

x x

i i

i6 =

j Y

⇒ λ (x) = =

j Q − −

x x x x

j i j i

i6 =

j 0≤i≤n, i6 =

j

Quindi  0

a = y 0

0

 0

n a = y

 1

 1

X 0

⇔ ∀i ⇔

p(x ) = y p(x ) = a δ = y = 0, 1, . . . , n .

i i i ij i

j ..

|{z} 

j=0 

(1.1) 

 0

a = y

 n

n

n n

X X

0

p(x) = a (x) = y λ (x)

λ j j j

j

|{z} j=0 j=0

(1) 0

∈ ∈ ∀i

Definition 5.5. Sia f C ([a, b]), x [a, b] = 0, 1, . . . , n definiamo la differenza divisa:

i

f [x ] = f (x ) ovvero la funzione f calcolata in x

0 0 0

f [x , . . . , x ] f [x , . . . , x ]

1 j 0 j−1

f [x , . . . , x ] = , per j > 0

0 j −

x x

j 0

33 n(n+1)

Proposition 5.6. Il seguente e’ un algoritmo che calcola f [x , . . . , x ] in flops:

0 n 2

∀i

a := f (x ) = 0, . . . , n

i i

∀k −

= n 1, . . . , 1, 0

∀i = 0, . . . , k −

a a

i+1 i

a :=

i −

x x

i+n−k i

Il risultato e’ in a

0

si puo’ visualizzare l’algoritmo cosi’:

f [x ] = f (x )

0 0 f [x , x ]

0 1

f [x ] = f (x ) f [x , x , x ]

1 1 0 1 2

f [x , x ] f [x , x , x , x ]

1 2 0 1 2 3

f [x ] = f (x ) f [x , x , x ]

2 2 1 2 3

f [x , x ]

2 3

f [x ] = f (x )

3 3

0

∈ ∈ ∀i

Theorem 5.7. Sia f C ([a, b]), x [a, b] = 0, 1, . . . , n. Il polinomio di Newton che

i

interpola i punti (x , f (x )) e’

i i n

X

p (x) = f [x , . . . , x ]w (x)

n 0 i i

i=0 ∀i

p (x ) = f (x ) = 0, 1, . . . , n dove

n i i

w (x) = 1

0 j−1

Y − ∀j

w (x) = (x x ) > 0

j i

i=0

p (x) ha il vantaggio di poter essere facilmente aggiornato per l’aggiunta di nuovi punti. Infatti,

n

dato x , basta calcolare

n+1 p (x) = p (x) + f [x , . . . , x ]w (x)

n+1 n 0 n+1 n+1

Inoltre, vale la seguente proposizione: ∀x ∈

f (x) = p (x) + R (x) [a, b]

n n

dove R (x) = f [x , . . . , x , x]w (x)

n 0 n n+1

Infine, per calcolare p (x) in x, si puo’ usare il seguente algoritmo che si basa sulla prop

n

[5.2,pg.32]: t := f [x , . . . , x ]

0 n

∀i = 1, . . . , n · −

t := f [x , . . . , x ] + t (x x )

0 n−i n−i 5

Dove f [x , . . . , x ] si puo’ calcolare con la prop [5.6,pg.33] (nota )

0 n−i

5 in realta’, non serve calcolare a ogni passo f [x , . . . , x ] usando [5.6,pg.33], perche’ alcuni calcoli vengono

0 n−i

ripetuti inutilmente. Si puo’ quindi pensare a un algoritmo piu’ ottimizzato.

34

Proof : Poniamo p (x) = f (x ), p (x) = 0

−1

0 0

h1i1. Dim. che k

X

∃A ∈ ∀k

, . . . , A : p (x) = p (x) + A w (x) = A w (x) = 0, . . . , n

R

1 n k k−1 k k i i

i=0

∀i

dove p (x) polinomio interpolatore di (x , f (x )) = 0, . . . , k.

k i i

Per Hp abbiamo: ∀i

(x ) = f (x ) = 0, 1, . . . , n

p n i i ≤

Usando un polinomio incognita g(x) di grado n + 1:

(x) = p (x) + g(x)

p n+1 n ∀i ⇒ ∀i

p (x ) = f (x ) = p (x ) = 0, . . . , n g(x ) = 0 = 0, . . . , n (1)

n+1 i i n i i

( ≤

deg p (x) n + 1

n+1 ⇒ ≤

deg g(x) n + 1 (2)

deg p (x) n

n n

Y

⇒ −

(1), (2) g(x) = A (x x ) = A w (x)

n+1 i n+1 n+1

|{z} i=0

Ruffini

(x) = p (x) + A w (x) (3)

p n+1 n n+1 n+1

In particolare abbiamo: 0

X

p (x) = f (x ) = A w (x)

0 0 i i

i=0 1

X

(x) = f (x ) + A (x x ) = A w (x)

p 1 0 1 0 i i

i=0 2

X

− − −

p (x) = f (x ) + A (x x ) + A (x x )(x x ) = A w (x)

2 0 1 0 2 0 1 i i

i=0

... n

X

p (x) = A w (x)

n i i

i=0

Adesso resta da determinare gli A .

i

h1i2. Dim. che vale la seguente implicazione: n

X

∀i ⇒

p (x) e’ il polinomio interpolatore di (x , f (x )) = 0, . . . , n p (x) = f [x , . . . , x ]w (x)

n i i n 0 i i

i=0

Procediamo per induzione su n.

h2i1. n = 0 (x) = f (x ) = f [x ]

p 0 0 0

h2i2. n = 1 − −

p (x) = p (x) + A (x x ) = f [x ] + A (x x )

1 0 1 0 0 1 0

|{z} |{z}

h1i1 h2i1

passo passo

Imponendo che p (x ) = f (x ), abbiamo

1 1 1 −

f (x ) f (x )

1 0 = F [x , x ]

A = 0 1

1 −

x x

1 0

h2i3. − ⇒

n 1 n 35

Sia h (x) il polinomio d’interpolazione relativo ai punti x , . . . , x , e g (x) quello rela-

n−1 1 n n−1

tivo ai punti x , . . . , x , ovvero

0 n−1 ∀i

h (x ) = f (x ) = 1, . . . , n (4)

n−1 i i ∀i −

g (x ) = f (x ) = 0, . . . , n 1 (4.1)

n−1 i i

Si ha −

− x x

x x n

0 h (x) + g (x) (5)

p (x) = n−1 n−1

n − −

x x x x

n 0 n 0

infatti, − −

x x x x

0 0 n 0

x = x h (x ) + g (x ) = g (x ) = f (x ) = p (x )

0 n−1 0 n−1 0 n−1 0 0 n 0

− −

x x x x

n 0 n 0 |{z}

(4.1)

− −

x x x x

n 0 n n

x = x h (x ) + g (x ) = h (x ) = f (x ) = p (x )

n n−1 n n−1 n n−1 n n n n

− −

x x x x

n 0 n 0 |{z}

(4)

− −

x x x x

i 0 n i

x = x , 0 < i < n h (x ) + g (x ) =

i n−1 i n−1 i

− −

x x x x

n 0 n 0

− −

x x x x

i 0 n i

+ f (x ) = f (x ) = p (x )

= i i n i

− −

x x x x

n 0 n 0

|{z}

h (x )=f (x )=g (x )

n−1 i i n−1 i si ha l’uguaglianza polinomiale nella (5).

e quindi per l’unicita’ di p (x) (vedi [5.3,pg.32])

n 6

Applicando l’ipotesi induttiva su h e g abbiamo

n−1 n−1

− − − · · ·

h (x) = f [x ] + f [x , x ](x x ) + f [x , x , x ](x x )(x x ) + +

n−1 1 1 2 1 1 2 3 1 2

− − −

+ f [x , x , . . . , x ](x x )(x x ) . . . (x x )

1 2 n 1 2 n−1

− − − · · ·

g (x) = f [x ] + f [x , x ](x x ) + f [x , x , x ](x x )(x x ) + +

n−1 1 1 2 0 0 1 2 0 1

− − −

+ f [x , x , . . . , x ](x x )(x x ) . . . (x x )

0 1 n−1 0 1 n−2

Indicando con M (h ) il coefficiente dell’indeterminata grado massimo di h , abbiamo

n−1 n−1

M (h ) = f [x , x , . . . , x ]

n−1 1 2 n

M (g ) = f [x , x , . . . , x ]

n−1 0 1 n−1 −

− f [x , x , . . . , x ] f [x , x , . . . , x ]

M (h ) M (g ) 1 2 n 0 1 n−1

n−1 n−1 = = f [x , . . . , x ] (5.1)

M (p ) = 0 n

n − −

x x x x

n 0 n 0 |{z}

|{z} per def.

(5)

D’altra parte,

p (x) = p (x) + A w (x)

n n−1 n n

 |{z}

 (3)

 ≤ − ⇒

deg p (x) n 1 M (p ) = M (A w (x)) = A = f [x , . . . , x ]

n−1 n n n n 0 n

|{z}

 ≤

w (x) n

deg A

 (5.1)

n n

w (x) e’ monico

n

p (x) = p (x) + f [x , . . . , x ]w (x) =

n n−1 0 n n

n−1 n

X X

= f [x , . . . , x ]w (x) + f [x , . . . , x ]w (x) = f [x , . . . , x ]w (x)

0 i i 0 n n 0 i i

|{z} i=0 i=0

Hp induttiva su p

n−1

h1i3. ∀x ∈

Dim. che f (x) = p (x) + R (x) [a, b]

n n

∈ ≡

Fissiamo x [a, b]. Se consideriamo (x , f (x )) (x, f (x)) come un nuovo punto di

n+1 n+1

6 se volessimo applicare l’Hp induttiva formalmente, ad esempio su h , basterebbe porre y = x , y =

n−1 0 1 1

x , . . . , y = x n

2 n−1 36

interpolazione, abbiamo

p (x) = p (x) + f [x , . . . , x ]w (x) = p (x) + R (x)

n+1 n 1 n+1 n+1 n n

|{z}

h1i2

passo

Proposition 5.8. Nel contesto di [5.7,pg.34], si ha (n+1)

f (ξ)

(n+1)

∃f ⇒ ∃ξ ∈

(x) in [a, b] [a, b] : R (x) = ω (x)

n n+1

(n + 1)!

(n+1)

dove f (x) e’ la derivata n + 1-esima di f .

Proposition 5.9. Sia −→

t (x) = cos nϕ(x) : [−1, 1] R

n −→

ϕ(x) = arcos x : [−1, 1] [0, π]

t (x) si puo’ esprimere come un polinomio, chiamato polinomio di Chebyshev:

n t (x) = 1

0

t (x) = x

1 − ∀n

t (x) = 2t (x)x t (x) > 1

n n−1 n−2

Valgono le seguenti proprieta’:

∀n ≥

deg t (x) = n 0

n n−1 ∀n ≥

c = 2 coeff. di grado massimo di t (x) 1

n n

π

π

+ k k = 0, . . . , n 1 sono tutte le radici di t (x)

α = cos n

k 2n n

2π k = 0, . . . , n sono tutti i punti di massimo

M = cos k

k n

π

= cos + k = 0, . . . , n sono tutti i punti di minimo

m k

k n n

Proof : −

t (x) = cos((n + 1)ϕ(x)) = cos(nϕ(x) + ϕ(x)) = cos(nϕ(x)) cos(ϕ(x)) sin(nϕ(x)) sin(ϕ(x))

n+1 − −

t (x) = cos((n 1)ϕ(x)) = cos(nϕ(x) ϕ(x)) = cos(nϕ(x)) cos(ϕ(x)) + sin(nϕ(x)) sin(ϕ(x))

n−1

Sommando membro a membro: ⇔ −

t (x) + t (x) = 2 cos(nϕ(x)) cos(ϕ(x)) = 2t (x)x t (x) = 2t (x)x t (x)

n+1 n−1 n n+1 n n−1

t (x) = 2t (x)x t (x)

n n−1 n−2

h1i1. deg t (x) = n

n

Per induzione forte su n: I casi n = 0, 1 sono veri.

deg t (x) = deg [2t (x)x t (x)]

n n−1 n−2 −

deg t (x)x = deg t (x) + 1 = (n 1) + 1 = n

n−1 n−1 |{z}

Hp induttiva

(x) = n 2

deg t n−2 |{z}

Hp induttiva

⇒ deg t (x) = n

n 37

n−1

h1i2. ≥

c = 2 , n 1

n 1−1

t (x) = x c = 1 = 2

1 1

 −

t (x) = 2t (x)x t (x)

n n−1 n−2

 n−2 n−1

⇒ c(t ) = c(2t (x)x) = 2c(t (x)x) = 2c(t ) = 22 = 2

(x)x = n

deg 2t n n−1 n−1 n−1

n−1 |{z}

 −

deg t = n 2 Hp induttiva

 n−2

dove con c(t ) abbiamo indicato il termine di grado massimo di t (x).

n n

Infine, per ricercare le radici, i massimi e i minimi, basta cercare i punti in cui la funzione cos nϕ

−1,

assume rispettivamente i valori 0, 1, e ricavarsi la x sapendo che x = cos(ϕ).

0

Proposition 5.10. Supponiamo di voler interpolare la funzione f (x) C ([a, b]) in n + 1

(x) il polinomio interpolatore nei punti

punti, che possiamo scegliere arbitrariamente. Sia p n

(x , f (x )), . . . , (x , f (x ))

0 0 n n |f −

Se [a, b] = [−1, 1], allora la scelta che minimizza l’errore (x) p (x)| e’

n

∀k

x = α = 0, . . . , n

k k

dove α sono le radici del polinomio di Chebyshev vedi [5.9,pg.37].

k

Example 5.11. Fenomeno di Runge.

Consideriamo la funzione 1 −→

: [−1, 1]

f (x) = R

2

1 + 25x

Se la interpoliamo in n + 1 punti distribuiti uniformemente in [−1, 1], otteniamo una pessima

∈ |f −

approssimazione di f (x), ovvero esistono alcuni punti x [−1, 1] dove (x) p (x)| e’ molto

n

→ ∞ ∈

alto. Al crescere di n la situazione peggiora: si dimostra che per n esistono degli x [−1, 1]

|f − → ∞.

t.c. (x) p (x)| Vedi figura:

n 38

f (x) e’ la linea in rosso.

Se invece scegliamo le radici di Chebyshev come punti per l’interpolazione, abbiamo:

La miglioria e’ dovuta al fatto che le radici di Chebyshev sono piu’ distribuite verso glii

estremi dell’intervallo [−1, 1].

I grafici sono stati prodotti con il seguente script octave:

clf

hold on

xf=[-1:0.002:1]; yf=1./(1+25*xf.^2);

plot(xf,yf,’;Funzione originale;’)

nodi=[5:4:16];

steps=2./(nodi-1);

function yi=polyinterp(x,y,xi)

n = length(x);

M = zeros(n);

for i=1:n

M(:,i) = (x.^(i-1))’;

endfor

a = M \ y’;

g = ones(1,length(xi)) * a(1);

for i=2:n

g = g + xi.^(i-1)*a(i);

endfor

yi=g;

endfunction 39


PAGINE

49

PESO

391.07 KB

AUTORE

Exxodus

PUBBLICATO

+1 anno fa


DETTAGLI
Corso di laurea: Corso di laurea in informatica
SSD:
Università: Catania - Unict
A.A.: 2009-2010

I contenuti di questa pagina costituiscono rielaborazioni personali del Publisher Exxodus di informazioni apprese con la frequenza delle lezioni di Analisi numerica e studio autonomo di eventuali libri di riferimento in preparazione dell'esame finale o della tesi. Non devono intendersi come materiale ufficiale dell'università Catania - Unict o del prof Russo Giovanni.

Acquista con carta o conto PayPal

Scarica il file tutte le volte che vuoi

Paga con un conto PayPal per usufruire della garanzia Soddisfatto o rimborsato

Recensioni
Ti è piaciuto questo appunto? Valutalo!

Altri appunti di Corso di laurea in informatica

Calcolo Combinatorio
Appunto
Tesi Marco Zennaro
Tesi
Matrici e sistemi lineari teoria ed esercizi parte 2
Appunto
Matrici e sistemi lineari teoria ed esercizi parte 1
Appunto