Élimination gaussienne à l'aide du pivotement en Python

L’élimination gaussienne est également connue sous le nom de méthode de réduction de ligne. C’est un algorithme couramment utilisé pour résoudre des problèmes linéaires. L’algorithme implique une série d’opérations de ligne sur une matrice de coefficients tirés des équations linéaires jusqu’à ce que la matrice soit réduite à la forme échelonnée.

Les trois opérations suivantes sont effectuées lors de la transformation d’une matrice en une forme échelonnée. Ceux-ci inclus:

• Multiplier une ligne par un scalaire et l’ajouter ou la soustraire d’une autre ligne.
• Échange de lignes.
• Multiplication de lignes avec un scalaire.

Une matrice est sous forme d’échelon de ligne si le premier élément de chaque ligne, également appelé entrée principale, est un zéro ; l’entrée principale de chaque ligne est une colonne à droite de l’entrée principale de la ligne précédente. De plus, les lignes avec des éléments nuls doivent être en dessous des éléments ayant des lignes non nulles.

Les opérations de ligne sont largement classées en deux. L’élimination directe consiste à réduire une matrice à la forme échelonnée pour déterminer si une solution viable existe et si elle est finie ou infinie. D’autre part, la rétrosubstitution réduit davantage la matrice en une forme d’échelon réduit en ligne.

Élimination gaussienne avec pivotement en Python

Le pivotement est l’échange de lignes et de colonnes pour obtenir l’élément pivot approprié. Un élément pivot approprié doit être à la fois différent de zéro et significativement grand mais plus petit par rapport aux autres entrées de ligne.

Le pivotement est classé en pivotement partiel et en pivotement complet. Dans l’algorithme de pivotement partiel, l’élément le plus grand est considéré comme l’élément pivot pour minimiser les erreurs d’arrondi.

D’autre part, le pivotement complet comprend l’échange de lignes et de colonnes pour obtenir le meilleur élément de pivot, augmentant ainsi la précision.

Un ensemble d’équations est considéré comme linéaire si aucune variable n’a plus d’un exposant. L’élimination gaussienne implique une série d’étapes ; la première étape consiste à créer une matrice de coefficients.

Une matrice de coefficients est simplement une matrice de coefficients tirés d’un ensemble d’équations linéaires. L’étape suivante consiste à créer une matrice augmentée qui est ensuite soumise à une série d’opérations qui la réduisent sous forme échelonnée.

Cependant, dans les cas où l’élément pivot est nul ou de très faible amplitude, nous devons échanger la rangée pivot avec une rangée inférieure.

Lors de la mise en œuvre de la méthode d’élimination gaussienne, nous devons noter les indices de tableau. Les itérables Python tels que les listes et les tableaux commencent souvent à l’index 0 et se terminent à l’index n-1.

On peut alors lire le contenu de la matrice augmentée, appliquer la méthode d’élimination, rétrosubstituer, et enfin afficher la réponse.

from numpy import array, zeros, fabs, linalg

a = array(
    [
        [0, 6, -1, 2, 2],
        [0, 3, 4, 1, 7],
        [5, 1, 0, 3, -1],
        [3, 1, 3, 0, 2],
        [4, 4, 1, -2, 1],
    ],
    float,
)
# the b matrix constant terms of the equations
b = array([5, 7, 2, 3, 4], float)

print("Solution by NumPy:")


print(linalg.solve(a, b))

n = len(b)
x = zeros(n, float)

# first loop specifys the fixed row
for k in range(n - 1):
    if fabs(a[k, k]) < 1.0e-12:

        for i in range(k + 1, n):
            if fabs(a[i, k]) > fabs(a[k, k]):
                a[[k, i]] = a[[i, k]]
                b[[k, i]] = b[[i, k]]
                break

    # applies the elimination below the fixed row

    for i in range(k + 1, n):
        if a[i, k] == 0:
            continue

        factor = a[k, k] / a[i, k]
        for j in range(k, n):
            a[i, j] = a[k, j] - a[i, j] * factor
            # we also calculate the b vector of each row
        b[i] = b[k] - b[i] * factor
print(a)
print(b)


x[n - 1] = b[n - 1] / a[n - 1, n - 1]
for i in range(n - 2, -1, -1):
    sum_ax = 0

    for j in range(i + 1, n):
        sum_ax += a[i, j] * x[j]

    x[i] = (b[i] - sum_ax) / a[i, i]

print("The solution of the system is:")
print(x)

Production :

Solution by NumPy:
[ 0.38947368  0.49473684 -0.10877193  0.12982456  0.83157895]
[[ 5.00000000e+00  1.00000000e+00  0.00000000e+00  3.00000000e+00  -1.00000000e+00]
 [ 0.00000000e+00  3.00000000e+00  4.00000000e+00  1.00000000e+00  7.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  4.50000000e+00  0.00000000e+00  6.00000000e+00]
 [ 0.00000000e+00  4.44089210e-16  0.00000000e+00  3.52702703e+00  2.95945946e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00  2.11354647e+00]]
[2.         7.         4.5        2.91891892 1.75758075]
The solution of the system is:
[ 0.38947368  0.49473684 -0.10877193  0.12982456  0.83157895]

Dans le code ci-dessus, nous avons utilisé la fonction array et la fonction fabs fournies par la bibliothèque NumPy pour créer une matrice et lire des valeurs absolues. Nous avons aussi utilisé Linalg ; une sous-bibliothèque NumPy utilisée pour effectuer des opérations telles que le calcul des valeurs propres, des vecteurs et des déterminants.

from numpy import array, zeros, fabs, linalg

La variable a représente les coefficients constants des variables, tandis que b représente la matrice des termes constants correspondant aux équations. Généralement, la matrice et le vecteur sont tirés de la représentation matricielle de quatre équations linéaires que l’on se propose de résoudre.

a = array(
    [
        [0, 6, -1, 2, 2],
        [0, 3, 4, 1, 7],
        [5, 1, 0, 3, -1],
        [3, 1, 3, 0, 2],
        [4, 4, 1, -2, 1],
    ],
    float,
)
# the b matrix constant terms of the equations
b = array([5, 7, 2, 3, 4], float)

Le but de l’élimination est de transformer la matrice pour avoir des zéros dans la diagonale principale.

Lors de l’exécution de l’élimination gaussienne, la première ligne doit rester fixe. Cependant, si l’élément pivot de la première ligne est égal à zéro, nous devons échanger la première ligne avec n’importe quelle ligne suivante sans zéro comme élément principal.

Nous avons utilisé la boucle k pour spécifier la ligne fixe dans le code ci-dessus. Dans ce cas, k représente l’indice de la ligne fixe qui va de 1 à n-1.

Les instructions if suivantes vérifient si l’élément pivot est égal à zéro et l’échangent avec la ligne appropriée. Le i utilisé représente l’indice des lignes en dessous de la ligne fixe.

for k in range(n - 1):
    if fabs(a[k, k]) < 1.0e-12:

        for i in range(k + 1, n):
            if fabs(a[i, k]) > fabs(a[k, k]):
                a[[k, i]] = a[[i, k]]
                b[[k, i]] = b[[i, k]]
                break

Dans le code ci-dessous, nous utilisons une boucle pour nous assurer que l’élimination n’est pas appliquée dans la première ligne.

Nous avons ensuite multiplié chaque ligne en dessous de la ligne fixe par un facteur. Un facteur est égal à l’élément diagonal correspondant à la ligne fixe divisé par le premier élément de la ligne correspondante.

Ensuite, nous soustrayons les deux lignes de la ligne fixe pour obtenir l’élément de la deuxième colonne sous la diagonale principale. Les lignes suivantes peuvent également être réduites dans la même procédure.

Comme nous l’avons mentionné, la première ligne restera toujours inchangée. Nous avons également calculé le facteur et le vecteur b pour chaque ligne, comme indiqué ci-dessous.

for i in range(k+1, n):
    # first row should remain intact
    if a[i, k] == 0:
        continue
    # calculating the factor
    factor = a[k, k]/a[i, k]
    for j in range(k, n):
        a[i, j] = a[k, j] - a[i, j]*factor
        # calculating the b vector
    b[i] = b[k] - b[i]*factor
print(a)
print(b)

Le terme rétrosubstitution fait référence au fait que les valeurs de x peuvent être calculées de la dernière équation à la première équation. Dans ce cas, nous exploitons les valeurs correspondantes du vecteur scalaire b multipliées par les coefficients de x. Le résultat est ensuite divisé par l’élément correspondant dans la diagonale principale.

Nous avons priorisé le calcul des valeurs de x à partir de n et les valeurs restantes de x de n-1 à la première. Par contre, la sommation devrait commencer à partir de j+1 puisque seules les valeurs de x nous sont connues des calculs précédents. De plus, la valeur de la somme est le résultat des valeurs de x multipliées par les valeurs de a. Nous pouvons maintenant calculer la valeur de x puisque nous avons la valeur de la somme.

# calculating the value of x beginning from n-1
x[n - 1] = b[n - 1] / a[n - 1, n - 1]
for i in range(n - 2, -1, -1):
    sum_ax = 0
    # calculating the value of the sum
    for j in range(i + 1, n):
        sum_ax += a[i, j] * x[j]
    # calculating the value of x
    x[i] = (b[i] - sum_ax) / a[i, i]

print(x)

Production :

[[ 2.00000000e+00  3.00000000e+00  4.00000000e+00  1.00000000e+00
   7.00000000e+00]
 [ 0.00000000e+00  7.00000000e+00 -1.00000000e+00  3.00000000e+00
   1.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00 -1.30000000e+01  2.00000000e+00
  -2.10000000e+01]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  2.81250000e+00
   5.81250000e+00]
 [ 0.00000000e+00  8.88178420e-16  0.00000000e+00 -4.44089210e-16
   4.95734797e+00]]
[  7.           5.         -14.           0.625        0.15371622]
The solution of the system is:
[0.02170543 0.79224806 1.05116279 0.15813953 0.03100775]

La fonction NumPy suivante devrait également aboutir à la même réponse.

print("The numpy solution:")
print(linalg.solve(a, b))

Production :

Solution by Nuumpy:
[0.02170543 0.79224806 1.05116279 0.15813953 0.03100775]