APPUNTI DI CALCOLO NUMERICO

[Questo file contiene le note da me redatte ad uso personale per un corso di Calcolo Numerico nell’ambito del programma di Dottorato di Ricerca a Firenze; il corso è stato tenuto dal Prof. Aldo Pasquali usando essenzialmente il libro di Monegato Fondamenti di calcolo numerico dal quale queste note sono prese; i programmi sono scritti in ANSI C].

§1.1. Eliminazione di Gauss

Supponiamo di voler risolvere il sistema lineare

Ax = b

ove A ∈ M_{n}(ℝ) è una matrice quadrata, b ∈ ℝ^{n} un vettore costante e x ∈ ℝ^{n} il vettore delle incognite del sistema; scriveremo sempre

A = ( \begin{matrix}a_{11} & a_{12} & ... & a_{1n} \\ a_{21} & a_{22} & ... & a_{2n} \\ ⋮ & ⋮ & ... & ⋮ \\ a_{n1} & a_{n2} & ... & a_{nn} \end{matrix} )

e b = (b_{1}, b_{2}, ..., b_{n})^{T}.

Supponiamo che A sia invertibile: questo implica, ad esempio, che non esiste una colonna di soli zeri in A, e quindi che, a meno di permutare le righe (il che vuol dire scambiare le equazioni del sistema), possiamo sempre supporre a_{ii} ≠ 0 per ogni i = 1, 2, ..., n.

Consideriamo per prima cosa il caso n = 1: questo è banale, poiché il sistema si riduce ad una sola equazione

a_{11}x_{1} = b_{1}

la cui soluzione è ovviamente x_{1} = b_{1}/a_{11} (a_{11} ≠ 0 perché la matrice (a_{11}) è appunto non singolare).

Nel caso n = 2 abbiamo

\{ \begin{matrix}a_{11}x_{1} + a_{12}x_{2} = b_{1} \\ a_{21}x_{1} + a_{22}x_{2} = b_{2} \end{matrix}

Allora, se a_{11} ≠ 0 possiamo sostituire la seconda equazione con una sua combinazione lineare del tipo:

a_{21}x_{1} + a_{22}x_{2} = b_{2} ↦ (a_{21}x_{1} + a_{22}x_{2} = b_{2}) − \frac{a_{21}}{a_{11}} (a_{11}x_{1} + a_{12}x_{2} = b_{1})

ottenendo cioè il sistema

\{ \begin{matrix}a_{11}x_{1} + a_{12}x_{2} = b_{1} \\ + (a_{22} − \frac{a_{21}}{a_{11}} a_{12}) x_{2} = b_{2} − \frac{a_{21}}{a_{11}} b_{1} \end{matrix}

2 §1.1. Eliminazione di Gauss

Ora la seconda equazione si risolve immediatamente come nel caso n = 1:

x_{2} = \frac{a_{11}b_{2} − a_{21}b_{1}}{a_{22}a_{11} − a_{21}a_{12}}

(di nuovo il denominatore è diverso da zero in quanto la matrice del sistema è invertibile: il denominatore è esattamente il determinante di questa matrice.)

Sostituendo questo valore di x_{2} nella prima equazione ci riconduciamo ad una equazione in una sola incognita x_{1} che sappiamo già risolvere.

Per n > 2 l’idea è di usare questo procedimento di eliminazione in modo ricorsivo: supponiamo cioè di saper risovere un sistema di n − 1 equazioni in n − 1 incognite con matrice invertibile: allora dato il sistema Ax = b possiamo trovare, dato che la prima colonna non può esser nulla in virtù dell’invertibilità della matrice A, un elemento a_{1i} ≠ 0 per qualche i ∈ {1, 2, ..., n}; scambiare di posto l’equazione i-esima con la prima non cambia ovviamente il sistema, e quindi possiamo farlo.

Ora consideriamo i numeri

m_{2} = −\frac{a_{21}}{a_{11}} , m_{3} = −\frac{a_{31}}{a_{11}} , ..., m_{n} = −\frac{a_{n1}}{a_{11}}

Ovviamente a_{i1} + m_{i}a_{11} = 0 e quindi sommando alla i-esima equazione la prima moltiplicata per m_{i} troviamo una equazione il cui coefficiente della variabile x_{1} è zero. Fatto questo il sistema diviene

\{ \begin{matrix}a_{11}x_{1} + a_{12}x_{2} + a_{13}x_{3} + ... + a_{1n}x_{n} = b_{1} \\ a′_{22}x_{2} + a′_{23}x_{3} + ... + a′_{2n}x_{n} = b′_{2} \\ a′_{32}x_{2} + a′_{33}x_{3} + ... + a′_{3n}x_{n} = b′_{3} \\ ................................................ \\ a′_{n2}x_{2} + a′_{n3}x_{3} + ... + a′_{nn}x_{n} = b′_{n} \end{matrix}

ove a′_{ij} = a_{ij} + m_{i}a_{1j} e b′_{i} = b_{i} + m_{i}.

Se ora escludiamo la prima equazione troviamo un sistema di n − 1 equazioni in n − 1 incognite, che abbiamo supposto di saper risolvere (la matrice di questo sistema è ancora non singolare perché il suo determinante non è cambiato se non nel segno, a causa dello scambio iniziale della riga i-sima con la prima, e per il fattore moltiplicativo a_{11}): dunque abbiamo la soluzione (ξ_{2}, ..., ξ_{n}) di questo sotto-sistema, e, sostituendo questi valori al posto delle incognite (x_{2}, ..., x_{n}) della prima equazione troviamo una equazione nella sola x_{1}, che possiamo risolvere come nel caso n = 1.

Infine, aggiungendo la soluzione di questa equazione alla soluzione del sotto-sistema (n − 1) × (n − 1) troviamo la soluzione del sistema di partenza.

Nello scambio di equazioni che si effettua per risolvere il sistema generale, l’idea è di sostituire uno zero con qualcosa di non nullo: sembra quindi naturale (ed in effetti garantisce maggior stabilità numerica) sostituire questo zero (cioè a_{11}) con il massimo numero (in valore assoluto) disponibile; poiché si tratta di scegliere un indice di riga i, sceglieremo l’equazione i-esima dove

|a_{i1}| = max_{1≤j≤n} |a_{j1}|

Questo modo di scegliere la riga da scambiare si chiama pivoting parziale.

Un altro metodo di effettuare la scelta è scegliere non solo per righe ma anche per colonne, cioè andare a cercare la coppia (i, j) tale che

|a_{ij}| = max_{1≤k,l≤n} |a_{kl}|

e scambiare la riga i-esima con la k-esima e la colonna j-esima con la l-esima: questo modo di procedere si chiama pivoting completo.

Volendo implementare l’algoritmo di eliminazione su una macchina numerica il pivoting completo, per quanto teoricamente migliore, può risultare assai più dispendioso del pivoting parziale.

Il condizionamento del problema della soluzione di un sistema di equazioni lineari dipende dalla matrice del sistema: supponiamo ad esempio di perturbare il vettore dei termini noti b di una quantità vettoriale δb e di risolvere il sistema Ax′ = b+δb in modo esatto: allora x′ = A^{−1}b + A^{−1}δb = x + δx, cioè δx = A^{−1}δb. Per stimare gli ordini di grandezza scegliamo una norma matriciale compatibile con una norma vettoriale: allora

∥δx∥ = ∥A^{−1}δb∥ ≤ ∥A^{−1}∥ ∥δb∥

cioè

\frac{∥δx∥}{∥δb∥} ≤ ∥A^{−1}∥

Ma ∥b∥ ≤ ∥A∥ ∥x∥, cioè (∥A∥ ∥x∥)^{−1} ≤ ∥b∥^{−1}, da cui

\frac{∥δx∥}{∥A∥ ∥x∥} ≤ ∥A^{−1}∥\frac{∥δb∥}{∥b∥}

Ne segue che

\frac{∥δx∥}{∥x∥} ≤ K(A)\frac{∥δb∥}{∥b∥}

ove K(A) = ∥A∥ ∥A^{−1}∥. Quindi il condizionamento del problema dipende dalla quantità K(A); supponendo di perturbare anche la matrice A, avremo

(A + δA)(x + δx) = b + δb ⟹ Aδx + δAx + δAδx = δb
⟹ (I + A^{−1}δA)δx = A^{−1}(δb − δAx)
⟹ δx = (I + A^{−1}δA)^{−1}A^{−1}(δb − δAx)

(supponendo che I + A^{−1}δA sia invertibile: poiché le matrici invertibili sono definite da un’equazione det A ≠ 0 sono un aperto denso nello spazio delle matrici, quindi esiste sempre una perturbazione dello stesso ordine di grandezza di δA che rende invertibile questa matrice.)

Ora vogliamo stimare

\frac{∥δx∥}{∥x∥} = \frac{∥(I + A^{−1}δA)^{−1}A^{−1}(δb − δAx)∥}{∥x∥} ≤ ∥(I + A^{−1}δA)^{−1}∥ ∥A^{−1}∥ (\frac{∥δb∥}{∥x∥} + ∥δA∥)

Supponiamo per questo che

∥δA∥ ≤ \frac{1}{∥A^{−1}∥}

4 §1.2. Decomposizione di Gauss

Allora

\frac{∥δx∥}{∥x∥} ≤\frac{∥A^{−1}∥}{∥I + A^{−1}δA∥} (\frac{∥δb∥}{∥x∥} + ∥δA∥) ≤ \frac{∥A^{−1}∥}{1 − ∥A^{−1}δA∥} (\frac{∥δb∥}{∥x∥} + ∥δA∥)
≤\frac{∥A^{−1}∥}{1 − ∥A^{−1}∥∥δA∥} (\frac{∥δb∥}{∥x∥} + ∥δA∥)

Ora teniamo conto che ∥b∥ ≤ ∥A∥∥x∥:

\frac{∥δx∥}{∥x∥} ≤\frac{∥A^{−1}∥}{1 − ∥A^{−1}∥∥δA∥} (∥A∥\frac{∥δb∥}{∥b∥} + ∥δA∥) ≤ \frac{∥A^{−1}∥∥A∥}{1 − ∥A^{−1}∥∥δA∥} (\frac{∥δb∥}{∥b∥} + \frac{∥δA∥}{∥A∥})
≤\frac{K(A)}{1 − K(A)\frac{∥δA∥}{∥A∥}} (\frac{∥δb∥}{∥b∥} + \frac{∥δA∥}{∥A∥})

Quindi il condizionamento del problema dipende essenzialmente da K(A): ad esempio vale il

1.1. Teorema. Se A è una matrice invertibile allora 1/K(A) è la distanza relativa (rispetto ad una norma compatibile con una vettoriale) di A dall’insieme delle matrici singolari, cioè

\frac{1}{K(A)} = min \{\frac{∥A − B∥}{∥A∥} | B singolare\}

§1.2. Decomposizione di Gauss

La tecnica dell’eliminazione consiste nel trasformare un sistema in un sistema la cui matrice sia triangolare: il legame fra la matrice di partenza e quella così ottenuta fa parte di un profondo risultato della teoria dei gruppi di Lie, noto come decomposizione di Iwasawa, ma che, nel caso da noi considerato di A ∈ GL_{n}(ℝ), risale appunto a Gauss.

L’idea è semplice: per risolvere il sistema Ax = b abbiamo manipolato la matrice A (ed il vettore b) sostituendo a delle sue righe combinazioni lineari di righe, e scambiando delle righe. Queste operazioni si possono esprimere in termini di moltiplicazione di matrici.

Ad esempio supponiamo di voler scambiare la riga i-esima di una matrice (quadrata) A con la riga j-esima di A: la matrice che ne risulta è P_{ij}A ove

P_{ij} = ( \begin{matrix}1 & 0 & 0 & ... & 0 & ... & 0 & ... & 0 \\ 0 & 1 & 0 & ... & 0 & ... & 0 & ... & 0 \\ 0 & 0 & 1 & ... & 0 & ... & 0 & ... & 0 \\ ⋮ & ⋮ & ⋮ & ... & ⋮ & ... & ⋮ & ... & ⋮ \\ 0 & 0 & 0 & ... & 0 & ... & 1 & ... & 0 \\ ⋮ & ⋮ & ⋮ & ... & ⋮ & ... & ⋮ & ... & ⋮ \\ 0 & 0 & 0 & ... & 1 & ... & 0 & ... & 0 \\ ⋮ & ⋮ & ⋮ & ... & ⋮ & ... & ⋮ & ... & ⋮ \\ 0 & 0 & 0 & ... & 0 & ... & 0 & ... & 1 \end{matrix} ) \begin{matrix}← i \\ ← j \end{matrix}

↑ ↑
i j

(le frecce indicano le colonne e le righe relative agli indici riportati.)

La proprietà che caratterizza la matrice P_{ij} è esattamente che quando si moltiplica a sinistra una matrice A per P_{ij} il risultato è la matrice P_{ij}A ottenuta da A scambiando la riga i-esima con la riga j-esima; se invece si moltiplica A per P_{ij} a destra, la matrice AP_{ij} che ne risulta è ottenuta da A scambiando la colonna i-esima con la colonna j-esima.

Allora l’operazione di scambiare l’equazione i-esima con la j-esima nel sistema equivale a considerare il nuovo sistema

P_{ij}A = P_{ij}b

le cui soluzioni coincidono con quelle di Ax = b.

Oltre a scambiare equazioni, abbiamo anche sostituito una equazione i-esima con se stessa sommata al multiplo scalare m di un’equazione j-esima: questa trasformazione si ottiene moltiplicando per la matrice

M_{ij}(m) = ( \begin{matrix}1 & 0 & 0 & ... & 0 & ... & 0 & ... & 0 \\ 0 & 1 & 0 & ... & 0 & ... & 0 & ... & 0 \\ 0 & 0 & 1 & ... & 0 & ... & 0 & ... & 0 \\ ⋮ & ⋮ & ⋮ & ... & ⋮ & ... & ⋮ & ... & ⋮ \\ 0 & 0 & 0 & ... & 1 & ... & 0 & ... & 0 \\ ⋮ & ⋮ & ⋮ & ... & ⋮ & ... & ⋮ & ... & ⋮ \\ 0 & 0 & 0 & ... & m & ... & 1 & ... & 0 \\ ⋮ & ⋮ & ⋮ & ... & ⋮ & ... & ⋮ & ... & ⋮ \\ 0 & 0 & 0 & ... & 0 & ... & 0 & ... & 1 \end{matrix} ) \begin{matrix}← i \\ ← j \end{matrix}

↑ ↑
i j

Questa matrice si può sinteticamente esprimere come

M_{ij}(m) = I + mE_{ij}

ove E_{ij} è la matrice che all’incrocio fra la riga i-esima e la riga j-esima ha un 1, mentre in tutte le altre entrate ha zero (cioè un elemento della base canonica dello spazio delle matrici).

Sia P_{ij} che M_{ij}(m) sono invertibili (per ogni i, j, m): infatti

P_{ij}P_{ij} = I_{n}

perché scambiando due volte i con j non si compie nessuna operazione, sicché P_{ij} = P_{ij}^{−1}; notiamo in particolare che P_{ii} = I_{n} per ogni i = 1, ..., n.

Inoltre

M_{ij}(m)M_{ij}(−m) = I_{n}

dato che sommando alla riga i-esima la j-esima moltiplicata per m e poi sommando a questo risultato la j-esima moltiplicata per −m otteniamo di nuovo la i-esima riga. Notiamo ancora che M_{ii}(0) = I_{n} per ciascun i = 1, ..., n.

Ora: il procedimento di Gauss può semplicemente interpretarsi come la moltiplicazione a sinistra di ambo i membri dell’equazione Ax = b per matrici P_{ij} e M_{ij}(m), ottenendo

GAx = Gb

ove abbiamo posto

G = M_{k_{n−1}l_{n−1}}(m_{k_{n−1}l_{n−1}})P_{(n−1)j_{n−1}}...M_{k_{2}l_{2}}(m_{k_{2}l_{2}})P_{2j_{2}}M_{k_{2}l_{2}}(m_{k_{2}l_{2}})P_{1j_{1}}

avendo compiuto n − 1 passi con opportuni indici i, j, k, l scelti col metodo del pivoting parziale ad esempio (ovviamente se non è necessario alcuno scambio al passo h poniamo ad esempio i_{h} = j_{h} = 1 in modo che P_{i_{h}j_{h}} = I_{n}).

La matrice T = GA è triangolare superiore per costruzione, ed è invertibile in quanto prodotto di matrici invertibili: la fattorizzazione A = G^{−1}T si dice decomposizione di Gauss.

Ovviamente per risolvere il sistema non è necessario determinare G ma solo le operazioni elementari cui corrispondono le matrici P_{ij} e M_{ij}(m) che, per moltiplicazione, danno luogo ad essa.

Stabiliamo ora un algoritmo che permetta di determinare le operazioni elementari da effettuare su una matrice A per ridurla in forma triangolare superiore T; non avremo bisogno di sapere quale sia G, ma semplicemente, ad ogni passo, ci basta sapere quale sia P_{ij} e quale sia M_{ij}(m).

Per determinare P_{ij} basta dare una coppia di numeri interi (i, j) e questo può farsi scrivendo un vettore p = (p_{1}, ..., p_{n−1}) con n − 1 componenti intere (tante quanti i passi da compiere), tale che p_{i} = j. Ad esempio, se al passo h-esimo si devono scambiare la riga h e la riga j allora p(h) = j.

Per determinare M_{ij}(m_{ij}) basta dare il numero m_{ij} e la sua posizione (i, j); poiché queste posizioni non possono coincidere per matrici diverse, possiamo metterli tutti in una matrice triangolare inferiore (i > j) nel modo seguente

( \begin{matrix}0 & 0 & 0 & ... & 0 \\ m_{21} & 0 & 0 & ... & 0 \\ m_{31} & m_{32} & 0 & ... & 0 \\ m_{41} & m_{42} & m_{43} & ... & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ m_{n1} & m_{n2} & m_{n3} & ... & 0 \end{matrix} )

Notiamo infine che possiamo calcolare il determinante di A durante il processo di decomposizione di Gauss, dato che

det G = det [M_{k_{n−1}l_{n−1}}(m_{k_{n−1}l_{n−1}})P_{(n−1)j_{n−1}}...M_{k_{2}l_{2}}(m_{k_{2}l_{2}})P_{2j_{2}}M_{k_{2}l_{2}}(m_{k_{2}l_{2}})P_{1j_{1}}]

Ma det M_{ij}(m) = 1 e det P_{ij} = ±1 ove il segno + vale se e solo se i = j; inoltre

det G det A = det T

e, dato che T è triangolare, il suo determinante è il prodotto degli elementi diagonali t_{11}, ..., t_{nn}; morale:

det A = (−1)^{s} ∏_{i=1}^{n} t_{ii}

ove s è il numero complessivo di scambi effettuati in tutto il processo.

Possiamo finalmente dare un algoritmo per la fattorizzazione di Gauss, che economizzerà gli spazi di memorizzazione nel modo seguente: la matrice A, che è data in input, verrà sovrascritta nel seguente modo

A ↦ ( \begin{matrix}t_{11} & t_{12} & t_{13} & ... & t_{1n} \\ m_{21} & t_{22} & t_{23} & ... & t_{2n} \\ m_{31} & m_{32} & t_{33} & ... & t_{3n} \\ m_{41} & m_{42} & m_{43} & ... & t_{4n} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ m_{n1} & m_{n2} & m_{n3} & ... & t_{nn} \end{matrix} )

ove m_{ij} sono i coefficienti delle matrici M_{ij}(m_{ij}) collocati nelle posizioni giuste e t_{ij} sono i coefficienti della matrice triangolare T.

In ANSI C il tutto si esprime come^{1}:

�define n /* numero di equazioni */ /*

Algoritmo di fattorizzazione di Gauss per una matrice A usando il pivoting parziale come scelta per gli scambi di riga; la matrice T ed i dati delle matrici M_{ij}(m) sono memorizzati al posto di A, ed il vettore dei pivot contiene gli scambi di riga effettuati nei vari passi; la funzione ritorna il determinante di A: se �e zero, il risultato in A �e da considerarsi indeterminato.

*/ float factor(

float A[n][n], /* matrice da fattorizzare */ int pivot[n-1] /* vettore degli scambi */ )

� float a; /* variabile ausiliaria */ float d=1; /* determinante */ int i, j, k; /* indici di ciclo */ int i0; /* indice per gli scambi */ for(k=0; k�n-1; k++)�

/* AbsMax calcola max_{k≤i≤n} |a_{ik}| e torna in i0 il pi�u piccolo degli indici per i quali il massimo �e assunto da |a_{ik}|

*/

a = AbsMax(A, &i0, k); /* i0 passato per indirizzo */ pivot[k] = i0; if(a==0) /* esce con un errore: matrice singolare */

return 0; if(i0!=k)�

/* scambia la riga k-esima con la riga i0-esima */ for(i=k; i�n; i++)� a = A[k][i]; A[k][i] = A[i0][i]; A[i0][i] = a;

1Si rammenti che: le matrici in C partono dall’indice 0 e non 1; è necessario dichiarare almeno la dimensione delle colonne di una matrice quadrata, se questa viene passata come parametro in una procedura; le matrici sono passate come parametri per mezzo del loro indirizzo, quindi se la funzione ne altera il contenuto questo cambiamento è permanente; per alterare il contenuto di un parametro intero o reale è necessario passarlo per indirizzo dichiarandolo come puntatore nell’intestazione della funzione.

�;

/* cambia segno al determinante perch�e ha

effettuato uno scambio */ d = - d;

�;

for(i=k+1; i�n; i++)� /* ora pone a_{ik} = −a_{ik}/a_{kk}(= m_{ik}) */ A[i][k] /= -A[k][k]; /* sarebbe m_{ik} */ /* e pone a_{ij} = a_{ij} + m_{ik}a_{kj} */ for(j=k+1; j�n; j++)

A[i][j] += A[i][k] * A[k][j]; �;

/* d = d · t_{kk} */

d *= A[k][k]; /* sarebbe t_{kk} */ �; return d * A[n-1][n-1];

�

Resta da scrivere la procedura per il calcolo del massimo valore assoluto degli elementi della colonna k-esima:

�include �math.h� /*

AbsMax calcola max_{k≤i≤n} |a_{ik}| e torna in i0 il pi�u piccolo degli indici per i quali il massimo �e assunto da |a_{ik}|

*/ float AbsMax(

float A[n][n], /* matrice */ int *i0, /* indice da trovare: passato per indirizzo */ int k /* riga da cui partire */ )

� float a = 0; /* conterr�a il risultato */ int i; /* indice di ciclo */ /* cerca il massimo */ for(i=k; i�n; i++)

if(a�fabs(A[i][k]))�/* fabs(x) = |x| */ a = fabs(A[i][k]); /* salva il risultato */ *i0 = i; /* annota l’indice */

�; return a; /* ritorna il massimo modulo */

�

Ovviamente per efficienza sarebbe bene includere questa funzione direttamente dentro factor. Si può usare questa routine in un programma del tipo:

�define n /* dimensione della matrice */

�include �stdio.h� �include �factor.h� void main()�

float A[n][n]= �...�; int pivot[n-1]; int i, j; printf(�Determinante =%g, matrice T=�n�, factor(A, p)); /* ora stampa la matrice per righe */ for(i=0; i�n; i++)�

printf(�� �); for(j=0; j�n; j++)

printf(�%g�t�, A[i][j]); printf(��b��n�);

�; �

ove il file factor.h contiene le dichiarazioni delle funzioni factor e AbsMax.

A questo punto, data una matrice A sappiamo trovare il suo fattore triangolare T, le matrici di scambio P (codificate nel vettore p e la matrice G, codificata nella matrice triangolare strettamente inferiore i cui coefficienti sono m_{ij}.

Per risolvere effettivamente il sistema Ax = b, dobbiamo considerare GAx = Gb cioè T x = Gb; si tratta quindi di calcolare il vettore Gb e di risolvere un sistema in forma triangolare.

Usando la funzione factor possiamo scrivere una funzione solve che svolga queste operazioni:

/* Algoritmo di soluzione di un sistema lineare che utilizza la fattorizzazione di Gauss: suppone di avere nella matrice A il risultato di factor() e nel vettore p il vettore dei pivot prodotto da factor; b deve contenere il vettore dei termini noti del sistema e, in uscita, conterr�a il vettore con la soluzione.

*/ void solve(

float A[n][n], /* matrice fattorizzata del sistema */ int pivot[n-1], /* vettore degli scambi */ float b[n] /* termine noto del sistema */ )

� float a; /* variabile ausiliaria */ int i, j, k; /* indici di ciclo */ /* ora calcola Gb e lo pone in b */ for(k=0; k�n-1; k++)�

j = pivot[k]; /* riga con cui va scambiata la k-esima */ if( j!=k )�

10 §1.3. Raffinamento iterativo

/* scambia b[j] con b[k] */ a = b[j]; b[j] = b[k]; b[k] = a; �; /* ora somma per i multipli degli m_{ik} */ for(i=k+1; i�n; i++)

b[i] += A[i][k] * b[k]; �; b[n-1] /= A[n-1][n-1]; /*

ora risolve ’a cascata’ le equazioni lineari risultanti partendo dall’ultima e sostituendo i valori trovati nelle incognite delle precedenti

*/

for(i=n-2; i�=0; i--)�

/* b_{i} = (b_{i} − ∑_{j=i+1}^{n} a_{ij}b_{j}) /a_{ii} */

a = b[i];

for(j=i+1; j�n; j++)

a -= A[i][j] * b[j]; b[i] = a / A[i][i];

�; �

Una ulteriore osservazione che possiamo fare sulla fattorizzazione di Gauss è che possiamo raggruppare le matrici di scambio P_{ij} e le matrici di somma M_{ij}(m) nel fattore G: precisamente, dato che P_{ij} = P_{ij}^{−1}, possiamo scrivere

P_{ij}X = (P_{ij}XP_{ij})P_{ij}

per ogni matrice X; questo consente di “spostare verso destra” le matrici di scambio nella fattorizzazione sostituendo alle M_{ij}(m) le loro coniugate per P_{kl}.

Questo non cambia la natura della fattorizzazione, dato che

P_{ij}M_{hk}(m)P_{ij} = P_{ij}(I + mE_{hk})P_{ij} = (P_{ij} + mP_{ij}E_{hk})P_{ij} = I + mP_{ij}E_{hk}P_{ij}

Ma la matrice P_{ij}E_{hk}P_{ij} è ancora della forma E_{rs} e quindi coniugare M_{hk}(m) per P_{ij} dà ancora luogo ad una matrice del tipo M_{rs}(m).

Dunque possiamo scrivere G nella forma

G = M_{n−1}(m_{n−1})...M_{2}(m_{2})M_{1}(m_{1})P_{n−1}...P_{2}P_{1}

ove {m_{1}, ..., m_{n−1}} = {m_{ij}}_{i>j} sono esattamente gli stessi coefficienti di prima, a meno di riordinamento; inoltre le matrici P_{i} sono esattamente le P_{ij} precedenti nello stesso ordine.

Possiamo quindi scrivere G = MP ed ottenere

P A = MT

La soluzione del sistema Ax = b equivale quindi a quella dei due sistemi

\{ \begin{matrix}My = P b \\ T x = y \end{matrix}

ambedue in forma triangolare (il primo inferiore, il secondo superiore).

§1.3. Raffinamento iterativo

Per limitare la propagazione dell’errore dovuta all’instabilità del problema nel metodo di Gauss, quando il problema non è eccessivamente mal condizionato, si può utilizzare l’algoritmo di raffinamento iterativo.

Se Ax = b è il sistema che vogliamo risolvere, possiamo codificare questo problema in \overline{A}\overline{x} = \overline{b}, ove \overline{A} rappresenta la matrice A scritta usando i numeri disponibili nella macchina; ora supponiamo che la soluzione di Ax = b fornita dall’algoritmo di Gauss sia x′, e che la soluzione esatta sia in realtà x′ + y′: allora

A(x′ + y′) = b

e quindi, per r′ = b − Ax′, il sistema

Ay′ = r′

è soddisfatto; possiamo risolverlo con Gauss (usando la stessa decomposizione GA = T che abbiamo adoperato per risolvere Ax = b) e troveremo non la soluzione y′ ma una sua approssimazione \overline{y′}.

Ora poniamo

x′′ = x′ + \overline{y′}

e ripetiamo il procedimento: cioè consideriamo la soluzione esatta x′′ + y′′ e così via.

In presenza di un buon condizionamento del problema, cioè se a una perturbazione dei dati iniziali A e b corrisponde una perturbazione della soluzione trovata dello stesso ordine di grandezza, il procedimento converge e consente di trovare approssimazioni migliori alla soluzione cercata.

Poiché il passaggio cruciale del metodo consiste nel calcolo di una differenza b−Ax^{(i)}, per limitare la cancellazione numerica valuteremo nell’algoritmo questa differenza in doppia precisione.

/* Raffinamento iterativo per migliorare l’approssimazione delle soluzioni di un sistema lineare risolto con factor/ solve; compie al pi�u kmax iterazioni e controlla se dopo queste la precisione della macchina �e stata raggiunta o meno; se questo non �e il caso, torna un valore non nullo.

*/

�define epsilon /* precisione della macchina */

int refine( float A[n][n], /* matrice di partenza */ float GT[n][n], /* matrice dopo la fattorizzazione */ float b[n], /* vettore dei termini costanti */ float pivot[n-1], /* vettore degli scambi */ float x[n], /* soluzione trovata da solve */ int kmax /* massimo numero di iterazioni */ )

� double d; /* usato nel calcolo in doppia precisione */ float y[n]; /* soluzione trovata da solve */

12 §1.3. Raffinamento iterativo

float in, xn, xn1, yn; /* norme di vettori */ int i, j, k; /* indici di loop */ /* nella prima iterazione usa x come dato iniziale: ora

ne calcola la norma iniziale */ in = norma(x); /* ciclo principale */ for(i=0; i�kmax; i++)�

/* ora calcola y = b − Ax (doppia precisione) */ for(j=0; j�n; j++)�

d = (double)b[j]; for(k=0; k�n; k++)

d -= (double)(A[j][k]*x[k]); /* ora pone y[i]=d (perdendo precisione) */ y[j] = (float)d;

�; /* ora risolve Ay = r */ solve(GT, pivot, y); /* calcola la norma della soluzione */ yn = norma(y); /* ora calcola ∥x∥ */ xn = norma(x); /* ora somma questa soluzione y a x ottenendo x’ */ for(j=0; j�; j++)

x[j] += y[j]; /* e ne calcola la norma */ xn1 = norma(x); /*

ora decide se continuare: se in∥x∥ ≤ ∥y∥ ∥x′∥ allora finisce con un errore; se ∥y∥ ≤ ε∥x′∥ allora finisce con un successo; altrimenti continua l’iterazione

*/ if( in*xn1 �= yn*xn ) return 1; if( yn �= epsilon*xn1 ) return 0; /* cambia la norma iniziale in quella di y */ in = yn;

�; /* se arriva qui ha superato il numero

di iterazioni consentite */ return 1;

�

Come norma possiamo ad esempio considerare quella euclidea:

�include �math.h�

/*

Calcola la norma euclidea di un vettore con n elementi */ float norma(

float v[n] /* vettore del quale calcolare la norma */

) �

float x=0; int i; for(i=0; i�n; i++)

x += v[i]*v[i]; x = sqrt(x);

�

Ad esempio su una macchina con numeri a 32 bit il sistema

( \begin{matrix}5 & 7 & 6 & 5 \\ 7 & 10 & 8 & 7 \\ 6 & 8 & 10 & 9 \\ 5 & 7 & 9 & 10 \end{matrix} ) ( \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \end{matrix} ) = ( \begin{matrix}23 \\ 32 \\ 33 \\ 31 \end{matrix} )

ha fornito col metodo di Gauss la soluzione

( \begin{matrix}0.999993 \\ 1 \\ 1 \\ 0.999999 \end{matrix} )

e, dopo un solo passo del raffinamento (con epsilon=.000001), la soluzione esatta

( \begin{matrix}1 \\ 1 \\ 1 \\ 1 \end{matrix} )

§1.4. Inverse di matrici

Nota la fattorizzazione di Gauss GA = T di una matrice invertibile A, possiamo calcolare A^{−1} come A^{−1} = T^{−1}G.

Questo fornisce una soluzione ovviamente approssimata che possiamo provare a raffinare se il problema non è mal condizionato: l’idea è che se B è una approssimazione di A^{−1}, se definiamo E = A^{−1} − B, allora A(B + E) = I e quindi, ponendo

R = AE = I − AB

otteniamo una migliore approssimazione

B′ = B + BR

Iterando questo procedimento su B′ si ottiene una successione {B^{(i)}} che converge alla soluzione esatta se

∥I − AB∥ < 1

Anche qui, nel calcolo di R, è dietro l’angolo il pericolo della cancellazione numerica.

14 §1.5. Metodi iterativi

§1.5. Metodi iterativi

Se esigenze di spazio non consentono una agevole manipolazione della matrice A, per risolvere il sistema lineare si può usare un metodo di “approssimazioni successive” che alteri ad ogni passo solo il vettore soluzione; questo si chiama metodo iterativo.

Una idea per applicare questo metodo è decomporre A in somma A = D + C; allora il sistema Ax = b diviene

Dx = −Cx + b

Allora, partendo da un dato iniziale x^{(0)} possiamo risolvere il sistema Dx^{(1)} = Cx^{(0)} +b, e poi il sistema Dx^{(2)} = Cx^{(1)} + b, e così via fino ad una soluzione x^{(n)} soddisfacente.

Ovviamente la scelta della decomposizione è cruciale: D sarà assunta con le seguenti proprietà:

1) det D ≠ 0

2) D è triangolare (o se possibile diagonale)

3) ∥I − D^{−1}A∥ < 1

In effetti 3) garantisce la convergenza della successione {x^{(i)}} alla soluzione del sistema.

Infatti abbiamo Dx^{(i+1)} = −Cx^{(i)} + b, cioè

x^{(i+1)} = −D^{−1}Cx^{(i)} + D^{−1}b

Poniamo ora B = −D^{−1}C = I − D^{−1}A: allora

x^{(i+1)} = Bx^{(i)} + D^{−1}b

Se e^{(i+1)} è l’errore assoluto relativo alla soluzione approssimata x^{(i+1)}, cioè

e^{(k+1)} = x − x^{(i+1)} = (Bx + D^{−1}b) − (Bx^{(i)} + D^{−1}b) = B(x − x^{(i)}) = Be^{(i)}

allora possiamo determinare la successione degli errori come

e^{(0)} = x − x^{(0)}
e^{(1)} = Be^{(0)}
e^{(2)} = Be^{(1)} = B^{2}e^{(0)}
e^{(3)} = Be^{(2)} = B^{3}e^{(0)}
...............................
e^{(i+1)} = Be^{(i)} = B^{i+1}e^{(0)}
...............................

La convergenza della successione {x^{(i)}} alla soluzione esatta equivale all’essere la successione {e^{(i)}} infinitesima, cioè che sia infinitesima la successione di matrici {B^{i}}: un criterio necessario e sufficiente per questo è che il raggio spettrale ρ(B) sia minore di 1.

Dunque il procedimento iterativo converge se e solo se

ρ(B) = ρ(−D^{−1}C) = ρ(I − D^{−1}A) < 1

Poiché una norma matriciale (compatibile con una norma vettoriale data) verifica la ρ(B) < ∥B∥, la 3) precedente implica la convergenza.

Una decomposizione A = D + C molto semplice è quella usata nel metodo di Jacobi: si sceglie come D la matrice diagonale i cui elementi diagonali sono quelli di A, e come C ovviamente A − D.

Ovviamente dobbiamo avere che a_{11} ≠ 0, a_{2} ≠ 0, ..., a_{nn} ≠ 0 altrimenti D risulta singolare.

Comunque se qualche elemento diagonale di A dovesse essere nullo, possiamo sempre scambiare due righe in A in modo che l’elemento diagonale in questione risulti diverso da zero (se questo non è possibile allora A possiede una colonna di zeri, e quindi è singolare, il che è escluso per ipotesi).

Quindi, a meno di scambiare le equazioni nel sistema, possiamo sempre supporre che la decomposizione di Jacobi sia possibile.

La convergenza è legata alla relazione 3): usando ad esempio la norma uniforme, questo vuol dire che:

∥I − D^{−1}A∥ < 1

Ma D^{−1} ha come elementi diagonali gli inversi degli elementi diagonali di A, e quindi

D^{−1}A = ( \begin{matrix}1 & a_{11}^{−1}a_{12} & a_{11}^{−1}a_{13} & ... & a_{11}^{−1}a_{1n} \\ a_{22}^{−1}a_{21} & 1 & a_{22}^{−1}a_{23} & ... & a_{22}^{−1}a_{2n} \\ a_{33}^{−1}a_{31} & a_{33}^{−1}a_{32} & 1 & ... & a_{33}^{−1}a_{3n} \\ ⋮ & ⋮ & ⋮ & ... & ⋮ \\ a_{nn}^{−1}a_{n1} & a_{nn}^{−1}a_{n2} & a_{nn}^{−1}a_{n3} & ... & 1 \end{matrix} )

sicché

∥I − D^{−1}A∥_{∞} = max_{1≤i≤n} ∑_{{j=1 \atop j≠−1}}^{n} |\frac{a_{ij}}{a_{ii}} |

e quindi la 3) si scrive semplicemente come

∑_{{j=1 \atop j≠−1}}^{n} |a_{ij}| < |a_{ii}|

per i = 1, ..., n.

Programmiamo ora il metodo di Jacobi: è molto semplice, dato che il sistema Ax = b si può scrivere come Dx = −Cx + b e D è facilmente invertibile essendo diagonale: x = −D^{−1}Cx + D^{−1}b equivale a

x_{i} = \frac{b_{i} − ∑_{{j=1 \atop j≠i}}^{n} a_{ij}x_{j}}{a_{ii}}

per i = 1, ..., n; partendo da una approssimazione x^{(i)} si calcola la x^{(i+1)} usando queste equazioni ove nei membri a destra si pongono le x^{(i)}.

/* Risoluzione del sistema Ax=b col metodo di Jacobi: richiede il numero massimo di iterazioni e restituisce nel vettore x una soluzione approssimata migliore di quella di partenza */ void jacobi(

float A[n][n], /* matrice del sistema */ float x[n], /* vettore delle incognite */ float b[n], /* vettore dei termini noti */ int kmax /* numero di iterazioni */ )

� float x1[n]; /* contiene copia dei vecchi dati */ float a; /* variabile ausiliaria */ int i, j, k; /* indici di ciclo */ /* per prima cosa ordina le equazioni in modo che

gli elementi diagonali di A siano tutti non nulli */ riordina(A, b); /* ora procede col metodo iterativo */ for(k=0; k�kmax; k++)�

/* copia x[] in x1[] */

for(i=0; i�n; i++) x1[i] = x[i]; for(i=0; i�n; i++)�

/* pone x_{i} = (b_{i} − ∑_{j≠i} a_{ij}x′_{j})/a_{ii} */

a = b[i];

for(j=0; j!=i && j�n; j++) a -= A[i][j]*x1[j];

x[i] = a/A[i][i]; �;

�; �

Resta solo da scrivere la funzione di scambio delle righe: si tratta di scambiare le righe della matrice A (e del vettore b) in modo da ottenere elementi diagonali tutti non nulli: una possibile soluzione (non molto elegante visto che utilizza una istruzione goto è la seguente funzione

/* Se non tutti gli elementi diagonali di A sono diversi da zero scambia le righe in A e b fino a soddisfare questo requisito; se non �e possibile (det A = 0) torna un numero diverso da zero

*/ int riordina(

float A[n][n], /* matrice da riordinare per righe */ float b[n] /* vettore dei termini noti */ )

� float a; /* variabile ausiliaria */ int i, j, k; /* indici di ciclo */ /* ciclo infinito */ for(i=0;; i++)�

/* cerca una riga il cui termine diagonale sia nullo */ for(; i�n; i++)

if(A[i][i]==0)�

/* ora i contiene l’indice della riga da scambiare: cerca una riga il cui termine A[i][i] sia non nullo fra quelle seguenti */

for(j=i+1; j�n; j++) if(A[j][i]!=0)

goto scambia; /* se arriva qui la riga cercata �e fra le precedenti: ne cerca una j-esima in modo che A[i][j]!=0 */

for(j=1; j�i; j++)

if(A[j][i]!=0 && A[i][j]!=0) goto scambia;

/* se arriva qui esiste una colonna di zeri */ return 1;

�; /* se arriva qui la matrice �e a posto */ return 0;

scambia: /* scambia la riga i con la riga j */ for(k=0; k�n; k++)�

a = A[j][k];

A[j][k] = A[i][k]; A[i][k] = a;

�; a = b[j]; b[j] = b[i]; b[i] = a;

�; /* ripete il ciclo infinito: ricomincia la ricerca */ �

La programmazione del metodo di Jacobi suggerisce un modo di migliorarlo: infatti abbiamo dovuto copiare il vettore x prima di riscriverlo, per poterne usare i valori precedenti; un miglioramento si avrebbe considerando i nuovi valori delle componenti di x appena calcolati: cioè usare la formula

x_{i} = \frac{b_{i} − ∑_{j=1}^{i−1} a_{ij}x_{j} − ∑_{j=i+1}^{n} a_{ij}x′_{j}}{a_{ii}}

ove x′ rappresenta la precedente approssimazione.

Questa modifica del metodo di Jacobi è nota come metodo di Gauss-Seidel; nel programmarlo inseriamo nell’algoritmo anche un controllo della precisione raggiunta nei passi stabiliti

/* Risoluzione del sistema Ax=b col metodo di Gauss-Seidel: richiede il numero massimo di iterazioni e la precisione desiderata nell’approssimazione della soluzione; restituisce nel vettore x una soluzione approssimata migliore di quella di partenza; se nel numero di passi richiesto �e raggiunta la precisione desiderata, torna 0, altrimenti torna 1

*/ int gauss seidel(

float A[n][n], /* matrice del sistema */ float x[n], /* vettore delle incognite */ float b[n], /* vettore dei termini noti */ int kmax, /* numero di iterazioni */ float prec /* precisione desiderata */ )

� float err; /* usato nel calcolo dell’errore */ float a; /* variabile ausiliaria */ int i, j, k; /* indici di ciclo */ /* per prima cosa ordina le equazioni in modo che

gli elementi diagonali di A siano tutti non nulli */ riordina(A, b); /* ora procede col metodo iterativo */ for(k=0; k�kmax; k++)�

/* calcola a = (b_{1} − ∑_{j=2}^{n} a_{1j}x′_{j})/a_{11} */

a = b[0];

for(j=1 /* parte dal 2o elemento */; j�n; j++)

a -= A[0][j]*x[j]; a /= A[0][0]; err = fabs(a-x[0]); x[0] = a; /* nuovo valore di x[0] */ for(i=1; i�n; i++)�

/* ora calcola a = (b_{i} − ∑_{j=1}^{i−1} a_{ij}x_{j} + ∑_{j=i+1}^{n} a_{ij}x′_{j})/a_{ii} */ a = b[i]; for(j=0; j!=i && j�n; j++) a -= A[i][j]*x[j];

a /= A[i][i];

/* sceglie l’errore pi�u grave commesso */

if(err�fabs(a-x[i]))

err = fabs(a-x[i]); x[i] = a; /* nuovo valore di x[i] */

�;

if(err�prec*norma(x)) return 0;

�; return 1;

�

Come si vede non c’è più bisogno di memorizzare il vettore x′ e quindi si ha anche un risparmio di spazio.

Il metodo di Gauss–Seidel corrisponde ad usare una decomposizione A = D + C ove D è la matrice triangolare inferiore estratta da A (supponendo che i suoi elementi diagonali siano tutti non nulli: altrimenti A andrebbe preventivamente riordinata per righe):

D = ( \begin{matrix}a_{11} & 0 & ... & 0 \\ a_{21} & a_{22} & ... & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ a_{n1} & a_{n2} & ... & a_{nn} \end{matrix} )

Di nuovo la condizione

∑_{{j=1 \atop j≠−1}}^{n} |a_{ij}| < |a_{ii}|

(per i = 1, ..., n) garantisce la convergenza dell’algoritmo.

Quindi per le matrici dominanti (per righe o per colonne) la convergenza è assicurata; anche nel caso delle matrici irriducibili abbiamo certamente convergenza dei metodi di Jacobi e Gauss–Seidel.

§1.6. Sovrarilassamento

Il metodo di Gauss–Seidel si riduce sostanzialmente alla formula iterativa

x_{i} = \frac{b_{i} − ∑_{j=1}^{i−1} a_{ij}x_{j} − ∑_{j=i+1}^{n} a_{ij}x′_{j}}{a_{ii}}

ove x′ rappresenta la precedente approssimazione; da questa, sottraendo ad ambo i membri il termine x′_{i}, definiamo

r′_{i} = x_{i} − x′_{i} = \frac{b_{i} − ∑_{j=1}^{i−1} a_{ij}x_{j} − ∑_{j=i}^{n} a_{ij}x′_{j}}{a_{ii}}

Questo valore r′_{i} può essere considerato una correzione di x′_{i} per ottenere x_{i}, dato che x_{i} = x′_{i} + r′_{i}.

Possiamo pensare di migliorare questa correzione riscalandola di un fattore ω:

x_{i} = x′_{i} + ωr′_{i}

Chiaramente il problema sta nello scegliere ω in modo da accelerare la convergenza della successione delle approssimanti alla soluzione esatta.

Nell’introduzione di ω consiste il metodo del rilassamento (di sottorilassamento se 0 < ω < 1 e di sovrarilassamento se 1 < ω); per ω = 1 ritroviamo il metodo di Gauss–Seidel.

Osserviamo che, spezzando A in somma L + D + U ove L è la matrice triangolare strettamente inferiore estratta da A, A la matrice diagonale estratta da A e U la matrice 20 §1.6. Sovrarilassamento

triangolare strettamente superiore estratta da A, troviamo che la formula per il metodo di Gauss–Seidel si scrive in forma vettoriale come

x = D^{−1}(b − Lx − Ux′)

Inoltre

r′ = D^{−1}(b − Lx − (D + U)x′)

e

(D + ωL)x = ((1 − ω)D − ωU)x′ + ωb

Dunque, la convergenza del metodo di rilassamento equivale alla condizione

ρ((D + ωL)^{−1}((1 − ω)D − ωU)) < 1

Il raggio spettrale misura in questo caso la rapidità della convergenza.

La ricerca del valore di ω ottimale rispetto alla convergenza dipende dal problema: per certe classi di matrici si riesce effettivamente a determinare il valore ottimo per ω.

Per le matrici simmetriche definite positive ogni valore ω ∈ (0, 2) il metodo converge; per ω ∉ (0, 2) il raggio spettrale precedente è almeno uno, e quindi il metodo certamente non converge.

La seguente funzione realizza il metodo di Gauss–Seidel tenendo conto di un parametro di rilassamento:

/* Risoluzione del sistema Ax=b col metodo di rilassamento: richiede il numero massimo di iterazioni e la precisione desiderata nell’approssimazione della soluzione; restituisce nel vettore x una soluzione approssimata migliore di quella di partenza; se nel numero di passi richiesto �e raggiunta la precisione desiderata, torna 0, altrimenti torna 1

*/ int gauss seidel(

float A[n][n], /* matrice del sistema */ float x[n], /* vettore delle incognite */ float b[n], /* vettore dei termini noti */ float omega, /* parametro di rilassamento */ int kmax, /* numero di iterazioni */ float prec /* precisione desiderata */ )

� float err; /* usato nel calcolo dell’errore */ float a; /* variabile ausiliaria */ int i, j, k; /* indici di ciclo */ /* per sicurezza controlla omega */ if(omega�=0 �� omega�=2)

omega = 1; /* ordina le equazioni in modo che gli elementi diagonali

di A siano tutti non nulli */ riordina(A, b); /* ora procede col metodo iterativo */ for(k=0; k�kmax; k++)�

/* calcola a = (b_{1} − ∑_{j=2}^{n} a_{1j}x′_{j})/a_{11} */ a = b[0]; for(j=0 /* parte dal 1o elemento */; j�n; j++)

a -= A[0][j]*x[j]; a *= omega/A[0][0]; err = fabs(a-x[0]); x[0] += a; /* nuovo valore di x[0] */ for(i=1; i�n; i++)�

/* ora calcola a = (b_{i} − ∑_{j=1}^{i−1} a_{ij}x_{j} + ∑_{j=i}^{n} a_{ij}x′_{j})/a_{ii} */ a = b[i]; for(j=0; j�n; j++)

a -= A[i][j]*x[j]; a *= omega/A[i][i]; /* sceglie l’errore pi�u grave commesso */ if(err�fabs(a-x[i]))

err = fabs(a-x[i]); x[i] += a; /* nuovo valore di x[i] */

�;

if(err�prec*norma(x)) return 0;

�;

return 1; �

§1.7. Metodo delle potenze

Supponiamo di avere una matrice reale A e di essere interessati ad i suoi autovalori λ_{1}, ..., λ_{n}; supporremo che questi siano ordinabili nel seguente modo:

|λ_{1}| > |λ_{2}| ≥ |λ_{3}| ≥ ⋯ ≥ |λ_{n}|

Supponiamo inoltre che esista una base di autovettori B = {x_{1}, ..., x_{n}} ove x_{i} è autovettore di autovalore λ_{i}.

Ora scegliamo un vettore v_{0} che sia tale che la sua decomposizione in coordinate rispetto alla base B

v_{0} = a_{1}x_{1} + ⋯ a_{n}x_{n}

abbia la prima coordinata non nulla: a_{1} ≠ 0.

Poiché A è reale, tale è il suo polinomio caratteristico, e le sue radici appaiono quindi o reali o in coppie complesse coniugate: dato che |λ_{1}| è più grande (e quindi diverso) da ogni altro modulo di autovalore, λ_{1} è reale, e quindi possiamo supporre che il suo autovettore x_{1} abbia componenti reali (rispetto alla base canonica).

22 §1.7. Metodo delle potenze

Ora consideriamo v_{1} = Av_{0}:

v_{1} =a_{1}Ax_{1} + a_{2}Ax_{2} + ⋯ + a_{n}Ax_{n} = a_{1}λ_{1}x_{1} + a_{2}λ_{2}x_{2} + ⋯ + a_{n}λ_{n}x_{n}
=λ_{1} (a_{1}x_{1} + \frac{λ_{2}}{λ_{1}} a_{2}x_{2} + ⋯ + \frac{λ_{n}}{λ_{1}} a_{n}x_{n})

Iteriamo questa operazione, considerando v_{2} = Av_{1} = A^{2}v_{0},..., v_{k} = A^{k}v_{0}:

v_{k} =a_{1}A^{k}x_{1} + a_{2}A^{k}x_{2} + ⋯ + a_{n}A^{k}x_{n} = a_{1}λ_{1}^{k}x_{1} + a_{2}λ_{2}^{k}x_{2} + ⋯ + a_{n}λ_{n}^{k}x_{n}
(*)  =λ_{1}^{k} (a_{1}x_{1} + (\frac{λ_{2}}{λ_{1}} )^{k} a_{2}x_{2} + ⋯ + (\frac{λ_{n}}{λ_{1}} )^{k} a_{n}x_{n})

Dato che |\frac{λ_{i}}{λ_{1}} | < 1 per i = 1, ..., n, abbiamo che

(**)  lim_{k→∞} \frac{1}{λ_{1}^{k}} v_{k} = a_{1}x_{1}

Ora ragioniamo sulle singole componenti di (v_{k})_{i} e (x_{1})_{i} rispetto alla base canonica; il limite appena scritto consente di calcolare, per i ∈ {1, ..., n} tale che (x_{1})_{i} ≠ 0:

lim_{k→∞} \frac{(v_{k+1})_{i}}{(v_{k})_{i}} = lim_{k→∞} \frac{(v_{k+1})_{i}}{(v_{k})_{i}}

= lim_{k→∞} \frac{λ_{1}^{k+1}}{(v_{k})_{i}} (a_{1}x_{1} + (\frac{λ_{2}}{λ_{1}} )^{k+1} a_{2}x_{2} + ⋯ + (\frac{λ_{n}}{λ_{1}} )^{k+1} a_{n}x_{n})_{i}

= lim_{k→∞} \frac{λ_{1}}{(a_{1}x_{1})_{i}} (a_{1}x_{1} + (\frac{λ_{2}}{λ_{1}} )^{k+1} a_{2}x_{2} + ⋯ + (\frac{λ_{n}}{λ_{1}} )^{k+1} a_{n}x_{n})_{i}

= \frac{λ_{1}}{(a_{1}x_{1})_{i}} (a_{1}x_{1})_{i}

= λ_{1}

(prima abbiamo usato la (*) e poi la (**)).

Dunque la successione (v_{k+1})_{i}/(v_{k})_{i} converge a λ_{1} e la velocità di questa convergenza è misurata dalle potenze |λ_{i}/λ_{1}|^{k}; naturalmente possiamo scegliere di normalizzare i vettori v_{m} rispetto a una certa norma, dato che ci interessa solo il rapporto delle loro compomenti; inoltre moltiplicare per un fattore non nullo non altera l’essere un vettore x_{k} un autovettore di autovalore λ_{k}.

Questo procedimento suggerisce un algoritmo iterativo per calcolare λ_{1} partendo da un valore iniziale arbitrario v_{0}.

/* Calcola un autovalore lambda1 di una matrice A (nelle ipotesi di crescenza degli autovalori) usando una successione di vettori normalizzati

*/

float lambda1( float A[n][n], /* matrice */ int kmax /* numero di iterazioni */ )

� float v[n]; /* vettore usato nel calcolo */ float w[n]; /* vettore usato nel calcolo */ float a; /* autovalore */ float f; /* variabile ausiliaria */ int i, j, k; /* indici di ciclo */ /* inizializza il vettore ad uno di norma (euclidea) 1 */ v[0] = 1; for(i=1; i�n; i++) v[i] = 0;

for(k=0; k�kmax; k++)� /* ora pone w = Av */ for(i=0; i�n; i++)� w[i] = 0; for(j=0; j�n; j++)

w[i] += A[i][j]*v[j]; �; /* sceglie un indice per cui v[i]!=0 */ for(i=0; v[i]==0 && i�n; i++); /* calcola l’autovalore */ a = w[i]/v[i]; /* ora assegna a v il normalizzato di w */ f = norma(w); for(i=0; i�n; i++) v[i] = w[i]/f;

�;

return a; �

Osserviamo che questa funzione calcola anche una approssimazione di un autovettore (di norma 1) relativo all’autovalore cercato; comunque la convergenza di questo algoritmo non segue esattamente dal ragionamento svolto in precedenza, perché la scelta dell’indice rispetto al quale considerare la componente del vettore v non è necessariamente la stessa ad ogni iterazione; in ogni caso la convergenza si può dimostrare anche in questo caso più generale. Inoltre questo metodo non funziona per autovalori complessi, perché abbiamo usato in modo essenziale la realtà dell’autovettore x_{1}.

Se la matrice A è simmetrica si può considerare un metodo alternativo che consiste nell’usare il quoziente di Rayleigh: ricordiamo che se v è un autovettore di una matrice, il suo autovalore è dato dal quoziente di Rayleigh (<, > denota il prodotto scalare canonico)

λ = \frac{< Av, v >}{< v, v >}

(infatti Av = λv e v ≠ 0, sicché < Av, v >= λ < v, v >).

Ora, se A è simmetrica, possiede una base ortonormale di autovettori {x_{1}, ..., x_{n}} e

v_{k} = \frac{1}{c_{k}} (a_{1}λ_{1}^{k}x_{1} + ... + a_{n}λ_{n}^{k}x_{n})

ove la costante c_{k} è il prodotto ∥w_{1}∥ · ∥w_{n}∥ delle norme delle approssimazioni degli autovettori (non normalizzati) trovati al passo precedente dell’iterazione; si pone allora

v_{k+1} = \frac{w_{k+1}}{∥w_{k+1}∥}

e si trova

< w_{k+1}, v_{k} >= \frac{< Av_{k}, v_{k} >}{< v_{k}, v_{k} >} = \frac{∑_{i=1}^{n} a_{i}^{2}λ_{i}^{2k+1}}{∑_{i=1}^{n} a_{i}^{2}λ_{i}^{2k}} = λ_{1} \frac{a_{1} + ∑_{i=2}^{n} a_{i}^{2} (\frac{λ_{i}}{λ_{1}} )^{2k+1}}{a_{1} + ∑_{i=2}^{n} a_{i}^{2} (\frac{λ_{i}}{λ_{1}} )^{2k}}

Allora, passando al limite per k → ∞, troviamo

lim_{k→∞} < w_{k+1}, v_{k} >= λ_{1}

con velocità di convergenza pari a quella con la quale |λ_{2}/λ_{1}|^{2m} tende a zero. Le opportune modifiche alla funzione precedente sono:

/* Calcola un autovalore lambda1 di una matrice SIMMETRICA A (nelle ipotesi di crescenza degli autovalori) usando una successione di vettori normalizzati

*/ float asimm(

float A[n][n], /* matrice */ int kmax /* numero di iterazioni */ )

� float v[n]; /* vettore usato nel calcolo */ float w[n]; /* vettore usato nel calcolo */ float a; /* autovalore */ float f; /* variabile ausiliaria */ int i, j, k; /* indici di ciclo */ /* inizializza il vettore ad uno di norma (euclidea) 1 */ v[0] = 1; for(i=1; i�n; i++)

v[i] = 0; for(k=0; k�kmax; k++)�

/* ora pone w = Av */

for(i=0; i�n; i++)� w[i] = 0; for(j=0; j�n; j++)

w[i] += A[i][j]*v[j]; �; /* calcola l’autovalore con Reyleigh: a =< v, w >,

cio�e a = v_{1}w_{1} + ... + v_{n}w_{n} */

a = w[0]*v[0]; for(i=1; i�n; i++)

a += w[i]*v[i]; /* ora assegna a v il normalizzato di w */ f = norma(w); for(i=0; i�n; i++) v[i] = w[i]/f;

�;

return a; �

Ad esempio sulla matrice

A = ( \begin{matrix}1 & 2 & 3 & 4 \\ 2 & 3 & 4 & 0 \\ 3 & 4 & 1 & 2 \\ 4 & 0 & 2 & 3 \end{matrix} )

questo metodo, su una macchina a 32 bit, ha prodotto la soluzione 9.5208 in 9 passi:

1.80562E-035 1 7.46667 9.31746 9.50276 9.51909 9.52062 9.52078 9.52079 9.5208

mentre il metodo delle potenze ne ha richiesti 16.

Supponiamo ora che gli autovalori siano ordinati come

|λ_{1}| ≥ |λ_{2}| ≥ ... ≥ |λ_{n}|

e di aver già calcolato λ_{1} ed un suo autovettore v_{1}; se vogliamo determinare λ_{2} ed un suo autovalore possiamo ancora ricorrere al metodo delle potenze, con la seguente modifica: supponiamo di poter determinare una matrice ortogonale U tale che

Uv_{1} = −σe_{1}

ove e_{1} è il primo vettore della base canonica e |σ| = ∥v_{1}∥_{2}; allora

(UAU^{T})Uv_{1} = λ_{1}Uv_{1} ⟹ UAU^{T} e_{1} = λ_{1}e_{1}

Ma A e UAU^{T} hanno gli stessi autovalori (U è ortogonale, quindi se λ è un autovalore di A e v un suo autovettore: UAU^{T} Uv = UAv = λUv); inoltre

UAU^{T} = ( \begin{matrix}λ_{1} & a^{T} \\ 0 & A_{1} \end{matrix} )

ove A_{1} è una matrice (n − 1) × (n − 1) con gli stessi autovalori di A (eccetto al più λ) e a ∈ ℝ^{n−1}. Questo si vede così: certamente

UAU^{T} = ( \begin{matrix}λ_{1} & a^{T} \\ b & A_{1} \end{matrix} )

con b ∈ ℝ^{n−1}; ma UAU^{T} e_{1} = λ_{1}Uv_{1} da cui

( \begin{matrix}λ_{1} \\ 0 \end{matrix} ) = ( \begin{matrix}λ_{1} & a^{T} \\ b & A_{1} \end{matrix} ) ( \begin{matrix}1 \\ 0 \end{matrix} ) = ( \begin{matrix}λ_{1} \\ b \end{matrix} )

Ora: poiché i rimanenti autovalori di A sono esattamente quelli di A_{1} possiamo applicare a A_{1} il metodo delle potenze per determinare il suo autovalore di modulo massimo, cioè λ_{2}; inoltre in questo modo abbiamo un autovettore w_{1} di A_{1} relativo a quell’autovalore.

Per trovare un autovettore di A relativo allo stesso autovalore basta trovare un x ∈ ℝ tale che

( \begin{matrix}λ_{1} & a^{T} \\ 0 & A_{1} \end{matrix} ) ( \begin{matrix}x \\ w_{1} \end{matrix} ) = λ_{2} ( \begin{matrix}x \\ w_{1} \end{matrix} )

cioè λ_{1}x + a^{T} w_{1} = λ_{2}x da cui

x = \frac{a^{T} w_{1}}{λ_{2} − λ_{1}}

Allora il vettore U^{T} ( \begin{matrix}x \\ w_{1} \end{matrix} ) è autovettore di A relativo all’autovalore λ_{2}.

Iterando questa tecnica in teoria potrebbero trovarsi tutti gli autovalori di A e per ciascuno di essi un autovettore non nullo; comunque la propagazione degli errori a ciascun passo è drammatica, perché si utilizzano soluzioni sempre più approssimate.

Una ulteriore modifica al metodo delle potenze potrebbe essere sfruttare la conoscenza di un valore approssimato per λ_{1} e cercare di migliorare questa approssimazione: questo è realizzato dal metodo delle potenze inverse.

Supponiamo che l sia una approssimazione di un autovalore λ di A: allora, se v è un autovettore, (A − lI)v = (λ − l)v e quindi è corretto affermare che (λ − l) è autovalore di (A − lI), sicché (λ − l)^{−1} è autovalore di (A − lI)^{−1}; se l’approssimazione l di λ è buona allora μ = (λ − l)^{−1} è l’autovalore di modulo massimo di (A − lI)^{−1}.

Dunque si può calcolare μ (e quindi λ = l +1/μ) usando il metodo delle potenze sulla matrice (A − lI)^{−1}.

/* Metodo delle potenze inverse: data la matrice A ed una approssimazione l dell’autovalore di modulo massimo torna una migliore approssimazione di l e calcola anche una approssimazione dell’autovettore relativo a l; richiede anche il massimo numero di iterazioni da effettuare ed un valore di precisione oltre il quale l’approssimazione �e ritenuta soddisfacente; torna 0 se �e possibile calcolare l con la precisione richiesta, torna 1 se la matrice non �e invertibile, torna 2 se non riesce a calcolare l con la precisione desiderata.

*/ int pot inv(

float A[n][n], /* matrice */ float x[n], /* autovettore: usato solo in output */ float *l, /* approssimazione iniziale e finale */ float prec, /* precisione desiderata */ int kmax /* massimo numero di iterazioni */ )

� float pivot[n]; /* usato nella chiamata a factor() */ float w[size]; /* approssimazione successiva di x[] */ float lambda; /* approssimazione successiva di l */ float a; /* variabile ausiliaria */ int i, j, k, k0; /* indici di ciclo */ /* per prima cosa cambia A in A − lI */ for(i=0; i�n; i++)

A[i][i] -= *l; if(factor(A, pivot)==0)

/* se il determinante di A �e zero... */ return 1;

/* inizializza x a (1,1,...,1) */ for(i=0; i�n; i++)

x[i] = 1; lambda = *l; for(k=0; k�kmax; k++)�

/* ora calcola w = Gx (G sta nella parte inferiore

di A): cio�e w_{i} = x_{i} + ∑_{j=1}^{i−1} a_{ij}x_{j} */ for(i=0; i�n; i++)�

w[i] = x[i]; for(j=0; j�i; j++)

w[i] += A[i][j]*x[j]; �; /* ora calcola a = ∥w∥_{∞} */ a = fabs(w[0]); for(i=1; i�n; i++)

if(fabs(w[i])�a)

a = fabs(w[i]); /* e trova l’indice k0 per il quale w[k0]

assume il valore massimo */ for(k0=0; ; k0++)

if(fabs(w[k0])==a) break;

/* ora mette in x il normalizzato di w */ for(i=0; i�n; i++)

x[i] = w[i]/a; /* e calcola una migliore approssimazione di l */ a = lambda; /* salva lambda */ lambda = *l + x[k0]/w[k0]; /* ora controlla se ha raggiunto la precisione voluta */ if(fabs(a-lambda)�=prec*fabs(lambda))

28 §1.8. Trasformazioni di Householder

return 0; �; *l = lambda; return 2;

�

§1.8. Trasformazioni di Householder

Partiamo da una matrice invertibile S: una similitudine è la trasformazione lineare nello spazio delle matrici data da

A ↦ SAS^{−1}

Se A e B sono tali che esiste S invertibile tale che B = SAS^{−1} allora A e B si dicono simili: la similitudine è iuna relazione di equivalenza.

Inoltre, se Av = λv sono un autovalore ed un suo autovettore di A allora

SAS^{−1}Sv = Sλv = λSv

sicché A e SAS^{−1} hanno gli stessi autovalori ed i relativi autovettori sono i trasformati di quelli di A per tramite dell’isomorfismo S.

Vale il

8.1. Teorema. Se A è una matrice reale di ordine n allora esiste una matrice ortogonale U tale che H = UAU^{T} è una matrice di Hessenberg, tridiagonale se A è simmetrica.

Chiaramente questo è un caso particolare del teorema spettrale: il vantaggio è che U può costruirsi con un numero finito di passi.

8.2. Definizione. Un riflettore elementare è una matrice del tipo

U = I − 2uu^{T}

ove u è un vettore di norma 1.

Chiaramente un riflettore elementare determina un isomorfismo che si chiama trasformazione di Householder.

8.3. Lemma. U è simmetrica, ortogonale ed idempotente.

(Infatti U^{T} = (I − 2uu^{T})^{T})I − 2u^{T T} u^{T} = U e U^{T} U = U^{2} = (I − 2uu^{T})(I − 2uu^{T}) = I − 2uu^{T} − 2uu^{T} + 4uu^{T} uu^{T} = I − 4uu^{T} + 4uu^{T} = I.)

Geometricamente, Ux è il vettore riflesso di x rispetto all’asse perpendicolare a u e passante per l’origine: infatti

Ux = x − 2(uu^{T})x = x − 2u(u^{T} x) = x − 2 < u, x > u

dato che < u, x > u è la proiezione di x su u.

I riflettori possono essere utilizzati per azzerare opportune coordinate di vettori:

8.4. Teorema. Se x è un vettore (non nullo), il riflettore

U = I − \frac{1}{p}uu^{T}

ove u = x + σe_{1}, essendo σ = ±∥x∥ e p = ∥u∥^{2}/2, è tale che

Ux = −σe_{1}

Notiamo intanto che

p =\frac{1}{2} < u, u >= \frac{1}{2} < x + σe_{1}, x + σe_{1} >= \frac{1}{2}(∥x∥^{2} + 2σx_{1} + σ^{2})
=\frac{1}{2}(∥x∥^{2} + 2σx_{1} + ∥x∥^{2}) = ∥x∥^{2} + σx_{1}

ove x_{1} è la prima coordinata del vettore x nella base canonica. Ora

Ux =x − \frac{1}{p}(uu^{T})x = x − \frac{1}{p}(x + σe_{1})(x + σe_{1})^{T} x
=x − \frac{1}{p}(x + σe_{1})(x^{T} + σe_{1}^{T})x
=x − \frac{1}{p}(xx^{T} + σxe_{1}^{T} + σe_{1}x^{T} + σ^{2}e_{1}e_{1}^{T})x
=x − \frac{1}{p}(xx^{T} x + σxe_{1}^{T} x + σe_{1}x^{T} x + σ^{2}e_{1}e_{1}^{T} x)
=x − \frac{1}{p}(∥x∥^{2}x + σx_{1}x + σ∥x∥^{2}e_{1} + σ^{2}x_{1}e_{1})
=x − \frac{∥x∥^{2} + σx_{1}}{p} x − \frac{σ∥x∥^{2} + σ^{2}x_{1}}{p} e_{1}
= − σe_{1}

(nell’ultimo passaggio abbiamo sfruttato la precedente p = ∥x∥^{2} + σx_{1}).

Per determinare completamente un riflettore elementare di questo tipo, basta quindi il numero p e il vettore u: costruiamo direttamente questi a partire da x:

�include �math.h�

/* Costruzione di un riflettore elementare in termini di un vettore non nullo x: restituisce in x stesso il vettore di riflessione e in come valore il termine p = \frac{1}{2}∥u∥^{2}

*/

float rifl elem( float x[n], /* vettore da riflettere */ )

30 §1.9. Fattorizzazione QR

� float sigma; /* coefficiente del riflettore */ float a; /* contiene la norma di x */ int i; /* indice di ciclo */ /* calcola a = ∥x∥_{∞} */ a = fabs(x[0]); for(i=1; i�n; i++)

if(fabs(x[i])�a) a = x[i];

/* ora inizia il ciclo principale */

sigma = 0;

for(i=0; i�n; i++)

/* se |x_{i}| > a\sqrt{ε} allora pone σ = σ + (x_{i}/a)^{2}

la condizione serve a trascurare quantit�a inutili */ if(fabs(x[i])�a*sqrt(epsilon))

/* normalizza per evitare over/under-flow */ sigma += pow(x[i]/a,2);

/* Il segno di σ �e arbitrario, quindi possiamo supporlo > 0 o < 0 a piacimento ed usare il segno che non d�a luogo a cancellazione numerica:

σ = sgn x_{1}\sqrt{σ}a */

sigma = sqrt(sigma)*a;

if(x[0]�0)

return sigma*(x[0]-=sigma); return sigma*(x[0]+=sigma);

�

Conoscendo p ed il vettore u possiamo calcolare z = Uy come segue: per prima cosa determiniamo

a = \frac{1}{p} ∑_{i=1}^{n} u_{i}y_{i}

quindi si calcola il risultato come

z = y − au

§1.9. Fattorizzazione QR

Generalizziamo la costruzione precedente dei riflettori elementari al seguente caso: abbiamo visto come, dato un vettore x ≠ 0, produrre un riflettore elementare che lo trasformi in un vettore con la sola prima coordinata non nulla (U = I − \frac{1}{p}uu^{T}); chiediamoci se non sia possibile trovare un riflettore che trasformi x in un vettore con le sole prime k coordinate non nulle.

L’idea è di considerare U^{(n−k+1)} = I− \frac{1}{p}uu^{T} (sempre con u = x±∥x∥e_{1} e p = ∥u∥^{2}/2) nello spazio ℝ^{n−k+1} e di usarlo per formare un operatore nello spazio ℝ^{n} definito come

U^{(n)} = I_{k−1} ⊕ U^{(n−k+1)} = I − \frac{1}{p_{k}} u_{k}u_{k}^{T}

ove

u_{k} = 0_{k−1} ⊕ u_{n−k+1}

e

p_{k} = p_{n−k+1} = \frac{1}{2}∥u_{n−k+1}∥^{2}

U^{(n−k+1)} è il riflettore elementare in ℝ^{n−k+1} tale che

U^{(n−k+1)}(x_{k}e_{k} + ... + x_{n}e_{n}) = −σ_{k}e_{k}

Avendo a disposizione per ciascun k = 1, ..., n riflettori U_{1}, ..., U_{n} di questo tipo potremmo, data una matrice A, ridurla in forma triangolare superiore per mezzo delle moltiplicazioni

U_{n−1}...U_{2}U_{1}A = R

in modo che, se Q = (U_{n−1}...U_{2}U_{1})^{T} = U_{1}^{T} U_{2}^{T} ...U_{n−1}^{T} allora

Q^{T} A = R ⟹ A = QR

Si noti che QQ è ortogonale essendo prodotto di matrici ortogonali.

Ciascun fattore U_{i} è determinato da un numero p_{i} e da un vettore u_{i} ∈ ℝ^{n−i+1}.

Applichiamo ora questa costruzione alla riduzione di una matrice simmetrica A in una matrice simile (quindi con gli stessi autovalori) che sia tridiagonale: l’idea è di azzerare successivamente gli elementi sotto la diagonale usando i riflettori U_{i}.

Si comincia con la prima prima colonna: consideriamo il vettore (a_{11}, a_{21}, ..., a_{n1})^{T} ed il riflettore U_{2} che azzeri gli elementi a_{31}, ..., a_{n1}: la matrice (i riflettori sono idempotenti ed ortogonali: U = U^{T} = U^{−1})

A_{2} = U_{2}AU_{2}

è simile ad A; ora scriviamo A a blocchi come

( \begin{matrix}a_{11} & v_{1}^{T} \\ v_{1} & A_{11} \end{matrix} )

e calcoliamo

U_{2}AU_{2} = ( \begin{matrix}1 & 0 \\ 0 & U^{(n−1)} \end{matrix} ) ( \begin{matrix}a_{11} & v_{1}^{T} \\ v_{1} & A_{11} \end{matrix} ) ( \begin{matrix}1 & 0 \\ 0 & U^{(n−1)} \end{matrix} )
= ( \begin{matrix}a_{11} & v_{1}^{T} \\ U^{(n−1)}v_{1} & U^{(n−1)}A_{11} \end{matrix} ) ( \begin{matrix}1 & 0 \\ 0 & U^{(n−1)} \end{matrix} )
= ( \begin{matrix}a_{11} & v_{1}^{T} U^{(n−1)} \\ U^{(n−1)}v_{1} & U^{(n−1)}A_{11}U^{(n−1)} \end{matrix} )

Ora osserviamo che U^{(n−1)} applicato ad un vettore di ℝ^{n−1} ne azzera tutte le componenti tranne la prima, sicché U^{(n−1)}v_{1} = (∗, 0, ..., 0)^{T} (con ∗ indichiamo un valore non nullo); inoltre v_{1}^{T} U^{(n−1)} = ((U^{(n−1)})T v_{1})^{T} = (U^{(n−1)}v_{1})^{T} = (∗, 0, ..., 0). Dunque A_{2} sarà una matrice (simmetrica) della forma

A_{2} = ( \begin{matrix}∗ & ∗ & 0 & ... & 0 \\ ∗ & ∗ & ∗ & ... & ∗ \\ 0 & ∗ & ∗ & ... & ∗ \\ ⋮ & ⋮ & ⋮ & ... & ⋮ \\ 0 & ∗ & ∗ & ... & ∗ \end{matrix} )

Ora iteriamo il procedimento considerando un riflettore U_{3} che applicato alla seconda colonna di questa matrice A_{2} ne azzeri gli ultimi n − 3 elementi, in modo che U_{3}A_{2}U_{3} sia una matrice

A_{2} = ( \begin{matrix}∗ & ∗ & 0 & 0 & ... & 0 \\ ∗ & ∗ & ∗ & 0 & ... & 0 \\ 0 & ∗ & ∗ & ∗ & ... & ∗ \\ 0 & 0 & ∗ & ∗ & ... & ∗ \\ ⋮ & ⋮ & ⋮ & ⋮ & ... & ⋮ \\ 0 & 0 & ∗ & ∗ & ... & ∗ \end{matrix} )

Proseguiamo in questo modo producendoa alla fine del processo una matrice simmetrica

A_{n−1} = U_{n−1}...U_{3}U_{2}AU_{2}U_{3}...U_{n−1}

tridiagonale, simile ad A.

Osserviamo inoltre che, se A non è simmetrica questo procedimento fornisce una matrice di Hessenberg simile ad A.

Diamo ora un algoritmo per questo procedimento: come al solito non sarà necessario costruire esplicitamente le matrici U_{i}, ma solo determinarle per mezzo delle coppie (p_{i}, u_{i}).

/* Data una matrice simmetrica A costruisce una tridiagonale a lei simile usando i riflettori elementari; poich�e il risultato �e determinato semplicemente dalla diagonale e dalla diagonale immediatamente superiore (quella inferiore �e uguale a questa dato che A �e simmetrica), il risultato consister�a in questa coppia di diagonali: la diagonale principale nel vettore d1 e la diagonale secondaria nel vettore d2; inoltre le coppie (p,u) che determinano i riflettori elementari sono poste nella parte triangolare strettamente inferiore di A come colonne il cui elemento sulla diagonale �e p (o zero se la trasformazione in quel passo �e l’identit�a) con sotto il vettore colonna u; la parte triangolare superiore di A �e lasciata intatta.

*/ void qr factor(

float A[n][n], /* matrice da trasformare */ float d1[n], /* conterr�a la diagonale principale */ float d2[n] /* conterr�a la codiagonale */ )

� float a; /* variabile ausiliaria */ float sigma; /* coefficiente del riflettore */ int i, j, k; /* indici di ciclo */ for(k=0; k�n-2; k++)�

/* ora decompone a blocchi la matrice */ d1[k] = A[k][k]; /*

ora determina il riflettore U_{k+1} */ /* calcola a = max_{i=k+1}^{n} |a_{ik}| */ a = fabs(A[k+1][k]); for(i=k+2; i�n; i++)

if(fabs(A[i][k])�a) a = A[i][k];

/* se la colonna �e di tutti zeri U_{k+1} = I */ if(a==0)�

A[k][k] = 0; /* p_{k+1} = 0 */

d2[k] = 0;

continue; /* ripete il ciclo for(k...) */ � /* ora calcola p_{k+1} */ sigma = 0; for(i=k+1; i�n; i++)

/* se |a_{ik}| > a\sqrt{ε} allora pone σ = σ + (a_{ik}/a)^{2}

serve a trascurare quantit�a inutili */ if(fabs(A[i][k])�a*sqrt(epsilon))

/* normalizza per evitare over/under-flow */ sigma += pow(A[i][k]/a,2);

/* Il segno di σ �e arbitrario, quindi possiamo supporlo > 0 o < 0 a piacimento ed usare il segno che non d�a luogo a cancellazione numerica:

σ = sgn a_{k+1,k}\sqrt{σ}a */

sigma = sqrt(sigma)*a; if(A[k+1][k]�0)

A[k+1][k] -=sigma; else

A[k+1][k] +=sigma; A[k][k] = sigma*A[k+1][k]; /* p_{k+1} */ d2[k] = -sigma; /* U_{k+1}v_{k} = −σe_{k} */ /*

ora calcola la matrice B = U_{k+1}A_{k}

*/ for(j=k+1; j�n; j++)�

/* a = \frac{1}{a_{kk}} ∑_{i=k+1}^{n} a_{ik}a_{ij} */

a = A[k+1][k]*A[k+1][j]; for(i=k+2; i�size; i++)

a += A[i][k]*A[i][j]; a /= A[k][k];

/* a_{ij} = a_{ij} − aa_{jk} */

for(i=k+1; i�n; i++)

A[i][j] -= a*A[i][k]; �; /*

ora calcola la matrice A_{k+1} = BU_{k+1}

*/ for(i=k+1; i�n; i++)�

/* a = \frac{1}{a_{kk}} ∑_{j=k+1}^{n} a_{ij}a_{jk} */

a = A[i][k+1]*A[k+1][k]; for(j=k+2; j�size; j++)

a += A[i][j]*A[j][k]; a /= A[k][k];

/* a_{ij} = a_{ij} − aa_{jk} */

for(j=k+1; j�n; j++)

A[i][j] -= a*A[j][k]; �;

�; /* ciclo for (k...) */ d1[n-2] = A[n-2][n-2]; d1[n-1] = A[n-1][n-1]; d2[n-2] = A[n][n-1];

�

Una applicazione di questo metodo è alla determinazione degli autovalori di una matrice simmetrica: per prima cosa la si riduce in forma tridiagonale

B = ( \begin{matrix}a_{1} & b_{1} & 0 & ... & 0 & 0 \\ b_{1} & a_{2} & b_{2} & ... & 0 & 0 \\ 0 & b_{2} & a_{3} & ... & 0 & 0 \\ ⋮ & ⋮ & ⋮ & ... & ⋮ & ⋮ \\ 0 & 0 & 0 & ... & a_{n−1} & b_{n−1} \\ 0 & 0 & 0 & ... & b_{n−1} & a_{n} \end{matrix} )

Evidentemente possiamo supporre che ciascun b_{i} sia diverso dallo zero, altrimenti la matrice si decompone in somma diretta di matrici tridiagonali con questa proprietà, e gli autovalori della matrice sono dati da quelli delle singole matrici addendi diretti.

Ora scriviamo λI − B il cui determinante è il polinomio caratteristico:

λI − B = ( \begin{matrix}λ − a_{1} & −b_{1} & 0 & ... & 0 & 0 \\ −b_{1} & λ − a_{2} & −b_{2} & ... & 0 & 0 \\ 0 & −b_{2} & λ − a_{3} & ... & 0 & 0 \\ ⋮ & ⋮ & ⋮ & ... & ⋮ & ⋮ \\ 0 & 0 & 0 & ... & λ − a_{n−1} & −b_{n−1} \\ 0 & 0 & 0 & ... & −b_{n−1} & λ − a_{n} \end{matrix} )

e consideriamo la successione di polinomi

\{ \begin{matrix}p_{0}(λ) = 1 \\ p_{1}(λ) = λ − a_{1} \\ ... \\ p_{i}(λ) = (λ − a_{1})p_{i−1}(λ) − b_{i−1}^{2}p_{i−2}(λ) \\ ... \end{matrix}

Si tratta dei determinanti delle sottomatrici principali di ordine i estratte dalla matrice λI − B; in particolare p_{n}(λ) è il polinomio caratteristico di B.

Si possono dimostrare le seguenti proprietà dei polinomi p_{i}:

1) Le radici di p_{i} sono reali e distinte (se i > 1);

2) Per ogni x ∈ ℝ tale che se p_{n}(x) ≠ 0, il numero di radici (cioè di autovalori) maggiori di x è dato dal numero di variazioni del segno nella successione di numeri reali

1, p_{1}(x), p_{2}(x), ..., 1, p_{n}(x)

dalla quale siano stati eliminati eventuali valori nulli.

Per localizzare le radici del polinomio caratteristico, cioè gli autovalori della matrice è anche utile il teorema di Gerschgorin:

9.1. Teorema. Se A è una matrice n × n, poniamo

r_{i} = ∑_{{j=1 \atop j≠i}}^{n} |a_{ij}|  e  c_{j} = ∑_{{i=1 \atop i≠j}}^{n} |a_{ij}|

Allora 1) Ogni autovalore di A appartiene all’insieme

R = ⋃_{i=1}^{n} R_{i}  ove  R_{i} = {z | |z − a_{ii}| < r_{i}}

2) Ogni autovalore di A appartiene all’insieme

C = ⋃_{i=j}^{n} C_{j}  ove  C_{j} = {z | |z − a_{jj}| < c_{j}}

3) Ogni componente connessa di R o C contiene tanti autovalori contati con le loro molteplicità quanti sono i cerchi della componente stessa.

Ora: supponiamo di conoscere un intervallo [a, b] che contenga tutti gli autovalori di B (ad esempio considerando il raggio spettrale della matrice o utilizzando il teorema di Gerschgorin); allora possiamo, considerato c = (b − a)/2, dire quanti autovalori si trovano nell’intervallo [c, b) e quanti nell’intervallo (a, c) usando la 2) con c = x; a loro volta suddividiamo questi intervalli a metà, considerando cioè c_{1} = (c − a)/2 e c_{2} = (b − c)/2 e di nuovo determiniamo il numero di autovalori che cadono in questi intervalli usando la 2).

Iterando queste bisezioni possiamo (dato che gli autovalori sono in numero finito, non pi di n) con un numero finito di bisezioni determinare una suddivisione di (a, b) ion intervalli finiti di eguale ampiezza tali che ciascuno di essi non contenga più di un autovalore; allora, se λ è contenuto nell’intervallo (c_{i}, c_{j}) possiamo considerare come valore iniziale per l’approssimazione di λ il punto medio (c_{j} −c_{i})/2; a questo punto, magari col metodo delle potenze inverse, possiamo provare a determinare migliori approssimazioni degli autovalori stessi.

36 §1.10. Approssimazione di funzioni

§1.10. Approssimazione di funzioni

Il problema dell’approssimazione consiste nel sostituire una funzione f con un’altra funzione che abbia un comportamento simile ma che consenta in pratica di effettuare i calcoli richiesti: ad esempio si potrebbe voler calcolare un integrale, o simili.

Le classi di funzioni alle quali “si attinge” per trovare candidati validi per queste sostituzioni sono:

1) Polinomi di grado n: ℙ_{n}

∑_{i=0}^{n} a_{n}x^{n}

ad esempio una funzione continua in un compatto si approssima sempre, con precisione desiderata ε > 0, con un poilinomio P_{ε};

2) Polinomi trigonometrici di grado n e periodo ω: 𝕋_{n}(ω)

a_{0} + ∑_{i=1}^{n}(a_{i} cos iωx + b_{i} sin iωx)

ad esempio una funzione continua e periodica si approssima in un compatto con polinomi trigonometrici;

3) Funzioni razionali quozienti di polinomi di gradi n e d: ℝ_{n,d}

\frac{P_{n}(x)}{Q_{d}(x)}

ad esempio funzioni che hanno poli o singolarità nell’intorno di un certo punto si lasciano approssimare in ℝ_{n,d}.

4) Somme esponenziali di ordine n: 𝔼_{n}

∑_{i=1}^{n} a_{i} exp (−b_{i}x)

ad esempio funzioni esponenziali possono approssimarsi in 𝔼_{n}.

5) Funzioni spline di ordine n: S_{n}; si tratta di funzioni C^{n} che a tratti sono polinomi di grado n.

Infine si possono moltiplicare fra loro elementi appartenenti a diverse fra le classi precedenti.

Supponiamo di voler approssimare una funzione f(x) della quale conosciamo il valore in un certo insieme finito {x_{0}, ..., x_{m}} di punti; lo faremo scegliendo l’approssimante f_{n}(x) in una classe fra le precedenti.

Questa scelta può avvenire ad esempio in base ad uno dei seguenti criteri:

a) Interpolazione: si sceglie f_{n} in modo che per ogni i = 0, ..., m si abbia f_{n}(x_{i}) = f(x_{i}); ovviamente la conoscenza degli f(x_{i}) deve essere sufficientemente precisa.

b) Minimi quadrati: si sceglie una f_{n} nella classe scelta “disti il meno possibile” da f nei punti nei quali questa è nota: cioè si sceglie la f_{n} che rende minimo il numero

∑_{i=0}^{m}(f_{n}(x_{i}) − f(x_{i}))^{2}

Talvolta è opportuno “pesare” la distanza come

∑_{i=0}^{m}(f_{n}(x_{i}) − f(x_{i}))^{2}p_{i}

ove p_{i} > 0 sono appunti dei pesi che sono tanto più grandi quanto più stimiamo essere accurata la conoscenza del valore della funzione f(x_{i}) la cui distanza stiamo pesando on p_{i}.

c) MiniMax: analogo al precedente ma si usa una distanza uniforme e non euclidea; cioè f_{n} è scelta in modo da minimizzare il valore

max_{0≤i≤m} |f_{n}(x_{i}) − f(x_{i})|

Poiché le funzioni f(x) che andremo ad approssimare appartengono sempre a qualche spazio vettoriale di funzioni, cercheremo di approssimarle con funzioni f_{n} che stanno in qualche sotto-spazio di dimensione finita: in questo modo, fissata una base di questo sotto-spazio, possiamo limitare il numero di funzioni approssimanti in considerazione; la scelta di questa base influenzerà in modo essenziale le prestazioni numeriche dell’approssimazione.

§1.11. Interpolazione di Lagrange

Proviamo ad approssimare funzioni definite in un intervallo [a, b] con polinomi usando il criterio di interpolazione.

Consideriamo dunque gli n + 1 punti {x_{0}, ..., x_{n}} nei quali delle approssimazioni y_{i} = f(x_{i}) del valore della funzione sono note (supporremo ovviamente distinti tali punti); consideriamo il polinomio di grado n

p_{n}(x) = a_{0} + a_{1}x + ... + a_{n}x^{n}

ed imponiamo i vincoli p_{n}(x_{i}) = y_{i}:

\{ \begin{matrix}a_{0} + a_{1}x_{0} + ... + a_{n}x_{0}^{n} = y_{0} \\ ............................ \\ a_{0} + a_{1}x_{n} + ... + a_{n}x_{n}^{n} = y_{n} \end{matrix}

Poiché vogliamo determinare p_{n}, le nostre incognite sono i numeri a_{i}, e questi soddisfano un sistema lineare la cui matrice dei coefficienti è quella di Vandermonde V_{n} ed il cui vettore dei termini noti è il vettore degli y_{i}:

( \begin{matrix}1 & x_{0} & x_{0}^{2} & ... & x_{0}^{n} \\ 1 & x_{1} & x_{1}^{2} & ... & x_{1}^{n} \\ 1 & x_{2} & x_{2}^{2} & ... & x_{2}^{n} \\ ⋮ & ⋮ & ⋮ & ... & ⋮ \\ 1 & x_{n} & x_{n}^{2} & ... & x_{n}^{n} \end{matrix} ) ( a_{0} ; a_{1} ; a_{2} ; ⋮ ; a_{n} ) = ( y_{0} ; y_{1} ; y_{2} ; ⋮ ; y_{n} )

38 §1.11. Interpolazione di Lagrange

Come noto il determinante di Vandermonde è

det V_{n} = ∏_{i>j}(x_{i} − x_{j})

ed è diverso da zero nell’ipotesi i ≠ j ⟹ x_{i} ≠ x_{j}: dunque il sistema lineare possiede un’unica soluzione, cioè esiste un unico polinomio p_{n} di grado n che soddisfa ai vincoli p_{n}(x_{i}) = y_{i} con i = 0, 1, ..., n.

Purtroppo la soluzione di un sistema lineare la cui matrice dei coefficienti è una di Vandermonde è mal condizionato, nel senso che a piccole variazioni dei dati corrispondono variazioni della soluzione approssimata non necessariamente dello stesso ordine di grandezza.

Dunque questo ragionamento vale solo teoricamente per accertare l’esistenza ed unicità del polinomio p_{n}, non per il suo calcolo effettivo.

Il metodo di Lagrange consiste in questo: fissiamo j ∈ {0, 1, ..., n} e consideriamo la funzione che nei punti x_{i} vale δ_{ij}; sia l_{j} il suo polinomio di interpolazione: siccome le x_{i} con i ≠ j (che sono n numeri distinti) debbono essere radici di l_{j} (che ha grado n), per il teorema fondamentale dell’algebra queste sono tutte e sole le radici di l_{j}(x), e dunque abbiamo che

l_{j}(x) = a_{n} ∏_{{i=0 \atop i≠j}}^{n}(x − x_{i})

ove a_{n} è il coefficiente del termine di grado n in l_{j}: questo si determina facilmente, dato che l_{j}(x_{j}) = 1 sicché

1 = a_{n} ∏_{{i=0 \atop i≠j}}^{n}(x_{j} − x_{i}) ⟹ a_{n} = \frac{1}{∏_{{i=0 \atop i≠j}}^{n}(x_{j} − x_{i})}

Dunque

l_{j}(x) = \frac{∏_{{i=0 \atop i≠j}}^{n}(x − x_{i})}{∏_{{i=0 \atop i≠j}}^{n}(x_{j} − x_{i})} = \frac{ω_{n+1}(x)}{(x − x_{j})ω′_{n+1}(x_{j})}

ove

ω_{n+1}(x) = ∏_{i=0}^{n}(x − x_{i})  e  ω′_{n+1}(x_{j}) = \frac{d}{dx}|_{x=x_{j}}ω_{n+1}(x) = ∏_{{i=0 \atop i≠j}}^{n}(x_{j} − x_{i})

A questo punto il polinomio di interpolazione di f(x) è combinazione lineare degli l_{j}(x) come

P_{n}(x) = ∑_{j=0}^{n} y_{j}l_{j}(x)

Infatti

P_{n}(x_{i}) = ∑_{j=0}^{n} y_{j}l_{j}(x_{i}) = y_{i}

I polinomi l_{j}(x) si dicono polinomi fondamentali di Lagrange.

Stimiamo ora l’errore commesso nel sostituire a f(x) il suo polinomio di interpolazione nel caso in cui f(x) presenti un comportamento regolare, cioè appartenga a qualche spazio C^{k}(a, b) con k > 0.

Precisamente supponiamo che f ∈ C^{n+1}(a, b) e stimiamo l’errore

E_{n}(x) = f(x) − P_{n}(x)

Poiché E(x_{i}) = 0 per i = 0, 1, ..., n, è naturale scrivere

E_{n}(x) = ∏_{i=0}^{n}(x − x_{i})R_{n}(x) = ω_{n+1}(x)R_{n}(x)

per qualche funzione R_{n}(x) che deve essere in C^{n+1}(a, b); dunque

f(x) = P_{n}(x) + ω_{n+1}(x)R_{n}(x)

Introduciamo ancora una funzione

G(t) = f(t) − P_{n}(t) − ω_{n+1}(t)R_{n}(x)

(con x ≠ x_{i}) che ovviamente si annulla in x, x_{0}, x_{1}, ..., x_{n}.

Poiché si tratta di una funzione derivabile n + 1 volte possiamo applicare il teorema di Rolle in ciascuno degli intervalli di estremi x, x_{0}, x_{1}, ..., x_{n}, ottenendo dei punti ξ_{1}^{(1)}, ξ_{2}^{(1)}, ...ξ_{n+1}^{(1)} nei quali la derivata di G si annulla; ripetiamo lo stesso ragionamento sulle derivate successive di G fino all’ordine n+1, troviamo senz’altro un punto ξ ∈ (a, b) tale che G^{(n+1)}(ξ) = 0, sicché

G^{(n+1)}(ξ) = f^{(n+1)}(ξ) − (n + 1)!R_{n}(x) = 0

Questo ξ dipenderà da x.

Comunque abbiamo determinato R_{n}(x) e quindi in qualche misura f:

R_{n}(x) = \frac{1}{(n + 1)!}f^{(n+1)}(ξ(x)) ⟹ f(x) = ∑_{j=0}^{n} f(x_{j})l_{j}(x) + \frac{ω_{n+1}(x)}{(n + 1)!}f^{(n+1)}(ξ)

Ad esempio, per f ≡ 1, troviamo

∑_{j=0}^{n} l_{j}(x) = 1

Il condizionamento dell’interpolazione di Lagrange può studiarsi come segue: supponiamo di perturbare i valori f(x_{0}), ..., f(x_{n}) come

\overline{y_{i}} = f(x_{i}) + ε_{i}

Allora possiamo considerare il polinomio

\overline{P_{n}}(x) = ∑_{j=0}^{n} \overline{y_{j}}l_{j}(x)

40 §1.12. Formule di Newton

e maggiorare l’errore

|P_{n}(x) − \overline{P_{n}}(x)| = |∑_{j=0}^{n} l_{j}(x)(f(x_{j}) − \overline{y_{j}})| ≤ (∑_{j=0}^{n} |l_{j}(x)|) max_{0≤j≤n} |f(x_{j}) − \overline{y_{j}}|

= (∑_{j=0}^{n} |l_{j}(x)|) max_{0≤j≤n} |ε_{j}| ≤ ||∑_{j=0}^{n} |l_{j}(x)|||_{∞} ∥ε∥_{∞}

ove ε = (ε_{0}, ..., ε_{n})^{T}); da questo, dato che ∥P_{n}∥_{∞} ≥ max_{j} |P_{n}(x_{j})| = max_{j} |f(x_{j})|, troviamo, per y = (f(x_{0}), ..., f(x_{n}))^{T}:

\frac{∥P_{n}(x) − \overline{P_{n}}(x)∥_{∞}}{∥P_{n}(x)∥_{∞}} ≤ Λ_{n} \frac{∥ε_{j}∥_{∞}}{∥f_{∞}∥}

avendo definito il numero di Lebesgue come

Λ_{n} = ||∑_{j=0}^{n} |l_{j}(x)|||_{∞}

11.1. Teorema. Se f ∈ C[a, b] e P_{n} è il polinomio di interpolazione di Lagrange relativo ai valori x_{0}, ..., x_{n} ∈ [a, b] allora

∥f − P_{n}∥_{∞} ≤ (1 + Λ_{n})E_{n}(f)

ove

E_{n}(f) = min_{Q∈ℙ_{n}} ∥f − Q∥_{∞}

(ove ℙ_{n} è lo spazio dei polinomi di grado non superiore a n).

L’idea è di considerare un polinomio Q ∈ ℙ_{n} che renda minima la ∥f − Q∥_{∞}: allora q(x) = ∑_{j} l_{j}(x)q(x_{j}) sicché

f(x) − P_{n}(x) =(f(x) − Q_{n}(x)) − (P_{n}(x) − Q_{n}(x))
=(f(x) − Q_{n}(x)) − ∑_{j=0}^{n} l_{j}(x)(f(x_{j}) − Q(x_{j}))

Fissato n, scelti x_{0}, ..., x_{n} in modo da minimizzare Λ_{n}, denotando con \overline{Λ_{n}} questo minimo si può dimostrare che

lim_{n→∞} \frac{\overline{Λ_{n}}}{log n} = \frac{2}{π}

§1.12. Formule di Newton

Consideriamo x_{0}, ..., x_{n} assegnati e distinti; la differenza divisa di ordine 1 di f è

f[x_{0}, x_{1}] = f[x_{1}, x_{0}] = \frac{f(x_{1}) − f(x_{0})}{x_{1} − x_{0}}

La differenza divisa di ordine 2 è

f[x_{0}, x_{1}, x_{2}] = \frac{f[x_{1}, x_{2}] − f[x_{0}, x_{1}]}{x_{2} − x_{0}}

ed in generale, quella di ordine n è

f[x_{0}, ..., x_{n}] = \frac{f[x_{1}, ...., x_{n}] − f[x_{0}, ..., x_{n−1}]}{x_{n} − x_{0}}

Per induzione si ha che

f[x_{0}, ..., x_{n}] = ∑_{i=0}^{n} \frac{f(x_{i})}{∏_{{j=0 \atop j≠i}}^{n}(x_{i} − x_{j})}

da cui si vede immediatamente che f[x_{0}, ..., x_{n}] è una funzione simmetrica (nel caso 1 lo è per definizione).

Vale il

12.1. Teorema. Se f ∈ C^{k}[a, b] e x0, ..., x_{n} ∈ [a, b] (non necessariamente distinti) allora esiste ξ ∈ (a, b) con min x_{i} ≤ ξ ≤ max x_{i} tale che

f[x_{0}, ..., x_{k}] = \frac{f^{(k)}(ξ)}{k!}

In particolare f[x, x, ..., x] = \frac{f^{(k)}(x)}{k!}.

Ora consideriamo

f[x, x_{0}] = \frac{f(x_{0}) − f(x)}{x_{0} − x} ⟹ f(x) = f(x_{0}) + (x − x_{0})f[x, x_{0}]

Questo definisce un polinomio interpolante di grado zero P_{0}(x) = f(x_{0}) e (x−x_{0})f[x, x_{0}] è l’errore f(x) − P_{0}(x); se aggiungiamo un valore x_{1} possiamo considerare

f[x, x_{0}, x_{1}] = \frac{f[x_{0}, x_{1}] − f[x, x_{0}]}{x1 − x} ⟹ f[x, x_{0}] = f[x_{0}, x_{1}] + (x − x_{1})f[x, x_{0}, x_{1}]

cioè

f(x) = f(x_{0}) + (x − x_{0})f[x_{0}, x_{1}] + (x − x_{0})(x − x_{1})f[x, x_{0}, x_{1}]

42 §1.13. Funzioni polinomiali a tratti

Abbiamo dunque un nuovo polinomio interpolante P_{1}(x) = f(x_{0}) + (x − x_{0})f[x_{0}, x_{1}]; proseguendo in questo modo troviamo una successione di polinomi di interpolazione per quanti sono i valori x_{i} a disposizione:

f(x) = P_{n}(x) + E_{n}(x)

ove

P_{n}(x) =f(x_{0}) + (x − x_{0})f[x_{0}, x_{1}] + (x − x_{0})(x − x_{1})f[x_{0}, x_{1}, x_{2}] + ...
+ (x − x_{0})...(x − x_{n−1})f[x_{0}, ..., x_{n}]

e

E_{n}(x) = (x − x_{0})...(x − x_{n})f[x, x_{0}, ..., x_{n}]

Questa è la formula di interpolazione di Newton.

Scriviamo un algoritmo che dati x_{0}, ..., x_{n} e f(x_{0}), ..., f(x_{n}) costruisca il valore del polinomio di interpolazione di Newton in un punto p.

�define n /* numero di punti x_{i} */

/* Dati i vettori x e y contenenti i punti ed i valori della funzione in questi punti, e dato un punto t restituisce il valore del polinomio di Newton in t

*/ float newton(

float x[n+1], /* x_{0}, ..., x_{n} */ float y[n+1], /* f(x_{0}), ..., f(x_{n}) */ float t /* valore nel quale calcolare il polinomio */ )

� float p; /* conterr�a il risultato */ int i, j; /* indici di ciclo */ /*

Calcola le colonne della tabella delle differenze divise e lascia in y[] la diagonale di questa tabella

*/ for(i=1; i�=n; i++)

for(j=n; j�i; j--)

y[j] = (y[j]-y[j-1])/(x[j]-x[j-i]); /*

Valuta il polinomio di Newton in t */ p = y[n]; for(i=n-1; i�=0; i--)

p = y[i]+(t-x[i])*p; return p;

�

§1.13. Funzioni polinomiali a tratti

Consideriamo la sequenza di punti

a = x_{0} < x_{1} < ... < x_{n−1} < x_{n} = b

Più n è grande e più la funzione polinomiale approssimante tende ad oscillare: può essere quindi conveniente considerare come approssimante non un singolo polinomio di grado n ma una sequenza di polinomi di grado d al più tre.

Ad esempio nel caso d = 1 questo si riduce ad approssimare il grafico della funzione con la spezzata che congiunge i punti del piano (x_{i}, f(x_{i})), e che ha equazione

x ∈ [x_{i}, x_{i+1}] ⟹ f_{n}(x) = \frac{(x_{i+1} − x)y_{i} + (x − x_{i})y_{i+1}}{x_{i+1} − x_{i}}

Tuttavia questa funzione, per quanto continua, non è sempre derivabile.

13.1. Definizione. Sia d ≥ 1: S_{d}(x) è la funzione spline di ordine d associata ai dati x_{0}, ..., x_{n} se

1) S_{d}(x) è un polinomio di grado d in ciascun intervallo [x_{i}, x_{i+1}];

2) S_{d}(x) ∈ C^{d−1}(a, b).

3) S_{d}(x_{i}) = f(x_{i}).

Ovviamente la derivata di una spline di ordine k è una spline di ordine k − 1 e la primitiva di una spline di ordine k è una spline di ordine k + 1.

Consideriamo ad esempio le spline cubiche: dati x_{0}, ..., x_{n} ed i valori assunti da f in questi punti: y_{0}, ..., y_{n} dobbiamo avere che, in ciascun tratto [x_{i}, x_{i+1}]

S_{3}(x) = a_{i} + b_{i}x + c_{i}x^{2} + d_{i}x^{3}

e che per ogni i = 0, ..., n:

lim_{x→x_{i}^{−}} \frac{S_{3}^{(k)}(x) − S_{3}^{(k)}(x)}{x − x_{i}} = lim_{x→x_{i}^{+}} \frac{S_{3}^{(k)}(x) − S_{3}^{(k)}(x)}{x − x_{i}}

ove k < 3 − 1 = 2. Ovviamente anche le S_{3}(x_{i}) = y_{i} vanno imposte.

Siano

M_{i} = S′′_{3}(x_{i})

In ciascun intervallo [x_{i}, x_{i+1}] S′′_{3} è una funzione lineare, precisamente

S′′_{3}(x) = \frac{(x_{i+1} − x)M_{i} − (x − x_{i})M_{i+1}}{x_{i+1} − x_{i}}

Infatti S′′_{3}(x) in [x_{i}, x_{i+1}] rappresenta la retta per i punti (x_{i}, M_{i}) e (x_{i+1}, M_{i+1}) cioè la retta di equazione

\frac{y − M_{i}}{M_{i+1} − M_{i}} = \frac{x − x_{i}}{x_{i+1} − x_{i}}

Quindi

y = \frac{(x − x_{i})(M_{i+1} − M_{i}) + M_{i}(x_{i+1} − x_{i})}{x_{i+1} − x_{i}} = \frac{(x − x_{i})M_{i+1} + M_{i}(x_{i+1} − x)}{x_{i+1} − x_{i}}

44 §1.13. Funzioni polinomiali a tratti

cioè l’equazione precedente, dato che M_{i} = S′′_{3}(x_{i}) e y = S′′_{3}(x).

A questo punto identifichiamo questa y con S′′_{3}(x) = 2c_{i} + 6d_{i}x ottenendo

c_{i} = \frac{1}{2} \frac{M_{i}x_{i+1} − M_{i+1}x_{i}}{x_{i+1} − x_{i}}

e

d_{i} = \frac{1}{6} \frac{M_{i+1} − M_{i}}{x_{i+1} − x_{i}}

da cui

S_{3}(x) = a_{i} + b_{i}x + \frac{1}{2} \frac{M_{i}x_{i+1} − M_{i+1}x_{i}}{x_{i+1} − x_{i}} x^{2} + \frac{1}{6} \frac{M_{i+1} − M_{i}}{x_{i+1} − x_{i}} x^{3}

Ora imponiamo che S_{3}(x_{i}) = y_{i} e S_{3}(x_{i+1}) = y_{i+1} ottenendo il sistema lineare

\{ \begin{matrix}y_{i} = a_{i} + b_{i}x_{i} + \frac{1}{2} \frac{M_{i}x_{i+1} − M_{i+1}x_{i}}{x_{i+1} − x_{i}} x_{i}^{2} + \frac{1}{6} \frac{M_{i+1} − M_{i}}{x_{i+1} − x_{i}} x_{i}^{3} \\ y_{i+1} = a_{i} + b_{i}x_{i+1} + \frac{1}{2} \frac{M_{i}x_{i+1} − M_{i+1}x_{i}}{x_{i+1} − x_{i}} x_{i+1}^{2} + \frac{1}{6} \frac{M_{i+1} − M_{i}}{x_{i+1} − x_{i}} x_{i+1}^{3} \end{matrix}

che consente di determinare a_{i} e b_{i}.

Infine si debbono trovare gli M_{i}, e questo si fa imponendo la continuità delle derivate prime delle S_{3}(x), cioè delle

S′_{3}(x) = b_{i} + \frac{M_{i}x_{i+1} − M_{i+1}x_{i}}{x_{i+1} − x_{i}} x + \frac{1}{2} \frac{M_{i+1} − M_{i}}{x_{i+1} − x_{i}} x^{2}

Si deve imporre che

lim_{x→x_{i}^{−}} \frac{S′_{3}(x) − S′_{3}(x)}{x − x_{i}} = lim_{x→x_{i}^{+}} \frac{S′_{3}(x) − S′_{3}(x)}{x − x_{i}}

cioè che

lim_{x→x_{i}^{−}} (b_{i−1} + \frac{M_{i−1}x_{i} − M_{i}x_{i−1}}{x_{i} − x_{i−1}} x + \frac{1}{2} \frac{M_{i} − M_{i−1}}{x_{i} − x_{i−1}} x^{2}) =

= lim_{x→x_{i}^{+}} (b_{i} + \frac{M_{i}x_{i+1} − M_{i+1}x_{i}}{x_{i+1} − x_{i}} x + \frac{1}{2} \frac{M_{i+1} − M_{i}}{x_{i+1} − x_{i}} x^{2})

Infine resteranno due parametri incogniti: M_{0} e M_{n}; per determinarli ci sono vari metodi.

1. Si può ad esempio imporre la coincidenza delle derivate prime agli estremi dell’intervallo:

S′_{3}(x_{0}) = f′(x_{0})  e  S′_{3}(x_{n}) = f′(x_{n})

A conti fatti si ottiene il sistema lineare tridiagonale simmetrico

( \begin{matrix}2(x_{1}−x_{0}) & x_{1}−x_{0} & 0 & ... & 0 & 0 & 0 \\ x_{1}−x_{0} & 2(x_{2}−x_{0}) & x_{2}−x_{1} & ... & 0 & 0 & 0 \\ 0 & x_{2}−x_{1} & 2(x_{3}−x_{2}) & ... & 0 & 0 & 0 \\ ⋮ & ⋮ & ⋮ & ... & ⋮ & ⋮ & ⋮ \\ 0 & 0 & 0 & ... & x_{n−1}−x_{n−2} & 2(x_{n−1}−x_{n−3}) & x_{n}−x_{n−1} \\ 0 & 0 & 0 & ... & 0 & x_{n}−x_{n−1} & 2(x_{n}−x_{n−1}) \end{matrix} ) ×

× ( M_{0} ; M_{1} ; M_{2} ; ⋮ ; M_{n−1} ; M_{n} ) = ( \frac{y_{1}−y_{0}}{x_{1}−x_{0}} −f′(x_{0}) ; \frac{y_{2}−y_{1}}{x_{2}−x_{1}} − \frac{y_{1}−y_{0}}{x_{1}−x_{0}} ; \frac{y_{3}−y_{2}}{x_{3}−x_{2}} − \frac{y_{2}−y_{1}}{x_{2}−x_{1}} ; ⋮ ; \frac{yn−y_{n−1}}{xn−x_{n−1}} − \frac{y_{n−1}−y_{n−2}}{x_{n−1}−x_{n−2}} ; f′(x_{n})− \frac{yn−y_{n−1}}{xn−x_{n−1}} )

Questo è un sistema non singolare a diagonale dominante che si può risolvere con Gauss. 2. Analogamente possiamo imporre la coincidenza delle derivate seconde:

M_{0} = S′′_{3}(x_{0}) = f′′(x_{0})  e  M_{n} = S′′_{3}(x_{n}) = f′′(x_{n})

Se si impone che

M_{0} = M_{n} = 0

le splines così ottenute si dicono naturali.

3. Se il fenomeno è periodico, cioè y_{0} = y_{n} si impone di solito

S′_{3}(x_{0}) = S′_{3}(x_{n})  e  S′′_{3}(x_{0}) = S′′_{3}(x_{n})

Fra tutte le funzioni f ∈ C^{2}(a, b) tali che

1) f(x_{i}) = y_{i} per i = 0, ..., n;

2) f′(x_{0}) = y′_{0} e f′(x_{n}) = y′_{n}, oppure f′′(x_{0}) = f′′(x_{n}) = 0 oppure f′(x_{0}) = f′(x_{n}) e f′′(x_{0}) = f′′(x_{n})

le splines cubiche sono quelle che minimizzano il funzionale

∫_{x_{0}}^{x_{n}}(f′′(x))^{2}dx

Le splines cubiche naturali minimizzano questo funzionale nella classe delle funzioni C^{2}(a, b) tali che f(x_{i}) = y_{i}.

§1.14. Minimi quadrati

Supponiamo di avere una quantità di m + 1 dati {x_{i}} (non necessariamente distinti) ed i corrispondenti valori {y_{i} = f(x_{i})}, e di voler approssimare f attingendo da un insieme di n + 1 funzioni {φ_{i}} con n ≪ m, per combinazione lineare

f_{n}(x) = ∑_{i=0}^{n} c_{i}φ_{i}(x)

a coefficienti costanti determinati col metodo dei minimi quadrati, cioè in modo che

ε^{2} = ∑_{i=0}^{m}(y_{i} − f_{n}(x_{i}))^{2}

46 §1.14. Minimi quadrati

sia minimo.

Se r_{i} = y_{i} − f_{n}(x_{i}) sono gli errori, i coefficienti c_{i} sono determinati come soluzioni del sistema

( \begin{matrix}φ_{0}(x_{0}) & φ_{1}(x_{0}) & ... & φ_{n}(x_{0}) \\ φ_{0}(x_{1}) & φ_{1}(x_{1}) & ... & φ_{n}(x_{1}) \\ ⋮ & ⋮ & ... & ⋮ \\ φ_{0}(x_{m}) & φ_{1}(x_{m}) & ... & φ_{n}(x_{m}) \end{matrix} ) ( c_{0} ; c_{1} ; ⋮ ; c_{n} ) + ( r_{0} ; r_{1} ; ⋮ ; r_{m} ) = ( y_{0} ; y_{1} ; ⋮ ; y_{m} )

in modo che la funzione \sqrt{∑_{i=0}^{m} r_{i}^{2}} sia minimizzata.

Si noti che m > n quindi ci sono più soluzioni che incognite e non si può sperare di trovare una soluzione con r = 0.

Abbiamo dunque una matrice A (n + 1) × (n + 1), un vettore b in ℝ^{n+1} e vogliamo determinare un vettore x ∈ ℝ^{n+1} tale che

ε^{2} = ∥b − Ax∥_{2}^{2}

sia minimo.

Calcoliamo allora i punti estremali di ε^{2} risolvendo il sistema

\frac{∂ε^{2}}{∂c_{i}} = 0

cioè

\{ \begin{matrix}∑_{j=0}^{n} ∑_{k=0}^{m} φ_{0}(x_{k})φ_{j}(x_{k})c_{j} = ∑_{k=0}^{m} y_{k}φ_{0}(x_{k}) \\ ∑_{j=0}^{n} ∑_{k=0}^{m} φ_{1}(x_{k})φ_{j}(x_{k})c_{j} = ∑_{k=0}^{m} y_{k}φ_{1}(x_{k}) \\ ............................................................. \\ ∑_{j=0}^{n} ∑_{k=0}^{m} φ_{n}(x_{k})φ_{j}(x_{k})c_{j} = ∑_{k=0}^{m} y_{k}φ_{n}(x_{k}) \end{matrix}

Notiamo che questo sistema ha come matrice dei coefficienti B = A^{T} A c come termine noto d = A^{T} b.

Dunque la matrice del sistema è simmetrica, e definita positiva se A ha le colonne linearmente indipendenti: in questo caso la soluzione è ovviamente

x = (A^{T} A)^{−1}A^{T} b

ed il residuo r = b − Ax_{0} di un’approssimazione ottima x_{0} è ortogonale allo spazio vettoriale generato dalle colonne di A:

A^{T} r = A^{T} b − A^{T} Ax_{0} = A^{T} b − A^{T} A(A^{T} A)^{−1}A^{T} b = 0

Vale il 14.1. Teorema. Il problema dei minimi quadrati ha sempre soluzione, e questa è unica se le colonne di A sono linearmente indipendenti.

Il problema è quindi trovare effettivamente una tale soluzione: ad esempio si potrebbe usare la decomposizione QR di A; supponiamo che le colonne di A siano linearmente indipendenti: al più con n + 1 trasformazioni di Householder possiamo trasformare A in una matrice

A′ = Q^{T} A = ( \begin{matrix}r_{11} & r_{12} & r_{13} & ... & r_{1n+1} \\ 0 & r_{22} & r_{23} & ... & r_{2n+1} \\ 0 & 0 & r_{33} & ... & r_{3+n1} \\ ⋮ & ⋮ & ⋮ & ... & ⋮ \\ 0 & 0 & 0 & ... & r_{n+1,n+1} \\ 0 & 0 & 0 & ... & 0 \\ ⋮ & ⋮ & ⋮ & ... & ⋮ \\ 0 & 0 & 0 & ... & 0 \end{matrix} )

ove la matrice triangolare superiore è non singolare se e solo se A ha le colonne linearmente indipendenti.

Dato che Q è ortogonale, abbiamo

∥b − Ax∥ = ∥Q^{T} (b − Ax)∥ = ∥b′ − A′x∥

che è minimizzata se Rx = c′ ove b′ − A′x = ( \begin{matrix}c′ − Rx \\ d′ \end{matrix} ) (essendo R la matrice trian-

golare superiore di cui sopra) ed assume in tal caso valore minimo ∥d′∥.