Appunti di Teoria di Calcolo Numerico

Floating Point: l'insieme dei numeri macchina(numeri a virgola mobile) è un sottoinsieme

proprio dei numeri reali ( ). Questo insieme ha la particolarità di avere un numero finito

di elementi. Dato la sua rappresentazione floating point è con

Troncamento: approssimazione che comporta nella perdita di alcune cifre decimali dopo una

certa posizione. Ad es.:

Arrotondamento: approssimazione che comporta nella perdita di alcune cifre decimali dopo

una certa posizione, con la particolarità di lasciare nell'ultima cifra rimanente un'informazione

sulla prima cifra sottoposta ad eliminazione: si aggiunge +0 all'ultima cifra se la prima cifra

eliminata mentre si aggiunge +1 all'ultima cifra se la prima cifra eliminata

Ad es.:

Tipologie di formato di dato

o format short: 4 cifre dopo la virgola

o format short e: una cifra prima la virgola e 4 cifre dopo la virgola, in esponenziale

o format long: 15 cifre dopo la virgola

o format long e: una cifra prima la virgola e 15 cifre dopo la virgola, in esponenziale

Rappresentazione floating point:

con : è l'ordine di grandezza, è il numero di cifre,

è la base, indica il segno (0 positivo, 1 negativo)

Mantissa: . La condizione serve ad avere una rappresentazione

univoca di un numero floating point(garantisce l'unicità della rappresentazione).

Il numero 0 sarebbe un eccezione e infatti (0 viene trattato in maniera differente)

Insieme floating point:

è la base, è il numero di cifre, è il valore minimo di mentre è il

valore massimo di . Il codice binario dei PC utilizza

Errore assoluto: differenza tra due valori, spesso il primo è ottenuto come approssimazione del

secondo e la loro differenza si chiama errore assoluto. Ad es.: è l'errore assoluto

tra un numero reale e il suo corrispondente floating point.

Errore relativo: rapporto tra l'errore assoluto e la quantità esatta. Questo errore è molto

significativo e affidabile in quanto non tiene conto dell'ordine di grandezza. Ad es.:

è l'errore relativo di un numero floating point rispetto al suo corrispondente reale.

Epsilon macchina: è l'errore dovuto al fatto che l'insieme ha un numero finito di cifre(16 in

decimale). tale che

.

L'epsilon macchina è la distanza tra e il suo successivo; inoltre all'aumentare della

dimensione del numero l' diventa sempre più grande mentre all'avvicinarsi a 0 diventa

sempre più piccolo.

Distribuzione floating point: l'insieme è una retta discontinua con distribuzione non

uniforme, bensì più densa vicino allo 0 e sempre meno densa man mano che si aumenta con

l'ordine di grandezza dei numeri.

Numero macchina minimo:

Numero macchina massimo:

Overflow:

Underflow:

Cardinalità di : quantità di elementi(numeri) nell'insieme floating point(con lo 0 aggiunto)

o è dovuto al segno

o è dovuto a tutte le possibili cifre della prima cifra della mantissa(tolto lo 0)

o è dovuto a tutte le possibili cifre delle altre cifre della mantissa

o tiene conto di tutti i valori assunti dall'esponente

o serve per includere lo 0 nel conteggio

Proprietà di :

o Commutativa:

o Non unicità dello zero: se succede che

o Associativa: purché oppure non vada in

overflow. Ad es.: se , posto , ma allora si

ha un fenomeno di overflow. Infatti non darà overflow in quanto non

viene mai superato in risultati intermedi, cosa che accade quando facciamo

o Cancellazione numerica: fenomeno dovuto all'approssimazione causata dalla differenza

tra la base con cui viene immesso un dato(decimale) e la base utilizzata dalla

macchina(binaria). Infatti se si calcola

con il calcolatore, posto si ha

come risultato e non come risulterebbe realmente.

Teorema di Abel:

Non esiste un numero finito di passi per trovare le radici di un polinomio di grado uguale o

superiore al 5°. In pratica non esistono formule chiuse per trovare le radici del polinomio.

Schema iterativo: metodo per il calcolo diretto delle radici di un polinomio mediante il quale,

grazie ad un punto di partenza, detto guess iniziale, si ripete un numero finito di passi

generando una successione di soluzioni che deve convergere alla soluzione con una certa

velocità e arrestandosi una volta raggiunte alcune prestabilite condizioni.

Ordine di convergenza: velocità impiegata da uno schema iterativo per raggiungere la

soluzione con la precisione richiesta. Maggiore è la velocità, migliore è lo schema utilizzato (in

termini economici di tempo).

Criterio di arresto: condizione con la quale decido se continuare a cercare la soluzione oppure

se fermarmi in quanto l'accuratezza della soluzione appena trovata mi soddisfa.

Stimatore: elemento con il quale impongo la precisione desiderata per la soluzione. Viene

solitamente relazionato all'errore, il quale viene maggiorato dallo stimatore.

Metodo di Bisezione:

e almeno un tale che (teorema Zeri)

Questo metodo diretto serve a trovare una radice di un polinomio, purché la soluzione sia

compresa tra e . Per farlo si sceglie l'intervallo , lo si divide a metà e si cerca in quale

semi-intervallo si trova la radice(sfruttando il teorema degli zeri). Si ripete questa procedura

volte fino ad avere il semi-intervallo -esimo minore della tolleranza richiesta sulla soluzione.

Il criterio di arresto è con tolleranza. L'intervallo -esimo è

poiché è

quello iniziale, e ad ogni iterata viene dimezzato. Il numero di iterazioni minimo al fine di

raggiungere la tolleranza richiesta si trova nella seguente maniera:

Essendo un numero non intero, se si esegue l'algoritmo di bisezione un numero di volte pari

all'intero successivo a si è sicuri di avere la precisione richiesta. Se il metodo è stato

eseguito correttamente si avrà che .

Criterio del residuo: criterio d'arresto che si basa sul residuo. Il residuo è il valore della funzione

nella posizione ovvero . Quando il valore assoluto di questa quantità è

minore della tolleranza allora si può arrestare il metodo diretto per trovare . Tuttavia se la

derivata della funzione in è molto maggiore di 1 si può avere una sovrastima mentre se la

derivata in è molto minore si può avere una sottostima. Lo stimatore che sfrutta il residuo è

buono quando la derivata in è circa 1.

Criterio dell'incremento: criterio d'arresto che si basa sull'incremento. L'incremento è la

differenza tra due iterate successive. Se questa differenza è minore della tolleranza si può dire

di essere abbastanza vicini alla soluzione.

Metodo di Newton

Questo metodo diretto per trovare le radici di una funzione sfrutta delle rette tangenti alla

funzione passanti per le varie .

tale che

Newton approssima Taylor, infatti la formula dell'iterata successiva si ricava da Taylor

troncata al 2° ordine:

Se trascuriamo l'O-grande del secondo ordine si ha che per sufficientemente grande

abbiamo come iterata finale e quindi

dalla quale si ricava l'espressione di

Convergenza del Metodo di Newton:

Condizioni per la convergenza del metodo di Newton:

1. deve essere sufficientemente vicino ad (cosa che si può fare facendo qualche

iterata con bisezione oppure guardando il grafico della funzione.

2. deve essere una radice semplice, ovvero

Se sono soddisfatte queste condizioni Newton converge. Con che velocità?

3. Se allora si può dire che l'ordine di convergenza è dato dal limite

Se il limite è verificato Newton converge quadraticamente.

4. Se è multipla allora, se Newton converge, lo fa linearmente.

Ordine di convergenza :

Uno schema iterativo ha convergenza se indipendente da tale che

Se il metodo converge linearmente mentre se il metodo converge

quadraticamente(e cosi via per valori maggiori). Essendo

si ha che per

è necessario che per avere convergenza, altrimenti l'errore non diminuisce ma

aumenta ad ogni iterata.

Metodo di Newton modificato

Il metodo di Newton modificato serve nel caso si abbia una radice multipla (ad esempio

ha radice multipla in e ). Se si ha radice multipla e se

questa radice ha molteplicità allora lo schema diventa

il quale converge più velocemente del metodo standard di Newton.

(esiste anche un ulteriore metodo, chiamato Newton Adattativo, che cambia ad ogni iterata)

Punto fisso: è punto fisso se

Trovare il punto fisso di una funzione equivale a trovare l'intersezione della funzione con

la bisettrice del primo e del terzo quadrante. Si possono avere più punti fissi( ), un solo

punto fisso( ) o nessun punto fisso( )

Metodo di Punto Fisso

Il metodo di punto fisso serve per trovare un punto tale che sfruttando una

funzione tale per cui . Esistono più adatte allo scopo, ad esempio

oppure

.

L'iterazione di punto fisso si basa su

Il metodo di Newton è un metodo di punto fisso con

Tipologie di convergenze a punto fisso:

è detto attrattore se si ha convergenza, repulsore se si ha divergenza

1. Convergenza non monotona(es.: e decrescente)

2. Convergenza monotona(es.: e decrescente)

3. Divergenza non monotona(es.: e crescente)

4. Divergenza monotona(es.: e decrescente)

Per monotona si intende (oppure a seconda

del caso) mentre per non monotona si intende che le varie iterate sono maggiori o minori

ad a seconda dell'iterata( ma ecc…)

Teorema Globale per la convergenza di Punto fisso:

1. Condizione sufficiente per l'esistenza di

2. Condizione sufficiente per l'esistenza di e per la convergenza di punto fisso

(richiesta lipschitziana)

Dimostrazione teorema globale per la convergenza di punto fisso

1. Prendiamo

Ciò vuol dire che, per il teorema degli zeri, almeno un

ovvero .

2. Supponiamo per assurdo che

entrambe punti fissi

Allora

Si ha quindi che però è assurdo perché (lipschitzianità)

Si conclude dunque che

Per dimostrare la convergenza si considera l'errore

per punto fisso(metodo e definizione)

Per lipschitzianità da cui

Di conseguenza

Posto l'errore si ha quindi che

dove è conosciuto e limitato mentre

tende a 0 per in quanto e perciò punto fisso converge.

Teorema locale di Ostrowski

Dimostrazione del teorema di Ostrowski

è un punto compreso tra e (per il teorema del valor medio).

Per ipotesi di continuità allora vale anche che, sostituendo, diventa

Fattore asintotico di convergenza: indica la velocità con la quale si raggiunge . La

velocità di convergenza della successione di iterate aumenta man mano che

Con la successione diverge, con non posso dire nulla sulla convergenza

mentre con la successione converge, come visto con Ostrowski. Nel teorema

globale il fattore asintotico è .

Generalizzazione di Ostrowski per convergenza quadratica

Dimostrazione di Ostrowski per convergenza quadratica

è un punto compreso tra e .

Essendo per ipotesi si ha dunque

Per ipotesi di continuità allora vale anche che diventa

Generalizzazione di Ostrowski per ordine di convergenza > 2

Dimostrazione di Ostrowski per ordine di convergenza > 2

è un punto compreso tra e .

Essendo per ipotesi si ha dunque

Per ipotesi di continuità allora vale anche che diventa

Per trovare basta cercare la prima derivata che non si annulla in quanto l'ordine di quella

derivata è l'ordine di convergenza.

Dimostrazione del fattore di convergenza di Newton

Newton converge con quindi

.

Valutando in (e ricordando che ) si ottiene che

(in Newton per con molteplicità si ha e

)

Legame tra errore e criterio d'arresto dell'incremento:

è l'incremento.

per la definizione di punto fisso.

Se allora compreso tra e

Da qui si nota che il caso ideale è . In questo caso l'incremento

diventa un buon stimatore dell'errore in quanto si ha . Se perdo

affidabilità nel criterio basato sull'incremento. per radici semplici, il che rende

questo criterio d'arresto favorito in questi casi.

Andamento errore/stimatore per radici con molteplicità diversa da 1

In Newton standard se ho con molteplicità si ha e

Questo comporta ad avere all'aumentare di e quindi a perdere affidabilità nel

criterio d'arresto dell'incremento. L'andamento dell'errore e dello stimatore(incremento) è

sempre parallelo per ogni ma la velocità di convergenza diminuisce all'aumentare di

mentre aumenta la distanza tra errore e stimatore.

Estensione di Newton a sistemi di equazioni non lineari

Dato un sistema di equazioni non lineari

Cercare le radici significa cercare dove ogni

elemento

Tuttavia il metodo di Newton prevede una divisione, cosa non possibile tra vettori: per questo

motivo si utilizza la Jacobiana dove indice equazione e indice incognita

Newton cosi diventa

Regola di Cramer per sistemi lineari

Dato il sistema lineare con soluzione (con ) la -esima soluzione è

con matrice ottenuta sostituendo la -esima colonna di con .

Il problema di questo metodo è che se sfrutta la regola di Laplace per il calcolo del

determinante ha un costo computazionale elevatissimo(circa )

Metodi di risoluzione:

1. Diretti: numero finito di passi

2. Iterativi: numero infinito di passi

Metodo della Fattorizzazione LU

Questo metodo per la risoluzione di sistemi lineari crea due matrici triangolari( e ) a partire

da tali per cui .

è una matrice triangolare inferiore(lower)

è una matrice triangolare superiore(upper)

Inoltre, essendo il prodotto di due matrici, : se il

significa che non ci sono elementi nulli sulle diagonali di e di e non ci sono autovalori nulli

(per le proprietà del determinante).

Il metodo consiste nel trovare e e risolvere il sistema con alcune sostituzioni:

e ponendo risolvere prima e poi . Essendo e

triangolari i due sistemi possono essere risolti con i due metodi seguenti.

Metodo delle sostituzioni in avanti

Si utilizza per risolvere un sistema dove è conosciuta e triangolare inferiore mentre

è vettore colonna di termini noti. L'elemento -esimo della soluzione si trova con la formula

Costo computazionale della risoluzione (con ordine della matrice )

Metodo delle sostituzioni all'indietro

Si utilizza per risolvere un sistema dove è conosciuta e triangolare superiore mentre

è vettore colonna di termini conosciuti. L'elemento -esimo della soluzione si trova con

Costo computazionale della risoluzione (con ordine della matrice U)

Fattorizzazione LU

Per trovare le matrici e bisognerebbe risolvere il sistema di equazioni(con

dimensione della matrice ) con incognite(le componenti di e ) che risulterebbe

sottodeterminato. A fronte di ciò si è deciso di porre uguale a ogni elemento della diagonale

della matrice . Si calcolano i moltiplicatori (che verranno posti sulla riga e sulla colonna

della matrice ) facendo scorrere Il numero in apice tra parentesi indica la

fase della fattorizzazione(oltre ad indicare la riga base di partenza) in quanto la matrice , una

volta svolta la fattorizzazione, diventerà (quindi ). Partendo da

Calcolati gli si procede mettendo in posizione

In questo modo si annullano tutti gli elementi della prima colonna sotto la diagonale di .

A questo punto si incrementa e si ripete il procedimento del calcolo degli e degli

.

Alla fine del processo, che viene eseguito per si ottiene una matrice

triangolare superiore che altri non è che la matrice .

Il costo computazionale della fattorizzazione è

Metodo di eliminazione di Gauss(MEG)

1. Metodo LU: si scompone , si risolve con e poi si risolve

al fine di ottenere la soluzione di . Questa strada si sceglie se si deve risolvere

un sistema del tipo

in quanto eseguo una volta sola la fattorizzazione.

Ha un costo computazionale di

2. A orlata b e sostituzione all'indietro: si orla con e si fattorizza la matrice ottenuta,

risolvendo poi con il metodo delle sostituzioni all'indietro. Ha un costo computazionale

di

Il MEG viene utilizzato per il calcolo dell'inversa, risolvendo sistemi che hanno per termine

noto le colonne della matrice identità e le cui soluzioni sono le colonne della matrice inversa.

Condizione necessaria e sufficiente per fattorizzazione LU

Affinché la fattorizzazione LU di una matrice esista e sia unica è necessario e

sufficiente che le sottomatrici principali di con siano non singolari. Le

sottomatrici principali sono matrici ottenute eliminando dalla matrice iniziale lo stesso numero

di righe e di colonne, partendo sia per le righe sia per le colonne dalla fine(dal basso per le

righe, da destra per le colonne). Nel caso qualche sottomatrice sia singolare si ha che le

fattorizzazioni LU possono essere o infinite o non esistenti per la matrice iniziale .

Condizione sufficiente per fattorizzazione LU

Affinché la fattorizzazione LU di una matrice esista e sia unica è sufficiente che

1. sia simmetrica definita positiva(sdp)

2. sia a dominanza diagonale stretta per righe, ovvero

3. sia a dominanza diagonale stretta per colonne, ovvero

Pivoting

Nella fattorizzazione LU gli elementi vengono ottenuti mediante una divisione il cui

denominatore è

. Tuttavia può capitare che questo elemento sia nullo, portando ad

un'operazione impossibile. A fronte di ciò si utilizza la tecnica del Pivoting che consiste nello

scambiare alcune righe(pivoting per righe) o alcune colonne(pivoting per colonne) della

matrice . Per farlo si utilizza una matrice ortogonale caratterizzata dall'avere ogni riga e

ogni colonna con tutti gli elementi nulli tranne uno pari ad 1. In questo modo la fattorizzazione

diventa e il sistema si risolve con e .

Questa tecnica, oltre ad evitare pivot nulli(elementi

), aiuta nell'accuratezza dei calcoli in

quanto evita l'utilizzo di numeri piccoli. Ad es.:

Sebbene le condizioni per attuare LU siano soddisfatte, senza pivoting si incorre in questo

errore quando gli

piccoli sono al denominatore, portando ad amplificazioni significative.

Il pivoting si dice parziale se fatto solo su una riga/colonna mentre si dice totale se viene

effettuato su tutta la matrice.

Matrice sparsa

Una matrice si dice sparsa quando il numero delle entrate non nulle è invece

di . In altre parole una matrice è considerata sparsa se il numero degli elementi nulli è

molto maggiore di quello degli elementi non nulli.

Matrice tridiagonale

Matrice con elementi non nulli su tre diagonali, solitamente sulla diagonale principale, sulla

prima sottodiagonale e sulla prima sopradiagonale.

Algoritmo di Thomas

Per risolvere un sistema con tridiagonale si può utilizzare l'algoritmo di Thomas.

Data

si fattorizza in

Il costo della fattorizzazione è di calcoli.

La risoluzione del sistema diventa quindi

Il costo per la risoluzione è di e quindi tutto l'algoritmo costa calcoli

Algoritmo di Cholesky

Per risolvere un sistema con sdp si può utilizzare l'algoritmo di Cholesky che prevede

una fattorizzazione di come dove e dove in questo caso gli

elementi diagonali di non valgono 1. Questo non comporta ad un sistema sotto determinato

in quanto le incognite sono solo . Gli elementi di si trovano quindi con l'algoritmo:

Efficacia della fattorizzazione LU

La fattorizzazione LU, per quanto possa essere buona, potrebbe non approssimare bene la

soluzione del sistema . Questo è dovuto al problema che si definisce mal

condizionato. Un esempio è la matrice di Hilbert che, per dimensioni di maggiori di 13,

presenta errori relativi maggiori del 1000%. Quando si risolve un sistema occorre tenere conto

di questi errori riscrivendo il problema come

considerando un problema ben condizionato se a piccoli e corrispondono piccoli .

Lemma 1 sulle matrici:

Data una matrice simmetrica definita positiva si ha che vale la relazione

Infatti se prendiamo autovalore di B e autovettore associato a abbiamo che

ed essendo risulta verificato il lemma.

Lemma 2 sulle matrici:

Data una matrice non singolare allora è simmetrica definita positiva.

Questo significa . Infatti se si prende allora

Posto si ha

Per ipotesi non è singolare quindi e dunque

Rapporto tra errore relativo di e

Dato un problema reale

Dimostrazione:

sottraiamo membro a membro ottenendo

dove abbiamo

Per il lemma 2 si ha che è sdp e quindi vale il lemma 1:

Moltiplicando per il denominatore si ha

Prendendo la prima disuguaglianza si ha

Sostituendo il risultato di prima

Per il lemma 1 si ha

Prendendo la seconda disuguaglianza

Per il lemma 2 e

Infine si ricavano, facendo la radice quadrata, le due disequazioni

Potendo essere moltiplicate tra loro membro a membro in quando ogni elemento è

strettamente positivo, si ha un'unica disequazione finale:

Questo risultato indica il legame tra gli errori relativi ed è utilizzata per il controllo tra gli errori

sui dati e quelli sulla soluzione. Il termine sotto radice dipende esclusivamente da A.

Condizionamento di una matrice

Numero che indica se un problema è condizionato bene(piccole variazioni sui dati causano

piccole variazioni sulla soluzione) o condizionato male. Il condizionamento di una matrice

viene definito come

Se è sdp la formula si riduce a

Questo numero dipende solo dal problema ed è sempre

Norme di matrici

1. Norma 1:

La norma 1 di una matrice è la massima somma ottenibile sommando gli elementi (in

valore assoluto) di una colonna della matrice.

2. Norma :

La norma 1 di una matrice è la massima somma ottenibile sommando gli elementi (in

valore assoluto) di una riga della matrice.

3. Norma 2:

Se

Raggio spettrale: è il valore assoluto(per numeri reali) o il modulo(per numeri complessi)

dell'autovalore massimo di una matrice. Si definisce che per una matrice sdp

diventa semplicemente . Infatti il condizionamento di una matrice si può scrivere

in quanto

Andamento degli errori relativi di un problema perturbato

1. Se

2. Se

Quando abbiamo che

Allora risulta, come da richiesta del denominatore, che

Per grande si ha un problema mal condizionato.

Per piccolo si ha un problema ben condizionato.

Per si ha un buon sistema, abbastanza ben condizionato.

Residuo in un sistema di equazioni lineari

Data una soluzione ad un sistema di equazioni lineari si ha che il residuo è

Perciò risulta, ricordando che

Il residuo piccolo comporta ad avere l'errore relativo sulla soluzione piccolo solo se anche il

condizionamento è piccolo.

Metodi iterativi

Un metodo iterativo è uno schema numerico per la risoluzione di un sistema dove il

vettore delle soluzioni per : ogni incognita viene iterata in una successione

dove

.

Lo schema iterativo viene attuato dalla formula . Ogni schema iterativo è

identificato dalla matrice quadrata (può indicare traslazione o rototraslazione) e dal vettore

(indica traslazione). La matrice viene detta matrice di iterazione. dipende da mentre

dipende da e da . Se è la soluzione esatta del sistema si ha che

ovvero

Proprietà degli schemi numerici

o Convergenza: proprietà dello schema ad approssimare in modo più corretto e preciso

all'aumentare delle iterazioni. Per si ha che soluzione esatta.

o Consistenza: proprietà che consiste nell'avere un problema con una approssimazione

che soddisfi la soluzione esatta, ovvero che la soluzione esatta sia soluzione anche del

sistema approssimato. Data soluzione esatta, se un problema è ben consistente allora

vale .

o Stabilità: uno schema numerico si dice stabile quando piccole perturbazioni sul

problema comportano piccole perturbazioni sui risultati(soluzioni)

Rapporto tra l'errore di due iterate consecutive

Posto che per la proprietà di convergenza e per la proprietà di

consistenza, si ha che se sottraggo la seconda equazione dalla prima ottengo

e considerando l'errore risulta che

Convergenza di un metodo iterativo

Condizione sufficiente affinché uno schema iterativo converga è che . Condizione

necessaria alla divergenza di un schema iterativo è che .

Queste condizioni valgono per ogni matrice .

Dimostrazione:

Da consegue che

dato che

Il valore massimo di e perciò l'errore si può scrivere

considerando quindi ( se è sdp)

Di conseguenza si crea una successione di disuguaglianze

Arrivando a concludere che con

Minore è il raggio spettrale e maggiore è la velocità di convergenza.

Dalla relazione si ricava il numero di iterate minimo a garantire l'abbattimento

dell'errore relativo rispetto alla tolleranza fissata.

Consistenza e convergenza di un metodo iterativo

Condizione necessaria e sufficiente affinché uno schema iterativo converga è che

e che lo schema sia consistente.

Precondizionatore

Matrice utilizzata nel metodo iterativo di Richardson che non deve essere singolare. Lo

scopo di questa matrice è di migliorare il condizionamento di A e di velocizzare la convergenza.

Per poterlo fare si utilizzano matrici particolari come diagonali, triangolari o tridiagonali.

Risoluzione di un problema con uno schema iterativo: metodo di Richardson

Si definisce il precondizionatore tale per cui . Partendo da si

sostituisce trovando . Associamo il primo membro all'iterata

mentre il secondo all'iterata in modo da avere la proprietà di convergenza; la proprietà di

consistenza è invece verificata per costruzione( è una forma consistente).

Troviamo cosi che e che

In alternativa si può risolvere lo schema utilizzando il residuo: si parte da

dove si identifica il residuo e quindi

Lo schema quindi vien attuato calcolando

dove viene detto residuo precondizionato.

Per trovarlo senza fare l'inversa di basta risolvere il sistema ad ogni iterata.

Per velocizzare lo schema e la convergenza si può introdurre un parametro reale tale che

e sviluppando come prima(passaggi analoghi al caso senza )

Uno schema di Richardson viene definito in due modi a seconda di

1. Stazionario: se costante

2. Dinamico: se

Metodo di Jacobi

Il metodo di Jacobi è uno schema di Richardson parallelizzabile in quanto il calcolo delle iterate

successive può essere fatto simultaneamente da più processori:

Dato il sistema

Si risolve isolando da ogni -esima equazione la componente dell'iterazione, legandola alla

nuova iterazione e legando le altre componenti di all'iterazione precedente

Ad esempio, nel caso di si ha che

Se si prende la matrice diagonale fatta dalla diagonale di la formula generale diventa

che, considerando tutti i valori di , può compattarsi in

L'ultima forma della formula indica che Jacobi è un metodo di Richardson stazionario con

e .

La matrice di iterazione risulta mentre si ha che . In

altre parole la matrice possiede tutti gli elementi di ,tranne la diagonale che è nulla,

cambiati di segno e divisi per con che indica la riga (

).

Condizione necessaria e sufficiente alla convergenza di Jacobi è che

Teorema sulla convergenza di Jacobi

Se la matrice è a dominanza diagonale stretta per righe(anche per colonne) allora il metodo

di Jacobi converge.

Dimostrazione:

La condizione richiesta da Jacobi è sempre verificata per matrici a dominanza

diagonale stretta per righe/colonne. Prendiamo come autovalore e autovettore ad esso

associato della matrice quindi . Considerando ogni riga si può riscrivere che

ovvero che

dove è la componente -esima

dell'autovettore e è l'elemento . Definiamo per semplicità che

( ) il quale può essere considerato vero in quanto ogni autovettore può essere

definito a meno di una costante moltiplicativa arbitraria.

Perciò si ha(considerando che gli elementi diagonali di sono nulli)

La disuguaglianza è dovuta a per definizione di che è la

componente massima. Per lo stesso motivo risulta che .

Dato che

e visto che per ipotesi di dominanza si ha

risulta che

la quale, ricollegata alla prima serie di uguaglianze/disuguaglianze, porta a

ovvero che e quindi alla convergenza del metodo.

Metodo di Gauss-Seidel

Il metodo di Gauss-Seidel è molto simile a quello di Jacobi, con la differenza che le iterate non

possono essere parallelizzabili in quanto ogni iterazione utilizza le componenti più aggiornate

della successione di iterazione . Questo vuol dire che, per il caso si ha

La formula generalizzata diventa quindi

Definendo gli elementi della matrice

e ricordando la matrice diagonale di (Jacobi) si può compattare la formula in

Perciò il precondizionatore che altri non è che la matrice triangolare inferiore

ottenuta rendendo triangolare (sostituendo con 0 ciò che non è sotto o nella diagonale)

Condizione sufficiente per convergenza di Gauss-Seidel è che sia sdp oppure che sia a

dominanza diagonale stretta per righe/colonne(anche entrambe le condizioni insieme).

Matrici tridiagonali e metodi di Jacobi/Gauss-Seidel

Se è tridiagonale, non singolare e se allora si possono avere due casi

distinti:

1. Entrambi gli schemi convergono con

Quindi Gauss-Seidel converge all'incirca il doppio velocemente rispetto a Jacobi.

2. Entrambi gli schemi divergono

Ad esempio, in ipotesi di convergenza, si ha

Quindi se

Essendo

I risultati ottenuti indicando che il numero di iterate minimo da eseguire per raggiungere la

tolleranza con Gauss-Seidel è circa la metà rispetto a Jacobi

Criteri d'arresto per schemi iterativi

Data la successione si arresta il metodo quando

1. Residuo:

Si definisce come il valore di

Il residuo iniziale è quindi . La relazione che lega il residuo al risultato è

Ricordando che

Risulta quindi che il residuo

è un buon stimatore per piccolo

2. Incremento:

Si definisce come il valore di

Per la disuguaglianza triangolare

Con la matrice sdp si ha che per convergenza. Dunque

Dall'ultima disuguaglianza si nota che, minore è è più affidabile diventa lo

stimatore basato sull'incremento. Da notare inoltre che, mentre nel caso del residuo si

controlla un errore relativo, qui si controlla un errore assoluto.

Richardson con convergenza ottimale

Se la matrice (del problema ) e la matrice sono sdp allora Richardson stazionario

converge se e solo se

dove è l'autovalore massimo della

matrice . Inoltre il fattore ottimale è

sempre riferito a

Dimostrazione:

Avendo ipotizzato convergenza si ha che dove è la matrice di

iterazione generica di Richardson. Lo schema diventa quindi

La matrice ha tutti gli autovalori strettamente positivi in quanto e sono sdp.

Definiamo come lo spettro degli autovalori di dove è l'autovalore

-esimo della matrice (ricordando che è l'autovalore -esimo di ).

Per l'ipotesi di convergenza si ha come conseguenza che

Essendo si può dedurre dalla prima disequazione che mentre dalla seconda

risulta che

.

Se per comodità prendiamo gli autovalori si ha che

è vera

solo se

in quanto

. Avendo preso si ha che

dimostrando la convergenza . Lo scopo è trovare tale che

sia più piccolo possibile.

Per dimostrare il valore

ottimale di si considera un

grafico con ascisse e

ordinate .

Ipotizzando sempre per

comodità

otteniamo la

figura qui a lato.

Dall'intersezione di

e rispettivamente

legati a e si trova il

valore minimo di possibile

al quale corrisponde in ascissa

al valore ottimale di .

Richardson stazionario

Matrice di iterazione

Autovalori di

Il raggio spettrale di indica la massima velocità di convergenza ottenibile in questo caso.

Richardson dinamico

Con e sdp si ha convergenza di Richardson se il parametro d'accelerazione è

Inoltre si ha che il rapporto tra l'errore iniziale e il -esimo è

La convenienza nell'utilizzo del metodo dinamico è che nel calcolo di compaiono termini

già calcolati in altri punti dell'algoritmo mentre nel caso stazionario bisogna calcolare anche gli

autovalori e cercare minimo e massimo, aumentando il costo computazionale di molto.

è la norma dell'energia e dipende dal problema(ovvero dalla matrice )

Metodo del gradiente precondizionato

L'algoritmo del gradiente precondizionato prevede i seguenti passi

1. Parto da e

2.

3. Trovo da

4. Calcolo

5.

6.

7.

8. Controllo con il criterio d'arresto scelto e riparto dal punto 3

Il residuo calcolato ad ogni iterata può essere scritto come

Per poter attuare questo metodo occorre avere una matrice non singolare, piuttosto

semplice(diagonale o simile) e il condizionamento deve essere abbastanza ridotto,

cosa spesso difficile da avere con matrici "semplici". Questo perché l'errore viene ridotto ad

ogni iterata dal fattore

che deve essere al fine di garantire convergenza.

Metodo del gradiente

Il metodo del gradiente (NON precondizionato) si utilizza solo per matrici sdp ed è simile a

quello precondizionato, con due differenze fondamentali:

1.

infatti si nota che non si utilizza più il residuo precondizionato

2. Il rapporto tra errore -esimo e quello iniziale diventa

Lo svantaggio dell'utilizzo di questo metodo è che risulta impossibile cambiare in quanto

dipende unicamente dal problema. Anche la mancanza di impedisce l'aumento della velocità

di convergenza, causando tempi più lunghi di risoluzione.

Metodo del gradiente e formula dell'energia

Questo metodo sfrutta una forma quadratica detta energia del sistema :

con il primo termine detto energia esterna e il secondo energia interna.

Il sistema con sdp allora se e solo se è soluzione del sistema ne consegue che

minimizza .

Essendo un paraboloide identifica il vertice più basso, ovvero il punto di minimo.

Dimostrazione

Definiamo il gradiente di

Se minimizza significa che

Se invece partiamo dall'ipotesi che possiamo scrivere che

e definendo

Ricordando la proprietà possiamo scrivere

Essendo

si ha che

rendendo cosi punto di minimo.

Norma dell'energia

Dato un vettore e una matrice si definisce norma dell'energia

Direzione del gradiente

Definiamo la direzione da prendere ad ogni valore di per raggiungere nel miglior

modo possibile(più in fretta) il minimo di un paraboloide partendo da un punto qualsiasi di

esso. Seguendo uno schema iterativo avremo che

ovvero

Come troviamo e ?

Essendo il gradiente la direzione di massima crescita allora la direzione sarà l'opposta al

gradiente. Perciò .

Per trovare ottimale occorre trovare un punto stazionario di rispetto alla variazione di e

quindi, dato che definiamo:

Ponendo

risulta che

ovvero

Isolando si raggiunge il risultato

Vettori A-ortogonali

Dati due vettori e una matrice sdp si dice che sono A-ortogonali se

risulta essere verificata.

Metodo del gradiente coniugato

Definiamo come direzione di discesa generica, dove quella iniziale è e dove

Il metodo prevede i seguenti passi

o

o

o

o

o

o

parametro che serve per dare la A-ortogonalità tra

o

o Si incrementa , si controlla con lo stimatore scelto se continuare. Se si, si riparte dal

punto 3 con la nuova iterata.

Questo metodo è molto efficace per 3 motivi:

1. Se è sdp allora il metodo converge in un numero finito di passi uguale o inferiore a

2. L'errore è ortogonale a ad ogni passo

3.

con

dove è il condizionamento fatto sulla norma 2

Approssimazione di funzioni e dati

Significa trovare funzioni semplici che approssimino delle funzioni reali o un insieme di dati in

maniera abbastanza precisa.

Sviluppo di Taylor

Interpolazione e approssimazione ai minimi quadrati

È un metodo che serve a costruire la funzione approssimata in modo che riproduca quei

valori sui quali è costruita. Si utilizza per quando si ha tanti dati e consiste nel mimare

l'andamento medio dei dati proposti(ogni volta minimizzando la somma degli scarti quadratici)

Interpolazione

Definiamo le coppie d'interpolazione dove sono i nodi di interpolazione(tutti

diversi tra loro, ovvero ) e sono i valori da interpolare. I valori

devono valere anche per la funzione quindi in ogni nodo.

La funzione si dice funzione interpolante e ci sono due tipologie di interpolazione:

o Interpolazione polinomiale: dove dipende dai

dati(specialmente dalla quantità di dati posseduti)

o Interpolazione trigonometrica: è una combinazione lineare di oggetti legati alla

trasformata di Fourier.

o Interpolazione razionale: è un rapporto tra polinomi

Unicità dell'interpolazione polinomiale

con distinti

dove viene chiamato polinomio interpolatore( dati si interpolano di grado ).

Considerando la funzione "reale" dell'andamento dei dati chiamata (continua nell'intervallo

dei dati considerati) si può definire il polinomio di interpolazione legato ad essa come

Dimostrazione:

Per assurdo poniamo che esistano due polinomi di interpolazione diversi per la stessa

funzione : e . Ne segue che

Definiamo un terzo polinomio con

Dall'ultima relazione si nota che il polinomio ha radici pari a zero. L'unico

polinomio di grado con radici pari a 0 è il polinomio .

Questo impone quindi che provando che il polinomio interpolatore è

unico.

Delta di Kronecker e funzione caratteristica

Si definisce Delta di Kronecker la funzione

In altre parole la funzione prende in ingresso una successione di valori e restituisce 0 per tutti i

valori tranne per un caso( ) in cui da risposta unitaria.

La funzione caratteristica è analoga al in quanto viene definita come

Restituisce 1 se l'oggetto appartiene all'insieme altrimenti l'output è nullo.

Costruzione del polinomio interpolatore

Prendiamo un insieme di nodi tali che siano legati ai

dati mentre per si abbia .

La funzione indica che con quindi si può scrivere

E in forma compatta diventa

viene detto polinomio caratteristico(in analogia alla funzione caratteristica)

Considerando ci si ritrova con un sistema a

equazioni e incognite. Tuttavia si nota che per in quanto i

termini diversi da sono tutti nulli. Si trova cosi che, per , ogni coefficiente

incognito vale esattamente arrivando a concludere

che in forma compatta diventa

Questo è il polinomio di interpolazione nella forma di Lagrange.

Errore tra polinomio reale e approssimato

L'errore viene definito come .

L'errore di interpolazione è indicato invece come .

Una proprietà dell'errore è la seguente:

Dato l'insieme aperto e le coppie con nodi distinti abbiamo che

Notare che la forma assomiglia molto allo sviluppo di Taylor, che è un polinomio di

grado e che .

Supponiamo un equispaziatura tra i vari nodi che valga

Quindi si ha che il passo è

rappresenta quindi una partizione uniforme di . Da questo possiamo dire che

Se sostituiamo questa disuguaglianza nell'espressione dell'errore di interpolazione abbiamo

Perciò, massimizzando su abbiamo che

L'ideale è avere un : ciò accade se

o in quanto diventa sempre più piccola all'aumentare di (a

maggior ragione quando )

o Se oppure con costante allora . Il problema sorge quando

, specie quando lo fa più velocemente rispetto alla rapidità di a tendere a

Fenomeno di Runge

Fenomeno di amplificazione dell'errore che causa oscillazioni di ampiezza piuttosto elevata

nella funzione interpolante rispetto a quella reale: questo è principalmente dovuto

all'equispaziatura dei nodi. Le oscillazioni si dicono spurie(non fisiche) e sono ben visibili nelle

zone estreme della funzione(ad esempio)

Infatti, per (e per questa funzione)si ha che

Nodi di Chebyshev

Uno dei metodi per risolvere il fenomeno di Runge è l'utilizzo di alcuni nodi speciali detti nodi

di Chebyshev. Dato un intervallo si ha che

Definiamo la funzione

.

La funzione indica il legame tra il dominio e il dominio dei nodi di interpolazione:

o

o

o

Per trovare i nodi si può costruire la

semicirconferenza di raggio unitario qui a lato.

La semicirconferenza viene divisa in

parti( nel disegno):

i punti scritti in nero sulla curva rappresentano

i nodi equispaziati(infatti gli archi di

circonferenza sono uguali) mentre i punti

cerchiati con indice sottolineato rappresentano

i valori di da sostituire in per poter ottenere

l'insieme dei nodi di Chebyshev.

Questo metodo, addensando i nodi agli estremi, può portare a oscillazioni considerevoli per

valori di piuttosto piccoli nelle zone poco "campionate" come per esempio quella centrale.

Tuttavia si ha che per un numero sufficientemente grande di queste oscillazioni si attenuano

molto e infatti la funzione

approssimata con ha un valore di .

La regolarità richiesta nell'utilizzo di Chebyshev è appena ovvero:

reale per (invece indica l'intervallo dei nodi senza

considerare quelli agli estremi e )

significa che si ha convergenza uniforme(mentre l'apice "c" indica Chebyshev)

Interpolazione composita

L'utilizzo dei nodi di Chebyshev è un buon metodo, tuttavia se si fanno delle misurazioni

sperimentali per un insieme di nodi e di dati non è detto che si trovino nodi di Chebyshev.

A fronte di ciò si utilizza un secondo metodo che previene il fenomeno di Runge: si chiama

interpolazione composita e utilizza una composizione di spezzate per approssimare una

funzione.

Interpolazione composita lineare

Definiamo le coppie

e gli intervalli

Applicando una restrizione su abbiamo l'espressione locale del polinomio che interpola le

due coppie di elementi e

La stima generica per l'errore del polinomio di grado è

Restringere sull' -esimo intervallo equivale porre , e perciò

Si conclude che se allora l'errore è maggiorato da

che tende

sempre a 0 per poiché è una costante e .

Spline: una spline è una polinomio interpolatore che ha la particolarità di essere più

regolare e liscia rispetto ad altri metodi di interpolazione già visti.

Spline cubica interpolatoria:

Definiamo le coppie di dati

. Le proprietà delle spline cubiche sono:

1.

2. e

3.

4.

5.

La prima proprietà indica condizione affinché la spline sia interpolatoria. Dalle proprietà

2,3 e 4 si hanno condizioni che impongono regolarità alla spline quindi .

L'ultima proprietà indica che ogni intervallo viene interpolato con un polinomio di terzo

grado e ciò comporta ad avere coefficienti incogniti al fine di definire la spline.

Tuttavia le condizioni totali a disposizioni(dalle prime quattro proprietà) sono e perciò

occorre imporre 2 condizioni arbitrariamente. Si può fare in due modi

1. Spline cubica interpolatoria naturale

Questa opzione consiste nell'imporre la derivata seconda del primo nodo e dell'ultimo

nulle, ovvero e

2. Spline cubica interpolatoria "not-a-knot"

Questa opzione consiste nell'imporre continuità della derivata terza nel secondo e nel

penultimo nodo:

e

La conseguenza è che a sinistra e a destra del secondo nodo(ciò accade anche al

penultimo) si ha la stessa cubica(ovvero stessi coefficienti) e perciò secondo e

penultimo nodo vengono definiti "not-a-knot"(non-nodi).

Errore sulle spline:

Si definisce errore il termine

Se allora

dove con passo -esimo tra due nodi(non è necessariamente

costante, quindi i nodi non sono necessariamente equispaziati)

La qualità dell'approssimazione delle derivate di considerando è solitamente migliore

dell'approssimazione della funzione nella spline . Tuttavia all'aumentare del grado delle

derivate diminuisce la qualità dell'approssimazione poiché H diminuisce di grado, ovvero:

Inoltre al crescere di cresce il parametro , rendendo l'approssimazione peggiore.

Il problema delle spline è che non conservano la monotonia della funzione di partenza, cosa

che si può risolvere con l'interpolazione di Hermite(che utilizza condizioni per costruire

l'approssimazione basate non solo sui valori della funzione ai nodi ma anche sui valori della

derivata prima sempre ai nodi).

Minimi quadrati

Estrapolazione che consiste nel trovare l' andamento di una funzione della quale si conosce

solo una nuvola di valori. Si differisce dall'interpolazione in quanto non è necessario che

l'approssimazione passi per i valori associati ai nodi. La funzione approssimata è un

polinomio di grado solitamente molto minore del numero di dati da utilizzare.

Si dice approssimazione ai minimi quadrati la

L'espressione indica lo scarto quadratico tra un generico polinomio di grado

e il generico dato . La disuguaglianza con le sommatorie indica che, tra tutti i polinomi di

grado , esiste un polinomio la cui somma di tutti gli scarti quadratici è minore(o al più

uguale) della somma di un qualunque altro polinomio dello stesso grado.

Criteri di scelta del grado del polinomio per i minimi quadrati

1. Minore è e più è veloce il metodo dei minimi quadrati

2. A seconda di come è disposta la nuvola di valori può essere più opportuno scegliere

una retta, una parabola …(guardando un grafico con i dati)

3. Mi faccio aiutare dalla fisica del problema(ad esempio per la forza di una molla mi

aspetto un andamento lineare rispetto l'allungamento e quindi pongo )

Coefficienti del polinomio per minimi quadrati

Definiamo come generico polinomio di grado .

Definiamo come il polinomio migliore per approssimare i dati.

Si dice funzione della somma degli scarti quadratici

Lo scopo è trovare i valori che corrispondono ai valori minimi di .

Per il caso si ha che e

L'ultima sommatoria rappresenta un paraboloide e perciò per trovare e devo utilizzare le

derivate parziali e uguagliarle a zero in e

Perciò basta risolvere il seguente sistema per ottenere e :

Per il caso generale il sistema da risolvere è

Notare che nel caso con numero di dati si ritorna ad un interpolazione(in quanto gli

scarti quadratici minimi corrispondono a zero).