Jacobi Risoluzione Sistemi Lineari (superiori)

Risoluzione di un sistema di equazioni lineari con Jacobi

Siamo partiti a risolvere un sistema lineare con il metodo della sostituzione adatto per piccoli sistemi e calcolato manualmente, siamo passati a risolverlo numericamente con Cramer e poi abbiamo visto una soluzione tramite Gauss. In tutti i casi precedenti abbiamo trovato l'esatta soluzione del problema. Ora ci spostiamo in una nuova dimensione delle tecniche risolutive e illustriamo il metodo di Carl Gustav Jacob Jacobi (Potsdam, 10 dicembre 1804 – Berlino, 18 febbraio 1851) . Pensiamo di essere nel 1825 e prendiamo carta e penna

Scriviamo il nostro problema:

5x1-x2+2x3=12
3x1+8x2-2x3=25
x1+x2+4x3=6

Ora isoliamo la x1 nella prima equazione, la x2 nella seconda e la x3 nella terza (facciamo finta per ora di non aver notato che le incognite isolate si trovano lungo la diagonale)

si ha :

x1 =          +1/5 x2 -2/5 x3 +12/5
x2 = -3/8 x1          +2/8 x3 +25/8
x3 = -1/4 x1  -1/4 x2         +6/4

ho spaziato per pensare successivamente in termini matriciali, ma la sostanza non cambia

ora ipotizzo che x1= 1 x2=1 x3=1 ( i valori attributi sono casuali) , ora se li sostituisco nel secondo membro delle equazioni, ottengo dei nuovi valori di x1 x2 x3
nel nostro caso ottengo

x1=2.2
x2=3
x3=1

(fate pure i calcoli a mano , io ho usato Octave con la formula

N=[0 1/5 -2/5 ; -3/8 0 2/8; -1/4 -1/4 0]
b=[12/5; 25/8; 6/4]
x=[1; 1; 1]
x=N*x+b

)

ora utilizzo questi nuovi valori per sotituirli nel secondo membro delle equazioni e ricavo delle nuove x1 x2 e x3 i cui valori sono

x1=2.60
x2=2.55
x3=0.20

( ho rieseguito la x=N*x+b )

ripeto la procedura iterativamente

iterazione3 x1=2.83 x2=2.20 x3=0.21250

iterazione4 x1=2.75500 x2=2.11687 x3=0.24250

iterazione5 x1=2.72638 x2=2.15250 x3=0.28203

iterazione6 x1=2.71769 x2=2.17312 x3=0.28028

iterazione7 x1=2.72251 x2=2.17594 x3=0.27730

iterazione8 x1=2.72427 x2=2.17338 x3=0.27539

iterazione 9 x1=2.72452 x2=2.17225 x3=0.27559

iterazione10 x1=2.72421 x2=2.17220 x3=0.27581

iterazione11 x1=2.72412 x2=2.17237 x3=0.27590

iterazione12 x1=2.72412 x2=2.17243 x3=0.27588

iterazione13 x1=2.72413 x2=2.17243 x3=0.27586

iterazione14 x1=2.72414 x2=2.17242 x3=0.27586

iterazione15 x1=2.72414 x2=2.17241 x3=0.27586

iterazione16 x1=2.72414 x2=2.17241 x3=0.27586

guardando i valori mi accorgo che c'e' stata una convergenza sui valori assunti dalle incognite , incredibile ma queste sono le soluzioni del sistema di equazioni iniziali , calcolandole con Cramer si ha infatti

A=[5 -1 2; 3 8 -2; 1 1 4]
b=[12; 25; 6]
x=inv(A)*b

x1= 2.27414
x2= 2.17241
x3= 0.27586

Ora Jacobi come metodo per risolvere il problema propone di scomporre la matrice A dei coefficienti delle incognite nella matrice D (che deve essere invertibile) che contiene solo i coefficienti sulla diagonale e nella matrice N in cui ci sono tutti i coefficienti della A invertiti di segno esclusi quelli della diagonale , cioe'

A = D - N

allora il sistema può essere scritto come

Ax=b
(D-N)x=b
Dx -Nx=b
Dx= Nx +b
x = inv(D)Nx + inv(D)b

ora pensiamo di calcolare la x a primo membro come soluzione del sistema al secondo membro , come se ci fossero due vettori x uno attuale e uno futuro, quello attuale lo indichiamo come x(k) e quello futuro come x(k+1) ora partendo da un generico x(0) ( i cui valori li decidiamo noi arbitrariamente ) possiamo ricavare x(1) poi da questo x(2) etc iterativamente (cioe' in modo ciclico)

ci fermiamo nel calcolo quando x(k) e' vicino a x(k+1) ammesso che ci sia la convergenza, naturalmente i matematici che sono precisi il vicino lo traducono con una espressione matematica che prende il nome di norma ( di norme ce ne sono più di una , vengono costruite in modo tale da garantire proprieta' particolari (funzione positiva, omogeneità, rispetto diseguaglianza triangolare) . Ora una norma (che prende il nome di norma euclidea) per stabilire quanto diversi/uguali sono due vettori x e y e' la seguente

$\left\|\mathbf {x} -y\right\|:={\sqrt {\sum _{i=1}^{n}(x_{i}-y_{i})^{2}}}$

la convergenza allora c'e' se norma $\left\|\mathbf {x} (k+1)-x(k)\right\|<{\text{quantità piccola che fissiamo noi}}$

Interessante data una matrice A è sapere prima di procedere con jacobi se ci sarà convergenza, ci sono alcune dimostrazioni che dicono c'è convergenza se la matrice A è una matrice diagonale dominante per righe, un'altra dimostrazione dice c'e' convergenza se l'autovalore max di inv(D)*N è minore in modulo di 1 (raggio spettrale minore di uno).

Altra cosa interessante da sapere e' la velocità di convergenza del metodo, cioe' quante iterazioni dobbiamo fare. Oltre alla tecnica di Jacobi ce ne sono altre similari (Gauss-Siedel ,SOR etc) che partono da una scomposizione della matrice A diversa da quella operata da Jacobi e che assicurano la convergenza in casi in cui non c'è con il metodo di Jacobi e altre volte garantiscono una velocità di convergenza maggiore di Jacobi. Non tutti i sistemi di equazioni lineari convergono con questi metodi iterativi.
In generale se la matrice è densa e di dimensione contenuta si usano i metodi diretti (Cramer, eliminazione di Gauss) che consentono la risoluzione esatta del sistema di equazioni e non presentano problemi di convergenza e terminano in un numero finito di passi; se la matrice e' sparsa o di grandi dimensioni si usano quelli indiretti (Jacobi, Gauss-Siedel, del Sovrarilassamento SOR etc.) che sono più veloci computazionalmente ma che possono presentare problemi di convergenza.

Ecco il codice che risolve iterativamente il nostro problema, le iterazioni si ripetono finche' la norma euclidea non diventa minore di 0.001 diag(diag(A)) genera una matrice nulla tranne sulla diagonale dove compaiono gli elementi della diagonale di A

A=[5 -1 2; 3 8 -2; 1 1 4]
b= [12; 25;6]
D=diag(diag(A));
N=-(A-D);

x=[1;1;1]
do
 y= inv(D)*N*x + inv(D)*b;
 norma=norm(y-x);
 x=y;
until ( norma<0.001)
 x