Mechanical Engineering - Metodi Analitici e Numerici per L'Ingegneria

Richiami di analisi e calcolo numerico

Other

Perr >0 numero reale ex0 2Rn ¯ssati, indicheremo conB r(x0 ) l'insieme dei puntix2Rn che distano dax0 meno dirnel senso della norma euclidea (si ricorda che perx= (x 1; x 2; :::; x n)2Rn la norma euclidea dix, indicata conjxjo conkxk µe data dal numero reale positivo¡ x2 1+x2 2+:::+x2 n¢ 1=2 e la distanza tra due punti xeyµekx¡yk). In simboli B r(x0 ) =n x2Rn¯ ¯ ¯kx¡x 0k< ro Indicheremo inoltre con¯ ¯ B r(x0 )¯ ¯ la misura diB r(x0 ), che in particolare per n= 1 µe 2r, pern= 2 µe l'area del cerchio di raggiore pern= 3 µe il volume della sfera di raggior. Teorema 1.1. ii) Sefµe continua inA£Binsieme con la derivata parziale@ f @ x iper qualche 1·i·n, al lora (2.1)@ F @ x i(x) =Z B@ f @ x i(x;y)dy: 3.Limiti, derivate e integrali di serie di funzioni Consideriamo una serie di funzioni+1 X n=0f n(x), conf n:A!R. Diremo che la serie soddisfa ilcriterio di WeierstrassinAse esistono numeri realia ntali che jf n(x)j< a nperx2Aen2N, e se+1 X n=0a nconverge. Una serie di funzioni che soddisfa il criterio di Weierstass inAin particolare converge per tutti glix2Ae quindi de¯nisce inAuna funzione sommaS(x). Teorema 3.1. Metodi Analitici e Numerici per l’ingegneria Cerutti Maria Cristina A cura di: Andrea Fuso, Gianpiero Gaeta Indice I Richiami di Analisi 2 1 Funzioni 2 2 Vettori 2 3 Matrici 3 3.1 Algebra delle matrici . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 4 Teorema spettrale 5 II Calcolo numerico 6 5 Risoluzione di sistemi lineari con metodi diretti 6 5.1 Metodo di Cramer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 5.2 Risoluzione sistemi triangolari . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 5.2.1 Backward substitution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 5.2.2 Forward substitution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 5.3 Metodo di eliminazione di Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 5.4 Fattorizzazione LU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 5.4.1 Pivoting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 5.5 Fattorizzazione di Cholesky . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 5.6 Errore e condizionamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 5.6.1 Numero di condizionamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 1 Parte I Richiami di Analisi 1 Funzioni Teorema 1 (Teorema di Weierstrass) .Ogni funzione continua f :K æ R,K µ Rn, insieme chiuso e limitato (compatto), ammette massimo e minimo. Teorema 2 (Teorema dei valori intermedi o teorema di Darboux) .Sia f:[a, b ]æ R una funzione continua su un intervallo chiuso e limitato [a, b ].Detti m eM rispettivamente il minimo e massimo valore di f(x)nell’intervallo [a, b ],risultache fassume tutti i valori compresi tra m eM (poichè è una funzione continua). ’cœ[a, b ]( Con m Æ cÆ M )÷x0œ[a, b ]t.c. f(x0)= c In generale per una funzione continua f:C æ R,C ™ Rn, insieme connesso, se x1,x2œC,fassume tutti i valori compresi tra f(x1)ef(x2). Teorema 3 (Teorema della media integrale) .Sia f:[a, b ]æ R una funzione continua (e quindi integrabile), allora esiste cœ[a, b ]tale che: 1 b≠a ⁄b a f(x)dx = f(c) ∆ ⁄b a f(x)dx =( b≠a)f(c) Osservazione 1.1 .Essendo l’integrale collegato al calcolo dell’area sottesa dall funzione, il teorema della media integrale a erma che esiste un punto cappartenente all’intervallo, tale per cui l’integrale è calcolabile come l’area di un rettangolo equivalente, avente per base la dimensione dell’intervallo e per altezza l’immagine di c. Teorema 4 (Teorema fondamentale dell’algebra) .Ogni polinomio a coe cienti reali o complessi, di grado maggiore o uguale a 1, ammette almeno una radice complessa. Quindi ogni polinomio a coe cienti reali o complessi di grado nØ 1ammente nel campo complesso C esattamente nsoluzioni (contate con la loro molteplicità). Definizione 1 (Distanza punto-retta) .Si definisce distanza di un punto P =( xP,y P)da una retta rdi equazione y= mx +q, la quantità. d(P, r )= |yP≠(mx P+q)| Ô1+ m2 oppure scrivendo la retta in forma implicita ax +by +c=0 si ha: d(P, r )= |ax P+byP+c| Ôa2+b2 2Vettori Teorema 5 (Disuguaglianza di Schwarz) .Per la disuguaglianza di Schwarz il valore assoluto del prodotto scalare di due vettori xeyè minore o uguale al prodotto delle loro norme: Îx·yÎÆÎ xηÎyÎ 2 3 Matrici 3.1 Algebra delle matrici Definizione 2 (Matrice diagonale) .Una matrice diagonale è una matrice diagonale D =[ dij]che ha tutti gli elementi al di fuori dalla diagonale principale nulli: dij=0 se i”= j Proposizione 3.1. Il prodotto di due matrici diagonali è ancora una matrice diagonale. Definizione 3 (Matrice inversa) .Siano A eB due matrici quadrate di ordine n.Sidiceche B è la matrica inversa di A se: AB = BA = I Dove Iè la matrice identità di ordine n. Se esiste una matrice inversa di A, allora si dice che A èinvertibile. Teorema 6 (Condizioni di invertibilità) .Per una matrice A quadrata di ordine n le seguenti condizioni sono equivalenti: 1. A ha rango massimo: r(A)= n 2. L’unica soluzione di Ax= 0èx= 0: Ker (A)= )0* 3. A èinvertibile. 4. A ha un’inversa sinistra. 5. A ha un’inversa destra. Definizione 4 (Matrice non singolare) .Una matrice A quadrata di ordine nche soddisfi le condizioni equivalenti del teorema 6 si dice non singolare . Quindi non singolare è sinonimo di invertibile, e singolare significa non invertibile. Proposizione 3.2 (Inversa di un prodotto) .Siano A eB due matrici quadrate di ordine n. La matrice prodotto AB è invertibile se e solo se A eB sono invertibili. In tal caso: (AB )≠1= B≠1A≠1 Definizione 5 (Matrice simmetrica) .Una matrice A si dice simmetrica se AT = A . Una matrice A si dice antissimetrica se AT= ≠A . Esempio 1. Per una matrice quadrata di ordine 3: Forma matrice simmetrica æ S U ab c bde cef T V 3 Forma matrice antissimetrica æ S U 0 bc ≠b 0 e ≠c ≠e 0 T V Esempio 2 (Matrice di Hilbert) .Una matrice H è una matrice di Hilbert se ha gli elementi hij pari a: hij=( i+j≠1)≠1 H = S WWWWWWWWWWWWWU 1 12 13 ... 1n 12 13 14 ... 1n+1 13 14 15 ... 1n+2 ... ... ... ... 1n 1n+1 1n+2 ... 12n≠1 T XXXXXXXXXXXXXV Proposizione 3.3. Una matrice simmetrica è sempre diagonalizzabile con autovalorei reali. Definizione 6 (Matrice definita positiva) .Sia A una matrice simmetrica di ordine na coe cienti reali. La matrice A si dice definita positiva se per ogni vettore x=[ x1,...,x n]œ Rn, il prodotto xA xTè maggiore di zero, ovvero: A definita positiva … xA xT> 0 ’x œ Rn Con x”= 0 Oppure: vedi teorema 7 Osservazione 3.1. Gli autovalori di una matrice definita positiva sono tutti positivi. Proposizione 3.4 (Proprietà matrice trasposta) . (A+B)T= AT+BT (tA )T= tA T !AT"T= A (AB )T= BTAT !AT"≠1= !A≠1"T Se A èsimmetrica æ AT= A Se A è antisimmetrica æ AT= ≠A Teorema 7 (Matrice definita positiva) .Per una matrice A simmetrica le seguenti condizioni sono equivalenti: 1. A è definita positiva. 2. I minori principali di nord-ovest di A sono positivi. 4 Definizione 7 (Matrice a dominanza diagonale per righe) .Una matrice A si dice che è dominante per righe se è una matrice quadrata di ordine ncon gli elementi sulla diagonale principale maggiori o uguali (in valore assoluto) della somma di tutti i restanti elementi della stessa riga (sempre in valore assoluto): Dominanza diagonale in senso debole æ |aii|Ø nÿ i=1 ,j”=i |aij| Dominanza diagonale in senso stretto æ |aii|> nÿ i=1 ,j”=i |aij| 4Teoremaspettrale Proposizione 4.1. Sia A una matrice reale simmetrica (Def: 5) di ordine nesia q(x)= xTAxallora: 1. A è diagonalizzabile con autovalori reali. 2. Se Av= ⁄v, allora q(v)= ⁄ÎvÎ2(quindi se vè autovettore della matrice A ). 3. Se ⁄min e⁄max sono il minimo e massimo degli autovalori di A, allora: ⁄min ÎxÎ2Æ xTAxÆ ⁄max ÎxÎ2 5 Parte I I Calcolo numerico 5 Risoluzione di sistemi lineari con metodi diretti 5.1 Metodo di Cramer Teorema 1. (Teo re m a d i C ra m e r ) Sia Ax= bun sistema lineare quadrato di n equazioni in n incognite. Se det (A)”=0 il sistema ammette un’unica soluzione v=[ x1,..., x n]Tdi componenti: xi= det (Ai) det (A) Dove Aiè ottenuta sostituendo la colonna di A con il termine noto b Esempio 1. I x+2 y=3 3x+7 y=2 ∆ A = 5 12 37 6 b= 5 3 2 6 x= det (Ax) det (A) = det 5 32 27 6 det 5 12 37 6= det 5 b1 2 b2 7 6 det 5 12 37 6 = 17 y= det (Ay) det (A) = det 5 13 32 6 det 5 12 37 6= det 5 1 b1 3 b2 6 det 5 12 37 6 = ≠7 Costo computazionale 1. Il metodo di Cramer richiede un costo computazionale pari a: #flops =2( n+ 1)! calcolando i determinanti con la regola di Laplace. 5.2 Risoluzione sistemi triangolari 5.2.1 Backward substitution Un sistema triangolare superiore, è un sistema che in forma matriciale si presenta nel seguente modo: S WWWU a1,1 a1,2 ... a 1,n 0 a2,2 ... a 2,n ... ... ... ... 00 ... a n,n T XXXV S WWWU x1 x2... xn T XXXV= S WWWU b1 b2... bn T XXXV Il termine xnè immediatamente calcolabile come: xn= bn an,n Per quanto riguarda il termine xn≠1l’operazione da eseguire corrisponde a: xn≠1=[ bn≠1≠an≠1,nxn] 1 an≠1,n≠1 Si può quindi estrapolare l’algoritmo per il calcolo del generico xi: 6 Proposizione 5.1. Data una matrice triangolare alta, è possibile trovare tutte le soluzioni attraverso il seguente algoritmo: xn= bn an,n xi= 1 ai,i S Ubi≠ nÿ j=i+1 ai,jxj T V i=1 ,...,n j= i+1 ,...n Questo algoritmo prende il nome di sostituzione all’indietro , conosciuto anche come backward substitu- tion . Costo computazionale 2. Per quanto riguarda il costo computazionale dell’algoritmo si distinguono due casi: il primo passaggio e i successivi. Per il primo passaggio (calcolo dell’elemento xn) si esegue una sola operazione (la divisione), mentre nel secondo passaggio le operazioni diventano tre: la divisione per ai,i , la moltiplicazione e la sottrazione all’interno della parentesi. Quindi il numero di operazioni elementari risulta essere: #flops = 1 + (1 + 2) + ··· +(1+2 i)+ ··· +(1+2( n≠1)) = nÿ i=0 (1 + 2 i) » #flops = n+2 n≠1 ÿ i=1 i= n+A2n(n≠1) A2 = n2 5.2.2 Forward substitution Un sistema triangolare inferiore , è un sistema che in forma matriciale si presenta nel seguente modo: S WWWU a1,1 0 ... 0 a1,2 a2,2 ... 0 ... ... ... ... an,1 an,2 ... a n,n T XXXV S WWWU x1 x2... xn T XXXV= S WWWU b1 b2... bn T XXXV Il termine x1è direttamente calcolabile come: x1= b1 a1,1 Il secondo termine risulta essere pari a: x2= 1 a2,2[b2≠a2,1x1] Algoritmo 1. L’algoritmo per il calcolo del la generica xirisulta quindi essere: xi= 1 ai,i S Ubi≠ i≠1ÿ j=1 ai,jxj T V j=1 ,...,n Proposizione 5.2. Questo algoritmo per calcolare le xiincognite prende il nome di algoritmo di sostituzione in avanti , conosciuto anche come forward substitution . Osservazione 5.1. Il costo computazionale per il forward substitution è perfettamente analogo a quello del backward substitution, è quindi pari a n2 Osservazione 5.2. Se A è triangolare non singolare (detA ”= 0) , gli algoritmi funzionano sempre poichè: detA ”=0 … ai,i ”=0 ’ i=1 ,...,n 7 5.3 Metodo di eliminazione di Gauss Esiste un algoritmo, detto algoritmo di eliminazione di Gauss, che riduce una matrice in una matrice a scala mediante un numero finito di operazioni elementari sulle righe. Il metodo di eliminazione gaussiana si basa sull’idea di ridurre il sistema Ax= bad un sistema equivalente (avente cioè la stessa soluzione) della forma Ux= ˆbdove U è triangolare superiore e ˆbè un nuovo termine noto. Quest’ultimo sistema potrà essere risolto con il metodo delle sostituzioni all’indietro ( backward substitution Prop. 5.1). Funzionamento del metodo Indicando il sistema originale come A(1) x= b(1) Si sostituisce la seconda riga con la di erenza tra la riga stessa e la prima riga moltiplicata per una costante non nulla, scelto in modo tale da rendere nullo il primo elemento della seconda riga; in questo modo si ottiene un sistema equivalente a quello di partenza æ cioè con la stessa soluzione. Tale costante è il motiplicatore m2,1: m2,1= a2,1 a1,1 b(2)2 = b(1)2 ≠m2,1b(1)1 tale procedimento viene ripetuto per tutte le righe successive. Facendo questo si ottiene, per la prima colonna, un vettore che ha come unico elemento non nullo quello sulla diagonale. S WWWU a1,1 a1,2 ... a 1,n 0 a2,2 ... a 2,n ... ... ... ... 0 an,2 ... a n,n --------- b1 b2... bn T XXXV∆ S WWWWU a(1)1,1 a(1)1,2 ... a (1)1,n 0 a(2)2,2 ... a (2)2,n ... ... ... ... 0 a(2)n,2 ... a (2)n,n T XXXXV S WWWU x1 x2... xn T XXXV= S WWWWU b(1)1 b(2)2... b(2)n T XXXXV Dove gli elementi indicati con: (2) corrispondono al risultato delle operazioni svolte sulle righe, dopo l’operazione di MEG sulla prima colonna. In generale si ha che: mi,k = a(k)i,k a(k)k,k i= k+1 ,...,n Dove gli elementi a(k)k,k sono detti pivot . Quindi nel caso generale, l’operazione eseguita è: a(k+1)i,j = a(k)i,j ≠mi,k a(k)k,j i, j = k+1 ,...,n b(k+1)i = b(k)i ≠m(k)i,k b(k)k i= k+1 ,...,n Osservazione 5.3. L’operazione di MEG sia arresta quando incontra un pivot nullo, per evitare questo si utilizza il pivoting (Thm 5.4.1). Esempio 3. Si consideri il seguente sistema: Y_] _[ 3x≠y+z=2 2x+y=1 ≠2x≠2y+z= ≠2 =∆ S WWWWWU 3 ≠11 210 ≠2 ≠21 ¸ ˚˙ ˝ A ------ 2 1 ≠1 ¸˚˙ ˝ ˛b T XXXXXV = B 8 Con i passaggi successivi si vuole rendere A (matrice dei coe cienti) triangolare alta, cioè tale da avere i termini sotto la diagonale principale tutti identicamente nulli, per risolvere così un sistema triangolare. Le operazioni da eseguire sono delle particolari combinazioni linari delle righe. A = S U a1 a2 a3 T V Con airighe di A B = Ë A ---˛b È = S U b1 b2 b3 T V Si costruiscono quindi successive matrici B(k), equivalenti a B, in cui gli elementi sotto la diagonale principale diventano progressivamente uguali a zero, e sia: B(1) = B B(2) è equivalente a B(1) ed ha gli elementi sottodiagonali della prima colonna nulli; in questo caso si cerca di avere a(2)21 = a(2)31 =0 ,quindi: b(2)1 = b(2)2 b(2)2 = b(1)2 ≠m21b(1)1 m21 è quel valore che permette di ottenere a(2)21 =0 : #2101 $≠ 2 3 #3 ≠112 $= #0 53 ≠23 ≠13 $ ∆ m21 = 2 3= a(1)21 a(1)11 pre rendere a(2)31 =0 : b(2)3 = b(1)3 ≠m31b(1)1 #≠2 ≠21 ≠2$+ 2 3 #3 ≠112 $= #0 ≠83 53 ≠23 $ ∆ m31 = ≠2 3= a(1)31 a(1)11 La matrice B(2) risulta quindi essere: B(2) = S WWWWU 3 ≠11 2 0 53 ≠23 ≠13 0 ≠83 53 ≠23 T XXXXV Ora si cerca la matrice B(3) tale da avere a(3)32 =0 ,quindi: b(3)1 = b(2)1 b(3)2 = b(2)2 b(3)3 = b(2)3 ≠m32b(2)2 ∆ #0 ≠83 53 ≠23 $+ 8 5 #0 53 ≠23 ≠13 $= #00 35 ≠65 $ 9 ∆ m32 = ≠8 5= a(2)32 a(2)22 La matrice B(3) risulta essere: B(3) = S WWWWU 3 ≠11 2 0 53 ≠23 ≠13 00 35 ≠65 T XXXXV=∆ Y______] ______[ 3x≠y+z=2 53y≠ 23z= ≠13 35z= ≠65 ∆ x= S U 1 ≠1 ≠2 T V Proposizione 5.3. Valgono le seguenti formule per la somma dei primi N interi, e per la somma dei quadrati dei primi N interi: Nÿ i=1 i= N (N + 1) 2 Nÿ i=1 i2= (2N + 1) ( N + 1) N 6 Costo computazionale 3. Contributo per il calcolo di m Essendo il generico mi,k pari a: mi,k = a(k)i,k a(k)k,k i= k+1 ,...,n Segue che si esegue una sola operazione (la divisione) per n≠(k≠1) volte, per ciascuna delle righe della matrice successiva alla prima ( k=2 ,...,n ). #flops =( n≠1) + ( n≠2) + ··· +1 = n≠1 ÿ i=1 i= (n≠1) ( n≠1 + 1) 2 = n(n≠1) 2 Contributo per il calcolo di ai,j Essendo il generico a(k+1)i,j = a(k)i,j ≠mi,k a(k)k,j , segue che per e ettuare tale calcolo si devono eseguire due operazioni elementari: il prodotto mi,k a(k)k,j e la sottrazione. Al primo passaggio bisogna calcolare gli elementi di una matrice (n≠1) ◊ (n≠1) (la prima riga della matrice A(2) , uguale alla prima riga della matrice A(1) , non si calcola così come la prima colonna sottodiagonale di A(2) , che ha elementi uguali a zero). Ad ogni passaggio bisogna calcolare gli elementi di una matrice quadrata (n≠(k≠1)) . Di conseguenza il numero di operazioni elementari totali (flops) corrisponde a: #flops =2 Ë (n≠1)2+( n≠2)2+··· +1 È =2 n≠1 ÿ i=1 i2=2 (2n≠2 + 1) ( n≠1 + 1) ( n≠1) 6 = (2n≠1) ( n≠1) n 3 10 Contributo per il calcolo di b Essendo il generico b(k+1)i = b(k)i ≠ m(k)i,k b(k)k , segue che per e ettuare tale calcolo si devono eseguire due operazioni elementari: il prodotto m(k)i,k b(k)k e la sottrazione, e ettuate per n≠ 1volte, cioè per gli n≠ 1elementi presenti nel vettore b(è il pedice ia scorrere); al passaggio successivo n≠2e così via. Di conseguenza il costo computazionale di tale operazione risulta essere: #flops =2[( n≠1) + ( n≠2) + ··· + 1] = 2 n≠1 ÿ i=1 i=2 n(n≠1) 2 = n(n≠1) Costo computazionale totale In definitiva il costo computazionale totale è pari a: #flops = n(n≠1) 2 + (2n≠1) ( n≠1) n 3 +2 n(n≠1) 2 = 3 2n(n≠1) + (2n≠1) ( n≠1) n 3 ƒ 3 2n2+ 2 3n3ƒ 2 3n3 (1) Osservazione 1.L’approssimazione è valida per valori di nmolto grandi. Osservazione 2.Il costo computazione per il MEG è molto minore di quello richiesto dal metodo di Cramer: MEG Cramer Costo computazionale 23n3 2( n+ 1)! Algoritmo 2. Sia k=1 ,...,n l’indice colonna, l’algoritmo del MEG conta tre passaggi: 1. Il calcolo del coe ciente mi,k : mi,k = a(k)i,k a(k)k,k i=1 ,...,n ≠1 2. Il calcolo del generico a(k+1)i,j : a(k+1)i,j = a(k)i,j ≠mi,k a(k)k,j i=1 ,...,n ≠1 3. Il calcolo del generico b(k+1)i : b(k+1)i = b(k)i ≠m(k)i,k b(k)k i=1 ,...,n ≠1 Osservazione 5.4. È ora possibile calcolare il costo computazionale corretto del MEG, al cui costo precedentemente ( Eq: 1) calcolato si deve aggiungere il risultato appena ottenuto: #flops tot = 3 2n(n≠1) + (2n≠1) ( n≠1) n 3 +n2ƒ 2 3n3 Di conseguenza per un numero di passaggi nmolto grande, il costo computazionale del MEG dipende sempre dal termine 23n3e il contributo dell’algoritmo di backward substitution risulta perciò trascurabile (di un ordine di grandezza inferiore). 5.4 Fattorizzazione LU Cosa fa 1. Fattorizza la matrice A in modo tale da avere A = LU dove L è una matrice di tipo triangolare bassa, e U una matrice triangolare alta. Osservazione 5.5. Se esiste una fattorizzazione LU di A, per risolvere il sistema Ax= b,sipuòscrivere LU x= b, che equivale a risolvere i sistemi: I Ly= b Ux= y di cui il primo è triangolare inferiore e il secondo è triangolare superiore (Sezione: 5.2). 11 Problema 5.1. Come si trovano le matrici Led U? Esempio 4. Considerando ad esempio una matrice 2◊2: 5 a1,1 a1,2 a2,1 a2,2 6 = 5 l1,1 0 l2,1 l2,2 65 u1,1 u1,2 0 u2,2 6 Il sistema presenta quattro equazioni e sei incognite (gli elementi delle matrici L eU), di conseguenza è un sistema di tipo sottodeterminato. Per risolvere questo problema, in questo caso si possono scegliere arbitrariamente due elementi, ad esempio l1,1el2,2(gli elementi sulla diagonale della matrice L) e porli uguali a uno, nell’esempio: 5 a1,1 a1,2 a2,1 a2,2 6 = 5 10 l2,1 1 65 u1,1 u1,2 0 u2,2 6 In questo modo il sistema diventa risolvibile: Y____] ____[ u1,1= a1,1 u1,2= a1,2 l2,1= a2,1a1,1(a1,1”= 0) u2,2= a2,2≠ a2,1a1,1a1,2= a2,2a1,1≠a2,1a1,2 a1,1 = detAa1,1 Osservazione 5.6. Se si fosse applicato il MEG alla matrice A: m2,1= a2,1 a1,1= l1,2 Inoltre la prima riga di U è uguale alla prima riga di A eu2,2= a2,2≠ m2,1a1,2, cioè la stessa operazione che si avrebbe ottenuto applicando il MEG; quindi U è la matrice che si ottiene alla fine del MEG. In generale per determinare Led U si devono trovare n2+nelementi con n2equazioni; si scelgono quindi gli n elementi della diagonale di L uguali ad 1 ( fattorizzazione LU di Gauss ). Teorema 8.Se tutti i pivots sono diversi da zero, la fattorizzazione LU di Gauss è unica, e si ha che: Læ Gli elementi sottodiagonali sono i moltiplicatori del MEG U æ È la matrice finale del MEG quindi condizione necessaria e su ciente per l’esistenza di questa fattorizzazione è che tutti i pivots a(k)k,k del MEG siano diversi da zero (condizione che si verifica quando la matrice non è singolare (Def: 4). Dimostrazione. Sia A(0) = A , si ossercva che A(1) , cioè la matrice trasformata dal primo passaggio iterativo del MEG che porta all’azzeramento degli elementi della prima colonna, si può ottenere come: A(1) = M1A(0) M1= S WWWWWWU 100 ... 0 ≠m2,1 10 ... 0 ≠m3,1 01 ... ... ... ... ... ... ... ≠mn,1 00 ... 1 T XXXXXXV • Ora per passare alla matrice A(2) il procedimento è analogo: A(2) = M2A(1) = M2M1A(0) : 12 M2= S WWWWWWU 100 ... 0 010 ... 0 0 ≠m3,2 1 ... ... ... ... ... ... ... 0 ≠mn,2 0 ... 1 T XXXXXXV • Si procede in questo modo fino ad arrivare alla matrice finale triangolare superiore del MEG, che indichiamo con U. Quindi la generica matrice Mkha una forma del tipo: Mk= S WWWWWWWWU 1 ... 00 ... 0 ... ... ... ... ... 0100 0 ≠mk+1 ,k 10 ... ... ... ... ... ... 0 ... ≠mn,k 0 ... 1 T XXXXXXXXV = In≠mkeTk dove ekè il k-esimo vettore della base canonica: eTk= #00 ··· 1 ··· 0$ mentre mkcorrisponde ad un vettore del tipo: mk= S WWWWWWWWU 0 0 ... mk+1 ,k... mn,k T XXXXXXXXV Esempio 5.Prendendo ad esempio una matrice 3◊ 3,incuil’elemento m3,2=4 , la matrice M2risulta essere: M2= S U 100 010 0 ≠41 T V= I3≠m2eTk= S U 100 010 001 T V≠ Q a S U 0 0 4 T V#010 $ R b = S U 100 010 001 T V≠ S U 000 000 040 T V • Quindi alla fine del processo si ottiene la matrice Akche è la matrice ridotta a scala, e che coincide quindi con U : A(k)= U = Mn≠1Mn≠2··· M2M1A • Ricordando la proprietà 3.2, è possibile scrivere: A = M ≠11 M ≠12 ··· M ≠1n≠2M ≠1n≠1U = LU • Le matrici Mksono matrici triangolari inferiori con elementi diagonali tutti pari ad uno e con inversa data da M ≠1k = In+mkeTk come si può verificare facilmente essendo !mieTi "! mjeTj "uguale alla matrice identicamente nulla per iÆ j. Di conseguenza si ha che: A = M ≠11 ··· M ≠1n≠1U = !In+m1eT1 "··· !In+mn≠1eTn≠1 "U = A In+ n≠1 ÿ i=1 mieTi B U 13 ∆ A = S WWWWWWWU 10 ... ... 0 m21 1 ... ... m32 ... ... ... ... ... 0 mn1 mn3 ... m n,n ≠1 1 T XXXXXXXV U • Di conseguenza risulta essere: M ≠11 M ≠12 ··· M ≠1n≠2M ≠1n≠1= L Proprietà 1. Sia A œ Rn◊n. La fattorizzazione LU di A con lii =1 per i=1 ,...,n (fattorizzazione LU di Gauss), esiste ed èunica se e solo se le sottomatrici principali di nord-ovest Aidi A di ordine i=1 ,...,n ≠1sono non singolari (Def: 4). Osservazione 5.7. Nel caso in cui la fattorizzazione LU esiste, il determinante di A si può calcolare come: detA = u11u22 ··· unn = nŸ i=1 uii Quando si può usare 1. Condizione necessaria e su ciente per fattorizzazione LU Una condizione necessaria e su ciente a nchè possa esistere la fattorizzazione LU è che i determinanti delle sottomatrici principali di nord-ovest siano ”=0 , cioè che non siano singolari (Def: 4). Condizioni su cienti per fattorizzazione LU Condizione su ciente per la fattorizzazione LU èche: 1. Che la matrice A sia simmetrica (Def: 5) e definita positiva (Def: 6). Siccome A èsimmetricaperla Proposizione 4.1 vale: xTAxØ ⁄min ÎxÎ2 2. Che la matrice A sia a dominanza diagonale stretta per le righe (Def: 7). |aii|> nÿ i=1 ,j”=i |aij| 5.4.1 Pivoting Se si incontra un pivot a(k)k,k =0 il procedimento si arresta. In questo caso siccome A è non singolare, per j>k ÷ a(k)j,k ”=0 . Se così non fosse il determinante della matrice sarebbe pari a zero, e dunque la matrice sarebbe singolare, dunque contraria all’ipotesi di partenza. Scambiando quindi la riga je la riga kè possibile procedere con la fattorizzazione. Questo procedimento viene detto pivoting per righe . Poiché un valore piccolo (in valore assoluto) del pivot a(k)k,k , può amplificare gli eventuali errori di arrotondamento presenti negli elementi a(k)k,j èutilee ettuare il pivoting per righe anche se a(k)k,k ”=0 utilizzando come pivot quello maggiore. Questa operazione viene e ettuata tramite la moltiplicazione a sinsitra per la matrice di permutazione P,cioè la matrice ottenuta dalla matrice identità in cui si sono scambiate le righe corrispondenti. 14 Esempio 6. Supponendo di avere una matrice A della forma: A = S U 500 100 300 T V e di voler scambiare le righe 2 e 3 di tale matrice per ottenere la matrice: B = S U 500 300 100 T V moltiplicando A per l’oppurtuna matrice P è posibile ottenere B: P = S U 100 001 010 T V B = PA = S U 100 001 010 T V S U 500 100 300 T V= S U 500 300 100 T V Se nell’operazione di fattorizzazione avviene il pivoting, alla fine del processo avrò L ed U e la matrice di permutazione P, prodotto di tutte le permutazioni e ettuate. In questo caso, per risolvere il sistema, dovrò risolvere: PA x= Pb ∆ LU y= Pb ∆ I Ly= Pb Ux= y È anche possibile e ettuate il pivoting totale, cioè cercare l’elemento maggiore all’interno della sottomatrice A(k)k = Ó a(k)i,j Ô i,j=k,...,n e scambiare le relative righe e colonne (si veda la Figura 1). Figura 1: Pivotazione per righe (a sinistra) o totale (a destra). In grigio più scuro, le aree interessate dalla ricerca dell’elemento pivotale Il vantaggio della fattorizzazione LU rispetto al MEG, per la risoluzione di sistemi è che L ed U dipendono dalla sola A e non dal termine noto, la stessa fattorizzazione può essere utilizzata per risolvere diversi sistemi lineari con la stessa matrice A dei coe cienti, ma con termine noto bvariabile. Matlab 1. In Matlab il comando per eseguire la fattorizzazione LU su una matrice A è: [L U P]=lu (A) y=f w s u b ( L , P b); x=b k s u b ( U , y ) ; 15 5.5 Fattorizzazione di Cholesky Rappresenta un metodo diretto alternativo a LU . Definizione 8 (Metodo Cholesky) .Permette di riscrivere la matrice A come: A = QTQ QTæ Triangolare inferiore Q æ Triangolare superiore Costo computazionale 4. Il costo computazionale associato a questo metodo è la metà di quello di LU : #flops = n3 3 5.6 Errore e condizionamento Definizione 9 (Errore relativo) .Sia ˆxla soluzione numerica calcolata e sia xla soluzione esatta, si definisce errore relativo la quantità: E = Îx≠ ˆxÎ ÎxÎ Definizione 10 (Residuo) .Sia ˆxla soluzione numerica calcolata per il sistema lineare del tipo Ax= b,sidefinisce residuo la quantità: r= b≠Aˆx= Ax≠Aˆx= A(x≠ ˆx) Allora: (x≠ ˆx)= A≠1r 5.6.1 Numero di condizionamento Definizione 11 (Norma di matrice) .La norma di una matrice A è definita come: ÎAÎ= maxxœRn ÎAxÎ ÎxÎ Poiché ’xœRnpuò essere scritto come x= Îxη xÎxÎ= Îxηy,dove yè il versore che indica la direzione di x, si ha: ÎAxÎ ÎxÎ = ÎA(Îxηy)Î ÎxÎ = ÎxÎÎAyÎ ÎxÎ da cui: ÎAÎ= max yœRn ÎyÎ=1 ÎAyÎ Dalla definizione di ÎAÎsi ha: ÎAÎØ ÎAxÎ ÎxÎ ≈∆ Î AxÎÆÎ AηÎxÎ Osservazione 5.8. Se A è una matrice simmetrica (Def: 5) e definita positiva (Def: 6) allora ÎAÎ= ⁄max cioè al massimo degli autovalori. 16 Dimostrazione. .... Definizione 12 (Numero di condizionamento) .Si definisce numero di condizionamento K (A)della matrice A ,la quantità: K (A)= ÎAηÎA≠1Î • Se il numero di condizionamento è molto grande, si dice che la matrice è malcondizionata . • Se il numero di condizionamento è piccolo (vicino al valore unitario), si dice che la matrice è bencondizionata . Matlab 2. Il comando Matlab per il numero di condizionamento di una matrice A è: cond (A) Osservazione 5.9. Nel caso ideale il numero di condizionamento K (A)di una matrice A dovrebbe essere pari a 1 (Def: 3), tuttavia gli errori di approssimazione che si generano in Matlab portano ad avere K (A)”=1 Proposizione 5.4. Esiste una relazione tra l’errore relativo E e il numero di condizionamento: E Æ K (A)ÎrÎ ÎbÎ Dimostrazione. Per definizione di errore relativo (Def: 9) si ha che: E = Îx≠ ˆxÎ ÎxÎ = ÎA≠1rÎ ÎxÎ • Per la definizione di norma di matrice (Def: 11) : E = ÎA≠1rÎ ÎxÎ Æ ÎA≠1ηÎrÎ ÎxÎ • Si può inoltre notare (sempre per la definizione di norma) che: ÎAxÎ ÎxÎ Æ ÎAηÎxÎ ÎxÎ ∆ xØ ÎAxÎ ÎAÎ = ÎbÎ ÎAÎ • Quindi si ottiene: E Æ ÎA≠1rÎ ÎbÎÎAÎ = ÎAηÎA≠1ÎÎrÎ ÎbÎ= K (A)ÎrÎ ÎbÎ (2) Osservazione 5.10. Se A è una matrice simmetrica e definita positiva, l’osservazione 5.8 implica che il numero di condizionamento è: K (A)= ⁄max ⁄min dove ⁄max e⁄min sono rispettivamente il massimo e il minimo degli autovalori di A. Osservazione 5.11. In genere il residuo ÎrÎsi mantiene piccolo, tuttavia può capitare che K (A)sia tanto grande da non essere compensato dal residuo. Esempio 7. Prendendo per esempio la matrice di Hilbert (Def: 2) di ordine 4: K (H4)> 1.5◊104 e Îr΃ 10≠16 • Passaggi tipici ÎAxÎ2=( Ax)T(Ax)= nÿ i=1 (xi)2(⁄i)2Æ ... ÎAxÎ2=( Ax)T(Ax) Ma ogni vettore xin funzione degli autovettori 17 6 Risoluzione di sistemi lineari con metodi iterativi I metodi iterativi per risolvere sistemi lineare del tipo Ax= bconsistono nel costruire una successione del tipo )x(k)* kœNtale che x(k)converga alla soluzione xper kche tende all’infinito, per qualunque vettore iniziale x(0) , scelto arbitrariamente, cioè: Ó x(k)Ô kœN t.c. x(k)kæŒ≠æ x limkæŒ x(k)= x Questi metodi si basano sostanzialmente su tre passi fondamentali: 1. Determinazione del metodo, cioè come determinare x(k). 2. Dimostrazione della convergenza del metodo. 3. Stima dell’errore per determinare il criterio di arresto. 6.1 Metodo del gradiente Quando si può usare 2. Quando A è simmetrica e definita positiva, si può usare perchè il metodo del gradiente converge. A cosa serve 1. Trova una soluzione ad un sistema lineare del tipo Ax=b Cosa fa 2. Tr o va i l p u n t o d i m i n i m o xdi una funzione partendo da un punto x(0) œRn. Vale infatti il seguente: Teorema 9. Si consideri (y)= 12yTAy≠ yTbcon :Rnæ R dove A eb sono la matrice dei coe cienti e il temine noto del sistema. Al lora xè soluzione di Ax= bse e solo se xè punto di minimo assoluto di . xsoluzione di Ax= b ≈∆ Ò (x)= 0 Osservazione 6.1. La funzione (y)= 12yTAy≠ yTbcon :Rnæ R nel caso n=1 rappresenta una parabola, e nel caso n=2 un paraboloide. Infatti: n=1 ∆ (y)= a 2y2≠by n=2 ∆ (y)= 1 2 #y1 y2$5 a11 a12 a12 a22 65 y1 y2 6 ≠#y1 y2$5 b1 b2 6 = 1 2 !a11y21+2 a12y1y2+a22y22 "≠y1b1≠y2b2 Dimostrazione. Si dimostrano le due implicazioni: 1. ∆ Si suppone che xsia soluzione di Ax= be si dimostra che sia anche il minimo della funzione : (x+v)= 1 2(x+v)TA(x+v)≠(x+v)Tb= = (x) ˙ ˝¸ ˚ 1 2xTAx≠xTb+1 2vTAx+ 1 2xTAv+ 1 2vTAv≠vTb= = (x)+ 1 2vTb+ 1 2xTAv+ 1 2vTAv≠vTb= Ricordando che: !xTAv"= !xTAv"T= vTAx= vTb 19 Poichè A è una matrice simmetrica (Def: 5 , Prop: 3.4): = (x)+ 1 2vTAvØ (x)+ 1 2⁄min ÎvÎ2> (x) ’vœRn,v”= 0 L’ultima disuguaglianza è valida poichè A è definita positiva (Prop: 4.1). Quindi ’v”=0 vale che (x+v)> (x),quindi xè il punto di minimo della funzione . 2. ≈ Si suppone che xsia punto di minimo della funzione e si dimostra che è soluzione di Ax= b: Se xè punto di minimo della funzione , allora per il teorema di Fermat Ò(x)= 0 (y1,...,y n)= 1 2 nÿ i,j=1 aijyiyj≠ nÿ i=1 biyi ˆ ˆyk= Se i=k ˙ ˝¸ ˚ 1 2 nÿ j=1 akjyj+ Se j=k ˙ ˝¸ ˚ 1 2 nÿ i=1 akiyi≠bk= 1 2 Sono uguali perchè Aèsimmetrica ˙ ˝¸ ˚ Q a1 2 nÿ j=1 akjyj+ 1 2 nÿ i=1 akiyi R b ≠bk= = nÿ i=1 akiyi≠bk= !Ay≠b" k Ovvero il k-esimo Ay≠b ∆Ò (y)= Ay≠b= ≠r ∆Ò (x)= ≠r= 0= ≠!Ax≠b" Si osservi che in realtà per questa seconda parte di dimostrazione non abbiamo utilizzato che A èdefinita positiva. Quindi l’essere punto di minimo è condizione su ciente per essere soluzione per una qualsiasi matrice A simmetrica. Siccome per il teorema di Fermat deve essere uguale a zero, implica che anche il residuo rsia pari a zero e che quindi xsia soluzione di Ax= b Quindi si tratta di costruire un algoritmo per trovare il punto di minimo di (y).Utilizzeremo il caso n=2 per farci guidare dall’intuizione geometrica. 20 Figura 2: Illustrazione metodo del gradiente. Per determinare xsi procede come segue: 1. Si sceglie un punto x(0) œRnarbitrario. 2. Si considera l’insieme di livello L(x(0)), che in questo caso è un’ellisse. L’obiettivo è quello di avvicinarsi al punto di minimo della funzione ;da x(0) la direzione di massima decrescita è ≠Ò !x(0) ". 3. Si considera la retta x(–)= x(0) +–r(0) ,da x(0) nella direzione opposta al gradiente e valuto lungo questa linea. In generale si ottiene una funzione in una sola variabile, di cui è facile trovare il minimo. Graficamente, nel caso n=2 , equivarrebbe a ”tagliare” il paraboloide con un piano verticale ottenendo una parabola. 4. Si trova il punto di minimo della funzione F(–)= !x(k)+–r(k)"(nel caso n=2 di una parabola, come si vede in Figura 2), che costituisce il punto x(1) da cui si parte per l’iterazione successiva. Figura 3: Le prime iterate generate dal metodo del gradiente sulle curve di livello di F 21 • In generale supponendo di avere x(k)si cerca x(k+1) della forma: x(k+1) = x≠– Direzione˙ ˝¸ ˚ Ò 1 x(k)2 = x(k)+–r(k) con r(k)= b≠Ax(k) minimizzando quindi F(–)= !x(k)+–r(k)".Quindi eguagliando a zero la derivata prima di F otteniamo: FÕ(–)= Ò 1 x(k+1) 2 = Ò 1 x(k)+–r(k)2 r(k)T= per definizione di derivata di una funzione composta in Rn. Ricordando la 2 si ha : = r(k)T1 A 1 x(k)+–r(k)2 ≠b 2 = r(k)T1 ≠r(k)+–Ar(k)2 = = ≠r(k)Tr(k)+–r(k)TAr(k)=0 FÕ(–)=0 ≈∆ –= –(k)= r(k)Tr(k) r(k)TAr(k) –(k)= r(k)Tr(k) r(k)TAr(k) (3) Algoritmo 3. 1. Si sceglie un x(0) arbitrario. 2. Si calcola r(k)(direzione di discesa): r(k)= b≠Ax(k) 3. Si calcola –(k)(lunghezza del passo): –(k)= r(k)Tr(k) r(k)TAr(k) 4. Si calcola x(k+1) : x(k+1) = x(k)+–(k)r(k) Osservazione 6.2. Il denominatore della 3 è diverso da zero perchè A è definita positiva (Def: 6) e simmetrica (Def: 5), e quando r(k)= 0significa che si è trovata la soluzione al passaggio precedente. Osservazione 6.3. Il metodo del gradiente fa parte dei metodi di Richardson dinamici, ovvero un metodo per cui l’elemento x(k+1) è ottenuto a partire da quello precedente x(k), a cui viene sommata un altro elemento (in questo caso pari a r(k)) moltiplicato per un coe ciente –, detto parametro di rilassamento (o di accelerazione): • Se –è costante (indipendente dall’iterazione k) si parla di metodo di Richardson stazionario. • Se –= –(k)(variabile e dipendente dall’iterazione k) si parla di metodo di Richardson dinamico. 6.2 Metodo del gradiente coniugato Quando si può usare 3. Quando A è simmetrica e definita positiva. A cosa serve 2. Serve per trovare la soluzione xdel sistema lineare del tipo Ax= b Cosa fa 3. Si cerca di migliorare la velocità di convergenza (Th: 12) del metodo del gradiente, scegliendo una direzione diversa da quella identificata da ≠Ò (data una generica direzione di discesa p(k), troveremo il valore –(k)come quel valore di –che rende minimo !x(k)+–p(k)"). 22 Figura 4: Direzione di p(0) rispetto a r(0) È possibile trovare delle direzioni p(k), lungo cui muoversi, migliori rispetto a r(k), rispetto cioè alla direzione definita dal gradiente ( r(k)= ≠Ò !x(k)"); si tratta cioè di una direzione per cui il metodo converge più velocemente alla soluzione x. Come si può osservare dalla figura 4, p(0) punta direttamente verso la soluzione x(che in questo caso coincide con il centro dell’ellisse, si sta sempre facendo riferimento all’esempio del paraboloide), e descrive quindi un angolo rispetto a r(0) che per costruzione deve essere minore di fi2(deve puntare verso l’interno). Figura 5: Le direzioni di discesa del gradiente coniugato (in linea tratteggiata e de- notate con CG) e quelle del gradiente (in linea continua ed indicate con G). Si noti come il metodo CG in due iterazioni abbia già raggiunto la soluzione. Osservazione 6.4. Dato p(k),l’ –(k)che minimizza !x(k)+–(k)p(k)"lo si trova eseguendo gli stessi passaggi che sono stati e ettuati per il calcolo nel metodo del gradiente (Eq: 3), ottenendo: –(k)= p(k)Tr(k) p(k)TAp(k) (4) Osservazione 6.5. Con questa scelta di –(k)si ottiene che il resituo al passaggio k+1 è ortogonale a p(k): r(k+1) ‹ p(k) Il problema a questo punto è determinare la direzione p(k). Determinazione di p(k) Adi erenza del metodo del gradiente, ogni passo deve essere ottimale tenendo conto anche dei passaggi precedenti: ad ogni passaggio si sceglie una direzione che sia ottimale non solo rispetto all’ultimo, ma a tutti i passaggi precedenti. 1. Si scelgono arbitrariamente x(0) (punto di partenza) e p(0) (ad esempio uguale a r(0) , anche se non è necessario). 23 2. Si trova il punto x(1) = x(0) +–(0) p(0) con –(0) che minimizza !x(0) +–p(0) ": –(0) = p(0) Tr(0) p(0) TAp(0) 3. Si trova poi il punto x(2) = x(1) +–(1) p(1) nella forma x(2) = x(0) +–(0) p(0) +–(1) p(1) e si cercano –(0) e–(1) tali che minimizzino !x(0) +–(0) p(0) +–(1) p(1) "= F!–(0) ,–(1) ", cioè una funzione a due variabili. x(2) è il punto di minimo di (v)tra tutti i vdella forma v= x(0) + –(0) p(0) + –(1) p(1) =∆ Si sta minimizzando lungo un piano (descritto dai vettori –(0) p(0) e–(1) p(1) ). Osservazione 6.6. Anzichè minimizzare solo l’ultimo passaggio, si scrive il vettore successivo come x(0) più una combinazione lineare di tutte le direzioni scelte in precedenza, e si fa in modo che x(k+1) sia minimo rispetto a tutte le precedenti direzioni. Cioè si minimizza prima lungo un’unica direzione (quindi lungo una retta), poi lungo due direzioni (quindi rispetto ad un piano) e così via. Dal teorema di Fermat nel punto di minimo deve essere ÒF!–(0) ,–(1) "= 0 4. Si calcola quindi ÒF!–(0) ,–(1) ": ˆF ˆ– (0) = p(0) TÒ 1 x(0) +–(0) p(0) +–(1) p(1) 2 = p(0) TË A 1 x(0) +–(0) p(0) +–(1) p(1) 2 ≠b È = = p(0) T1 ≠r(0) +–(0) Ap(0) +–(1) Ap(1) 2 =0 = ≠p(0) Tr(0) +–(0) p(0) TAp(0) +–(1) p(0) TAp(1) =0 (Per i passaggi vedi quelli e ettuati per il calcolo della 3) ˆF ˆ– (1) = p(1) TÒ 1 x(0) +–(0) p(0) +–(1) p(1) 2 = p(1) TË A 1 x(0) +–(0) p(0) +–(1) p(1) 2 ≠b È = = p(1) T1 ≠r(0) +–(0) Ap(0) +–(1) Ap(1) 2 = = ≠p(1) Tr(0) +–(0) p(1) TAp(0) +–(1) p(1) TAp(1) =0 Osservazione 6.7. Osserviamo che se nella prima equazione si ottiene lo stesso –(0) calcolato al primo passo se –(1) p(0) TAp(1) =0 5. Scegliamo p(1) tale che p(0) TAp(1) =0 ,sidiceche p(1) èA-ortogonale ap(0) , in questo modo: –(0) = p(0) Tr(0) p(0) TAp(0) –(1) = p(1) Tr(1) p(1) TAp(1) Osservazione 6.8. Dalla scelta di p(1) A-ortogonale a p(0) , posso mettere r(1) nell’espressione di –(k)al posto di r(0) perchè p(1) Tr(0) = p(1) Tr(1) Osservazione 6.9. A corregge il vettore per far sì che punti verso il centro. 24 Quindi nel caso generale si cerca: x(k+1) = x(0) + Combinazione lineare ˙ ˝¸ ˚ kÿ j=0 –(j)p(j) t.c. x(k+1) = Punto di minimo di ˙ ˝¸ ˚ argmin (v) Tr a t u t t i i v= x(0) + kÿ j=0 —(j)p(j) Si tratta cioè del minimo sullo spazio di kdimensioni. Con conti analoghi a quelli precedenti si trova che la scelta ottimale di p(k)è che sia A-ortogonale ai p(j)con j=1 ,...,k ≠1, tale che sia cioè ortogonali a tutti. p(k)TAp(k)=0 p(k)TA Ë r(k+1) ≠—(k)p(k)È = p(k)TAr(k+1) ≠—(k)p(k)TAp(k)=0 —(k)= p(k)TAr(k+1) p(k)TAp(k) –(k)= p(k)Tr(k) p(k)TAp(k) Teorema 10. Con questa scelta p(k+1) è A-ortogonale non solo a p(k)ma a tutti i p(j)(con j=0 ,...,k )esiccome la matrice A è simmetrica (Def: 5) e definita positiva (Def: 6), i vettori sono lineramente indipendenti. p(k)TAp(k)=0 Dimostrazione. L’unica combinazione lineare che possa dare il vettore nullo è quella con tutti i coe cienti nulli, prendendo una combinazione lineare: v= 0= kÿ j=0 v(j)p(j) p(i)TAv= p(i)TA Q a kÿ j=0 v(j)p(j) R b = kÿ j=0 v(j) ú ˙ ˝¸ ˚ p(i)TAp(j)= v(i)p(i)Ap(i)”=0 úæ È uguale a zero tranne per i= j (Per A-ortogonalità) Siccome A è definita positiva (Def: 6): p(i)TAp(i)Ø ⁄min Îp(i)Î2> 0 questo per il teorema th: 4.1, ∆ v(j)=0 ’i Tale espressione di —(k)deriva dal fatto che p(k+1) debba essere A-ortogonale a p(k),cioè: p(k)TAp(k)=0 p(k)TA Ë r(k+1) ≠—(k)p(k)È = p(k)TAr(k+1) ≠—(k)p(k)TAp(k)=0 Quindi vœRnex(n)èilminimoditutto Rn∆ È IL minimo. 25 Corollario 1. Il metodo è esatto dopo al più npa s s i x(k+1) = argmin (v) con v= x(0) + kÿ j=0 —(j)p(j) —(k)= p(k)TAr(k+1) p(k)TAp(k) Dove —(k)è l’entità del la correzione (del lo spostamento) tra l’elemento p(k+1) e quel lo p(k): p(k+1) = r(k+1) ≠—(k)p(k) Osservazione 6.10. Questa conclusione è vera in aritmetica esatta, in Matlab ci sono gli errori di calcolo ed è quindi necessario un criterio di arresto. Algoritmo 4. L’algoritmo si articola nei seguenti passi: 1. Si scelgono arbitrariamente x(0) (punto di partenza) e p(0) 2. Si calcola –(k): –(k)= p(k)Tr(k) p(k)TAp(k) 3. Si calcola x(k+1) : x(k+1) = x(k)+–(k)p(k) 4. Si calcola r(k+1) : r(k+1) = b≠Ax(k+1) = r(k)≠–(k)p(k) 5. Si calcola p(k+1) : p(k+1) = r(k+1) ≠—(k)p(k) 6. Si calcola —(k)in modo tale che p(k+1) sia A-ortogonale a p(k): —(k)= p(k)TAr(k+1) p(k)TAp(k) 6.3 Velocità di convergenza Definizione 13. Si definisce norma A (simmetrica e definita positiva) di xla quantità: ÎxÎ2A= xTAx Osservazione 6.11. Le proprietà di questa norma sono: • Se A è la matrice identità si ottiene la norma del vettore x • È sempre maggiore o uguale di zero. • È uguale a zero solo se x= 0 • Vale omogeneità: ΖvÎ= |–|·ÎvÎ • Disuguaglianza triangolare: Îx+yÎÆÎ xÎ+ÎyÎ 26 Teorema 11 (Velocità di convergenza metodo gradiente) .’x(0) œRn: Îx≠x(k)ÎAÆ K (A)≠1 K (A)+1 Îx≠x(k≠1)ÎA Îx≠x(k)ÎAÆ 3K (A)≠1 K (A)+1 4k Îx≠x(0) ÎA Osservazione 6.12. Se K (A)∫ 1 ∆ K(A)≠1 K(A)+1 ƒ 1e quindi si ha una convergenza molto lenta. Teorema 12 (Velocità di convergenza metodo gradiente coniugato) .’x(0) œRn p(0) œRn: Îx≠x(k)ÎÆ 2 AK (A)≠1 K (A)+1 Bk Îx≠x(0) ÎA Osservazione 6.13. Per la presenza del termine K (A)la situazione migliora rispetto al caso del metodo del gradiente, però se K(A)∫ 1si ha che 3ÔK(A)≠1 ÔK(A)+1 4 ƒ 1e la convergenza è lenta. 6.4 Criterio di arresto I criteri di arresto più comuni sono due: Criterio di arresto sul residuo normalizzato Sul residuo normalizzato si sceglie una tolleranza Á> 0e ci si arresta al primo k= kmin , ovvero il primo valore tale per cui: Îr(kmin )Î ÎbÎ Æ Á 1 r(k)= b≠Ax(k)2 Ponendo ˆx= x(kmin )e facendo gli stessi conti e ettuati per il caso dei metodi diretti (Eq: 2) si ottiene: Îx≠ ˆxÎ ÎxÎ Æ K (A)Îr(kmin )Î ÎbÎ Æ ÁK (A) Osservazione 6.14. Siccome di solito il numero di condizionamento K (A)è un numero molto grande, per avere una buona precisione è necessario scegliere un numero Ámolto piccolo. Di conseguenza il numero di condizionamento influenza anche il numero di iterazioni (quello minimo per avere una soluzione con una precisione accettabile). Criterio di arresto sull’incremento Si basa sulla di erenza tra una iterazione e la successiva: ”(k)= x(k+1) ≠x(k) Quindi il criterio d’arresto, preso un Á> 0, corrisponde al primo k= kmin tale per cui: Δ(kmin )ÎÆ ÁÎbÎ 6.5 Precondizionamento Idea: anzichè risolvere il sistema Ax= bpossiamo risolvere: P≠1Ax= P≠1b e scegliere una matrice P abbastanza vicina ad A tale che: P ¥ A t.c.=∆ P≠1A ¥ I (Matrice identità )∆ K !P≠1A"¥ 1 27 in modo tale da abbattere il numero di condizionamento, per velocizzare il metodo. Idealmente, se si prendesse P = A si avrebbe K !P≠1A"= K !A≠1A"= K (I)=1 . Tuttavia calolare A≠1 equivale a risolvere il sistema Ax= b, di conseguenza è necessario trovare una matrice P che sia facile da invertire da un lato, e che tale da rendere K !P≠1A"vicino a 1. Procederemo in realtà in questo modo: supponendo che P sia una matrice simmetrica e definita positiva, allora: ÷P12t.c. P12P12= P Infatti una matrice simmetrica e definita positiva è diagonalizzabile con tutti gli autovalori positivi e con matrice V degli autovettori ortogonale, e può quindi essere scritta come: P = VDV T VTV = Iæ Perchè V è una matrice ortogonale D æ Matrice diagonale con gli autovalori sulla diagonale principale Di coseguenza risulta essere: P12= VD 12VT dove D12è la matrice diagonale che sulla diagonale principale ha le radici quadrate degli autovalori Ô⁄i(ciò è possibile perchè A è definita positiva). Infatti: P12P12= 1 VD 12VT21 VD 12VT2 ú= VDV T= P Quindi si può scrivere: Ax= b≈∆ P≠12A x ˙˝¸ ˚ P≠12y= P≠12b dove: P≠12= 1 P122≠1= 1 VD 12VT2≠1ú= 1 VD ≠12VT2 (5) l’uguaglianza * è possibile grazie alle proprietà della matrice trasposta (Prop: 3.4). D≠12è una matrice diagonale che sulla diagonale principale ha i reciproci delle radici degli autovalori; si trova ye poi si calcola xcome x= P≠12y P≠12AP ≠12æ Simmetrica e definita positiva • Per verificare la simmetria: 1 P≠12AP ≠122T= 1 P≠122TA 1 P≠122T essendo A simmetrica per ipotesi (ricordando che per applicare il metodo del gradiente A deve essere simme- trica definita positiva), quindi AT= A Inoltre (considerando l’Eq: 5): 1 P≠122T= 1 VD ≠12VT2T= VD ≠12VT= P≠12 ∆ 1 P≠12AP ≠122T= P≠12AP ≠12= AP • Per verificare che è definita positiva: yTAPy= xT ˙ ˝¸ ˚ yTP≠12A x ˙˝¸ ˚ P≠12y= xTAx> 0 tranne per x= 0… y= 0 3 Essendo A definita positiva 4 28 Osservazione 6.15. Si possono quindi applicare il metodo del gradiente o il metodo del gradiente coniugato a P≠12AP ≠12y= P≠12b, con x= P≠12y Esiste una vasta letteratura sulle possibili P facili da invertire e tali che K 1 P≠12A 2 ¥ 1, un esempio è scegliere P come matrice diagonale avente gli elementi diagonali di A sulla diagonale: {aii}ni=1 Quando si può usare 4. Si usa nei metodi iterativi (ad esempio nel metodo del gradiente). A cosa serve 3. Serve a ridurre il numero di condizionamento della matrice. Algoritmo 5. Precondizionatori diagonali: in generale la scelta di P come la diagonale di A si può rivelare e ca ce nel caso in cui A sia una matrice simmetrica definita positiva. In generale la scelta di P come la diagonale di A si può rivelare e cace nel caso in cui A sia una matrice simmetrica definita positiva. Matlab 3 (Metodo del gradiente precondizionato) . X= p c g ( A , b ) ; 29 Metodi Analitici e Numerici per l’ingegneria Cerutti Maria Cristina A cura di: Andrea Fuso, Gianpiero Gaeta Indice I Richiami di Analisi 3 1 Funzioni 3 2 Vettori 4 3 Matrici 4 3.1 Algebra delle matrici . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 4 Teorema spettrale 6 II Calcolo numerico 7 5 Risoluzione di sistemi lineari con metodi diretti 7 5.1 Metodo di Cramer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 5.2 Risoluzione sistemi triangolari . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 5.2.1 Backward substitution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 5.2.2 Forward substitution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 5.3 Metodo di eliminazione di Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 5.4 Fattorizzazione LU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 5.4.1 Pivoting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 5.5 Fattorizzazione di Cholesky . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 5.6 Errore e condizionamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 6 Risoluzione di sistemi lineari con metodi iterativi 20 6.1 Metodo del gradiente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 6.2 Metodo del gradiente coniugato . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 6.3 Velocità di convergenza . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 6.4 Criterio di arresto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 6.5 Precondizionamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 7 Approssimazione di autovalori e autovettori 31 7.1 Metodo delle potenze . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 7.2 Metodo delle potenze generalizzato . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 8 Equazioni e sistemi non lineari 34 8.1 Metodo di bisezione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 8.2 Metodo di Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 8.2.1 Problemi di punto fisso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 8.2.2 Formula di Newton modificata . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 8.3 Criteri di arresto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 8.4 Soluzione di sistemi di equazioni non lineari . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 1 7 Approssimazione di autovalori e autovettori 7.1 Metodo delle potenze Quando si può usare 5. Si usa quando A è diagonalizzabile, che significa che la molteplicità geometrica di tutti gli autovalori è uguale a quella algebrica, o alternativamente quando esiste una base di autovettori. A cosa serve 4. Il metodo delle potenze serve per il calcolo dell’autovalore di modulo massimo, ⁄1, sotto opportune ipotesi. Funzionamento del metodo Euristicamente procediamo in questo modo: • Dato un vettore iniziale x(0) , il metodo consente di calcolare la successione di vettori, il cui k≠esimo è: x(k)= A(k)x(0) • Supponendo che A sia una matrice diagonalizzabile ( prima ipotesi), esiste una base di autovettori, di conse- guenza il vettore iniziale può essere scritto come una loro combinazione lineare: x(0) = nÿ i=1 –(i)v(i) • Ricordando poi che Av(i)= ⁄iv(i)(con i=1 ,...,n ) e sfruttando la precedente relazione si ottiene: x(k)= A(k)x(0) = nÿ i=1 –(i)A(k)v(i)= nÿ i=1 –(i)⁄kv(i) • Come seconda ipotesi si assume –(1) ”=0 e si ha: x(k)= –(k)⁄k1v(1) + nÿ i=2 –(i)⁄kiv(i)= –(1) ⁄k1 A v(1) + nÿ i=2 3–(i) –(1) 43 ⁄i ⁄1 4k v(i) B Osservazione 7.1. Il vettore iniziale x(0) ha una componente non nulla lungo l’autovettore v(1) , corrispondente al’autovalore ⁄1. • Supponendo che esiste un autovalore ⁄1strettamente maggiore degli altri: |⁄1|> |⁄2|Ø ··· Ø ··· > |⁄n|(terza ipotesi), i rapporti ---⁄i⁄1 ---< 1equindi 1⁄i⁄1 2kconvergono a zero per kche tende a +Œ , quindi si potrebbe dire che il vettore x(k)tende ad essere parallelo all’autovettore v(1) . Tale procedimento in realtà non converge ad un vettore in quanto: Underflow æ Quando |⁄1|< 1,⁄k1converge a zero Overflow æ Quando |⁄1|> 1,⁄k1diverge a +Œ Per evitare questo tipo di problemi si normalizza il vettore x(k)ad ogni iterazione, per evitare che la sua norma cresca o decresca eccessivamente: y(k)= x(k) Îx(k)Î x(k)= Ay(k≠1) 31 ∆ —(k)= 1 nr i=0 Îx(k)Î ∆ x(k)= —(k)⁄k1 A –(1) v(1) + nÿ i=2 –(i)3⁄i ⁄1 4k v(i) B In questo modo —(k)compensa ⁄k1: Non diverge se |⁄1|∫ 1 Non converge a zero se |⁄1|< 1 Infatti: —(k)= 1 ÎA(k)x(0) Î= 1 Î⁄k1 !–(1) v(1) +q(k)"Î q(k)= nÿ i=2 –(i)3⁄i ⁄1 4k v(i)≠æ 0 (Per kæ +Œ) ∆ —(k)|⁄1|k≠æ 1 |–(1) |·Îv(1) Î (Per kæ +Œ) ∆ x(k)≠æ 1 |–(1) |·Îv(1) Ζ(1) v(1) (Per kæ +Œ) Teorema 13. Sia A œCn◊nuna matrice diagonalizzabile i cui autovalori soddisfano la seguente relazione: |⁄1|> |⁄2|Ø ··· Ø ··· > |⁄n|,con ⁄1autovalore di modulo massimo. Assumendo –(1)”=0 , esiste C> 0tale che: Θx≠v(1)Î= C ----⁄2 ⁄1 ---- k +o A----⁄2 ⁄1 ---- kB kØ 1 dove: ˜x= x(k)ÎA(k)x(0)Î –(1) = v(1)+ nÿ i=1 –(i) –(1) 3⁄i ⁄1 4k v(i) k=1 ,2,... Osservazione 7.2. Questo metodo funziona se –(1) ”=0 , cioè se si è scelto un vettore x(0) con componente lungo v(1) diversa da zero. In questo caso però, nell’implementazione in Matlab, non bisoggna preoccuparsi di tale problematica, poichè in questo caso gli errori di approssimazione giocano a proprio favore: se si avesse scelto x(0) con –(0) =0 le iterazioni successive avrebbero fatto comparire uan componente parallela a v(1) . Osservazione 7.3. La velocità di convergenza di tale metodo è proporzionale a: Velocità di convergenza ƒ ---- ⁄2 ⁄1 ---- k cioè a quel termine che converge a zero più lentamente di tutti gli altri. 32 7.2 Metodo delle potenze generalizzato Considerando l’inversa di A:A≠1.Se A è diagonalizzabile, lo è anche l’inversa, quindi se |⁄n≠1|> |⁄n|(autovalore di A ), allora è possibile applicare il metodo delle potenze a A≠1per trovare ⁄n, reciproco dell’autovalore di modulo massimo per A≠1. Algoritmo 6 (Metodo delle potenze inverso) .L’algoritmo si articola nei seguenti passi: 1. Si sceglie un x(0) arbitrario (punto di partenza). 2. Si calcola y(0) come: y(0) = x(0) Îx(0) Î 3. L’elemento x(k+1) è poi calcolato come: x(k+1) = A≠1y(k) Osservazione 7.4. y(k)è conosciuto al passaggio precedente. 4. Una volta calcolato x(k+1) è possibile calcolare y(k+1) : y(k+1) = x(k+1) Îx(k+1) Î Osservazione 7.5. Ad ogni passo si ha un sistema lineare di cui cambia il termine noto (in questo caso y(k)), di conseguenza viene molto comodo usare il metodo LU (Sezione: 5.4). È possibile un’ulteriore generalizzazione: se esiste un autovalore di A vicino ad un certo µœC allora è possibile applicare il metodo delle potenze inverse per determinare l’autovalore ⁄, che ha distanza minima da µ, applicando alla matrice A≠µI Criterio di arresto Îy(k+1) ≠yÎÆ Á 33 8 Equazioni e sistemi non lineari 8.1 Metodo di bisezione Teorema 14 (Teorema degli zeri) .Sia f:[a, b ]æ R una funzione continua e sia il prodotto f(a)f(b)< 0al lora: ÷– œ (a, b )t.c. f(–)=0 Dimostrazione. La dimostrazione è e ettuata mediant eil metodo di bisezione: • Si suppone (ad esempio) che f(a)< 0eche f(b)> 0 • Si prende in considerazione punto medio cpreso a metà dell’intervallo (a, b ): a(0) = a ; b(0) = b ; c(0) = a(0) +b(0) 2 • A questo punto si calcola f!c(0) "e si hanno tre possibilità: 1. f!c(0) "=0 ∆ il metodo si arresta perchè è stato trovato quel valore di –tale che f(–)=0 2. f!c(0) "< 0 ∆ si considera il nuovo intervallo #a(1) ,b (1) $= #c(0) ,b (0) $ 3. f!c(0) "> 0 ∆ si considera il nuovo intervallo #a(1) ,b (1) $= #a(0) ,c (0) $ • Il nuovo intervallo #a(1) ,b (1) $èinclusonelprecedente( #a(1) ,b (1) $µ #a(0) ,b (0) $) e di ampiezza: b(1) ≠a(1) = b(0) ≠a(0) 2 = b≠a 2 • Quindi generalizzando al k-esimo passaggio si ha: Ë a(k),b (k)È ; c(k)= a(k)+b(k) 2 = b≠a 2k • Si valuta il generico c(k)nella funzione f: 1. Se f!c(k)"=0 ∆ il metodo si arresta perchè è stato trovato quel valore di –tale che f(–)=0 2. Se f!c(k)"< 0 ∆ si considera il nuovo intervallo #a(k+1) ,b (k+1) $= #c(k),b (k)$ 3. Se f!c(k)"> 0 ∆ si considera il nuovo intervallo #a(k+1) ,b (k+1) $= #a(k),c (k)$ • Si ottiene quindi l’intervallo #a(k+1) ,b (k+1) $tale che f!a(k+1) "< 0ef!b(k+1) "> 0e ampiezza pari a : b(k+1) ≠a(k+1) = b(k)≠a(k) 2 = b(0) ≠a(0) 2k+1 = b≠a 2k+1 • Si costruiscono ora le successioni )a(k)* kœNe)b(k)* kœNtali che: f 1 a(k)2 < 0; f 1 b(k)2 > 0 b(k)≠a(k)= b≠a 2k ≠æ 0( Per kæ +Œ) )a(k)* è non decrescente )b(k)* è non crescente 6 æ Sono limitate ∆ Teorema esistenza del limite per le successioni monotone Quindi a(k)¬ a(a(k)tende per difetto a a)e b(k)√ b(b(k)tende per eccesso a b), ma b≠a=lim kæ+Œb(k)≠a(k)=lim kæ+Œ b≠a 2k =0 ∆ b= a= – (Candidata radice dell’equazione ) 34 f(–)= lim kæ+Œf!a(k)" Æ 0 øø Per Teorema continuità permanenza segno ¿¿ f(–)= lim kæ+Œf!b(k)" Ø 0 T XXXXXXXV ∆ f(–)=0 Implementazione metodo di bisezione Algoritmo 7. L’algoritmo per il metodo di bisezione si articola nei seguenti passi: 1. Scelta degli estremi del l’interval lo a(0) = aeb(0) = bin modo tale che f(a)f(b)< 0: ∆ x(0) = b≠a 2 dove x(0) è l’ampiezza del l’interval lo iniziale. 2. Per i passaggi kØ 1: (a) Se f!x(k≠1)"=0 ∆ –= x(k≠1) (b) Se f!x(k≠1)"< 0 ∆ a(k)= x(k≠1) eb(k)= b(k≠1) (c) Se f!x(k≠1)"> 0 ∆ a(k)= a(k≠1) eb(k≠1) = x(k≠1) 3. Il nuovo valore: x(k)= a(k≠1)+b(k≠1) 2 4. Si calcola l’errore: e(k)= ---x(k)≠– ---< 1 2 ---I(k)---= 31 2 4k+1 < Á dove --I(k)--è la misura del l’interval lo I(k)= #a(k),b (k)$ 5. Si determina il criterio di arresto, ci si arresta a kmin ,primo kper cui: 31 2 4k+1 (b≠a)< Á 31 2 4k+1 < Á b≠a ≈∆ 2k+1 > b≠a Á kmin > log 2 3b≠a Á 4 ≠1 Osservazione 8.1. È un criterio molto semplice ma converge mlto lentamente e non tiene conto di come è fatta la funzione. 35 8.2 Metodo di Newton Sia f :I=[ a, b ]æ R tale da essere fœC1(I)(1), con fÕ(–)”=0 , allora il metodo di Newton permette di trovare quel valore –tale per cui f(–)=0 attraverso i seguenti passaggi: 1. Si parte da un x(0) arbitrario e si scrive l’equazione della retta tangente al grafico di fnel punto !x(0) ,f !x(0) "": y= f 1 x(0) 2 +fÕ1 x(0) 2 · 1 x≠x(0) 2 Osservazione 8.2. Corrisponde al polinomio di Taylor della funzione f del primo ordine, che determina appunto la retta (quella tangente) che meglio approssima la funzione fnell’intorno di x(0) 2. Si trova il punto di intersezione x(1) della retta ycon l’asse delle ascisse, si trova cioè la radice della retta tangente: f 1 x(0) 2 +fÕ1 x(0) 2 · 1 x(1) ≠x(0) 2 =0 ∆ x(1) = x(0) ≠ f!x(0) " fÕ!x(0) " 3. In generale alla k-esima iterazione la retta tangente al punto x(k)risulta essere: y= f 1 x(k)2 +fÕ1 x(k)2 · 1 x≠x(k)2 cioè quella retta tangente al punto !x(k),f !x(k)"" . 4. Si trova il punto x(k+1) come l’intersezione della retta yprecedentemente calcolata e l’asse delle x: x(k+1) = x(k)≠ f!x(k