Informatica 5 Liceo Scientifico Scienze Applicate/Jacobi Risoluzione Sistemi Lineari

Wikibooks, manuali e libri di testo liberi.
Indice del libro

Risoluzione di un sistema di equazioni lineari con Jacobi[modifica]

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 puo' 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 piu' di una, vengono costruite in modo tale da garantire proprieta' particolari (funzione positiva, omogeneita', 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


la convergenza allora c'e' se norma


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

Altra cosa interessante da sapere e' la velocita' 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'e con il metodo di Jacobi e altre volte garantiscono una velocita' di convergenza maggiore di Jacobi. Non tutti i sistemi di equazioni lineari convergono con questi metodi iterativi.
In generale se la matrice e' 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 piu' 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