Risolvere equazioni lineari utilizzando le matrici e Python

Mostreremo un esempio di risoluzione di equazioni lineari Ax = B usando Python, in particolare quando la dimensione di A rende computazionalmente costoso calcolare l'inversa.

Risoluzione di equazioni lineari utilizzando matrici e Python

In un precedente articolo, abbiamo esaminato la risoluzione di un problema LP, cioè un sistema di equazioni lineari con vincoli di disuguaglianza. Se il nostro insieme di equazioni lineari ha vincoli deterministici, possiamo rappresentare il problema come matrici e applicare l’algebra matriciale. I metodi a matrice rappresentano più equazioni lineari in modo compatto mantenendo l’utilizzo delle funzioni della libreria di matrici esistenti.

Useremo NumPy (un buon tutorial qui) e SciPy (una guida di riferimento qui). Per l’installazione di questi moduli esistono centinaia di risorse nel web, ne indicherò solo una Installare lo stack SciPy.

Un esempio: sistema di equazioni lineari

Come nostra pratica, procederemo con un esempio, prima scrivendo il modello a matrice e poi usando Numpy per le soluzioni.

    \[\begin{array}{lcl} 2a+b+c & = & 4 \\ a+3b+2c & = & 5 \\ a&=&6 \end{array}\]

Ora possiamo formalizzare il problema con le matrici:

    \[A=\begin{bmatrix} 2 &1 &1 \\ 1 &3 &2 \\ 1 &0 &0 \end{bmatrix}, x=\begin{bmatrix} a\\ b\\ c \end{bmatrix}, B=\begin{bmatrix} 4\\ 5\\ 6 \end{bmatrix}\]

Quindi, le equazioni lineari possono essere scritte in questo modo:

    \[Ax=B\]

per risolvere per il vettore x, dobbiamo prendere l’inversa della matrice A e l’equazione viene scritta come segue:

    \[x=A^{-1}B\]

Utilizzo di numpy per risolvere il sistema

import numpy as np
# definiamo la matrice A usando Numpy
A = np.array([[2, 1, 1],
 [1, 3, 2],
 [1, 0, 0]]) 

#definiamo la matrice B
B = np.array([4, 5, 6]) 

# linalg.solve è la funzione di NumPy per risolvere un sistema di equazioni lineari
print "Soluzioni:\n",np.linalg.solve(A, B ) 

Soluzioni:
[  6.  15. -23.]

Dunque le soluzioni del sistema di equazioni lineari sono: a=6, b=15, c=-23

Quando le matrici crescono

All’aumentare del numero di variabili, aumenta anche la dimensione della matrice A e diventa computazionalmente costoso ottenere l’inversa della matrice A. Tra i vari metodi, considereremo una tra le procedure più utilizzate per ottenere la matrice A fattorizzata in matrici più semplici: la decomposizione LU.

La decomposizione LU

La decomposizione LU, nota anche come fattorizzazione inferiore superiore, è uno dei metodi per risolvere sistemi quadrati di equazioni lineari. Come suggerisce il nome, la fattorizzazione LU scompone la matrice A nel prodotto di due matrici: una matrice triangolare inferiore L e una matrice triangolare superiore U. La decomposizione può essere rappresentata come segue:

    \[LU=A\]

    \[\begin{bmatrix} l_{11} &0 &0 \\ l_{21} &l_{22} &0 \\ l_{31} &l_{32} &l_{33} \end{bmatrix}\begin{bmatrix} u_{11} &u_{12} &u_{13} \\ 0 &u_{22} &u_{23} \\ 0 &0 &u_{33} \end{bmatrix} =\begin{bmatrix} a_{11} &a_{12} &a_{13} \\ a_{21} &a_{22} &a_{23} \\ a_{31} &a_{32} &a_{33} \end{bmatrix}\]

    \[ \begin{bmatrix} l_{11}u_{11} &l_{11}u_{12} &l_{11}u_{13} \\ l_{21}u_{11} &l_{21}u_{12}+l_{22}u_{22} &l_{21}u_{13}+l_{22}u_{23} \\ l_{31}u_{11} &l_{31}u_{12}+l_{32}u_{22} &l_{31}u_{13}+l_{32}u_{23}+l_{33}u_{33} \end{bmatrix} =\begin{bmatrix} a_{11} &a_{12} &a_{13} \\ a_{21} &a_{22} &a_{23} \\ a_{31} &a_{32} &a_{33} \end{bmatrix} \]

La decomposizione LU con SciPy

#importo il pacchetto linalg di SciPy per la decomposizione LU 
import scipy.linalg as linalg 
#import NumPy
import numpy as np
#definisco A come prima
A = np.array([[2., 1., 1.],
 [1., 3., 2.],
 [1., 0., 0.]]) 

#definisco B
B = np.array([4., 5., 6.]) 

#invochiamo la funzione di fattorizzazione LU
LU = linalg.lu_factor(A) 

#risolviamo date LU and B
x = linalg.lu_solve(LU, B) 
print ("Soluzioni:\n",x)

#vediamo come è stata fattorizzata A, dove P è la cosiddetta matrice di permutazione
P, L, U = linalg.lu(A) 
print("P=\n",P,"\n","U=\n",U,"\n","L=\n",L)
print (np.dot(L,U))

Si ottiene:


[evaluate linalg.py]
Soluzioni:
 [  6.  15. -23.]
P=
 [[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]] 
 U=
 [[ 2.   1.   1. ]
 [ 0.   2.5  1.5]
 [ 0.   0.  -0.2]] 
 L=
 [[ 1.   0.   0. ]
 [ 0.5  1.   0. ]
 [ 0.5 -0.2  1. ]]
[[ 2.00000000e+00  1.00000000e+00  1.00000000e+00]
 [ 1.00000000e+00  3.00000000e+00  2.00000000e+00]
 [ 1.00000000e+00 -2.77555756e-17  2.77555756e-17]]


Date L e U, possiamo rappresentare la matrice fattorizata A come:

    \[\begin{bmatrix} 1 &0 &0 \\ 0.5 &1 &0 \\ 0.5 &-0.2 &1 \end{bmatrix}\begin{bmatrix} 2 &1 &1 \\ 0 &2.5 &1.5 \\ 0 &0 &-0.2 \end{bmatrix} =A=\begin{bmatrix} 2 &1 &1 \\ 1 &3 &2 \\ 1 &0 &0 \end{bmatrix}\]

Informazioni su Carlo Bazzo 9 Articoli
Sysadmin & network eng. @Epysoft, editor @TheTechGoggler, CTO @HDEMO. Ex Developer @MSFT, @GOOG. Ex Marketing Ops @XRX @HPQ. LinkedIn: it.linkedin.com/in/carlobazzo