Obiettivo: Vogliamo calcolare, utilizzando Python, tutte le soluzioni per un sistema lineare di primo grado utilizzando le matrici e il metodo di Cramer. Per farlo non utilizzeremo librerie aggiuntive, se non quanto disponibile nativamente dentro Python.
Ovviamente lo scopo dell’esercizio è l’utilizzo di Python, visto che si potrebbe risolvere agilmente con librerie come SciPy e NumPy.
Supponiamo di avere un sistema lineare come quello seguente (l’ho scopiazzato spudoratamente da YouMath, dove c’è anche un’esaustiva spiegazione sul metodo di Cramer):
La matrice associata ai coefficienti ed il relativo determinante sarebbe questo:
Mentre la matrice dei termini noti sarebbe:
Giusto per curiosità la soluzione immediata con NumPy sarebbe la seguente:
|
1 2 3 4 5 6 |
import numpy as np a = np.array([[2,1,1],[4,-1,1],[-1,1,2]]) b = np.array([1,-5,5]) x = np.linalg.solve(a, b) print x |
Il risultato nel caso specifico darà: [-1. 2. 1.]
Noi vogliamo costruire qualcosa di simile in Python 2.7.
Cominciamo creando la base della nostra classe perché possa accettare i medesimi argomenti:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
class SistemaLineare(): __coeff = [] __tnoti = [] def __init__(s,coeff,tnoti): s.__coeff = coeff s.__tnoti = tnoti def solve(s): return [] sl = SistemaLineare([[2,1,1],[4,-1,1],[-1,1,2]],[1,-5,5]) x = sl.solve() print x |
Fatto questo cominciamo a creare i metodi di cui avremo bisogno. Il metodo di Creamer con il calcolo del determinante implica, per matrici di dimensione superiore a 2×2, che si debba poter ridurre la matrice in sottomatrici.
Giusto per chiarezza, nel caso del determinante D della suddetta matrice, il calcolo si svolgerebbe in questo modo:
Ogni sottomatrice prende tutti i valori della matrice corrente eccetto quelli nella riga e nella colonna della posizione attuale. Scriviamo quindi la nostra funzione nella maniera seguente:
|
1 2 3 4 5 6 7 8 9 10 11 12 |
def __subMatrice(s,r,c,matrice): nuova = [] nr = len(matrice) nc = len(matrice[0]) for i in range(nr): if i == r: continue riga = [] for j in range(nc): if j == c: continue riga.append(matrice[i][j]) nuova.append(riga) return nuova |
Per calcolare il determinante vogliamo fare due operazioni distinte: il calcolo per matrici 2×2 e il calcolo per riduzione delle matrici di dimensione superiore. Iniziamo scrivendo l’operazione per quelle 2×2:
|
1 2 |
def __determ2x2(s,matrice): return matrice[0][0]*matrice[1][1]-matrice[0][1]*matrice[1][0] |
Per calcolare gli altri determinanti procederemo invece con:
|
1 2 3 4 5 6 7 8 9 10 |
def __determ(s,matrice): nr = len(matrice) nc = len(matrice[0]) if nr == 2 and nc == 2: return s.__determ2x2(matrice) else: risultato = 0 for c in range(nc): risultato += (1-2*(c%2))*matrice[0][c]*s.__determ(s.__subMatrice(0,c,matrice)) return risultato |
In particolare con (1-2*(c%2)) calcoliamo il segno per ogni elemento della prima riga. In particolare ricordiamo che gli elementi saranno calcolati in ordine da 0, 1, 2, 3, 4 ecc. Se calcoliamo il resto per ognuno di essi avremo 0, 1, 0, 1, 0 ecc. Eseguendo la suddetta operazione avremo quindi 1-0, 1-2, 1-0, 1-2, 1-0 ecc. Questo ci darà i segni come 1, -1, 1, -1, 1 ecc.
A questo punto ci manca solo un metodo che ci permetta di creare nuove matrici per calcolare i vari determinanti per x, y, z ecc. Ricordiamoci che i vari risultati saranno dati da:
Dobbiamo fare attenzione ad un dettaglio importante. In Python le liste vengono passate per riferimento, questo vuol dire che modificando ogni lista “copiata” si modifica in realtà la lista originale all’indirizzo di memoria. Per copiare una matrice multidimensionale è necessario quindi copiare i singoli valori nella nuova matrice, nella maniera seguente darebbe errore:
|
1 |
tmp = [c for c in [r for r in matrice]] |
Quindi procediamo in questo modo:
|
1 2 3 4 5 6 7 |
nr = len(matrice) nc = len(matrice[0]) tmp = [] for i in range(nr): tmp.append([]) for j in range(nc): tmp[i].append(matrice[i][j]) |
La nostra funzione di sostituzione sostituirà quindi i termini noti nella colonna preselezionata.
|
1 2 3 4 5 6 |
def __sostitMatr(matrice,tnoti,c): nr = len(matrice) tmp = [c for c in [r for r in matrice]] for r in range(nr): tmp[r][c] = tnoti[r] return tmp |
A questo punto possiamo scrivere il metodo che calcolerà la soluzione, alla maniera seguente:
|
1 2 3 4 5 6 7 8 9 10 11 |
def __soluzione(s,matrice,tnoti): nc = len(matrice[0]) d = s.__determ(matrice) if d != 0: x = [] for c in range(nc): dx = s.__determ(s.__sostitMatr(matrice,tnoti,c)) x.append(1.0*dx/d) return x else: return [] |
Infine implementiamo il metodo solve aggiungendo un paio di chiamate di errore, nel caso in cui la matrice non sia adatta.
|
1 2 3 4 5 6 7 8 9 |
def solve(s): nr = len(s.__coeff) if nr == 0: raise Exception("Matrice con 0 righe") else: nc = len(s.__coeff[0]) if nc != nr: raise Exception("Matrice non quadrata") else: return s.__soluzione(s.__coeff,s.__tnoti) return [] |
Nel suo complesso la classe che avremo creato sarà la seguente:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 |
class SistemaLineare(): __coeff = [] __tnoti = [] def __init__(s,coeff,tnoti): s.__coeff = coeff s.__tnoti = tnoti def solve(s): nr = len(s.__coeff) if nr == 0: raise Exception("Matrice con 0 righe") else: nc = len(s.__coeff[0]) if nc != nr: raise Exception("Matrice non quadrata") else: return s.__soluzione(s.__coeff,s.__tnoti) return [] def __subMatrice(s,r,c,matrice): nuova = [] nr = len(matrice) nc = len(matrice[0]) for i in range(nr): if i == r: continue riga = [] for j in range(nc): if j == c: continue riga.append(matrice[i][j]) nuova.append(riga) return nuova def __determ2x2(s,matrice): return matrice[0][0]*matrice[1][1]-matrice[0][1]*matrice[1][0] def __determ(s,matrice): nr = len(matrice) nc = len(matrice[0]) if nr == 2 and nc == 2: return s.__determ2x2(matrice) else: risultato = 0 for c in range(nc): risultato += (1-2*(c%2))*matrice[0][c]*s.__determ(s.__subMatrice(0,c,matrice)) return risultato def __sostitMatr(s,matrice,tnoti,c): nr = len(matrice) nc = len(matrice[0]) tmp = [] for i in range(nr): tmp.append([]) for j in range(nc): tmp[i].append(matrice[i][j]) for r in range(nr): tmp[r][c] = tnoti[r] return tmp def __soluzione(s,matrice,tnoti): nc = len(matrice[0]) d = s.__determ(matrice) if d != 0: x = [] for c in range(nc): dx = s.__determ(s.__sostitMatr(matrice,tnoti,c)) x.append(1.0*dx/d) return x else: return [] |
Per utilizzarla sarà sufficiente scrivere:
|
1 2 3 |
sl = SistemaLineare([[2,1,1],[4,-1,1],[-1,1,2]],[1,-5,5]) x = sl.solve() print x |
Se abbiamo fatto tutto correttamente il risultato che otterremo sarà: [-1.0, 2.0, 1.0]