Teoria calcolo

Modellistica matematica e problemi numerici

Partiamo sempre da un problema fisico. Attraverso la modellistica matematica arrivo al problema matematico, che sarebbe solo la traduzione del problema fisico tramite equazioni, (d sono i dati). Passo attraverso l'approssimazione numerica e arrivo al problema numerico, che è un'approssimazione del problema matematico. Poi all'algoritmo, che è una sequenza finita di operazioni. Attraverso il linguaggio MATLAB e Octave ho una soluzione computazionale. Il calcolatore lavora con un'aritmetica finita, che ha una memoria finita. Aritmetica floating point, per certi numeri ci sarà un'approssimazione.

Devo verificare che la soluzione numerica Xn esiste ed è unica. Verificare che sia stabile, cioè la soluzione numerica deve dipendere dai dati in una maniera limitata. Verificare che sia convergente, ovvero che |Xn-Xm|, che sarebbe l'accuratezza, voglio la differenza sia piccola.

Problemi matematici

• Radici di equazioni non lineari
• Sistemi lineari
• Interpolazione → fornire una funzione che passa in più punti
• Integrazione
• Equazioni differenziali ordinarie (EDO) → che sono i problemi di Cauchy

Problema di variazione batterica

Se la derivata è uguale a 0, la variazione nel tempo del numero di batteri è nulla. Quindi io cresco ma arrivo a una certa che non cresco più. Questo era il problema matematico. Per arrivare al problema numerico uso le approssimazioni. L'idea è di discretizzare l'intervallo di tempo T in intervalli di ampiezza h. Tra un intervallo e l'altro dico che c'è un tempo tk. Devo soddisfare l'equazione iniziale per ogni tempo, provo così a soddisfarla per i tempi discreti che io ho generato.

Ricerca di radici di funzioni non lineari

Radici ed equazioni non lineari

Immagino di investire una quota v=quota annuale, per un numero di anni n. R= tasso di rendimento. La quota che io ottengo dopo tot anni → Quanto deve valere r per accumulare M dopo n anni, dato l'investimento annuale v? Obiettivo: data f(x), f : I → R Voglio trovare numero reale α ∈ R tale che f(α)=0. α è la radice.

Esempio: Ho una funzione f(x) Voglio trovare α. Ma questa era un’equazione lineare. Spesso non è possibile determinare α in maniera analitica, quindi uso metodi numerici.

Metodi iterativi locali

L'idea è di costruire una successione di approssimazioni di α, nello specifico avrò:

Esempio: Voglio costruire una successione di approssimazioni che tendono ad andare verso la mia radice.

Convergenza

NB non posso fare infinite iterazioni.

Metodo di Newton

• Suppongo di essere in un punto E di aver già trovato un'approssimazione al passo k-esimo. Valuti f al passo k-esimo e costruisci la retta tangente. L'iterata di Newton è tale che: è l'intersezione tra Yk(x) con l'asse delle ascisse. L'iterata di Newton è:

Metodo di Newton (algoritmo)

• Scelgo iterata iniziale
• Per k=0,1,2...

Osservazione: sulla carta il metodo converge in un infinito numero di iterazioni. L'algoritmo viene in realtà arrestato quando un certo criterio di arresto è soddisfatto. Mi voglio fermare quando:

Esempio: I metodi iterativi locali, come Newton, convergono se l'iterata iniziale è sufficientemente vicino ad alfa. Sufficientemente vicino vuol dire che:

Esempio: Mi sono allontanato dalla mia radice. Questa è una divergenza. Questo succede perché il mio metodo genera una successione che non converge, perché l'iterata iniziale non era sufficientemente vicino ad alfa.

Osservazione: questo comportamento e la convergenza di Newton si può spiegare interpretando il metodo come metodo delle iterazioni di punto fisso.

Osservazione: per usare Newton devo calcolare la derivata prima della funzione, questa può essere difficoltosa. Si possono usare metodi alternativi a Newton.

Metodo delle corde

• Q approssima la derivata prima. Usare questo metodo vuol dire costruire una pendenza che è sempre fissa per ogni iterata, è data dalla retta passante per il punto iniziale e finale dell'intervallo.

Metodo delle secanti

• Dice che: per k=1,2,...

Osservazione: funzionano se la prima iterata è sufficientemente vicino a alfa 24/2/21.

Iterazioni di punto fisso

Osservazione:

Definizione: detta phi(x) la funzione di iterazione, alfa è punto fisso se alfa=phi(alfa).

Osservazione: se phi(x)= f(x) + x allora alfa è radice di f <-> alfa è punto fisso di phi.

Interpretazione grafica del punto fisso

Traccio bisettrice I-III quadrante, e identifico l'intersezione.

Osservazione: voglio approssimare il punto fisso alfa di phi. C'è un algoritmo → delle iterazioni di algoritmo punto fisso.

Osservazione:

Esempio: La derivata prima è minore di 1 ma maggiore di -1, se è negativa converge a ragnatela, se no convergo sempre dallo stesso lato.

Proposizione: convergenza globale in un intervallo

• Se phi è continua in [a,b] e phi(x) appartiene ad [a,b] per ogni x appartenente ad [a,b], allora esiste almeno un punto fisso di phi, alfa appartenente a [a,b]
• Se oltre all'ipotesi 1, esiste una costante L tale che 0<=L<1, e per cui per ogni x1,x2 appartenenti a [a,b], allora il punto fisso è unico e l'algoritmo delle iterazioni di punto fisso converge ad alfa per ogni scelta di (posso partire da qualsiasi punto nell'intervallo per arrivare poi a convergere)

Per avere convergenza richiedo che tutti i rapporti incrementali in [a,b] siano minori di 1 in modulo.

Proposizione: convergenza globale in un intervallo

• Se phi è 1) continua in [a,b] e phi(x) appartiene ad [a,b] per ogni x appartenente ad [a,b], allora esiste almeno un punto fisso alfa in [a,b]
• Se phi è derivabile con continuità in [a,b] allora il punto fisso è unico in [a,b] e le iterazioni di punto fisso convergono ad alfa, ovvero:

Osservazione:

Osservazione: L è il valore massimo che assume la derivata prima di phi in [a,b], se phi è derivabile con continuità nell'intervallo.

Osservazione: se L è vicino a 0, cioè se la derivata prima di phi è circa zero, allora l'errore si abbatte in poche iterazioni. Se L è circa 1, ma minore di 1, servono più iterazioni per ridurre l'errore. (Se ho funzione di iterazione più piatta possibile allora mi servono poche iterazioni).

La convergenza globale spiega come convergono le iterazioni di punto fisso per ogni scelta di un’iterata iniziale in [a,b].

Teorema: convergenza locale/di Ostrowski

• (Ovvero in un intorno del punto fisso) Ad ogni iterata l'errore si riduce. Se phi’ è 0.99 l’errore si riduce di poco.

Osservazione:

Proposizione (di convergenza locale)

Allora il metodo delle iterazioni di punto fisso converge per sufficientemente vicino ad alfa con ordine p=2, (meglio convergere con ordine 2 perché l’errore sta scendendo in modo più rapido.). Convergenza rapida per funzioni che tagliano la bisettrice in modo tipo piatto, cioè la derivata prima è zero. E si dice convergere con ordine 2.

Osservazione: Criterio di arresto del metodo di iterazioni di punto fisso

Non posso iterare fino ad infinito criterio differenza tra iterate successive.

Introduco un indicatore dell’errore → → mi fermo a k+1 tale che:

Il criterio è soddisfacente se

Il criterio non è soddisfacente se. Perché l’errore è sottostimato dalla differenza di iterate successive.

• Se errore sovrastimato dalla differenza, criterio è comunque soddisfacente anche se faccio più del necessario.
• Metodo di Newton come metodo di iterazioni di punto fisso. Se scelgo, allora posso interpretare il metodo di Newton come metodo di iterazioni di punto fisso.

Definizione: la radice alfa di f è semplice se f(alfa)=0 ma f’(alfa) diverso da 0

Definizione: la radice di f ha molteplicità m se f(alfa)=0

Proposizione

Se alfa è radice semplice (quindi m =1) allora il metodo di Newton converge con ordine p=2, infatti ho che:

Se alfa è radice multipla (m>1) allora il metodo di Newton converge con ordine p=1

Metodo di Newton modificato

• L’iterata k-esima dipende da m. Nota m posso ripristinare l’ordine di convergenza p=2 anche se m>1.

Criteri di arresto per il metodo di Newton

• Differenza tra iterate successive: arresto le iterazioni quando la differenza tra due iterate successive scende sotto una certa tolleranza (scelta a priori). Questo è un buon criterio quando:

Siccome il metodo di Newton si può interpretare come metodo delle iterazioni di punto fisso, è noto che il criterio sia soddisfacente se

Se alfa è semplice allora

Se alfa è multiplo il criterio non è soddisfacente perché l’errore è sottostimato dalla differenza tra iterate successive

• Criterio del residuo: mi fermo quando il valore del residuo al passo k che è. Deve essere minore di una certa tolleranza.

Sappiamo che f(alfa)=0. Quando il criterio è soddisfacente? Cioè. Soddisfacente se f’(x)=1 per ogni x appartenente all’interno di alfa.

Se invece la pendenza è maggiore di 45 →. Il criterio è insoddisfacente perché l’errore vero è sovrastimato dal residuo. Vuol dire che mi fermo quando residuo è sotto una certa tolleranza ma il mio errore è molto più piccolo. Se invece la pendenza è molto bassa l'errore è molto maggiore del residuo. Il criterio è quindi insoddisfacente perché l’errore è sottostimato dal residuo.

Metodo delle corde

• Tale metodo può essere interpretato come metodo delle iterazioni di punto fisso con iterata: Il metodo converge se:

Osservazione: il metodo delle secanti non è interpretabile come metodo delle iterazioni di punto fisso.

Sistemi lineari

Risolvere un sistema lineare Ax=b, dove l’elemento della matrice A

Definizione: la matrice A è non singolare <-> det(A) diverso da 0

Proposizione: se A è non singolare, allora esiste un’unica soluzione x del sistema lineare.

Osservazione: in generale, dobbiamo risolvere sistemi lineari di dimensione n grande. Devo usare un calcolatore. Devo usare metodi per risolvere il sistema lineare che si traducano in algoritmi computazionalmente efficienti.

Esempio: Emessi dei raggi X e una parte viene assorbita dal corpo in maniera diversificata a seconda delle proprietà del materiale (la densità). L’idea è che il tubo che emette i raggi e il rilevatore devono girare. Il nostro corpo venga diviso in piccole dimensioni, se la densità di queste porzioni è grande allora arriveranno meno raggi X al rilevatore.

In bidimensionale: bi = porzione di raggi X assorbita lungo la direzione i-esima. Questo può essere visto come un sistema Ax=b (A)_ij = 1 se il pixel j-esimo è attraversato dal raggio i-esimo.

Osservazione: come risolvere sistema lineare?

Il numero di operazioni da effettuare con il metodo di Cramer è (3(n+1)!), un calcolatore ha una CPU che è in grado di effettuare 1 miliardo di operazioni al secondo (pari a 1GHz). Tempo CPU= n operazioni che farei a meno/potenza in GHz.

Osservazione: se il determinante di A è diverso da 0 esiste. Con il calcolatore non lo faccio mai perché costa molto di più e richiede molta memoria al computer per eseguire l’operazione. Gli errori computazionali inoltre si amplificano.

Metodi diretti

• Sono metodi numerici per risolvere il sistema lineare in cui si effettuano un numero finito di passaggi noto a priori.

Osservazione: i metodi diretti riconducono la soluzione di Ax=b a quella di sistemi lineari più “semplici” da risolvere. Vediamo i metodi semplici:

• Matrice triangolare inferiore

Algoritmo delle sostituzioni in avanti: ovvero parto dalla prima riga trovando X1, lo metto nella seconda riga e trovo X2 ecc. In particolare l’algoritmo è:

Quante operazioni: divisioni → (una sulla prima riga e una sulle righe successive) n divisioni. Prodotti → Somme e sottrazioni → La somma di queste operazioni fa n²

Osservazione: det(L) = prodotto degli elementi che stanno sulla diagonale principale

• Matrice triangolare superiore

Algoritmo delle sostituzioni all’indietro N² operazioni

Metodo della fattorizzazione LU

Def: data una matrice non singolare, la sua fattorizzazione LU consiste nel trovare una matrice triangolare inferiore L e una matrice triangolare superiore U, tale che A=LU (se esistono). Oss: voglio risolvere Ax=b. Se esistono L e U tale che A=LU allora LUx=b Introduco vettore ausiliario y=Ux Il primo sistema è triangolare inferiore, il secondo triangolare superiore. Li risolvo con gli algoritmi trovati sopra

Def: il metodo della fattorizzazione LU per risolvere il sistema lineare Ax=b consiste in:

• Determinare se esiste la fattorizzazione LU di A
• Risolvere sistema triangolare inferiore (Ly=b) con l’algoritmo delle sostituzioni in avanti
• Risolvere sistema triangolare superiore (Ux=y) con l’algoritmo delle sostituzioni all’indietro

Oss: n=2 Problema sottodeterminato perché ho più coefficienti che equazioni. In generale ho n+n² coefficienti e solo n² equazioni (vincoli) → per convenzione fisso n elementi di L pari a 1 sulla diagonale principale.

Metodo di eliminazione di Gauss (MEG)

È un metodo che fornisce i fattori L e U di una matrice A non singolare, ammesso che esistano. → algoritmo MEG

Esempio: Come faccio a fare venire lo 0 lì? Prendo la seconda riga di A la devo sottrarre a un coefficiente e moltiplico con la prima.

Oss: stavamo vedendo le soluzioni del sistema lineare con metodi diretti, in particolare con il metodo della fattorizzazione LU

• Determino L e U tale che A=LU tramite MEG
• Soluzione Ly=b con algoritmo sostituzioni in avanti
• Soluzione Ux=y con algoritmo sostituzione indietro

L’algoritmo trasforma A in U e crea la matrice L dei moltiplicatori, sfruttando gli elementi pivotali che compaiono

Esempio: Oss il MEG comporta un costo computazionale di (2n³/3) operazioni

• MEG 2n³/3
• Ly=b n²
• Ux=y n²

Per n grande prevale MEG che è comunque