Simultaneous Equations

RISOLUZIONE DEI SISTEMI LINEARI

I sistemi lineari, se possiedono una soluzione determinata, possono essere risolti da *metodi diretti* o da metodi iterativi. I metodi diretti permettono di trovare la soluzione esatta con un numero di operazioni finito e predicibile che dipende dall'ordine del sistema n. I metodi iterativi producono invece soluzioni approssimate ma con un numero di operazioni che è approssimativamente indipendente dall'ordine del sistema. Per piccoli sistemi (n<100) i metodi diretti sono più veloci mentre per sistemi di ordine più elevato è preferibile usare metodi iterativi.
METODI DIRETTI
Con il metodo di Cramer si risolve direttamente un sistema di n equazioni a n incognite mediante il calcolo di n+1 determinanti di ordine n. Tale metodo tuttavia è impraticabile perchè il numero delle operazioni richieste è elevato. I metodi di eliminazione invece offrono un numero relativamente basso di operazioni. Essi sono basati nel produrre combinazioni lineari delle equazioni in modo opportuno in modo da pervenire ad una equazione con una sola incognita. Poi, ripercorrendo i passi percorsi all'indietro, si trovano via via tutte le altre incognite . Se si scrive il sistema di equazioni in forma matriciale si può osservare che, quando si finiscono le eliminazioni, si deve ottenere una matrice di coefficienti di tipo triangolare . I sistemi di questo tipo si chiamano sistemi lineari triangolari.
RISOLUZIONE DEI SISTEMI LINEARI TRIANGOLARI. I sistemi lineari del tipo: $\begin{pmatrix} a_{11} & a_{12} & …& …& a_{1n}\\ 0& a_{22} &⋯&⋯&a_{2n}\\ 0& 0 &a_{33}&⋯&a_{3n}\\\,&\,&⋱&⋱&⋮\\\,& \, &\,&0&a_{nn}\end{pmatrix}\begin{pmatrix} x_1 \\ x_2\\x_3\\ \\⋮\\x_n \end{pmatrix}=\begin{pmatrix} b_1 \\ b_2\\b_3\\ \\⋮\\b_n \end{pmatrix}\qquad oppure \qquad\begin{pmatrix} \\ a_{11}&0&0&&0\\a_{21}& a_{22} &0&&0\\a_{31}& a_{32} &a_{33}&⋯&0\\ ⋮& ⋮ & ⋮&⋱&⋮\\a_{n1} & a_{n2} &a_{n3}& …& a_{nn}\end{pmatrix}\begin{pmatrix} x_1 \\ x_2\\x_3\\ \\⋮\\x_n \\\end{pmatrix}=\begin{pmatrix} b_1 \\ b_2\\b_3\\ \\⋮\\b_n \end{pmatrix}$ sono detti sistemi triangolari. Il primo è triangolare superiore ed è risolvibile con il metodo di sostituzione all'indietro (back substitution), il secondo è triangolare inferiore ed è risolvibile con il metodo di sostituzione in avanti (forward substitution). In ambedue i casi i sistemi sono risolvibili solo se tutti gli elementi diagonali della matrice dei coefficienti sono diversi d un a 0. È possibile ricavare delle espressioni generali per le sostituzioni ALGORITMO PER LA FORWARD SUBSTITUTION. $\,\quad\qquad\qquad \begin{cases} x_1=\frac{b_1}{a_11} \\ x_i=\frac{1}{a_{ii}}\left(b_i-{\sum_{1}^{i-1}}_ja_{ij}x_j\right) \qquad i=2,3,⋯,n \end{cases}$ ALGORITMO PER LA BACK SUBSTITUTION . $\;\qquad\qquad\qquad \begin{cases} x_n=\frac{b_n}{a_{nn}} \\ x_i=\frac{1}{a_{ii}}\left(b_i-{\sum_{i+1}^{n}}_ja_{ij}x_j\right) \qquad i=n-1,n-2, ⋯,1 \end{cases}$ Con le espressioni degli algoritmi è possibile scrivere dei segmenti di codice che li implementano. ESEMPIO: CODICE IN FORTRAN PER LA BACK SUBSTITUTION: x(nd)=b(nd)/a(nd,nd) backsubstitution: DO i=nd-1,1, -1 sum=0. DO j=i+1,nd sum=sum+a(i,j)*x(j) END DO x(i)=(b(i)-sum)/a(i,i) END DO backsubstitution
ELIMINAZIONE DI GAUSS. Il metodo di eliminazione successiva delle incognite è detto *metodo di Gauss. Il metodo di eliminazione di Gauss consiste quindi nel produrre ripetutamente combinazioni lineari delle equazioni in modo da trasformare il sistema assegnato in un sistema a matrice triangolare superiore. Considderiamo un sistema di n equazioni a n incognite. Il primo passo consiste nell'eliminare gli elementi della prima colonna al di sotto della diagonale principale. Per questo fine: sottraiamo la prima equazione moltiplicata con m₁₂= a₂₁/a₁₁ dalla seconda equazione : $\begin{split} a_21x_1+a_22x_2+a_23x_3+⋯+a_{2n}x_n=b_2\qquad - \\\frac{a_21}{a_11}\left(a_11x_1+a_12x_2+a_13x_3+⋯+a_{1n}x_n\right)=\frac{a_21}{a_11}b_1\end{split}$ ________________________________________________________ $\left(a_22-\frac{a_21}{a_11}a_12\right)x_2+\left(a_23-\frac{a_21}{a_11}a_13\right)x_3+⋯+\left(a_{2n}-\frac{a_21}{a_11}a_{1n}\right)x_n=b_2-\frac{a_21}{a_11}b_1$ Abbiamo ottenuto un'equazione in cui x₁ è stata eliminata. ora sottraiamo la prima equazione moltiplicata con m₁₃= a₃₁/a₁₁ dalla terza equazione: $\begin{split} a_31x_1+a_32x_2+a_33x_3+⋯+a_{3n}x_n=b_3\qquad - \\\frac{a_31}{a_11}\left(a_11x_1+a_12x_2+a_13x_3+⋯+a_{1n}x_n\right)=\frac{a_31}{a_11}b_1\end{split}$ ________________________________________________________ $\left(a_32-\frac{a_31}{a_11}a_12\right)x_2+\left(a_23-\frac{a_31}{a_11}a_13\right)x_3+⋯+\left(a_{3n}-\frac{a_21}{a_11}a_{1n}\right)x_n=b_3-\frac{a_31}{a_11}b_1$ Ecco una seconda equazione dove manca x₁. ........ ........ sottraiamo la prima equazione moltiplicata con m_1n=a_n1/a₁₁ dalla n-esima equazione. Ecco una n-1-esima equazione dove manca x₁. Come risultato si avrà la nuova matrice con elementi nulli sulla prima colonna, dalla seconda riga in poi. Il sistema diventa: $\begin{pmatrix}a_11& a_{12} & …& …& a_{1n} \\ 0&a_{22}^{(1)} &⋯&⋯&a_{2n}^{(1)}\\0& a_23^{(1)} &a_{33}^{(1)}&⋯&a_{3n}^{(1)}\\\,&\,&⋮&⋱&⋮\\\,& \, &\,&a_{n-1n}^{(1)}&a_{nn}^{(1)}\end{pmatrix}\begin{pmatrix} x_1\\x_2\\x_3\\ \\⋮\\x_n \end{pmatrix}=\begin{pmatrix} b_1\\b_2^{(1)}\\b_3^{(1)}\\ \\⋮\\b_n^{(1)} \end{pmatrix}$ Ora consideriamo il sistema ridotto (n-1)x(n-1) che si ottiene scartando la prima riga e la prima colonna: $\begin{pmatrix} \\ a_{22}^{(1)} &⋯&⋯&a_{2n}^{(1)}\\ a_23^{(1)} &a_{33}^{(1)}&⋯&a_{3n}^{(1)}\\⋮&⋮&⋱\\\,& \, &a_{n-1n}^{(1)}&a_{nn}^{(1)}\end{pmatrix}\begin{pmatrix} x_2\\x_3\\ \\⋮\\x_n \end{pmatrix}=\begin{pmatrix} b_2^{(1)}\\b_3^{(1)}\\ \\⋮\\b_n^{(1)} \end{pmatrix}$ A questa sottomatrice applichiamo lo stesso procedimento di prima, sottraendo la prima equazione della sottomatrice moltiplicata con $m_23=a_32^{(1)}/a_22^{(1)}$ dalla seconda equazione della sottomatrice, e così via come prima fino ad eliminare la n-1 esima incognita. Dopo tutti questi passi, rimettendo insieme le prime righe di tutti i sottosistemi che si formano, si ottiene il sistema triangolare superiore cercato : $\begin{pmatrix}a_11& a_{12} & …& …& a_{1n} \\ 0&a_{22}^{(1)} &⋯&⋯&a_{2n}^{(1)}\\0& 0 &a_{33}^{(2)}&⋯&a_{3n}^{(2)}\\⋮&⋮&⋮&⋱&⋮\\0 &0 &0&0&a_{nn}^{(n-1)}\end{pmatrix}\begin{pmatrix} x_1\\x_2\\x_3\\ \\⋮\\x_n \end{pmatrix}=\begin{pmatrix} b_1\\b_2^{(1)}\\b_3^{(2)}\\ \\⋮\\b_n^{(n-1)} \end{pmatrix}$ CODICE IN FORTRAN PER L'ELIMINAZIONE DI GAUSS: eliminavariabili: DO i=1,nd-1 righe: DO j=i+1,nd factor=a(j,i)/a(i,i) colonne: DO k=1,nd a(j,k)=a(j,k)-factora(i,k) END DO colonne b(j)=b(j)-factor*b(i) END DO righe END DO eliminavariabili
L'eliminazione di Gauss è un algoritmo che può essere implementato solo se gli elementi sulla diagonale principale sono diversi da zeri, altrimenti nel calcolo dei coefficienti accade di dividere con zero. Se si incontra questo problema si può cerca scambiare righe sottostanti alla riga presa in considerazione in modo da evitare elementi diagonali nulli. Questa operazione di scambio è nota come tecnica di pivoting. Esempio: $\begin{cases} x_1+2x_2+3x_3=6 \\ 2x_1+4x_2-x_3=5\\x_1-2x_2+2x_3=1 \end{cases}⇒\begin{cases} x_1+2x_2+3x_3=6 \\ 2x_1+4x_2-x_3=5→m_12=2/1\\x_1-2x_2+2x_3=1→m_13=1/1 \end{cases}⇒\begin{cases} x_1+2x_2+3x_3=6 \\ 0x_2-7x_3=-7\\-4x_2-x_3=-5→m_23=-4/0 \end{cases}$ Sembra che il metodo di eliminazione non funzioni. In realtà il sistema ha soluzione se scambiamo la seconda e la terza riga: $\begin{cases} x_1+2x_2+3x_3=6 \\ -4x_2-x_3=-5\\-7x_3=-7 \end{cases}$ Anzi si ottiene subito il sistema triangolare.
LA TECNICA DEL PIVOTING. In generale, se si ottiene un elemento pivotale k-esimo nullo, a_kk = 0, con la tecnica del pivoting si cerca una riga i sotto la riga k-esima che abbia un elemento non nullo nella colonna k e poi si scambiano queste due righe e se la matrice dei coefficienti è invertibile questa riga sicuramente c'è. La scelta è sempre tra le equazioni che seguono la k-esima. Tuttavia la strategia di permutare equazioni/righe viene usata anche se si incontrano elementi pivot diversi da zero per aumentare la stabilità dell’eliminazione gaussiana. Si può infatti dimostrare che per aumentare la stabilità è opportuno evitare moltiplicatori troppo grandi, poiché amplificano gli errori di arrotondamento, con conseguente instabilità dell’algoritmo. Per prevenire la crescita dei moltiplicatori è necessario prevenire elementi pivot troppo piccoli (perchè poi, nei moltiplicatori, sono un quoziente). Per raggiungere questo obiettivo possiamo applicare la strategia del *pivoting parziale* ammettendo scambi di righe anche se al passo k l’elemento pivot a_kk non è nullo. Più precisamente, al passo k si determina l’elemento di modulo massimo fra tutti gli elementi a_rk , con r ≥ k e, una volta individuato l’indice r corrispondente all’elemento di modulo massimo, si scambiano fra loro la riga k-esima e r -esima. Se l’elemento più grande in modulo si cerca tanto nelle colonne quanto nelle righe allora il metodo prende il nome di pivoting totale. Un codice che implementi l'eliminazione di Gauss deve contenere tre moduli: eliminazione in avanti pivoting sostituzione all’indietro. Ma prima si deve monitorizzare il termine diagonale durante la fase di pivoting per rilevare la presenza di uno zero (o di un numero minore di una tolleranza predefinita) e quindi avvertire che il sistema è singolare. Se ciò avviene la ricerca delle soluzioni termina (det=0)
METODO DI GAUSS-JORDAN. Il metodo di Gauss-Jordan è molto simile all'eliminazione di Gauss per le matrici quadrate, ed è un metodo che tramite operazioni elementari (quali lo scambio di righe, combinazioni lineari, ecc...) opera la seguente trasformazione della matrice completa associata : $A_c=\begin{pmatrix} a_11 & a_12&⋯&a_{1i}&⋯&a_{1n}&b_1 \\ a_21 & a_22&⋯&a_{2i}&⋯&a_{2n}&b_2\\ ⋯ & ⋯&⋱&⋯&⋯&⋯&⋯ \\a_{i1} & a_{i2}&⋯&a_{ii}&⋯&a_{in}&b_i\\⋯ & ⋯&⋯&⋯&⋱&⋯&⋯\\a_{n1} & a_{n2}&⋯&a_{ni}&⋯&a_{nn}&b_n\end{pmatrix} ⇒A'_c=\begin{pmatrix} 1 & 0&⋯&0&⋯&0&b'_1 \\ 0 & 1&⋯&0&⋯&0&b'_2\\ ⋯ & ⋯&⋱&⋯&⋯&⋯&⋯ \\0 & 0&⋯&1&⋯&0&b'_i\\⋯ & ⋯&⋯&⋯&⋱&⋯&⋯\\0 & 0&⋯&0&⋯&1&b'_m\end{pmatrix}$ Una volta eseguita la trasformazione, la soluzione del sistema è rappresentata dal vettore colonna dei termini noti. La trasformazione della matrice è basata sul fatto che se si applicano a un sistema lineare una o più delle tre operazioni elementari seguenti: Scambiare fra loro due equazioni. Moltiplicare ambo i membri di una equazione per uno stesso numero diverso da zero. Sostituire a una equazione la sua somma con una qualsiasi combinazione lineare di alcune altre. (Nel caso più semplice: sommare a una equazione un’altra equazione moltiplicata per un numero) si ottiene un sistema equivalente al precedente, (cioè che ammette tutte e sole le soluzioni di quello di partenza). Un algoritmo che, passo per passo, dalla prima riga alla n-esima, trasforma la matrice incompleta in matrice unitaria è : $\begin{cases} a^k_{kj}=\frac{a^{k-1}_{kj}}{a^{k-1}_{kk}}\quad \;\; \qquad\qquad\quad se \; i=k \\ a^k_{ij}=a^{k-1}_{ij}-a^{k-1}_{ij}\cdota^k_{kj}\quad se \; i≠k \end{cases}$ Qui k non indica un esponente ma l'indice di ogni passo del processo di trasformazione. $\begin{cases} k=1,⋯,n \qquad (servono\; n\; passi) \\ i=1,⋯,n\qquad (per\;tutte\;le\;n\;righe)\\ j=k,⋯,n+1\qquad (per\;k,...,n+1\;colonne) \end{cases}$ Esempio. Consideriamo il sistema lineare: $\begin{aligned} 4x_1+x_2+2x_3 &= 16 \\ x_1+3x_2+x_3 &= 10 \\x_1+2x_2+5x_3&=12\end{aligned}$ La matrice completa è: $A_c=\begin{pmatrix} 4 & 1&2&16 \\ 1 & 3&1&10\\1&2&5&12 \end{pmatrix}$ e k=0. Si inizia con k=1 (primo passo), i=1,2,3 e j=1,2,3,4 $A^1_c=\begin{pmatrix} a^1_11=\frac{a^0_11}{a^0_11} & a^1_12=\frac{a^0_12}{a^0_11} & a^1_13=\frac{a^0_13}{a^0_11}& a^1_14=\frac{a^0_14}{a^0_11}\\ a^1_21=a^0_21-a^0_21a^1_11 \quad& a^1_22=a^0_22-a^0_21a^1_12\quad &a^1_23=a^0_23-a^0_21a^1_13 \quad&a^1_24=a^0_24-a^0_21a^1_14\\ a^1_31=a^0_31-a^0_31a^1_11 \quad& a^1_32=a^0_32-a^0_31a^1_12\quad &a^1_33=a^0_33-a^0_31a^1_13 \quad&a^1_34=a^0_34-a^0_31a^1_14 \end{pmatrix}⇒A^1_c=\begin{pmatrix} 1 & \frac{1}{4}&\frac{1}{2}&4 \\ 0&\frac{11}{4}&\frac{1}{2}& 6\\0&\frac{7}{4}&\frac{9}{2}&8 \end{pmatrix}$ Secondo passo: k=2, i=1,2,3, j=2,3,4 $A^2_c=\begin{pmatrix} a^1_11=\frac{a^0_11}{a^0_11} & a^2_12=a^1_12-a^1_12a^2_22\quad &a^2_13=a^1_13-a^1_12a^2_23 \quad&a^2_14=a^1_14-a^1_12a^2_24 \\ a^1_21=a^0_21-a^0_21a^1_11 \quad& a^2_22=\frac{a^1_22}{a^1_22}\quad &a^2_23=\frac{a^1_23}{a^1_22} \quad&a^2_24=\frac{a^1_24}{a^1_22}\\ a^1_31=a^0_31-a^0_31a^1_31 \quad& a^2_32=a^1_32-a^1_32a^2_22\quad &a^2_33=a^1_33-a^1_32a^2_23 \quad&a^2_34=a^1_34-a^1_32a^2_24 \end{pmatrix}⇒ A^2_c=\begin{pmatrix} 1 &0&\frac{5}{11}&\frac{38}{11}\\ 0&1&\frac{2}{11}& \frac{24}{11}\\0&0&\frac{46}{11}&\frac{46}{11} \end{pmatrix}$ Terzo passo: k=3, i=1,2,3, j=3,4 $A^3_c=\begin{pmatrix} a^1_11=\frac{a^0_11}{a^0_11} & a^2_12=a^1_12-a^1_12a^2_22 & a^3_13=a^2_13-a^2_13a^3_33 \quad&a^3_14=a^2_14-a^2_13a^3_34\\ a^1_21=a^0_21-a^0_21a^1_11 \quad& a^2_22=\frac{a^1_22}{a^1_22}\quad &a^3_23=a^2_23-a^2_23a^3_33 \quad&a^3_24=a^2_24-a^2_23a^3_34\\ a^1_31=a^0_31-a^0_31a^1_31 \quad& a^2_32=a^1_32-a^1_32a^2_22\quad &a^3_33=\frac{a^2_33}{a^2_33} \quad&a^3_34=\frac{a^2_34}{a^2_33} \end{pmatrix}⇒ A^3_c=\begin{pmatrix} 1 & 0&0&3\\ 0&1&0& 2\\0&0&1&1\end{pmatri$ Il metodo di Gauss –Jordan e quello di Gauss possono sembrare simili, ma quello di Gauss-Jordan richiede un pò di più tempo di calcolo. Si può dimostrare che il numero di operazioni con Gauss-Jordan è dato da $\frac{n^3}{2}+O(n^2)$ mentre con Gauss è $\frac{n^3}{3}+O(n^2)$ che aumenta con n ma Gauss-Jordan ne ha il 33% in più. Tuttavia il metodo di Gauss-Jordan ha il vantaggio di una maggiore semplicità di programmazione.
SISTEMI TRIDIAGONALI. Per ottenere le soluzioni numeriche delle equazioni differenziali ordinarie e alle derivate parziali e anche in certi metodi di interpolazione (es: spline) si ha la necessità di risolvere il sistema lineare della forma generale: $\begin{matrix} a_11x_1 & +a_12x_2 & ⋯ & ⋯ & ⋯ & ⋯ & ⋯& =a_{1,n+1}\\ a_21x_1 & +a_22x_2 & +a_23x_3 & ⋯ & ⋯ & ⋯ & ⋯& =a_{2,n+1} \\⋯&a_32x_2 & +a_33x_3 & +a_34x_4 & ⋯ & ⋯ & ⋯& =a_{3,n+1}\\&&&⋮&&&\\&&a_{i,i-1}x_{i-1}&+a_{ii}x_i&+a_{i,i+1}x_{i+1}&&&=a_{i,n+1}\\&&&&&&⋮\\\\&&&&&a_{n,n-1}x_{n-1} & +a_{nn}x_n &=a_{n,n+1}\\ \end{matrix}$ nel quale in ogni equazione tutti i coefficienti sono nulli eccetto quelli nella diagonale principale e nelle diagonali adiacenti. Questi sistemi si dicono tridiagonali e possono essere risolti mediante algoritmi di eliminazione. Una versione semplificata dell'eliminazione di Gauss per i sistemi tridiagonali è l'algoritmo di Thomas .
METODI ITERATIVI
I metodi iterativi non sono solo usati quando il sistema è di ordine elevato ma anche quando si ha un sistema con un buon numero di coefficienti uguali a zero, come, per esempio, nella matrice completa: $\begin{pmatrix} 1 & 3&0&0&2 \\0 & 2&1&0&5 \\1 & 0&3&4&-2 \\6 & 0&0&5&4\end{pmatrix}$ (queste matrici si dicono "sparse"). Un metodo iterativo, tenendo conto dall'inizio degli zeri, può dare risultati più rapidi rispetto ad un metodo diretto. In un metodo iterativo si trasforma il sistema iniziale: $\begin{cases} a_11x_1+a_12x_2+⋯+a_{1n}x_n=b_1 \\ a_12x_1+a_22x_2+⋯+a_{2n}x_n=b_2\\⋯\quad⋯\quad⋯\quad⋯\quad⋯\quad⋯\quad⋯\quad⋯\\a_n1x_1+a_{n2}x_2+⋯+a_{nn}x_n=b_n \end{cases}⇒\begin{pmatrix} a_{11} & a_{12} & …& …& a_{1n}\\ a_21& a_{22} &⋯&⋯&a_{2n}\\ a_31& a_32 &a_{33}&⋯&a_{3n}\\⋯&⋯&⋯&⋱&⋯\\a_{n1}& a_{n2} &a_{n3}&⋯&a_{nn}\end{pmatrix}\begin{pmatrix} x_1 \\ x_2\\x_3\\ \\⋮\\x_n \end{pmatrix}=\begin{pmatrix} b_1 \\ b_2\\b_3\\ \\⋮\\b_n \end{pmatrix}⇒A\cdotx=b$ nella forma equivalente: $\widevec{x}^{k+1}=\widevec{h}+C\cdot\widevec{x}^k$ con $\widevec{x}$ e $\widevec{h}$ vettore colonna nx1 e C una matrice nxn. Si inizia l'iterazione fissando ad arbitrio un valore per l'approssimazione iniziale x⁰, prendendo, per esempio, per approssimazione iniziale la colonna dei termini noti: $\widevec{x}^0=\widevec{h}=\widevec{b}$ Successivamente costruiscono le matrici colonne: $\begin{matrix} \widevec{x}^1=\widevec{h}+C\;\widevec{x}^0 \qquad (prima\;\;approssimazione)\\ \widevec{x}^2=\widevec{h}+C\;\widevec{x}^1 \qquad (seconda\;\;approssimazione)\end{matrix}$ e così via, l'approssimazione si calcola in generale con la formula: $\widevec{x}^{i+1}=\widevec{h}+C\;\widevec{x}^{i}$ Se la successione delle approssimazioni $\widevec{x}^0, \widevec{x}^1, \widevec{x}^2 ,⋯⋯$ possiede un limite : $\widevec{x}=\underset{i⇒\infinity}{lim}\widevec{x}^i$ questo limite è la soluzione del sistema. Esistono diverse formulazioni nei vari metodi iterativi tutti che si differenziano nella scelta della matrice C. I più importanti sono il metodo di Jacobi e il metodo di Gauss-Siedel.
METODO DI JACOBI. Il metodo funziona se $a_{ii}\ne 0$ . Dato il sistema $A\cdot\widevec{x}=\widevec{b}$ , consideriamo le tre seguenti matrici: $D=\begin{pmatrix} a_11 & 0&⋯&0 \\ 0&a_22&⋯&0\\⋯&⋯&⋯ & ⋯\\0&0&⋯&a_{nn} \end{pmatrix}\qquad E=\begin{pmatrix} 0 & 0&0&⋯&0&0 \\ a_21&0&0&⋯&0&0\\a_31&a_32&0&⋯&0&0\\⋯&⋯&⋯&⋯&⋯&⋯\\a_{n-1,1}&a_{n-1,2}&a_{n-1,3}&⋯&0&0\\a_{n1}&a_{n2}&a_{n3}&&a_{n,n-1}&0\end{pmatrix}\qquad F=\begin{pmatrix} 0 & a_12&a_13&⋯&a_{1n} \\ 0& 0&a_23&⋯&a_{2n}\\⋯&⋯&⋯&⋯&⋯\\ 0&0&0&⋯&a_{n-1,n}\\0&0&0&⋯&0\end{pmatrix}$ La matrice D è una matrice diagonale avente, nell'ordine, sulla diagonale gli elementi della diagonale principale di, le matrici E ed F sono rispettivamente la triangolare inferiore e la triangolare superiore. Si ha evidentemente che : A= D + E + F. Sostituendo questa scomposizione di A nel sistema $A\cdot\widevec{x}=\widevec{b}$ si ottiene $\left(D+E+F\right)\cdot \widevec{x}=\widevec{b}$ da cui: $D\cdot\widevec{x}=\widevec{b}-\left(E+F\right)\cdot\widevec{x}$ ed essendo $a_{ii}\ne 0$ esiste l'inversa $D^{-1}$ , per cui moltiplicando a sinistra per $D^{-1}$ , ambo i membri si ottiene: $\widevec{x}=D^{-1}\widevec{b}-D^{-1}\left(E+F\right)\widevec{x}$ che fornisce il metodo iterativo di Jacobi nel quale ogni i+1 esima approssimazione si calcola con la formula: $\widevec{x}^{i+1}=D^{-1}\widevec{b}-D^{-1}\left(E+F\right)\widevec{x}^$
METODO DI GAUSS-SIEDEL. Con le stesse posizioni del metodo di Jacobi possiamo scrivere la $\left(D+E+F\right)\cdot \widevec{x}=\widevec{b}$ nella forma: $\left(D+E\right)\widevec{x}=\widevec{b}-F\widevec{x}$ Essendo $a_{ii}\ne 0$ la matrice D + E non è singolare e quindi si può ricavare direttamente il vettore $\widevec{x}$ : $\widevec{x}=\left(D+E\right)^{-1}\widevec{b}-\left(D+E\right)^{-1}\cdotF\widevec{x}$ che fornisce il metodo iterativo di Gauss-Siedel nel quale ogni i+1 esima approssimazione si calcola con la formula: $\widevec{x}^{i+1}=(D+E)^{-1}\widevec{b}-(D+E)^{-1}\cdotF\widevec{x}^i$