Pour tester

A = [[2, 3], [5, -2]]
B = [[5], [-16]]
I = [[1, 0], [0, 1]]

import numpy as np
print(np.linalg.solve(A, B))
print()
print(np.linalg.inv(A))

[[-2.]
 [ 3.]]

[[ 0.10526316  0.15789474]
 [ 0.26315789 -0.10526316]]

Outils

def dimensions(M):
    return len(M), len(M[0])

def matrice_nulle(lignes, colonnes):
    """ Surtout pas [[0]*p]*n """
    return [[0] * colonnes for _ in range(lignes)]

def copie_matrice(M):
    lignes, colonnes = dimensions(M)
    R = matrice_nulle(lignes, colonnes)
    for i in range(lignes):
        for j in range(colonnes):
            R[i][j] = M[i][j]
    return R

def identite(n):
    M = matrice_nulle(n, n)
    for i in range(n):
        M[i][i]=1
    return M

def echange_ligne(A, i, j):
    lignes, colonnes = dimensions(A)
    for k in range(colonnes):
        A[i][k], A[j][k] = A[j][k], A[i][k]

def transvection(A, i1, i2, mu):
    """ Li1 <- Li1 + mu * Li2 """
    _, colonnes = dimensions(A)
    for k in range(colonnes):
        A[i1][k] = A[i1][k] + mu *  A[i2][k]

def dilatation(A, i, coeff):
    """ Li <- coeff * Li """
    _, colonnes = dimensions(A)
    for k in range(colonnes):
        A[i][k] = A[i][k] * coeff

Pivot de Gauss

# Attention, il manque le pivot partiel. Ne fonctionne pas tout le temps !

def gauss_bof(A0, B0):
    """ Pivot de Gauss : on se ramène à la matrice identité
        A.X=B <=> I.X=solution  """
    A, B = copie_matrice(A0), copie_matrice(B0)
    lignes, colonnes = dimensions(A)
    for j in range(colonnes):
        for i in range(lignes):
            if i != j:
                mu = -A[i][j]/A[j][j]     # Ajj : pivot
                transvection(A, i, j, mu)
                transvection(B, i, j, mu)
    # ici, A est une matrice diagonale
    for i in range(lignes):
        coeff = 1/A[i][i]
        dilatation(A, i, coeff) # Pour avoir des 1 sur la diagonale
        dilatation(B, i, coeff)
    return A, B                  # A est la matrice identité ; B est la solution

print(gauss_bof(A, B))
print(gauss_bof(A, I))

([[1.0, 0.0], [-0.0, 1.0]], [[-2.0], [3.0]])
([[1.0, 0.0], [-0.0, 1.0]], [[0.10526315789473684, 0.15789473684210525], [0.2631578947368421, -0.10526315789473684]])

def pivot_partiel(M, j):
    lignes, colonnes = dimensions(M)
    p = j
    for i in range(j, lignes):
        if abs(M[i][j]) > abs(M[p][j]):
            p = i
    return p

def gauss_ok(A0, B0):
    """ Pivot de Gauss : on se ramène à la matrice identité
        A.X=B <=> I.X=solution  """
    A, B = copie_matrice(A0), copie_matrice(B0)
    lignes, colonnes = dimensions(A)
    for j in range(colonnes):
        p = pivot_partiel(A, j)
        echange_ligne(A, j, p)
        echange_ligne(B, j, p) # Avec la même ligne p
        for i in range(lignes):
            if i != j:
                mu = -A[i][j]/A[j][j]     # Ajj : pivot
                transvection(A, i, j, mu)
                transvection(B, i, j, mu)
    for i in range(lignes):
        coeff = 1/A[i][i]
        dilatation(A, i, coeff) # Pour avoir des 1 sur la diagonale
        dilatation(B, i, coeff)
    return A, B

T = [[1, 1, 1], [2, 2, 0], [1, 2, 3]]
gauss_ok(T, identite(3))

([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]],
 [[3.0, -0.5, -1.0], [-3.0, 1.0, 1.0], [1.0, -0.5, 0.0]])

np.linalg.solve(T, identite(3))

array([[ 3. , -0.5, -1. ],
       [-3. ,  1. ,  1. ],
       [ 1. , -0.5,  0. ]])

%xmode plain
gauss_bof(T, identite(3)) # Ne fonctionnait pas dans certains cas !

Exception reporting mode: Plain

Traceback (most recent call last):

  File "<ipython-input-48-2d5b9a1656a9>", line 2, in <module>
    gauss(T, identite(3))

  File "<ipython-input-32-f1b5582cddc6>", line 11, in gauss
    mu = -A[i][j]/A[j][j]     # Ajj : pivot

ZeroDivisionError: float division by zero

Autre méthode

Avec résolution d'équation à la fin... vraiment moins bien

$$x_i = \frac{b_i - \sum_{k=i+1}^{n}{a_{i,k}\cdot x_k}}{a_{i,i}}$$

def gauss_resolution(A0, B0):
    """ B : matrice colonne uniquement ! """
    A, B = copie_matrice(A0), copie_matrice(B0)
    lignes, colonnes = dimensions(A)
    for j in range(colonnes):
        p = pivot_partiel(A, j)
        echange_ligne(A, j, p)
        echange_ligne(B, j, p)
        for i in range(lignes):
            if i > j:           # Ici, seulement triangulaire supérieure
                mu = -A[i][j]/A[j][j]
                transvection(A, i, j, mu)
                transvection(B, i, j, mu)
    # A est triangulaire supérieure, il faut résoudre des équations
    print(A, B)
    X = matrice_nulle(lignes, 1)
    for i in range(lignes-1, -1, -1):
        # C'est compliqué hein
        X[i][0] = 1/A[i][i]*(B[i][0]-sum([A[i][k]*X[k][0] for k in range(i+1, colonnes)]))
    return X

gauss_resolution(A, B)

[[5, -2], [0.0, 3.8]] [[-16], [11.4]]

[[-2.0], [3.0]]

np.linalg.solve(A, B)

array([[-2.],
       [ 3.]])