Metodi iterativi per la soluzione di sistemi lineari

Alvise Sommariva

Università degli Studi di Padova

Dipartimento di Matematica

8 aprile 2014

Metodi iterativi e diretti

Supponiamo che siano A ∈ ℝ^{n×n} (matrice non singolare), b ∈ ℝ^{n} (vettore colonna) e di dover risolvere il problema Ax = b avente sol. unica x^{∗}.

A tal proposito si può utilizzare la fatt. LU con pivoting. Il costo computazionale è in generale di O(n^{3}/3) operazioni moltiplicative. Questo diventa proibitivo se n è particolarmente elevato.

L’idea dei metodi iterativi è quello di ottenere una successione di vettori x^{(k)} → x^{∗} cosicchè per k̄ ≪ n sia x^{(k̄)} ≈ x^{∗}.

In generale, la soluzione non è ottenuta esattamente come nei metodi diretti in un numero finito di operazioni (in aritmetica esatta), ma quale limite, accontentandosi di poche iterazioni ognuna dal costo quadratico. Quindi il costo totale sarà di ordine O(k̄ · n^{2}).

Metodi iterativi stazionari

Supponiamo che la matrice non singolare A ∈ ℝ^{n×n} sia tale che

A = M − N, con M non singolare.

Allora Mx − Nx = Ax = b e quindi Mx = Nx + b. Moltiplicando ambo i membri per M^{−1} e posto ϕ(x) = M^{−1}Nx + b abbiamo x = M^{−1}Nx + M^{−1}b = ϕ(x). Viene quindi naturale utilizzare la succ. del metodo di punto fisso

x^{(k+1)} = ϕ(x^{(k)}) = M^{−1}Nx^{(k)} + M^{−1}b.

La matrice P = M^{−1}N si dice di iterazione e non dipende, come pure b dall’indice di iterazione k. Per questo motivo tali metodi si chiamano iterativi stazionari. Quale utile notazione, sia inoltre A = D − E − F con D la matrice diagonale estratta da A, E, F rispettivamente triangolari inferiore e superiore.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  3/ 81

Esempio A, E, F

>> A=[1 2 3 4 ; 5 6 7 2 ; 8 9 1 2 ; 3 4 5 1 ]

A =

1  2  3  4
5  6  7  2
8  9  1  2
3  4  5  1
>> E=−( t r i l ( A )−d i a g ( d i a g ( A ) ) )

E =

0  0  0  0
−5  0  0  0
−8  −9  0  0
−3  −4  −5  0
>> F=−( t r i u ( A )−d i a g ( d i a g ( A ) ) )

F =

0  −2  −3  −4
0  0  −7  −2
0  0  0  −2
0  0  0  0
>> % A=d i a g ( d i a g (A) )−E−F .

Metodo di Jacobi

Nel caso del metodo di Jacobi si ha

M = D, N = E + F  (1)

e quindi

P = M^{−1}N = D^{−1}(E + F ) = D^{−1}(D − D + E + F )
= D^{−1}(D − A) = I − D^{−1}A  (2)

Si osservi che se D è singolare allora il metodo di Jacobi, almeno in questa versione di base, non può essere utilizzato visto che in (5) non ha senso la scrittura D^{−1}. Qualora sia a_{ii} ≠ 0 per ogni i = 1, …, n, il metodo di Jacobi può essere descritto come metodo delle sostituzioni simultanee

x_{i}^{(k+1)} = (b_{i} − ∑_{j=1}^{i−1}a_{ij}x_{j}^{(k)} − ∑_{j=i+1}^{n}a_{ij}x_{j}^{(k)})/a_{ii}, i = 1, …, n. (3)

Metodo di Jacobi

Vediamo perchè questa derivazione. Da Ax = b, con A ∈ ℝ^{n×n} abbiamo

∑_{j<i}a_{i,j}x_{j} + a_{i,i}x_{i} + ∑_{j>i}a_{i,j}x_{j} = ∑_{j}a_{i,j}x_{j} = b_{i}, i = 1, …, n

e quindi evidenziando x_{i}, supposto a_{i,i} ≠ 0 deduciamo

x_{i} = \frac{1}{a_{i,i}} (b_{i} − (∑_{j<i}a_{i,j}x_{j} + ∑_{j>i}a_{i,j}x_{j} )), i = 1, …, n

Quindi note delle approssimazioni della soluzione x^{(k)} = {x_{j}^{(k)}} per j = 1, …, n è naturale introdurre il metodo iterativo

x_{i}^{(k+1)} = \frac{1}{a_{i,i}} (b_{i} − (∑_{j<i}a_{i,j}x_{j}^{(k)} + ∑_{j>i}a_{i,j}x_{j}^{(k)})), i = 1, …, n

Metodo di Gauss-Seidel

Il metodo di Gauss-Seidel è definito quale metodo staz. in cui

M = D − E, N = F  (4)

e quindi

P = M^{−1}N = (D − E)^{−1}F  (5)

Similmente al metodo di Jacobi, possiamo riscrivere più semplicemente anche Gauss-Seidel come

x_{i}^{(k+1)} = (b_{i} − ∑_{j=1}^{i−1}a_{ij}x_{j}^{(k+1)} − ∑_{j=i+1}^{n}a_{ij}x_{j}^{(k)})/a_{ii}.  (6)

Da (6) si capisce perchè tale metodo è noto anche come metodo delle sostituzioni successive.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  7/ 81

Metodo di Gauss-Seidel

Vediamo perchè questa derivazione. Da Ax = b, con A ∈ ℝ^{n×n} abbiamo

∑_{j<i}a_{i,j}x_{j} + a_{i,i}x_{i} + ∑_{j>i}a_{i,j}x_{j} = ∑_{j}a_{i,j}x_{j} = b_{i}, i = 1, …, n

e quindi evidenziando x_{i}, supposto a_{i,i} ≠ 0 deduciamo

x_{i} = \frac{1}{a_{i,i}} (b_{i} − (∑_{j<i}a_{i,j}x_{j} + ∑_{j>i}a_{i,j}x_{j} )), i = 1, …, n

Quindi note delle approssimazioni della soluzione x^{(k)} = {x_{j}^{(k)}} per

j > i, e già calcolate x^{(k+1)} = {x_{j}^{(k+1)}} per j < i è naturale introdurre il metodo iterativo

x_{i}^{(k+1)} = \frac{1}{a_{i,i}} (b_{i} − (∑_{j<i}a_{i,j}x_{j}^{(k+1)} + ∑_{j>i}a_{i,j}x_{j}^{(k)})), i = 1, …, n.

SOR

La versione di Gauss-Seidel con la scelta del parametro ω è nota in letteratura come successive over relaxation (SOR)

x^{(k+1)} = (\frac{D}{ω} − E)^{−1} ((\frac{1}{ω} − 1) D + F )x^{(k)} + (\frac{D}{ω} − E)^{−1} b.

(7)

dove D, E, F sono al solito tre matrici tali cha A = D − E − F con D la matrice diagonale estratta da A, E, F rispettivamente triangolari inferiore e superiore.

SOR

Il metodo SOR deriva dallo descrivere esplicitamente una iter. di Gauss-Seidel

\overline{x}_{i}^{(k+1)} = \frac{1}{a_{ii}} [b_{i} − ∑_{j=1}^{i−1}a_{ij}x_{j}^{(k+1)} − ∑_{j=i+1}^{n}a_{ij}x_{j}^{(k)}]

e sostituirla con la combinazione convessa

x_{i}^{(k+1)} = ω\overline{x}_{i}^{(k+1)} + (1 − ω)x^{(k)}.

cioè

x_{i}^{(k+1)} = \frac{ω}{a_{ii}} [b_{i} − ∑_{j=1}^{i−1}a_{ij}x_{j}^{(k+1)} − ∑_{j=i+1}^{n}a_{ij}x_{j}^{(k)}] + (1 − ω)x^{(k)}.

Gauss-Seidel è SOR per ω = 1.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  10/ 81

SOR

Notiamo che le iterazioni di SOR verificano l’uguaglianza

x^{(k+1)} = M^{−1}Nx^{k} + M^{−1}b

con

M = \frac{D}{ω} − E, N = (\frac{1}{ω} − 1) D + F

ed è

M − N = (\frac{D}{ω} − E) − ((\frac{1}{ω} − 1) D + F ) = −E + D − F = A,

per cui SOR è un metodo iterativo stazionario.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  11/ 81

I metodi di Richardson stazionari

Definizione

Fissato α, un metodo di Richardson stazionario, con matrice di precondizionamento P, verifica

P(x^{(k+1)} − x^{(k)}) = αr^{(k)}.  (8)

dove

r^{(k)} = b − Ax^{(k)}  (9)

è il residuo alla k-sima iterazione.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  12/ 81

I metodi di Richardson non stazionari

Definizione

Fissati α_{k} dipendenti dalla iterazione k, un metodo di Richardson non stazionario, con matrice di precondizionamento P, verifica

P(x^{(k+1)} − x^{(k)}) = α_{k}r^{(k)}.  (10)

Si osservi che se α_{k} = α per ogni k, allora il metodo di Richardson non stazionario diventa stazionario.

I metodi di Richardson

I metodi di Jacobi e di Gauss-Seidel, SOR, sono metodi iterativi del tipo

Mx^{(k+1)} = Nx^{(k)} + b,  (11)

per opportune scelte delle matrici M (che dev’essere invertibile), N tali che

A = M − N.  (12)

Essendo r^{(k)} = b − Ax^{(k)},

M(x^{(k+1)} − x^{(k)}) = Nx^{(k)} + b − Mx^{(k)} = b − Ax^{(k)} = r^{(k)}. (13)

I metodi di Richardson

Ne consegue che i metodi di Jacobi e di Gauss-Seidel, SOR, verificano

M(x^{(k+1)} − x^{(k)}) = r^{(k)}  (14)

In altri termini sono dei metodi di Richardson sono metodi di Richardson stazionari, con α = 1 e matrice di precondizionamento P = M.

Per quanto riguarda i metodi di Richardson precondizionati e non stazionari, un classico esempio è il metodo del gradiente classico che vedremo in seguito.

Norma di matrici

Sia ρ(A) il massimo degli autovalori in modulo della matrice A (il cosidetto raggio spettrale).

Sia ∥ · ∥ : ℝ^{n} → ℝ_{+} una norma vettoriale. Definiamo norma naturale (in alcuni testi norma indotta) di una matrice A ∈ R^{n×n} la quantità

∥A∥ := sup_{x∈ℝ^{n},x≠0} \frac{∥Ax∥}{∥x∥}.

Si nota subito che questa definizione coincide con quella di norma di un operatore lineare e continuo in spazi normati.

Norma di matrici

Vediamo alcuni esempi (cf. [4, p.24]). Sia x un arbitrario elemento di R^{n}, A ∈ R^{n×n}.

▶ Si definisce ∥x∥_{1} := ∑_{k=1}^{n} |x_{k}| e si dimostra che la norma naturale corrispondente è (cf. [4, p.26])

∥A∥_{1} = max_{j} ∑_{i=1}^{n}|a_{i,j}|.

▶ Si definisce ∥x∥_{∞} := max_{k} |x_{k}| e si dimostra che la norma naturale corrispondente è (cf. [4, p.26])

∥A∥_{∞} = max_{i} ∑_{j=1}^{n}|a_{i,j}|.

▶ Si definisce ∥x∥_{2} := (∑_{k=1}^{n} |x_{k}|^{2})^{1/2} e si dimostra che la norma naturale corrispondente è (cf. [4, p.27])

∥A∥_{2} = ρ^{1/2}(A^{T}A).

Per quanto riguarda un esempio chiarificatore in Matlab/Octave

>> A=[1 5 ; 7 1 3 ]

A =

1  5
7  13
>> norm ( A , 1 )
ans =

18

>> norm ( A , inf )
ans =

20

>> norm ( A , 2 )
ans =

15.5563

>> e i g ( A∗A ’ )
ans =

2 242

>> s q r t (242)
ans =

15.5563

>> raggio_spettrale_A=max( abs ( e i g ( A ) ) )
raggio_spettrale_A =

15.4261 >>

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  17/ 81

Norma di matrici

Si dimostra che (cf. [4, p.28])

Teorema

Per ogni norma naturale ∥ · ∥ e ogni matrice quadrata A si ha ρ(A) ≤ ∥A∥. Inoltre per ogni matrice A di ordine n e per ogni ϵ > 0 esiste una norma naturale ∥ · ∥ tale che

ρ(A) ≤ ∥A∥ ≤ ρ(A) + ϵ.

e inoltre (cf. [4, p.29], [3, p.232])

Teorema

Fissata una norma naturale ∥ · ∥, i seguenti asserti sono equivalenti 1. A^{m} → 0;

2. ∥A^{m}∥ → 0;

3. ρ(A) < 1.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  18/ 81

Sul raggio spettrale

Ricordiamo che il raggio spettrale

ρ(A) = max_{k}(|λ_{k}|)

(con {λ_{k}}_{k=1,...,n} autovalori della matrice A ∈ ℝ^{n×n}) non è una norma. Infatti la matrice

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

ha raggio spettrale nullo, ma non è la matrice nulla. Osserviamo che dagli esempi il raggio spettrale di una matrice A non coincide in generale con la norma 1, 2, ∞, ma che a volte ρ(A) = ∥A∥_{2} come nel caso di una matrice diagonale A (essendo gli autovalori di una matrice diagonale, proprio i suoi elementi diagonali).

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  19/ 81

Teorema di Hensel

Teorema

Se P è diagonalizzabile allora un metodo iterativo stazionario consistente x^{(k+1)} = Px^{(k)} + c converge per ogni vettore iniziale x_{0} se e solo se ρ(P) < 1.

Dimostrazione.

Ricordiamo che A ∈ ℝ^{n×n} è diagonalizzabile se e solo ℝ^{n} possiede una base composta di autovettori di A (cf. [3, p.57]). Consideriamo un metodo iterativo stazionario x^{(k+1)} = Px^{(k)} + c in cui scelto x^{(0)} si abbia

x^{(0)} − x^{∗} = ∑_{s=1}^{n}c_{s}u_{s}

dove {u_{k}}_{k} è una base di autovettori di P avente autovalori {λ_{k}}_{k}.

Teorema di Hensel

Dimostrazione.

Supponiamo |λ_{s}| < 1 per s = 1, …, n. Se il metodo è consistente, cioè

x^{∗} = Px^{∗} + c

abbiamo

x^{(k)} − x^{∗} = x^{(k−1)} + c − (Px^{∗} + c)
= P(x^{(k−1)} − x^{∗}) = … = P^{k}(x^{0} − x^{∗}
= ∑_{s=1}^{n}c_{s}P^{k}u_{s}
= ∑_{s=1}^{n}c_{s}λ_{s}^{k}u_{s}  (15)
Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  21/ 81

Dimostrazione.

Quindi se |λ_{s}^{k}| < 1 per ogni s = 1, …, n e k = 1, 2, …, abbiamo

∥x^{(k)} − x^{∗}∥ = ∥ ∑_{s=1}^{n}c_{s}λ_{s}^{k}u_{s}∥ ≤ ∑_{s=1}^{n}|c_{s}||λ_{s}^{k}|∥u_{s}∥ → 0.

Se invece per qualche k si ha |λ^{k}| ≥ 1 e c_{k} ≠ 0 allora ∥x^{(k)} − x^{∗}∥ non converge a 0 al crescere di k. Infatti, se λ_{l} ≥ 1 è l’autovalore di massimo modulo, posto

x^{(0)} − x^{∗} = u_{l}

abbiamo da (15)

∥x^{(k)} − x^{∗}∥ = ∥λ_{l}^{k}u_{l}∥ = |λ_{l}^{k}|∥u_{l}∥

che non tende a 0. Di conseguenza non è vero che il metodo è convergente per qualsiasi scelta del vettore x^{(0)}.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  22/ 81

Dimostriamo ora una sua generalizzazione, scoperta da Hensel nel 1926.

Teorema

[Hensel] Un metodo iterativo stazionario consistente x^{(k+1)} = Px^{(k)} + c converge per ogni vettore iniziale x_{0} se e solo se ρ(P) < 1.

Si noti che il teorema riguarda la convergenza per ogni vettore iniziale x_{0} ed è quindi di convergenza globale. Inoltre non si richiede che la matrice P sia diagonalizzabile.

Dimostrazione.

La dimostrazione è tratta da [3, p.236]. ⇐ Se ρ(P) < 1, allora x = Px + c ha una e una sola sol. x^{∗}. Infatti,

x = Px + c ⇔ (I − P)x = c

e la matrice I − P ha autovalori 1 − λ_{k} con k = 1, …, n tali che

0 < |1 − |λ_{k}|_{ℂ}|_{ℝ} ≤ |1 − λ_{k}|_{ℂ},

poichè |λ_{k}|_{ℂ} ≤ ρ(P) < 1 e quindi

det(I − P) = ∏_{k=1}^{n}(1 − λ_{k}) ≠ 0,

per cui la matrice I − P è invertibile e il sistema (I − P)x = c ha una e una sola soluzione x^{∗}.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  24/ 81

Dimostrazione.

Sia e(k) = x^{(k)} − x^{∗}. Come stabilito dal Teorema 0.1, sia inoltre una norma naturale ∥ · ∥ tale che

ρ(P) ≤ ∥P∥ = ρ(P) + (1 − ρ(P))/2 < 1.

Essendo x^{(k+1)} = Px^{(k)} + c e x^{∗} = Px^{∗} + c, sottraendo membro a membro le equazioni si ottiene

e^{(k+1)} = Pe^{(k)} = P^{k}e^{(0)}

da cui essendo ∥ · ∥ una norma naturale

∥e^{(k+1)}∥ = ∥Pe^{(k)}∥ = ∥P^{k}e^{(0)}∥ ≤ ∥P^{k}∥∥e^{(0)}∥.  (16)

Poichè il raggio spettrale è minore di 1 dal Teorema 0.2 abbiamo che ∥P^{k}∥ → 0 da cui per (16) necessariamente ∥e^{(k+1)}∥ → 0 e quindi per le proprietà delle norme e^{(k+1)} → 0 cioè x^{(k)} → x^{∗}.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  25/ 81

Dimostrazione.

⇒ Supponiamo che la successione x^{(k+1)} = Px^{(k)} + c converga a x^{∗} per qualsiasi x^{(0)} ∈ ℝ^{n} ma che sia ρ(P) ≥ 1. Sia λmax il massimo autovalore in modulo di P e scegliamo x^{(0)} tale che e^{(0)} = x^{(0)} − x^{∗} sia autovettore di P relativamente all’autovalore λmax. Essendo Pe^{(0)} = λmaxe^{(0)} e e^{(k+1)} = P^{k}e^{(0)} abbiamo che

e^{(k+1)} = λmax^{k}e^{(0)}

da cui, qualsiasi sia la norma ∥ · ∥, per ogni k = 1, 2, … si ha

∥e^{(k+1)}∥ = |λmax^{k}|_{ℂ}∥e^{(0)}∥ ≥ ∥e^{(0)}∥

il che comporta che la successione non è convergente (altrimenti per qualche k sarebbe e^{(k)} < e^{(0)}).

Velocità di convergenza

L’analisi che faremo in questa sezione non è rigorosamente matematica, ciò nonostante permette di capire il legame tra il raggio spettrale della matrice di iterazione P e la riduzione dell’errore.

Si vede facilmente che se intendiamo calcolare x^{∗} tale che Ax^{∗} = b con un metodo stazionario x^{(k+1)} = Px^{(k)} + c, posto e^{(k)} = x^{(k)} − x^{∗} si ha, supposto il metodo stazionario sia consistente, cioè x^{∗} = Px^{∗} + c

e^{(k)} = x^{(k)} − x^{∗}
= (Px^{(k−1)} + c) − (Px^{∗} + c)
= Pe^{(k−1)} = … = P^{k}∥e^{(0)}∥  (17)

e quindi

∥e^{(k)}∥ ≤ ∥P^{k}∥∥e^{(0)}∥,  (18)

Se e^{(k−1)} ≠ 0, la quantità ∥e^{(k)}∥/∥e^{(k−1)}∥ esprime la riduzione dell’errore al k-simo passo e

σ_{k} = (\frac{∥e^{(k)}∥}{∥e^{(k−1)}∥} … \frac{∥e^{(1)}∥}{∥e^{(0)}∥})^{\frac{1}{k}}

la riduzione media per passo dell’errore relativo ai primi k passi (cf. [3, p.239]).

Velocità di convergenza

Si dimostra che

Teorema

Sia A ∈ ℂ^{n×n} e ∥ · ∥ una norma naturale. Allora lim_{k}∥A^{k}∥^{\frac{1}{k}} = ρ(A)

Quindi per k sufficientemente grande si ha

∥P^{k}∥^{1/k} ≈ ρ(P).

Sotto queste ipotesi, se ∥e^{(k+m)}∥ ≈ ∥P^{m}∥∥e^{(k)}∥ abbiamo

∥e^{(k+m)}∥ ≈ ∥P^{m}∥|∥e^{(k)}∥ ≈ ρ^{m}(P)∥e^{(k)}∥  (19)
Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  29/ 81

per cui affinchè

∥e^{(k+m)}∥/∥e^{(k)}∥ ≈ ρ^{m}(P) ≈ ϵ

applicando il logaritmo naturale ad ambo i membri, si vede serve sia,

m log (ρ(P)) ≈ log ϵ

e quindi

m ≈ \frac{log ϵ}{log (ρ(P))}.

Se

R(P) = − log(ρ(P))

è la cosidetta velocità di convergenza asintotica del metodo iterativo relativo a P, si può così stimare che il numero di iterazioni m necessarie per ridurre l’errore di un fattore ϵ relativamente alla k-sima iterazione, cioè affinchè

∥e^{(k+m)}∥/∥e^{(k)}∥ = ϵ.

Si vede facilmente che è circa

m ≈ ⌈\frac{− log(ϵ)}{R(P)}⌉ .

Conseguentemente, da

m ≈ ⌈\frac{− log(ϵ)}{R(P)}⌉

e

R(P) = − log(ρ(P))

se P è la matrice d’iterazione di un metodo stazionario convergente (e consistente), essendo ρ(P) < 1, minore è ρ(P) necessariamente è maggiore R(P) e si può stimare il numero di iterazioni per ridurre l’errore di un fattore. ϵ.

Si desidera quindi cercare metodi con ρ(P) più piccolo possibile.

Alcuni teoremi di convergenza

Definizione

Una matrice A si dice tridiagonale se A_{i,j} = 0 qualora |i − j| > 1. Esempio:

( \begin{matrix}1 & 2 & 0 & 0 & 0 \\ 1 & 5 & 7 & 0 & 0 \\ 0 & 2 & 0 & 9 & 0 \\ 0 & 0 & 4 & 4 & 2 \\ 0 & 0 & 0 & 5 & 3 \end{matrix} )

Alcuni teoremi di convergenza

Teorema

Per matrici tridiagonali A = (a_{i,j}) con componenti diagonali non nulle, i metodi di Jacobi e Gauss-Seidel sono o entrambi convergenti o divergenti e il tasso di convergenza del metodo di Gauss-Seidel è il doppio di quello del metodo di Jacobi.

Il che vuol dire che asintoticamente sono necessarie metà iterazioni del metodo di Gauss-Seidel per ottenere la stessa precisione del metodo di Jacobi.

Nota.

Il teorema si estende a matrici tridiagonali a blocchi, con blocchi diagonali non singolari.

Alcune definizioni

Definizione

La matrice A è a predominanza diagonale (per righe) se per ogni i = 1, …, n risulta

|a_{i,i}| ≥ ∑_{j=1,j≠i}^{n} |a_{i,j}|

e per almeno un indice s si abbia

|a_{s,s}| > ∑_{j=1,j≠s}^{n} |a_{s,j}|.

La matrice A è a predominanza diagonale per colonne se A^{T} a predominanza diagonale per righe.

Alcune definizioni

Definizione

Se

|a_{s,s}| > ∑_{j=1,j≠s}^{n} |a_{s,j}|, s = 1, …, n

allora la matrice A si dice a predominanza diagonale (per righe) in senso stretto.

La matrice A è a predominanza diagonale per colonne in senso stretto se A^{T} a predominanza diagonale per righe in senso stretto.

Ad esempio la matrice

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

è a predominanza diagonale per righe (non stretta).

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  36/ 81

Teoremi di convergenza

Teorema

Sia A una matrice quadrata a predominanza diagonale per righe in senso stretto. Allora il metodo di Jacobi converge alla soluzione di Ax = b, qualsiasi sia il punto x^{(0)} iniziale.

Teorema

Sia A è a predominanza diagonale per righe in senso stretto. Allora il metodo di Gauss-Seidel converge, qualsiasi sia il punto x^{(0)} iniziale.

Tali teoremi valgono anche se A è a predominanza diagonale per colonne in senso stretto.

Alcune definizioni

La matrice A è simmetrica se A = A^{T}. Una matrice A è definita positiva se ha tutti gli autovalori positivi. Ad esempio la matrice

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

è simmetrica e definita positiva:

>> A=[4 −1 0 ; −1 4 −1; 0 −1 4 ] ;
>> e i g ( A ) % AUTOVALORI DI A .
ans =

2.5858 4.0000 5.4142

>>

Alcuni teoremi di convergenza

Ricordiamo che equivalentemente una matrice A è definita positiva se

x^{T}Ax > 0, per ogni x ≠ 0.

Teorema

Sia A simmetrica con elementi diagonali positivi. Allora il metodo di Gauss-Seidel converge se e solo A è definita positiva.

Più in generale,

Teorema

Sia A simmetrica con elementi diagonali positivi. Allora il metodo SOR converge se e solo se 0 < w < 2 e A è definita positiva.

Per una dimostrazione si veda [6, p.215].

Test di arresto

Consideriamo il sistema lineare Ax = b avente un’unica soluzione x^{∗} e supponiamo di risolverlo numericamente con un metodo iterativo stazionario del tipo

x^{(k+1)} = Px^{(k)} + c,

che sia consistente cioè

x^{∗} = Px^{∗} + c.

Desideriamo introdurre un test di arresto che interrompa le iterazioni, qualora una certa quantità relativa al sistema lineare Ax = b e alle iterazioni eseguite, sia al di sotto di una tolleranza ϵ > 0 fissata dall’utente.

Test di arresto: criterio dello step

Posto Δ^{(k)} := x^{(k+1)} − x^{(k)} e e^{(k)} = x^{∗} − x^{(k)}, essendo

e^{(k)} = x^{∗} − x^{(k)} = (Px^{∗} + c) − (Px^{(k−1)} + c)
= P(x^{∗} − x^{(k−1)}) = Pe^{(k−1)}  (20)

abbiamo

∥e^{(k)}∥_{2} = ∥x^{∗} − x^{(k)}∥_{2} = ∥(x^{∗} − x^{(k+1)}) + (x^{(k+1)} − x^{(k)})∥_{2}
= ∥e^{(k+1)} + Δ^{(k)}∥_{2} = ∥Pe^{(k)} + Δ^{(k)}∥_{2}
≤ ∥P∥_{2} · ∥e^{(k)}∥_{2} + ∥Δ^{(k)}∥_{2}  (21)
Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  41/ 81

Test di arresto: criterio dello step

Fissata dall’utente una tolleranza tol, si desidera interrompere il processo iterativo quando ∥x^{∗} − x^{(k)}∥ ≤ tol.

Non disponendo di x^{∗}, il criterio dello step, consiste nell’interrompere il metodo iterativo alla k + 1-sima iterazione qualora ∥x^{(k+1)} − x^{(k)}∥ ≤ tol.

Di seguito desideriamo vedere quando tale criterio risulti attendibile cioè

|x^{(k+1)} − x^{(k)}| ≈ |x^{∗} − x^{(k)}|

Se P è simmetrica, allora esistono una matrice ortogonale U, cioè tale che U^{T} = U^{−1}, e una matrice diagonale a coefficienti reali Λ per cui

P = UΛU^{T}

ed essendo P e Λ simili hanno gli stessi autovalori {λ_{k}}_{k}.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  42/ 81

Di conseguenza, se P è simmetrica

∥P∥_{2} = \sqrt{ρ(PP^{T})} = \sqrt{ρ(UΛU^{T}(UΛU^{T})^{T})}
= \sqrt{ρ(UΛ^{2}U^{T})}  (22)

Essendo UΛ^{2}U^{T} simile a Λ^{2}, UΛ^{2}U^{T} e Λ^{2} hanno gli stessi autovalori uguali a {λ_{k}^{2}}_{k} e di conseguenza lo stesso raggio spettrale, da cui

ρ(UΛ^{2}U^{T}) = ρ(Λ^{2})

e quindi ricaviamo

∥P∥_{2} = \sqrt{ρ(Λ^{2})} = \sqrt{max_{k}|λ_{k}^{2}|}
= \sqrt{(max_{k}|λ_{k}|^{2})} = \sqrt{(max_{k}|λ_{k}|)^{2}}
= max_{k}|λ_{k}| = ρ(P)  (23)

Quindi da (21)

∥e^{(k)}∥_{2} ≤ ∥P∥_{2} · ∥e^{(k)}∥_{2} + ∥Δ^{(k)}∥_{2}
= ρ(P) · ∥e^{(k)}∥_{2} + ∥Δ^{(k)}∥_{2}  (24)

e se ρ(P) < 1, cioè il metodo iterativo stazionario converge per qualsiasi scelta del vettore iniziale, portando ρ(P) · ∥e^{(k)}∥_{2} a primo membro e dividendo per 1 − ρ(P) deduciamo

∥x^{∗} − x^{(k)}∥_{2} = ∥e^{(k)}∥_{2} =≤ \frac{∥Δ^{(k)}∥_{2}}{1 − ρ(P)} = \frac{∥x^{(k+1)} − x^{(k)}∥_{2}}{1 − ρ(P)}

da cui se P è simmetrica allora il criterio dello step è affidabile se ρ(P) è piccolo.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  44/ 81

Test di arresto: criterio del residuo

Si definisce residuo alla k-sima iterazione la quantità

r^{(k)} := b − Ax^{(k)}

ed essendo b = Ax^{∗} abbiamo per e^{(k)} = x^{∗} − x^{(k)}

b − Ax^{(k)} = Ax^{∗} − Ax^{(k)} = A(x^{∗} − x^{(k)}) = Ae^{(k)} da cui

r^{(k)} = Ae^{(k)}.

Interromperemo il processo iterativo quando

∥r^{(k)}∥ ≤ tol ,

desiderando pure

\frac{∥x^{(k)} − x^{∗}∥}{∥x^{∗}∥} ≤ tol

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  45/ 81

Test di arresto: criterio del residuo

Osserviamo che

▶ essendo A invertibile e r^{(k)} = Ae^{(k)} ricaviamo e^{(k)} = A^{−1}r^{(k)}

da cui

∥e^{(k)}∥ = ∥A^{−1}r^{(k)}∥ ≤ ∥A^{−1}∥∥r^{(k)}∥;

▶ poichè b = Ax^{∗} abbiamo ∥b∥ ≤ ∥A∥∥x^{∗}∥ e quindi

\frac{1}{∥x^{∗}∥} ≤ \frac{∥A∥}{∥b∥}.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  46/ 81

Test di arresto: criterio del residuo

Di conseguenza, denotato con κ(A) = ∥A∥∥A^{−1}∥ il numero di condizionamento (necessariamente maggiore o uguale a 1), se x^{∗} ≠ 0 abbiamo

\frac{∥e^{(k)}∥}{∥x^{∗}∥} ≤ \frac{∥A∥}{∥b∥}∥e^{(k)}∥ ≤ \frac{∥A∥}{∥b∥} · ∥A^{−1}∥∥r^{(k)}∥ ≤ κ(A)\frac{∥r^{(k)}∥}{∥b∥}

Quindi

\frac{∥x^{(k)} − x^{∗}∥}{∥x^{∗}∥} = \frac{∥e^{(k)}∥}{∥x^{∗}∥} ≤ κ(A)\frac{∥r^{(k)}∥}{∥b∥}.

Il criterio d’arresto \frac{∥r^{(k)}∥}{∥b∥} ≤ tol è quindi inappropriato quando κ(A) ≫ 1 in quanto non implica necessariamente che

\frac{∥x^{(k)} − x^{∗}∥}{∥x^{∗}∥} ≤ tol.

I metodi di discesa

Sia A una matrice simmetrica definita positiva.

Si osserva che se x^{∗} è l’unica soluzione di Ax = b allora è pure il minimo del funzionale dell’energia

ϕ(x) = \frac{1}{2}x^{T}Ax − b^{T}x, x ∈ ℝ^{n}

in quanto grad(ϕ(x)) = Ax − b = 0 ⇔ Ax = b.

Quindi invece di calcolare la soluzione del sistema lineare, intendiamo calcolare il minimo del funzionale ϕ.

Un generico metodo di discesa consiste nel generare una successione

x^{(k+1)} = x^{(k)} + α_{k}p^{(k)}

dove p^{(k)} è una direzione fissata secondo qualche criterio. Lo scopo ovviamente è che

ϕ(x^{(k+1)}) < ϕ(x^{(k)}),

e che il punto x^{∗}, in cui si ha il minimo di ϕ, venga calcolato rapidamente.

Il metodo del gradiente classico consiste nello scegliere α_{k} e p^{(k)} così da ottenere la massima riduzione del funzionale dell’energia a partire dal punto x^{(k)}.

Differenziando

ϕ(x) = \frac{1}{2}x^{T}Ax − b^{T}x, x ∈ ℝ^{n}

si vede che tale scelta coincide con lo scegliere

p^{(k)} = r^{(k)}
α_{k} = \frac{∥r^{(k)}∥_{2}^{2}}{(r^{(k)})^{T}Ar^{(k)}}.  (25)

Essendo x^{(k+1)} = x^{(k)} + α_{k}r^{(k)} è un metodo di Richardson non stazionario con P = I e parametro α_{k} definito da (25).

Il metodo del gradiente coniugato (1952).

Supponiamo di dover risolvere il sistema lineare Ax = b. Con il termine residuo in x^{(k)} si intende la quantità r^{(k)} = b − Ax^{(k)}. La succ. delle iterazioni del gradiente coniugato è quella propria dei metodi di discesa,

x^{(k+1)} = x^{(k)} + α_{k}p^{(k)}, α_{k} = \frac{(r^{(k)})^{T}r^{(k)}}{(p^{(k)})^{T}Ap^{(k)}}

dove p^{(0)} = r^{(0)} e

p^{(k)} = r^{(k)} + β_{k}p^{(k−1)}, β_{k} = \frac{(r^{(k)})^{T}r^{(k)}}{(r^{(k−1)})^{T}r^{(k−1)}}.

Con questa scelta si prova che p^{(k)} e p^{(k−1)} sono A-coniugati.

(p^{(k)})^{T}Ap^{(k−1)} = 0.

Il metodo del gradiente coniugato: proprietà

Il metodo del gradiente coniugato ha molte proprietà particolari. Ne citiamo alcune.

Teorema

Se A è una matrice simmetrica e definita positiva di ordine n, allora il metodo del gradiente coniugato è convergente e fornisce in aritmetica esatta la soluzione del sistema Ax = b in al massimo n iterazioni.

Questo teorema tradisce un po’ le attese, sia perchè in generale i calcoli non sono compiuti in aritmetica esatta, sia perchè in molti casi della modellistica matematica n risulta essere molto alto.

Il metodo del gradiente coniugato: proprietà

Definizione

Lo spazio

K_{k} = span(r^{(0)}, Ar^{(0)}, …, A^{k−1}r^{(0)})

per k ≥ 1 si dice spazio di Krylov.

Teorema

Sia

K_{k} = span(r^{(0)}, Ar^{(0)}, …, A^{k−1}r^{(0)})

per k ≥ 1. Allora la k-sima iterata dal metodo del gradiente coniugato, minimizza il funzionale ϕ nell’insieme x^{(0)} + K_{k}.

Per una dimostrazione si veda [7, p.12]. Si osservi che essendo la k-sima iterazione del gradiente classico pure in x^{(0)} + K_{k}, il gradiente classico non minimizza in generale il funzionale ϕ in x^{(0)} + K_{k}.

Il metodo del gradiente coniugato: proprietà

Si può dimostrare che se A è simmetrica e definita positiva,

∥x∥_{A} = \sqrt{x^{T}Ax}

e

e_{k} = x^{∗} − x^{(k)}

allora

∥e_{k}∥_{A} ≤ (\frac{\sqrt{K_{2}(A)} − 1}{\sqrt{K_{2}(A)} + 1})^{2k} ∥e_{0}∥_{A}.

Questo risultato stabilisce che la convergenza del gradiente coniugato è lenta qualora si abbiano alti numeri di condizionamento

K_{2}(A) := ∥A∥_{2}∥A^{−1}∥_{2} = \frac{max_{i} |λ_{i}|}{min_{j} |λ_{j}|}

(ove al solito {λ_{i}} sono gli autovalori di A).

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  54/ 81

Metodo di Jacobi in Matlab

Un codice gratuito del metodo di Jacobi, è jacobi.m tratto da Netlib:

f u n c t i o n [ x , e r r o r , iter , f l a g ]= jacobi ( A , x , b , max_it , tol )
% i n p u t
% A, REAL m a t r i x
% x , REAL i n i t i a l g u e s s v e c t o r
% b , REAL r i g h t hand s i d e v e c t o r
% max it , INTEGER maximum number o f i t e r a t i o n s
% t o l , REAL e r r o r t o l e r a n c e

%

% output
% x , REAL s o l u t i o n v e c t o r
% e r r o r , REAL e r r o r norm
% i t e r , INTEGER number o f i t e r a t i o n s performed
% f l a g , INTEGER : 0 = s o l u t i o n found to t o l e r a n c e
% 1 = no c o n v e r g e n c e g i v e n m a x i t

Nota.

Osserviamo che jacobi inizialmente usa un test del residuo (relativo), per poi usarne uno basato sullo step (relativo).

Metodo di Jacobi in Matlab

iter = 0 ; % i n i t i a l i z a t i o n
f l a g = 0 ;
bnrm2 = norm ( b ) ;
i f ( bnrm2 = 0 . 0 ) , bnrm2 = 1 . 0 ; end

r = b − A∗x ;

e r r o r = norm ( r ) / bnrm2 ;
i f ( e r r o r < tol ) r e t u r n , end
[ m , n]= s i z e ( A ) ;
[ M , N ] = split ( A , b , 1 . 0 , 1 ) ; % m a t r i x s p l i t t i n g
f o r iter = 1 : max_it , % b e g i n i t e r a t i o n

  x_1 = x ;

x = M \ ( N∗x + b ) ; % update a p p r o x i m a t i o n
e r r o r = norm ( x − x_1 ) / norm ( x ) ; % compute e r r o r
i f ( e r r o r <= tol ) , break , end % check c o n v e r g e n c e

end

i f ( e r r o r > tol ) f l a g = 1 ; end % no c o n v e r g e n c e

Metodo di SOR in Matlab

f u n c t i o n [ x , e r r o r , iter , f l a g ] = sor ( A , x , b , w , max_it , tol )
% s o r .m s o l v e s the l i n e a r system Ax=b u s i n g the
% S u c c e s s i v e Over−R e l a x a t i o n Method ( Gauss−S e i d e l method when omega = 1 ) .
% i n p u t
% A, REAL m a t r i x
% x REAL i n i t i a l g u e s s v e c t o r
% b REAL r i g h t hand s i d e v e c t o r
% w REAL r e l a x a t i o n s c a l a r
% m a x i t INTEGER maximum number o f i t e r a t i o n s
% t o l REAL e r r o r t o l e r a n c e
% output
% x REAL s o l u t i o n v e c t o r
% e r r o r REAL e r r o r norm
% i t e r INTEGER number o f i t e r a t i o n s performed
% f l a g INTEGER : 0 = s o l u t i o n found to t o l e r a n c e
%  1 = no c o n v e r g e n c e g i v e n m a x i t

Metodo di SOR in Matlab

f l a g = 0 ; % i n i t i a l i z a t i o n

iter = 0 ;

bnrm2 = norm ( b ) ;
i f ( bnrm2 = 0 . 0 ) , bnrm2 = 1 . 0 ; end

r = b − A∗x ;

e r r o r = norm ( r ) / bnrm2 ;
i f ( e r r o r < tol ) r e t u r n , end
[ M , N , b ] = split ( A , b , w , 2 ) ; % m a t r i x s p l i t t i n g
f o r iter = 1 : max_it % b e g i n i t e r a t i o n

 x_1 = x ;

x = M \ ( N∗x + b ) ; % update a p p r o x i m a t i o n
e r r o r = norm ( x − x_1 ) / norm ( x ) ; % compute e r r o r
i f ( e r r o r <= tol ) , break , end % check c o n v e r g e n c e

end

b = b / w ; % r e s t o r e r h s
i f ( e r r o r > tol ) f l a g = 1 ; end ; % no c o n v e r g e n c e

Routine split

La routine split chiamata da jacobi e sor, tratta da Netlib

f u n c t i o n [ M , N , b ] = split ( A , b , w , f l a g )
% s p l i t .m s e t s up the m a t r i x s p l i t t i n g f o r the s t a t .
% i t e r a t i v e methods : j a c o b i and s o r ( gauss−s e i d e l , w=1)
% i n p u t
% A DOUBLE PRECISION m a t r i x
% b DOUBLE PRECISION r i g h t hand s i d e v e c t o r ( f o r SOR)
% w DOUBLE PRECISION r e l a x a t i o n s c a l a r
% f l a g INTEGER f l a g f o r method : 1 = j a c o b i 2 = s o r .
% output
% M DOUBLE PRECISION m a t r i x
% N DOUBLE PRECISION m a t r i x such t h a t A = M − N
% b DOUBLE PRECISION r h s v e c t o r ( a l t e r e d f o r SOR )
[ m , n ] = s i z e ( A ) ;
i f ( f l a g = 1 ) , % j a c o b i s p l i t t i n g
M = d i a g ( d i a g ( A ) ) ; N = d i a g ( d i a g ( A ) ) − A ;
e l s e i f ( f l a g = 2 ) , % s o r / gauss−s e i d e l s p l i t t i n g

   b = w ∗ b ;

M = w ∗ t r i l ( A , −1 ) + d i a g ( d i a g ( A ) ) ;
N = −w ∗ t r i u ( A , 1 ) + ( 1 . 0 − w ) ∗ d i a g ( d i a g ( A ) ) ;
end ;

Gradiente coniugato in Matlab

Un codice gratuito del Grad. Coniugato, è cg.m tratto da Netlib:

f u n c t i o n [ x , e r r o r , iter , f l a g ]=cg ( A , x , b , M , max_it , tol )
f l a g = 0 ; iter = 0 ; bnrm2 = norm ( b ) ;
i f ( bnrm2 = 0 . 0 ) , bnrm2 = 1 . 0 ; end
r = b − A∗x ; e r r o r = norm ( r ) / bnrm2 ;
i f ( e r r o r < tol ) r e t u r n , end

f o r iter = 1 : max_it

z = M \ r ; rho = ( r ’ ∗ z ) ;
i f ( iter > 1 )
beta = rho / rho_1 ; p = z + beta ∗p ;

e l s e

  p = z ;

end

q = A∗p ; alpha = rho / ( p ’ ∗ q ) ; x = x + alpha ∗ p ;

r = r − alpha∗q ; e r r o r = norm ( r ) / bnrm2 ;
i f ( e r r o r <= tol ) , break , end

rho_1 = rho ;

end

i f ( e r r o r > tol ) f l a g = 1 ; end

Esercizio in Matlab: matrice di Poisson

Consideriamo il sistema lineare Ax = b dove A è la matrice tridiagonale a blocchi (di Poisson)

A = ( \begin{matrix}B & −I & 0 & … & 0 \\ −I & B & −I & … & 0 \\ 0 & −I & B & … & … \\ 0 & … & … & … & −I \\ 0 & 0 & … & −I & B \end{matrix} )

con

B = ( \begin{matrix}4 & −1 & 0 & … & 0 \\ −1 & 4 & −1 & … & 0 \\ 0 & −1 & 4 & … & … \\ 0 & … & … & … & −1 \\ 0 & 0 & … & −1 & 4 \end{matrix} )

Esercizio in Matlab: matrice di Poisson

La matrice A è facilmente disponibile, con il comando gallery di Matlab. Vediamo un esempio:

>> A=g a l l e r y ( ’ p o i s s o n ’ , 3 ) ; % A s p a r s e .
>> A= f u l l ( A ) ; % A p i e n a .
>> A

A =

4  −1  0  −1  0  0  0  0  0
−1  4  −1  0  −1  0  0  0  0
0  −1  4  0  0  −1  0  0  0
−1  0  0  4  −1  0  −1  0  0
0  −1  0  −1  4  −1  0  −1  0
0  0  −1  0  −1  4  0  0  −1
0  0  0  −1  0  0  4  −1  0
0  0  0  0  −1  0  −1  4  −1
0  0  0  0  0  −1  0  −1  4

>>

Esercizio in Matlab: matrice di Poisson

Evidentemente A è una matrice di Poisson con B matrice quadrata di ordine 3, dove

A = ( \begin{matrix}B & −I & 0 & … & 0 \\ −I & B & −I & … & 0 \\ 0 & −I & B & … & … \\ 0 & … & … & … & −I \\ 0 & 0 & … & −I & B \end{matrix} )

in cui

B = ( \begin{matrix}4 & −1 & 0 \\ −1 & 4 & −1 \\ 0 & −1 & 4 \end{matrix} )

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  63/ 81

Esercizio in Matlab: matrice di Poisson

Per ulteriori dettagli sulle origini della matrice di Poisson, si considerino ad esempio [1, p. 557], [3, p. 283], [4, p. 334]. Le matrici di Poisson sono

▶ simmetriche; ▶ tridiagonali a blocchi; ▶ diagonalmente dominanti; ▶ non singolari (deriva dal primo e dal secondo teorema di

Gerschgorin [3, p. 76-80], [4, p. 955]); ▶ definite positive.

Per accertarsene, calcoliamo il minimo autovalore della matrice di Poisson con B ∈ M_{5}, semplicemente digitando sulla shell di Matlab-Octave

>> A=makefish ( 5 ) ;

>> m=min ( e i g ( A ) )

m =

0.5359 >>

Tale matrice di Poisson non è malcondizionata essendo

>> A=makefish ( 5 ) ;

>> cond ( A )
ans =

13.9282 >>

Poniamo ora

b=ones ( s i z e ( A , 1 ) , 1 ) ;

e risolviamo il sistema Ax = b digitando

x_sol=A\b ;

Nota la soluzione esatta confrontiamo i vari metodi risolvendo il sistema lineare con un numero massimo di iterazioni maxit e una tolleranza tol come segue

maxit =200; tol=10ˆ (−8) ;

Esercizio in Matlab: matrice di Poisson

A tal proposito consideriamo l’m-file demo algebra lineare.m, contenente il codice

maxit =200; tol=10ˆ (−8) ;

siz =5;

A = makefish ( siz ) ;  % MATRICE DI POISSON .
b=ones ( s i z e ( A , 1 ) , 1 ) ; % TERMINE NOTO.
x_sol=A\b ;  % SOLUZIONE ESATTA . METODO LU .

norm_x_sol=norm ( x_sol ) ;

i f norm ( x_sol ) = 0

     norm_x_sol =1;

end

x=z e r o s ( s i z e ( b ) ) ;  % VALORE INIZIALE .

% JACOBI .

[ x_j , error_j , iter_j , flag_j ]= jacobi ( A , x , b , maxit , tol ) ;

f p r i n t f ( ’\t \n [ JACOBI ] [ STEP REL . , NORMA 2 ] : % 2 . 2 e [ REL . ERR . ] : % 2 . 2 e ’ , error_j ,

      norm ( x_j−x_sol ) / norm_x_sol ) ;

f p r i n t f ( ’\t \n [ ITER . ] : % 3 . 0 f [ FLAG] : % 1 . 0 f \n ’ , iter_j , flag_j ) ;
% GAUSS−SEIDEL .
w=1;

[ x_gs , error_gs , iter_gs , flag_gs ]=sor ( A , x , b , w , maxit , tol ) ;

f p r i n t f ( ’\t \n [ GS ] [ STEP REL . , NORMA 2 ] : % 2 . 2 e [ REL . ERR . ] : % 2 . 2 e ’ , error_gs , norm (

      x_gs−x_sol ) / norm_x_sol ) ;

f p r i n t f ( ’\t \n [ ITER . ] : % 3 . 0 f [ FLAG] : % 1 . 0 f \n ’ , iter_gs , flag_gs ) ;
% SOR .
w_vett = 0 . 8 : 0 . 0 2 5 : 2 ;
f o r index =1: l e n g t h ( w_vett )

     w=w_vett ( index ) ;
    [ x_sor , error_sor ( index ) , iter_sor ( index ) , flag_sor ( index ) ] = sor ( A , x , b , w , maxit
          , tol ) ;
   relerr ( index )=norm ( x_sor−x_sol ) / norm_x_sol ;

end

[ min_iter_sor , min_index ]=min ( iter_sor ) ;

f p r i n t f ( ’\t \n [SOR OTT . ] [ STEP REL . ,NORMA 2 ] : % 2 . 2 e [ REL . ERR . ] : % 2 . 2 e ’ , error_sor (

      min_index ) , relerr ( min_index ) ) ;

f p r i n t f ( ’\t \n [ ITER . ] : % 3 . 0 f [ FLAG] : % 1 . 0 f [w] : % 2 . 3 f \n ’ , min_iter_sor , flag_sor (

      min_index ) , w_vett ( min_index ) ) ;
p l o t ( w_vett , iter_sor , ’ r−’ ) ;

% GRADIENTE CONIUGATO.
M=eye ( s i z e ( A ) ) ;

[ x_gc , error_gc , iter_gc , flag_gc ]=cg ( A , x , b , M , maxit , tol ) ;

f p r i n t f ( ’\t \n [GC ] [ STEP REL . , NORMA 2 ] : % 2 . 2 e [ REL . ERR . ] : % 2 . 2 e ’ , error_gc , norm (

      x_gc−x_sol ) / norm_x_sol ) ;

f p r i n t f ( ’\t \n [ ITER . ] : % 3 . 0 f [ FLAG] : % 1 . 0 f \n ’ , iter_gc , flag_gc ) ;

Esercizio in Matlab: matrice di Poisson

Lanciamo la demo nella shell di Matlab-Octave e otteniamo

>> demo_algebra_lineare

[ JACOBI ] [ STEP  REL . , NORMA 2 ] : 8 . 7 3 e−009 [ REL . ERR . ] : 5 . 6 5 e−008
[ ITER . ] : 116 [ FLAG ] : 0
[ GAU . SEI . ] [ STEP  REL . , NORMA 2 ] : 9 . 2 2 e−009 [ REL . ERR . ] : 2 . 7 6 e−008
[ ITER . ] : 61 [ FLAG ] : 0
[ SOR  OTT . ] [ STEP  REL . , NORMA 2 ] : 2 . 3 1 e−009 [ REL . ERR . ] : 1 . 1 0 e−009
[ ITER . ] : 21 [ FLAG ] : 0 [ w ] : 1.350
[ GC ] [ STEP  REL . , NORMA 2 ] : 4 . 4 1 e−017 [ REL . ERR . ] : 2 . 2 1 e−016
[ ITER . ] :  5 [ FLAG ] : 0

>>

Una breve analisi ci dice che

▶ Come previsto dalla teoria, il metodo di Gauss-Seidel converge

in approssimativamente metà iterazioni di Jacobi; ▶ Il metodo SOR ha quale costante quasi ottimale w = 1.350; ▶ Il metodo del gradiente coniugato converge in meno iterazioni

rispetto agli altri metodi (solo 5 iterazioni, ma si osservi il test d’arresto differente). Essendo la matrice di Poisson di ordine 25, in effetti ciò accade in meno di 25 iterazioni come previsto. Vediamo cosa succede dopo 25 iterazioni:

>> A=g a l l e r y ( ’ p o i s s o n ’ , 5 ) ;
>> A= f u l l ( A ) ; b=ones ( s i z e ( A , 1 ) , 1 ) ;
>> maxit =25; tol =0;
>> x=z e r o s ( s i z e ( b ) ) ; M=eye ( s i z e ( A ) ) ;

>> [ x_gc , error_gc , iter_gc , flag_gc ]=cg ( A , x , b , M , maxit , tol ) ;
>> error_gc
error_gc =

8.3759 e−39

Il residuo relativo, seppur non nullo è molto piccolo.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  71/ 81

Un punto delicato riguarda la scelta del parametro ω ottimale (cioè minimizzante il raggio spettrale di SOR). Sia questo valore uguale a ω^{∗}. Nel nostro codice abbiamo calcolato per forza bruta ω^{+}, tra i numeri reali ω^{+} ≤ 2 del tipo w_{j} = 0.8 + j · 0.025 quello per cui venivano compiute meno iterazioni.

E’ possibile calcolare ω^{∗} matematicamente? Nel caso della matrice di Poisson la risposta è affermativa. Da [4, Teor.5.10, p.333]

ω^{∗} = \frac{2}{1 + \sqrt{1 − ρ^{2}(B_{J})}}

dove ρ(B_{J}) è il massimo degli autovalori in modulo della matrice B_{J} la matrice di iterazione del metodo di Jacobi.

Il raggio spettrale della matrice di iterazione di SOR ott. vale ω^{∗} − 1.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  72/ 81

Esercizio in Matlab: matrice di Poisson

Vediamo di calcolare questo valore nel caso della sopracitata matrice di Poisson. Dalla teoria, con ovvie notazioni,

B_{J} = I − D^{−1}A

e quindi

>> format long ;
>> D=d i a g ( d i a g ( A ) ) ;
>> BJ=eye ( s i z e ( A ) )−i n v ( D ) ∗A ;
>> s=e i g ( BJ ) ;
>> s_abs=abs ( s ) ;
>> rho=max( s_abs ) ;
>> w=2/(1+ s q r t (1−rho ˆ 2) )

 w =

1.33333333333333

>> maxit =50; tol=10ˆ (−8) ;
>> b=ones ( s i z e ( A , 1 ) , 1 ) ;

>> [ x_sor , error_sor , iter_sor , flag_sor ] = sor ( A , x , b , w , maxit , tol ) ;
>> iter_sor
 iter_sor =

22 >>

Si rimane un po’ sorpresi dal fatto che per w = 1.350 il numero di iterazioni fosse inferiore di quello fornito dal valore ottimale teorico w^{∗} = 1.333 . . ..

Il fatto è che questo è ottenuto cercando di massimizzare la velocità asintotica di convergenza. Purtroppo questo minimizza una stima del numero di iterazioni k minime da compiere e non quello effettivo.

Abbiamo detto che un punto chiave è la grandezza del raggio spettrale delle matrici di iterazione e che è desiderabile che questo numero oltre ad essere strettamente minore di uno sia il più piccolo possibile. Vediamo i raggi spettrali dei metodi esposti.

Esercizio in Matlab: matrice di Poisson

Salviamo in raggispettrali.m il seguente programma principale

maxit =50; tol =0;
siz =5;

A = makefish ( siz ) ;  % MATRICE DI POISSON .
b=ones ( s i z e ( A , 1 ) , 1 ) ; % TERMINE NOTO.
[ M , N ] = split ( A , b , 1 . 0 , 1 ) ; % JACOBI .
P=i n v ( M ) ∗N ;
rho_J=max( abs ( e i g ( P ) ) ) ;
f p r i n t f ( ’\n \t [ RAGGIO SPETTRALE ] [ JACOBI ] : %2.15 f ’ , rho_J ) ;
[ M , N , b ] = split ( A , b , 1 , 2 ) ; % GS .
P=i n v ( M ) ∗N ;
rho_gs=max( abs ( e i g ( P ) ) ) ;
f p r i n t f ( ’\n \t [ RAGGIO SPETTRALE ] [ GAUSS−SEIDEL ] : %2.15 f ’ , rho_gs ) ;

D=d i a g ( d i a g ( A ) ) ; E=−( t r i l ( A )−D ) ; F=−( t r i u ( A )−D ) ;
w =1.350;
M=D/w−E ; N=(1/w−1)∗D+F ; P=i n v ( M ) ∗N ;
rho_sor=max( abs ( e i g ( P ) ) ) ;
f p r i n t f ( ’\n \t [ RAGGIO SPETTRALE ] [ SOR BEST] : % 2 . 1 5 f ’ , rho_sor ) ;
w =1.33333333333333;
[ M , N , b ] = split ( A , b , w , 2 ) ; % SOR OPT.
M=D/w−E ; N=(1/w−1)∗D+F ; P=i n v ( M ) ∗N ;
rho_sor_opt=max( abs ( e i g ( P ) ) ) ;
f p r i n t f ( ’\n \t [ RAGGIO SPETTRALE ] [ SOR OPT ] : %2.15 f ’ , rho_sor_opt ) ;

Di seguito:

>> raggispettrali

[ RAGGIO  SPETTRALE ] [ JACOBI ] : 0.866025403784438

    [ RAGGIO  SPETTRALE ] [ GAUSS−SEIDEL ] : 0.750000000000000
    [ RAGGIO  SPETTRALE ] [ SOR  BEST ] : 0.350000000000001

[ RAGGIO  SPETTRALE ] [ SOR  OPT ] : 0.333333380707781

>>

Il valore del raggio spettrale della matrice di iterazione del metodo SOR per parametro ottimale, per quanto visto anticipatamente vale ω^{∗} − 1, e l’esperimento numerico lo conferma. Abbiamo poi osservato che in questo caso la velocità di convergenza del metodo di Gauss-Seidel è il doppio di quella di Jacobi. Poste B_{GS}, B_{J} le rispettive matrici di iterazione, e detta R la velocità di convergenza, osserviamo che da

R(B_{J}) := − ln (ρ(B_{J}))  (26)
R(B_{GS}) := − ln (ρ(B_{GS}))  (27)
R(B_{GS}) := 2R(B_{J})  (28)

si ha

− ln (ρ(B_{GS})) = R(B_{GS}) = 2R(B_{J}) = −2 ln (ρ(B_{J})) = − ln (ρ(B_{J}))^{2}

da cui essendo il logaritmo una funzione invertibile

ρ(B_{GS}) = (ρ(B_{J}))^{2}.
Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  78/ 81

Il raggio spettrale della matrice di iterazione di Gauss-Seidel coincide quindi col quadrato di quella di Jacobi ed infatti come è facile verificare

>> 0.866025403784438ˆ 2
ans =

0.75000000000000 >>

Esercizio in Matlab II

▶ Si calcoli in Matlab la matrice di Poisson P_{20} di ordine 20 aiutandosi con

>> h e l p g a l l e r y

In Octave si usi la routine makefish scaricabile da internet. ▶ Sia b il vettore composto di componenti uguali a 1, avente lo

stesso numero di righe di P_{20}.

1. Si risolva con Gauss-Seidel e Jacobi, il problema P_{20}x = b, con tolleranza di 10^{(−12)}, partendo da x^{0} = [0 … 0] e b = [1, …, 1]. Quante iterazioni servono?

2. Si risolva col gradiente coniugato il problema P_{20}x = b, con tolleranza di 10^{(−12)}, partendo da x^{0} = [0 … 0]. Quante iterazioni servono?

Bibliografia

K. Atkinson, Introduction to Numerical Analysis, Wiley, 1989. K. Atkinson e W. Han, Theoretical Numerical Analysis, Springer, 2001. D. Bini, M. Capovani e O. Menchi, Metodi numerici per l’algebra lineare, Zanichelli, 1988. V. Comincioli, Analisi Numerica, metodi modelli applicazioni, Mc Graw-Hill, 1990. S.D. Conte e C. de Boor, Elementary Numerical Analysis, 3rd Edition, Mc Graw-Hill, 1980. L.A. Hageman e D.M. Young Applied Iterative Methods, Dover, 2004. C.T. Kelley, Iterative Methods for Linear and Nonlinear Equations, SIAM, 1995. The MathWorks Inc., Numerical Computing with Matlab, http://www.mathworks.com/moler. Netlib, http://www.netlib.org/templates/matlab/. A. Quarteroni e F. Saleri, Introduzione al calcolo scientifico, Springer Verlag, 2006. A. Suli e D. Mayers, An Introduction to Numerical Analysis, Cambridge University Press, 2003. Wikipedia (Metodo di Gauss-Seidel), http://it.wikipedia.org/wiki/Metodo di Gauss-Seidel. Wikipedia (Metodo del Gradiente Coniugato), http://it.wikipedia.org/wiki/Metodo del gradiente coniugato. Wikipedia (Metodo di Jacobi), http://it.wikipedia.org/wiki/Metodo di Jacobi. Wikipedia (Successive Over Relaxation), http://it.wikipedia.org/wiki/Successive Over Relaxation.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  81/ 81