Numerical Analysis: 2017

Costi computazionali

Formule sommatorie
$\sum\limits_{k=1}^nk^2 = \frac{n(n+1)(n+2)}{6}$
$\sum\limits_{k=1}^nk = \frac{n(n+1)}{2}$

Laplace: $O(n!)$
Cramer: $O((n+1)!)$
Risoluzione sistema $Tx=b$: $O(\frac{n^2}{2})$, $T$ è triangolare
Confronti nel Pivoting parziale: $O(\frac{n^2}{2})$
Confronti nel pivoting totale: $O(\frac{n^3}{3})$
Tecniche compatte per determinare la fattorizzazione LU: $O(\frac{n^3}{3})$
Risolvere un sistema AX=B con fattorizzazione LU: $O(\frac{n^3}{3}+rn^2)$
Metodo di Gauss (triangolarizzazione): $O(\frac{n^3}{3})$
Diagonalizzare una matrice con la Variante di Gauss Jordan del metodo di eliminazione di Gauss: $O(\frac{n^3}{2})$
Metodo di Gauss-Jordan per il calcolo dell'inversa: $O(n^3)$
Metodi iterativi per sistemi lineari: $O(Kn^2)$, dove $K$ è il numero di iterazioni
Costruzionee tabella delle differenze divise: $O(\frac{n^2}{2})$
Polinomio di Lagrange: $O(2n^2+2n)$

Interpolazione di Hermite

E' un caso particolare di polinomio osculatore.
Dati $n+1$ nodi $x_0, x_1, ...,x_n$, il polinomio interpolante di Hermite di grado $2n+1$ è tale che soddisfa:
$p_H=f(x_i), i=0,...,n$ ($n+1$ condizioni)
$p^{'}_H=f^{'}(x_i), i=0,...,n$ ($n+1$ condizioni)
Abbiamo $2n+2$ condizioni.

Il polinomio di Hermite può essere rappresentato usando gli $L_j(x)$ di Lagrange:
$p_H(x)=\sum\limits_{j=0}^n[U_j(x)f(x_j)+V_j(x)f^{'}(x)]$ dove
per $j =0,...,n$

$U_j(x)=L^2_j(x)[1-2L^{'}_j(x_j)(x-x_j)]$
$V_j(x)=L^2_j(x)(x-x_j)$

Resto nell'interpolazione di Hermite

Siano $a=\min\limits_{0 \leq i \leq n}\{x_i\}, b=\max\limits_{0\leq i\leq n}\{x_i\}$ e sia $f(x) \in C^{2n+2}([a,b])$.
Allora $\exists\xi\in]a,b[$ tale che il resto del polinomio di Hermite è dato da $r(x)=\pi_{n+1}^2(x)\frac{f^{2n+2}(\xi)}{(2n+2)!}$

Nodi di Chebyshev

Osservando i polinomi interpolanti la funzione di Runge su nodi equidistanti si vede che l'errore di interpolazione è minore al centro dell'intervallo [-5,5], mentre cresce, al crescere di n,agli estremi dell'intervallo [-5, 5]. Ci suggerisce di interpolare con una distribuzione di nodi non equidistanti. Una buona scelta è quella dei nodi di Chebyshev:
$x_i=\frac{a+b}{2}-\frac{b-a}{2}cos[\frac{(2i+1)}{2(n+1)}\pi], i=0,1,...,n$

Questa scelta è quella che nell'espressione del resto $r(x)$, minimizza il $\pi_{n+1}(x)$ .
Per concludere, la "bontà" dell'interpolazione, dipende da:

il numero dei nodi (non ha senso costruire polinomi interpolatori per $n \geq 7$)
la distribuzione dei nodi

Polinomi osculatori

// $f^{(0)}(x_i) = {p^{(0)}}_n(x_i), i=0,...,n$

I polinomi osculatori sono polinomi che nei nodi $x_i, i=0,...,n$ soddisfano condizioni che coinvolgono la funzione interpolanda $f(x)$ e le sue derivate

$p^{(k)}(x_i)=f^{(k)}(x_i), i=0,...n, k=0,...n$

Prendere n ci garantisce che $\exists!$ il polinomio osculatore.

Il grado del polinomio osculatore è pari al numero delle condizioni di interpolazione meno uno.

Condizione sufficiente per l'$\exists!$ del polinomio osculatore

Il polinomio osculatore $\exists!$ quando solo se, quando si impone una condizione di interpolazione sulla $k$-esima derivata in un nodo $x_j$ fissato, $p^{(k)}(x_j)=f^{(k)}(x_j)$ allora vengono imposte anche tutte le condizioni di interpolazione sule derivata di ordine inferiore a $k$, cioè $0,1,...,k-1$ sullo stesso nodo $x_j$.

Metodo di Gauss Seidel

$M = D-B$
$N = C$
$B_{GS} = M^{-1}N=(D-B)^{-1}C$
$d=M^{-1}b=(D-B)^{-1}b$

Poiché in generale, non è sempre agevole e può risultare computazionalmente oneroso valutare l'inversa della matrice $(D-B)$ per valutare la matrice di iterazione $B_{GS}=(D-B)^{-1}C$ conviene considerare la seguente forma equivalente:

$x^{(k)} = [(D-B)^{-1}C]x^{(k-1)}+b$
$Dx^{(k)}-Bx^{(k)}=Cx^{(k-1)}+b$
$Dx^{(k)}=Bx^{(k)}+Cx^{(k-1)}+b$
$x^{(k)}=D^{-1}Bx^{(k)}+D^{-1}Cx^{(k-1)}+D^{-1}b$ (molto simile a Jacobi)

Tale espressione può essere interpretata come l'iterazione di Jacobi nella quale allo step $k$-esimo, per calcolare la componente $x_i^{(k)}$ si usano le componenti già calcolate nello stesso step, $x_j^{(k)}, j=1,...,i-1$ e le componenti $x_j^{(k-1)}, j=i+1,...,n$ dello step precedente cioè il $(k-1)$-simo. Per tale motivo il metodo di Gauss-Seidel è detto "Metodo degli spostamenti successivi".

[... mancano i calcoli per ricavarsi le componenti ]

Quindi, in termini di componenti:

$x_i^{(k)} = \frac{1}{a_{ii}}[b_i-(\sum\limits_{j=1}^{i-1}a_{ij}x_j^{(k)}+\sum\limits_{j=i+1}^{n}a_{ij}x_j^{(k-1)})]$ per $i=1,2,...,n$ e $k=1,..,n$

Metodo di Cholesky

Se A è simmetrica e definita positiva esiste ed è unica la fattorizzazione $A=LL^T$ con $L$ matrice triangolare inferiore.
La tecnica compatta che va sotto il nome di Metodo di Cholesky, usa le seguenti relazioni per determinare gli elementi $l_{ij}$ della matrice $L$.

$j=1,2,...,n$ $\begin{cases}l_{jj}=\sqrt{a_{jj}-\sum\limits_{k=1}^{j-1}l_{jk}^2}, \\l_{ij}=\frac{1}{l_{jj}}[a_{ij}-\sum\limits_{k=1}^{j-1}l_{ik}l_{jk}], & \mbox{if } i=j+1,...,n, j \neq n\end{cases}$

Resto nell'interpolazione

Viste le condizioni di interpolazione $f(x_i)=p(x_i), i=0,...,n$, per "costruzione" (e prescindendo dagli errori sui dati e dagli errori di arrotondamento) l'operazione di interpolazione restituisce una risposta "esatta" sui nodi per cui se la funzione $f(x)$

Insieme dei numeri macchina

Il calcolatore, essendo una macchina finita, consente la rappresentazione di un sottoinsieme finito di reali.

Definizione.

Si definisce insieme dei numeri macchina con $t>0$ cifre significative, base $\beta \geq 2$ e rango $(1,U)$ il sottoinsieme di $\mathbb{R}$ così definito:
$\mathbb{F}(t, \beta, L, U)=\{0\} \cup \{ x \in \mathbb{R} / x = segno(x) \beta^p \sum \limits_{i=1}^t d_i \beta^{-1} \}$ con $0 \leq d_i \leq \beta -1, i=1,2...,t$
$d_1 \neq 0$
$L \leq p \leq U$ dove $L$ intero negativo, $U$ intero positivo.
$x=\pm(.x_1x_2...x_t)\beta^p$

Il calcolatore usa solo un sottoinsieme finito di $\mathbb{R}$ e ogni volta che deve rappresentare un numero $x \in \mathbb{R}$ tale che $x \notin \mathbb{F}$ deve associargli un opportuno $\hat{x}\in \mathbb{F}$ tale che $x \approx \hat{x}$. Se $x \in \mathbb{R}, x\neq0, x\notin \mathbb{F}(t, \beta, L, U)$ si procede nel modo seguente:

se $p<L \Rightarrow x \rightarrow \hat{x}=0$ e indicazione di underflow
se $p>U \Rightarrow x \rightarrow$ nessun $\hat{x}$ ("Inf") e indicazione di overflow
se $L \leq p \leq U$ ma $x \notin \mathbb{F}$ perché le sue cifre con $i>t$ non sono tutte nulle, si procede per troncamento o per arrotondamento:
- se $d_{t+1} < \frac{\beta}{2}$
  $\Rightarrow x \rightarrow \hat{x}=trn(x)=\beta^p\sum\limits_{i=1}^td_i\beta^{-1}$
  $\Rightarrow \hat{x}=(.d_1d_2...d_t)\beta^p$
- se $d_{t+1} \geq \frac{\beta}{2}$
  $\Rightarrow x \rightarrow \hat{x}=arr(x)=trn(x)\beta^{p-t}$
  $\Rightarrow \hat{x}=(.d_1d_2...(d_t)+1)\beta^p$

Rappresentazione in macchina

Un calcolatore può rappresentare solo un numero finito di cifre, ciò implica che i numeri reali introdotti nel calcolatore vengono approssimati così come vengono approssimati i risultati delle operazioni elementari su tali numeri.
Quindi quando eseguiamo un algoritmo su un calcolatore, si ha una generazione e propagazione di errori.
Per esaminare la generazione degli errori inerenti bisognare parlare di rappresentazione in base di un numero. Dato un intero $\beta>1$ (detto base della rappresentazione), un nnumero $x \in \mathbb{R}, x \neq 0$, si può scrivere/esprimere in modo univoco nella forma:
$x=segno(x)[d_1\beta^{-1}+d_2\beta^{-2}+...]\beta^p=\pm\beta^p\sum\limits_{i=1}^{+\infty}d_i\beta^{-1}$

$x=segno(x)=\begin{cases}1, & \mbox{if } x>0\\-1, & \mbox{if } x<0\end{cases}$
$p \in \mathbb{Z}$ detto caratteristica (o esponente) della rappresentazione di $x$.
$d_i \in \mathbb{Z}, i=1,2,...$, sono dette cifre della rappresentazione di $x$ e verificano le condizioni

$0 \leq d_i \leq \beta - 1, i=1,2,...$
$d_1 \neq 0$ e $d_i$ non sono mai definitivamente uguali a $\beta-1$

Notazione posizionale.

$x=\pm(.d_1d_2...)\beta^p$ $d_1 \neq 0 \rightarrow$ notazione normalizzata.

La normalizzazione oltre ad essere necessaria per l'unicità, si rivela vantaggiosa quando si rappresenta un numero solo con un numero finito di cifre.

Esempio:
$x=frac{1}{7000}$ può essere rappresentato in $\beta=10$ nei due modi seguenti:

$x_1=(.142857142857...)10^{-3}$ ovvero in rappresentazione normalizzata
$x_2=(.000142857142857...)$ rappresentazione non normalizzata.

per un numero massimo di $t=12$ cifre di mantissa

$x_1=(.142857142857...)10^{-3}$
$x_2=(.000142857142...)$

$\Rightarrow x_1$ è una rappresentazione di $x$ migliore di $x_2$ perché guadagno 3 numeri di mantissa. Ovvero fornisce un'approssimazione migliore di $x$ rispetto ad $x_2$, perché l'informazione che in $x_2$ è contenuta negli zeri viene spostata all'esponente.

Condizionamento di un problema

Stima degli errori sui dati input: quando si realizza una computazione per via automatica, qualunque metodo si usi, non si può prescindere dall'errore di rappresentazione dei dati (errore inerente) e dall'errore di arrotondamento sulle singole operazioni (errore algoritmico), che dipende dal metodo usato ed è indipendente dall'errore sui dati.
L'errore inerente da una stima della sensibilità della soluzione agli errori sui dati.

Definizione.

Un problema si dice malcondizionato se a "piccole" perturbazioni/errori sui dati iniziali corrispondono "grandi" perturbazioni sui dati finali. Viceversa si dice ben condizionato. Chiaramente i termini "piccole" e "grandi" dipendono dal contesto.

Matrici definite in segno

Sia $A \in \mathbb{R}^{nxn}$ una matrice simmetrica ($A=A^T, a_{ij}=a_{ji}$) e $x \in \mathbb{R}^n$, allora $\alpha = x^TAx$
$\alpha = x^TAx$ è uno scalare, cioè $\alpha \in \mathbb{R}$.

Definizione.

Se $A \in \mathbb{R}^{nxn}$ è una matrice simmetrica e $\forall x \in \mathbb{R}^n, x \neq 0$ il numero reale $\alpha = x^TAx$ mantiene lo stesso segno, la matrice $A$ si dice definita in segno.
In particolare:

se $\alpha = x^TAx > 0$, $A$ è detta definita positiva
se $\alpha = x^TAx \geq 0$, $A$ è detta semidefinita positiva
se $\alpha = x^TAx < 0$, $A$ è detta definita negativa
se $\alpha = x^TAx \leq 0$, $A$ è detta semidefinita negativa

Matrice a predominanza diagonale

Definizione.

Una matrice $A \in \mathbb{R}^{nxn}$ si dice:

a predominanza diagonale in senso stretto (o forte) per righe se:
$|a_{ii}|>\sum\limits_{j=1, i \neq j}^n|a_{ij}|, \forall i=1,...,n$
a predominanza diagonale debole per righe se:
$|a_{ii}|\geq \sum\limits_{j=1, i \neq j}^n|a_{ij}|, \forall i=1,...,n$
a predominanza diagonale in senso stretto per colonne se:
$|a_{jj}|> \sum\limits_{i=1, i \neq j}^n|a_{ij}|, \forall j=1,...,n$
a predominanza diagonale debole per colonne se:

Proprietà delle matrici a dominanza diagonale

Se $A$ è simmetrica $\Rightarrow$ pdf per righe $\equiv$ pdf per colonne
Se $A$ è a pdf per colonne $\Rightarrow$ nel procedimento di eliminazione di Gauss non è necessario fare pivoting parziale.
L'algoritmo di Gauss è sempre stabile.
Se $A$ è a pdf:

$A$ non è singolare, cioè $|A| \neq 0$
Tutti i minori principali di testa di $A$ sono non nulli
$A$ è fattorizzabile in forma LU

Fattorizzazioni di Matrici "Speciali"

Pro e contro della strategia di pivoting

pro:

rende più stabile l'algoritmo di fattorizzazione

contro:

modificare la struttura della matrice (ad esempio si può avere perdita di simmetria: una matrice simmetrica occupa memoria solo $\frac{n^2}{2}$ posizioni anziché $n^2$.

Ci chiediamo se esistono classi di matrici per cui la condizione di pivoting è automaticamente soddisfatta.

Sì, cioè:

Matrici a diagonale dominante (o "a predominanza diagonale",..)
Matrici simmetriche definite in segno

Metodi iterativi per sistemi lineari

Prerequisiti

Definizione.

Una matrice $A$ di ordine $n, n \geq 2$, si dice riducibile se $\exists$ una matrice di permutazione $P$ e un intero $k, 0<k<n$ tale che:

$ B=PAP^T = \left( \begin{array}{cc} A_{11} & A_{12} \\ 0 & A_{22} \end{array} \right)$ in cui $A_{11} \in \mathbb{R}^{kxk}$ e $A_{22} \in \mathbb{R}^{(n-k)x(n-k)}$
Se la matrice A non è riducibile allora si dice irriducibile.
N.B.: Se A è riducibile, possono esistere differenti matrici di permutazione $P$ che consentono di trasformare $A$ nella forma $B$.
$A = \left( \begin{array}{ccc} 1 & 3 & 0 \\ 0 & 2 & -1 \\ -1 & 0 & 2 \end{array} \right) n=3 \Rightarrow$ avremo un grafo con 3 nodi $P_1,P_2,P_3$. Se $a_{ij} \neq 0 \Rightarrow \exists$ un arco uscente da $P_i$ entrante in $P_j$. [! Aggiungere immagine grafo] [P1 .> P1; P1 .> P2; P2 .> P2; P2 .> P3; P3 .> P3; P3 .> P1] Per determinare se $A$ di ordine $n$ è riducibile si può usare il grafo orientato associato ad $A$, cioè un grafo che ha $n$ nodi $p_i , i=1,...,n$ e un arco orientato uscente da $P_i$ ed entrante in $P_j$ se $a_{ij} \neq 0, \forall i,k$
Se il nodo di arrivo di un arco è il nodo di partenza di un altro arco $\Rightarrow$ questi sono contigui.

Definizione.

Due archi di un grafo si dicono contigui se il nodo di arrivo del primo arco è anche il nodo di partenza del secondo arco.

Definizione.

Una successione di archi contigui si dice cammino orientato

Definizione.

Un grafo orientato si dice fortemente connesso se per ogni coppia di indici $(i,j), 1 \leq i,j \leq n, i \neq j$, esiste un cammino orientato che parte da $P_i$ e arriva in $P_j$

Teorema
(CNS)

Una matrice $A$ è riducibile $\Leftrightarrow$ il suo grafo orientato NON è fortemente connesso.

Segue che C.N.:

se $A$ è irriducibile $\Rightarrow \exists$ un cammino orientato chiuso che tocca tutti i nodi

Metodo di Jacobi

Il metodo di Jacobi è anche detto "Metodo degli spostamenti simultanei", perché le componenti del vettore soluzione al passo k-esimo, $x^{(k)}$, sostituiscono simultaneamente, alla fine dell'iterata, le componenti del vettore soluzione al passo precedente, $x^{(k-1)}$. Inoltre il metodo di Jacobi risolve un sistema lineare, equivalente ad $Ax=b$, in cui la i-esima equazione, $i=1,...,m$ viene "isolata" l'incogninita i-esima, $i=1,...,n$
La matrice di iterazione di Jacobi è:
$B_j=M^{-1}N=D^{-1}(B+C) = \left( \begin{array}{ccccc} 0 & -\frac{a_{12}}{a_{11}} & -\frac{a_{13}}{a_{11}} & ... & -\frac{a_{1n}}{a_{11}} \\ -\frac{a_{21}}{a_{22}} & 0 & ... & ... & -\frac{a_{2n}}{a_{22}} \\ . & . & & & . \\ . & & . & & . \\ . & & & . & . \\ . & & & . & . \\ -\frac{a_{n1}}{a_{nn}} & -\frac{a_{n2}}{a_{nn}} & ... & ... & 0 \end{array} \right)$

Sinossi Interpolazione

Quale forma del polinomio interpolante si deve scegliere?

Dipende, oltre che dalla distribuzione dei nodi (equidistanti, Chebychev, etc.), dalla posizione del punto $\bar{x}$ in cui può essere richiesta una valutazione del $p_n(x)$:

se $\bar{x}$ è prossimo all'estremo $a$ (si costruisce solo su nodi equidistanti) $\rightarrow$ Polinomio di Newton-Gregory alle Differenze Finite in Avanti $(x_0, \Delta)$
se $\bar{x}$ è prossimo all'estremo $b$ (si costruisce solo su nodi equidistanti) $\rightarrow$ Polinomio di Newton-Gregory alle Differenze Finite all'Indietro $(x_n, \nabla)$
se $\bar{x}$ è prossimo al punto medio di $[a,b]\rightarrow$ Polinomio di Lagrange, Polinomio di Newton alle Differenze Divise

Poiché ${e_{ALG}}_{n \rightarrow +\infty} \rightarrow +\infty$, non è opportuno costruire polinomi interpolanti di grado $n \geq7$

Formula interpolante di Lagrange

Il polinomio interpolante di Lagrange che si ottiene scegliendo come funzione base $\varphi_j(x),j=0,...,n$ le seguenti funzioni:
$L_j(x)=\frac{(x-x_0)(x-x_1)...(x-x_{i-1})(x-x_{i+1})...(x-x_n)}{(x_j-x_0)(x_j-x_1)...(x_j-x_{i-1})(x_j-x_{i+1})...(x_j-x_n)}=\prod\limits_{i=0,i \neq j}^n\frac{x-x_i}{x_j-x_i}$ per $j=0,1,...,n$ tali che $L_j(x_i)=\delta_{ji}=\begin{cases}1, & \mbox{if } i=j\\0, & \mbox{if } i\neq j\end{cases}$
Allora $p_n(x)=\sum\limits_{j=0}^n y_jL_j(x)$ dove $y_j=f(x_j)$ per $j=0,...,n$ e soddisfa le condizioni di interpolazione $p_n(x_i)=\sum\limits_{j=0}^n y_jL_j(x_i)=y_i, i=0,...,n$

Osservazioni.

Ciascun $L_j(x)$ coinvolto nella forma $p_n(x)=\sum\limits_{j=0}^n y_jL_j(x)$ è un polinomio di grano $n$.
al variare di $j$, gli $L_j(x)$ sono funzioni linearmente indipendenti
Le funzioni vanno scelte sempre diverse l'una dall'altra e su uno stesso intervallo $[a,b]$

Costo computazionale del polinomio di Lagrange.

Consideriamo il caso di 4 nodi $x_0, x_1, x_2, x_3$
$p_3(x)=\sum\limits_{j=0}^3y_jL_j(x)=y_0L_0(x)+y_1L_1(x)+y_2L_2(x)+y_3L_3(x)$
$4 = n+1$ moltiplicazioni

Stimiamo il costo del calcolo di un $L_j(x)$

$L_0(x)=\frac{(x-x_1)(x-x_2)(x-x_3)}{(x_0-x_1)(x_0-x_2)(x_0-x_3)}$

Sia al numeratore che al denominatore, contiamo $n-1=2$ moltiplicazioni. A queste $2n-2$ moltiplicazioni, sommiamo anche la divisione del numeratore per il denominatore.

Allora il costo di un polinomio $L_j(x)$ è di $2n-1$.

Dovendo valutare $4=n+1$ polinomi $L_j(x), j=0,1,2,3 \Rightarrow$ il costo totale di tutti gli $L_j(x)$ è:
$(n+1)(2n-1)=2n^2+n-1$

A queste $2n^2+n-1$ operazioni sommiamo gli $(n+1)$ prodotti tra ogni polinomio $L_j(x)$ e il suo relativo $y_j$, ottenendo quindi:
$2n^2+n-1+n+1 = 2n^2+2n \Rightarrow$ il costo computazionale del polinomio interpolante di Lagrange è dell'$O(2n^2)$

Interpolazione con polinomi algebrici

Si considera come base interpolante (o modello interpolante) l'insieme dei monomi $\{\varphi_j(x)=x^t\}_{j=0}^n$

Teorema (fondamentale dell'algebra)
Se $(x_i,y_i), i=0,...,n$ sono $n+1$ punti tale che $x_i \neq x_j,i \neq j$
$\exists!$ il polinomio $p_n(x)$ di grado (al più) $n$ tale che:
$p_n(x_i)=y_i,i=0,...,n$ (detta condizione di interpolazione)
Tale polinomio è detto Polinomio Interpolante.

[... ! Aggiungere la dimostrazione ]

Il polinomio interpolante pur essendo unico può essere espresso in diversi polinomi, ci chiediamo quale sia il più conveniente.

Il polinomio interpolante, pur essendo unico, può essere rappresentato in forme diverse (ossia mediante algoritmi diversi= più convenienti sotto due punti di vista:

il costo computazionale;
la stabilità numerica del metodo.

Interpolazione polinomiale

Il problema dell'interpolazione

Data una funzione $f(x)$ definitiva su un intervallo $[a,b]$ di cui sono noti solo i valori $y_i=f(x_i)$ che essa assume in corrispondenza di punti (detti nodi dell'interpolazione) $x_i\in[a,b],i=0,1,...,n$ e fissato un insieme di funzioni $\varphi_j(x),j=0,1,...,n$ definita su $[a,b]$ e ivi linearmente indipendenti, il problema dell'interpolazione consiste nel determinare la funzione

$\hat{f}(x)=\sum\limits_{j=0}^n\alpha_j\varphi_j(x)$ (detta funzione interpolante) tale che

$\hat{f}(x_i)=f(x_i),i=0,1,...,n$.
E' importante scegliere la classe di funzioni interpolanti $\{\varphi_j(x)\}_{j=0}^n$ più adeguato al problema in esame. Le classi di funzioni più usate sono:

$\mathbb{P}_n=\{p_n(x)=\sum\limits_{i=0}^na_ix^i=a_0+a_1x^1+a_2x^2+...+a_nx^n\}$
ovvero la classe dei polinomi algebrici di grado n. Tale classe è opportuna per approssimare funzioni continue e limitate in intervalli chiusi. La determinazione di un elemento $p_n(x)\in \mathbb{P}_n$ richiede di fissare $n+1$ coefficienti $a_0, a_1, ..., a_n$
$\mathbb{T}_n=\{t_n(x)=a_0+\sum[a_k(cos(kx))+b_k(sen(kx))]\}$
classe di polinomi trigonometrici di grado $n$. Tale classe è valida per $f(x)$ periodiche e continue, per determinare un $t_n \in \mathbb{T}_n$ è necessario determinare $2n+1$ coefficienti $a_0, a_1,..., a_n, b_1, b_2,..., b_n$
$\mathbb{R}_{n, d}=\{p_n(x)/p_d(x), p_n \in\mathbb{P}_n, p_d\in\mathbb{P}_d\}$
classe di funzioni razionali. Indicata nel caso in cui $f(x)$ presenti dello singolarità o sia non periodica in intervalli infiniti
$\mathbb{E}_n$
$\mathbb{S}_n$

Criterio di Sylvester

$A \in \mathbb{R}^{nxn}$ è simmetria, $A$ è definita positiva $\Leftrightarrow$ tutti i minori principali di testa compreso $det(A)$ sono positivi.
Proprietà delle Matrici Definite Positive

$det(A) \neq 0$
Se $A$ è definita positiva $\Rightarrow$ tutte le sottomatrici sono definite positive
Se $A$ è simmetrica con autovalori $\lambda_1, \lambda_2, ..., \lambda_n$ è definita positiva $\Leftrightarrow \lambda_i > 0, \forall i =1, ..., n$
Se $A$ è definita positiva, poiché $det(A)=\lambda_1, \lambda_2, ..., \lambda_n \Rightarrow det(A)>0$

Metodo di Eliminazione di Gauss (Triangolarizzazione)

Premessa

Osserviamo preliminarmente che la soluzione di un sistema non singolare di forma triangolare, superiore o inferiore, è quasi immediata.

Infatti nel caso di un sistema triangolare superiore:

$a_{ii} \neq 0, \forall i = 1, ..., n$

la soluzione si ottiene con
$O( \frac{n^2}{2})$ operazioni moltiplicative con il seguente algoritmo: $x_i = \frac{b_i - \sum \limits_{j=i+1}^n a_{ij} x_j }{a_{ii}}, i=n-1, n-2, ..., 1$

Questo ci suggerisce che triangolarizzare un qualunque sistema con matrice quadrata è una utile strategia risolutiva.

E' noto che quando ad un'equazione di un sistema sostituiamo una combinazione lineare con un'altra equazione del sistema, il nuovo sistema risulta equivalente a quello originario, cioè pur avendo forma diversa, la soluzione è uguale.

Col metodo di Gaus è possibile trasformare in n-1 passi, ed eventuali permutazioni di righe e/o di colonne (pivoting parziale e pivoting totale), un generico sistema in uno equivalente ma triangolare.

Calcolo del costo computazionale del metodo di Gauss
Ci chiediamo:

quanto costa un passo?

Al passo k-esimo si effettuano $n-k$ divisioni

infatti si valuta

$m_{ik}=\frac{a_{ik}^{(k)}}{a_{kk}^{(k)}}$ per $i=k+1, ..., n$;
e $(n-k)(n-k+1)$ moltiplicazioni per le entrate $a_{ij}^{(k+1)}$ e $b_{i}^{(k+1)}$ per $i=k+1, ..., n$ e $j=k+1, ..., n$

quanti passi facciamo?

Poiché realizziamo $n-1$ passi allora avremo

Metodi per la risoluzione dei sistemi lineari

Due classi di metodi numeri per la risoluzione dei sistemi lineari sono:

Metodi diretti (es: Gauss, Gauss-Jordan, fattorizzazione $LU$, fattorizzazione $LL^t$ (Cholesky)
Si utilizzano su matrici dense, le quali vengono modificate passe dopo passo (fenomeno fill-in);
di dimensioni $n \leq 10^3$;
numero di passi finito: $O(n)$.
Metodi iterativi (es: Jacobi, Gauss-Seidel)
Si utilizzano su matrici sparse, che restano invariate poiché il metodo lavora su una opportuna matrice di iterazione.
di dimensioni $n > 10^3$
numero di passi (potenzialmente) infinito.

Definizione.

A è sparsa se il numero dei suoi elementi $a_{ij} \neq 0$ è dell'$O(n)$.

Viceversa si dice densa o piena.

Fenomeno del Fill-In.
Abbiamo n matrici intermedie che ci portano da A al risultato.

Complessità computazionale

L'efficienza di un metodo dipende, oltre che dalla sua stabilità numerica, anche dal suo costo computazionale (c.c.), cioè dal numero di operazioni aritmetiche elementari che esso richiede.
Di seguito considereremo il numero delle operazioni moltiplicative (cosiddette pesanti),

Prodotto matrice-vettore: $A \in \mathbb{R}^{nxm}, x \in \mathbb{R}^n \Rightarrow y = Ax, y \in \mathbb{R}^n $
$y_i = a_{i1}x_1+a_{i2}x_2+...+a_{in}x_n = \sum\limits_{j=1}^n a_{ij}x_j, \forall i = 1, ..., n$
Ciascuna componente $y_i$ del vettore $y \in \mathbb{R}^n$ costa $n$ prodotti, le componenti sono $n$
$ \Rightarrow $ il costo computazionale del prodotto matrice vettore è dell' $O(n^2)$
Prodotto righe per colonne: $A, B \in \mathbb{R}^{nxn} \Rightarrow \mathbb{C} = \mathbb{r}^{nxn}$
$c_{ij} = \sum\limits_{k=1}^n a_{ik}b_{kj}=a_{i1}b_{1j}+a_{i2}b_{2j}+...+a_{in}b_{nj}$
Abbiamo quindi $n$ prodotti e $n-1$ addizioni (ma come detto prima ci occuperemo solo delle moltiplicazioni).
Il generico elemento $c_{ij}$ ha un costo computazionale dell'$O(n)$ prodotti. Gli elementi $c_{ij}$ della matrice $C \in \mathbb{R}^{nxn}$ sono in numero $n^2$.
$\Rightarrow$ il costo computazionale del prodotto righe-colonne tra $A, B \in \mathbb{R}^{nxn}$ è dell'$O(n^3)$.

Se $A \in \mathbb{R}^{mxn}, B \in \mathbb{R}^{nxp} \Rightarrow C = AB \in \mathbb{R}^{mxp}$, il quale c.c. è dell'$O(mnp)$

Norme Matriciali

Definizione.

Una norma matriciale è una funzione $||.||: M_n(\mathbb{C}) \rightarrow \mathbb{R}_+ \cup \{ 0 \}$ che soddisfa le seguenti proprietà:

$||A|| \geq 0$ e $||A||=0 \Leftrightarrow A=\Omega, \forall A \in M_n$
$|| \alpha A || = | \alpha | ||A||, \forall a \in \mathbb{C}$ e $\forall A \in M_n $
$||A+B|| \leq ||A|| + ||B||, \forall A,B \in M_n$
$||AB|| \leq ||A||||B||, \forall A,B \in M_n$

Definizione.

La norma matriciale definita da $||A||_M = \max_{||x||=1} ||A X ||_v$ viene detta Norma Indotta.
Le 3 norme vettoriali $||.||_1, ||.||_2, ||.||_\infty$ inducono delle corrispondenti norme matriciali

norma 1 (somma più alta tra le somme degli elementi di una colonna): $||A||_1 = \max\limits_{1 \leq j \leq n } \sum\limits_{i=1}^n |a_{ij}|$
norma 2: $||A||_2= \sqrt{ \rho (A^t A)}$
norma $\infty$ (somma più alta tra le somme degli elementi di una riga): $||A||_\infty = \max\limits_{1 \leq i \leq n} \{ \sum\limits_{j=1}^n |a_{ij}| \}$

Spazi vettoriali normati

Definizione.

Una norma è una funzione $||.||: V \rightarrow \Re_+ U \{ 0 \}$ che soddisfa le seguenti proprietà:

$||x|| \geqslant 0$ e $||x||=0 \Leftrightarrow x = 0, x \in V$
$|| \alpha || = | \alpha | * || x ||, \forall \alpha \in \Re$
$|| x+y || \leq || x || + ||y||, \forall x, y \in V$

Le norme vettoriali più comunemente usate in $\mathbb{R}^n(\mathbb{C}^n)$:

norma 1: $||x||_1 = \sum\limits_{i=1}^n |x_i| $ (Array-sum)
norma 2: $||x||_2 = \sqrt{x^t x} = \sum\limits_{i=1}^n |x_i|^2 $
norma $\infty$: $||x||_\infty = \max_{1 \leq i \leq n}\{ |x_i| \}$ maxAbs(array a)

Algoritmo per la determinazione della EPS (Precisione di macchina)

Algorithm pseudo code (Cleve Moler)

a = 4/3

b = a - 1

c = b + b + b

e = 1 - c

eps = e/2

Algoritmo per la determinazione della EPS (Precisione di macchina)

La eps di macchina è caratteristica di ogni calcolatore ed è la più piccola quantità positiva significativa (ossia percepita) in macchina, tale che sommata a 1 fornisce un valore 1+eps>1.

Algorithm pseudo code

eps = 1

repeat

eps = eps/2

until eps + 1 = 1

eps = eps*2