\documentclass[10pt,landscape]{article} \usepackage{amssymb,amsmath,amsthm,amsfonts} \usepackage{multicol,multirow} \usepackage{marvosym} \usepackage{calc} \usepackage{ifthen} \usepackage[landscape]{geometry} \usepackage[colorlinks=true,citecolor=blue,linkcolor=blue]{hyperref} \usepackage{notes_2023} \usepackage{quiver} \setlength{\extrarowheight}{0pt} \renewcommand{\sp}{\operatorname{sp}} \ifthenelse{\lengthtest { \paperwidth = 11in}} { \geometry{top=.5in,left=.5in,right=.5in,bottom=.5in} } {\ifthenelse{ \lengthtest{ \paperwidth = 297mm}} {\geometry{top=1cm,left=1cm,right=1cm,bottom=1cm} } {\geometry{top=1cm,left=1cm,right=1cm,bottom=1cm} } } \DeclareMathOperator{\DFT}{DFT} \DeclareMathOperator{\IDFT}{IDFT} \newcommand{\ff}{\mathcal{F}} \newcommand{\ein}{\varepsilon_{\text{in}}} \newcommand{\ealg}{\varepsilon_{\text{alg}}} \newcommand{\ean}{\varepsilon_{\text{an}}} \newcommand{\etot}{\varepsilon_{\text{tot}}} \renewcommand{\eps}{\varepsilon} %\pagestyle{empty} \newcommand{\tx}{\Tilde{x}} \newcommand{\ty}{\Tilde{y}} \newcommand{\tz}{\Tilde{z}} \makeatletter \renewcommand{\section}{\@startsection{section}{1}{0mm}% {-1ex plus -.5ex minus -.2ex}% {0.5ex plus .2ex}%x {\normalfont\large\bfseries}} \renewcommand{\subsection}{\@startsection{subsection}{2}{0mm}% {-1explus -.5ex minus -.2ex}% {0.5ex plus .2ex}% {\normalfont\normalsize\bfseries}} \renewcommand{\subsubsection}{\@startsection{subsubsection}{3}{0mm}% {-1ex plus -.5ex minus -.2ex}% {1ex plus .2ex}% {\normalfont\small\bfseries}} \makeatother \setcounter{secnumdepth}{0} \setlength{\parindent}{0pt} \setlength{\parskip}{0pt plus 0.5ex} % ----------------------------------------------------------------------- \title{Schede riassuntive di Analisi Numerica} \begin{document} \parskip=0.7ex \raggedright \footnotesize \begin{center} \Large{\textbf{Schede riassuntive di Analisi Numerica}} \\ \end{center} \begin{multicols}{3} \setlength{\premulticols}{1pt} \setlength{\postmulticols}{1pt} \setlength{\multicolsep}{1pt} \setlength{\columnsep}{2pt} \section{Analisi dell'errore} Se $\tilde{x}$ è un'approssimazione di $x$ definiamo $\tilde{x}-x$ come l'\textit{errore assoluto}, $\eps = (\tx-x)/x$ come l'\textit{errore relativo} (rispetto a $x$) e $\eta = (x-\tx)/\tx$ come l'\textit{errore relativo} (rispetto a $\tx$). In particolare vale che $\tilde{x} = x(1+\eps) = x/(1+\eta)$. \subsection{Rappresentazione in base} Sia $B\ge2$ un numero intero, vale allora il seguente risultato: \begin{theorem}[di rappresentazione in base] Per ogni numero reale $x\neq0$ esistono e sono unici un intero $p$ ed una successione $\{d_i\}_{i\ge1}$ con le seguenti proprietà: \begin{enumerate}[(1.)] \item $0\le d_i\le B-1$ \item $d_1\neq0$ \item per ogni $k>0$ esiste un $j\ge k$ tale che $d_j\neq B-1$ (ossia $\{d_i\}_{i \geq 1}$ è frequentemente diversa da $B-1$) \item $x$ si scrive come: $$x=\text{sgn}(x)B^p\sum_{i\ge1}d_iB^{-i}$$ \end{enumerate} \end{theorem} L'intero $B$ è detto \textit{base della rappresentazione}, gli interi $d_i$ sono dette \textit{cifre della rappresentazione} ed il numero $\sum_{i\ge1}d_iB^{-i}$ viene chiamato \textit{mantissa}. Si scrive $\pp(x)$ per indicare l'esponente $p$ relativo a $x$. \medskip La condizione ($2$.) è detta di \textit{normalizzazione} e serve a garantire l'unicità della rappresentazione e a memorizzare il numero in maniera più efficiente, mentre la condizione $3$ esclude rappresentazioni con cifre uguali a $B-1$ da un certo punto in poi (per esempio la rappresentazione di $1$ come $0.\overline{9}$ è esclusa). \subsection{Numeri \textit{floating point}} Dati gli interi $B\ge2$, $t\ge1$ ed $M$, $m>0$, si definisce l'insieme $\mathcal{F}(t,B,m,M)$ dei \textbf{numeri di macchina} o \textbf{dei numeri in virgola mobile} o ancora dei \textbf{numeri \textit{floating point}} come l'insieme: $$\zeroset\cup\{\pm B^p\sum_i^td_iB^{-i}|d_1\neq0,0\le d_i\le B-1,-m\le p\le M\},$$ in particolare $B$ rappresenta la base della rappresentazione, $t$ la lunghezza della mantissa, $B^{-m}$ il minimo numero che moltiplica la mantissa e $B^M$ il massimo. \subsubsection{Troncamento su \texorpdfstring{$\RR$}{ℝ} e \texorpdfstring{$\RR^n$}{ℝⁿ} e numero di macchina \texorpdfstring{$u$}{u}} Sia $x$ un numero reale. Se $-m \leq \pp(x) \leq M$, allora $x$ viene ben rappresentato in $\mathcal{F}$ dal numero $\tx=\fl(x)$, dove: \[ x=\sgn(x)B^p\sum_{i\ge1}d_iB^{-i} \implies \tx=\sgn(x)B^p\sum_{i=1}^t d_iB^{-i}, \] ossia $\tx$ è ottenuto da $x$ troncandone la mantissa al $t$-esimo termine (\textbf{troncamento}). Facendo così si ottiene un errore relativo di rappresentazione tale per cui: \[ \abs{\frac{\tx-x}{x}}M$ il numero non è rappresentabile e ci ritroviamo rispettivamente in una situazione di \textit{underflow} o di \textit{overflow}. \medskip Dato $x\in\RR^n$, si definisce la rappresentazione $\tx$ rispetto a $x$ in $\mathcal{F}$ come: \[ \tx=(\tx_i)=(\fl(x_i)=(x_i(1+\varepsilon_i)), \quad \abs{\varepsilon_i}(t_2-1)\frac{\log B_2}{\log B_1}+1 \] Si può definire in modo analogo l'approssimazione per arrotondamento, osservando che in questo caso l'errore relativo di rappresentazione è limitato da $\frac{1}{2}B^{1-t} = \frac{1}2 u$. \subsubsection{Aritmetica di macchina} Siano $a$, $b\in\mathcal{F}(t,B,m,M)$ e sia $\op$ una delle quattro operazioni aritmetiche ($+$, $-$, $\cdot$ e $/$). Si consideri $c=a \bop b$. Allora la macchina calcola $c$ con l'approssimazione $\hat{c} := \fl(c)$. Vale pertanto che: \[ \hat{c}=\fl(a \bop b)=c(1+\delta)=c/(1+\eta)\;\con \abs\delta, \abs\etar_i = \sum_{j=1, \, j\neq i}^n \abs{a_{ij}}$ $\forall i=1$, ..., $n$. Per il primo teorema di Gershgorin, $A$ non può essere singolare ($0$ non appartiene per ipotesi a nessun cerchio). Ogni sottomatrice principale di una matrice fortemente diagonale in senso stretto è allo stesso modo fortemente diagonale. In particolare, se $A$ è fortemente dominante diagonale, allora $a_{ii} \neq 0$. \medskip Si possono definire in modo analogo le matrici dominanti diagonali (in senso debole) cambiando $>$ in $\geq$, ma non è detto che tali matrici siano non singolari. Un controesempio è infatti la matrice: \[\begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix}. \] \subsubsection{Secondo teorema di Gershgorin} \begin{theorem}[Secondo teorema di Gershgorin] Sia $A\in \CC^{n \times n}$ e sia $K=\bigcup_{i=1}^n K_i$ l'unione dei cerchi $K_i$ relativi ad $A$. Sia inoltre \[ K=M_1\cup M_2\;\con M_1\cap M_2=\emptyset, \] dove $M_1$ è costituito da $n_1$ cerchi e $M_2$ è costituito da $n_2$ cerchi. Allora $M_1$ contiene $n_1$ autovalori e $M_2$ ne contiene $n_2$. \end{theorem} \begin{note} La dimostrazione di questo teorema sfrutta un argomento di continuità riguarda al segmento $A(t) = D + t(A-D)$ con $t \in [0, 1]$ e $D = \diag(A)$. Un simile argomento di continuità può risultare utile in altri contesti. \end{note} Quest'ultimo teorema risulta utile per dimostrare che una matrice $A \in \RR^{n \times n}$ ha un autovalore reale. Se infatti $K_i$ è un cerchio di $A$ disgiunto dagli altri, allora $K_i$ contiene un unico autovalore, che deve essere dunque necessariamente reale (altrimenti, siccome $A$ è reale, vi sarebbe anche il suo coniugato, e quindi vi sarebbero almeno $2$ autovalori, assurdo). \subsubsection{Terzo teorema di Gershgorin} \begin{theorem}[Terzo teorema di Gershgorin] Sia $\lambda$ un autovalore di $A$. Si ipotizzi che val le seguenti condizioni \begin{enumerate}[(1.)] \item $A$ è irriducibile, \item $\lambda \in K_i \implies \lambda \in \partial K_i$ (dove $\partial K_i$ è la frontiera di $K_i$). \end{enumerate} Allora $\lambda \in \bigcap \partial K_i$. \end{theorem} \subsubsection{Matrici irriducibilmente dominanti diagonali (i.d.d.)} Una matrice $A\in \CC^{n \times n}$ si dice \textbf{irriducibilmente dominante diagonale} se \begin{enumerate}[(i)] \item $A$ è irriducibile, \item $A$ è dominante diagonale \underline{in senso debole}, \item $\exists h$ t.c.~$\abs{a_{hh}}>r_h = \sum_{j=1, \, j\neq h}^n{a_{hj}}$ (ossia con $0 \notin K_h$). \end{enumerate} Per il terzo teorema di Gershgorin le matrici i.d.d.~sono non singolari. Inoltre, se $A$ è i.d.d., allora $a_{ii} \neq 0$ per ogni $i$. %\\Sia $p(x)=x^n+\sum_{i=0}^{n-1}a_ix^i$ un polinomio a coefficienti in $\CC$. Sia $C$ la matrice compagna di $p(x)$ allora trovare gli zeri di $p(x)$ corrispondono agli autovalori di $C$. (si mostra per induzione che $\det(xI-C)=p(x)$). Si possono dunque localizzare gli zeri di un polinomio tramite i teoremi di Gershgorin. Se $a_0\neq0$ si può utilizzare lo stesso risultato a $q(x)=(1/a_0)x^np(x\inv)$ che fornisce inclusioni per i reciproci degli zeri. \\Un'altra localizzazione degli zeri consiste nel costruire la matrice $A=D-ue^t$ dove $D=diag(x_1,\dots,x_n),e^t=(1,\dots,1)\E u=(u_i)\;\con u_i=p(x_i)/\prod_{j=1,j\neq i}(x_i-x_j)$. Il teorema di Gershgorin applicato ad A fornisce delle stime dell'errore con cui i valori $x_i$ approssimano gli zeri. \subsection{Forma normale di Schur} \begin{theorem} Data $A\in \CC^{n \times n}$, esiste una matrice $U\in \CC^{n \times n}$ unitaria, ovverosia tale per cui $U^HU=UU^H=I_n$, tale che: \[ U^HAU=T\in T_+(n,\CC), \] ovverosia con $U^H A U$ triangolare superiore. Si dice in tal caso che $T$ è una \textbf{forma normale di Schur} di $A$. \end{theorem} Generalmente esistono più forme normali di Schur per una matrice, e possono essere ottenute riordinando gli autovalori e/o le basi scelte. \subsubsection{Teorema spettrale} \begin{theorem}[spettrale] Se $A$ è hermitiana, allora una sua forma normale di Schur è sempre una matrice diagonale con elementi reali, ovverosia $A$ è ortogonalmente diagonalizzabile con autovalori reali. Se invece $A$ è anti-hermitiana ($A^H = -A$), allora una sua forma normale di Schur è una matrice diagonale con elementi immaginari puri, ovverosia $A$ è ortogonalmente diagonalizzabile con autovalori immaginari puri. \end{theorem} \subsubsection{Caratterizzazione delle matrici normali} \begin{definition} Una matrice $A \in \CC^{n \times n}$ si dice \textbf{normale} se $A A^H = A^H A$. \end{definition} Si enuncia una caratterizzazione delle matrici triangolari normali: \begin{proposition} Una matrice triangolare $T$ è normale se e solo se è diagonale. \end{proposition} Dalla precedente proposizione si ottiene facilmente la seguente fondamentale altra caratterizzazione: \begin{theorem} Una matrice $A \in \CC^{n \times n}$ è normale se e solo se la sua forma normale di Schur è diagonale (i.e.~è ortogonalmente diagonalizzabile). \end{theorem} \subsubsection{Autovalori di una matrice unitaria} Gli autovalori di una matrice unitaria $A$ hanno modulo $1$. Se infatti $v$ è un autovettore relativo all'autovalore $\lambda$, vale che: \[ A v = \lambda v \implies v^H A^H = \overline\lambda v^H, \] da cui: \[ v^H v = v^H A^H A v = \overline\lambda v^H A v = \overline\lambda \lambda v^H v \implies \abs{\lambda} = 1. \] In particolare, se $A$ è reale, $\lambda \in \{\pm1\}$. \subsubsection{Forma normale di Schur reale} \begin{definition}[Matrice quasi-triangolare superiore] Una matrice $T\in \RR^{n \times n}$ si dice \textbf{quasi-triangolare superiore} se si può scrivere nella forma: \[ \begin{pmatrix} T_{11} & \dots & T_{1n} \\ & \ddots & \vdots \\ & & T_{mm} \end{pmatrix}, \] dove i blocchi matriciali $T_{ii}$ possono essere matrici $2 \times 2$ oppure matrici $1 \times 1$, ossia numeri reali. \end{definition} Si definisce analogamente la nozione di matrice quasi-triangolare \textit{inferiore}. \smallskip Gli autovalori di una matrice quasi-triangolare $T$ sono gli autovalori delle sottomatrici $T_{ii}$. Il polinomio caratteristico di $T$ è il prodotto dei polinomi dei blocchi $T_{ii}$. Il determinante di $T$ è il prodotto dei determinanti dei blocchi $T_{ii}$, la traccia è la somma delle tracce dei blocchi $T_{ii}$. \smallskip Usando la nozione di matrice quasi-triangolare si può enunciare un teorema analogo a quello della forma normale di Schur, ma riadattato per le matrici reali: \begin{theorem}[Forma normale di Schur reale] Per ogni matrice $A\in \RR^{n \times n}$ esistono $T\in \RR^{n \times n}$ quasi-triangolare superiore e $Q \in O(n)$ tali per cui: \[ Q^tAQ=T. \] \end{theorem} \subsection{Norme di vettori} Una \textbf{norma} (vettoriale) su $\CC^n$ è un'applicazione \[ \norm{\cdot}:\CC^n\longrightarrow\RR \] tale per cui: \begin{enumerate}[(1.)] \item $\forall x\in\CC^n$, $\norm{x} \geq 0$ e $\norm{x}=0\iff x=0$ (definitezza positiva). \item $\forall\alpha\in\CC$, $\forall x\in\CC^n$, $\norm{\alpha x}=\abs\alpha\norm{x}$ (omogeneità). \item $\forall x$, $y\in\CC^n$, $\norm{x+y}\le\norm{x}+\norm{y}$ (disuguaglianza triangolare). \end{enumerate} Ogni prodotto hermitiano $\langle \cdot, \cdot \rangle$ definito positivo di $\CC^n$ induce una norma ponendo $\norm{v} = \sqrt{\langle v, v \rangle}$. \medskip Per $p \in [1, \infty)$ si dice \textbf{norma di Hölder} di ordine $p$ la norma: \[ \norm{x}_p=\left(\sum_{i=1}^n|x_i|^p\right)^{\frac{1}{p}}. \] Alcuni esempi noti di norme di Hölder sono: \begin{itemize} \item $\norm{x}_1=\sum_{i=1}^n \abs{x_i}$, la norma $1$, \item $\norm{x}_2=\left(\sum_{i=1}^nx_i^2\right)^{\frac{1}{2}} = \sqrt{x^H x}$, detta anche \textit{norma euclidea}, indotta dal prodotto hermitiano standard di $\CC^n$, \item $\norm{x}_{\infty}:=\lim_{p \to \infty} \norm{x}_p = \max_{i}\abs{x_i}$. \end{itemize} Per una norma l'insieme $S = \{x\in\CC^n \mid \norm{x}\le1\}$ è un insieme convesso. Da ciò si deduce che per $0 0$, esiste $\delta > 0$ tale per cui: \[ \forall x,y\in\CC^n, \norm{x-y}_2 < \delta \implies \abs{\norm{x}-\norm{y}} < \varepsilon, \] espressione che segue dal fatto che una norma è Lipschitziana (deriva dalla \textit{disuguaglianza triangolare}). Equivalentemente vale che per ogni $\eps > 0$ esiste $\delta > 0$ tale per cui: \[ \forall x,y\in\CC^n, \abs{x_i-y_i} < \delta \implies \abs{\norm{x}-\norm{y}} \leq \varepsilon. \] A partire da questo risultato si può dimostrare il seguente fondamentale teorema: \begin{theorem}[Equivalenza tra norme] Per ogni coppia di norme $\norm{\cdot}'$ e $\norm{\cdot}''$ su $\CC^n$, esistono due costanti positive $\alpha$ e $\beta$ (dipendenti da $n$) tali per cui: \[ \alpha\norm{x}'\le\norm{x}''\le\beta\norm{x}' \quad \forall x\in\CC^n. \] In altre parole, le norme su $\CC^n$ sono bi-lipschitziane tra loro. \end{theorem} Infatti $S = \{ x \in \CC^n \mid \norm{x}_2 = 1 \}$ è un compatto euclideo (essendo chiuso e limitato, per il teorema di Heine-Borel), e $\norm{\cdot}$ è continua sulla topologia euclidea. Pertanto tutte le norme di $\CC^n$ inducono la stessa topologia, ossia quella indotta da $\norm{\cdot}_2$, e quindi coincidono tutti gli aperti e i chiusi. \smallskip In particolare per le norme $1$, $2$ e $\infty$ vale che: \begin{itemize} \item $\norm{x}_\infty \leq \norm{x}_1 \leq n \norm{x}_\infty$, \item $\norm{x}_2 \leq \norm{x}_1 \leq \sqrt{n} \norm{x}_2$, \item $\norm{x}_\infty \leq \norm{x}_2 \leq \sqrt{n} \norm{x}_\infty$. \end{itemize} Se $U$ è una matrice unitaria, allora $U$ agisce lasciando invariata la norma $2$, ovverosia $\norm{Ux}_2 = \norm{x}_2$. \subsection{Norme di matrici} Una \textbf{norma matriciale} è un'applicazione \[ \norm{\cdot}:\CC^{n \times n}\to\RR \] per la quale valgono le seguenti quattro proprietà: \begin{enumerate}[(1.)] \item $\forall A\in \CC^{n \times n}$, $\norm{A}\ge0$ e $\norm{A}=0\iff A=0$ (definitezza positiva), \item $\forall \lambda\in\CC$, $\forall A\in \CC^{n \times n}$, $\norm{\lambda A}=\abs\lambda\norm{A}$ (omogeneità), \item $\norm{A+B}\le\norm{A}+\norm{B}$ $\forall A$,$B\in \CC^{n \times n}$ (disuguauglianza triangolare), \item $\norm{AB}\le\norm{A}\norm{B}$ $\forall A$, $B\in \CC^{n \times n}$ (proprietà submoltiplicativa). \end{enumerate} \subsubsection{Norma di Frobenius} Si definisce la \textbf{norma di Frobenius} l'applicazione che agisce su $A \in \CC^{n \times n}$ in modo tale che: \[ \norm{A}_F=\left(\sum_{i=1}^n\sum_{j=1}^n\abs{a_{ij}}^2\right)^{\frac{1}{2}}=\sqrt{\tr(A^HA)}. \] In particolare la norma di Frobenius è esattamente la norma euclidea dello spazio $\CC^{n \times n}$ immerso in $\CC^{n^2}$ a cui è naturalmente isomorfo associando ad $A$ il vettore $(A_i^\top)_i$ ottenuto trasponendo le righe e sovrapponendole ordinatamente. \subsubsection{Norme matriciali indotte} \begin{note} Si ricorda che, data una norma $\norm{\cdot}$, l'insieme $S=\{x\in\CC^n \tc \norm{x}=1\}$ è chiuso e limitato (in tutte le topologie indotte da norme, essendo equivalenti). Allora, per il teorema di Heine-Borel, $S$ è compatto, e dunque una funzione continua $f : S \to \RR$ ammette massimo per il teorema di Weierstrass. \end{note} \begin{note} Una mappa lineare da $\CC^n$ in sé è continua. Dunque è possibile applicare il teorema di Weierstrass alla funzione $\norm{ \cdot } \circ f_A$ ristretta su $S$, dove $f_A$ è l'app.~indotta da una matrice $A$. \end{note} \begin{definition} Data una norma vettoriale $\norm{\cdot}:\CC^n\to\RR$, si definisce la \textbf{norma matriciale indotta} (o la corrispondente \textit{norma operatore}), come la norma che agisce su $A \in \CC^{n \times n}$ in modo tale che: \[ \norm{A}:=\max_{\norm{x}=1} \norm{Ax}=\max_{x \in S} \norm{Ax}=\max_{x \in \CC^n \setminuszero} \frac{\norm{Ax}}{\norm{x}}. \] \end{definition} Per ogni norma operatore valgono le due seguenti aggiuntive proprietà, \textit{oltre quelle caratterizzanti una norma}: \begin{enumerate}[(1.)] \item $\norm{Ax}\le\norm{A}\norm{x}$, $\forall A\in \CC^{n \times n}$, $\forall x\in \CC^n$, \item $\norm{I}=1$. \end{enumerate} Poiché $\norm{I}_F = \sqrt{n}$, la norma di Frobenius non è in generale indotta da alcuna norma vettoriale. \smallskip Per le norme $1$, $2$, $\infty$ valgono le seguenti identità: \begin{itemize} \item $\displaystyle \norm{A}_1=\max_{j=1,\dots,n}\sum_{i=1}^n \abs{a_{ij}} = \max_{j=1,\ldots,n} \norm{A^j}_1$, \item $\displaystyle \norm{A}_\infty=\max_{i=1,\dots,n}\sum_{j=1}^n \abs{a_{ij}} = \max_{i=1,\ldots,n} \norm{A_i}_1$, \item $\norm{A}_2=\sqrt{\rho(A^HA)}$, dove $\rho(A^H A)$ è l'autovalore di valore assoluto maggiore di $A^H A$ (\textit{raggio spettrale}). \end{itemize} Si osserva immediatamente che la norma matriciale $2$ è computazionalmente più difficile da calcolare rispetto alle norme $1$ e $\infty$. Inoltre, se $A$ è simmetrica, le norme $1$ e $\infty$ coincidono. \smallskip Poiché $\norm{A}_F = \sqrt{\tr(A^H A)}$ è la radice della somma degli autovalori di $A^H A$, vale in particolare che $\norm{A}_F \geq \norm{A}_2$. Si osserva che $A^H A$ è semidefinita positiva, e dunque i suoi autovalori sono non negativi. \smallskip Se $U$ è unitaria e $B = UA$, vale che: \[ B^H B = (UA)^H UA = A^H U^H U A = A^H A, \] e quindi $B$ e $A$ condividono la stessa norma di Frobenius e la stessa norma $2$. Analogamente si vede che $AU$ e $A$ condividono le stesse due norme. Equivalentemente, la moltiplicazione per matrice unitaria (sia a destra che a sinistra) lascia invariata sia la norma $2$ che quella di Frobenius. In particolare, se $U$ e $V$ sono unitarie, $\norm{UAV}_2 = \norm{A}_2$ e $\norm{UAV}_F = \norm{A}_F$. \smallskip Se $A$ è normale, allora gli autovalori di $A^H A$ sono i moduli quadrati degli autovalori di $A$. Pertanto, per $A$ normale, $\norm{A}_2 = \max_{\lambda \in \sp(A)} \abs{\lambda} = \rho(A)$ e $\norm{A}_F = \sqrt{\sum_{\lambda \in \sp(A)} \abs{\lambda}^2}$. Si è usato che la forma normale di Schur di $A$ è diagonale e che le trasformazioni unitarie non variano né $\norm{A}_2$ né $\norm{A}_F$. \smallskip Valgono inoltre le seguenti altre due disuguaglianze: \begin{itemize} \item $\norm{A}_F \leq \sqrt{r} \norm{A}_2 \leq \sqrt{n} \norm{A}_2$, dove $r = \rg(A)$, \item $\norm{A}_2^2 \leq \norm{A}_1 \cdot \norm{A}_\infty$. \end{itemize} \subsubsection{Norma indotta da una matrice non singolare \texorpdfstring{$S$}{S}} Sia $\norm{\cdot}$ una norma vettoriale e $S\in \CC^{n \times n}$ una matrice non singolare. Allora, data la norma vettoriale $\norm{x}_S := \norm{Sx}$, vale che: \[ \norm{A}_S=\max_{\norm{x}_S=1}\norm{Ax}_S=\max_{\norm{x}_S=1}\norm{SAx}= \max_{\norm{Sx}=1}{\norm{SAS\inv Sx}}, \] e quindi, sfruttando che $S$ è non singolare: \[ \norm{A}_S=\max_{\norm{y}=1}{\norm{SAS\inv y}}, \] ossia vale che: \[ \norm{A}_S = \norm{SAS\inv}. \] \subsubsection{Norme e raggi spettrali} \begin{definition} Data $A\in \CC^{n \times n}$, si definisce \textbf{raggio spettrale} $\rho(A)$ il modulo dell'autovalore massimo di $A$, ossia: \[ \rho(A)=\max\{\abs\lambda \mid \lambda\text{ autovalore di } A\}. \] \end{definition} Se $x$ è un autovettore relativo a $\rho(A)$ di modulo unitario rispetto a una norma $\norm{\cdot}$, allora vale che: \[ \norm{Ax} = \norm{\rho(A) x} = \rho(A) \norm{x} = \rho(A), \] e dunque $\rho(A) \leq \norm{A}$ per ogni $A$. In particolare $\norm{A} \geq \abs{\lambda}$ per ogni autovalore $\lambda$. \smallskip Si enunciano inoltre i seguenti teoremi: \begin{theorem} Sia $A \in \CC^{n \times n}$. Allora per ogni $\varepsilon > 0$ esiste una norma indotta $\norm{\cdot}$ tale per cui: \[ \rho(A) \leq \norm{A} \leq \rho(A) + \varepsilon. \] Inoltre, se gli autovalori di modulo massimo di $A$ appartengono solo a blocchi di Jordan di taglia $1$, allora esiste una norma per cui $\norm{A} = \rho(A)$. \end{theorem} \begin{theorem} Sia $\norm{\cdot}$ una norma matriciale. Allora vale la seguente identità: \[ \lim_{k\to+\infty}\norm{A^k}^{\frac{1}{k}}=\rho(A). \] \end{theorem} Se $A\in \CC^{n \times n}$ è tale per cui $\norm{A}<1$ dove $\norm{\cdot}$ è una norma matriciale indotta, allora $1$ non può essere autovalore di $A$, e dunque $I-A$ è invertibile. Inoltre vale che:\[ \norm{(I-A)\inv}\le\frac{1}{1-\norm{A}}. \] \subsection{Condizionamento di un sistema lineare e numero di condizionamento} Data $A\in \CC^{n \times n}$ non singolare e dato un vettore $b\in\CC^n\setminuszero$ si vuole studiare il condizionamento del problema $Ax=b$, ossia di un sistema lineare. \medskip Consideriamo il problema $(A+\delta_A)y=b+\delta_b$, dove perturbiamo il sistema originale mediante dei parametri $\delta_A$ e $\delta_B$ di cui conosciamo i rapporti $\frac{\norm{\delta_b}}{\norm{b}}$ e $\frac{\norm{\delta_A}}{\norm{A}}$, cercando di ottenere informazioni riguardo $\frac{\norm{\delta_x}}{\norm{x}}$, sostituendo $y=x+\delta_x$. \medskip \begin{definition} Si dice \textbf{numero di condizionamento} $\mu(A)$ di $A$ nella norma $\norm{\cdot}$ il valore: \[ \mu(A)=\norm{A}\norm{A\inv}. \] Si scrive $\mu_p(A)$ per intendere $\norm{A}_p \norm{A\inv}_p$. \end{definition} Studiamo in particolare la perturbazione di $Ax=b$ nel caso in cui $\delta_A=0$. \begin{proposition} Se $\delta_A = 0$, allora vale la seguente disuguaglianza: \[ \frac{\norm{\delta_x}}{\norm{x}}\le\norm{A}\norm{A\inv}\frac{\norm{\delta_b}}{\norm{b}} = \mu(A) \frac{\norm{\delta_b}}{\norm{b}}. \] In generale, per $\delta_A \neq 0$ e $A + \delta_A$ invertibile, vale che: \[ \varepsilon_x\le\frac{\norm{A}\norm{A\inv}}{1-\varepsilon_A\norm{A}\norm{A\inv}}(\varepsilon_A+\varepsilon_B), \] dove $\eps_t = \norm{\delta_t}/\norm{t}$. \end{proposition} Pertanto il sistema è \textbf{ben condizionato} se $\mu(A)$ è relativamente piccolo. \smallskip Per il numero di condizionamento valgono le seguenti proprietà: \begin{itemize} \item $\mu(A) \geq \norm{I}$ per la \textit{proprietà submoltiplicativa}, \item $\mu(A) \geq 1$ per $\norm{\cdot}$ norma operatore ($\norm{I} = 1$), \item $\mu_2(U)=1$ per $U$ unitaria. \end{itemize} Per $A$ normale vale la seguente identità: \[\mu_2(A) = \frac{\max_{\lambda \in \sp(A)} \abs{\lambda}}{\min_{\lambda \in \sp(A)} \abs{\lambda}}. \] Inoltre vale la seguente disuguaglianza: \[ \mu_2(A) \leq \mu(A), \] dove $\mu(A)$ è riferito a qualsiasi altra norma operatore. \section{Metodi diretti per sistemi lineari} Ci si propone di risolvere il sistema $Ax=b$ scrivendo $A$ come prodotto di matrici ``consone'' e facilmente invertibili. Se infatti $A=PQ$, allora il sistema $PQx=b$ può essere risolto come: \[ \begin{cases} Py=b \\ Qx=y \end{cases} \] risolvendo dunque prima $Py=b$ e poi $Qx=y$. \subsection{Risoluzione di \texorpdfstring{$Ax=b$}{Ax=b} per \texorpdfstring{$A$}{A} triangolare o unitaria} Sono presentati di seguito tre tipi di matrice per le quali l'invertibilità è garantita e per cui il sistema $Ax=b$ è facilmente risolvibile. \smallskip Se $A$ è \textbf{triangolare inferiore} con $a_{ii}\neq0$ $\forall i$, allora, per risolvere $Ax=b$, è possibile applicare il \textit{metodo di sostituzione in avanti} ponendo: \[ x_1=\frac{b_1}{a_1}, \quad x_i=\frac{1}{a_{ii}}\left(b_i-\sum_{j=1}^{i-1}a_{ij}x_j\right)\;\,\con i\geq 2. \] Se $A$ è \textbf{triangolare superiore} con $a_{ii}\neq0$ $\forall i$, allora, per risolvere $Ax=b$, è possibile applicare il \textit{metodo di sostituzione all'indietro} ponendo: \[ x_n=\frac{b_n}{a_n}, \quad x_i=\frac{1}{a_{ii}}\left(b_i-\sum_{j=i+1}^n a_{ij} x_j\right)\;\,\con i < n.\] I due algoritmi presentati hanno un costo computazionale di $n^2$ flops e sono entrambi numericamente stabili all'indietro. \smallskip Se $A$ è \textbf{unitaria}, $Ax=b \implies x=A^Hb$, e dunque si verifica che: \[ x_j=\sum_{i=1}^n \overline{a_{ij}} b_i. \] Questo algoritmo richiede un costo computazionale di $2n^2-n$ flops ed è ancora numericamente stabile all'indietro. \subsection{Fattorizzazione classiche di matrici} In letteratura esistono $4$ fattorizzazioni classiche: \begin{enumerate} \item $A=LU$ con $L$ triangolare inferiore con tutti $1$ sulla diagonale e $U$ triangolare superiore (possibile solo per alcune classi di matrici); \item $A=PLU$ con $L$, $U$ come sopra e $P$ matrice di permutazione (sempre possibile); \item $A=P_1LUP_2$ con $L$, $U$ come sopra e $P_1$, $P_2$ matrici di permutazione (sempre possibile); \item $A=QR$ con $Q$ unitaria e $R$ triangolare superiore (sempre possibile). \end{enumerate} \subsubsection{Condizioni per l'esistenza e l'unicità di una fattorizzazione LU} \begin{definition} Si dicono \textbf{sottomatrici principali di testa} le sottomatrici di $A$ di cui prendiamo le prime $k$ righe e colonne. Quando si utilizza l'aggettivo \textit{proprie}, si esclude $A$ stessa. \end{definition} Vale la seguente condizione sufficiente per l'esistenza e l'unicità di una fattorizzazione LU: \begin{proposition} Se tutte le sottomatrici principali di testa proprie sono non singolari allora esiste ed è unica la fattorizzazione $LU$ di $A$. \end{proposition} Se non sono verificate le ipotesi può comunque esistere una fattorizzazione $LU$. Per esempio: \[ \begin{pmatrix} 0 & 1 \\ 0 & 2 \end{pmatrix}=\begin{pmatrix} 1 & 0 \\ 1 & 1 \end{pmatrix}\begin{pmatrix} 0 & 1 \\ 0 & 1 \end{pmatrix}. \] Se $A$ è invertibile, allora la non singolarità delle sottomatrici principali di testa diventa una condizione necessaria oltre che sufficiente per l'esistenza e l'unicità della fattorizzazione LU. \smallskip Le seguenti classi di matrici ammettono sempre un'unica fattorizzazione LU: \begin{itemize} \item Matrici fortemente dominanti diagonali, le cui sottomatrici principali di testa sono fortemente dominanti diagonali e dunque invertibili, \item Matrici hermitiane definite positive (criterio di Sylvester), \item Matrici hermitiane definite positive su $\Span(e_1, \ldots, e_{n-1})$ e semidefinite positive su $\CC^n$ (metodo di Jacobi, criterio di Sylvester). \end{itemize} \subsection{Matrici elementari} Dati $\sigma\in\CC$, $u$, $v\in\CC^n$ si dice \textbf{matrice elementare} una matrice $M$ del tipo: \[M=I-\sigma uv^H.\] Si osserva che $uv^H$ è della seguente forma: \[ uv^H = (u_i \overline{v_j})_{ij}. \] Inoltre $\rg(uv^H) \leq \rg(u) \leq 1$. Se $u$ o $v$ sono nulli $uv^H$ è necessariamente nullo. \underline{Non} si deve confondere $uv^H$ con $u^H v$, che invece è il prodotto hermitiano complesso computato su $u$ e $v$. \medskip Un vettore $x$ è autovettore di $M$ se: \begin{itemize} \item $x=u \implies Mu=(1-\sigma(v^Hu))u$, e dunque $u$ è relativo all'autovalore $1-\sigma(v^Hx)$; \item $x$ t.c.~$v^Hx=0$ ($x$ è ortogonale a $v$) $\implies Mx=x$, il cui relativo autovalore è $1$. \end{itemize} Supponiamo che $\sigma(v^H u)$ sia diverso da $0$. La traccia di $M$ è: \[ \tr(M)=\tr(I)-\sigma \cdot \tr(uv^H))=(n-1) + (1-\sigma v^Hu). \] Dal momento che $v^\perp$ ha dimensione $n-1$ -- il prodotto hermitiano è positivo definito --, allora $\mu_{g}(1) \geq n-1$. \smallskip Osservando allora che $\tr(M)$ è la somma degli autovalori di $M$ e che $\mu_{g}(1-\sigma(v^H u)) \geq 1$, si conclude che gli unici autovalori di $M$ sono proprio $1$, con molteplicità algebrica e geometrica $n-1$, e $1-\sigma(v^H u)$, con molteplicità $1$. In particolare $M$ è sempre diagonalizzabile e vale la seguente proposizione: \begin{proposition} Se $\sigma(v^H u) \neq 0$, gli unici autovettori di $M$ sono i vettori ortogonali a $v$, che sono punti fissi, e i multipli di $u$. \end{proposition} $M$ è non singolare se e solo se $0$ non è autovalore, ossia se e solo se: \[ 1-\sigma(v^H u) \neq 0 \iff \sigma(v^H u) \neq 1. \] Se $M$ è non singolare, allora vale che: \[ M\inv=I-\tau uv^H\;\con \tau=\frac{-\sigma}{1-\sigma v^Hu}, \] e dunque anche $M\inv$ è una matrice elementare. \begin{proposition} Per ogni $x$ e $y \in \CC^n\setminuszero$ esiste $M$ matrice elementare con $\det M \neq0$ tale che $Mx=y$. In particolare è sufficiente scegliere $v$ non ortogonale sia a $x$ che a $y$ e porre $\sigma u=\frac{x-y}{v^Hx}$. \end{proposition} \subsubsection{Matrici elementari di Gauss e applicazione alla fattorizzazione LU} \begin{definition} Dato $x\in\CC^n$ con $x_1\neq0$ si definisce la relativa \textbf{matrice elementare di Gauss} come la matrice elementare: \[ M=I-ue_1^T, \quad u^T=\frac{1}{x_1}(0, x_2, \ldots, x_n). \] \end{definition} Vale sempre $Mx=x_1 \, e_1$. Inoltre una matrice elementare di Gauss è sempre invertibile, e $M^{-1}=I+ue_1^T$. Vale sempre $\det(M) = 1$. Una matrice elementare di Gauss corrisponde a uno step di annichilimento degli elementi sotto il \textit{pivot} dell'algoritmo di Gauss. \smallskip In particolare: \[ M=\begin{pmatrix} 1 & 0 & \dots & 0 \\ -u_2 & 1 & 0 \\ \vdots & 0 &\ddots \\ -u_n & 0 & \dots & 1 \end{pmatrix}\E M\inv=\begin{pmatrix} 1 & 0 & \dots & 0 \\ u_2 & 1 & 0 \\ \vdots & 0 &\ddots \\ u_n & 0 & \dots & 1 \end{pmatrix}. \] \begin{algorithm}[Fattorizzazione LU con le matrici di Gauss] Sia $A\in \CC^{n \times n}$ che soddisfa l'esistenza e unicità della fattorizzazione $LU$. Se $M_1$ è la matrice elementare di Gauss associata alla prima colonna di $A$, allora: \[ M_1A=\begin{pmatrix} a_{11} & \dots & a_{1n} \\ 0 \\ \vdots & {A_1} \\ 0 \end{pmatrix}. \] Se $A_1=(a_{ij}^{(1)})_{i,j>1}$, allora $a_{ij}^{(1)}=a_{ij}-\frac{a_{i1}}{a_{11}} a_{1j}$. Se $\hat{M}_2$ è la matrice elementare di Gauss che annulla la prima colonna di $A_1$, si definisce $M_2$ in modo tale che: \[ M_2=\begin{pmatrix} 1 & 0 \\ 0 & \hat{M}_2 \end{pmatrix}. \] Iterando questo procedimento si ottiene il seguente prodotto: \[ M_{n-1}\dots M_1 A=U, \quad M_1\inv\dots M_{n-1}\inv=L, \] dove $L$ è triangolare inferiore con $1$ sulla diagonale e $U$ è triangolare superiore. \smallskip In particolare vale che $\diag(U)=(a_{11},a_{22}^{(1)},\dots,a_{nn}^{(n-1)})$ e la $j$-esima colonna di $L=(\ell_{ij})$ è la $j$-esima colonna di $M_j$ cambiata di segno. Sono espresse di seguite le relazioni di ricorrenza: \[ \begin{cases} a_{ij}^{(k+1)}=a_{ij}^{(k)}+m_{ij}^{(k)}a_{jk}^{(k)} & i,j=k+1,\dots,n \\ m_{ik}^{(k)}=-\frac{a_{ik}^{(k)}}{a_{kk}^{(k)}} & i=k+1,\dots,n \\ \ell_{ij} = -m_{ij}^{(j)} & i \geq j \end{cases} \] Nel caso della risoluzione di un sistema lineare $Ax=b$, l'algoritmo può essere esteso aggiornando $b$ ad ogni step; in particolare $b^{(k+1)} = M_k b^{(k)}$ dove $M_k$ è la $k$-esima matrice di Gauss. In tal caso vale che: \[ b_i^{(k+1)} = b_i^{(k)} + m_{ik}^{(k)} b_k^{(k)}, \quad i = k+1, \ldots, n. \] Il costo totale dell'algoritmo è di $\frac{2}{3}n^3+O(n^2)$ operazioni aritmetiche. Se $A$ è anche tridiagonale il costo è $O(n)$. \end{algorithm} Riguardo al precedente algoritmo vale la seguente proposizione: \begin{proposition} Sia $A \in \CC^{n \times n}$ una matrice con sottomatrici principali proprie di testa non singolari per cui esiste ed è unica la fattorizzazione $A = LU$. Se $\tilde{L}$ e $\tilde{U}$ i valori effettivamente calcolati di $L$ e $U$ tramite il precedente algoritmo in aritmetica di macchina e sia $\Delta_A=A-\tilde{L}\tilde{U}$. Allora vale elemento per elemento di $\Delta_A$ la seguente disuguaglianza: \[ \abs{\Delta_A}\le 2nu\left(\abs{A}+|\Tilde{L}|\,|\Tilde{U}|\right)+O(u^2), \] dove $|T| := (|t_{ij}|)$. \smallskip Se $\tilde{y}$ è la soluzione del sistema $\tilde{L}y=b$ calcolato realmente in aritmetica di macchina mediante l'algoritmo di sostituzione in avanti e $\tx$ è il vettore effettivamente calcolato risolvendo $\tilde{U}x = \tilde{y}$ mediante sostituzione all'indietro, allora vale che: \[ (A + \hat{\Delta}_A) \tx = b, \] con \[ |\hat{\Delta}_A| \leq 4nu\left(\abs{A}+|\Tilde{L}|\,|\Tilde{U}|\right)+O(u^2), \] dove $|T| := (|t_{ij}|)$. \end{proposition} In particolare questa proposizione mostra che i valori $\tilde{L}$ e $\tilde{U}$ formano una decomposizione $LU$ di una perturbazione di $A$, e dunque è possibile effettuare un'analisi all'indietro dell'algoritmo di Gauss. Se i valori in modulo di $A$ sono troppo grandi, ci si può aspettare un mal funzionamento in senso numerico del metodo di eliminazione gaussiana (anche per il calcolo della soluzione $\tx$). \subsubsection{Strategia di \textit{pivoting} parziale per la decomposizione PLU} Permutando i pivot è sempre possibile fornire una fattorizzazione del tipo $PLU$. La \textbf{strategia del \textit{pivoting parziale}} (\textit{massimo pivot parziale}) consiste nel scegliere al passo $k$ come pivot il termine $a^{(k)}_{hk}$ con $h\ge k$ tale che $|a_{hk}^{(k)}|\ge|a_{ik}^{(k)}|$ per $i=k,\dots,n$. In questo modo se $a^{(k)}_{hk}=0$, allora tutta la parte di colonna è nulla e si può procedere allo step successivo dell'algoritmo; altrimenti si può scambiare la riga $h$-esima con la riga $k$-esima e applicare poi l'algoritmo di Gauss. \smallskip In questo modo si sta moltiplicando per una matrice di permutazione $P$ relativa alla trasposizione $\tau=(h, k)$. Nel caso in cui $a^{(k)}_{hk}\neq 0$ vale allora che: \[ A_{k+1}=M_k(P_kA_k). \] In particolare vale sempre $m_{ij}^{k} \leq 1$ per $i\geq j$. \smallskip Si osserva inoltre che: \[ M_k P_k = P_k \tilde{M_k}, \] dove, se $M_k = I - u e_q^T$, allora $\tilde{M_k} = I - P_k u e_q^T$ (infatti $\tilde{M_k} = P_k^T M_k P_k$). Applicando questa strategia ad ogni passo si ottiene allora una fattorizzazione $PA=LU$ dove $P$ è un'opportuna matrice di permutazione, e dunque si ottiene una fattorizzazione $PLU$. \subsubsection{Matrici elementari di Householder e applicazioni alla fattorizzazione QR} \begin{definition} Si definiscono \textbf{matrici elementari di Householder} le matrici elementari della forma: \[ M=I-\beta uu^H, \] con $u\in\CC^n$ e $\beta\in\RR$. \end{definition} Se $u \neq 0$ e $\beta = 0$ o $\beta = 2/(u^Hu) = 2/\langle u, u \rangle$, allora una matrice di Householder è unitaria. Inoltre una matrice di Householder è sempre hermitiana. \begin{proposition} Per ogni $x\in\CC^n\setminuszero$ esiste $M$ matrice elementare di Householder tale per cui $Mx=\alpha e_1$ per un dato $\alpha = \theta \norm{x}_2$, dove: \[ \theta=\begin{cases} \pm 1 & \text{se } x_1=0, \\ \pm \frac{x_1}{|x_1|} & \text{se } x_1\neq0. \end{cases}\] Tale matrice $M$ si ottiene ponendo $u^t=x-\alpha e_1 = (x_1-\alpha,x_2,\dots,x_n)$ e $\beta=\frac{2}{u^Hu}$. \end{proposition} Per evitare le cancellazioni nell'implementazione del calcolo di una matrice di Householder, e dunque migliorare la stabilità numerica, è consigliato scegliere sempre $\theta=-\frac{x_1}{\abs{x_1}}$ per $x_1 \neq 0$ (nel caso in cui $x_1 = 0$, la scelta è indiffernete). \medskip \begin{algorithm} Applicando la stessa filosofia dell'algoritmo di Gauss si possono utilizzare le matrici di Householder per calcolare la fattorizzazione $QR$ di una matrice $A$. Sia $u^{(k)}$ il vettore relativo alla $k$-esima matrice di Householder. Allora vale che: \[ u_i^{(k)} = \begin{cases} 0 & \se i < k, \\ a_{kk}^{(k)}\left(1 + \frac{\sqrt{\sum_{i=k}^n \abs{a_{ik}^{(k)}}^2}}{\abs{a_{kk}^{(k)}}}\right) & \se i = k, \\ a_{ik}^{(k)} & \altrimenti. \end{cases} \] mentre il parametro $\beta^{(k)}$ della stessa matrice è così dato: \[ \beta^{(k)} = 2\bigg/\sum_{i=k}^n \abs{u_i^{(k)}}^2 \] A partire da questi, si ottengono i termini di $A_{k+1}$: \[ a_{i, j}^{(k+1)} = \begin{cases} a_{ij}^{(k)} & \se j < k, \se i \leq k, \\ 0 & \se j = k, \E i > k, \\ a_{ij}^{(k)} - \beta^{(k)} u_i^{(k)} \sum_{r=k}^n \overline{u}_r^{(k)} a_{rj}^{(k)} & \se i \geq k, j > k. \end{cases} \] Nel caso della risoluzione di un sistema lineare $Ax=b$, l'algoritmo può essere esteso aggiornando $b$ ad ogni step; in particolare $b^{(k+1)} = M_k b^{(k)}$ dove $M_k$ è la $k$-esima matrice di Householder. In tal caso vale che: \[ b_i^{(k+1)} = b_i^{(k)} - \beta^{(k)} u_i^{(k)} \sum_{r=k}^n \overline{u}_r^{(k)} b_r^{(k)}, \quad i=k, \ldots, n. \] Il costo totale è di $\frac{4}{3}n^3 + O(n^2)$ operazioni, il \textit{doppio} di quello dell'algoritmo gaussiano. \end{algorithm} Riguardo a quest'ultimo algoritmo vale la seguente proposizione: \begin{proposition} Sia $A \in \CC^{n \times n}$ e sia $\tilde{R}$ la matrice triangolare superiore ottenuta in aritmetica di macchina applicando il metodo di Householder con l'algoritmo dato in precedenza. Sia inoltre $\tx$ la soluzione ottenuta risolvendo in aritmetica di macchina il sistema $Ax=b$ con il metodo di Householder, aggiornando le entrate di $b$. Allora $\tx$ risolve il sistema $(A + \Delta_A) \tx = b + \delta_b$, dove: \[ \norm{\Delta_A}_F \leq u(\gamma n^2 \norm{A}_F + n \norm{\tilde{R}}_F) + O(u^2), \] \[ \norm{\delta_b}_2 \leq \gamma n^2 u \norm{b}_2 + O(u^2), \] dove $\gamma$ è una costante positiva. Inoltre $\norm{\tilde{R}}_F$ differisce da $\norm{\tilde{R}}_F = \norm{\tilde{A}}_F$ ($A = QR$) per un termine proporzionale a $u$, e dunque: \[ \norm{\Delta_A}_F \leq \gamma u(n^2 + n) \norm{A}_F + O(u^2). \] \end{proposition} In particolare, il metodo di Householder è numericamente stabile. \section{Metodi stazionari iterativi per sistemi lineari} \begin{definition} Dato il sistema lineare $Ax=b$ con $A\in \CC^{n \times n}$ si definisce un generico \textbf{partizionamento additivo di $A$} una partizione della forma $A=M-N$ con $\det M\neq0$. \end{definition} A partire da ciò possiamo riscrivere equivalentemente il sistema come $Mx=Nx+b$, il quale riconduce alla seguente scrittura equivalente del sistema originale: \[ x=Px+q, \quad P=M\inv N, \, q=M\inv b. \] Questa formulazione del sistema viene detta \textbf{problema di punto fisso}. Una volta fissato $x^{(0)}\in\CC^n$ si può generare in modo naturale una successione di vettori: \[ x^{(k+1)}=Px^{(k)}+q. \] Se tale successione ammette limite $x^*$, allora il punto limite $x^*$ è un punto fisso di $Px+q$, ed è dunque l'unica soluzione di $Ax=b$. Questo metodo è detto \textbf{metodo iterativo stazionario}. La matrice $P$ ed il vettore $q$ sono indipendenti da $k$, e quindi la matrice $P$ viene detta \textbf{matrice di iterazione} del metodo. \medskip Si dice che il metodo iterativo è \textbf{convergente} se per ogni scelta del vettore $x^{(0)}$ la successione $x^{(k)}$ converge alla soluzione $x^*$ del sistema. \smallskip Se $A$ è invertibile, definiamo la successione $e^{(k)}=x^{(k)}-x^*$ dove $x^*$ è soluzione del sistema. Allora $e^{(k+1)} = Px^{(k)} + Q - (Px^* + Q) = P e^{(k)}$. A partire da questa osservazione si dimostra il seguente teorema: \begin{theorem} Se esiste una norma di matrice indotta $\norm{\cdot}$ tale che $\norm{P}<1$, allora la matrice $A$ è invertibile. In tal caso il metodo iterativo è convergente ($\lim_{k \to \infty} \norm{e^{(k)}} \leq \lim_{k \to \infty} \norm{P}^k \norm{e^{(0)}} = 0$). \end{theorem} \begin{theorem} Il metodo iterativo è convergente e $\det A \neq 0$ se e solo se $\rho(P)<1$. \end{theorem} Il quoziente $\norm{e^{(k)}}/\norm{e^{(k-1)}}$ rappresenta la riduzione dell'errore al passo $k$-esimo del metodo iterativo rispetto alla norma scelta. Se si considera la media geometrica $\theta_k\left(e^{(0)}\right)$ delle prime $k$ riduzioni, si ricava che: \[ \theta_k\left(e^{(0)}\right)=\left(\frac{\norm{e^{(1)}}}{\norm{e^{(0)}}}\dots\frac{\norm{e^{(k)}}}{\norm{e^{(k-1)}}}\right)^{\frac{1}{k}}=\left(\frac{\norm{P^ke^{(0)}}}{\norm{e^{(0)}}}\right)^{\frac{1}{k}}\le\norm{P^k}^{\frac{1}{k}}. \] \begin{definition} Si definisce \textbf{riduzione asintotica media per passo} di un metodo iterativo con errore iniziale $e^{(0)}$ il valore: \[ \theta(e^{(0)})=\lim_{k \to \infty}\theta_k(e^{(0)}). \] Tale riduzione rappresenta la velocità di convergenza del metodo: più è piccolo il valore, più veloce è il metodo. \end{definition} In particolare vale che: \[ \theta(e^{(0)}) \leq \lim_{k \to \infty}\norm{P^k}^{\frac{1}{k}}=\rho(P). \] L'uguaglianza è raggiunta per $e^{(0)}$ pari ad un autovettore relativo al raggio spettrale di $P$. Se $P$ è nilpotente, allora il metodo converge in un numero finito di passi. \medskip Un esempio di metodo iterativo è dato dal \textbf{metodo di Richardson}, che pone $P = I - \alpha A$, $q = M\inv b = \alpha A$. In particolare tale metodo genera la successione $x^{(k+1)}=x^{(k)}-\alpha(Ax^{(k)}-b)$, dove $\alpha \neq 0$ è un opportuno parametro e $Ax_n - b$ è detto \textit{errore residuo}. \smallskip Poiché il metodo converge se $\rho(P) < 1$, l'idea naturale è quella di scegliere $\alpha$ in modo tale che $\rho(P)$ sia minimo, ricordandosi che $\lambda$ è autovalore di $P$ se e solo se $\alpha (\lambda + 1)$ è autovalore di $A$. Infatti: \[ \det(\lambda I- P) = \det(\lambda I + I - \alpha A) = \frac{1}{\alpha^n} \det(\alpha(\lambda + 1)I - A).\] Se la matrice è hermitiana e definita positiva, allora $\alpha=1/\norm{A}$, dove $\norm{\cdot}$ è una norma indotta, garantisce la convergenza del metodo. \subsubsection{I metodi di Jacobi e Gauss-Seidel} Sia $A$ tale per cui $a_{ii} \neq 0$. Si può allora decomporre $A$ nel seguente modo: \[ A=D-B-C, \] dove: \begin{itemize} \item $D$ è diagonale con $d_{ii}=a_{ii}\neq 0$ ($D = \diag(A)$, $D$ invertibile); \item $B$ è strettamente triangolare inferiore con $b_{ij}=-a_{ij}$ con $i>j$ ($B = -\tril(A)$); \item $C$ è strettamente triangolare superiore con $c_{ij}=-a_{ij}$ con $ij\E C$ la matrice triangolare superiore a blocchi tale che $C_{ij}=-A_{ij}\;\con i 0$ (e dunque $f(a) < 0$); altrimenti è sufficiente applicare l'algoritmo al contrario a $-f(x)$. \smallskip Al passo $(k+1)$-esimo si considera $c_k = (a_k + b_k)/2$, dove $a_0 = a$ e $b_0 = b$. Se $f(c_k) = 0$, l'algoritmo termina; altrimenti, se $f(c_k) > 0$, $b_{k+1} = c_k$ e $a_{k+1} = a_k$, e se $f(c_k) < 0$, $b_{k+1} = b_k$ e $a_{k+1} = c_k$. \smallskip Dal momento che gli intervalli $I_k = [a_k, b_k]$ hanno lunghezza $(b-a)/2^k$, per $k \to \infty$ la successione $x_k$ ha limite $\alpha$, che è tale per cui $f(\alpha) = 0$. Inoltre l'intervallo si dimezza a ogni step, e dunque il metodo converge esponenzialmente (non per questo è veloce, anzi: come vedremo, ci sono metodi estremamenti più veloci che convergono in modo \textit{doppiamente esponenziale}). \subsection{Metodi del punto fisso} I metodi del punto fisso si ottengono trasformando il problema $f(x)=0$ in un problema del tipo $g(x)=x$. Questa trasformazione si può ottenere in infiniti modi diversi. Ad esempio, data una qualsiasi funzione $h(x)$ (diversa da zero nei punti del dominio di $f(x)$), si può porre: \[ g(x)=x-\frac{f(x)}{h(x)}, \] i cui punti fissi sono gli zeri di $f(x)$. \smallskip Se $g(x)$ è continua, la successione $x_{k+1} = g(x_k)$, con $x_0 \in \RR$, se convergente, converge ad un punto fisso di $g(x)$ (e dunque ad uno zero di $f(x)$). Questo metodo è detto \textbf{metodo del punto fisso} (o \textit{di iterazione funzionale}) associato a $g(x)$. \begin{theorem}[del punto fisso] Sia $\mathcal{I}=[\alpha-\rho,\alpha+\rho]$ e $g(x)\in C^1(\mathcal{I})$ dove $\alpha=g(\alpha)$ e $\rho>0$. Si denoti con $\lambda=\max_{x \in \mathcal{I}}\lvert g'(x) \rvert$. Se $\lambda<1$ allora per ogni $x_0\in\mathcal{I}$, posto $x_{k+1}=g(x_k)$, vale che: \[ \lvert x_k-\alpha \rvert \le\lambda^k\rho. \] Pertanto $x_k\in\mathcal{I}$ e $\lim_k x_k=\alpha$. Inoltre $\alpha$ è l'unico punto fisso di $g(x)$ in $\mathcal{I}$. \end{theorem} Supponiamo di avere un intervallo $[a,b]$ in cui è presente un punto fisso $\alpha$ e che $\lvert g'(x) \rvert<1$ su tale intervallo. Sotto queste ipotesi siamo certi che almeno con una delle due scelte $x_0=a$ o $x_0=b$ la successione è ben definita e converge al punto fisso $\alpha$. Infatti basta prendere rispettivamente $\rho=\alpha-a$ o $\rho=b-\alpha$ ed applicare il teorema. Se la successione $x_k$ cade fuori dall'intervallo, allora si arrestano le iterazioni e si riparte con $x_0$ uguale all'altro estremo. \smallskip In questo caso la convergenza vale in un intorno opportuno del punto fisso $\alpha$. Denotiamo questo fatto con l'espressione \textbf{convergenza locale}. Si parla di \textbf{convergenza globale} qualora la convergenza vi sia per ogni scelta iniziale. \smallskip Nell'aritmetica a virgola mobile il teorema diventa: \begin{theorem} Nelle ipotesi del teorema del punto fisso sia: \[ \Tilde{x}_{k+1}=g(\Tilde{x}_k)+\delta_k, \] dove $\lvert \delta_k \rvert \le\delta$ è l'errore commesso nel calcolo di $g(\Tilde{x}_k)$ e $\delta$ è una quantità nota. Posto $\sigma=\delta/(1-\lambda),\se\sigma<\rho$ allora: \[ \lvert\Tilde{x}_k-\alpha\rvert \le(\rho-\sigma)\lambda^k+\sigma. \] \end{theorem} Questo teorema ci dice che la distanza di $\Tilde{x}_k$ da $\alpha$ è limitata dalla somma di due parti. La prima converge a zero in modo esponenziale su base $\lambda$. La seconda è costante e rappresenta l'intervallo di incertezza sotto il quale non è consentito andare. \smallskip Per la funzione $g(x)=x-f(x)$ l'intervallo di incertezza è: \[ I=\left[\alpha-\frac{\delta}{\lvert f'(\alpha) \rvert},\alpha+\frac{\delta}{\lvert f'(\alpha) \rvert}\right]. \] Se $g'(x)>0$, allora $I\subseteq [\alpha-\sigma,\alpha+\sigma]$. \smallskip Si ricorda che: \[ x_{k+1}-\alpha=g'(\xi_k)(x_k-\alpha), \] dove $\xi_k \in (\alpha, x_k)$ viene dall'applicazione del teorema di Lagrange. Ciò implica che se $0\alpha\Rightarrow x_k>\alpha \; \forall k, \quad \alpha \alpha$ la sottosuccessione $\{x_{2k}\}$ cresce mentre $\{x_{2k+1}\}$ decresce (entrambe al punto $\alpha$). \subsubsection{Velocità di convergenza} \begin{definition} Sia $\{x_k\}$ una successione tale che $\lim_kx_k=\alpha$. Si ponga allora $e_k = x_k - \alpha$ come l'errore assoluto al passo $k$-esimo. Supponiamo esista il limite della riduzione dell'errore al passo $k$-esimo: \[ \gamma=\lim_k \abs{\frac{x_{k+1}-\alpha}{x_k-\alpha}} = \lim_k \abs{\frac{e_{k+1}}{e_k}}. \] La convergenza di $\{x_k\}$ a $\alpha$ è detta allora: \begin{itemize} \item lineare (o geometrica) se $0<\gamma<1$; \item sublineare se $\gamma=1$; \item superlineare se $\gamma=0$. \end{itemize} \end{definition} Nel caso di convergenza superlineare, se $p>1$ ed esiste il limite: \[ \lim_k \abs{\frac{x_{k+1}-\alpha}{(x_k-\alpha)^p}} =\sigma, \quad 0<\sigma<\infty \] si dice che la successione \textbf{converge con ordine $p$}. \begin{remark} L'ordine di convergenza non è obbligatoriamente intero. \end{remark} \begin{theorem} Sia $g(x)\in C^1([a,b])$ e $\alpha\in(a,b)$ t.c.~$g(\alpha)=\alpha$. Se esiste un $x_0\in[a,b]$ tale che la successione $x_{k+1}=g(x_k)$ converge linearmente ad $\alpha$ con fattore $\gamma$ allora: \[ \lvert g'(\alpha)\rvert=\gamma. \] Viceversa, se $0<\lvert g'(\alpha)\rvert<1$ allora esiste un intorno $I$ di $\alpha$ contenuto in $[a,b]$ tale che per ogni $x_0\in I$ la successione $\{x_k\}$ converge ad $\alpha$ in modo lineare con fattore $\gamma=\lvert g'(\alpha) \rvert$. \end{theorem} \begin{theorem} Sia $g(x)\in C^1([a,b])$ e $\alpha\in (a,b)$ t.c.~$g(\alpha)=\alpha$. Se esiste un $x_0\in[a,b]$ tale che la successione $x_{k+1}=g(x_k)$ converge sublinearmente ad $\alpha$ allora: \[ \abs{g'(\alpha)}=1. \] Viceversa, se $\abs{g'(\alpha)}=1$ esiste un intorno $I$ di $\alpha$ contenuto in $[a,b]$ tale che per ogni $x\in I$, $x\neq\alpha$ vale $\abs{g'(x)}<1$ e $g'(x)$ non cambia segno su $I$ allora tutte le successioni $\{x_k\}$ con $x_0\in I$ convergono ad $\alpha$ in modo sublineare. \end{theorem} \begin{theorem} Sia $g(x)\in C^p([a,b])$ con $p>1$ intero e $\alpha\in(a,b)$ t.c.~$g(\alpha)=\alpha$. Se esiste un $x_0\in[a,b]$ tale che la successione $x_{k+1}=g(x_k)$ converge superlinearmente ad $\alpha$ con ordine di convergenza $p$ allora: \[ \abs{g^{(k)}(\alpha)}=0 \;\;\;\, 1 \leq k \leq p-1, \quad g^{(p)}(\alpha)\neq0. \] Viceversa se $\abs{g^{(k)}(\alpha)}=0$ per $k=1,\dots,p-1\E g^{(p)}(\alpha)\neq0$ allora esiste un intorno $I$ di $\alpha$ tale che per ogni $x_0\in I$ tutte le successioni $\{x_k\}$ convergono ad $\alpha$ in modo superlineare con ordine $p$. \end{theorem} \begin{definition} La successione $\{x_k\}$ converge ad $\alpha$ \textbf{con ordine almeno $p$} se esiste una costante $\beta$ tale che \[ \abs{x_{k+1}-\alpha}\le\beta \, \abs{x_k-\alpha}^p. \] \end{definition} Se una successione converge con ordine $q\ge p$ allora converge anche con ordine almeno $p$. \smallskip Se una successione $x_k$ converge ad $\alpha$ in modo che l'errore relativo $\eps_k$ al passo $k$-esimo è limitato nel seguente modo: \[ \varepsilon_k=\abs{x_k-\alpha/\alpha}\le\beta\gamma^{p^k}\] allora il numero di cifre significative ($1+\log_2(\varepsilon_k\inv)$) è tale che \[ 1+\log_2 \, \varepsilon_k\inv \ge1+\log_2 \, \beta\inv+p^k\log_2 \,\gamma\inv. \] \subsubsection{Confronto tra metodi} Siano dati due metodi iterativi del punto fisso definiti da due funzioni $g_1(x)$ e $g_2(x)$. Si supponga che siano entrambi siano o a convergenza lineare o a convergenza superlineare. Denotiamo con $c_1$, $c_2$ il numero di flops per passo dei due metodi. Se siamo nel caso di convergenza lineare, detti $\gamma_1$, $\gamma_2$ i due fattori di convergenza allora il primo metodo risulta più conveniente se \[ \frac{c_1}{c_2}<\frac{\log{\gamma_1}}{\log{\gamma_2}}. \] Nel caso di convergenza superlineare, se $p_1$ e $p_2$ sono i due ordini di convergenza, il primo metodo è più efficiente del secondo se \[\frac{c_1}{c_2}<\frac{\log{p_1}}{\log{p_2}}.\] \subsection{Alcuni metodi del punto fisso} \subsubsection{Metodo delle secanti} Sia $f(x)\in C^1([a,b])\E\alpha\in[a,b]\tc f(\alpha)=0$. Il metodo definito dalla funzione $$g(x)=x-\frac{f(x)}{m}$$ dove $m$ è un'opportuna costante è detto \textbf{metodo delle costanti}. \smallskip Il metodo consiste nel tracciare la retta passante per il punto $(x_k,f(x_k))$ di coefficiente angolare $m$ e considerare come $x_{k+1}$ l'ascissa dell'intersezione di tale retta con l'asse delle ascisse, reiterando. \smallskip Una condizione sufficiente di convergenza è che valga $\abs{1-f'(x)/m}<1$ in un intorno di circolare di $\alpha$. È sufficiente quindi scegliere $m$ in modo che abbia lo stesso segno di $f'(x)$ e $\abs{m}>\frac{1}{2}\abs{f'(x)}$. Se $f'(\alpha)$ fosse nota, la scelta $m=f'(\alpha)$ darebbe una convergenza superlineare. \subsubsection{Metodo delle tangenti di Newton} La funzione $g(x)$ del \textbf{metodo di Newton} è definita nel modo seguente: \[ g(x)=x-\frac{f(x)}{f'(x)}. \] Per tale metodo vale: \[ g'(x)=1-\frac{f'(x)^2 - f''(x)f(x)}{f'(x)^2} = \frac{f''(x) f(x)}{f'(x)^2}. \] Il metodo consiste nel tracciare la retta passante per il punto $(x_k,f(x_k))$ e tangente a esso e considerare come $x_{k+1}$ l'ascissa dell'intersezione di tale retta con l'asse delle ascisse, reiterando. \begin{theorem} Sia $f(x)\in C^2([a,b])$ e sia $\alpha\in(a,b)$ t.c.~$f(\alpha)=0$. Se $f'(\alpha)\neq0$, allora esiste un intorno $I=[\alpha-\rho,\alpha+\rho]\subset[a,b]$ tale che per ogni $x_0\in I$ la successione generata dal metodo di Newton converge ad $\alpha$. Inoltre se $f''(\alpha)\neq0$ la convergenza è superlineare di ordine $2$, se $f''(\alpha)=0$ la convergenza è di ordine almeno $2$. \end{theorem} \begin{theorem} Sia $f(x)\in C^p([a,b])$ con $p>2$ e $\alpha\in(a,b)$ t.c.~$f(\alpha)=0$. Se $f'(\alpha)=\dots=f^{(p-1)}(\alpha)=0$, $f^{(p)}(\alpha)\neq0$ e $f'(x)\neq0$ per $x\neq\alpha$, allora esiste un intorno $I=[\alpha-\rho,\alpha+\rho]\subset[a,b]$ in cui $f'(x)\neq0$ per $x\in I,x\neq\alpha$ e tale che per ogni $x_0\in I$, la successione generata dal metodo di Newton converge ad $\alpha$. Inoltre $\alpha$ è l'unico zero di $f(x)$ in $I$. La convergenza è lineare con fattore di convergenza $\gamma = 1-1/p$. \end{theorem} \begin{theorem} Se la funzione $f(x)$ è di classe $C^2$ sull'intervallo $I=[\alpha,\alpha+\rho]$ ed è tale che $f'(x)f''(x)>0$ per $x\in I$ allora per ogni $x_0\in I$ la successione generata dal metodo di Newton applicato ad $f(x)$ converge decrescendo ad $\alpha$ con ordine $2$. Un risultato analogo vale su intervalli del tipo $[\alpha-\rho,\alpha]$, su cui invece la successione cresce. \end{theorem} \begin{example}[Calcolo del reciproco di un numero $a$] Per calcolare $1/a$ si può porre $f(x)=a-1/x$ e poi applicare il metodo di Newton. \end{example} \begin{example}[Calcolo della radice $p$-esima di un numero $a$] Per calcolare $\sqrt[p]{a}$ è sufficiente applicare il metodo di Newton con $f(x) = (x^p-a)x^{-q}$, dove $q\geq 0$. \end{example} \begin{example}[Metodo di Aberth] Se si hanno delle approssimazione $t_1$, ..., $t_n$ di zeri di $p(x)$, si può applicare il metodo di Newton su $f(x) = p(x) \big/ \prod_{i=1}^n (x-t_i)$, in modo tale che $f(t_i)$ sia sempre più grande, diminuendo gli errori e migliorando le approssimazioni. \end{example} \section{Interpolazione di funzioni} Siano $\varphi_0(x)$, ..., $\varphi_n(x) : [a,b] \to \RR$ funzioni linearmente indipendenti. Siano $(x_i,y_i)$ $n+1$ valori assegnati (\textit{nodi}) tali che $x_i\in[a,b]$ e $x_i\neq x_j$ se $i\neq j$. \smallskip Il problema dell'interpolazione consiste nel determinare i coefficienti $a_i \in \RR$ tali per cui: \[ f(x)=\sum_{i=0}^na_i \, \varphi_i(x), \qquad f(x_i)=y_i\;\forall i. \] Tali condizioni sono dette \textbf{condizioni di interpolazione}. \subsection{Interpolazione polinomiale} L'interpolazione polinomiale ricerca il polinomio $p$ di grado $n+2$ che soddisfi le condizioni di interpolazione di $n+1$ nodi. Sia pertanto $\varphi_i(x)=x^i$. Siano $(x_i, y_i)$ per $i = 0$, ..., $n+1$ i nodi dell'interpolazione. \subsubsection{Interpolazione monomiale} \begin{definition} Si dice \textbf{matrice di Vandermonde} nei nodi $x_0$, ..., $x_n$ la matrice di taglia $n+2$: \[ V_n=\begin{pmatrix} 1 & x_0 & \dots & x_0^n \\ 1 & x_1 & \dots & x_1^n \\ \vdots & \vdots & &\vdots \\ 1 & x_n & \dots & x_n^n \end{pmatrix} \] \end{definition} \begin{theorem} Vale $\det V_n=\prod_{j 0$, ovverosia: \[ \Pi_n = \begin{pmatrix} 1 & 0 \\ 0 & P_{n-1} \end{pmatrix}, \] dove $P_{n-1}$ è la matrice identità specchiata orizzontalmente (ossia con soli $1$ sull'antidiagonale), che è tale per cui $P_{n-1}^2 = I$. \end{itemize} \end{proposition} Come corollario della precedente proposizione, si ottiene che $\Omega_n^{-1} = \frac{1}{n} \Omega_n^H$ e che $\Omega_n^{-1} = \frac{1}{n} \Pi_n \Omega_n = \frac{1}{n} \Omega_n \Pi_n$. Dunque il problema dell'interpolazione si risolve ponendo $z= \Omega_n^{-1} y$ e applicando una delle formule proposte precedentemente. \begin{remark} La matrice $F_n=\Omega_n/\sqrt{n}$ è unitaria, e dunque $\norm{F_n}_2=1$. Ciò implica che $\norm{\Omega_n}_2=\sqrt{n}$, dunque il numero di condizionamento di $\Omega_n$ è $1$, ovverosia il problema dell'interpolazione ai nodi di Fourier è ben condizionato. \end{remark} \subsection{Calcolo di (I)DFT: \textit{Fast Fourier Transform} (FFT)} Consideriamo il caso del calcolo di $\IDFT(y)$, dove abbiamo i coefficienti $z_0$, ..., $z_{n-1}$ e vogliamo trovare i valori $y_i$ tali per cui $\Omega_n (z_i) = (y_i)$. Il calcolo di $\DFT(z)$ si può poi fare eseguendo $\IDFT(z)$ e applicando le relazioni introdotte nella proposizione precedente, aggiungendo $n$ divisioni e permutando gli indici. \medskip Consideriamo il caso in cui $n=2^q$. Allora: \[ y_i = \sum_{j=0}^{n-1} \omega_n^{ij} z_j. \] Ricordando che $\omega_n^2 = \omega_{\frac{n}{2}}$, allora, separando gli indici pari da quelli dispari, vale che: \[ y_i =\sum_{j=0}^{\frac{n}{2}-1}\omega_{\frac{n}{2}}^{ij}z_{2j}+\omega_n^i \sum_{j=0}^{\frac{n}{2}-1}\omega_{\frac{n}{2}}^{ij}z_{2j+1}. \] Pertanto, detti $Y = \left(y_0,\dots,y_{\frac{n}{2}-1}\right)^\top$, $Y' = \left(y_{\frac{n}{2}},\dots,y_{n-1}\right)^\top$ e $D_n = \diag\left(1,\omega_n,\dots,\omega_n^{\frac{n}{2}-1}\right)$, valgono le seguenti due identità: \[ Y=\IDFT_{\frac{n}{2}}(z_{\text{pari}})+D_n \, \IDFT_{\frac{n}{2}}(z_\text{disp}), \] \[ Y'=\IDFT_{\frac{n}{2}}(z_\text{pari})-D_n \, \IDFT_{\frac{n}{2}}(z_\text{disp}). \] In forma matriciale le due identità si scrivono infine come: \[ y=\begin{pmatrix} \Omega_{\frac{n}{2}}z_\text{pari} \\ \Omega_{\frac{n}{2}}z_\text{pari} \end{pmatrix}+\begin{pmatrix} D_n\Omega_{\frac{n}{2}}z_\text{disp} \\ -D_n\Omega_{\frac{n}{2}}z_\text{disp} \end{pmatrix}. \] Poiché $n$ è potenza di $2$ possiamo ripetere questa strategia per calcolare le due trasformate di ordine $\frac{n}{2}$ mediante quattro trasformate di ordine $\frac{n}{4}$, ecc... fino ad avere trasformate di ordine $1$ che non richiedono alcuna operazione. \medskip Il costo $c(n)$ di flops di una IDFT su $n$ nodi con questo metodo è ricorsivamente: \[ c(n)=2 \, c\!\left(\frac{n}{2}\right)+\frac{n}{2}+n=2 \, c\!\left(\frac{n}{2}\right)+\frac{3}{2}n, \] dove $n/2$ sono le moltiplicazioni effettuate e $n$ le addizioni. \smallskip Dacché $c(1)=0$ e $n=2^q$ vale che: \[ c(n)=\frac{3}{2}n\log_2n=3q2^{q-1}. \] In generale, l'algoritmo ha dunque complessità $O(n \log_2(n))$. L'algoritmo per il calcolo della IDFT che si ottiene in questo modo è noto come \textbf{algoritmo di Cooley-Tukey}. \subsubsection{Applicazioni della FFT: calcolo del prodotto di polinomi} Dati $a(t)$, $b(t)\in\CC[x]$ allora $c(t)=a(t)b(t)$ si può calcolare trovando i coefficienti $c_i=\sum_{j+k=i}a_jb_k$. Se $a(t)$, $b(t)$ hanno grado rispettivamente $p$, $q$, questo algoritmo ha complessità $O((p+q)^2)$. \medskip Sia $n = 2^q \geq p+q$. Con la FFT posso valutare $a(t)$, $b(t)$ nelle radici $n$-esime dell'unità. Dopo aver moltiplicato i valori ottenuti, si può applicare una DFT, ricavando i coefficienti di $c(t)$. Questo algoritmo utilizza $2$ IDFT, $n$ moltiplicazioni ed una DFT; pertanto ha un costo di $O(n\log_2(n))$ flops. \section{Metodi dell'integrazione approssimata} In questa sezione si illustrano alcuni metodi per approssimare gli integrali su un intervallo. \begin{definition} Data $f : [a, b] \to \RR$ continua e $n+1$ nodi $x_0$, ..., $x_n$ si definiscono le due quantità: \[ S[f]=\int_a^bf(x)\,dx, \quad S_{n+1}[f]=\sum_{i=0}^nw_if(x_i), \] dove i termini $w_i$ sono positivi reali detti \textit{pesi}, possibilmente variabili. Una scelta di pesi corrisponde a una \textbf{formula di integrazione approssimata}. \end{definition} \begin{definition} Si dice \textbf{resto} il valore \[ r_{n+1}=S[f]-S_{n+1}[f]. \] \end{definition} \begin{definition} Si dice che una formula di integrazione approssimata ha \textbf{grado di precisione} (massimo) $k\ge0$ se vale \[ r_{n+1}=0 \impliedby f(x)=x^{j}, \quad 0 \leq j \leq k \] e se \[ r_{n+1}\neq0 \impliedby f(x)=x^{k+1}. \] \end{definition} \subsection{Formule di integrazione dell'interpolazione di Lagrange} Si può approssimare un integrale utilizzando l'interpolazione di Lagrange sui nodi, ponendo \[ S_{n+1}[f] = \int_a^b \tilde{f}(x)\,dx = \sum_{i=0}^n \left[ f(x_i) \int_a^b L_i(x)\,dx\right], \] dove $\tilde{f}(x)$ è il polinomio ottenuto applicando l'interpolazione di Lagrange su $f$ dati i nodi $(x_0, f(x_0))$, ..., $(x_n, f(x_n))$ e i pesi sono tali per cui $w_i = \int_a^b L_i(x)\,dx$. \smallskip Le formule di integrazione interpolatorie hanno grado di precisione almeno $n$ (infatti se $j \leq n$, $x^j$ è il polinomio che interpola i nodi dati) ed hanno grado di precisione massimo $2n+1$. Se una formula di integrazione ha grado di precisione almeno $n$, allora è interpolatoria. \smallskip Se $f\in C^{(n+1)}([a,b])$ e $\abs{f^{(n+1)}(x)}\le M$ allora \[ \abs{r_{n+1}}\leq\frac{M}{(n+1)!}\abs{\int_a^b\prod_{j=0}^n(x-x_j) \, dx}. \] \subsection{Formule di Newton-Cotes (semplici)} Scegliendo nodi equispaziati, con $h=\frac{b-a}{n}$, si ottengono le formule di Newton-Cotes. nel caso $n=1$ vale \[ w_0=\int_a^bL_0(x)\,dx=\frac{h}{2}=w_1\Rightarrow S_2[f]=\frac{h}{2}(f(x_0)+f(x_1)), \] ovverosia corrisponde all'area del trapezio con di base $\abs{x_0}$ e $\abs{x_1}$ e con altezza $h$. Il resto $r_2$ corrisponde invece a \[ \int_a^bf(x)\,dx-S_2[f]=-\frac{1}{12}h^3f''(\xi), \quad \xi\in(a,b). \] Nel caso $n=2$ invece \[ w_0=w_2=\frac{h}{3}, \quad w_1=\frac{4}{3}h \] ed il resto è dunque \[ \int_a^bf(x)\,dx-S_3[f]=-\frac{1}{90}h^5f^{(4)}(\xi),\quad \xi\in(a,b) \] \subsection{Formule di Newton-Cotes composte (trapezi, Cavalieri-Simpson)} Si può rendere più precisa la formula di integrazione introducendo altri nodi equispaziati $z_0$, ..., $z_N$ e risolvendo \[ \int_a^bf(x)\,dx=\sum_{i=0}^{N-1}\int_{z_i}^{z_{i+1}}f(x)\,dx, \] approssimando ciascuno degli integrali con le formule di Newton-Cotes semplici per $n=1$ o $2$. \smallskip Scegliendo sempre $n=1$ si ottiene la \textbf{formula dei trapezi}, per la quale: \[ J_2^{(n)}[f]=\sum_{i=0}^{n-1}S_2^{(i)}[f]=\frac{b-a}{2n}\left[f(z_0)+f(z_n)+2\sum_{i=1}^{n-1}f(z_i)\right] \] Inoltre vale che: \[ S[f]-J_2^{(nN)}[f]=-\frac{(b-a)h^2}{12}f''(\xi). \] Scegliendo invece sempre $n=2$ si ottiene la \textbf{formula di Cavalieri-Simpson}, per la quale: \begin{multline*} J_3^{(n)}[f] = \frac{b-a}{6n} \!\left[f(z_0) + f(z_n) + 2 \sum_{k=1}^{n-1} f(z_k) + \right. \\ \left.+ 4 \sum_{k=0}^{n-1} f\!\left(\frac{z_k + z_{k+1}}{2}\right) \right]. \end{multline*} \section{Risoluzione di problemi di Cauchy con metodi a un passo} Sia $I \subseteq \RR$ un intervallo e consideriamo una funzione $f(t, y) : I \times \RR \to \RR$ continua. Dati i parametri iniziali $t_0 \in I$, $y_0 \in \RR$ si associa a questi il seguente \textit{problema di Cauchy}: \[ \begin{cases} y'(t) = f(t, y(t)), \\ y(t_0) = y_0, \end{cases} \] su cui si ipotizza che l'incognita $y$ è $C^1(I, \RR)$. \smallskip Se $f(t, y)$ è Lipschitziana rispetto a $y$, ovverosia esiste $L \in \RR$ t.c.~$\abs{f(t, y_1) - f(t, y_2)} \leq L \abs{y_1 - y_2}$ per ogni $t \in I$, $y_1$, $y_2 \in \RR$, allora esiste ed è unica la soluzione $y$ al problema di Cauchy proposto. \begin{notation}[Intervallo destro di $t_0$ e successione per passi] Per intervallo destro di un $t_0 \in \RR$ intendiamo un intervallo della forma $[t_0, t_0 + T]$ con $t > 0$. \smallskip Dato $h > 0$, definiamo $t_n := t_0 + nh$ per $n = 0$, ..., $N_h$, dove $N_h$ è il massimo naturale $p$ per cui $ph \leq T$ (ossia il massimo naturale per cui la successione definita sta dentro l'intervallo destro definito). \end{notation} \begin{notation}[Approssimazioni] Scriveremo $y_n$ per indicare $y(t_n)$ e $u_n$ per indicare l'approssimazione ricavata di $y(t_n)$, mentre $f_n := f(t_n, u_n)$. \end{notation} \begin{definition}[Metodo a un passo] Si dice \textbf{metodo a un passo} un metodo di risoluzione approssimato di un problema di Cauchy che imposti una sequenza della forma: \[ \begin{cases} u_{n+1} = u_n + h \phi(t_n, u_n, f_n, h) & 0 \leq n \leq N_h - 1, \\ u_0 = y_0. \end{cases} \] \end{definition} \begin{definition} Sia $\eps_{n+1} := y_{n+1} - (y_n + h \phi(t_n, y_n, f(t_n, y_n), h))$ la differenza tra il valore esatto $y_{n+1} = y(t_{n+1})$ e quello dato dal metodo applicato al valore esatto $y_{n+1} = y(t_{n+1})$. Si definisce allora \textbf{errore locale di troncamento} al nodo $(n+1)$-esimo il valore: \[ \tau_{n+1}(h) = \frac{\eps_{n+1}}{h}. \] Si definisce inoltre l'\textbf{errore globale di troncamento} come massimo del modulo degli errori locali: \[ \tau(h) = \max_{n=0,\ldots,N_h - 1} \abs{\tau_{n+1}(h)}. \] \end{definition} \begin{definition}[Consistenza del metodo] Si dice che un metodo è consistente di ordine $p$ per il problema se $\tau(h) = o(h)$ ($\lim_{h \to 0} \tau(h) = 0$) e se $\tau(h) = O(h^p)$ per $h \to 0$. \end{definition} \begin{definition}[Convergenza del metodo] Posto $e_n = y_n - u_n$, detto \textit{errore globale}, un metodo si dice \textbf{convergente} di ordine $p$ se esiste $C(h) = o(h)$ ($\lim_{h \to 0} C(h) = 0$), $C(h) = O(h^p)$ per $h \to 0$ tale per cui $\abs{y_n - u_n} \leq C(h)$ per ogni $n$ della successione. \end{definition} \subsection{Metodo di Eulero} Si definisce \textbf{metodo di Eulero} il metodo a un passo che si ottiene imponendo $\phi(t_n, u_n, f_n, h) = f_n$, ovverosia: \[ \begin{cases} u_{n+1} = u_n + h f_n & 0 \leq n \leq N_h - 1, \\ u_0 = y_0. \end{cases} \] L'idea del metodo di Eulero è quella di approssimare innanzitutto i reali con una successione $t_n$ prendendo $h \ll 1$ (in questo modo $N_h \gg 1$ e si ottiene una sequenza arbitrariamente lunga di reali). In questo modo l'equazione \[ y'(t) = f(t, y(t)) \] può essere approssimata ponendo l'equazione vera sui $t_n$ con $y'(t) \approx \frac{u_{n+1} - u_n}{h}$ e $f(t, y(t)) \approx f_n$. La seconda equazione $y(t_0) = y_0$ si riconduce invece facilmente a $u_0 = y_0$. \smallskip Per ottenere dunque un'approssimazione è dunque necessario risolvere la relazione di ricorrenza tra $u_{n+1}$ e $u_n$, ponendo $u_0 = y_0$, o calcolare direttamente la sequenza. Per il metodo di Eulero vale che: \[ \eps_{n+1} = y_{n+1} - (y_n + h f(t_n, y_n)). \] Il metodo di Eulero è sempre consistente di ordine $1$ ed è sempre convergente, ancora di ordine $1$. \subsubsection{Metodo di Eulero implicito} Si definisce \textbf{metodo di Eulero implicito} il metodo che si ottiene sostituendo a $\phi(t_n, u_n, f_n, h)$ il valore $f_{n+1}$, ovverosia: \[ \begin{cases} u_{n+1} = u_n + h f_{n+1} & 0 \leq n \leq N_h - 1, \\ u_0 = y_0. \end{cases} \] Tale metodo è consistente di ordine $1$ e convergente, ancora di ordine $1$. \begin{example}[Problema test di Dahlquist] Si consideri il seguente problema di Cauchy: \[ \begin{cases} y' = \lambda y, \\ y(0) = 1, \end{cases} \] con $I = [0, \infty)$. Applicando il metodo di Eulero implicito si deve risolvere il sistema: \[ \begin{cases} u_{n+1} = u_n + \lambda u_{n+1}, u_0 = 1, \end{cases} \] che ha come soluzione $u_n = \left(\frac{1}{1-h\lambda}\right)^n$. \end{example} \vfill \hrule ~\\ Opera originale di Mario Zito, modifiche e aggiornamenti a cura di Gabriel Antonio Videtta. ~\\Reperibile su \url{https://notes.hearot.it}. \end{multicols} \end{document}