Gaußsche Eliminierung mit Pivotierung in Python

Die Gaußsche Elimination ist auch als Zeilenreduktionsverfahren bekannt. Es ist ein Algorithmus, der häufig verwendet wird, um lineare Probleme zu lösen. Der Algorithmus umfasst eine Reihe von Zeilenoperationen an einer Matrix von Koeffizienten, die aus den linearen Gleichungen gezogen werden, bis die Matrix auf eine Stufenform reduziert wird.

Die folgenden drei Operationen werden durchgeführt, wenn eine Matrix in eine Stufenform transformiert wird. Diese schließen ein:

• Multiplizieren einer Zeile mit einem Skalar und Addieren oder Subtrahieren von einer anderen Zeile.
• Zeilen vertauschen.
• Multiplizieren von Zeilen mit einem Skalar.

Eine Matrix ist in Zeilenstufenform, wenn das erste Element in jeder Zeile, auch bekannt als führender Eintrag, eine Null ist; der führende Eintrag in jeder Zeile befindet sich eine Spalte rechts vom führenden Eintrag in der vorherigen Zeile. Darüber hinaus sollten Zeilen mit Null-Elementen unter Elementen mit Nicht-Null-Zeilen liegen.

Zeilenoperationen werden weitgehend in zwei Klassen eingeteilt. Bei der Vorwärtseliminierung wird eine Matrix in eine Stufenform reduziert, um zu bestimmen, ob eine brauchbare Lösung existiert und ob sie endlich oder unendlich ist. Andererseits reduziert die Rücksubstitution die Matrix weiter auf eine zeilenreduzierte Stufenform.

Gaußsche Eliminierung mit Pivotierung in Python

Pivotieren ist das Vertauschen von Zeilen und Spalten, um das passende Pivot-Element zu erhalten. Ein geeignetes Pivot-Element sollte sowohl ungleich Null als auch signifikant groß, aber kleiner im Vergleich zu den anderen Zeileneinträgen sein.

Schwenken wird in teilweises Schwenken und vollständiges Schwenken eingeteilt. Beim partiellen Pivot-Algorithmus wird das größte Element als Pivot-Element betrachtet, um Rundungsfehler zu minimieren.

Andererseits beinhaltet das vollständige Pivotieren das Vertauschen von Zeilen und Spalten, um das beste Pivot-Element zu erhalten, wodurch die Genauigkeit erhöht wird.

Ein Satz von Gleichungen wird als linear betrachtet, wenn keine Variable einen Exponenten von mehr als eins hat. Die Gaußsche Eliminierung umfasst eine Reihe von Schritten; Der erste Schritt besteht darin, eine Koeffizientenmatrix zu erstellen.

Eine Koeffizientenmatrix ist einfach eine Matrix von Koeffizienten, die aus einem Satz linearer Gleichungen gezogen werden. Der nächste Schritt besteht darin, eine erweiterte Matrix zu erstellen, die dann einer Reihe von Operationen unterzogen wird, die sie in eine Stufenform reduzieren.

In Fällen, in denen das Pivot-Element jedoch null oder sehr klein ist, müssen wir die Pivot-Zeile mit einer niedrigeren Zeile vertauschen.

Bei der Implementierung der Gaußschen Eliminationsmethode müssen wir die Array-Indizes beachten. Python-Iterables wie Listen und Arrays beginnen oft bei Index 0 und enden bei Index n-1.

Wir können dann den Inhalt der erweiterten Matrix lesen, die Eliminationsmethode anwenden, die Substitution rückgängig machen und schließlich die Antwort anzeigen.

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)

Ausgabe:

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]

Im obigen Code haben wir die Funktion array und die Funktion verwendet, die von der NumPy-Bibliothek bereitgestellt werden, um eine Matrix zu erstellen und absolute Werte zu lesen. Wir haben auch Linalg verwendet; eine NumPy-Unterbibliothek, die zum Ausführen von Operationen wie der Berechnung von Eigenwerten und Vektoren und Determinanten verwendet wird.

from numpy import array, zeros, fabs, linalg

Die Variable a stellt die konstanten Koeffizienten der Variablen dar, während b die Matrix der konstanten Terme darstellt, die den Gleichungen entsprechen. Im Allgemeinen werden die Matrix und der Vektor aus der Matrixdarstellung von vier linearen Gleichungen gezogen, die wir lösen wollen.

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)

Das Ziel der Eliminierung besteht darin, die Matrix so zu transformieren, dass sie Nullen in der führenden Diagonale hat.

Bei der Gaußschen Elimination sollte die erste Reihe unverändert bleiben. Wenn das Pivot-Element in der ersten Zeile jedoch gleich Null ist, müssen wir die erste Zeile mit jeder nachfolgenden Zeile ohne Null als führendes Element tauschen.

Wir haben die k-Schleife verwendet, um die feste Zeile im obigen Code anzugeben. In diesem Fall stellt k den Index der festen Reihe dar, der von 1 bis n-1 reicht.

Die nachfolgenden if-Anweisungen prüfen, ob das Pivot-Element Null ist und tauschen es mit der entsprechenden Zeile aus. Das verwendete i repräsentiert den Index der Zeilen unterhalb der festen Zeile.

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

Im folgenden Code verwenden wir eine Schleife, um sicherzustellen, dass die Eliminierung nicht in der ersten Zeile angewendet wird.

Wir haben dann jede Zeile unterhalb der festen Zeile mit einem Faktor multipliziert. Ein Faktor ist gleich dem diagonalen Element, das der festen Reihe entspricht, dividiert durch das erste Element in der entsprechenden Reihe.

Als nächstes subtrahieren wir die beiden Zeilen von der festen Zeile, um das Element in der zweiten Spalte unter der Hauptdiagonalen zu erhalten. In gleicher Weise können auch die nachfolgenden Reihen gekürzt werden.

Wie bereits erwähnt, bleibt die erste Zeile immer unverändert. Wir haben auch den Faktor und den b-Vektor für jede Zeile wie unten gezeigt berechnet.

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)

Der Begriff Rücksubstitution bezieht sich darauf, dass die Werte von x von der letzten Gleichung zur ersten Gleichung berechnet werden können. In diesem Fall nutzen wir die entsprechenden Werte des Skalarvektors b multipliziert mit den Koeffizienten von x. Das Ergebnis wird dann durch das entsprechende Element in der Hauptdiagonale dividiert.

Wir haben die Berechnung der x-Werte beginnend mit n und der verbleibenden x-Werte von n-1 bis zum ersten priorisiert. Andererseits sollte die Summation bei j+1 beginnen, da uns aus den vorherigen Berechnungen nur die Werte von x bekannt sind. Außerdem ergibt sich der Wert der Summe aus den Werten von x multipliziert mit den Werten von a. Wir können jetzt den Wert von x berechnen, da wir den Wert der Summe haben.

# 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)

Ausgabe:

[[ 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]

Die folgende NumPy-Funktion sollte ebenfalls zu derselben Antwort führen.

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

Ausgabe:

Solution by Nuumpy:
[0.02170543 0.79224806 1.05116279 0.15813953 0.03100775]