\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}}<B^{1-t}, \quad \abs{\frac{\tx-x}{\tx}}<B^{1-t}, \]
 e si definisce pertanto $u := B^{1-t}$ come la \textbf{precisione di macchina}. In generale, se $y$ è un'approssimazione di $x \in \RR$
 tale per cui $\abs{(y-x)/x}$ è minore di $B^{1-c}$ con $c$ intero, si dice che $y$ ha $c$ \textit{cifre significative}. \medskip

 
  Se invece $\pp(x)<-m$ o $\pp(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}<u. \]

  \subsubsection{Precisione di un insieme \texorpdfstring{$\ff$}{ℱ}}
  
  Dati due sistemi floating point in base $B_1$ e $B_2$ di $t_1$ e $t_2$ cifre, il primo è più preciso del secondo se e solo se vale:
  \[ B_1^{1-t_1} < B_2^{1-t_2} \iff t_1>(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\eta<u. \]
  Gli errori relativi $\delta$ e $\eta$ sono detti \textbf{errori locali} generati dall'operazione $\op$. Un'operazione del tipo $\fl \circ \op$ è detta
  \textbf{flop}.

  
  \subsection{Errori nel calcolo di una funzione}

  Si utilizza il simbolo $\doteq$ per indicare l'uguaglianza di termini
  al \textit{primo ordine}. Si scrive $\varepsilon(x)$ per valutare
  l'errore nella variabile $x$.

  \subsubsection{Errore inerente}
  
  L'errore relativo di $f(\tx)$ rispetto $f(x)$ viene detto \textbf{errore inerente} e vale:
  \[ \ein=\frac{f(\tx)-f(x)}{f(x)}. \]
  Tale errore misura il \textit{condizionamento} di un problema, ossia
  quanto varia il risultato perturbando l'input. \smallskip

    Nel caso di funzioni $f:\RR^n \to \RR$ sufficientemente regolari vale che:
    \[ \ein 
\doteq \sum_{i=1}^n C_i \delta_i, \] 
    dove:
    \[ \delta_i:=\eps_{x_i}=\frac{\tx_i-x_i}{x_i}, \quad C_i=x_i \frac{\partial f/\partial x_i(x)}{f(x)}. \]
    I termini $C_i$ sono detti \textbf{coefficienti di amplificazione} rispetto alla variabile $x_i$. \smallskip

  Un problema si dice \textbf{ben condizionato} se una piccola variazione nelle condizioni iniziali produce una piccola variazione in output e si dice \textbf{mal condizionato}
  altrimenti. Il condizionamento di un problema varia a seconda della grandezza in modulo dei coefficienti di amplificazione (più sono grandi, più il problema è mal condizionato).

  Per le $4$ operazioni aritmetiche elementari $x \bop y$ sono presentati i relativi coeff.~di amplificazione:
    \begin{center}
    \begin{tabular}{|c|c|c|}
    \hline
     Operazione ($\op$) & $C_1$ & $C_2$ \\ \hline
     addizione ($+$) & $x/(x+y)$ & $y/(x+y)$ \\ \hline
     sottrazione ($-$) & $x/(x-y)$ & $-y/(x-y)$ \\ \hline
     moltiplicazione ($\cdot$) & 1 & 1 \\ \hline
     divisione ($/$) &  1 & -1 \\ \hline
    \end{tabular}
    \end{center}
  
  \subsubsection{Funzioni razionali ed errore algoritmico}
  Si ricorda che una funzione $f$ si dice \textit{razionale} se
  si può esprimere tramite una formula \textit{finita} che comprende le quattro
  operazioni aritmetiche. \smallskip
  
  Se la funzione $f : \RR^n \to \RR$ è
  razionale, si può agevolmente calcolarne l'approssimazione
  $\varphi = \tilde{f}$ tramite troncamento. Generalmente esistono più modi di implementare la priorità delle operazioni, e una tale scelta di priorità $\varphi$
  è detta \textit{algoritmo}. \smallskip

  Il discostamento relativo di $\varphi(\tx)$ da $f(\tx)$ (il valore effettivamente
  calcolato sull'input perturbato)
  viene detto \textbf{errore algoritmico}, che pertanto è così definito:
  \[ \ealg=\frac{\varphi(\tx)-f(\tx)}{f(\tx)}. \]

  Un problema si dice \textbf{numericamente stabile} (in avanti) se $\abs{\ealg} < ku$ per $k$ costante, e \textbf{numericamente instabile} (in avanti) se non lo è.

  \subsubsection{Funzioni regolari non razionali ed errore analitico}

  Se $f : \RR^n \to \RR$ non è
  razionale, ma è comunque sufficientemente regolare, si può approssimare $f$ tramite una funzione
  razionale $g$. Chiamiamo \textbf{errore analitico} il discostamento relativo
  di $g$ da $f$ sull'input non perturbato (ossia a priori dell'uso di
  aritmetica di macchina):
  \[ \ean=\frac{g(x)-f(x)}{f(x)}. \]
  Questo errore può essere studiato attraverso gli strumenti
  dell'analisi matematica e della teoria dell'approssimazione
  delle funzioni (e.g.~polinomi di Taylor e resti di Peano, Lagrange, Cauchy o resto integrale).

  \subsubsection{Errore totale}
  
  Definiamo l'\textbf{errore totale} come il discostamento relativo dell'algoritmo $\varphi$ valutato sull'input perturbato dal valore esatto $f(x)$: \[ \etot=\frac{\varphi(\tx)-f(x)}{f(x)}. \]
  
  \begin{proposition}
      A meno di termini non lineari vale sempre che:
      \[
        \etot \doteq \ein + \ealg + \ean(\tx).
      \]
  \end{proposition}

  Pertanto, se la funzione $f$ è razionale,
  $\etot \doteq \ein + \ealg$. Questa proposizione giustifica
  lo studio separato dei tre errori.
  
  \subsubsection{Calcolo dell'errore algoritmico e dell'errore totale (per funzioni razionali)}

  L'unica differenza tra lo studio dell'errore algoritmico e dell'errore totale (per funzioni razionali)
  secondo il metodo presentato sta nell'errore $\eps_{x_i}$ che viene dato alle variabili
  di input. Infatti, nel caso dell'errore algoritmico, le variabili sono
  già ben rappresentate come numeri di macchina e quindi non hanno errore,
  mentre nel caso dell'errore totale lo hanno. Pertanto, se una traccia
  richiede di calcolare un errore totale su dei numeri di macchina, ciò è
  equivalente a calcolarne l'errore algoritmico. \smallskip

  Per procedere allo studio dell'errore algoritmico conviene utilizzare dei grafi orientati con le op.~aritmetiche segnando ogni flop come un nuovo nodo,
  a cui sono rivolte le frecce da tutte le variabili impiegate per portare
  a termine l'op.~aritmetica. Il più
  semplice di questi grafi, che rappresenta $S = A \bop B$, è presentato
  di seguito:

    % https://q.uiver.app/#q=WzAsNixbMSwwLCJBIl0sWzMsMCwiQiJdLFsyLDEsIlMiXSxbNCwwXSxbMCwwXSxbMiwyXSxbMCwyLCJDX0EiLDIseyJsYWJlbF9wb3NpdGlvbiI6NjB9XSxbMSwyLCJDX0IiLDAseyJsYWJlbF9wb3NpdGlvbiI6NjB9XSxbNSwyLCJcXGRlbHRhIiwxLHsibGFiZWxfcG9zaXRpb24iOjEwMCwib2Zmc2V0IjoyLCJzdHlsZSI6eyJib2R5Ijp7Im5hbWUiOiJub25lIn0sImhlYWQiOnsibmFtZSI6Im5vbmUifX19XSxbNCwwLCJcXHZhcmVwc2lsb25fQSIsMCx7ImxhYmVsX3Bvc2l0aW9uIjo5MCwib2Zmc2V0Ijo1LCJzdHlsZSI6eyJib2R5Ijp7Im5hbWUiOiJub25lIn0sImhlYWQiOnsibmFtZSI6Im5vbmUifX19XSxbMSwzLCJcXHZhcmVwc2lsb25fQiIsMCx7ImxhYmVsX3Bvc2l0aW9uIjoyMCwib2Zmc2V0Ijo1LCJzdHlsZSI6eyJib2R5Ijp7Im5hbWUiOiJub25lIn0sImhlYWQiOnsibmFtZSI6Im5vbmUifX19XV0=
    \[\begin{tikzcd}
    	{} & A && B & {} \\
    	&& S \\
    	&& {}
    	\arrow["{\varepsilon_A}"{pos=0.9}, shift right=5, draw=none, from=1-1, to=1-2]
    	\arrow["{C_A}"'{pos=0.6}, from=1-2, to=2-3]
    	\arrow["{\varepsilon_B}"{pos=0.2}, shift right=5, draw=none, from=1-4, to=1-5]
    	\arrow["{C_B}"{pos=0.6}, from=1-4, to=2-3]
    	\arrow["\delta"{description, pos=1}, shift right=2, draw=none, from=3-3, to=2-3]
    \end{tikzcd}\]
    \vskip -0.1in
    In tal caso vale che:
    \[ \eps_S=\delta+C_A\eps_A+C_B\eps_B. \]
    Una volta ricavato l'errore, è sufficiente maggiorarlo in modulo con una
    certa funzione di $u$ per determinare se l'algoritmo è stabile o meno.

    \subsubsection{Cancellazione numerica}

    In generale per la scelta di un algoritmo è consigliabile \underline{evitare} di \textit{sottrarre numeri dello stesso segno} o di \textit{sommare numeri discordi}. Infatti, se la somma o la differenza di tali numeri è molto vicina
    a $0$, i coefficienti di amplificazione hanno modulo molto maggiore
    di $1$, e dunque il problema diviene \textit{numericamente instabile}.

    \subsubsection{Analisi dell'errore all'indietro}

    Sia $f(x_1,\dots,x_n)$ una funzione razionale e si denoti con $\varphi(x_1,\dots,x_n)$ l'algoritmo definito su $\mathcal{F}^n$ i cui valori
    sono ottenuti calcolando $f(x_1,\dots,x_n)$ con l'aritmetica di macchina.
    
    Nell'analisi \textit{all'indietro} dell'errore si cercano delle perturbazioni $\delta_i$ tali che denotando $\hat{x}_i=x_i(1+\delta_i)$ risulti:
    \[ \varphi(x_1,\dots,x_n)=f(\hat{x}_1,\dots,\hat{x}_n), \]
    ovverosia valutare l'algoritmo $\varphi$ su $(x_1, \ldots, x_n)$ diviene equivalente
    a calcolare il valore esatto di $f(\hat{x}_1, \ldots, \hat{x}_n)$. \smallskip
    
    Agendo in questo modo per calcolare l'errore algoritmico $\ealg$ è
    sufficiente calcolare in realtà l'errore inerente $\ein$ sulle perturbazioni $\delta_i$. \smallskip

    Si riporta lo schema risolutivo di un tipico esercizio sull'analisi all'indietro:
    \begin{itemize}
        \item Si calcoli $f(\hat{x}_1, \ldots, \hat{x}_n)$, dove $\hat{x}_i = x(1+\delta_i)$ e i vari $\delta_i$ sono variabili libere.
        \item Si calcoli $\varphi(x_1, \ldots, x_n)$ sostituendo a ogni $j$-esima flop
        $\fl(\tilde{s_1} \bop \tilde{s_2})$ il valore esatto $(\tilde{s_1} \bop \tilde{s_2})(1+\eps_j)$ dove $\abs{\eps_j} < u$, procedendo poi a ritroso
        nelle flop di $\tilde{s_1}$ e di $\tilde{s_2}$.
        \item Si risolva il sistema $f(\hat{x}_1, \ldots, \hat{x}_n) = \varphi(x_1, \ldots, x_n)$ nei $\delta_i$ in funzione degli $\eps_j$, maggiorandoli alla fine.
    \end{itemize}
    
    Un algoritmo si dice \textbf{numericamente stabile all'indietro} se è
    possibile effettuare un analisi all'indietro per cui $\ealg < ku$ con
    $k$ costante.

    \section{Algebra lineare numerica}

    \subsection{Matrici di permutazione e matrici riducibili}

    \begin{definition}
        Si dice che $P \in \CC^{n \times n}$ è una \textbf{matrice di permutazione}
        se esiste $\sigma \in S_n$ tale per cui $P^j = e_{\sigma(j)}$; e in
        tal caso si scrive $P_\sigma := P$ per rimarcare la permutazione a cui
        $P$ è associata.
    \end{definition}

    Si osserva in particolare che $\sigma \mapsto P_\sigma$ è un omomorfismo
    da $S_n$ alle matrici di permutazioni, che dunque formano un gruppo moltiplicativo. Per le matrici di permutazioni valgono dunque le seguenti
    proprietà:

    \begin{itemize}
        \item $P_\sigma \in \RR^{n \times n}$,
        \item $P_\sigma P_\tau = P_{\sigma \circ \tau}$,
        \item $P_\sigma^T P_\sigma = P_\sigma P_\sigma^T = I$ ($P_\sigma$ è ortogonale e unitaria),
        \item $P_\sigma^{-1} = P_\sigma^T = P_{\sigma^{-1}}$,
        \item $\det(P_\sigma) = \sgn(\sigma)$,
        \item $P_\sigma$ rappresenta la matrice di cambio di base da quella
        canonica a quella canonica permutata secondo $\sigma$, ovverosia
        da $(e_1, \ldots, e_n)$ a $\left(e_{\sigma(1)}, \ldots, e_{\sigma(n)}\right)$,
        \item Se $A \in \CC^{n \times n}$ e $B = P_\sigma^T A P_\sigma$, allora
        $b_{ij} = a_{\sigma(i)\sigma(j)}$.
        \item Analogamente, se $A \in \CC^{n \times n}$ e $B = P_\sigma A P_\sigma^T = P_{\sigma^{-1}}^T A P_{\sigma^{-1}}$, allora
        $b_{ij} = a_{\sigma^{-1}(i)\sigma^{-1}(j)}$, e dunque
        $b_{\sigma(i)\sigma(j)} = a_{ij}$.
    \end{itemize}

    \begin{definition}[Matrice riducibile]
        Una matrice $A\in \CC^{n \times n}$ si dice \textbf{riducibile} se esiste una matrice di permutazione $P$ per cui:
        \[ PAP^\top=PAP\inv=\begin{pmatrix}
            A_{11} & A_{12} \\
            0 & A_{22} \\
        \end{pmatrix} \]
    con $A_{11}$ e $A_{22}$ matrici quadrate, ovverosia se $PAP^\top$ è triangolare a blocchi. \smallskip

    Una matrice si dice \textbf{irriducibile} se non è riducibile.
    \end{definition}

    \vskip -0.1in
    
    Equivalentemente una matrice è riducibile se esiste un sottoinsieme proprio
    non vuoto $I$ di $\{e_1, \ldots, e_n\}$ tale per cui $A\cdot I \subseteq \Span(I)$, ovverosia se $I$ è $A$-invariante.

    \subsubsection{Grafo associato a una matrice}
    
    Data $A\in \CC^{n \times n}$ si definisce \textbf{grafo associato} $G[A]$ ad $A$ come
    il grafo \underline{orientato} con le seguenti proprietà:
    
    \begin{itemize}
        \item $G[A]$ ha come nodi $\{1, \ldots, n\}$,
        \item Esiste un arco orientato da $i$ a $j$ se e solo se $a_{ij} \neq 0$. 
    \end{itemize}

    Un cammino è una sequenza di archi $i_1 \to i_2 \to \cdots i_m$.
    Un grafo si dice \textbf{fortemente connesso} se $\forall i$, $j$ esiste
    un cammino da $i$ a $j$.

    Se $P$ è una matrice di permutazione, $B := PAP^\top$ ha lo stesso
    grafo di $A$ a meno di rinominare i nodi secondo $\sigma$, ovverosia
    esiste un arco da $\sigma(i)$ a $\sigma(j)$ in $G[B]$ se e solo se
    esiste un arco da $i$ a $j$ in $G[A]$. Pertanto
    il grafo di $A$ è fortemente connesso se e solo se
    quello di $PAP^\top$ lo è. Questo fatto risulta fondamentale nel
    dimostrare il seguente teorema:

    \begin{theorem}
    $A$ è irriducibile se e solo se il suo grafo associato
    è fortemente connesso. Equivalentemente, $A$ è riducibile se e solo se
    il suo grafo associato \underline{non} è fortemente connesso.
    \end{theorem}
    
    \subsection{Teoremi di Gershgorin e applicazioni}

    Sia $A\in \CC^{n \times n}$. Allora si definisce l'$i$-esimo \textbf{cerchio di
    Gershgorin} come l'insieme $K_i \subset \CC$ tale per cui:
    \[
        K_i=\left\{z\in\CC:\abs{z-a_{ii}} \leq r_i\right\},
    \]
    dove
    \[
        r_i = \sum_{j=1, \, j\neq i}^n\abs{a_{ij}}
    \]
    è detto $i$-esimo \textbf{raggio di Gershgorin}. In particolare
    $K_i$ è un cerchio in $\CC \cong \RR^2$ di centro $a_{ii}$ e raggio
    $r_i$.

    \subsubsection{Primo teorema di Gershgorin}

    \begin{theorem}[Gershgorin I]
    Data $A\in \CC^{n \times n}$, gli autovalori di $A$ appartengono all'insieme
    $\bigcup_{i=1}^n K_i$, ossia per ogni autovalore $\lambda$ di $A$ esiste
    un $i$ tale per cui $\lambda$ appartiene all'$i$-esimo cerchio $K_i$ di
    Gershgorin. \smallskip

    Inoltre, se $v$ è un autovettore relativo all'autovalore $\lambda$,
    $\lambda \in K_h$, dove $h = \argmax  \abs{v_i}$, ovverosia
    $h$ è l'indice tra gli $i$ per cui $\abs{v_i}$ è massimo. 
    \end{theorem}


    Dal momento che $A^\top$ è simile ad $A$, $A^\top$ e $A$ condividono gli stessi autovalori. Si può dunque rafforzare la tesi del
    teorema cercando gli autovalori nell'intersezione delle unioni dei cerchi relativi ad $A$ e a $A^\top$, ovverosia, se $H_i$ rappresenta l'$i$-esimo cerchio di Gershgorin di $A^\top$, ogni autovalore di $A$ appartiene a:
    \[ \left(\bigcup_{i=1}^n K_i\right) \cap \left(\bigcup_{i=1}^n H_i\right). \]
    Per lo stesso motivo è utile coniugare $A$ per similitudine, specie se si coniuga $A$ utilizzando matrici diagonali, dacché il coniugio per matrice diagonale, oltre a non alterare gli autovalori, lascia invariati i centri dei cerchi e riscalare i raggi. \smallskip

    In particolare, se $D = \diag(d_1, \ldots, d_n)$, allora
    $DAD\inv = (d_i a_{ij} d_j^{-1})$ e $D\inv A D = (d_i^{-1} a_{ij} d_j)$.
    Pertanto vale il seguente teorema:

    \begin{theorem}
        Sia $A \in \CC^{n \times n}$. Allora gli autovalori di $A$ appartengono
        all'insieme:
        \[
            \bigcap_{d \in D} \bigcup_{i=1}^n \left\{ z \in \CC : \abs{z - a_{ii}} \leq \frac{1}{\abs{d_i}} \sum_{j=1, \, j\neq i}^n \abs{a_{ij}} \abs{d_j} \right\},
        \]
        dove $D$ è un sottoinsieme non vuoto di $\CC^n$.
    \end{theorem}

    \subsubsection{Matrici fortemente dominanti diagonali}

    Una matrice $A$ si dice \textbf{fortemente dominante diagonale} se $\abs{a_{ii}}>r_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<p<1$ l'espressione (generalizzata) $\norm{x}_p$ non è una norma. \smallskip


    Data una norma $\norm{\cdot}$ su $\CC^n$ e data $S\in \CC^{n \times n}$ non singolare, allora esiste la norma (indotta da) $S$, definita in
    modo tale che:
    \[ \norm{\cdot}_S:x\mapsto\norm{Sx}. \]

    Inoltre vale la seguente caratterizzazione:
    \begin{proposition}
         Una norma $\norm{\cdot}$ è indotta da un prodotto scalare se e solo se vale la \textit{legge del parallelogramma},
        ossia se e solo se:
        \[ \norm{x+y}^2+\norm{x-y}^2=2(\norm{x}^2+\norm{y}^2) \]
    \end{proposition}

    A partire da questo risultato, si deduce che le norme $1$ e $\infty$ non sono indotte da nessun prodotto scalare. \medskip

    Ogni norma è una \textbf{funzione uniformemente continua}, cioè vale che, $\forall \varepsilon > 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 $i<j$ ($C = -\triu(A)$).
    \end{itemize}
    Il metodo iterativo ottenuto con il partizionamento additivo $M=D$ e $N=B+C$ è detto \textbf{metodo di Jacobi}, mentre
    ponendo $M=D-B$ e $N=C$ si applica il \textbf{metodo di Gauss-Seidel}. \smallskip
    
    Le matrici di iterazione $P=M\inv N$ dei due metodi illustrati sono rispettivamente:
    \[ J=D\inv(B+C), \quad G=(D-B)\inv C. \]

    \begin{theorem}
        Se vale una delle seguenti condizioni allora $\rho(J)$ e $\rho(G)<1$ (e quindi $A$ è invertibile e il metodo iterativo è convergente):
        \begin{enumerate}
            \item $A$ è fortemente dominante diagonale;
            \item $A^T$ è fortemente dominante diagonale;
            \item $A$ è irriducibilmente dominante diagonale (i.d.d.);
            \item $A^T$ è irriducibilmente dominante diagonale (i.d.d.).
        \end{enumerate}
    \end{theorem}
    
    L'iterazione del metodo di Jacobi si scrive come:
    \[x^{(k+1)}=D\inv((B+C)x^{(k)}+b),\]
    che diventa:
    \[x_i^{(k+1)}=\frac{1}{a_{ii}}\left(b_i-\sum_{j=1, \,j\neq i}^na_{ij}x_j^{(k)}\right),\]
    richiedendo $O(n^2)$ flops. \smallskip

    L'iterazione del metodo di Gauss-Seidel si scrive invece come
    \[ x^{(k+1)}=(D-B)\inv(Cx^{(k)}+b),\] \[ x^{(k+1)}=D\inv(Bx^{(k+1)}+Cx^{(k)}+b), \]
    che diventa:
    \[ x_i^{(k+1)}=\frac{1}{a_{ii}}\left(b_i-\sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)}-\sum_{j=i+1}^na_{ij}x_j^{(k)}\right), \]
    richiedendo sempre $O(n^2)$ flops. \smallskip

    Il metodo di Jacobi esegue le operazioni su $x^{(k)}$ mentre quello di Gauss-Siedel usa anche le componenti già aggiornate di $x^{(k+1)}$. In generale dunque, il metodo di Gauss-Siedel ha migliori proprietà di convergenza, anche se non è sempre così.


    \begin{theorem}[di Stein-Rosenberg]
        Se $a_{ii}\neq0$ $\forall i$ e $J$ ha elementi non negativi allora vale una sola delle seguenti proprietà:
        
        \begin{itemize}
                \item $\rho(J)=\rho(G)=0$;
                \item $0<\rho(G)<\rho(J)<1$;
                \item $\rho(J)=\rho(G)=1$;
                \item $1<\rho(J)<\rho(G)$.
        \end{itemize}
    \end{theorem}

    \begin{theorem}
        Se $A$ è una matrice tridiagonale con elementi diagonali non nulli, allora per ogni autovalore $\lambda$ di $J$ esiste un autovalore $\mu$ di $G$ tale che $\mu=\lambda^2$. 
        Per ogni autovalore non nullo di $G$ esiste un autovalore $\lambda$ di $J$ tale che $\mu=\lambda^2$. In particolare $\rho(G)=\rho(J)^2$.
    \end{theorem}

    \subsubsection{Metodi a blocchi}

    Se $A\in M_{nm}(\CC)$ è partizionata in $n^2$ blocchi $m\times n$ $A=(A_{ij})$ possiamo considerare una decomposizione additiva $A=D-B-C$ dove $D$ matrice diagonale a blocchi con blocchi diagonali uguali a $A_{ii}$, $B$ la matrice triangolare inferiore a blocchi tale che $B_{ij}=-A_{ij}\;\con i>j\E C$ la matrice triangolare superiore a blocchi tale che $C_{ij}=-A_{ij}\;\con i<j$. \medskip

    
    In questo modo se i blocchi diagonali $A_{ii}$ sono non singolari possiamo considerare i metodi iterativi con $M=D$ e $N=B+C$, e con $M=D-B$ e $N=C$. Il primo è detto \textbf{metodo di Jacobi a blocchi}, mentre il secondo è detto \textbf{metodo di Gauss-Siedel a blocchi}.


    \begin{theorem}
        Sia $A=(A_{ij})$ una matrice tridiagonale a blocchi, cioè tale che $A_{ij}=0\;\se|i-j|\ge2$ con blocchi diagonali $A_{ii}$ non singolari. Siano $J_B$ e $G_B$ le matrici di iterazione dei metodi di Jacobi a blocchi e di Gauss-Siedel a blocchi. \nl Allora per ogni autovalore $\lambda$ di $J_B$ esiste un autovalore $\mu$ di $G_B$ tale che $\mu=\lambda^2$. Per ogni autovalore non nullo $\mu$ di $G_B$ esiste un autovalore $\lambda$ di $J_B$ tale che $\mu=\lambda^2$. In particolare vale che $\rho(G_B)=\rho(J_B)^2$.
    \end{theorem}

    \section{Calcolo di zeri di funzioni continue su \texorpdfstring{$\RR$}{ℝ}}
    
    Sia $f: [a,b] \to \RR$ continua tale per cui $f(a)f(b)<0$ (i.e.~$f(a)$ ed $f(b)$ sono discordi). Allora, per il teorema degli zeri (o di Bolzano), esiste un punto $\alpha\in[a,b]$ tale
    per cui $f(\alpha)=0$, ossia esiste uno zero in $(a, b)$. Lo scopo di questa
    sezione è illustrare i metodi
    che permettono di generare una successione $\{x_k\}$ che converga ad un
    tale $\alpha$.

    \subsection{Metodo di bisezione (o \textit{dicotomico})}

    Uno dei metodi più semplici, e il più costruttivo per la dimostrazione
    del teorema degli zeri, è il metodo di \textbf{bisezione} (o \textit{dicotomico}). Innanzitutto
    si presuppone $f(b) > 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<g'(x)<1$ con $x\in[\alpha-\rho,\alpha+\rho]$ allora $x_{k+1}-\alpha$ ha lo stesso segno di $x_k-\alpha$. Vale
    quindi che:
    \[ x_0>\alpha\Rightarrow x_k>\alpha \; \forall k, \quad \alpha<x_{k+1}<x_k, \]
    cioè la successione è decrescente. Analogamente se $x_0<\alpha$ la successione è crescente. \medskip

    Se invece $-1<g'(x)<0$, la differenza di un punto con $\alpha$ alterna segno,
    in particolare per $x_0 > \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<i}(x_i-x_j)$.
    \end{theorem}

    Se $a=(a_i)$ è il vettore dei coefficienti di $p(x)$ e $y=(y_i)$ è il vettore
    delle immagini nei nodi, il problema dell'interpolazione polinomiale si
    riduce a risolvere il sistema lineare $V_n a = y$. \smallskip

    Il problema ha dunque un costo computazionale di $\frac{2}{3} (n+2)^3$ (ossia
    quello di risoluzione di un sistema lineare), ma risulta essere
    numericamente instabile.

    \subsubsection{Interpolazione di Lagrange}

    \begin{definition}[polinomio di Lagrange]
        Si definisce l'$i$-esimo polinomio di Lagrange $L_i(x)$ come:
        \[
            L_i(x) = \frac{\prod_{j=0,j\neq i}^n(x-x_j)}{\prod_{j=0,j\neq i}^n(x_i-x_j)}.
        \]
        Questo polinomio è tale per cui $L_i(x_i) = 1$ e $L_i(x_j) = 0$ per
        $j \neq i$.
    \end{definition}

    Considerando dunque come base dell'interpolazione $\phi_i = L_i$,
    si ottiene che:
    \[
        p(x) = \sum_{i=0}^n y_i L_i(x).
    \]
    Infatti, dacché $V_n$ è invertibile ($\det V_n \neq 0$), allora
    esiste un solo polinomio di grado al più $n+2$ che interpola i nodi
    $(x_i, y_i)$, e $p(x)$ soddisfa le condizioni di interpolazione.

    Inoltre, per $x \neq x_i$, vale la seguente riscrittura:
    \[
        p(x) = \left[ \prod_{i=0}^n (x-x_i) \right] \left[ \sum_{i=0}^n \frac{y_i/\theta_i}{x-x_i} \right], \quad \theta_i = \prod_{\substack{j=0 \\ j \neq i}}^n (x_i - x_j).
    \]

    Tramite questa riscrittura il calcolo del primo valore di $p(x)$ impiega $O(n^2)$ flops, ed il costo della valutazione in un nuovo valore costa solo $O(n)$ flops (dacché $\sigma_i$ è già calcolato).

    \subsection{Resto dell'interpolazione polinomiale}

    Se $f(x)$ è sufficientemente regolare, allora definiamo il \textbf{resto dell'interpolazione} come:
    \[ r_n(x)=f(x)-p_n(x), \]
    dove $p_n(x)$ è il polinomio di interpolazione di $f(x)$ relativo ai nodi $x_0<\dots<x_n$.

    \begin{theorem}
        Sia $f(x)\in C^{n+1}([a,b])$ e sia $p_n(x)$ il polinomio di interpolazione di $f(x)$ relativo ai nodi $a\le x_0<\dots<x_n\le b$.
        Allora per ogni $x\in[a,b]$ esiste $\xi\in(a,b)$ tale per cui: \[r_n(x)=\left[\prod_{i=0}^n(x-x_i)\right]\frac{f^{(n+1)}(\xi)}{(n+1)!}.\]
    \end{theorem}
    
    \theorem{Sia $f(x)\in C^{n+1}[a,b]\E p(x)$ il polinomio di interpolazione di $f(x)$ relativo ai nodi $a\le x_0<\dots<x_n\le b$. Per ogni $x\in[a,b]$ esiste $\xi\in(a,b)$ tale che $$r_n(x)=\prod_{i=0}^n(x-x_i)\frac{f^{(n+1)}(\xi)}{(n+1)!}$$}

    \begin{remark}
        Avvicinandosi con tutti i nodi a $x_i$ si ottiene un'espressione simile
        a quella del resto di Lagrange (ponendo infatti $x_i = x_j$ per ogni $i$, $j$, la formula è proprio quella del resto di Lagrange).
    \end{remark}

    \subsection{Interpolazione polinomiale nelle radici \texorpdfstring{$n$}{n}-esime dell'unità}

    \begin{definition}
        Si dice \textbf{radice $n$-esima dell'unità} una radice
        di $x^n - 1$ su $\CC$. Si dice \textbf{radice primitiva $n$-esima
        dell'unità} il numero complesso $\omega_n$ tale per cui:
        \[
            \omega_n=\cos{\frac{2\pi}{n}}+i\sin{\frac{2\pi}{n}}.
        \]
    \end{definition}

    Si verifica facilmente che le radici $n$-esime dell'unità formano
    un gruppo moltiplicativo generato da $\omega_n$. Le radici dell'unità
    soddisfano la seguente proposizione:

    \begin{proposition}
        Vale l'identità:
        \[ \sum_{i=0}^{n-1}\omega_n^{ki}=\begin{cases}
            n & \se k\equiv0\!\!\!\pmod{n}, \\
            0 & \altrimenti\end{cases} \]
    \end{proposition}

    Scegliamo come nodi dell'interpolazione le radici $n$-esime dell'unità, detti
    \textit{nodi di Fourier}.
    Dato il vettore $y=(y_0,\dots,y_{n-1})$, il nostro obiettivo è dunque 
    trovare i coefficienti $z=(z_0,\dots,z_{n-1})$ del polinomio 
    $p(t)=\sum_{j=0}^{n-1}z_jt^j$ che soddisfa $p(\omega_n^i)=y_i$ $\forall i$.

    \begin{definition}
        Si definisce $z = \DFT_n(y) = \DFT(y)$ come la
        \textbf{trasformata discreta di Fourier}, e
        $y = \IDFT_n(z) = \IDFT(z)$ come la \textbf{trasformata discreta
        inversa di Fourier} (l'inverso è ben definito
        dal momento che le radici $n$-esime sono distinte).
    \end{definition}

    \begin{definition}
        Si definisce \textbf{matrice di Fourier} di taglia $n-1$ la
        matrice di Vandermonde nei nodi di Fourier:
        \[
            \Omega_n = (\omega_n^{ij}) = (\omega_n^{ij \bmod{n}}).
        \]
        $\Omega_n$ soddisfa la relazione $\DFT(y) = \Omega_n^{-1} y$ e
        $\IDFT(z) = \Omega_n z$.
    \end{definition}

    \begin{proposition}
        Valgono le seguenti proprietà per la matrice di Fourier:
        \begin{itemize}
            \item $\Omega_n=\Omega_n^T$;
            \item $\Omega_n^H\Omega_n=nI$ (pertanto $F_n := \Omega_n/\sqrt{n}$ è unitaria);
            \item $\Omega_n^2=n \, \Pi_n$, dove $\Pi_n$ è la matrice di permutazione corrispondente alla permutazione $\sigma$ su $\{0, \ldots, n-1\}$ tale per cui $\sigma(0)=0$ e $\sigma(j)=n-j$ per $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}
