Pivot de Gauss en Python

[lesurveilleur]

Bonsoir,

Je ne comprends pas bien la méthode derrière ce pivot de Gauss.

import numpy as np

def lu_decomposition_pivot(A):
    n = A.shape[0]
    L = np.zeros((n, n))
    U = np.zeros((n, n))
    P = np.identity(n)
    
    
    for i in range(n):
        # Pivotage partiel
        max_index = np.argmax(np.abs(A[i:, i])) + i
        if A[max_index, i] == 0:
            raise ValueError("La matrice est singulière")

        # Échange des lignes
        A[[i, max_index]] = A[[max_index, i]]
        P[[i, max_index]] = P[[max_index, i]]
        if i > 0:
            L[[i, max_index], :i] = L[[max_index, i], :i]

        L[i, i] = 1
        for j in range(i, n):
            U[i, j] = A[i, j] - np.sum(L[i, :i] * U[:i, j])
            if i != j:
                L[j, i] = (A[j, i] - np.sum(L[j, :i] * U[:i, i])) / U[i, i]

    return P, L, U
     

# Exemple d'utilisation
A = np.array([[2, 0, 2, 0.6],
              [3, 3, 4, -0.1],
              [1, 1, 1, 2],
              [1, 3, 3, 1]])

P, L, U = lu_decomposition_pivot(A)

print("Matrice P :")
print(P)

print("Matrice L :")
print(L)

print("Matrice U :")
print(U)

print("Vérification PA = LU :")
print(np.dot(P, A) - np.dot(L, U))

Je comprends mathématiquement cette méthode, mais pas numériquement.

Dans l'attente de vous lire
-----

[Biname]

Salut,
Mon assistant 4o trouve ceci :
------------------------------
Il semble que le problème vient de la manière dont la matrice A est modifiée directement lors du processus de décomposition, ce qui perturbe les calculs ultérieurs des matrices L et U. En modifiant A directement, les éléments de L et U calculés se basent sur une matrice qui change au fil des itérations, causant des incohérences.

Pour mettre en évidence l'erreur, regardons les points spécifiques où la matrice A est modifiée :

Échange des lignes dans A :

A[[i, max_index]] = A[[max_index, i]]

Cet échange affecte la matrice A et par conséquent, peut fausser les calculs de L et U.

Calcul des éléments de U :

U[i, j] = A[i, j] - np.sum(L[i, :i] * U[:i, j])

À ce stade, A a déjà été modifiée par les échanges de lignes, ce qui peut entraîner des erreurs dans les calculs de U.

Calcul des éléments de L :

    L[j, i] = (A[j, i] - np.sum(L[j, :i] * U[:i, i])) / U[i, i]

De même, les modifications précédentes de A peuvent induire des erreurs dans les calculs de L.

La solution consiste à utiliser une copie de la matrice A pour effectuer les opérations de pivotage et les calculs de U et L. Cela permet de préserver l'intégrité de la matrice originale et d'éviter les incohérences.

Voici le correctif appliqué :

 Cliquez pour afficher

import numpy as np

def lu_decomposition_pivot(A):
    n = A.shape[0]
    L = np.zeros((n, n))
    U = np.zeros((n, n))
    P = np.identity(n)
    A_copy = A.copy()  # Travail sur une copie de A pour éviter les modifications en place

    for i in range(n):
        # Pivotage partiel
        max_index = np.argmax(np.abs(A_copy[i:, i])) + i
        if A_copy[max_index, i] == 0:
            raise ValueError("La matrice est singulière")

        # Échange des lignes dans A_copy et P
        A_copy[[i, max_index]] = A_copy[[max_index, i]]
        P[[i, max_index]] = P[[max_index, i]]
        if i > 0:
            L[[i, max_index], :i] = L[[max_index, i], :i]

        # Définir la diagonale de L à 1
        L[i, i] = 1

        # Calcul des éléments de U
        for j in range(i, n):
            U[i, j] = A_copy[i, j] - np.sum(L[i, :i] * U[:i, j])

        # Calcul des éléments de L
        for j in range(i + 1, n):
            L[j, i] = (A_copy[j, i] - np.sum(L[j, :i] * U[:i, i])) / U[i, i]

    return P, L, U

# Exemple d'utilisation
A = np.array([[2, 0, 2, 0.6],
              [3, 3, 4, -0.1],
              [1, 1, 1, 2],
              [1, 3, 3, 1]])

P, L, U = lu_decomposition_pivot(A)

print("Matrice P :")
print(P)

print("Matrice L :")
print(L)

print("Matrice U :")
print(U)

print("Vérification PA = LU :")
print(np.dot(P, A) - np.dot(L, U))

En utilisant A_copy pour les opérations de décomposition, nous évitons les problèmes liés aux modifications en place de A, assurant ainsi la précision des calculs de L et U.
------------------------------
Est-ce que ceci répond à ta question ?