Equazioni lineari con Python: metodo di eliminazione di Gauss

E' il metodo più familiare per la soluzione dei sistemi di equazioni lineari. Esso si compone di due fasi: la fase di eliminazione e quella di sostituzione all'indietro.

Risolvere i sistemi lineari con Python: metodo di eliminazione di Gauss

Le equazioni lineari e algebriche si possono trovare in quasi tutti i rami dell’analisi numerica. Ma la loro applicazione più naturale nell’ingegneria è nell’analisi dei sistemi lineari. L’ampia classe dei sistemi lineari include strutture, solidi elastici, flussi di calore, campi elettromagnetici, circuiti elettrici e molto altro.

La modellazione di sistemi lineari dà origine a equazioni lineari della forma Ax = b, dove b è il vettore degli ingressi e il vettore x rappresenta la risposta del sistema. La matrice di coefficienti A, che riflette le caratteristiche intrinseche del sistema, è indipendente dall’ingresso. Ovvero, se modifichiamo l’ingresso, il sistema di equazioni lineari da risolvere avrà un vettore b diverso ma la stessa matrice dei coefficienti A.

Metodi di risoluzione dei sistemi di equazioni lineari

Oltre ai metodi iterativi, dei quali non parleremo in questa sede, esistono i metodi cosiddetti diretti. La loro caratteristica comune è che essi cercano di trasformare le equazioni originali in un sistema equivalente di più facile soluzione.

Tale trasformazione può essere ottenuta applicando 3 operazioni fondamentali:

  1. scambiare due righe di A (il determinante di A cambia segno);
  2. moltiplicare una riga di A per uno scalare λ non nullo (il determinante di A viene moltiplicato per lo stesso scalare);
  3. sostituire una riga di A con quella che si ottiene sommando tale riga con un’altra riga moltiplicata per uno scalare (lascia invariato il determinante di A);

Naturalmente queste operazioni non influenzano le soluzioni del sistema che rimangono le stesse, ma possono influenzare i coefficienti di A ed il suo determinante.

La tabella seguente, riassume i 3 principali metodi diretti di risoluzione:

Metodo Forma iniziale Forma finale
Eliminazione di Gauss Ax=b Ux=c
Decomposizione LU Ax=b LUx=b
Eliminazione di Gauss-Jordan Ax=b Ix=c

In tabella, nella colonna della forma finale delle equazioni, appaiono 3 matrici standard, che sono:

U=\begin{bmatrix} U_{11}&U_{12}&U_{13} \\ 0&U_{22}&U_{23} \\ 0&0&U_{33} \\\end{bmatrix} detta triangolare superiore;

L=\begin{bmatrix} L_{11}&0&0 \\ L_{21}&L_{22}&0 \\ L_{31}&L_{32}&L_{33} \\\end{bmatrix} detta triangolare inferiore;

I=\begin{bmatrix} 1&0&0 \\ 0&1&0 \\ 0&0&1 \\\end{bmatrix} detta matrice identità.

Il motivo per il quale i metodi diretti di risoluzione cercano di trasformare la matrice A in una matrice triangolare, risulta subito evidente quando si scriva, ad esempio Ux=c, che significa risolvere il sistema di equazioni lineari seguente:

\left\{\begin{matrix} U_{11}x_{1}+U_{12}x_{2}+U_{13}x_{3}=c_{1}\\ U_{22}x_{2}+U_{23}x_{3}=c_{2}\\U_{33}x_{3}=c_{3}\end{matrix}\right.

come si vede esso risulta di facile risoluzione partendo dall’ultima equazione.

Il metodo di eliminazione di Gauss

E’ il metodo più familiare per la soluzione dei sistemi di equazioni lineari. Esso si compone di due fasi: la fase di eliminazione e quella di sostituzione all’indietro. La prima fase ha lo scopo, come indicato nella precedente tabella, di trasformare le equazioni dalla forma Ax=b a quella di immediata soluzione Ux=c.

La fase di eliminazione procede eseguendo ripetutamente la 3^ operazione fondamentale dell’elenco precedentemente discusso. Senza perderci nelle complicazioni della notazione per il caso generale, procederemo con un esempio.

Risolveremo il sistema di equazioni lineari seguente:

\left\{\begin{matrix} x+y+z=1\\x-y-z=1 \\x-2y+3z=-5\end{matrix}\right.

Soluzione manuale con il metodo di Gauss

Prima di vedere l’algoritmo codificato in Python, procederemo manualmente con il metodo di eliminazione di Gauss.

La matrice A del sistema è:

A=\begin{bmatrix} 1&1 &1 \\ 1&-1 &-1 \\ 1&-2 &3 \\\end{bmatrix}

e il vettore dei coefficienti b è:

b=\begin{bmatrix} 1\\1 \\-5\end{bmatrix}

Come detto, lo scopo della prima fase del metodo di Gauss è quello di trasformare la matrice A in una matrice traingolare superiore U. Come si vede da A, il primo passo sarà quello di rendere nulli i termini a21 e a31. Per fare questo dobbiamo individuare una costante λ tale che moltiplicata per la riga R1 e sommata alla riga R2 e successivamente R3, rendano nulli i coefficienti desiderati.

1. Per annullare a21 moltiplicheremo R1 per -1 e poi la sommeremo a R2:

-1\cdot R_{1}|b_{1}+R_{2}|b_{2}=\begin{bmatrix} 0&-2 &-2 &|0\\\end{bmatrix}

Da notare che la stessa operazione viene anche eseguita sui corrispondenti elementi del vettore b.

Analogamente per l’annullamento del coefficiente a31:

-1\cdot R_{1}|b_{1}+R_{3}|b_{3}=\begin{bmatrix} 0&-3 &2 &|-6\\\end{bmatrix}

Sostituendo le nuove righe trasformate nella matrice A, abbiamo:

A=\begin{bmatrix} 1&1 &1 \\ 0&-2 &-2 \\ 0&-3 &2 \\\end{bmatrix}

e il vettore b trasformato adesso è:

b=\begin{bmatrix} 1\\0 \\-6\end{bmatrix}

2. Procediamo adesso ad annullare il coefficiente a32.

-\frac{3}{2}\cdot R_{2}|b_{2}+R_{3}|b_{3}=\begin{bmatrix} 0&0 &5 &|-6\\\end{bmatrix}

Sostituendo la nuova riga così ottenuta:

A=\begin{bmatrix} 1&1 &1 \\ 0&-2 &-2 \\ 0&0 &5 \\\end{bmatrix}

e il vettore b, che resta immutato:

b=\begin{bmatrix} 1\\0 \\-6\end{bmatrix}

Il procedimento di eliminazione è concluso, come si può vedere la matrice A è stata trasformata in una matrice triangolare superiore U e il vettore b è pure trasformato in un nuovo vettore che nella tabella abbiamo chiamato c. A questo punto, partendo dal basso, è agevole trovare le soluzioni del sistema lineare:

5z=-6; z=\frac{-6}{5}
-2y-2z=0; y=-z=\frac{6}{5}
x=1-y-z=1

Le soluzioni allora sono:

x=1; y=\frac{6}{5}; z=-\frac{6}{5}

Algoritmo di eliminazione di Gauss in Python

Vi sono due variazioni rispetto al procedimento manuale: la prima è che le righe vengono trasformate per sottrazione invece che per somma; la seconda è che le righe trasformate non vengono sostituite alle righe originali della matrice A, solo gli elementi propri della matrice triangolare superiore vengono sostituiti. Infatti, gli elementi non appartenenti a U (matrice trasformata) sono irrilevanti per il calcolo delle soluzioni.


import numpy as np
def gaussElim(a,b):
    n = len(b)
    # Fase di eliminazione
    for k in range(0,n-1):
        for i in range(k+1,n):
            if a[i,k] != 0.0:
                #se elemento non nullo calcolo lambda opportuno
                lam = a [i,k]/a[k,k]
                #calcolo la nuova riga della matrice
                a[i,k+1:n] = a[i,k+1:n] - lam*a[k,k+1:n]
                #aggiorno il vettore b
                b[i] = b[i] - lam*b[k]
                # sostituzione all'indietro
    for k in range(n-1,-1,-1):
        b[k] = (b[k] - np.dot(a[k,k+1:n],b[k+1:n]))/a[k,k]
    
    return b

#inserisco i coefficienti iniziali
a=np.array([[1.0,1.0,1.0],[1.0,-1.0,-1.0],[1.0,-2.0,3.0]])
b=np.array([1.0,1.0,-5.0])
aOrig = a.copy() # Salvo matrice originale
bOrig = b.copy() #salvo vettore b iniziale
x = gaussElim(a,b)
#stampo la matrice A per verifica
print(a)
print("x =\n",x)
#det = np.prod(np.diagonal(a)) #eventuale calcolo determinante
#print("\ndet =",det)
#eventuale calcolo precisione del risultato
print("\nCheck risultato: [a]{x} - b =\n",np.dot(aOrig,x) - bOrig)

Esecuzione:


Python 3.8.8 
[evaluate ex-gauss.py]
A=
 [[ 1.  1.  1.]
 [ 1. -2. -2.]
 [ 1. -3.  5.]]
x =
 [ 1.   1.2 -1.2]

Check risultato: [a]{x} - b =
 [2.22044605e-16 0.00000000e+00 0.00000000e+00]


Come si diceva, la matrice A non appare triangolarizzata ma solo perchè gli elementi al di sotto della diagonale non vengono sostituiti effettivamente essendo irrilevanti per il calcolo delle soluzioni. Se guardiamo invece gli elementi diversi da zero di U:

\begin{Bmatrix} a_{11}&a_{12} &a_{13} &a_{22} &a_{23} &a_{33} \\\end{Bmatrix}

essi corrispondono esattamente a quelli ottenuti nel nostro procedimento manuale.

servizialleimprese
Avatar photo
Informazioni su Carlo Bazzo 16 Articoli
Sysadmin & network eng. @Epysoft, editor @TheTechGoggler, CTO @HDEMO. Former developer @MSFT @GOOG. Former MOps consultant @XRX @HPQ. LinkedIn: it.linkedin.com/in/carlobazzo