Risoluzione sistemi

# Codice Python per risoluzione di sistemi
# (non ottimizzato)

import numpy as np

# Input: matrice completa A|B
# Output: matrice ridotta
def solveGauss(m,verbose=True):
    n=m.copy()
    rows, cols=n.shape
    print("-------------------")
    print("Matrice incompleta:")
    print( m[:,range(cols-1)] )
    print("Termini noti:")
    print( m[:,cols-1] )
    print("\nRiduzione")
    # itera su tutte le posizioni
    for i in range(rows):
        # Riduzione righe
        k=np.nonzero(n[i,:])[1]
        if len(k)==0:
            continue
        k=k[0]
        n[i,:]/=n[i,k]
        for j in range(rows):
            if i==j:
                continue
            n[j,:]-=n[i,:]*n[j,k]
        if verbose:
            print(n,"\n")
    return n

# Stampa soluzione in forma polinomiale
def getSol(n):
    rows, cols = n.shape
    fix = []
    sol=np.zeros(cols-1)
    gen=np.matrix(np.zeros((cols-1,cols)))
    for i in range(rows):
        k=np.nonzero(n[i,:])[1]
        if len(k)==0:
            continue
        if k[0]==cols-1:
            print("Sistema non compatibile")
            return False
        fix.append(k[0])
        print ("x(",k[0]+1,")=",n[i,cols-1],sep='',end='')
        sol[k[0]]=n[i,cols-1]
        for j in k:
            if j>k[0] and j<cols-1:
                print("+(",-n[i,j],")*x(",j+1,")",sep='',end='')
                gen[k[0],j]=-n[i,j]
        print("")
    for i in range(cols-1):
        if not i in fix:
            print("x(",i+1,")=x(",i+1,")",sep='')
            gen[i,i]=1.
    # restituisce solo le colonne non nulle (e traspone)
    gen = [ x for x in np.array(gen.transpose()) if x.any() ]
    print("")
    return np.matrix(sol), np.matrix(gen)

def checkSol(m,s,gen):
    rows, cols=m.shape
    print("Matrice incompleta:")
    print( m[:,range(cols-1)] )
    print("Termini noti:")
    print( m[:,cols-1] )
    print("Soluzione particolare: ")
    print(s)
    print("Generatori soluzioni sistema omogeneo associato:")
    print(gen)
    print(" --- verifica AS^t-B=0 ---")
    res=m[:,range(cols-1)]*s.transpose()-m[:,cols-1]
    print(res)
    print(" --- verifica AZ^t=0 --- ")
    res2=m[:,range(cols-1)]*gen.transpose()
    print(res2)
    return res,res2
# Esempi

AB=np.matrix('1,2,3,0,5,6,2; 2,4,6,0,10,12,4; 0,3,2,1,0,0,-1; 1,2,5,3,2,1,2; 0,0,0,1,2,0,1')*1.

R=solveGauss(AB)
s,gen=getSol(R)
print(checkSol(AB,s,gen))

AB2=np.matrix(np.random.randint(0,high=10,size=(6,10)))*1.
R2=solveGauss(AB2)
s2,gen2=getSol(R2)
print(checkSol(AB2,s2,gen2))