Calcolo numerico

Quando il numero di condizionamento diverge, la soluzione calcolata in macchina può essere

molto diversa dalla soluzione reale.

Osservazione 6. La soluzione in macchina del sistema reale = , in realtà, si riconduce alla

soluzione di un sistema perturbato, definito come

( + )( + ) = + (46)

La matrice e il termine noto sono alterati dall’approssimazione dovuta alla necessità di

rappresentare i numeri nel calcolatore. I termini e sono piccoli, generalmente dell’ordine

di 10 , ma vanno considerati, in quanto dipendono dagli errori di approssimazione.

-16 15

piccole

Se le perturbazioni sui dati e sono , il metodo numerico utilizzato è stabile se anche

la perturbazione sulla soluzione ha lo stesso ordine di grandezza. Il teorema di stabilità

afferma rigorosamente quanto accennato.

Teorema 8. Teorema di stabilità.

Data la matrice , non singolare, sia la perturbazione sui dati contenuta, ovvero tale che valga

−1

‖ ‖ < 1 , allora l’errore sulla soluzione è limitato da un numero che modula le

perturbazioni sui dati. ()

2

≤( )( + ) (47)

()

‖‖ ‖‖ ‖‖

()

1 − 2 ‖‖

Nel caso particolare in cui = 0, ovvero se la matrice A è rappresentata in modo esatto, si ha

che l’errore sulla soluzione dipende linearmente dal numero di condizionamento di A.

≤ () (48)

2

‖‖ ‖‖

Tanto più aumenta l’ordine di () e tanto peggiore è la soluzione, in quanto ad ogni aumento

2

di ordine del numero di condizionamento si perde una cifra significativa della soluzione.

Questo teorema è utilizzato per stimare l’errore dato dai criteri d’arresto negli algoritmi basati

su metodi iterativi. Si definisce come criterio d’arresto per metodi iterativi la condizione sul

residuo normalizzato ‖ ‖

≤

‖‖

Sfruttando il teorema di stabilità, si ha che l’errore al passo è dato dalla differenza fra la

() ()

soluzione esatta e la soluzione al passo , ovvero = − . L’errore soddisfa la condizione

‖ ‖ ‖ ‖

≤ () ∙ (49)

2

‖‖ ‖‖

La tolleranza è scelta arbitrariamente, in base alla complessità del problema. Il criterio

d’arresto, di per sé, garantisce che l’errore sia piccolo nel momento in cui la tolleranza è di ordine

infinitesimo; se la matrice A è malcondizionata, però, l’ordine di grandezza del residuo può

divergere e fornire una soluzione inesatta. 16

Dimostrazione del teorema di stabilità nel caso particolare = 0. Il sistema perturbato ha

equazione ( + ) = + (50)

−1

Essendo = e moltiplicando a sinistra per , si ottiene

= (51)

−1

= (52)

Passando alla norma, ovvero scrivendo l’uguaglianza sulle norme, sfruttando il fatto che la

norma del prodotto fra una matrice e vettore è sempre minore del prodotto fra le relative norme,

si può scrivere l’equazione (53). −1 −1

‖‖ ‖ ‖ ‖‖‖

= ‖ ≤ (53)

‖‖ ‖‖ ‖‖‖‖

Divideno a sinistra e destra per la norma di x e osservando che = ≤ si può

scrivere −1

‖‖ ‖ ‖‖‖

≤ (54)

‖‖ ‖‖ −1

‖‖ ‖‖ ‖‖‖ ‖‖‖

1 ≤ → ≤ (55)

‖‖ ‖‖ ‖‖ ‖‖

Grazie alla sostituzione, nella (55) compare il numero di condizionamento della matrice A. Si

ottiene infine cvd

≤ () ( + )

2

‖‖ ‖‖ ‖‖

−1 −1

Osservazione 7. Nella pratica, si risolve il sistema = e quindi si può scegliere la

−1

(

matrice di precondizionamento P in modo tale che valga ) ≪ (). Questa operazione,

2 2

in realtà, ha un costo e non prevede il calcolo di una matrice inversa, ma si risolve mediante il

sistema lineare definito dal precondizionatore e dal residuo, ovvero = .

17

3. Problema agli autovalori

Data una matrice ∈ ℝ , si devono determinare lo scalare ∈ ℂ e il vettore ∈ ℝ tali che

l’autovettore sia parallelo al vettore , ovvero sia = .

Il problema agli autovalori è un problema tipico non solo di sistemi lineari, ma risolve anche

equazioni differenziali.

Per il problema agli autovalori esiste la soluzione analitica, che prevede la risoluzione del

polinomio caratteristico associato alla matrice , di cui l’autovalore è radice.

()

= det( − ) = 0

Il polinomio carattistico () è un polinomio di ordine a coefficienti reali, ovvero definito da

=1

() ∑

= , ove ∈ ℝ, ∀. Da ciò, dato che la matrice A è reale, deriva che il polinomio

caratteristico ha esattamente soluzioni, reali o complesse coniugate. Tuttavia, ciò non esclude

l’esistenza di autovalori coincidenti, ovvero aventi moltepliticità algebrica maggiore di uno.

Una volta determinati gli autovalori, bisogna determinare gli autovettori, che non esistono

base

univoci: generalmente, si determinano gli autovalori a norma unitaria, definendo così una .

La soluzione analitica del problema agli autovalori comporta il calcolo del determinante della

matrice − , operazione il cui tempo di calcolo incrementa notevolmente (~! operazioni) al

crescere dell’ordine . Il polinomio caratteristico, inoltre, ha soluzione esatta solamente per

< 5, mentre per ordini maggiori bisogna ricorrere ad una soluzione numerica. Per questi

motivi, sono stati implementati metodi numerici per la risoluzione in macchina di problemi agli

autovalori.

Casi particolari. Esistono matrici particolari per cui gli autovalori sono facilmente ricavabili,

anche per ordine elevato. Infatti, nel caso in cui la matrice sia diagonale o triangolare, gli

autovalori sono semplicemente dati dai termini sulla diagonale di . Una matrice , inoltre, è

−1

diagonalizzabile nel momento in cui esiste una matrice U tale che = , ove la matrice

= ( , … , ) è una matrice diagonale che contiene gli autovalori di .

1

3.1 Metodo delle potenze

Il metodo delle potenze fornisce una stima dell’autovalore a modulo massimo della matrice .

Questo metodo è basato su due ipotesi di partenza: gli autovalori sono ordinati in modulo,

| | | |

ovvero vale > ≥ ⋯ ≥ | |, e inoltre la molteplicità di deve essere unitaria; da qui

1 2 1

deriva che la prima disuguaglianza è scritta con il segno di maggiore strettamente. Il metodo

propone un algoritmo iterativo, definito di seguito.

18

(0) (0)

Algoritmo di calcolo. Dato il vettore ∈ ℝ , sia l’autovettore a norma unitaria associato.

λ 2

Ogni iterazione costa circa operazioni e, se il numero massimo di iterazioni, , è minore

di , si tratta di un costo accettabile per l’esecuzione dell’algoritmo.

(0)

Si osserva che il vettore è descrivibile come combinazione lineari degli autovettori.

(0)

= ∑ (1)

=1

1

(0)

= ∑ (2)

(0)

‖ ‖ =1

(1) (0)

Alla prima iterazioni, ovvero per = 1, = , ove y è il vettore definito dalla (2), a cui

(0)

può essere applicata la definizione di autovettore e dunque si ottiene

1 1

(1) (0)

= = ∑ = ∑ (3)

(0) (0)

‖ ‖ ‖ ‖

=1 =1

(1)

(1)

= (4)

(1)

‖ ‖

Al generico passo si può dedurre che il vettore normalizzato sarà dato dalla formula

1

()

= ∑ (5)

=1 ()

∏ ‖ ‖ =1

−1

() =1 ()

∏

Sia [ ] = ‖ ‖

. Riscrivendo la (5) e considerando il limite → +∞, si ottiene

()

() ()

= + ∑ ( )

[ ] (6)

1 1

1

=2 1

Dal momento che, per ipotesi, l’autovalore è quello a modulo massimo, la sommatoria

1 ()

converge a zero. Per → ∞, inoltre, si ha che → , ovvero al crescere del numero di

1

()

iterazioni, l’autovettore normalizzato si allinea all’autovettore , relativo all’autovalore .

1 1

19

Si dimostra che l’errore di convergenza dipende dal rapporto fra e , sempre minore di uno.

1 2

Se è molto maggiore di tutti gli altri autovalori allora l’algoritmo converge velocemente,

1 ( )

invece se = non converge mai, da qui deve essere valida l’ipotesi = 1.

1 2 1

2

() (−1)

‖ − ‖ ≤ ‖ − ‖

| | (7)

1 1

1

Il criterio d’arresto viene scelto, invece, valutando l’errore normalizzato rispetto ad una

() (1)

tolleranza arbitraria sull’autovalore , che converge a .

() (−1)

| − | (8)

<

()

3.2 Metodo delle potenze inverse −1

Il metodo applicato alla matrice inversa permette di stimare l’autovalore a modulo minimo

di . ()

L’algoritmo per questo metodo è analogo, fatta eccezione per la definizione del vettore , che

−1

prevede il calcolo della matrice inversa , operazione non effettuabile in macchina.

Applicando i metodi per la risoluzione di sistemi lineari, è possibile ottenere una

()

approssimazione numerica della soluzione. Il vettore , dunque, deriva dalla soluzione del

() (−1)

sistema lineare = . Si n