2  Calcolo Numerico

Autore/Autrice
Affiliazione

Enrico Bertolazzi

University of Trento, Department of Industrial Engineering

Formulario del corso. Le sezioni seguono l’ordine dei capitoli: ognuna rimanda al capitolo in cui l’argomento è svolto per esteso, con le ipotesi e le dimostrazioni. Qui ci sono solo le formule da avere sotto gli occhi.

2.1 Sistemi lineari: metodi diretti

Capitolo di riferimento: Metodo di eliminazione di Gauss, Fattorizzazione di Cholesky

Matrice di Frobenius

\begin{pmatrix} a_1 \\ a_2 \\ \vdots \\ a_n \end{pmatrix} \quad \implies \quad \bm{F}_1 = \bm{I} - \begin{pmatrix} 0 \\ \bm{v}_1 \end{pmatrix} \bm{e}_1^\top \qquad \bm{v}_1= \begin{pmatrix} a_2/a_1 \\ \vdots \\ a_n/a_1 \end{pmatrix}

per cui

\bm{F}_1 \begin{pmatrix} a_1 \\ a_2 \\ \vdots \\ a_n \end{pmatrix} = \begin{pmatrix} a_1 \\ 0 \\ \vdots \\ 0 \end{pmatrix}

ed

\begin{pmatrix} a_1 \\ \vdots\\ a_k \\ a_{k+1} \\ \vdots \\ a_n \end{pmatrix} \quad \implies \quad \bm{F}_k = \bm{I} - \begin{pmatrix} 0 \\ \vdots\\ 0 \\ \bm{v}_k \end{pmatrix} \bm{e}_k^\top \qquad \bm{v}_k= \begin{pmatrix} a_{k+1}/a_k \\ \vdots \\ a_n/a_k \end{pmatrix}

per cui

\bm{F}_k \begin{pmatrix} a_1 \\ \vdots\\ a_k \\ a_{k+1} \\ \vdots \\ a_n \end{pmatrix} = \begin{pmatrix} a_1 \\ \vdots\\ a_k \\ 0 \\ \vdots \\ 0 \end{pmatrix}

inoltre

\bm{L} = \begin{pmatrix} 1 & 0 & \cdots & & 0 \\ & 1 & & & \vdots \\ & & \ddots & & \vdots \\ & & & 1 & \vdots \\ \bm{v}_1 & \bm{v}_2 & & \bm{v}_{n-1} & 1 \\ \end{pmatrix}

Cholesky: passo generico

Per \bm{A} simmetrica definita positiva, \bm{A}=\bm{R}^\top\bm{R} con \bm{R} triangolare superiore. Al primo passo

\bm{A} = \begin{pmatrix} \alpha & \bm{a}^\top \\ \bm{a} & \tilde{\bm{A}} \end{pmatrix} = \bm{R}_1^\top \begin{pmatrix} 1 & \bm{0}^\top \\ \bm{0} & \bm{A}' \end{pmatrix} \bm{R}_1, \qquad \bm{R}_1 = \begin{pmatrix} \sqrt{\alpha} & \dfrac{1}{\sqrt{\alpha}}\bm{a}^\top \\[4pt] \bm{0} & \bm{I}_{n-1} \end{pmatrix}

e si itera su \bm{A}'. Costo \sim n^3/3 flop, metà di Gauss1; nessun pivoting necessario.

2.2 Fattorizzazione QR

Capitolo di riferimento: Metodo QR per sistemi lineari, Costruzione della fattorizzazione QR, La pseudo inversa di Moore-Penrose

Rotazione di Givens

\begin{pmatrix} a \\ b \end{pmatrix} \quad \implies \quad \begin{pmatrix} \sqrt{a^2+b^2} \\ 0 \end{pmatrix}

\bm{G} = \begin{pmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{pmatrix} \qquad \cos\theta = \dfrac{a}{\sqrt{a^2+b^2}},\quad \sin\theta = \dfrac{b}{\sqrt{a^2+b^2}}

Riflessione di Householder

\bm{x}= \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix} \quad \implies \quad \begin{pmatrix} \|\bm{x}\|_2 \\ 0 \\ \vdots \\ 0 \end{pmatrix}

\bm{H} = \bm{I} -2\dfrac{\bm{v}\bm{v}^\top}{\bm{v}\cdot\bm{v}}, \qquad \bm{v} = \bm{x}+s\|\bm{x}\|_2 \bm{e}_1 = \begin{pmatrix} x_1+s\|\bm{x}\|_2\\ x_2 \\ \vdots \\ x_n \end{pmatrix} \qquad s = \mathrm{sign}(x_1)

Pseudo-inversa di Moore-Penrose

\bm{A}^+ è l’unica matrice che soddisfa le quattro condizioni di Penrose2

\bm{A}\bm{A}^+\bm{A}=\bm{A}, \qquad \bm{A}^+\bm{A}\bm{A}^+=\bm{A}^+, \qquad (\bm{A}^+\bm{A})^\top=\bm{A}^+\bm{A}, \qquad (\bm{A}\bm{A}^+)^\top=\bm{A}\bm{A}^+

Se \bm{A} è quadrata non singolare allora \bm{A}^+=\bm{A}^{-1}. Per un vettore, \bm{a}^+=(\bm{a}\cdot\bm{a})^{-1}\bm{a}^\top.

2.3 Interpolazione polinomiale

Capitolo di riferimento: Interpolazione Polinomiale, Interpolazione di Lagrange, Interpolazione di Newton, Esempio di Runge

Errore di interpolazione

f(x) - p(x) = \dfrac{f^{(n+1)}(\zeta)}{(n+1)!}(x-x_0)(x-x_1)\cdots(x-x_n)

con \zeta\in I(x,x_0,x_1,\ldots,x_n)

Interpolazione di Lagrange

L_{i}(x) = \dfrac{\ell_{i}(x)}{\ell_{i}(x_i)}, \qquad \ell_{i}(x) = \prod_{{j=0}\atop{j\neq i}}^n (x-x_j)

p(x) = y_0 L_0(x) + y_1 L_1(x)+\cdots+y_n L_n(x)

Interpolazione di Newton

\begin{aligned} p(x) & = f[x_0]\omega_0(x) \\ & + f[x_0,x_1] \omega_1(x) \\ & + f[x_0,x_1,x_2] \omega_2(x) \\ & \vdots \\ & + f[x_0,x_1,x_2,\ldots,x_n] \omega_n(x) \end{aligned}

base polinomiale per il metodo di Newton

\begin{aligned} \omega_0(x) & = 1 \\ \omega_1(x) & = x-x_0 \\ \omega_2(x) & = (x-x_0)(x-x_1) = \omega_1(x)(x-x_1) \\ & \vdots \\ \omega_n(x) & = (x-x_0)(x-x_1)\cdots (x-x_{n-1}) = \omega_{n-1}(x)(x-x_{n-1}) \end{aligned}

Differenze divise

f[a,\bullet,b] = \dfrac{f[a,\bullet]-f[\bullet,b]}{a-b}

Scelta dei nodi: Chebyshev

Su nodi equispaziati la costante di Lebesgue3 cresce come \Lambda_n \approx \dfrac{2^{n+1}}{n e(\ln n+\gamma)}: da qui il fenomeno di Runge4. I nodi di Chebyshev5 la riportano a crescita logaritmica:

x_k=\dfrac{a+b}{2}+\dfrac{b-a}{2}\cos\left(\pi\dfrac{2k+1}{2n+2}\right), \qquad x_k=\dfrac{a+b}{2}+\dfrac{b-a}{2}\cos\left(\dfrac{k\pi}{n}\right)

(Chebyshev e Chebyshev-Lobatto). L’errore di interpolazione resta entro un fattore 1+\Lambda_n da quello della miglior approssimazione.

2.4 Integrazione numerica

Capitolo di riferimento: Integrazione numerica, Integrazione composita

Trapezi

\int_a^b f(x) dx = \underbrace{\dfrac{h}{2}\left[f(a) + f(b) + 2\sum_{k=1}^{n-1} f(x_i) \right]}_{\text{formula quadratura}} \quad \underbrace{-\dfrac{b-a}{12}h^2 f''(\zeta)}_{\text{termine di errore}}

Simpson

\int_a^b f(x) dx = \dfrac{h}{3}\left[ f(x_0) + 4f(x_1) + 2f(x_2)+4f(x_3) + \cdots + 4 f(x_{n-1}) + f(x_n) \right]-\dfrac{b-a}{180}h^4 f^{(4)}(\zeta)

Errore integrazione per intervallo [-1,1]

\begin{aligned} \mathrm{Trapezi: } \quad c & = -\dfrac{2}{3} g''(z)\\[1em] \mathrm{Simpson: } \quad c & = -\dfrac{1}{90} g''''(z)\\[1em] {\mathrm{Gauss}\atop\mathrm{Legendre}}\quad c & = \dfrac{2^{2m+1}(m!)^4}{(2m+1)((2m)!)^3} g^{(2m)}(z)\quad \end{aligned}

Trasformazione per l’intervallo [a,b]:

N \geq \left(\dfrac{b-a}{2}\right)\left(\dfrac{|c|(b-a)}{2\epsilon}M_p\right)^{1/p}

2.5 Zeri di funzione

Capitolo di riferimento: Zeri di funzione, Iterazioni di punto fisso, Zeri di polinomi

Newton secanti Halley

\begin{aligned} x_{k+1} &= x_k - \dfrac{f(x_k)}{f'(x_k)},\\[1em] x_{k+1} &= x_k - \dfrac{(x_k-x_{k-1})f(x_k)}{f(x_k)-f(x_{k-1})},\\[1em] x_{k+1} &= x_k - \dfrac{f(x_k)f'(x_k)}{f'(x_k)^2-\dfrac{1}{2}f(x_k)f''(x_k)}, \end{aligned}

Verifica numerica dell’ordine con radice \alpha nota

\boxed{ p \approx\dfrac{\log\dfrac{\left\lvert x_{k+1}-\alpha\right\rvert}{\left\lvert x_{k+2}-\alpha\right\rvert} }{\log\dfrac{\left\lvert x_{k}-\alpha\right\rvert}{\left\lvert x_{k+1}-\alpha\right\rvert}}, }

Punto fisso: convergenza e ordine

\alpha=g(\alpha) e g lipschitziana su I=[\alpha-\rho,\alpha+\rho] con costante L<1: le iterate x_{k+1}=g(x_k) restano in I, convergono ad \alpha ed \alpha è l’unico punto fisso in I.

\left\lvert g(x)-g(x')\right\rvert \leq L\left\lvert x-x'\right\rvert, \qquad 0 \leq L < 1

L’ordine si legge dalle derivate in \alpha: se g'(\alpha)=\cdots=g^{(p-1)}(\alpha)=0 e g^{(p)}(\alpha)\neq 0 il metodo ha ordine p. Newton6 su radice semplice ha p=2, sulle secanti p=(1+\sqrt{5})/2\approx 1.618, Halley7 p=3.

2.6 Sistemi lineari: metodi iterativi

Capitolo di riferimento: Metodi Iterativi per Sistemi Lineari

Schema generale \bm{P}\bm{x}_{k+1} = \bm{b}+\bm{Q}\bm{x}_k

nota: mai invertire \bm{P}.

Tabella 2.1: Quadro riassuntivo metodi iterativi classici
Jacobi Gauss-Seidel SOR
\bm{P} \bm{D} \bm{D}-\bm{L} \dfrac{\bm{D}}{\omega}-\bm{L}
\bm{Q} \bm{L}+\bm{U} \bm{U} \dfrac{1-\omega}{\omega}\bm{D}+\bm{U}
\bm{I}-\bm{P}^{-1}\bm{A} \bm{D}^{-1}(\bm{L}+\bm{U}) (\bm{D}-\bm{L})^{-1}\bm{U} \left(\dfrac{\bm{D}}{\omega}-\bm{L}\right)^{-1}\left(\dfrac{1-\omega}{\omega}\bm{D}+\bm{U} \right)

Per il calcolo del raggio spettrale evitando l’inversione della matrice \bm{P}

0=\underbrace{|\bm{P}^{-1}\bm{Q}-\lambda\bm{I}|}_{p(\lambda)} =\underbrace{|\bm{P}^{-1}|}_{\neq 0}\underbrace{|\bm{Q}-\lambda\bm{P}|}_{q(\lambda)}

Criterio di convergenza

Lo schema \bm{P}\bm{x}_{k+1}=\bm{Q}\bm{x}_k+\bm{b} converge per ogni dato iniziale se e solo se

\varrho(\bm{P}^{-1}\bm{Q}) < 1 \qquad\Longleftrightarrow\qquad \varrho(\bm{Q}\bm{P}^{-1}) < 1

Condizione necessaria per SOR: \omega\in(0,2).

2.7 Equazioni differenziali ordinarie: Runge-Kutta

Capitolo di riferimento: Equazioni Differenziali Ordinarie, Metodi di Runge-Kutta

y_{k+1} = y_k + \sum_{i=1}^s b_i K_i, \qquad K_i = h f\Big( t_k+c_i h, y_k + \sum_{j=1}^s a_{ij} K_j \Big), \qquad \def\arraystretch{1.5} \begin{array}{c|c} \bm{c} & \bm{A} \\ \hline & \bm{b}^T \end{array}

\large \color{magenta} \left. \begin{aligned} &\color{DarkGreen}\underbrace{ \color{DarkRed}\overbrace{\color{blue}\overbrace{ \underbrace{\color{blue}\bm{A}\bm{1}=\bm{c}}_{\mathrm{consistenza}},\quad \color{blue}\bm{b}\cdot\bm{1}=1, }^{\mathrm{ordine 1}}\quad \color{DarkRed}\bm{b}\cdot\bm{c} = \dfrac{1}{2}}^{\mathrm{ordine 2}},\quad \color{DarkGreen}\bm{b}\cdot\bm{c}^2 = \dfrac{1}{3},\quad \bm{b}\cdot\bm{d} = \dfrac{1}{6} }_{\mathrm{ordine 3}} \\ & \color{magenta}\bm{b}\cdot\bm{c}^3 = \dfrac{1}{4},\quad (\bm{b}\otimes\bm{c})\cdot\bm{d} = \dfrac{1}{8},\quad \bm{e}\cdot\bm{c}^2 = \dfrac{1}{12},\quad \bm{e}\cdot\bm{d} = \dfrac{1}{24} \quad \end{aligned} \right\} _{\color{magenta}\mathrm{ordine 4}}

\bm{d} = \bm{A}\bm{c}, \qquad \bm{e} = \bm{A}^T\bm{b}, \qquad \bm{c}^2 = \bm{c}\otimes\bm{c},\qquad \bm{c}^3 = \bm{c}\otimes\bm{c}\otimes\bm{c}

Ordine di convergenza

Errore asintotico

  • y(x) soluzione analitica
  • y_k soluzione numerica in x_k

E_h = \max_{k=0}^n |y_k-y(x_k)| \approx C h^p.

stima dell’errore

\boxed{ \quad p \approx \dfrac{\log\dfrac{E_{h_1}}{E_{h_2}}} {\log\dfrac{h_1}{h_2}},\quad }

Stabilità assoluta

Capitolo di riferimento: Stabilità assoluta

Applicando un metodo a un passo all’equazione test z'=\lambda z si ottiene z_{n+1}=R(q)z_n con q=h\lambda; per un metodo di Runge-Kutta

R(q) = 1 + q\,\bm{b}^T(\bm{I}-q\bm{A})^{-1}\bm{1}, \qquad \mathcal{S} = \big\lbrace q : |R(q)|\leq 1 \big\rbrace

  • metodo esplicito \Longrightarrow R polinomio \Longrightarrow \mathcal{S} limitata
  • A-stabile: \mathrm{Re}(q)\leq 0 \Rightarrow q\in\mathcal{S} (solo metodi impliciti)
  • L-stabile: A-stabile e R(q)\to 0 per |q|\to\infty
Tabella 2.2: Funzioni di stabilità dei metodi a un passo.
R(q) \mathbb{R}^{-} A L
Eulero esplicito 1+q [-2,0] no no
RK2 1+q+\dfrac{q^2}{2} [-2,0] no no
RK4 \sum_{k\leq 4}\dfrac{q^k}{k!} \approx[-2.785,0] no no
Eulero implicito \dfrac{1}{1-q} (-\infty,0] sì sì
Trapezi \dfrac{1+q/2}{1-q/2} (-\infty,0] sì no

Problema stiff: rapporto \max|\mathrm{Re}\lambda_j| / \min|\mathrm{Re}\lambda_j| grande. Il passo di un metodo esplicito è allora dettato dalla stabilità (h\lambda_j\in\mathcal{S} per ogni j) e non dall’accuratezza.

  1. Carl Friedrich Gauss (1777-1855)↩︎

  2. Roger Penrose (1931)↩︎

  3. Henri Léon Lebesgue (1875-1941)↩︎

  4. Carl David Tolmé Runge (1856-1927)↩︎

  5. Pafnutij L’vovič Čebyšëv (1821-1894)↩︎

  6. Isaac Newton (1642-1727)↩︎

  7. Edmond Halley (1656-1742)↩︎