Introduzione

Questi appunti sono un estratto del corso di Analisi Numerica tenuto dalla Professoressa R. M. Pidatella (anno 2021/22) al corso magistrale di informatica dell'Università degli Studi di Catania. Gli appunti non sostituiscono il materiale didattico e consiglio fortemente di approfondire gli argomenti trattati a lezione con gli adeguati libri di testo.

La repository è open su GitHub a questo link. Se trovi errori all'interno del testo, segnalali su GitHub attraverso una issue. Se gli appunti sono stati d'aiuto, valuta di lasciare una star alla repository. Per maggiori informazioni, contattami via email.

Numan - Implementazioni Python

Alcuni degli algoritmi sono implementati in Python nella repository integrativa Numan.

1. Teoria degli errori

Si vuole stimare l'influenza degli errori sulla soluzione, ovvero di quanto la soluzione calcolata $\hat{y}$ si discosti da quella reale $y$.

Errore assoluto

Calcola la bontà dell'approssimazione in termini di cifre decimali. $$E_a=|y-\hat{y}|$$

Errore relativo

Calcola la bontà dell'approssimazione in termini di cifre significative. È spesso rappresentato in percentuale ed utilizzato per avere una rapida intuizione della grandezza dell'errore. $$E_r=\left|\frac{y-\hat{y}}{y}\right|$$

L'errore relativo è una migliore stima della bontà dell'approssimazione, tranne quando $y$ è prossimo a zero. Solitamente è più importante, poiché tiene in considerazione la grandezza dei numeri rispetto all'errore.

Approssimazioni corrette

Mentre con le cifre decimali ci riferiamo a tutte le cifre dopo la virgola, con cifre significative indichiamo le cifre necessarie ad esprimere il risultato di una misura con la precisione con cui è stata fatta.

• Diciamo che $\hat{y}$ è una corretta approssimazione di $y$ a $p$ cifre decimali se $E_a<10^{-p}$

  • $x=624.428731,y=624.420711,E_a=0.8\times 10^{-2}<10^{-2}\Rightarrow$ prime 2 cifre decimali uguali

• Diciamo che $\hat{y}$ è una corretta approssimazione di $y$ a $p$ cifre significative se $E_r<10^{-p+1}$

  • $x=624.428731,y=624.420711,E_r=1.3\times 10^{-5}<10^{-4}\Rightarrow$ prime 5 cifre significative uguali

Stime errore relativo

Il numero di cifre significative dà una stima dell'errore relativo. Consideriamo due numeri 990 e 110 con 2 cifre significative. L'errore sarà di $\pm 5$, quindi avremo che $E_r(110)\approx 5\%$ $E_r(990)\approx 0.5\%$, quindi il range in cui varia l'errore relativo è $.5÷5\%$. Si può stilare una tabella con la corrispondenza cifre significative - errore relativo.

Metodi iterativi

Lavoreremo spesso con metodi iterativi, che permettono di determinare una successione di soluzioni approssimate $\{x_n{\}}_n$ che converge alla soluzione esatta del problema $x^{\ast }$. Per ciascuna iterazione possiamo definire un errore come segue: $$e_n=|x_n-x^{\ast }|$$ Il metodo si dice convergente se si ha che: $$\underset{n\to \infty }{\lim }|e_n|=0\text{oppure}\underset{n\to \infty }{\lim }|x_n|=x^{\ast }$$

Ordine di convergenza

Dopo aver definito la nozione di convergenza di un metodo iterativo (def. dipendente dalla norma $|\cdot |$ utilizzata), si passa a definire l'ordine di convergenza. Diciamo che il metodo ha ordine di convergenza $p$ se: $$\underset{n\to \infty }{\lim }\frac{|e_{n+1}|}{|e_n{|}^p}=\text{costante}\text{oppure}|e_{n+1}|\le c|e_n{|}^p$$ La seconda vale definitivamente a partire da un $n$ sufficientemente grande e per qualche costante $c$. L'ordine di convergenza esprime il numero di cifre decimali che il metodo guadagna, ad ogni iterazione, rispetto alla soluzione esatta.

Altro sui metodi iterativi

Il metodo potrebbe iterare all'infinito, per cui bisogna stabilire un criterio di arresto. Per ogni tipo di problema è possibile progettare diversi algoritmi, che verranno valutati in base a:

• Stabilità
• Efficienza - Complessità computazionale
• Occupazione dello spazio

Esempio: Metodo di Gauss

Nel caso di un sistema lineare, il numero caratteristico del problema è il numero di equazioni del sistema che, in un sistema con matrice dei coefficienti quadrata, è uguale al numero delle incognite $n$. Nel caso del metodo di Gauss, che esamineremo più avanti, la complessità computazionale è $\frac{4}{3}{n}^3$.

Esempio: Schema di Horner

Lo schema di Horner serve ad abbassare il costo computazionale del calcolo di un polinomio. Si vuole eseguire il calcolo di: $$p(x)=a{x}^3+b{x}^2+cx+d$$ Sono necessarie 3 somme e 6 moltiplicazioni per un totale di 9 operazioni floating point (flops). In generale, per calcolare: $$p(x)=a_n{x}^n+a_{n-1}{x}^{n-1}+\cdots +a_{1}x+a_0$$ Occorrono $n$ somme ed $\frac{n(n+1)}{2}$ moltiplicazioni, per un totale di $$n+\frac{n(n+1)}{n}=\frac{n(n+3)}{n}\sim n^2\text{ flopls}$$ Se invece si applica lo schema di Horner: $$p(x)=a_0+x(a_1+x(a_2+x(\dots a_{n-1}+x{a}_n)\dots ))$$ Otteniamo $n$ somme ed $n$ moltiplicazioni per un totale di $2n$ flops.

Sistemi di numerazione

I sistemi di rappresentazione numerica sono posizionali, ovvero ogni cifra occupa una posizione corrispondente ad una potenza della base del sistema adottato. Una fonte di errore è data dal passaggio da un sistema di numerazione all'altro. Scelta una base $b$, ogni reale $a\in \mathbb{R}$ può essere scritto come $$a=\pm (a_m{b}^m+\cdots +a_0{b}^0+a_{-1}{b}^{-1}+\dots )$$ Dove $0\le a_i\le b-1$, $a_i\in \mathbb{N}$. La rappresentazione è univoca a meno che il numero non necessiti di infinite cifre consecutive $a_{-k}=b-1$ (in decimale ad esempio: 0.1299999...), allora una rappresentazione equivalente consiste nel sopprimere la successione aggiungendo un'unità all'ultima cifra rimasta (dall'esempio precedente: 0.13).

Teorema

Sia $b\in \mathbb{N}$ e $b\ge 2$ la base, sia $x\in \mathbb{R}\setminus \{0\}$ un numero reale, allora esiste unico $e\in \mathbb{Z}$ ed una successione $\{a_i\}$ di numeri naturali $a_i\in \mathbb{N}$ tale che $$x=\pm \left(\sum _{i=1}^{\infty }{a}_i{b}^{-i}\right)\cdot b^e$$ Dove

• $0\le a_i\le b-1$
• $a_1\ne 0$
• $a_i\ne b-1$ definitivamente

Piccolo hint: sembra somigliare alla notazione scientifica.

Rappresentazione numerica in un calcolatore

In un calcolatore lo spazio di memoria è finito, e la sommatoria precedente si può estendere fino ad un $t$ numero finito.

Rappresentazione in virgola fissa

In sintesi $N_1$ cifre prima della virgola, $N_2$ cifre decimali, $N=N_1+N_2$ cifre totali. Un bit $s$ è solitamente conservato per il segno.

Rappresentazione in virgola mobile

La posizione della virgola non è fissa, ma è data dall'esponente. Un numero di cifre $t$ è riservato alla mantissa $m$, un numero $N_e$ è riservato all'esponente $e$. Un bit $s$ è solitamente conservato per il segno. I bit della mantissa $m$ determinano la precisione con cui viene rappresentato un numero, mentre i bit dell'esponente $e$ determinano il massimo ed il minimo numero rappresentabili.

Insieme dei numeri macchina

L'insieme dei numeri macchina $F$ è l'insieme dei numeri reali che sono rappresentabili attraverso una mantissa da $t$ cifre e il cui esponente sta tra due interi $L$ (lower) ed $U$ (upper) definiti nel calcolatore. $$F(t,b,L,U)=\{0\}\cup \{x\in \mathbb{R}:x=\pm b^e\sum _{i=1}^t{a}_i{b}^{-i}\}$$ Lo 0 è rappresentato a parte poiché ha una rappresentazione particolare. L'insieme $F$ è finito e numerabile ed ha la seguente cardinalità: $$|F|=1+2(b-1)b^{t-1}(U-L+1)$$ La rappresentazione non è unica, ma si dice normalizzata se $a_1\ne 0$ e la mantissa $m\ge b^{t-1}$ (in pratica le regole utilizzate nella notazione scientifica).

Memorizzazione esponente

Tenendo conto che il lower bound $L$ dell'esponente è un numero negativo (minimo numero rappresentabile) allora per conservare un esponente $e$ si memorizza $e^{\ast }=e-L$, e dato che $e\ge L$ allora $e^{\ast }\ge 0$.

Conversione decimale-binario

La parte intera si divide per due: se c'è il resto si mette 1, altrimenti si mette 0 e si assegnano potenze di due crescenti.

La parte decimale si moltiplica per 2 (divide per 1/2), se il prodotto è minore di 1, allora si ha 0, altrimenti si ha 1. Quando il numero diventa maggiore di 1 si sottrae 1 procedendo come prima.

0.2 | 2 => 0 x 2^-1
0.4 | 2 => 0 x 2^-2 
0.8 | 2 => 1 x 2^-3
1.6 | 2 => 1 x 2^-4  
... (periodico)
0.2 = .0011

In questo caso abbiamo una rappresentazione finita in decimale ed una rappresentazione approssimata in binario, e ciò da luogo ad un errore di arrotondamento (round-off).

Conversione decimale a base qualunque

Sia $n$ un numero in base 10 che vogliamo convertire in base $b\in \mathbb{N}$. Sia $n=(a_j{a}_{j-1}\dots a_1{a}_0{)}_b$ la conversione attesa. Dividendo $n$ per $b$ sia ha: $$\frac{n}{b}=a_j\cdot b^{j-1}+\dots +a_1\cdot b^0+\frac{a_0}{b}⟹n=b{n}_0+a_0$$ Dove $n_0<n$. Quindi l'ultima cifra della rappresentazione $a_0$ non è che il resto intero di $n/b$. Per ottenere la penultima cifra si divide $n_0/b$ ed il procedimento si arresta ad $a_j$ tale che $n_j=0$.

Scegliere la rappresentazione in floating point

Dato $x\in \mathbb{R}$ come scelgo $fl(x)\in F$? Si hanno i seguenti casi:

• Underflow (esponente troppo piccolo) $e<L$
  • Si emette un warning e si pone $fl(x)=0$
• Overflow (esponente troppo grande) $e>U$
  • Segnale di errore e arresto del programma
• Se $L\le x\le U$, allora si procede:
  • Se $a_k=0$ per $k\ge t+1$ allora il numero $x$ è in $F$ e $x=fl(x)$
  • Altrimenti si ha:
    • (Chopping) Si esclude la parte destra della $t$-esima cifra
    • (Rounding) Si arrotonda in base alla $t$-esima cifra

Epsilon macchina

Tutti i floating point che hanno lo stesso esponente $e$ hanno una spaziatura $b^{e-t+1)}$, questo implica che due numeri con lo stesso esponente e consecutivi $p_1$ e $p_2$ sono legati dalla relazione $p_2=p_1+b^{e-t+1}$. La spaziatura cambia prima e dopo le potenze esatte della base. L'epsilon macchina è un upper bound dell'errore relativo che si ha arrotondando o troncando la rappresentazione reale per ottenere la rappresentazione in floating point. Dato che parliamo di errore relativo, possiamo studiare l'epsilon macchina limitandoci ad $e=0$ e ai numeri positivi.

Utilizzando il rounding, l'errore assoluto più grande ottenibile è $b^{e-t+1}/2$, basta porre un numero reale al centro dello spacing. Il denominatore nell'errore relativo è il numero che stiamo approssimando. Per trovare un upper bound dell'errore relativo, questo numero deve essere il più piccolo possibile. I peggiori errori relativi capitano arrotondando numeri del tipo $1+a$ dove $a\in [0,b^{-t+1}/2]$, nello specifico il peggior caso è $a=b^{-t+1}/2$. Calcolando l'errore relativo, il denominatore è $1+a\approx 1$, quindi abbiamo un errore relativo pari a $b^{-t+1}/2$. Questo upper bound è proprio l'epsilon macchina. Utilizzando il chopping, il ragionamento si ripete ma considerando $a\in [0,b^{-t+1}]$, quindi l'epsilon macchina è proprio pari a $b^{-t+1}$.

Aritmetica Standard IEEE

Wobbling precision

I numeri in floating point non sono equispaziati, ma si addensano in prossimità del più piccolo numero rappresentabile. All'interno dell'intervallo $[b^e,b^{e+1}]$ i punti sono equispaziati e la loro distanza è $b^{e-t}$. Ogni volta che si diminuisce o aumenta l'esponente $e$, diminuisce o aumenta la spaziatura. Il fenomeno è detto wobbling precision ed ha un andamento oscillatorio

Condizionamento e stabilità

Def. Consideriamo il problema: trovare $x$ tale che $F(x)=d$, dove $d$ è il dato (o i dati) da cui dipende la soluzione $x$ ed $F$ è la relazione funzionale che lega $x$ e $d$​. Diremo che tale problema è ben posto, se per un certo dato, la soluzione esiste, è unica e dipende con continuità dai dati.

L'unicità e l'esistenza della soluzione sono problemi analitici, mentre la dipendenza è un problema numerico. La dipendenza continua dei dati significa che piccole perturbazioni dei dati danno luogo a piccole variazioni della soluzione, dove "piccolo" può essere inteso in senso relativo o assoluto. Nascono due problemi, ci si chiede:

1. Alterando i dati del problema, di quanto si altera la soluzione?
2. Come si propagano gli errori?

Numero di condizionamento

Il primo problema è connesso con la dipendenza continua dai dati della soluzione e può essere stimato con il numero di condizionamento del problema, numero che non dipende dall'uso dell'aritmetica finita del calcolatore, ma dal tipo di problema.

Def. Sia $\delta d$ la perturbazione applicata ai dati, che scatena una variazione $\delta x$ nella soluzione, quindi $d+\delta d\to x+\delta x$. Diremo $K$ numero di condizionamento relativo il valore ottenuto come segue: $$K=\frac{\Vert \delta x\Vert /\Vert x\Vert }{\Vert \delta d\Vert /\Vert d\Vert }$$ Se $x=0,d=0$ si calcola $K_{ass}$ il numero di condizionamento assoluto: $$K_{ass}=\frac{\Vert \delta x\Vert }{\Vert \delta d\Vert }$$

Se $K$ è grande allora il problema è mal condizionato. Se un problema è ben posto ma $K$ è grande basta riformulare il problema.

Propagazione dell'errore

Il problema (2) dipende dalla stabilità dell'algoritmo. Ad ogni problema numerico si possono associare più algoritmi.

Un algoritmo è stabile se la propagazione degli errori dovuti all'aritmetica di macchina è limitata.

Un algoritmo è più stabile di un altro se in esso l'influenza degli errori è minore.

Errori dovuti a operazioni aritmetiche

In generale, i risultati di operazioni aritmetiche tra numeri macchina non sono numeri macchina. Indichiamo con $a\oplus b$ il valore reale di un calcolo, dove con $\oplus$ indichiamo una tra le 4 operazioni. Quando tale risultato non appartiene ai numeri macchina $F$, indichiamo la sua approssimazione con $fl(a\oplus b)$. Chiamiamo round-off error la differenza $$a\oplus b-fl(a\oplus b)$$ Tale errore si può minimizzare con degli accorgimenti.

Esempio con somma. Quando un calcolatore somma due numeri, porta il secondo numero allo stesso ordine del primo (stesso esponente) e somma la mantissa (da $m$ cifre) in un accumulatore da $2m$ cifre, dopodiché aggiusta l'esponente. Questo può portare alla non validità della proprietà commutativa, distributiva ed associativa della somma e della moltiplicazione.

Esempio numerico. Supponiamo $m=4$ e di voler eseguire la somma ${\sum }_{i=1}^{11}{x}_i$, dove $x_1=0.5055\times 10^4$ e per $i=2,\dots ,11$ abbiamo $x_i=0.4000\times 10^0$. Nella prima somma avremo $x_1+x_2$, quindi portiamo tutto all'ordine di $x_1$: $(0.5055+0.00004)\times 10^4$, ma essendo che la mantissa ha 4 cifre, la somma che viene calcolata è $(0.5055+0.0000)\times 10^4$ che rimane $x_1$. Questo non avviene calcolando prima la somma tra i termini $x_2,\dots ,x_{1}1$ e dopodiché $x_1$, quindi non vale la proprietà commutativa. Per minimizzare l'errore conviene sommare i numeri in base alla loro precisione, dal più piccolo al più grande (in valore assoluto).

Se $fl(a\oplus b)$ è arrotondato correttamente (rounding, come nel sistema IEEE), allora a meno di underflow o overflow si ha che: $$fl(a\oplus b)=(a\oplus b)(1+\delta )\text{con }|\delta |<{ϵ}_M$$ In realtà, per evitare errori dovuti ad underflow si aggiunge una piccolissima quantità ($10^{-45}$) che da luogo ad una "eccezione" che viene segnalata da un flag.

Errori nelle 4 operazioni

Consideriamo l'errore nella forma precedente, per l'errore assoluto si ha:

• Somma algebrica: $x_1(1+{ϵ}_1)\pm x_2(1+{ϵ}_2)$
• Prodotto: $x_1(1+{ϵ}_1)\cdot x_2(1+{ϵ}_2)$
• Divisione: $x_1(1+{ϵ}_1)/x_2(1+{ϵ}_2)$

Per l'errore relativo si ha:

Prodotto: $$\frac{x_1{x}_2-x_1(1+{ϵ}_1)x_2(1+{ϵ}_2)}{x_1{x}_2}=-{ϵ}_1-{ϵ}_2-{ϵ}_1{ϵ}_2\approx -{ϵ}_1-{ϵ}_2$$ Divisione: $$\frac{\frac{x_1}{x_2}-\frac{x_1(1+{ϵ}_1)}{x_2(1+{ϵ}_2)}}{\frac{x_1}{x_2}}=\frac{-{ϵ}_1+{ϵ}_2}{1+{ϵ}_2}\approx -{ϵ}_1+{ϵ}_2$$ Somma algebrica: $$\frac{(x_1\pm x_2)-[x_1(1+{ϵ}_1)+x_2(1+{ϵ}_2)]}{x_1\pm x_2}=\frac{-x_1{ϵ}_1-x_2{ϵ}_2}{x_1\pm x_2}$$ Pertanto abbiamo che:

• La sottrazione va bene per l'errore assoluto, ma si ha un errore relativo grande se $x_1\approx x_2$.
• Il prodotto va bene per l'errore relativo, l'errore assoluto dipende dall'ordine di grandezza dei fattori.
• La divisione va bene per l'errore relativo, l'errore assoluto è grande se $x_2\approx 0$.

2. Richiami sulle matrici

Inserirò solo teoremi e definizioni non di base.

Traccia

La traccia di una matrice è definita solo per matrici quadrate ed è la somma degli elementi della diagonale principale.

$$tr(A)=\sum _i{a}_{ii}$$

Proprietà sulla traccia

Prop 1. Calcolare la traccia è una trasformazione lineare:

$$tr(A+B)=tr(A)+tr(B)tr(cA)=c\times tr(A)$$

Prop 2. Una matrice e la sua trasposta hanno la stessa traccia

$$tr(A)=tr(A^T)$$

Prop 3. Data una matrice $A$ di dimensione $m\times n$ ed una matrice $B$ di dimensione $n\times m$ allora

$$tr(AB)=tr(BA)$$

Prop 4. La traccia è pari alla somma degli autovalori della matrice.

Vedi tutte le propietà.

Tipi di matrici

Esploriamo alcune definizioni particolari.

Matrice trasposta coniugata

La matrice A \in \C^{m\times n} con elementi complessi ha una trasposta coniugata $C=A^{\ast }$, ovvero una semplice trasposta dove ogni elemento è il coniugato dell'elemento corrispondente in $A$: $c_{ij}={\bar{a}}_{ij}$.

Matrice di Hessemberg

Una matrice di Hessemberg è una matrice triangolare superiore in cui è aggiunta una subdiagonale inferiore. La regola è che l'elemento $a_{ij}=0$ se $i>j+1$ oppure $j>i+1$. $$\left[\begin{array}{cccc}1& 4& 2& 3\\ 2& 2& 8& 9\\ 0& 7& 2& 4\\ 0& 0& 7& 2\end{array}\right]$$

Matrice definita positiva e negativa

Una matrice simmetrica si dice definita positiva se $\forall x\in {\mathbb{R}}^n$ e $x\ne 0$ si ha $x^{T}Ax>0$. Si dice definita negativa il viceversa con ($<$). Se si usa il minore-maggiore uguale, allora anziché definita la matrice è semidefinita.

Matrice diagonalmente dominante

Una matrice è strettamente diagonalmente dominante se per ogni riga, la somma dei valori assoluti degli elementi non in diagonale è minore stretta del valore assoluto dell'elemento in diagonale. Se anziché essere minore stretta è minore o uguale, allora la matrice si dice debolmente diagonalmente dominante.

Matrici unitarie, ortogonali e normali

• Matrice unitaria: $U^{\ast }U=U{U}^{\ast }=I$ da cui $U^{-1}=U^{\ast }$
• Matrice ortogonale: $U^{T}U=U{U}^T=I$ da cui $U^{-1}=U^T$
• Matrice normale: $U^{T}U=U{U}^T$

Teoremi sul determinante

Regola di Laplace per il determinante

Sia $A_{ij}$ la matrice privata dell'$i$-esima riga e della $j$-esima colonna. Si definisce complemento algebrico dell'elemento $a_{ij}$ il numero definito come segue: $$c_{ij}=(-1{)}^{i+j}\cdot \det (A_{ji})$$ Si definisce il determinante di $A$ sviluppato rispetto all'$i$-esima riga come segue: $$\det (A)=\sum _{j=1}^n{a}_{ij}\cdot c_{ij}$$

Teorema di Binet

Date due matrici moltiplicabili $A,B$ si ha che $\det (AB)=\det (A)\cdot \det (B)$.

Teorema di Sylvester

Indichiamo con $A_k$ la matrice $A$ formata dalle prime $k$ righe e $k$ colonne. Una matrice $A$ simmetrica si dice definita positiva se e solo se $\det (A_k)>0$ per $k=1,\dots ,n$.

Corollario del teorema di Sylvester

Sia $A$ una matrice simmetrica, diagonalmente dominante a diagonale positiva, allora la matrice è definita positiva.

Matrice non degenere

Una matrice non degenere è una matrice con determinante non nullo.

Teoremi sulla matrice inversa

Costruzione della matrice inversa

Sia $\hat{A}$ la matrice trasposta dei complementi algebrici, tale matrice gode della seguente proprietà: $$A\hat{A}=\hat{A}A=\det (A)\cdot I_n$$ Questo implica che per ottenere l'inversa $A^{-1}$ di $A$ basta calcolare $$\frac{\hat{A}}{\det (A)}=A^{-1}$$ Per calcolare direttamente ogni elemento dell'inversa $b_{ij}$ basterà calcolare come segue: $$b_{ij}=(-1{)}^{i+j}\frac{\det (A_{ji})}{\det (A)}$$

Unicità dell'inversa

La matrice inversa quando esiste è unica.

Dimostrazione (per assurdo). Supponiamo che esista una matrice $B$ diversa da $A^{-1}$ e tale che $BA=I$. Allora si ha $BA{A}^{-1}=I{A}^{-1}⟹B=A^{-1}$, assurdo $\perp$.

Inversa del prodotto tra matrici

Siano $A,B$ non degeneri e sia $C=AB$, che per il teorema di Binet è non degenere. Allora si ha che $C^{-1}=B^{-1}{A}^{-1}$ (si noti l'inversione).

Norme vettoriali e matriciali

Norma p o norma Holderiana

Sia $1\le p\le \infty$, la norma $p$ è definita come segue:

|x|p = \root{p}\of{\sum{i} | x_i |^p}

Definiamo le norme più comuni come segue:

\begin{split}

& |x|1 = & \sum_i |x_i| \hspace{1cm} & \text{norma 1} \ & |x|2 = & \sqrt{\sum_i x_i^2} \hspace{1cm} & \text{norma 2} \ & |x|\infty = & \max{i} |x_i| \hspace{1cm} & \text{norma }\infty \

\end{split}

Norme equivalenti in spazi discreti

In ${\mathbb{R}}^n$ le norme $1,2,\infty$ sono equivalenti, ovvero esistono $\alpha ,\beta \in \mathbb{R}$ con $0<\alpha <\beta$ tali che: $$\alpha \Vert x{\Vert }_a\le \Vert x{\Vert }_b\le \beta \Vert x{\Vert }_a$$ Dove $a$ e $b$ possono essere due norme distinte tra le sopracitate.

Norme matriciali indotte

Credits: Advanced LAFF, Robert van de Geijn

Supponiamo di avere una trasformazione lineare L: \C^n \to \C^m e supponiamo di voler misurare quanto il vettore x \in \C^n venga magnificato dalla trasformazione. Per avere una sorta di magnitudo del vettore utilizziamo la norma vettoriale $||x||$. Per capire il grado di magnificazione dato dalla trasformazione $L$ potremmo considerare il rapporto tra la norma del vettore di output $y=L(x)$ e quella del vettore di input, come segue: $$\frac{||L(x)||}{||x||}$$ Supponiamo ancora che per descrivere $L$ vogliamo trovare $x$ che ci dia il massimo rapporto: $$\underset{x\ne 0}{\sup }\frac{||L(x)||}{||x||}$$ Abbiamo quindi una definizione della magnificazione data da $L$, quindi potremmo costruire la norma di $L$ basandoci su questo concetto: $$||L||=\underset{x\ne 0}{\sup }\frac{||L(x)||}{||x||}$$ Se abbiamo una matrice $A$ che definisce la trasformazione lineare $L$, ovvero $L(x)=Ax$, allora possiamo definire la norma matriciale di $A$ come: $$||A||=\underset{x\ne 0}{\sup }\frac{||Ax||}{||x||}$$ La norma matriciale di $A$ si dice indotta dalla norma vettoriale utilizzata. Ritornando alla trasformazione, se effettuiamo prima lo scaling e dopodiché trasformiamo, il risultato non cambia. Questo vuol dire che le due seguenti espressioni sono equivalenti: $$||L||=\underset{x\ne 0}{\sup }\frac{||L(x)||}{||x||}=\underset{x\ne 0}{\sup }\frac{||L(\frac{x}{||x||})||}{||\frac{x}{||x||}||}=\underset{x\ne 0}{\sup }\frac{||L(\frac{x}{||x||})||}{\frac{||x||}{||x||}}=\underset{x\ne 0}{\sup }||L(\frac{x}{||x||})||$$ Ma $\frac{x}{||x||}$ è un vettore unitario, quindi possiamo scrivere: $$||L||=\underset{||x||=1}{\sup }||normL(x)||$$ L'insieme $\{x\mid ||x||=1\}$ ha il massimo ed il minimo inclusi (compact set), quindi l'estremo superiore combacia con il massimo. $$||L||=\underset{||x||=1}{\max }||L(x)||$$ Possiamo effettuare esattamente gli stessi passaggi con la norma matriciale indotta, quindi: $$||A||=\underset{||x||=1}{\max }||Ax||$$ Dalle norme vettoriali $1,2,\infty$ vengono indotte le omonime norme matriciali:

\begin{split}

||{x}||1 = \sum_i \abs{x_i} & \to & ||{A}||1 = \max{1 \le j \le n} \sum{i=1}^n |a_{ij}|\ ||{x}||2 = \sqrt{\sum_i x_i^2} & \to & ||{A}||2 = \sqrt{\rho(A^* A)} \ ||{x}||{\infty} = \max{i} \abs{x_i} & \to & ||{A}||\infty =\max{1 \le i \le n} \sum_{j=1}^n |a_{ij}|\

\end{split}

Nota: nella norma 1 si prende la somma in valore maggiore assoluto della colonna, mentre nella norma infinito si prende la stessa somma, ma della riga.

Autovalori e autovettori

Definizione classica

Sia $A\in {\mathbb{R}}^{n\times n}$ diciamo che $\lambda$ è autovalore di $A$ se \exist \bar x \in \C^n che chiameremo autovettore, $\bar{x}\ne 0$, tale che $(A-\lambda I)\bar{x}=0$ o analogamente $A\bar{x}=\lambda \bar{x}$. Definiamo $\sigma (A)$ l'insieme degli autovalori di $A$.

Autovalore nullo

Una matrice contiene un autovalore nullo se e solo se è singolare ($\det (A)=0$).

Raggio spettrale

Il raggio spettrale è l'autovalore più grande in valore assoluto: $\rho ={\max }_{\lambda \in \sigma (A)}|\lambda |$.

Trasformata per contragradienza (simili)

Siano $A,B\in {\mathbb{R}}^{n\times n}$ ed esista $B^{-1}$ e sia $C=B^{-1}AB$, allora $C$ si dice trasformata per contragradienza di $A$ mediante $B$. Due matrici trasformate per contragradienza l'una dall'altra si dicono simili.

Proprietà matrici simili

Se $A$ e $C$ sono simili allora:

• $A^s$ e $C^s$ per $s\in \mathbb{N}$.
• $A^{-1}$ e $C^{-1}$ sono simili.
• se $\det (A)\ne 0$ allora anche $\det (C)\ne 0$

Matrice diagonalizzabile

$A\in {\mathbb{R}}^{n\times n}$ si dice diagonalizzabile se è simile ad una matrice diagonale.

6 rapide proprietà sugli autovalori

1. La matrice e la sua trasposta hanno gli stessi autovalori.
2. Due matrici simili hanno gli stessi autovalori.
3. Una matrice degenere ha almeno un autovalore nullo.
4. Se $\lambda$ è autovalore di $A$ ed esiste l'inversa, allora ${\lambda }^{-1}$ è autovalore di $A^{-1}$.
5. Se $\lambda$ è autovalore di $A$ con autovettore $x$, ${\lambda }^s$ è autovalore di $A^s$ con lo stesso autovettore ($s\in \mathbb{N}$)
6. Se $\lambda$ è autovalore di $A$, allora ${\lambda }^s$ è autovalore di $A^s$ con $s\in \mathbb{Z}$.

Molteplicità degli autovalori

Definiamo:

• Molteplicità algebrica: molteplicità come radice del polinomio caratteristico.
• Molteplicità geometrica: numero di vettori linearmente indipendenti associati ad esso.

Teorema dei cerchi di Gershgorin

Credits: Mu Prime Math - Learning Linear Algebra

Sia $A$ una matrice quadrata $n\times n$. Il teorema di Gerschgorin ci fornisce esattamente $n$ dischi nel piano complesso (dato che gli autovalori possono anche essere complessi) all'interno dei quali risiede un (vedasi note) autovalore. Viene creato un disco per ogni riga $i$ della matrice, centrato nell'elemento $a_{ii}$ e di raggio ${\rho }_i$ calcolato come segue:

\rho_i = \sum_{j=1, j \ne i}^{n} \abs{a_{ij}} \hspace{1cm} i=1, \dots, n

Possiamo definire ognuno dei dischi matematicamente come segue:

\gamma_i = \big{ z \in \C : \abs{z - {a_{ii}}} \le \rho_i \big} \hspace{1cm} i=1, \dots, n

Note In realtà il teorema non ci dimostra che in un cerchio risieda un solo autovalore, o che ce ne sia almeno uno. Una estensione del teorema dimostra che se il disco non interseca con altri dischi, allora contiene uno ed un solo autovalore.

Possiamo concludere che, definendo l'insieme dei dischi $\gamma ={\bigcup }_{i=1}^n{\gamma }_i$ allora: $$\lambda \in \sigma (A)⟹\lambda \in \gamma$$ Dove $\sigma (A)$ è lo spettro di $A$, ovvero l'insieme dei suoi autovalori. La dimostrazione è contenuta nel documento della professoressa ed è molto simile a quella indicata nel video.

Conseguenze del teorema

Una matrice strettamente diagonalmente dominante è non degenere.

Dimostrazione (per assurdo). Supponiamo che $A$, strettamente diagonalmente dominante, sia degenere, quindi $\det (A)=0$. Questo implica che $0\in \sigma (A)$ per i teoremi precedenti, e quindi per il teorema di Gershgorin $0$ è all'interno di uno dei dischi! Ovvero esiste il disco i-esimo tale che:

|0 - a_{ii}| \le \rho_i \Longrightarrow |a_{ii}| \le \sum_{j=1, j \ne i}^{n} \abs{a_{ij}}

Il che è assurdo per ipotesi di matrice strettamente diagonalmente dominante.

Teorema di Hermite

Se la matrice $A$ è hermitiana, ovvero $A=A^{\ast }$ (uguale alla sua trasposta coniugata) allora gli autovalori sono tutti reali. Il teorema vale anche per matrici simmetriche, quindi se $A=A^T$, a patto che sia anche definita positiva.

Teorema sulle matrici diagonalizzabili

Una matrice $A$ è diagonalizzabile se e solo se ha $n$ autovettori linearmente indipendenti.

Teorema di Shur

Sia A \in \C^{n\times n} allora esiste una matrice unitaria $U$ tale che $T=U^{\ast }AU$ dove $T$ è triangolare superiore. Se $A$ è reale, allora $U$ è ortogonale.

Matrice convergente

Una matrice $A$ si dice convergente se ${\lim }_{n\to \infty }{A}^s=0$ (matrice zero).

Altri teoremi

• Una matrice hermitiana è diagonalizzabile
• Una matrice è convergente se il raggio spettrale è inferiore a 1
• Il raggio spettrale è sempre minore o uguale alla norma della matrice
• Una matrice è convergente se la norma dell'esponenziazione della matrice tende a 0

3. Risoluzione di Sistemi Lineari

Definizione problema

Sia $A\in {\mathbb{R}}^{m\times n}$ e $b\in {\mathbb{R}}^m$ vogliamo trovare un vettore $x\in {\mathbb{R}}^n$ che soddisfi l'equazione $Ax=b$.

Tratteremo solo sistemi quadrati, in cui esiste un'unica soluzione $x$ se e solo se:

• esiste l'inversa $A^{-1}$ della matrice $A$
• oppure il rango della matrice $A$ è $n$
• oppure $Ax=0$ e in tal caso $x=0$

Teorema di Cramer

Se il determinante della matrice è non nullo, allora esiste unica una soluzione del sistema, ed è data da: $$x_i=\frac{\det ({\Delta }_i)}{\det (A)}$$ Dove ${\Delta }_i$ è $A$ a cui sostituiamo l'$i$-esima riga con il vettore dei termini noti.

Numero di condizionamento nei sistemi lineari

Sia $A$ una matrice quadrata, il numero di condizionamento $k(A)$ della matrice è un valore che misura la sensibilità della soluzione di $Ax=b$ alle perturbazioni nei dati. Idealmente vorremmo lavorare con matrici con numero di condizionamento basso. Supponendo l'esistenza dell'inversa, $k(A)$ si calcola come segue: $$k(A)=\Vert A\Vert \cdot \Vert A^{-1}\Vert$$

La matrice identità ha il più basso numero di condizionamento (1) ed è l'esempio perfetto di matrice ben condizionata nel caso della risoluzione di un sistema lineare associato.

In generale, in numero di condizionamento è correlato all'errore sui dati come segue: $$k(A)\ge \frac{\text{errore sul risultato}}{\text{errore sui dati}}$$

Problemi di calcolo del numero di condizionamento

La complessità del calcolo della matrice inversa cresce insieme alla crescita del numero di condizionamento. Ma per calcolare il numero di condizionamento è necessaria la matrice inversa! Affermiamo che:

1. Se $k(A)$ è vicino ad 1 allora $k(A)$ è facilmente calcolabile
2. Se $k(A)$ è grande, allora è difficilmente ricavabile

Correlazione tra numero di condizionamento e raggio spettrale

Supponendo l'esistenza dell'inversa, è presente la seguente relazione tra i due valori: $$k(A)\ge \rho (A)\rho (A^{-1})$$ Questo implica che: $$k(A)\ge \frac{\max \lambda }{\min \lambda }$$ Banalmente al numeratore è presente il raggio spettrale di $A$. Sappiamo che se esiste $\lambda$ autovalore di $A$, allora ${\lambda }^{-1}$ è autovalore di $A^{-1}$. Dato che invertiamo tutti gli autovalori, l'autovalore più piccolo di $A$, una volta invertito, diventa l'autovalore più grande di $A^{-1}$: $$\rho (A)\rho (A^{-1})=\max \lambda \cdot \frac{1}{\min \lambda }=\frac{\max \lambda }{\min \lambda }$$

Matrice di Hilbert

La matrice di Hilbert è un esempio di matrice mal condizionata per i sistemi lineari (ma è ben condizionata per altri problemi). In generale, l'elemento $h_{ij}=(i+j-1{)}^{-1}$, approfondire qui.

Metodi diretti

I metodi diretti trovano la soluzione del sistema lineare in un numero finito di passi. Sono adatti a sistemi con matrici piene, essendo che tendono a rimpire gli zeri della matrice.

Metodo delle sostituzioni in avanti

Supponiamo che la matrice $A$ del sistema sia triangolare inferiore. Il metodo delle sostituzioni in avanti risolve il sistema una soluzione per volta. Partendo da $x_1$ come segue: $$x_1=\frac{b_1}{l_{11}}$$ Si calcola il generico $x_i$ $$x_i=\frac{b_i-{\sum }_{j=1}^{i-1}{l}_{ij}{x}_j}{l_{ii}}$$

Metodo delle sostituzioni indietro

Supponiamo che la matrice $A$ del sistema sia triangolare superiore. Il metodo delle sostituzioni indietro risolve il sistema una soluzione per volta, partendo dal basso. Partendo da $x_n$ come segue: $$x_n=\frac{b_n}{u_{nn}}$$ Si calcola il generico $x_i$ $$x_i=\frac{b_i-{\sum }_{j=i+1}^n{u}_{ij}{x}_j}{u_{ii}}$$

Complessità metodi di sostituzione e considerazioni

I metodi di sostituzione presentati eseguono circa $\frac{n(n+1)}{2}$ divisioni e moltiplicazioni e $\frac{n(n-1)}{2}$ somme algebriche, per un costo computazionale $\Theta (n^2)$. Per derivare la complessità basta osservare che nella prima iterazione facciamo 1 div/mul, nella seconda 2, nella terza 3 e così via. Avendo $n$ righe totali basta calcolare la somma dei primi $n$ numeri naturali, che corrisponde proprio alla prima formula mostrata.

Questi metodi sono molto leggeri, ma necessitano di matrici triangolari. Gli altri metodi presentati cercheranno di triangolarizzare la matrice, per poi utilizzare i metodi di sostituzione per risolvere il sistema lineare. Se dell'istanza del problema cambia solo il vettore dei coefficienti, è possibile precalcolarsi alcuni dei parametri per rendere il metodo più rapido.

Metodo di eliminazione naive di Gauss (MEG)

Supponendo che la matrice $A$ sia non degenere, e sia $a_{11}\ne 0$ (altrimenti scambia righe), allora è possibile applicare il metodo di Gauss (a meno di un'altra condizione che enunceremo dopo). Alla prima iterazione si calcolano i moltiplicatori: $$m_{i1}^{(1)}=-\frac{a_{i1}}{a_{11}}i=2,\dots ,n$$ Aggiungiamo alla $i$-esima equazione la prima equazione moltiplicata per $m_{i1}$. Così facendo andremo ad annullare tutti gli elementi della prima colonna meno che il primo. In generale, durante la $i$-esima iterazione lo scopo è annullare tutti gli elementi della $i$-esima colonna al di sotto dell'$i$-esimo, quindi si calcolano i moltiplicatori tramite l'elemento $a_{ii}$, che prende il nome di pivot: $$m_{ji}^{(i)}=-\frac{a_{ji}}{a_{ii}}j=i+1,\dots ,n$$ E si aggiunge l'$i$-esima equazione alle successive moltiplicata per il rispettivo moltiplicatore. Al passo $n-1$ si ottiene un sistema triangolare superiore che può essere risolto con la sostituzione all'indietro. Il costo computazionale del metodo è circa $\frac{4}{3}{n}^3$.

Affinché il metodo funzioni è necessario che gli elementi della diagonale $a_{ii}$ siano non nulli ad ogni iterazione. Questo è garantito se tutti i minori principali di $A$ sono non nulli.

MEG con Pivot Parziale

La tecnica del pivot parziale evita le divisioni per zero o per numeri prossimi allo zero nel calcolo dei moltiplicatori. Alla $i$-esima iterazione si cerca la riga $i\le k\le n$ con l'elemento maggiore nella colonna $i$. $$a_{ki}=\underset{i\le s\le n}{\max }|a_{si}|$$ Si sostituisce la riga $k$-esima con la riga $i$-esima. Il beneficio sta nel fatto che viene minimizzata la grandezza del moltiplicatore, minimizzando anche l'errore amplificato durante le moltiplicazioni. La complessità totale del metodo è $O(n^2)$.

MEG con Pivot Totale

Il pivot totale ha gli stessi benefici del pivot parziale, ma amplificati: il metodo può scegliere tra tutti gli elementi della matrice (incompleta) e non solo quelli della colonna analizzata. Il drawback sta nella complessità implementativa e nel cambio di variabile necessario.

Supponiamo che il pivot corrente sia $a_{ii}$ e che l'elemento più grande della matrice sia $a_{rs}$. Allora prima si effettua una permutazione delle colonne $i\leftrightarrow s$ e si memorizza il cambio di variabile, e dopodiché si scambiano le righe $i\leftrightarrow r$. Quando si ottiene la soluzione, per ottenere la soluzione rispetto al problema originale bisogna applicare i cambi di variabile a ritroso alla soluzione.

Fattorizzazione LU

Sia $Ax=b$ il sistema lineare da risolvere, un metodo di fattorizzazione matriciale consiste nei seguenti passi:

1. Si trova una matrice $S$ non singolare tale che $SAx=Sb$ e $SA=U$ triangolare superiore.
2. Se $S$ è triangolare inferiore, lo sarà anche $S^{-1}$ e poniamo $L=S^{-1}$
3. Osserviamo che $A=LU=S^{-1}SA=A$ quindi $LU$ è una fattorizzazione di $A$

La fattorizzazione non dipende dai termini noti, quindi se nel sistema lineare in analisi variano solo i termini noti, precalcolando la fattorizzazione si ha un risparmio in efficienza. Attraverso la fattorizzazione LU riformuliamo il metodo di eliminazione di Gauss. Siano: $$A=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& & \ddots & \\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right]L_1=\left[\begin{array}{cccc}1& 0& \cdots & 0\\ m_{21}& 1& \cdots & 0\\ ⋮& & \ddots & \\ m_{n1}& 0& \cdots & 1\end{array}\right]$$ dove $m_{i1}$ per $i=2,\dots ,n$ sono i moltiplicatori mostrati in precedenza. Il prodotto $L_{1}A$ equivale al primo passo di Gauss. In generale, la matrice $L_i$ è definita come segue: $$L_i=\left[\begin{array}{ccccc}1& \dots & 0& \dots & 0\\ ⋮& & ⋮& & ⋮\\ 0& \dots & 1& \dots & 0\\ 0& \dots & m_{i+1,i}& \dots & 0\\ ⋮& & ⋮& & ⋮\\ 0& \dots & m_{n,i}& \dots & 1\end{array}\right]$$ Alla fine si ha: $$U=L_{n-1}{L}_{n-2}\dots L_2{L}_{1}A$$ Poniamo $L=L_{n-1}{L}_{n-2}\dots L_2{L}_1$, quindi $U=LA$. Poniamo $L=L^{-}1$ e osserviamo che: $$LU=L^{-1}U=L^{-1}LA=A$$ Adesso, la soluzione del sistema lineare $Ax=b$ si risolve in due passaggi:

1. $Ly=b$ e si risolve per $y$
2. $Ux=y$ e si risolve per $x$

Questo poiché con il passo (1) si trova $y=L^{-1}b$, mentre con il passo 2 si trova $x=U^{-1}{L}^{-1}b$, che è sicuramente soluzione del sistema $Ax=LUx=b$, ed entrambi possono essere risolti con sostituzioni in avanti ed indietro, essendo matrici triangolari inferiori (L) e superiori (U).

Non sempre esiste una fattorizzazione LU della matrice. Condizione necessaria di esistenza è che la matrice sia non degenere, quindi che abbia il determinante non nullo. Se vale tale condizione, allora esiste sicuramente una matrice di permutazione $P$ tale che $PA=LU$. L'esistenza è garantita per le matrici diagonalmente dominanti e simmetriche definite positive, senza dover permutare la matrice.

Matrici di permutazione

Le tecniche del pivot parziale e totale possono essere implementate nella forma matriciale del metodo di eliminazione di Gauss attraverso delle matrici di permutazione. Una matrice di permutazione è una matrice ottenuta scambiando righe o colonne della matrice identità. In particolare, scambiando la riga $i$ con la riga $j$ di $I$ e

• premoltiplicandola per $A$ si scambiano le righe
• postmoltiplicandola per $A$ si scambiano le colonne

Esempio: scambiamo la riga 1 e la riga 3 della matrice identità e vediamo cosa succede se la premoltiplichiamo ad una matrice qualunque: $$\left[\begin{array}{ccc}0& 0& 1\\ 0& 1& 0\\ 1& 0& 0\end{array}\right]\cdot \left[\begin{array}{ccc}4& 2& 7\\ 6& 2& 9\\ 7& 3& 5\end{array}\right]=\left[\begin{array}{ccc}7& 3& 5\\ 6& 2& 9\\ 4& 2& 7\end{array}\right]$$ L'effetto è stato uno scambio di righe. Se invece postmoltiplichiamo la matrice di permutazione: $$\left[\begin{array}{ccc}4& 2& 7\\ 6& 2& 9\\ 7& 3& 5\end{array}\right]\cdot \left[\begin{array}{ccc}0& 0& 1\\ 0& 1& 0\\ 1& 0& 0\end{array}\right]=\left[\begin{array}{ccc}7& 2& 4\\ 9& 2& 6\\ 5& 3& 7\end{array}\right]$$ L'effetto è stato uno scambio di colonne.

Metodo di Doolittle e di Croud

Se esplicitassimo il sistema lineare $LU=A$ avremmo $n^2$ equazioni (una per ogni elemento di $A$) ed $n^2+n$ incognite (elementi di $L$ ed $U$ supponendo che siano rispettivamente triangolari inferiore e superiore), questo implica che abbiamo $n$ gradi di libertà, ovvero la fattorizzazione $LU$ non è unica.

I seguenti metodi eliminano i gradi di liberta imponendo dei vincoli:

• Doolittle: si fissa $l_{ii}=1$ in $L$ (equivalente a eliminazione gaussiana senza pivot)
• Crout: si fissa $u_{ii}=1$ in $U$

Effetto Fill-in

I metodi di fattorizzazione modificano la matrice iniziale causando un effetto fill-in, ovvero gli zeri diventano elementi non nulli. Se la matrice è inizialmente sparsa, conviene optare per metodi iterativi.

Metodo di Cholesky

Sia $A\in {\mathbb{R}}^{n\times n}$ una matrice simmetrica ($A=A^T$) e definita positiva, allora esiste almeno una matrice $L$ triangolare inferiore tale che $A=L{L}^T$. Inoltre, se si impone che $l_{ii}>0$ la fattorizzazione è unica.

Dimostrazione (costruttiva).

Per il criterio di Sylvester, la matrice ha determinante strettamente positivo, e quindi esiste una fattorizzazione LU. Proviamo a costruire la fattorizzazione $U^{T}U$ come segue: $$U^{T}U=\left[\begin{array}{ccc}{u}_{11}& & 0\\ ⋮& \ddots & \\ u_{1n}& \dots & u_{nn}\end{array}\right]\left[\begin{array}{ccc}{u}_{11}& \dots & u_{1n}\\ ⋮& \ddots & \\ 0& \dots & u_{nn}\end{array}\right]=\left[\begin{array}{ccc}{a}_{11}& \dots & a_{1n}\\ ⋮& \ddots & \\ a_{n1}& \dots & a_{nn}\end{array}\right]=A$$ Dato che le incognite sono gli elementi di $U$, che sono $n(n+1)/2$, possiamo strutturare $n^2$ equazioni (non lineari) tramite il prodotto riga-colonna.

Esplicitiamo il calcolo dell'elemento $a_{kj}$ supponendo che $k>j$ (le moltiplicazioni per indici maggiori di $j$ sono nulle per la triangolarità delle matrici), quindi abbiamo: $$a_{kj}=\sum _{i=1}^j{u}_{ki}{u}_{ij}=u_{kj}{u}_{jj}+\sum _{i=1}^{j-1}{u}_{ki}{u}_{ij}$$ Estraiamo l'elemento $u_{kj}$ dall'equazione: $$u_{kj}=\frac{a_{kj}-{\sum }_{i=1}^{j-1}{u}_{ki}{u}_{ij}}{u_{jj}}$$ Vediamo invece come calcolare un elemento della diagonale $u_{kk}$. Tramite la seguente equazione cerchiamo di calcolare l'elemento $a_{kk}$: $$a_{kk}=\sum _{p=1}^k{u}_{pk}^2=u_{kk}^2+\sum _{p=1}^{k-1}{u}_{pk}^2$$ Oss. la sommatoria termina a $k$ e non ad $n$ perché per la struttura triangolare della matrice gli elementi $a_{k,k+1},\dots ,a_{k,n}$ sono nulli. Ora possiamo estrapolare l'elemento $u_{kk}$ come: $$u_{kk}=\sqrt{a_{kk}-\sum _{p=1}^{k-1}{u}_{pk}^2}$$ Calcolando $u_{11}$ la sommatoria verrebbe annullata, quindi $u_{11}=\sqrt{a_{11}}$. Bisogna seguire un certo andamento per poter calcolare tutti gli elementi di $U$, dato che un calcolo potrebbe richiedere dei calcoli di altri elementi della matrice. Nella seguente foto è illustrato tale andamento su una matrice $3\times 3$:

Algoritmo di Thomas

L'algoritmo di Thomas è un TDMA (tridiagonal matrix algorithm), ovvero un algoritmo che opera su matrici tridiagonali. Una matrice tridiagonale di dimensione $n\times n$ ha $3n-2$ elementi. Scriviamola come prodotto di due matrici particolari, in cui le incognite sono ${\alpha }_i$ per $i=1,\dots ,n$ e ${\gamma }_i$ per $i=1,\dots ,n-1$.

I coefficienti sono determinabili come segue:

• ${\alpha }_1=a_1$ e ${\alpha }_i=a_i-b_i{\gamma }_{i-1}$ per $i=2,\dots ,n$
• ${\gamma }_1=c_1/{\alpha }_1$ e ${\alpha }_i=c_i/{\alpha }_i$ per $i=2,\dots ,n-1$

Per una trattazione completa consultare il seguente documento.

Metodi iterativi

Un metodo iterativo produce una successione di soluzioni che, sotto opportune condizioni, converge alla soluzione reale. Comunemente si parte da una soluzione randomica $x^{(0)}$, in quanto si dimostra che il metodo converge comunque.

Descrizione generale di un algoritmo iterativo

Sia $A\in {\mathbb{R}}^{n\times n}$ una matrice non degenere, si impone che: $$\{\begin{array}{l}Ax=b\\ A=M-N\end{array}$$ Quindi si ha: $$\begin{array}{rl}Ax& =b\\ (M-N)x& =b\\ Mx-Nx& =b\\ M{x}^{(k+1)}& =N{x}^{(k)}+b\end{array}$$ E quindi si determina la soluzione al passo $k+1$ dalla soluzione al passo precedente: $$x^{(k+1)}=M^{-1}N{x}^{(k)}+M^{-1}b$$ L'equazione $A=M-N$ prende il nome di decomposizione, e si dice regolare se $M$ non è degenere e gli elementi dell'inversa $M^{-1}$ e di $N$ sono tutti maggiori di zero. Spesso la matrice $M^{-1}N$ viene chiamata "matrice del metodo" a cui si fa riferimento, poiché studiandola è possibile studiarne la convergenza.

Criteri d'arresto

Fissata una tolleranza $ϵ$, allora un criterio d'arresto standard utilizzabile per fermare un algoritmo iterativo è il seguente: $$\frac{\Vert x^{(k+1)}-x^{(k)}{\Vert }_{\infty }}{\Vert x^{(k)}{\Vert }_{\infty }}$$

Metodo iterativo convergente

Un metodo iterativo si dice convergente se per qualunque vettore iniziale $x^{(0)}$ il metodo sarà convergente.

Teorema sulla convergenza (1)

Sia $A=M-N$ una decomposizione regolare di $A$ e sia $\Vert M^{-1}N\Vert \le \lambda <1$. Allora:

1. $A$ non è singolare
2. Il metodo iterativo associato alla decomposizione è convergente
3. $\Vert x^{(k)}-x\Vert$ ci dà un limite all'errore commesso

Teorema sulla convergenza (2)

Condizione necessaria e sufficiente affinché il metodo iterativo sia convergente è che il raggio spettrale $\rho$ della matrice $M^{-1}N$ sia minore di 1: $$\rho (M^{-1}N)<1$$

Corollario sulla convergenza (3)

Dato che il determinante di una matrice è anche definito come il prodotto dei suoi autovalori, allora condizione necessaria affinché il metodo converga è che $$|\det (M^{-1}N)|<1$$ Se il determinante è maggiore di 1 allora esiste almeno un autovalore $\ge 1$, e quindi il metodo non può convergere per il teorema precedente.

Corollario sulla convergenza (4)

La traccia di una matrice è definita come la somma dei suoi autovalori. Condizione necessaria affinché il metodo converga è che: $$|t_r(M^{-1}N)|<n$$ Se vale il contrario, allora almeno uno degli autovalori è $\ge 1$, quindi il metodo non può convergere per il teorema precedente.

Correlazione tra raggio spettrale e convergenza

Il raggio spettrale $\rho$ relativo alla matrice di uno specifico metodo è correlato alla convergenza dello stesso come segue:

• Se $\rho \ge 1$ il metodo diverge
• Se $\rho <1$ il metodo converge
• Tanto più $\rho$ si avvicina ad 1, tanto più il metodo convergerà lentamente

Velocità di convergenza di un metodo iterativo

Tramite la seguente valutazione è possibile ricavare il numero di iterazioni $k$ affinché l'errore $e$ generato da un certo metodo scenda al di sotto di un errore assegnato $m$. Siano $e^{(k)}$ ed $x^{(k)}$ rispettivamente l'errore e la soluzione generati al passo $k$. L'errore è ricavato come segue: $$e^{(k+1)}=x^{(k+1)}-x^{(k)}$$ Che deriva dalla formula di approssimazione dell'errore, dove l'$x$ soluzione reale viene approssimato con la soluzione calcolata allo step successivo. È possibile calcolare un lower bound del numero di iterazioni $k$ necessarie a far scendere l'errore sotto $m$, come segue: $$k\ge -\frac{m}{{\log }_{10}\rho (B)}$$ Dove $B$ è la matrice del metodo. Sia $R$ il fattore di convergenza del metodo iterativo, definito come segue: $$R={\log }_{10}\rho (B)$$ Allora si avrà: $$k\ge -\frac{m}{R}$$ Più grande è $R$, più la convergenza sarà veloce. Da questo deriva che più $\rho (B)$ è vicino ad 1, più l'intero denominatore si annulla e più le esecuzioni necessarie crescono. Se $0<\rho (B)<1$ il logaritmo assumerà un valore negativo che è più grande tanto più è vicino allo zero.

Metodo di Jacobi

Partendo da una soluzione iniziale $x^{(0)}$, la soluzione $x^{(k)}$ si determina dalla soluzione $x^{(k-1)}$ come segue:

Condizione sufficiente di convergenza per Jacobi

Sia $Ax=b$ il sistema da risolvere, condizione sufficiente affinché Jacobi converga è che $A$ sia strettamente diagonalmente dominante.

Formulazione matriciale

La decomposizione presa in considerazione nel metodo è la seguente: $$A=D-E-F$$ Se si pone $M=D$ ed $N=E+F$ allora si ottiene la decomposizione classica: $$A=M-N$$ Le tre matrici rappresentano:

• $D$ la diagonale principale di $A$
• $E$ gli elementi al di sotto della diagonale principale di $A$, cambiati di segno.
• $F$ gli elementi al di sopra della diagonale principale di $A$, cambiati di segno.

Dalle matrici, ritroviamo l'equazione di aggiornamento del metodo di Jacobi: $$\begin{array}{rl}Ax& =b\\ (D-E-F)x& =b\\ Dx-(E+F)x& =b\\ Dx& =(E+F)x+b\\ x& =D^{-1}(E+F)x+D^{-1}b\end{array}$$ Applichiamo gli indici all'ultimo passaggio per indicare l'aggiornamento: $$x^{(k+1)}=D^{-1}(E+F)x^{(k)}+D^{-1}b$$

Ricordiamo che l'inversa di una matrice diagonale è pari al reciproco dei singoli elementi non nulli.

La matrice del metodo è $D^{-1}(E+F)$ e studiandola è possibile studiare la convergenza.

Metodo di Gauss-Seidel (GS)

È una variante del metodo di Jacobi che sfrutta il calcolo delle componenti precedenti della soluzione corrente per calcolare le successive. Si vede empiricamente che il metodo converge spesso più velocemente rispetto a quello di Jacobi (quindi il raggio spettrale di GS è minore). Nell'esempio sottostante, gli apici rossi sono quelli di Gauss-Seidel, mentre quelli neri sono di Jacobi.

Condizione sufficiente di convergenza per GS

Sia $Ax=b$ il sistema da risolvere, condizione sufficiente affinché Jacobi converga è che $A$ sia strettamente diagonalmente dominante o che $A$ sia simmetrica definita positiva.

Formulazione matriciale

La formulazione è analoga a quella di Jacobi, ma viene isolata la matrice $(D-E)$ anziché isolare solo la matrice $D$. Se volessimo affibiare le matrici alla trattazione generale, avremmo che $M=(D-E)$ e $N=F$. Vediamo i passaggi: $$\begin{array}{rl}Ax& =b\\ (D-E-F)x& =b\\ (D-E)x-Fx& =b\\ (D-E)x& =Fx+b\\ x& =(D-E{)}^{-1}Fx+(D-E{)}^{-1}b\end{array}$$ Supponendo l'esistenza dell'inversa. Applichiamo gli indici all'ultimo passaggio per indicare l'aggiornamento: $$x^{(k+1)}=(D-E{)}^{-1}F{x}^{(k)}+(D-E{)}^{-1}b$$ La matrice del metodo in questo caso è $(D-E{)}^{-1}F$ e studiandola è possibile studiare la convergenza.

Metodo Successive Over-Relaxation (SOR)

SOR è un metodo basato su GS che applica un rilassamento al calcolo della soluzione. Si fissa un certo parametro $\omega$ e si calcola la soluzione utilizzata come la media pesata tra la soluzione di GS e la soluzione precedente: $$x^{(k+1)}=\omega x^{(k+1)}+(1-\omega )x^{(k)}$$

Formulazione matriciale

Ricaviamo l'equazione di aggiornamento del metodo, partendo dal sistema lineare: $$\begin{array}{rl}Ax& =b\\ \omega Ax& =\omega b\\ \omega (D-E-F)x& =\omega b\\ Dx+\omega (D-E-F)x& =Dx+\omega b\\ Dx+\omega Dx-\omega Ex-\omega Fx& =Dx+\omega b\\ (D-\omega E)x& =\omega Fx+Dx-\omega Dx+\omega b\\ (D-\omega E)x& =[\omega F+(1-\omega )D]x+\omega b\\ x& =(D-\omega E{)}^{-1}[\omega F+(1-\omega )D]x+(D-\omega E{)}^{-1}\omega b\end{array}$$ Introduciamo la matrice $H(w)=(D-\omega E{)}^{-1}[\omega F+(1-\omega )D]$, sostituiamo nell'ultima equazione, e introduciamo gli apici del metodo iterativo: $$x^{(k+1)}=H(\omega )x^{(k)}+(D-\omega E{)}^{-1}\omega b$$ Poniamo: $$\begin{gathered} L=D^{-1}E \\ R=D^{-1}F \end{gathered}$$ E manipoliamo la matrice $H(\omega )$ come segue: $$\begin{array}{rl}H(\omega )& =(D-\omega E{)}^{-1}[(1-\omega )D+\omega F]\\ & =[D(I-\omega D^{-1}E){]}^{-1}D[(1-\omega )I+\omega D^{-1}F]\\ & =[D(I-\omega L){]}^{-1}D[(1-\omega )I+\omega R]\\ & =(I-\omega L{)}^{-1}{D}^{-1}D[(1-\omega )I+\omega R]\\ & =(I-\omega L{)}^{-1}[(1-\omega )I+\omega R]\end{array}$$ La matrice finale è data da: $$H(\omega )=(I-\omega L{)}^{-1}[(1-\omega )I+\omega R]$$ E studiando $\rho (H(\omega ))$ non solo è possibile studiare la convergenza del metodo, ma dato il parametro $\omega$ arbitrario sarà possibile aggiustare la matrice per manipolare la sua convergenza.

Teorema di condizionamento del parametro di SOR

Si ha la seguente relazione sempre verificata $$\rho (H(\omega ))\ge |\omega -1|\forall \omega \in \mathbb{R}$$ Da cui osserviamo che:

\begin{cases} \text{SOR diverge} & \omega \le 0 \or \omega \ge 2 \ \text{SOR converge} & 0 < \omega < 2 \end{cases}

Fissato $\omega$ nel range convergente, è possibile muoversi all'interno del range di convergenza per tarare l'algoritmo e renderlo più rapido.

Variante SSOR

Symmetric SOR (SSOR) è una variante di SOR in cui viene preservata la possibile simmetria della matrice.

Metodo di discesa del gradiente

Quando un sistema $Ax=b$ è molto grande, è possibile risolverlo riconducendolo ad un problema di minimizzazione. Si considera la forma quadratica del problema: $$\varphi (x)=\frac{1}{2}{x}^{T}Ax+x^{T}b$$ Si calcola il gradiente della funzione: $$\Delta \varphi =\left(\frac{\partial \varphi }{\partial x_1},\dots ,\frac{\partial \varphi }{\partial x_n}\right)$$ Utilizziamo il l'opposto del gradiente come direzione di discesa: $$d^{(k)}=-\Delta \varphi (x^{(k)})$$ Sia $a_k$ la lunghezza del passo di aggiornamento all'$k$-esimo passo, aggiorniamo la soluzione: $$x^{(k+1)}=x^{(k)}+a_k{d}^{(k)}$$

4. Interpolazione

Definizione del problema

È un problema di approssimazione di una funzione o di un insieme di dati con una funzione che sia più semplice e che abbia buone proprietà di regolarità. È utilizzata quando i dati sono forniti con precisione. La condizione di interpolazione è che la funzione interpolante e la funzione interpolata coincidano nei nodi che sono le ascisse dei dati.

Interpolazione polinomiale

Nell'interpolazione polinomiale la funzione interpolata sarà un polinomio. Con $P_n$ indicheremo l'insieme dei polinomi di grado minore o uguale ad $n$.

Teorema di Vandermonde

Siano dato gli $n+1$ nodi $x_0,\dots ,x_n$ distinti ed i dati in corrispondenza $y_0,\dots ,y_n$, allora esiste sempre uno ed un solo polinomio interpolante di grado $n$ passante per tutti i punti $x_i$. $${\exists }_{1}p\in P_n:p(x_i)=y_{i}i=0,1,\dots ,n$$

Metodo dei coefficienti indeterminati

Anche detto metodo di Vandermonde, consiste nel ricavare il polinomio risolvendo un sistema lineare. Un polinomio $p\in P_n$ ha $n+1$ coefficienti, che sarebbero le incognite. Se abbiamo $n+1$ nodi allora possiamo facilmente costruire un sistema lineare quadrato per ottenere i coefficienti del polinomioì, e quindi un polinomio che rispetti la condizione di interpolazione. $$\{\begin{array}{l}{a}_0+a_1{x}_0^1+\cdots +a_n{x}_0^n=y_0\\ a_0+a_1{x}_1^1+\cdots +a_n{x}_1^n=y_1\\ ⋮\\ a_0+a_1{x}_n^1+\cdots +a_n{x}_n^n=y_n\end{array}$$ Purtroppo il metodo genera una matrice mal condizionata per il problema dei sistemi lineari al crescere di $n$. Esistono metodi migliori per il problema dell'interpolazione.

Metodo dei polinomi di Lagrange

Si pone il polinomio $p\in P_n$ alla seguente espressione: $$p(x)=\sum _{i=0}^n{y}_i{L}_i(x)$$ Dove gli $L_i$ sono i polinomi di Lagrange e sono definibili come una base dello spazio dei polinomi $P_n$, mentre nei casi precedenti viene usata la base canonica ($1,x,x^2,\dots ,x^n$). La costruzione avviene come segue: $$L_i(x)=\prod _{j\ne i,j=0}^n\frac{x-x_j}{x_i-x_j}$$ Il comportamento di $L_i$ con i nodi è quello di un interruttore: $$Li(x_j)=\{\begin{array}{ll}0& x\ne x_i\\ 1& x=x_i\end{array}$$ Quindi è ovvio che il polinomio definito precedentemente interpola i nodi, essendo che quando è sottomesso a $x_i$ si attiva $L_i(x_i)=1$ e si annullano tutti gli altri, quindi $p(x_i)=y_i{L}_i(x_i)=y_i$.

Metodo delle differenze divise di Newton

Anche in questo metodo si utilizza una base differente, e si crea $p\in P_n$ come segue: $$p(x)=b_0+b_1(x-x_0)+b_2(x-x_0)(x-x_1)+\cdots +b_n\prod _{i=0}^{n-1}(x-x_i)$$ Osserviamo che:

1. $p(x_0)=b_0$
2. $p(x_1)=b_0+b_1(x-x_0)$
3. etc.