Calcolo_Numerico

8/13/2019 Calcolo_Numerico

1/264

CALCOLO NUMERICOLaurea Specialistica Ingegneria Elettronica, Gestionale, Informatica

Prof. M.L. Lo Cascio

(A.A. 2006-2007)

Lezione del 15.01.2007

Argomenti

Il calcolo numerico

Il sistema dei numeri macchina

Errore di arrotondamento

Errore assoluto, errore relativo

Esempi sugli effetti della propagazione dellerrore

Stabilita di un algoritmo

1

Cose il CALCOLO NUMERICO ?

IlCalcolooAnalisi Numerica e larea della matematica che

- crea, analizza i metodi per fornire una soluzione nume-

rica dei problemi della matematica continua

- studia gli algoritmi, implementamentabili su un calcola-tore, dedotti dai metodi stessi.

2


2/264

Problemi della matematica continua(Punto di partenza): Problemi

matematici in cui sono coinvolte variabili appartenenti a domini del

campo reale R o complesso C (coppie di numeri reali).

Sono i problemi, costruiti dalle scienze applicate, quali lingegneria o

la medicina, formalizzati attraverso modelli matematici nel continuo -

per esempio, integrali, equazioni differenziali, sistemi lineari.

Soluzione numerica (Punto di arrivo): e sempre costituita da un in-

sieme discreto di numeri, che rappresenta la soluzione del problema

ottenuta mediante luso di un calcolatore sul quale si e implementato

lalgoritmo.

3

Schema di procedura

problema matematico nel continuo

scelta del metodo/algoritmo numerico nel discreto.La scelta e un arte!!

implementazione dellalgoritmo in un linguaggio di programmazione

verifica del programma su problemi test

esecuzione del programma sul problema posto

analisi dei risultati !!!

La soluzione numerica e accettabile solo se si sa stimare lattendibilitadei risultati

perche nei vari passaggi si introducono inevitabilmente degli errori

4


3/264

Quali errori?

Errori di arrotondamento: dovuti essenzialmente al fatto che un in un

calcolatore si possono rappresentare solo numeri di lungheza finita

Errori di discretizzazione o troncamento: legati alla trasformazione di

un problema nel continuo in un problema nel discreto

5

Il sistema di numeri macchina

Linsieme F dei numeri disponibile su un calcolatore e detto sistema

di numeri macchina

F e un sottoinsieme di R discreto

F contiene numeri costituiti da un numero finito di cifre [M, m] {0} [m, M], m > 0, M > 0 dipendenti dal tipo di calcolatore e/o dilinguaggio usato.

Analisi MatematicaGeometria, Agebra

R: Numeri reali

insieme infinitonumeri di lunghezza infinita

Analisi Numerica

F: Numeri macchina

insieme finitonumeri di lunghezza finita

6


4/264

Errori di arrotondamento

Esempi

3.1415926535897932385 e 13711.704699910719625109

x R x F10 x F15 3.141592654 3.14159265358979137 11.70469991 11.7046999107196

Definizione

x R valore esatto, x valore approssimato

Errore assoluto e(x) =x x

Errore relativo er(x) =x x

x(x= 0)

7

Tabella degli errori sui valori ,

137

F10 F15v R e(v) er(v) e(v) er(v)

-.41e-9 -.13e-9 .32e-14 .10e-14

137 .72e-9 .61e-10 .25e-13 .21e-14

Lerrore assoluto misura il numero di cifre decimali esatte

Lerrore relativo misura il numero di cifre significative esatte

x e il valore di x arrotondato sulla r-sima cifra decimale e sulla t-simacifra significativa se r e t sono il massimo intero positivo per cui

|e(x)

| 0.5

10r ,

|er(x)

| 0.5

10t+1

8


5/264

Osservazione

- Lerrore assoluto si considera (spesso) dotato di segno:

errore assoluto positivo x approssimazione per difetto

errore assoluto negativo x approssimazione per eccesso

- Generalmente, non essendo disponibile x, si valutera lerrore relativo

con x, disponibile, al posto di x.

er(x) =x x

x (x= 0)

9

Quando si introducono gli errori di arrotondamento?

- sui dati di inputquando / F.

Esempio: (0.1)10= (0.0001100)2,

100k=1

0.1 = 9.99999999999998= 10 (Matlab).

- su ogni risultato che eccede la rappresentabilita di memoria

Esempio: x= e15 = 3269017.372472110639301855....

Matlab fornisce: x= 3269017.372472111

Larrotondamento sui numeri e la prima fonte di errore.

10


6/264

Errori di arrotondamento: alcuni danni

Esempio 1 - Luguaglianza fra numeri reali

q1(x) = (x 1)7 q2(x) =x

7 7x6 + 21x5 35x4 + 35x321x2 + 7x 1

Per lalgebra q1(x) e q2(x) sono identichex.

Tabuliamo q1 e q2 per x

[0.9998, 1.0002] con MATLAB.

11

0.9998 0.9999 1 1.0001 1.00021.5

1

0.5

0

0.5

1

1.5x 10

26

0.9998 0.9999 1 1.0001 1.00021.5

1

0.5

0

0.5

1

1.5x 10

14

: q1 : q2

I valori assunti da q1 e q2 non solo non coincidono, ma hanno ordine

di grandezza sensibilmente differente:

Luguaglianza fra numeri reali puo non avere significato!

12


7/264

Esempio 2 - La cancellazione numerica

Consideriamo lequazione di secondo grado

q(x) =ax2 + bx + c= 0

e la formula classica che fornisce le soluzioni

x1=b+

2a x2=b

2a

=b2 4ac

13

a= 1, b=124.81, c= 1.3432;

Soluzionex1= 124.79923711369504356, x2 = 0.010762886304956441142

n. cifre x1 x2 er(x1) er(x2)

7 124.7992 .01075000 .30e-6 .12e-2

12 124.799237114 .0107628865000 -.24e-11 -.18e-7

15 124.799237113695 .107628863050000e-1 .35e-15 -.40e-11

Le soluzioni dellequazione di secondo grado fornite dalla formula classica non siottengono con la stessa accuratezza

Perche ?b,

esatti nella precisione fissata

|b| nel calcolo di x2 si sottraggono numeri molto vicini si perdono cifresignificative, ovvero viene amplificato lerrore relativo. (Cancellazione numerica)

14


8/264

Bisogna evitare la sottrazione.

Ricordando che x1x2 = c/a un altro algoritmo di calcolo e il seguente

x1=b + sign(b)

2ax2=

c

a x1

La radice di massimo modulo resta inalterata e per la radice x2 si

ottiene

n. cifre x2 er(x1)

7 .01076289 -.34e-6

12 .0107628863049 .52e-11

15 .0107628863049564 .38e-14

15

Stabilita e condizionamento

Al problema della propagazione degli errori di arrotondamento sono

associati il concetto di

stabilita di un algoritmo

condizionamento di un problema

Un algoritmo si dice (numericamente) stabile se loutput subisce solo

piccole variazioni quando i dati di input subiscono piccole variazioni;

si dira instabilese non verifica tale condizione.

16


9/264

Stabilita di un algoritmo: esempio

Modello matematico: In= 1e

1

0

xnex dx

Integrando per parti si ottiene

In= 1e

e

10

nxn1ex dx

= 1 nIn1

Algoritmo di calcolo

I0= 1

1e

In= 1 nIn1, n1 (legge ricorsiva)

17

Osservazioni

- xnex >0 per x(0, 1]In>0n finito

- xn < xn1 per x(0, 1]In< In1

ovvero{In} e unasuccessione monotona decrescente di valori positivi.

Applichiamo la ricorrenza: In= 1 nIn1, n1

I00.63212055882856 =I0 numero macchinaI1= 1 I01 I0 = 0.36787944117144 =I1I2= 1

2I1

1

2I1= 0.36787944117144 =I2

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

18


10/264

n In0 0.632120558828561 0.367879441171442 0.264241117657123 0.20727664702865

9 0.09161229299417

16 0.0554593017295717 0.0571918705973118 -0.02945367075154

qualcosa non va!negativo? impossibile !!!

E proseguendo si ha

I19=1.559619744279I20=30.192394885584I21=635.040292597259

19

Tutta colpa della propagazione dellerrore iniziale!

Algoritmo applicato

In= 1 nIn1 valore iniziale I0, numero di macchina

Algoritmo esatto

In= 1 nIn1 valore iniziale I0 R

Ad ogni passo si commette un errore

en=In In= 1 nIn1 (1 nIn1) =nen1, n1

e1

=

e0

, e2

=

2e1

= 2e0

, e3

=

3e2

=

2

3e0

20


11/264

per induzioneen =(1)nn!

e0

Coeff. di amplificazioneIl valore di innesco, I0, e arrotondato sulla 14a cifra decimale

|e0| 12 1014, ma lerrore iniziale cresce in modo proporzionale

al fattoriale e diventa dominante rispetto ai valori esatti In.

|e20| 20! 12 1014 0.12e5

L algoritmo non e stabile!

21

Un nuovo algoritmo

- la In= 1 nIn1 puo essere letta anche

In1= 1n(1 In) (ricorsione allindietro)

- poiche In

0, per n sufficientemente grande si puo assumere

In1In In1 nIn In 1n+1.

Da I35= 0.02777777777778 applicando la ricorsione allindietro si haI20= 0.04554488407582, che e il valore esatto nella precisione usata.

Verificate che in questo caso lerrore inizialee35viene progressivamente

abbattuto e pertanto

L algoritmo e stabile!

22


12/264



(A.A. 2006-2007)


Argomenti

Condizionamento di un problema

Valutazione di una funzione

Errore di toncamento o discretizzazione

Formule ricorsiveIl metodo del punto unito

Teorema di esistenza del punto unito

Esempi

1

Condizionamento di un problema

In alcuni casi gli effetti devastanti della propagazione dellerrore di

arrotondamento hanno unorigine sotto certi aspetti piu grave, perche

dipendono dal problema stesso.

Rappresentazione schematica di un problema

ID P ODID: dati di input, per esempio la matrice e i termini noti di un sistema

lineare

OD: valori di output, per esempio la soluzione del sistema

Quanto e sensibile loutput in presenza di perturbazioni sullinput, in-

dipendentemente da come si eseguono i calcoli?

Un problema si dice ben condizionato se la soluzione subisce solo pic-

cole variazioni quando i parametri da cui essa dipende subiscono piccole

variazioni; si dira mal condizionato se non verifica tale condizione.

2


13/264

Esempio 1 - Zeri di un polinomio

p(x) =x3 + a1x2 + a2x + a3

a1 a2 a3 1 2 3

-9.5 20.0 12.5 -0.5 5.0 5.0

-9.498 20.001 12.5 -0.5 4.999+0.099995i 4.999-0.099995i

Errore percentuale sui dati (in modulo):

er(a1) = 0.02%, er(a2) = 0.005%, er(a3) = 0

Errore sulla soluzione: non valutabile:

cambia lo spazio cui appartengono gli zeri di p

3

Esempio 2: sistemi lineari

Si considerino i sistemi Ax= b con

I. A= 5 70.7 1 , b= 121.7 , x= 11

II. A=

5 70.7 1.02

, b=

12.02

1.7

, x=

2.720.2

III. A=

5 70.7 1

, b=

121.69

, x=

1.70.5

4


14/264

Errori sui dati e sulle soluzioni

Considerano i Casi II e III come perturbazioni del Caso I.

Caso Errore sui dati Errore sulla soluzione

II. er(a22) = 2% er(x1) = 172%

er(b1) = 0.2% er(x2) = 120%

III. er(b2) = 0.6% er(x1) = 70%

er(x2) = 50%

Gli errori sui risultati non sono proprio comparabili con le perturbazioni

introdotte sui dati!

5

Esempio 3 - Valutazione di una funzione in un punto

p(x) = 19x5 135x4 + 76x3 540x2 + 57x 405 = 0Zero di p (Valore vero) Valore approssimato

= 135/197.105263157894736842 = 7.10526

Corrispondenti valori di p

p() = 0 p() =0.1652204

Quale la relazione fra lerrore relativo sui dati (variabile indipendente)

e lerrore relativo sui risultati (valori della funzione)?

6


15/264

Indice di condizionamento nel calcolo di una y =f(x)

Sia xR, valore esatto; x , valore approssimato

Sviluppo in serie di Taylor:

y =y y= f(x) f(x) =f(x)(x x) + ......

Errore relativo:

yy

f(x)f(x)

|x|=

f(x)x

f(x)

xx

Indice di condizionamento del problema: CP := f(x)xf(x)

7

Osservazioni

- Nella pratica si sostituisce il valore esatto x (in genere non disponibile)

con il valore approssimato xCP :=

f(x)xf(x)

- CP e il fattore di amplificazione dellerrore relativo sui dati: maggiore

e Cp piu il problema e mal-condizionato

- Il calcolo di una funzione in uno zero e un problema mal-condizionato

- Lindice di condizionamento cresce al crescere della derivata

- Se lerrore relativonon e definito, si utilizza lerrore assoluto

Esempio Nel caso del polinomio: p(x)5.23e + 004,

Cp() =p()0.37e06

8


16/264

... concludendo

instabilita , mal-condizionamento: piccole variazioni percentuali sui

dati producono variazioni molto maggiori sui risultati

Linstabilita di un algoritmo e dovuta al tipo di operazioni che si fanno si puo ovviare cambiando algoritmo

Nel mal-condizionamento di un problema lattenzione e focalizzata su

l istante iniziale e l istante finale, nel senso che il fenomeno e

intrinseco al problema e non dipende dalle operazione che si eseguono

per risolverlo- uno stesso problema puo essere mal-condizionato su alcuni dati

e ben-condizionato su altri- ogni algoritmo applicato ad un problema mal-condizionato risulta

instabile

9

Lerrore di troncamento o discretizzazione

Esempio 1: le successioni

an, an= 1 +1

2+ +1

n log(n)

Soluzione matematica: limn an= =0.57721566490152

(costante di Eulero)

Soluzione numerica: aN, N finito

Errore di troncamento: EN = aN.

Per N= 1000: a1000 = 0.57771558156821, E1000 0.5e 3.

10


17/264

Esempio 2: le serie

Problema xy+ y+ xy = 0, y(0) = 1

Soluzione matematica: J0(x) = k=0

(1)k x2k

22k(k!)2, x0

(funzione di Bessel)

Soluzione numerica: SN(x) =N

k=0

ckx2k P2N, N finito

Errore di troncamento o discretizzazione:

|E(x)

|=

|J0(x)

SN(x)

| x2(N+1)

22(N+1)

((N+1)!)2

x= 0.5, N = 5 |J0(0.5) SN(0.5)| 0.11e 12.

11

Errore di troncamento o discretizzazione:

-E presente ogni volta che la soluzione del problema matematico

richiede un numero infinito di operazioni e si costruisce una soluzione

numerica con un numero finito di operazioni.

-E un errore indipendente dallerrore di arrotondamento.

- Nella pratica lerrore di discretizzazione e , in generale, piu rilevante

dellerrore di arrotondamento.

12


18/264

Le formule ricorsive

Una formula ricorsiva e una formula che calcola lelemento ak di una

sequenza in funzione di m

1 elementi precedenti.

m= 1 ak =g(ak1), k1 f. r. ad un passo

m >1 ak =g(ak1, ak2, , akm), km f . r . a m passi

La sequenza e univocamente definita

dal valore dinnesco o valore iniziale a0

dai valore dinnesco a0, a1,

, am

1

Esempio In= 1 nIn1 (v. Lezione precedente - stabilita )

13

Esempio 1. Il fattoriale: fn =n! =n

j=1

j.

Formula ricorsiva: fk =kfk1 k = 1, 2 , n, f 0 = 1

Per n grande fk /F.

32 bits ( 6-7 cifre decimali significative) max fn= 13!

64 bits (15-16 cifre decimali significative) max fn= 21!.

Per valori di n maggiori si perdono cifre significative e solo lordine di

grandezza risuta corretto.

Per esempio con Matlab (64 bits): e(22!) = 32000, e(30!) 0.6e17;e(22!)0.28e 15, e(22!)0.22e 15.

14


19/264

Esempio 2. Algoritmo di corner cutting Dati: N punti P(0)

j R2

Formula ricorsiva:

P

(k+1)

2j =

1

4P

(k)

j +

3

4P

(k)

j+1 P

(k+1)

2j1 = 3

4P

(k)

j +

1

4P

(k)

j+1

spezzata al passo 0

spezzata al passo 1

spezzata al passo 0

spezzata al passo 1

spezzata al passo 2

15

Esempio 3. Algoritmo di di Horner Si vuole calcolare

p(z) =n

k=0

akznk =a0zn + a1zn1 + + an1z+ an.

in un fissato punto z =t.

Algoritmo basato sullespressione di p

z


20/264

Algoritmo di Horner

b0=a0 bk =ak+ t bk1, k = 1, 2, . . . , n

fornisce p(t) =bn.

Costo computazionale: n prodotti.

Algoritmo di Horner generalizzato: si itera la stessa legge ricorsiva:

b(j)0 =b

(j1)0 b

(j)k =b

(j1)k + t b

(j)k1, k= 1, 2, . . . , n j, j = 1, 2, ...n j,

b(0)k =bk, k= 0, 1, . . . , n

b(j)nj = pj

(t)j!

,j = 0, 1, , n 1.

17

Esempio 4.

yn= [(2n)!!]2

[(2n 1)!!]2(2n + 1)/4 (n )

(molto lentamente!) Con il calcolo esplicito dei fattoriali e poi del

quoziente si va in overflaw ( per n= 151 con Matlab).

Formula ricorsiva: posto yn= yn2n + 1

, yn= [(2n)!!]2[(2n 1)!!]2

yn= 2n2n 1

2 yn1, n2, y1 = 4.non ci sono problemi di overflaw.

18


21/264

Le formule ricorsive sono una tecnica di programmazione elegante che

consente di

implementare codici di calcolo molto semplici

ridurre i costi computazionali

proseguire nei calcoli, in casi in cui si avrebbero fenomeni di overflaw

(underflaw)

MA...

bisogna sapere quando fermare il procedimento ricorsivo

se non funziona correttamente si hanno risultati privi di significato

(si veda lesempio sulla stabilita e del fattoriale)

19

I metodi iterativi

Caratteristica comune

cercano la soluzione numerica del problema calcolando successive ap-

prossimazioni a partire da una approssimazione iniziale.

Ad ogni passo del metodo iterativo si vuole ottenere una stima della

soluzione piu accurata di quella calcolata nel passo precedente.

Molti metodi iterativi si basano su formule ricorsive: la conoscenza

delle proprieta e il controllo della ricorsione diventano un problema

cruciale.

20


22/264

il metodo del punto unito in R1

Formula ricorsiva ad 1 passo:

xk+1=(xk), k0, x0 dato: RRLa f. r. rappresenta una trasformazione della retta in se.

Si dice punto unito o punto fisso di ogni R t.c.

= ()

Il metodo iterativo (basato sulla f.r.) viene detto metodo del puntofisso con funzione di iterazione .

21

Non sempre esiste un punto unito (x) =x2 + 1

La curva y =x2 +1non ha punti di intersezione con la bisettricey =x:

la trasformazione x= x2 + 1 non ha punti fissi in R.

y=(x)=x2+1

y=x

x

y

22


23/264

Teorema di esistenza del punto fisso

Teorema 1 Se la successionexk+1=(xk), k0, x0 dato, convergead e la funzione e continua in , allora e punto unito di .

Dim. passando al limite su ambo i membri della f.r.

= limk

xk+1= limk

(xk) =()

per la convergenza per la continuita

In pratica: Se si applica il metodo e si verifica la convergenza i valoriottenuti sono approssimazioni di un punto fisso di .

23

Se ammette un punto unito,

xk =(xk1)

converge?

Non e detto!!

Esempio. 1(x) =ex22 ha un punto fisso in [1, 2]. Infatti

x= 1(x)g(x) =x 1(x) = 0

(un punto fisso di 1(x) e uno zero di g)g(1) = 1 e1 >0, g(2) = 2 e2


24/264

Se si assume x0= 1 si ottengono i valori xk:

k xk k xk1 0.36787944117144 7 0.13793482701375

2 0.15494815269479 8 0.13793482562036

3 0.13862385819704 9 0.13793482556734

4 0.13796111264843 10 0.13793482556532

5 0.13793582594112 11 0.13793482556525

6 0.13793486363172 12 0.13793482556524

quindi la successione converge. A cosa converge?

La convergenza del metodo del punto unito dipende dalla funzione di

iterazione, ma anche dalla scelta dellapprossimazione iniziale.

25

Teorema 2 Sia t. c.

i1 C1[a, b]

i2 (x)[a, b], x[a, b]

i3 |(x)| , 0< 1

26


25/264



(A.A. 2006-2007)


Argomenti

Teorema del punto unito

Criteri darresto

Esempi

Ordine di convergenza

Esercizi

Zeri di equazioni: introduzione

1

Il metodo del punto unito in R1

Il problema: Considerato il metodo iterativo

xk+1=(xk), k 0, x0 dato( metodo del punto fisso con funzione di iterazione )

sotto quali ipotesi e possibile assicurare che

- esiste almeno un punto unito (punto fisso) R t.c.

= ()

- il metodo iterativo converge ad ?

2


26/264

Teorema 2 Sia t. c.

i1 C1[a, b]

i2 (x) [a, b], x [a, b]

i3 |(x)| , 0< 1

3

Dim.

t1 (x) [a, b] (per i2) (a) a, (b) b.Se (a) =a o (b) =b il teor. e dimostarto. Altrimenti

(a) a >0

(b) b


27/264

t3 Segue immediatamente per induzione da i2.

t4 Per il Teorema di Lagrange

xk

= (xk1

)

() =(

k) xk1 (*)

con k nellintervallo di estremi xk e , per t3, k [a, b].Posto ej =xj , la (*) si scrive

ek =(k)ek1

|ek| |ek1| e, per induzione|ek| k|e0|.

Essendo

(0, 1),

lim

k |ek

| |e0

| lim

kk = 0

e per il Teorema del confronto limk

|ek| = 0 .

5

t5 Osserviamo poi che

|xk+1 xk| = |(xk) (xk1)| |xk xk1| e, per induzione

|xk+1 xk| k|x1 x0|

|xk | |xk xk+1| + |xk+1 | |xk+1 xk| + |xk |

(1 )|xk | |xk+1 xk| k|x1 x0| .

6


28/264

Alcune conseguenze pratiche della dimostrazione

|ek| k|e0| k(b a)consente di stimare il numero di iterate sufficienti ad assicurare una

prefissata accuratezza :

|ek| k

(b a) da cui k log

(b

a)

log

(1 )|xk | |xk+1 xk| |xk xk1| da cui

|xk | 1|xk xk1|suggerisce di fermare il calcolo delle iterate controllando|xk xk1|:|xk xk1| , |xk | 1.

ek =(k)ek1 implica che se(x) >0, x [a, b] xk e monotona(x)


29/264

Esempio. 1(x) =

log x + 2, x 1.

1(x) = 1

2x

log x+2>0, x 1 1 crescente (>0) , per x 1.

1(1) 1(x) 1(2), calcolando

1(x) [1.42, 1.65] [1, 2]

ammette (almeno) un punto unito [1, 2].

1(x)= 1

2x

log x+2 1

2

2 0.35 =

quindi e unico e il metodo del punto unito converge ad .

Stima delle iterate: = 0.5e7 (arrotondamento sulla 7a cifra di partedecimale), b a= 1

k log0.5e7log0.15 8.9

9

Esempio. 2(x) = log x + 2

xha gli stessi punti fissi di 1(x) =

log x + 2. Infatti

x= 2(x) =log x + 2

x x2 = log x + 2 x= 1(x) =

log x + 2.

2(x) = log x + 1x2 0, x [1, 2] 2(x) crescente

|2(x)| max|2(1)|, |2(2)|

= 1

NON sono verificate le condizioni del teorema!

x0= 1, condizioni darresto: |xk xk1| 0.5e 7.

10


30/264

k |2(xk)| |1(xk)|1 2.00000000000000 1.414213562373102 1.34657359027997 1.531852992385363 1.70622927805222 1.557715670057004 1.48531376629473 1.563080431084675 1.61287540365716 1.564179820079056 1.53639800411291 1.564404554451687 1.58125740330165 1.564450470677778 1.55459847984262 1.564459850997149 1.57032013836520 1.5644617672809110 1.56100622317187 1.5644621587523911 1.5665092117402112 1.5632526725662813 1.5651780035462414 1.5640390755957915 1.5647125857613716 1.5643142249440317 1.5645498161472618 1.5644104776568719 1.5644928849484420 1.56444414663178

21 1.5644729716301722 1.5644559236992623 1.5644660062849524 1.5644600431665625 1.5644635699126526 1.5644614840995127 1.5644627177051428 1.5644619881175529 1.56446241961521

11

0.8 1 1.2 1.4

y=(log(x)+2)1/2

y=x

x0

x1 x2

0.8 1 2 2.2

y=(log(x)+2)/x y=x

x0

x1

x3

x2

F. di iterazione 1 F. di iterazione 2

Il metodo del punto unito con funzione di iterazione 2 converge

(lecondizioni del teorema sono solo sufficienti!)

Il numero di iterate richiesto dalle due funzioni di iterazione e

molto diverso!

Come misurare della velocita di convergenza delle successioni?

12


31/264

Ordine di convergenza

{xk, k 0} convergente ad ha ordine di convergenza p 1 se esisteuna costante C >0, finita tale che

limk

|xk+1 |

|xk

|p

=C

C costante asintotica dellerrore

Per estensione: il metodo iterativo e di ordine p se la successione che

genera ha ordine di c. p.

La definizione implica: |ek+1|=C|ek|p (uguaglianza asintotica)

p= 1 convergenza lineare - deve essere C


32/264

Esempio Per 1, 2 si ha i(x) = 0, x [1, 2]

il metodo del punto unito con funzione di iterazione i haconvergenza lineare.

Costante asintotica: min(|(x)|) C max(|(x)|)

1: 0.15 C 0.35

2: 0.42 C


33/264

1. (x) =12

x

x [2, 1] 2.5 (x) 1.5 il m. del punto unito non puoconvergere ad 1 perche ad ogni passo si ha

|ek| = |(xk) (1)| =()|ek1| > |ek1|

x [0, 1] 0.5 (x) 0.5 il m. del punto unito convergead 2.

(x) = 0 x = 12, ma (12)= 12 qundi x = 12 non e punto fisso (2)= 0 il metodo ha ordine di convergenza p= 1 e costanteasintotica C 0.5.

(x) non e a segno costante in [0, 1], non si puo dir nulla a priorisullandamento della successione.

Stima delle iterate: = 0.5 (da k(b a) )

k log(0.5e 6)log(0.5)

20.9

17

x0= 0.5

k x0= 0.5 x0= 1 x0 = 1.51 0.37500000000000 -0.75000000000000 -1.625000002 0.36718750000000 -0.40625000000000 -1.88281250

3 0.36618041992188 -0.03564453125000 -2.463897714 0.36604615999386 0.23154246807098 -4.017344805 0.36602818437380 0.33896527677529 -9.828202036 0.36602577630891 0.36203390895797 -52.960878647 0.36602545369318 0.36548267886129 -1428.657772358 0.36595254515683 -1021245.594138469 0.3660156399250310 0.3660240956276511 0.36602522852381

12 0.36602538030395

18


34/264

Esercizio (Proposto)

Si considerino, per x [3, 4] le funzioni1(x) =

4(x1)2 + 3, 2(x) = 1 +

4+3(x1)2

x , 3(x) = 1 +

4

x3

1. Verificare che le funzioni date ammetto tutte lo stesso punto unito

in [3, 4].2. Stabilire per ciascuna delle funzioni se il procedimento iterativo

xk =i(xk1) converge ad .

3. Nei casi in cui il procedimento converge si precisino le proprieta

della successione.

19

Equazioni non lineari

Il problema: determinare le soluzioni o radici, dellequazione non li-

neare

f(x) = 0

cioe quei valori tali che f() = 0 (zeri della funzione ).

Idea base dei metodi numerici: costruire una successione{xk} tale chelim

kxk =

Si assume come approssimazione di il valore xN, per un certo indice

N, finito, opportunamente individuato.

Generalmente per costruire

{xk

} si usano tecniche ricorsive.

Una possibilita :xk =(xk1), k 1, x0 dato

cioe il metodo del punto unito.

20


35/264

I passi per risolvere numericamente il problema

I. Separazione delle radici: per gestire la legge che genera la succes-sione{xn} occorre sapere (in maniera piu o meno accurata) la dislo-cazione degli zeri di f

II. Costruzione della successione{xn}: si tratta di scegliere la leggeper generare la successione in modo da garantire la convergenza a ;

in pratica si vorrebbe anche una convergenza veloce

III. Scelta dei criteri darresto: poiche lalgoritmo iterativo e infinito,

si deve troncarlo, opportunamente.

21


36/264



(A.A. 2006-2007)


Argomenti

Separazione degli zeri

Esempi

Cerchio di inclusione degli zeri di un polinomioMetodo di bisezione

Radici multiple: definizioni

1

Separazione delle radici

Ipotesi di base: lequazione f(x) = 0 ha una (almeno) soluzione

separare una radice= determinare un intervallo [a,b] tale che:

- f C[a, b]

- [a, b]

- e lunica radice dellequazione in [a, b]

2


37/264

Determinare [a, b]

per via grafica

a) si individuano gli intervalli in cui il grafico della curva y = f(x)

attraversa lasse delle x

b) si pone f(x) = 0 nella forma equivalente f1(x) = f2(x), si indi-

viduano i punti di intersezione delle curve y =f1(x) e y =f2(x) e di

conseguenza [a, b].

per tabulazione della funzione: [a, b] e tale che f(a)f(b)0,x possono esistere radici solo per|x|> 1

f1(1) =e > f2(1) = 0,

ma l esponenziale tende allinfinito piu velocemente di ogni potenza

lequazione non ha radici positive

4


38/264


39/264

Se si vuole restringere lintervallo si puo procedere per tabulazione

f(1.587).55e 3 > 0 (lo zero e molto prossimo a 1.587) f(/2).12< 0[/2, 34]

Si osservi che quando si parla di separare gli zeri di f non si richiede

che lintervallo sia piccolo; tuttavia in molte circostanze e opportuno

o necessario affinare lintervallo di separazione.

7

Separazione delle radici: esempio 3 f(x) = cos x xex = 0

g(x) = cos x= xex =h(x)

Analisi preliminare:

| cos x| 1 eventuali radici si hanno per x tale che|xex| 1

Per x > 0 h(x) > 0 e crescente; h(x) > 1, x > 1 possono esistereradici positive solo in [0, 1]; se esiste una radice e unica perche g(x) e

decrescente in [0, 1]

Per tabulazione:

f(0) = 1> 0, f(1) =

2.18< 0

!

[0, 1]

8


40/264

Per x < 0 h(x) < 0 e h(x) 0, x (ex e un infinitesimo diordine superiore alla potenza) esistono infinite radici negative. Leradici devono cadere negli intervalli in cui cos x


41/264

Separazione delle radici: esempio 4

Si studi il numero e la dislocazione delle radici dellequazione

f(x) =ex x= 0

al variare del parametro reale .

Si tratta di individuare le intersezioni fra la curva y =ex e le rette del

fascio per lorigine y =x

Per = 0 lequazione diventa ex = 0 e quindi non ha soluzioni.

11

3 2.5 2 1.5 1 0.5 0 0.5 15

0

5

10

15

= 5

= 1 y=ex

= 0.5

1 0.5 0 0.5 1 1.5 2 2.5 35

0

5

10

15

20

25

=4 =2

=1

=0.5

y=ex

Per 0 esistono 2 radici 0< 1< , < 2


42/264

Chi sono e ?

e il valore di per cui la retta del fascio e tangente a y =ex e e il

punto di tangenza. Sono quindi individuati dal sistema

ex =x uguaglianza delle funzioniex = uguaglianza delle derivate

da cui si ricava = e, = 1

Riassumendo


43/264

Equazioni algebriche: esempio 1

Separare gli zeri reali di p(x) =x3 9.5x2 + 5x + 2.3

p di grado dispari, a coefficienti reali R (almeno)

a = max(2.3, 10.5)b= max(1, 16.8), c= 2 max(9.5,

5, 3

2.3)

= 10.5 zeri reali [10.5, 10.5].

Punti estremali di p: p(x) = 0x1=0.3, x2=6.1

p(x1)> 0, p(x2)< 0 p(x)< 0, x a, p(x)> 0, xa

Zeri di p 1[10.5, 0.3], 21[0.3, 6.1], 3[6.1, 10.5]

15

Equazioni algebriche: esempio 2

Stabilire se lequazione p(x) = 31x4 8x3 17x 20 = 0 ammette zerireali.

p di grado pari zeri reali in numetro pari (0).

a = max

20

31, 1 +

17

31

, b = max

1,

1

31( 8 + 1 7 + 2 0 )

,

c = 2 max 8

31, 317

31, 420

31 = b = 1.45

x

yy=p(x)

16


44/264

Esercizi

I. Separare le radici delle seguenti equazioni:

1. f(x) =x(x 1) x 0.5 = 02. f(x) = 1 x cosh x= 0

II. Data lequazione f(x) =x cos(x + ) = 0, [0, ] individuarei valori di per cui lequazione ammette una radice [0, /2].

III.Stimare il raggio del cerchio di inclusione degli zeri dei polinomi

p(x) = 2x3

+

23

5x2

18

7p(x) = 256x4 26x3 + 226x 120

17

Costruzione della successione: il metodo di bisezione

Il metodo di bisezione o dicotomico si basa sul Teorema degli zeri.

Ipotesi di applicabilita :i1 E dato un intervallo, I= [a, b], di separazione della la radice

di f(x) = 0

i2 f(a)f(b)< 0.

i1 f C[a, b] e e lunica radice in [a, b]

i2

ha molteplicita dispari.

18


45/264

Step 1. Calcolare, x= (a + b)/2, (punto medio dellintervallo)

Step 2. Se f(x)= 0 ( f(x) = 0x= soluzione ), allora

[a1, b1]=

[x, b1] se f(a)f(x0)> 0

[a, x] se f(a)f(x0)< 0

I passi precedenti si applicano quindi di nuovo allintervallo [a1, b1] e

poi iterativamente.

19

a0

b1=b

0

f(x)

x

x1

a2=a

1

x2

b2

b1 a1=b a

2 , b2 a2=

b1 a12

=b a

22 ,

20


46/264

Il metodo costruisce quindi le successioni

a0 =a, b0=b {ak} ,{bk} , k >0Proprieta delle successioni

p1 akak1, bkbk1,kp2 akbk kp3 bk ak =

b a2k

{ak} ,{bk} sono monotone, limitate convergono a .

I valori ak, bk forniscono approssimazioni rispettivamente per difetto eper eccesso della radice .

21

Convergenza del metodo di bisezione

Errore di troncamento: ek =xk

|ek

|= xk

< bk

ak =

b a

2k

|ek+1|= 12|ek| (asintoticamente)

Il metodo ha ordine p= 1 e costante asintotica C= 12

La convergenza e lenta, in quanto ad ogni passo lerrore viene dimez-

zato, = ad ogni passo si guadagna una cifra binaria

per guadagnare una cifra decimaleservono 3-4 iterazioni.

22


47/264

Stima del numero di iterate

|ek| b

a

2k

|ek| se (C.S.)b a

2k k log(ba)log()log2

Condizione darresto bk

ak < , tolleranza prefissata .

23

Esempio f(x) =x3 x 1, [a, b] = [1, 2]

Stima delle iterate: = 108:

2k

108

k

26.6

a26= 1.3247179 5380116<


48/264

Esercizio

Data lequazione f(x) =x3/2 2 log(x2 + 0.6) = 0I. Separare le radici con un intervallo di ampiezza non superiore a 0.5

II.Stimare quante iterate sono sufficienti per ottenere unapprossimazione

alla 10 cifra decimale della radice maggiore con il metodo di bisezioni.

I. f(x) = 0f1(x) =x3/2 =2 log(x2 + 0.6) =f2(x).f1, f2 entrambe crescenti, f1>0, x >0, f2(0)< 0Landamento tipo puo essere uno dei seguenti

25

0 0.5 1 1.5 2 2.5 32

0

2

4

6

8

10

12

14

16

Caso a:la potenza si mantiene sempre aldi sopra del logaritmo

0 0.5 1 1.5 2 2.5 32

1

0

1

2

3

4

5

6

Caso b: se la potenza i ncontra il logaritmo in un punto,allora gli due zeri sono necessariamente due

f1 e sempre sopra f2: non ci sonoradici

f1 incontra f2 in un punto: allora cisono necessariamente 2 radici

Pertabulazione si individuano le variazioni di segno:

f(0) =2log0.5> 0, f(1.5) 0.2< 0, f(3)0.7> 0

quindi esiste una radice 1[0, 1.5] e una radice 2[1.5, 3]

26


49/264

Affinamento dellintervallo(ampiezza 0.5)

- Radice 1: f(0)> 0, f(1)> 0 1[1, 1.5]

- Radice 2: applichiamo bisezioni (f(1.5)< 0, f(3)>0)

n ak() bk(+) x f(x)

0 1.5 3 2.25 -.09

1 2.25 3 2.625 0.23

Segue 2[2.25, 2.625] = [a2, b2] e b2 a2= 0.375< 0.5.

II. Stima del numero di iterateba2k

121010 klog2log

2 0.375 1010

n32.8

27

Alcuni richiami

si diceradicedi f(x) = 0 ( zero di f) ,di molteplicita 1 se esiste

limx

f(x)

(x ) == 0, finito ()

Se f Cr(), 1 r 1 allora, e zero di molteplicita rper f(r) e

limx

f(r)(x)

(x )r = !

( r)!

Se f C() la condizione (*) equivale a

f(j)() = 0, j = 0, 1, . . . ,

1, f()()

= 0

e =f()()

! .

28


50/264

CALCOLO NUMERICO

Laurea Specialistica Ingegneria Elettronica, Gestionale, Informatica


(A.A. 2006-2007)


Agomanti

Metodo di Newton

Interpretazione geometrica

Teoremi di convergenza

Esempi

1

Costruzione della successione: metodo di Newton

Si basa sullosservazione che, in un intorno I di , la funzione f puo

essere approssimata con la sua tangente in un punto appartenente ad

I.

Equazione della tangente t a f nel punto (x, f(x)):

t: y f(x) =f(x)(x x)

Lintersezione di f con lasse delle x ( y = 0) viene approssimata con

lintersezione di t con lasse delle x:

x =x f(x)

f(x)

2


51/264

Interpretazione geometrica del metodo di Newton

Prima iterazione:

Approssimazione iniziale: x0Intersezione di t0 (tangente a f in (x0, f(x0))) con y = 0

x1 =x0 f(x0)f(x0) nuova approssimazione di

x

a b

f(x)

x0

t0

(x0

,f(x0

)) .

x1

3

Seconda iterazione:

Approssimazione iniziale: x1Intersezione di t1 (tangente a f in (x1, f(x1))) con y = 0

x2 =x1 f(x1)f(x1)

nuova approssimazione di

x

a b

f(x)

x0

t0

(x0,f(x0)) .

x1

t1

.(x

1,f(x

1))

4


52/264

La funzione diterazione del metodo di Newton

In generale, se f almeno in un intorno di

x0 assegnato

xk =xk1 f(xk1)

f(xk1), k = 1, 2, . . .

Metodo del punto unito con funzione diterazione: N R(x) =x f(x)f(x)

f() = 0 =N R() qualunque sia la molteplicita . Infatti:

se f()= 0 e immediato.

se ha molteplicita

N R(x) =x (x ) f(x)

(x )(x )1

f(x)

(da limxf(r)(x)

(x)r =r, r0 (v. lezione precedente))

limx

N R(x) = N R e prolungabile con continuita in con

=N R().

5

Il metodo di Newton: un esempio

f(x) = 4x3 5x Vogliamo approssimare = 0

a) x0= 0.5

0.5 0.4 0.3 0.2 0.1 0 0.1 0.2 0.3 0.4 0.52

1.5

1

0.5

0

0.5

1

1.5

2

f(x)

.

x

t0

a b0.5 0.4 0.3 0.2 0.1 0 0.1 0.2 0.3 0.4 0.52

1.5

1

0.5

0

0.5

1

1.5

2

f(x)

.

x

t0

t1

a b

x2k = 0.5 x2k+1=0.5

6


53/264

b) x0= 0.35

k xk |xk xk1| f(xk)

1 -0.09716713881020 0.447e+00 0.482166

2 0.00150186993989 0.987e-01 0.750934e-02

3 -0.00000000542025 0.150e-02 0.2710125e-07

4 0 0.542025e-08 0

7

c) x0 = 0.65

k xk |xk xk1| f(xk)

1 31.385714286 3.138e+01 1.2351071e+06

2 20.9326637216 1.046e+01 3.6584134e+04

3 13.9683918420 6.964e+00 1.0831983e+04

4 9.33218995276 4.636e+00 3.2042921e+03

5 6.251368612684 3.081e+00 9.4594733e+02

6 4.212492657810 2.039e+00 2.7794186e+02

7 2.87585544490976 1.337e+00 8.0760285+01

8 2.01895088696529 0.857e+00 2.8235346e+00

9 1.49921778432550 0.520e+00 5.9828023e+00

14 1.11803398874990 0.647e-07 6.0e-14

15 1.11803398874989 0.1e-13 0

8


54/264

Cosa succede dal punto di vista geometrico

0.5 0 0.5 1 1.5 210

5

0

5

10

15

20

25

t0

x0

=0.35

Lapprosssimazione iniziale e vicina ad un punto estremale: si esce

dallintorno della radice cercata.

La successione e convergente e quindi (Teorema 1) converge comunque

ad uno zero di f

9

Metodo di Newton: convergenza

xk =N R(xk1) =xk1 f(xk1)

f(xk1), k = 1, 2, . . . ; x0dato

Teorema di convergenza 1. Sia lunica radice di f(x) = 0 in [a, b],

di molteplicita 1

i1 fC2[a, b]

i2 f(x)= 0,x [a, b], x=.

t1 un intorno di , I() [a, b] tale che x0 I() il metodo di

Newton converge.

t2 Se = 1 allora p 2

t3

Se >1 allora p= 1, C= 1 1

10


55/264


56/264

Sia radice di molteplicita : per x= si puo scrivere

N R(x) = f(x)(x )

f(x)(x )2

(x )1

f(x)

2

e dalla limxf(r)(x)

(x)r = !

(r)!,= 0 passando al limite, si ottiene

limx

N R(x) = !

( 2)!

( 1)!

!

2=1

1

= 0

N RC0() p = 1, C= 1 1

(teorema sullordine di convergenza).

13

Definizione. Sia f C2[a, b] tale che f(a)f(b) < 0, f(x) = 0, x

[a, b]. Si dice estremo di Fourier di f quello dei due estremi a, b per il

quale risulta f()f()> 0.

a

b

a

b

a

b

a

b

I quattro casi possibili di estremo di Fourier

14


57/264

Teorema di convergenza 2. Sia data lequazione f(x) = 0 avente una

radice [a, b] e siano verificate le seguenti ipotesi

i1 fC

2

[a, b]i2 f(a)f(b)< 0

i3 f(x)= 0,x [a, b].

allora

t1 e lunico zero di f in [a, b]

t2 sex0 coincide con con lestremo di Fourier dif, {xk} generata dalmetodo di Newton e monotona e convergente a con ordine p= 2.

15

Esempio f(x) = 4x3 5x= 0 (v. prima) Si ha

f(x) = 24x >0 x >0 e f(x)> 0 x >2

ogni valore di x >2 e estremo di Fourier per la radice maggiore

Nel caso x0 = 0.65 si ha x1 31.4 = lestremo di Fourier Il m. diN. converge alla radice maggiore (Teor. di convergenza 2)

Se si sceglie per esempio x0 [0.4, 0.4] il m. di N. converge a = 0

(Teor. di convergenza 1)

NOTA: A priori e generalmenete impossibile sapere quanto x0 deve

esserevicino alla radice.

16


58/264

Esempio Si consideri lequazione f(x) = cos x xex = 0 (v. Lezione

precedente).

Si vuole approssimare [,/2].

f(x) = cos x ex(x + 2) non e a segno costante in [, /2].

Si applica il teor. di convergenza 1.

x0 = 3/4

3a iterate: 6 cifre decimali

4a 14 cifre decimali (si conquistano tutte le cifre con cui si lavora):

=1.86399519240253

Verificate che per approssimare la radice positiva si puo appliacare il

teor. di convergenza 2.

17

Esempio p(x) = x3 9.5x2 + 20x+ 12.5 = 0: si vuole verificare il

comportamento del metodo di Newton quando si approssima la radice

doppia.

Il polinomio (v. L02) ha uno zero negativo e uno zero doppio = 5.

Imponendo unaccuratezza pari a 1010 si ottiene

n. iterate (N) x0 xN p(xN) p(xN)

27 8.9 5.00000006170276 2.842e-014 6.787e-00726 4.5 5.00000004576961 0 5.035e-007

Ordine di convergenza: p= 1, C= 12 (Teorema di convergenza 1)

E evidente limpossibilita di migliorare laccuratezza raggiunta (mal-

condizionamento)

18


59/264



(A.A. 2006-2007)


Argomenti

Condizionamento del calcolo degli zeri di polinomi

Efficienza computazionale

Il metodo delel secanti

Esempi ed esercizi

1

Condizionamento del calcolo degli zeri di un polinomio

p(x) =n

j=0

ajxnj = 0, t.c. p() = 0

Fissato i aiai+ i + ()

Condizionamneto: determinare il fattore che lega

er(ai) =i

aie er() =

()

Dalla formula di Taylor:

2


60/264

= 0p (a0, , ai+ i, , an, + ())=

p(a0, , an, ) + p

ai

a,

i + p

a,

() + R2

= 0

=p()

p

ai=

ai

nj=0

ajxnj =xni

0nii+p()() |()|

ni

p()

|i|

Da questa relazione sugli errori assoluti si ricava quella sugli errori

relativi

()

aini1

p()

iai

3

Alcune osservazioni

La determinazione di una radice multipla di un polinomio e un

problema malcondizionato

In generale se si applica Newton per approssimare una radicemultipla di f, possono sorgere problemi di propagazione degli errori di

arrotondamento nel calcolo di f(xk)

f(xk)

In presenza di radici di molteplicita si puo ristabilire la conver-

genza quadratica del metodo di Newton assumendo come funzione di

iterazione

(x) =x f(x)

f(x). Nonsi evitano pero i problemi di calcolo

detti prima.

4


61/264

Esercizio 1 Lequazione f(x) = log xx

+x2 5.4 = 0 ha una radice[2, 3].

I. E possibile approssimare con il metodo di Newton?

II. In caso affermativo si precisi la scelta dellapprossimazione iniziale

e lordine di convergenza del metodo.

Soluzione Conviene vedere se e applicabile il Teorema 2 (funzioni a

concavita fissa. Si ha f C[2, 3]. Occorre verificare se

- f e a segno costante- la radice e semplice (f a segno costante)

5

Posto f(x) =g(x) + x2 5.4, g(x) =x1/2 log x si ha

g(x) =x3/2 12x3/2 log xf

(x) =x

3/2

1 1

2log x + 2x

g(x) =32x5/2

1 12 log x

12x5/2 f(x) =14x5/2 (8 3log x) + 2

ed infine g(x) =f(x)= 58x7/2 (8 3log x) + 34x7/2 >0 segue che

f e monotona (quindi assume i valori estremali agli estremi)f(2)1.7, f(3)1.9f(x)> 0, x[2, 3]

6


62/264

f >0 implica f crescente:f(2)4.2, f(3)6.1f(x)> 0, x[2, 3] la radice e semplice.

Approssimazione iniziale nellestremo di Fourier:f(2) 0.91, f(3)4.2x0= 3Il metodo di Newton converge con ordine p= 2.

n xn n xn0 3.00000000000000 3 2.20619119992521

1 2.30434241770304 4 2.20619053704335

2 2.20802828937289 5 2.20619053704326

7

Esercizio 2 Si consideri lequazione x log(2x + 1) = 1

I. Quante radici ammette?

II. A quale radice converge il metodo di Newton se si assume come

approssimazione iniziale x0 = 0.5?

Soluzione f e definita in D = (12, +)lim

x1/2+f(x) = +, f(0) = 1, lim

x f(x) = + e anche ungrafico sommario delle curve

f1(x) = log(2x + 1) , f2(x) =1

xper x >12

conferma che lequazione ammette due radici una positiva e una ne-

gativa.

8


63/264

x0= 0.5 e estremo di Fourier?

- f(0.5) 0.65< 0 (la radice positiva verifica >0.5)

- f C(D)- f(x) = log(2x + 1 ) + 2x2x+1, sign f

(x) = sign x, xD

- f(x) = 4 x+1(2x+1)2

>0,x 12 (funzione a concavita fissax)

x0= 0.5 non e un estremo di Fourier.

Calcoliamo x1 = 0.5

f(0.5)

f1(0.5)

= 1.047649460

Ma f(x1).18> 0 x1 e estremo di Fourier per la radice positiva.

9

Efficienza computazionale

Lefficienza di un metodo iterativo tiene conto sia dellordine di con-

vergenza p che del costo computazionale, cioe della quantita di calcoli

richiesta ad ogni passo, valutato come numero di valutazioni funzionali

richieste ad ogni passo, r

Efficienza: E=p1/r

Esempio Metodo di bisezione:

p= 1, r = 1 (valutazione nel punto medio), E= 1

In generale, per il metodo del punto unito aconvergenza lineare: E= 1

10


64/264

Schema riassuntivo sul metodo di Newton

f(x) = 0 [a, b] radice f C2[a, b]

1. f()= 0 (radice semplice): x0I(), p2, E

2

2. f()= 0, f= 0, x[a, b] (radice semplice, concavita fissa):

x0 = estremo di Fourier, p= 2, C=

f()

2f()

, E=

2

3. f(k)() = 0, k = 0, , 1, f()()= 0 (radice di molteplicita ):

x0I(), p= 1, C= 1 1

, E= 1,

11

Il metodo delle secanti

Il metodo di Newton richiede 2 valutazioni funzionali.

E possibile costruire un metodo veloceche richieda solo una valu-

tazione funzionale?

Si sostituisce alla derivata prima un rapporto incrementale occorronodue punti

xk+1=xkf(xk)(xk xk1)

f(xk) f(xk1), k1, x0, x1 dati

E una f. ricorsiva a 2 passi non e un metodo del punto unito deltipo visto.

E un metodo di linearizzazione.

12


65/264

Interpretazione geometrica del m. delle secanti

Step 1 s0 retta per i punti (x0, f(x0)) , (x1, f(x1)),

x2=s0 {y = 0}

Step 2 s1 retta per i punti (x1, f(x1)) , (x2, f(x2)),

x3=s1 {y = 0

}

1.5 2 2.5 3 3.53

2

1

0

1

2

3

4

5

x0

x1

x2

s0

s1

13

Teorema di convergenza Data f(x) = 0, sia [a, b] una sua radice.Se

i1 f C2[a, b]

i2 f(a)f(b)< 0

i3 f()= 0

i4 f(x)= 0,x[a, b]

allora

t1 la successione generata dal metodo delle secanti e convergente a

comunque si scelgano x0, x1

[a, b] tali che f(

)f(

)> 0

14


66/264

Ordine di convergenza Se converge a , radice semplice di f(x) = 0,

allora il m. delle secanti ha ordine di convergenza

p=1 +

5

2 1.62 con C=

f()2f()

1/p

Dim. Posto ej =xj dalla formula iterativa sottraendo ad ambo imembri

ek+1=ek1f(xk) ekf(xk1)

f(xk) f(xk1)Applicando Taylor

ek1f(xk) ekf(xk1)= ek1

ekf() + f

()2 e

2k+ R1

ek ek1f() + f()2 e2k1+ R2= ekek1(ek ek1)f()2 + Rf(xk)f(xk1)= f()ek +R2

f()ek1+ R3

= f()(ek ek1)+R

15

Sostituendo e trascurando i resti

ek+1=ekek1

f()2f()

(*)

Dalla definizione di ordine di convergenza:

ek=C epk1ek1=e

1/pk C

1/p sostituendo in (*):

ek+1=e1+1/p

k C1/pf

()2f()

Ma deve anche essere ek+1= Cepk e, uguaglaindo membro a membro

p=1

p+ 1, C=C1/pf

()2f()

C1p +1 =Cp =

f()2f()

p =

1 +

5

2

.

Ad ogni passo (escluso il primo) si deve eseguire solo una nuova valu-

tazione di funzione (r = 1) E=p1.6216


67/264

Esempio f(x) = ex = x, per = 3 (vedi L04). Nella tabella sono

riportate le approssimazioni fornite dal m. di Newton (x0 = 2) e dal

m. delle secanti (x0= 2, x1= 1.9).

k Newton Secanti

1 1.68351826278341 1.655459273962712 1.54348197207520 1.563650594911693 1.51348862277163 1.521063585702424 1.51213725012625 1.512769802200435 1.51213455166859 1.512142834968896 1.51213455165784 1.512134559420427 1.512134551657948 1.51213455165784

Costo computazionale totale

Newton: rtot= 12

Secanti: rtot= 9

17

Esercizi

I. Stabilire se e possibile applicare il metodo delle secanti in condizioni

di convergenza per approssimare lo zero di f(x) = 1 x cosh(x), pre-cisando, in caso affermativo, i valori di innesco.

II. Si consideri lequazione f(x) = (x + 2)log(x) + 1 = 0.

II.1 Si localizzi la radice dellequazione.

II.2 Si studi se e possibile approssimare con il metodo di Newton,

precisando, una volta fissata lapprossimazione iniziale, le proprieta

della successione.

II.3 Verificare se il metodo del punto fisso con funzione di iterazione

(x) = e1/(x+2)

e convergente a e, in caso affermativo, stimare ilnumero di iterate sufficienti ad approssimare la radice alla 10a cifra

decimale.

18


68/264

Note allestensione nel campo complesso

La ricerca di radici C e un problema prevalentemente associato alcaso delle equazioni algebriche.

- La separazione degli zeri per tabulazione o per via grafica non e

applicabile in C

- Il metodo di bisezioni non e applicabile in C

- Il m. di Newton/delle secantipuo essere applicato in C emantiene le

proprieta viste (in particolare la convergenza quadratica/superlineare,

se la radice e semplice).

- Nel caso di equazione algebrica, se ak R, k e si assume x0 R /x0, x1 R il m. di Newton/delle secanti non puo convergere a unaradice complessa ma solo ad una radice reale (se esiste).

19

Criteri di arresto

Oltre ai criteri (v. L03) dellerrore assoluto (errore relativo) e del

numero massimo di iterate, si puo pensare di utilizzare

|f(xn)|

basato sul concetto di continuita .

Ma f(xn) non sempre e significativo

a b

f(x)

x

xk

a

f(x)

x b

xk

f(xn) e grande anche se xn

e vicino a

f(xn) e piccolo anche se xn

e lontano da

20


69/264



(A.A. 2006-2007)


Argomenti

Norma di vettoreDistanza, Intorno, convergenzaNorme holderiane

Norma di matriceNorma sub-moltiplicative, compatibile, indotta

Il metodo del punto unito in Rn

Convergenza

1

Esempio Sistema lineare dato

A=

8 1 11 5 1

1 1 5

x1x2

x3

=

267

7

.

Soluzione: xv = (3, 1, 1)T.

Approssimazioni: v1= (2.99976750, 1.00023650, 1.00030050)T

v2 = (3.00027250, 0.99961650, 0.99955250)T

Come possiamo stabilire quale delle due approssimazioni e piu vicina

alla soluzione del sistema?

Vettore dellerrore:

e1

= xv

v1

= (0.2325e

3,

0.2365e

3,

0.3005e

3)T

e2 = xv v2 = (0.2725e 3, 0.3835e 3, 0.4475e 3)T

quale differisce meno dal vettore nullo?

2


70/264


71/264

Le norme holderiane

Norma holderiane di parametro p o norma p, p 1:

x = n

j=1 |xj|p

1/p

I casi piu usati:

p= 1 x =n

j=1

|xj| (norma 1)

p= 2

x

=

n

j=1 |xj|2 (norma 2 o euclidea)

p= x = max1jn

|xj|

(norma infinito o uniforme)

5

Lintorno unitario in norma 1, 2, infinito

n = 2, e quindi un vettore x = (x1, x2)T R2 e rappresentabile come

un punto del piano cartesiano.

Intorno unitario: I1(0, p) = x R2 t.c. xp 1Frontiera dellintorno:

x R2 t.c. xp = 1

p= 2 (e il caso della distanza euclidea)

x2=

|x1|2 + |x2|2

= 1 |x1|2 + |x2|2 = 1

ovvero I1(0,

2) = cerchio di centro lorigine e raggio 1

6


72/264

p=

x= max (|x1|, |x2|)= 1,|xi| 1 x1= 1, x2= 1

x

y

(1,0) (1,0)

(0,1)

(0,1)

I1(0, ) = quadrato con lati agli assi

7

p= 1

x1= |x1| + |x2| = 1,|xi| 1

x1+ x2= 1, x1+ x2= 1, x1 x2= 1, x1 x2 = 1

x

y

(1,0) (1,0)

(0,1)

(0,1)

I1(0, 1) = quadrato con lati alle bisettrici8


73/264

norma uniforme

norma 1

norma 2

O x

y

(0,1)

(0,1)

(1,0)(1,0)

I(x; 1) I(x; 2) I(x; )

una condizione del tipo: x(k) xv

e piu restrittiva se si usa la norma 1, e piu debole se si usa la norma

uniforme.

9

Esempio Riprendiamo lesempio dato allinizio: si vuole una approssi-

mazione di x tale chexv x 0.5e 3. Calcoliamo :

e11 = 0.2325e 3 + 0.2365e 3 + 0.3005e 3 0.77e 3,

e12 = (0.2325e 3)2 + (0.2365e 3)2 + (0.3005e 3)2 0.45e 3,

e1 = max(0.2325e 3, 0.2365e 3, 0.3005e 3) 0.30e 3

e21 0.11e 2; e22 0.65e 3; e2 0.45e 3

1 2

v1 no si si

v1 no no si

10


74/264

Norma di matrice

Si definisce norma di una matrice A Rmn una funzione

: Rmn

R+ {

0

}1.A = 0 A= 0

2.A = || A, R

3.A + B A + B (disuguaglianza triangolare)

NOTA: La definizione si estende al campo complesso

Si estendono le definizioni di distanza, intorno e convergenza in modo

formalmente identico al caso dei vettori.

11

Un esempio: la norma di Frobenius A= ((aij)) Rmn

||A||F := m

i=1

nj=1

|aij|21/2

Attenzione: La formula individua una classe di norme.

Esempio: ||A||F, A R23 =||A||F, A R44

sono diversi gli spazi su cui sono definite le funzioni|| ||F.

Quale relazione esiste fra le norme se si eseguono operazioni (di prodotto)

fra matrici di spazi diversi?

12


75/264

Norme sub-moltiplicative

Una norma di matrice si dice sub-moltiplicativa se

||A B| | | |A| | | |B||

quando il prodotto e definito.

| | | |F e sub-moltiplicativa

Non tutte le norme sono sub-moltiplicative

EsempioA = max1i,jn |aij|, A R

nn

(e immediato verificare che la funzione definita e una norma)

Cosideriamo per n= 2

A= B =

1 10 1

A= B

1 20 1

A B = 2> A B=1

13

Norme compatibili

Una norma di matrice si dice compatibile rispetto ad una norma di

vettore se

norma di matrice

||Ax||

||A| | | |x||

, A Rmn, x Rn

norma di vettore

Come ottenere norme che siano contemporaneamente sub-moltiplicative

e compatibili?

14


76/264

Norme indotte

Si definisce la norma indotta (da una fissata norma di vettore) di

matrice

A = maxx=1

Ax (norma indotta)

Si dimostra che questa definizione ha senso ( il max).

Proprieta delle norme indotte

Ogni norma indotta e una norma sub-moltiplicativa.

Ogni norma indotta e compatibile (rispetto alla norma di vettore

da cui e definita)

In ogni norma idottaI = 1

Infatti:I = maxx=1

Ix = maxx=1

x = 1

15

Norme di matrice dalle norme holderiane

Norma uno: A1:= max1jn

n

i=1

|aij| (per colonne)

Norma infinito: A:= max1in

n

j=1

|aij| (per righe)

Norma 2 o di Hilbert: A2 =

(AHA)

(M) eil raggio spettrale= autovalore di modulo massimodella matrice

M.

||A||F non e unanorma indotta(IF =

n) ma ecompatibilerispetto

alla norma euclidea di vettore:

Ax

2 A

F x2

Chiameremo norma canonica ogni norma di matrice che sia sub-

moltiplicativa e compatibile.

16


77/264

Esercizio Sia data la matrice, dipendente dal parametro reale

A=

3 0

3 2

0.5

0 0.5 1

Studiare il valore diA al variare di in ognuna delle norme canoniche1, e di Frobenius.

a)A2F =2 + 1 8 + 4 + 2 ( 0.5)2 + 1 =32 2 + 23.5

parabola con minimo in = 1/3.

AF = 32 2 + 23.5 decresce da a 4.813 e poi cresce, pervalori di crescenti, e tende a

17

A= AT A1= A= A

A = max (|| + 3, 5 + | 0.5|, 1 + | 0.5|)

A

= max(

|

|+ 3, 5 +

|

0.5

|)

0A = max( + 3, 5 + 0.5) = 5.5 + ||

0 0.5A = max( + 3, 5 + 0.5)

3< + 3 3.5, 5< 5.5 5.5 A = 5.5

>0

A

= max( + 3, 5 +

0.5) = 4.5 +

18


78/264

Alcune proprieta delle norme

Ogni norma e una funzione continua dei suoi argomenti.

Se a, b sono norme introdotte in uno spazio di dimen-

sione finita allora lenormesonoequivalentinel senso che esistono due

costanti c1, c2>0 tali che

c1 a b c2 a

Assicura che la convergenza e indipendente dalla norma.

Casi particolari A Rmn

A2 AF

nA21

nA A2

mA

1m

A1 A2

nA1

19

Trasormazioni e punto fisso in Rn

Trasformazione in Rn:

x1 = 1(x1, x2, . . . , xn)

x

i =

i(x

1, x

2, . . . , x

n)

xn = n(x1, x2, . . . , xn)

Notazione maticiale: x= (xi)T Rn e = (1, , n)

x= (x)

Rn si chiama punto unito della trasformazione se

= ()

20


79/264

Esempio n= 2x1 =1(x1, x2) =

2x1 x22

x2 =2(x1, x2) = 5 log x1

x

y

0

y=5log x

x2+y

22x=0

21

Esistono 2 punti uniti:

1 D1= {[1, 2] [0, 1]}

2 D2= {[0, 1] [1, 0]}

21

Il m. del punto fisso in Rn

x(k) = (x(k1)), k 1, x(0) dato ()

x(j) Rn , = (1, , n)

Teoremi di convergenza

Teorema 1 Se

i1 C0(D)i2 lim

kx(k) =

allora

t1 = ()

Come in R questo teorema ha una validita a posteriori: se si ottiene

una successione convergente, allora il punto limite e un punto fisso.

22


80/264

Teorema 2 Sia D un sottoinsieme chiuso e limitato diRn oRn stesso.

Se e tale che

i1 C0(D),

i2

(D)

D,

i3 esiste una costante 0<


81/264



(A.A. 2006-2007)


Argomenti

Metodo del punto unito -Esercizio

Soluzione di sistemi lineari

Classificazione dei metodi

Metodi iterativi - Condizioni di convergenza

Metodo di Jacobi

1

Esercizio Consideriamo la trsformazione = (1, 2)

x1 = 1(x1, x2) = 0.5

2(3x2+ 15)

x2 = 2(x1, x2) = x1+ log x1

per x D = [3, 4] [2, 3].1. Verificare che ammette un unico punto unito in D.

2. Stabilire se il metodo del punto unito x(k) = (x(k1)) converge ad

3. In caso affermativo stimare il numero di iterate che consente di

ottenere una approssimazione con 7 cifre decimali.

2


82/264

1. 1, 2 C(D), monotone

assumo i valori estremali agli estremi degli intervalli [2, 3] e [3, 4] rispet-

tivamente:

3.25 x1 = 1(x1, x2) 3.47, 2.03 x2=2(x1, x2) 2.33

(1(D), 2(D)) D ! D punto unito di .

2. Proviamo ad utilizzare la C.S. di contrazione. Occorre

a) Calcolare le derivate parziali di 1, 2

b) Calcolare le maggioranti del modulo di ogni derivata: i

xj mijc) Posto M = ((mij)), verificare cheM


83/264

3. Fissiamo x(0) = (3.5, 2.5).

Si ricava x

(1)

= (3.354101966, 2.180083248) x(1) x(0) = (.145898034, .319916752).

Dalla condizione

0.33k 0.5e 7 0.67(.145898034, .319916752)si ricava

k 14.5 se si adotta la norma, k 14.85 se si adotta la norma 1.

5

Soluzione numerica dei sistemi lineari

Dati A= ((aij)) Rnn, b= (bj)T Rn

determinare la soluzione xv Rn din

j=1

aijxj =bi, i= 1, 2, , n (Sistema lineare di ordine n)

In forma matriciale Ax= b

Il Teorema di Cramer assicura che il sistema ammette una e una sola

soluzione il determinante di A e non nullo xv =A1b.

Supporremo pertanto che sia detA= 0

6


84/264

I metodi numerici per la soluzione di sistemi lineari si dividono in due

classi:

Metodi iterativi Metodi diretti

Metodi iterativi

Costruiscono una successione di vettori

x(k)

, k 0 t.c.

limk

x(k) =xv

Si assume come soluzione numerica (approssimazione di xv) un ele-

mento della successione xM

, M finito.

La soluzione numerica e affetta da un errore di troncamento

7

Metodi diretti

Si basano sul

Teorema A Se sul sistema Ax = b si esegue una delle seguenti oper-

azioni:

scambio di due righe moltiplicazione (divisione) di una riga per uno scalare non nullo sostituzione di una riga con una combinazione della riga stessa conaltre righe del sistema

allora

il nuovo sistema ammette soluzione se e solo se Ax= b ammettesoluzione

le soluzioni dei due sistemi coincidono (i due sistemi sono equi-

valenti)

8


85/264

Idea base dei metodi diretti:

trasformare Ax=b in un sistema equivalente di cui sia direttamente

calcolabile la soluzione, applicando le operazioni del Teorema A unnumero finito di volte.

La matrice dei coefficienti del nuovo sistema avra una struttura parti-

colare.

In algebra infinita un metodo diretto fornisce la soluzione esatta

Nelsistema di numeri macchinaun metodo diretto fornisce unasoluzione

(numerica) affetta da errori di arrotondamento

9

Metodi iterativi per sistemi lineari

Passo 1. Ax= b x= Cx + q t.c.

Axv b xv Cxv + q

Tecnica di base

A= M+ N (splitting della matrice) t.c. M1

(M+ N)x= b Mx= Nx + b x= M1N =C

x + M1b =q

Passo 2. Generare la successione con il metodo del punto unito

x(k) =Cx(k1) + q, k= 1, 2, . . .

x(0) Rn datoC : matrice di iterazione

10


86/264

Convergenza dei metodi iterativi per sistemi lineari

Chi e lo jacobiano del sistema?

i(x) =n

j=1

cijxj+ qi ixj

=cij

indipendente da x

La C.S. per la convergenza del m. del punto unito diventa

Teorema 1 Se in una qualche norma canonica risultaC


87/264

Teorema 2C.N e S. affinche il procedimento iterativox(k) =Cx(k1) +q, k = 1, 2, . . . converga e che sia (C)< 1 e la convergenza e indipen-

dente dalla scelta dellapprossimazione iniziale

((C) = raggio spettrale= autovalore di modulo massimo di C)

Relazione fra i due teoremi

B Rnn : Bx= x, x =0 (in una norma compatibile )

x = || x = Bx B x || B,

In particolare (B) B,B Rnn quindi C


88/264

Esempio C=

0.25 0.4 00.34 0.35 0.01

0 0.4 0.3

, q=

10

2

Si vuole stimare il numero di iterate sufficienti ad ottenere unappros-

simazione con una accuratezza = 108.

Condizione da imporre Ck1 Cq 10

8

Innorma e norma di Frobenius e verificata la C.S. di convergenza;si puo usare la norma.

C

= max(0.65, 0.7, 0.7) = 0.7,

q

= 2

0.7k

0.3 2 108 k log0.7 log(0.15 108) k 57

15

Costruzione di C: il metodo di Jacobi

Dalla i-sima riga: ai1x1+ + aiixi+ ainxn = bi

Se aii = 0, si ricava xi:

xi = ai1aii

x1 ai,i1aii

xi1 ai,i+1aii

xi+1 ainaii

xn+ biaii

Ipotesi: aii = 0, i= 1, 2, , n

Ax= b x= Cx + q, C= ((cij)), q= (qi) t.c.

cij =

0 i= j

aijaii i =j , qi=

biaii

16


89/264

La matrice di iterazione di Jacobi dallo splitting di A

A= L + D+ U

A=

a11 a12 a1na21 a22 a2na31 a32 a33 a3n

an1 ann

C=

D1(L + U), q= D1b

17

Condizioni sufficienti di convergenza del metodo di Jacobi

A diagonalmente dominante per righe :

|aii| >k=i

|aik| Ck=i

|aki| C1


90/264

EsempioSi consideri il sistema

[A|b] = 1 2 2 01 1 1 1

2 2 1 1

La soluzione del sistema e x= (2, 6, 7)T

[CJ|q] =

0 2 2 01 0 1 12 2 0 1

19

Si assuma x(0) =q = (0, 1, 1)T Prima iterata

x(1) =

0 2 21 0 12 2 0

01

1

+

01

1

=

401

Seconda iterata

x(2) =

0 2 21 0 1

2 2 0

40

1

+

01

1

=

26

7

Il metodo di Jacobi non solo converge, ma con la scelta fatta per x(0)

raggiunge la soluzione in un numero finito di iterate.

20


91/264

In nessuna delle norme canoniche viste (norma 1,...,F) risulta verificata

CJ


92/264



(A.A. 2006-2007)


Argomenti

Metodo di Gauss-Seidel

Matrici definite positive e C.S. di convergenza

Esercizi sui metodi iterativi (Jacobi, Gauss-Seidel)

Metodi diretti.

Algoritmi per i metodi di sostituzione

Il metodo di eliminazioni di Gauss

Tecnica del pivot

1

Il metodo di Gauss SeidelStruttura iterativa del metodo di Jacobi:

x(k+1)i =

ai1

aiix

(k)1

ai,i1aii

x(k)i1

ai,i+1aii

x(k)i+1 ainaiix(k)n + bi

aii, i= 1, 2, , n

Perche non usare via via le componenti gia aggiornate?

x(k+1)i =

ai1

aiix

(k+1)1

ai,i1aii

x((k+1)i1

ai,i+1

aii

x(k)i+1

ain

aii

x(k)n +

bi

aii

, i= 1, 2,

, n

2


93/264

In forma matriciale: A= L + D+ U

xk+1=D1

(L)xk+1+ D1

(U)xk+ D1

b

da cui

xk+1= (D+ L)1Uxk+ (D+ L)1b

CGS= (D+ L)1U

Il nuovo metodo e diverso dal metodo di Jacobi enon si puo considerare

un suo miglioramento.

3

Esempio Sistema (v. lezione precedente - m. di Jacobi)

[A|b] = 1 2 2 01 1 1 1

2 2 1 1

Utilizziamo la C.N.e S. (CGS)< 1

D+ L=

1 0 01 1 02 2 1

(D+ L)1 = 1 0 01 1 0

4 2 1

CGS=

1 0 01 1 04 2 1

0 2 20 0 10 0 0

=

0 2 20 2 10

8 6

4


94/264

Gli autovalori della matrice di iterazione di Gauss-Seidel sono

1= 0, 2= 2

2, 1 = 2 +

2 (CGS) = 2 +

2> 1

Il m. di Jacobi converge

Il metodo di Gauss Seidel non converge

5

Matrici definite positive

Definizione: A Rnn e definita positiva se:

1. A= AT (A simmetrica)

2. xTAx= i,j aijxixj >0, x

= 0

Teorema Criterio di Sylvester: A simmetrica e definita positiva Mk >0, k = 1, , n

Mk = det[aij], i , j = 1, 2, , k minori principali di testa di A:

Teorema Proprieta degli autovalori: A simmetrica e definita positiva

i suoi autovalori sono tutti reali e positivi.

una matrice definita positiva e regolare (invertibile)

6


95/264

C. S. di convergenza del metodo di Gauss-Seidel

Teorema 1 Se A e (simmetrica) definita postiva, il metodo di Gauss-

Seidel converge per ogni scelta di x0

Teorema 2SeAe diagonalmente dominante(per righee/oper colonne),

il metodo di Gauss-Seidel converge per ogni scelta di x0

ovvero:

le condizioni di diagonale dominanza sono sufficienti per la conver-

genza sia del m. di Jacobi che del m. di Gauss-Seidel

7

Esercizio Sia assegnato il sistema lineare Ax= b dipendente dal para-

metro reale

A=

3 1 01 1 20 2 2

1. Per quali valori di sussistono le condizioni sufficienti di convergenza

del metodo iterativo di Gauss-Seidel?2. Per i valori di determinati al punto precedente sono verificate il

metodo di Jacobi risulta convergente?

Soluzione 1. A = AT, non diagonalmente dominante (vedi ultima

riga/colonna).

Dal criterio di Sylvester:

M1= 3> 0, M2 = 3(1 ) 1> 0 0 < 4

3

8


96/264

Il metodo di Gauss-Seidel converge per < 43

(C.S.)

2. Occorre verificare se (CJ

)< 1

CJ =

0 13 0

11 0 210 1 0

det(CJ I) = 3 + 73(1)Autovalori: 0,

7

3(1 )

(CJ) = 73(1 )


97/264

Nota C.N.e S. di convergenza al m. di Gauss-Seidel

(D+ L)1

U =

13 0 0

1

3(1)1

1 0

13(1)

11

12

0 1 0

0 0 2

0 0 0

CGS=

0 13 0

0 13(1) 21

0 13(1

)

21

autovalori 0, 0, 73(1)

Pertanto linsiemi degli per cui il m. di Gauss-Seidel converge coin-

cide con linsieme di convergenza del m. di Jacobi.

11

I metodi diretti

Trasformiano, con un numero finito di passi, Ax= b in

Ax=bsistema equivalente

A strutturata (triangolare, diagonale, prodotto di matrici triango-lari).

12


98/264

Soluzione di sistemi triangolari

Sistemi triangolari superiori: Ux=b

u11x1 + u12x2 + + u1nxn = b1

ukkxk + + uknxn = bk

unnxn = bnAlgoritmo di sostituzione allindietro:

xn= bnunn

xk = 1ukk

bk ni=k+1

ukixi , k = n 1, n 2, ..., 2, 1

13

Sistemi triangolari inferiori: Lx=b

11x1 = b1

11

x1

+

+kk

xk

= bk

n1x1 + n2x2 + + nnxn = bnAlgoritmo di sostituzione in avanti:

x1= b111

xk = 1kk

bk

k1i=1

kixi

, k = 2, , n

14


99/264

Metodo di eliminazione di GaussA= U triangolare superiorePasso ksimo A(k1) A(k)

u11 u12 unn

0 u22 u2n0 0

. . . u3n

0 . . . 0 0 kk kn0 0 k+1,k k+1,n

0 0 ... 0 0 n,k n,n

u11 u12 unn

0 u22 u2n0 0

. . . u3n

0 . . . 0 0 ukk uk,k+1 kn0 0 0

k+1,k+1

k+1,n

0 0 ... 0 0 0

n,k+1 n,n

riga i-ma riga i-ma - ikkk

riga k-ma, i= k + 1, , n

15

Ipotesi: kk= 0 ad ogni passo k( kk valore attuale dellelemento della matrice)

Se kk = 0 si deve fare uno scambio di righe

(jk = 0, j =k n A non e regolare).Teorema SeMk= 0 k = 1, 2, , n, Mk i minori principali di testadella matriceA Rnn, allora il m. di eliminazioni di Gauss puo essereeseguito senza scambi di righe.

Teorema Se A e regolare e sempre possibile con scambi di righe su A

ottenere una matrice che soddisfa le ipotesi precedenti.

16


100/264

Ridurre la propagazione degli errori: la tecnica del pivot

Ad ogni passo, simassimizzail denominatore kk (kk elemento attuale

della matrice)

Ricerca del pivot al passo k

si determina lelemento rk t.c. |rk| = maxi=k, ,n

|ik|

a) se|rk| det(A) 0 e il sistema non e numericamenterisolubile

b)|rk| > si scambia la riga k-ma con la riga r-ma

17

Il passo ksimo in simboli matricialiA(k) =L(k)A(k1)

L(k) =

1 0 0

0 1 0

0 . . . 0

0 1 0 00 k+1,kk,k 1 0

0 ... 0 0

0 n,kk,k

0 1

Le eliminazioni in simboli matriciali

U =L(n1) L(2)L(1)A=A(Servono n 1 passi) triangolare inferiore a diagonale unitaria det = 1

18


101/264


Prof. M.L. Lo Cascio(A.A. 2006-2007)


Argomenti

Metodi diretti per sistemi lineari

Fattorizzazione LU

Fattorizzazione di Choleski

Sistemi simultanei, calcolo del determinante e dellinversa

Il condizionamento del problema

Indice di condizionamento

Determinante normalizzato

Esercizi

1

Fattorizzazione LU

Il metodo di eliminazioni di Gauss in simboli matriciali

U =L(n1) L(2)L(1)A=A

triangolare inferiore a diagonale unitaria det = 1U triangolare superiore a diagonale ukk= 0

1 =L

Teorema Una matrice regolare A ammette una sola fattorizzazione

A = LU nel prodotto di una matrice triangolare inferiore a diagonale

unitaria e di una matrice triangolare superiore ( Fattorizzazione diBanachiewicz-Doolittle )

2


102/264

Teorema

i1 A= AT definita positiva

t1 ! L triangolare inferiore t.c.A = L LT (fattorizzazione di Cholesky)

Dim. Per il teorema precedente A= LU

Siano: D = diag(u11, , unn), ukk elementi diagonali di U, B t.c.U =D B (B a diagonale unitaria).

A= LDB, AT =B TDTLT, ma (i1) A= AT, D=DT

A= LDB =BTDLT L= B T (unicita della fattorizzazione)

Osserviamo che:

Mk =

k1 ujj

(i2) ujj >0

j (criterio di Sylvester),

si puo porre D =

D

D.

A= LDDLT = (LD)(LD)T.

3

Il caso delle matrici tridiagonali

A=

d1 s1 0 0a2 d2 s2 0

. . . . . . . . .

0 0 dn1 sn10 0 an dn

=

=

1 0 0 02 1 0 00 2 1 0... ... . . . . . . ...0 0

n 1

u1 v1 0 00 u2 v2 0... ... . . . . . . ...0 0 un1 vn10 0

0 un

=LU

4


103/264

Algoritmo di Thomas

vi=si i= 1, 2, . . . , n 1

u1=d1

i=ai/ui1 i= 2, 3, . . . , n

ui=di ivi1 i= 2, 3, . . . , n

Algoritmi di sostituzione

y1= b1 yi= b

i

iy

i1 i= 2, 3, . . . , n

xn=yn/un xi= (yi vixi+1)/ui i= n 1, . . . , 1

5

Costo computazionale di un algoritmo

Algoritmi finiti, in particolare nel caso dei m. diretti per la soluzione

numerica del sistema Ax= b si introduce

Costo computazionale:

Cc = numero di moltiplicazioni e divisioni

Cc dellalgoritmo di sostituzione allindietro

Ad ogni passo ci sono

n k moltiplicazioni (n. degli addendi), k = 1, , n 1 divisioneCc =

n

k=1(n k+ 1) =n

k=1 k =(n + 1)

2

n= n2

2

Verificare che lalgoritmo delle sostituzioni in avanti ha lo stesso Cc

6


104/264

Costi computazionali

Eliminazione di Gauss e Fattorizzazione Cc=

n3

3 + O(n2

).

Fattorizzzazione di Choleski Cc =n3

6 + O(n2).

Sostituzioni (in avanti o allindietro) Cc=n2

2 + O(n).

Esercizio: Calcolare il costo computazionale dellalgoritmo di Thomas

7

Soluzione di sistemi lineari simultanei

sistemi lineari simultanei= sistemi aventi stessa matricedei coefficienti

e diversi vettori termini noti

Axi = bi i= 1, 2, . . . , r

Si fattorizza A= LU

Si risolvono, per ogni vettore bi, i due sistemi triangolari

Lyi =bii= 1, 2, . . . , r

Uxi=yi

Costo computazionale

Cc= n3

3 + O(n2)

1 fattorizzazione+r 2

n2

2 + O(n)

sostituzioni

Cc n3

3 + r n2

8


105/264

Calcolo del determinante di A

detA= det(L U) = detL =1

detU =n

k=1

ukk

Se durante la fattorizzazione sono stati fatti s scambi di righe, allora

detA= (1)sn

k=1

ukk

Il calcolo del determinante ha il costo computazionale della fattoriz-

zazione (o delle eliminazioni di Gauss): Cc n3

3

9

Calcolo dellinversa di A

A e regolare, la matrice inversa A1 e per definizioneA A1 =I I : matrice identita

I:=

1 0 00 1 0

0 0 1

= [e1 e2 en]

ei = [0, 0, , 1i

, 0, , 0]T Vettori della base canonica

A A1 =I A xi =ei i= 1, 2, , nA1 = [x1 x2 xn]

Il calcolo dellinversa equivale alla soluzione di n sistemi simultanei

Cc n3

3 + n n2 4n

3

3

10


106/264

Perche no il metodo di Cramer

Il m. di Cramer e un metodo diretto

xi = DidetA

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

n + 1 determinanti + n divisioni (trascurabili)

Cc = (n + 1)n3

3 n

4

3 + O(n3)

Luso del m. di Cramer equivale a risolveren +1 sistemi (indipendenti)

11

Indice di condizionamento di una matrice

Sistema assegnato: Ax= b, xv soluzione

Sistema perturbato: Ax= b + b, xv + x soluzione

A(xv + x) =b + b A x= b x= A1 b

x A1 bAxv =b b A xv 1xv A

1

b

e infine xxv A A1

bb

Indice di condizionamento

12


107/264

Proprieta

cond(A) = A A1 [1, +)

A A1 AA1 = I = 1 (norma indotta)Se detA 0 (A non invertibile), gli elementi di A1 e quindiA1tendono ad.

Norma di Hilbert e matrici simmetriche A= AT A1 =

A1T

A2=

(ATA) =

(A2) =

2(A) =(A) =max

A1

2= (A1)2= 2(A1) =(A1) =

1min

min autovalore di modulo minimo: cond2(A) = maxmin

13

Determinante normalizzato di una matrice

DN = detA

a12 an2ai = (ai1, , ain)T (riga i-ma di A)

Proprieta

E invariante rispetto alle moltiplicazioni delle righe di A per costantinon nulle

DN [1, 1]

DN = 0 (condizionamento pessimo) detA= 0

|DN| = 1 (condizionamento ottimo) per esempio quando i vettori aisono ortogonali.

14


108/264


109/264

a) cond(A) = A A1 in una norma compatibile.

Norma 1 (Norma)A = 1 + 12+ 13, A1 = 36 + 192 + 180

cond(A) = 748

Norma FA2F= 1 + 2/4 + 2/9 + 2/1 6 + 1/25

A12F = 9 + 1922 + 1802 + 2 362 + 2 302 + 2 1802

cond(A) = 526.02

b)DN(A) = det(A)a12a22a32a12 = 1 + 14+ 19, a22 = 14+ 19+ 116 a32 = 19+ 116+ 125DN(A) = 0.13e 2; 1DN = 758.05

17

Esercizio 2 Sia assegnata la matrice

A=

1 11 2 1

1 1 2

, R.

a)Individuare per quali valori di la matrice ammette la fattorizzazione

di Cholesky.

b) Per i valori di diversi da quelli individuati al punto a) la matrice

ammette sempre una fattorizzazzione?

c)Per quali valori di A puo essere fattorizzata senza scambi di righe?

18


110/264

a) C.N. per applicare Cholesky: A definita positiva. A e simmetrica;

C.N.e S. di Sylvester:

M1 = > 0, M2 = 2 1 > 0 > 1/2, M3 = 3 2 > 0 >2/3 A e definita positiva e quindi ammette la fattorizzazione diCholesky >2/3.

b) A e fattorizzabile, in generale nella forma LU se det(A)= 0 = 2/3.

c)C.S. per eseguire la fattorizzazzione senza scambi di righe eMk= 0e quindi (v. punto a)) oltre a

= 2/3 (matrice singolare)

= 0, = 1/2

19

Esercizio 3 Si consideri il sistema Ax= b con

A=

1 1 11 1.01 1

1.01 0.99 1

, b= b1= (2, 2, 2)T, b= b2= (2.01, 2, 2)T

la cui soluzione e x1 = (0, 0, 2)T

e x2 = (2, 1, 5.01)T

, rispettiva-mente

a)Si stimi lindice di condizionamento di A utilizzando le informazioni

disponibili.

b) Sapendo che A1 =

200 100 100100 100 0

301 200 100

stimare lindice di con-

dizionamento di A mediante la definizione e confrontare il risultato con

quello ottenuto al punto a).

20


111/264

a) Daxx cond(A)

bb cond(A)

x/xb/b

b= (0.01, 0, 0)T

, x= (2, 1, 3.01)T

Norma cond(A) 3.012

20.01

cond(A) 301

b)A= 3.01;A1= 601 cond(A) = 1809.01La stima ottenuta al punto a) suggerisce che la matrice e mal con-

dizionata, pur fornendo un valore molto inferiore al valor

Date post:	04-Jun-2018
Category:	Documents
Upload:	sascia1775
View:	216 times
Download:	1 times

Calcolo_Numerico

Documents