6    Gleichungssysteme

In diesem Kapitel lernen Sie, wie Sie mit dem Gauß-Algorithmus, mit der cramerschen Regel und mit Matrizenoperationen lineare Gleichungssysteme lösen können. Außerdem werden mit NumPy-Funktionen und SymPy-Methoden sowohl lineare als auch nichtlineare Gleichungssysteme gelöst.

In den Ingenieurwissenschaften kommen Aufgabenstellungen vor, die die Lösung linearer Gleichungssysteme mit bis zu mehreren Tausend Unbekannten verlangen. Solche komplexen Aufgabenstellungen sind zwangsläufig nur mit Rechenmaschinen zu bewältigen. Daher ist es nicht überraschend, dass der erste Computer von einem Ingenieur konstruiert wurde. Der Bauingenieur Konrad Zuse (1910–1995) löste bereits mit seinem elektromechanischen Relaisrechner lineare Gleichungssysteme, um die Kräfteverteilung in Fachwerken zu ermitteln. Auch in elektrischen Netzwerken werden Ströme und Spannungen mithilfe linearer Gleichungssysteme berechnet. Jeweils ein Anwendungsbeispiel zu Fachwerken und elektrischen Netzwerken finden Sie in [Arens: 514 und 518].

Auch bei Problemstellungen aus der Mathematik, z. B. bei Partialbruchzerlegungen, Interpolationen und Approximationen von Kurven sowie beim Lösen von Differenzialgleichungen, müssen Gleichungssysteme gelöst werden.

Neben den linearen Gleichungssystemen gibt es auch nichtlineare Gleichungssysteme. In diesen Gleichungssystemen kommen Produkte oder höhere Potenzen als 1 der Unbekannten vor. Das ist z. B. der Fall, wenn die Schnittpunkte von zwei Parabeln berechnet werden sollen. Schauen wir uns aber zunächst die linearen Gleichungssysteme an.

6.1    Lineare Gleichungssysteme

Lineare Gleichungssysteme (LGS) haben die allgemeine Form:

Das System hat genauso viele Gleichungen (Zeilen) wie Unbekannte, nämlich x₁ bis x_n. In der ersten Spalte steht die Unbekannte x₁, in der zweiten Spalte steht die Unbekannte x₂, und in der n-ten Spalte steht die Unbekannte x_n. Die Faktoren a₁₁ bis a_nn werden als Koeffizienten bezeichnet. Der erste Index gibt an, in welcher Zeile der Koeffizient steht, und der zweite Index kennzeichnet die Position des Koeffizienten in einer Spalte.

Ein LGS zu lösen heißt, für die Unbekannten x₁ bis x_n Werte zu finden, die für jede Zeile die Gleichung erfüllen. Es handelt sich also um eine sehr anspruchsvolle Problemstellung, die einige Vorarbeiten verlangt.

In Matrizenschreibweise können LGS übersichtlicher dargestellt werden:

LGS können entweder mit einer direkten Methode, z. B. mit dem Gauß-Algorithmus, oder mit iterativen Methoden, z. B. mit dem Jacobi- oder Gauß-Seidel-Algorithmus, gelöst werden.

Bevor Gleichungssysteme mit n Unbekannten behandelt werden, möchte ich Ihnen einige wichtige Begriffe, wie lineare Abhängigkeit bzw. Unabhängigkeit, Lösungsmenge und Lösbarkeit, am Beispiel eines Gleichungssystems mit zwei Unbekannten zeigen. Danach zeige ich Ihnen ab Abschnitt 6.1.2 einige Algorithmen, mit denen Sie LGS lösen können.

6.1.1    Lösbarkeit eines linearen Gleichungssystems

Betrachten wir folgendes Gleichungssystem:

Die Gleichungen sind miteinander gekoppelt, weil die Unbekannten in beiden Gleichungen vorkommen. Die Konstanten a bis d sind gegeben, gesucht werden die Werte für die Unbekannten x und y. Die Frage lautet nun: Welche Bedingungen müssen erfüllt sein, damit die Gleichung lösbar ist?

Beide Gleichungen können nach y aufgelöst werden:

Sie können erkennen, dass es sich um zwei Gleichungen handelt, die jeweils eine Gerade beschreiben. Die Parameter a und c repräsentieren die Steigungen und die Parameter b und d die Achsenabschnitte auf der y-Achse. Eine Lösung (x, y) dieses Gleichungssystems kann also geometrisch als Schnittpunkt dieser zwei Geraden interpretiert werden.

Mit dem Matplotlib-Programm aus Abbildung 6.1 können Sie die Positionen beider Geraden verändern. Mit den Slidern a und c werden die Steigungen und mit den Slidern b und d die Achsenabschnitte auf der y-Achse gewählt.

Abbildung 6.1     Schnittpunkte zweier linearer Funktionen

Beim Programmstart schneiden sich die Geraden bei x = 2 und y = 3. Das ist genau die Lösung des Gleichungssystems. Oder mengentheoretisch formuliert: Die Lösungsmenge beträgt: L = {(2, 3)}. Die Lösung kann auch als Zeilenvektor (2, 3) oder als Spaltenvektor (2, 3)^T dargestellt werden. Das hochgestellte T (Transponieren) bedeutet, dass die Zeile des Vektors in eine Spalte umgewandelt werden muss.

Wenn die beiden Geraden sich schneiden, dann sind beide Funktionsgleichungen linear unabhängig, und es existiert eine Lösung. Wenn die beiden Geraden dagegen parallel verlaufen, dann nennt man die beiden Funktionsgleichungen linear abhängig, und es existiert keine Lösung.

Lösbarkeit

1. Bedingung: Ein lineares Gleichungssystem mit zwei Unbekannten ist genau dann lösbar, wenn die Geraden nicht parallel verlaufen. Sie dürfen also nicht die gleiche Steigung haben. Anders formuliert: Beide Gleichungen müssen linear unabhängig voneinander sein.

2. Bedingung: Die Anzahl der Gleichungen muss mit der Anzahl der Unbekannten übereinstimmen. Für n Unbekannte müssen n Gleichungen gegeben sein.

Die Sonderfälle von unterbestimmtem oder überbestimmtem Gleichungssystem werden hier nicht behandelt.

Diese Bedingungen sind allerdings noch unvollständig. Sie werden weiter unten noch präzisiert.

Mit Listing 6.1 können Sie beide Fälle (Lösung, keine Lösung) simulieren, indem Sie die Slider-Knöpfe entsprechend verschieben.

01  #01_sld_schnittpunkt.py
02  import numpy as np
03  import matplotlib.pyplot as plt
04  from matplotlib.widgets import Slider
05  #1. Funktion
06  def f1(x,a=1,b=1):
07      return a*x+b
08  #2. Funktion
09  def f2(x,c=-1,d=5):
10      return c*x+d
11  #Slider-Einstellungen erfassen
12  def update(val):
13      a = sldA.val
14      b = sldB.val
15      c = sldC.val
16      d = sldD.val
17      y1.set_data(x,f1(x,a,b))
18      y2.set_data(x,f2(x,c,d))
19      fig.canvas.draw_idle()
20  #Grafikbereich
21  fig, ax = plt.subplots(figsize=(6,6))
22  x = np.arange(-10,10,0.001)
23  y1,y2 = ax.plot(x,f1(x,1,1),x,f2(x,-1,5),lw=1.5)
24  ax.set_title(r'$y_{1}=ax+b,\  y_{2}=cx+d$')
25  ax.set_xlabel('x')
26  ax.set_ylabel('y',rotation=True)
27  ax.grid(True)
28  ax.set_xlim(-10,10)
29  ax.set_ylim(-10,10)
30  ax.set_xticks(np.arange(-10, 11, 1))
31  ax.set_yticks(np.arange(-10, 11, 1))
32  fig.subplots_adjust(left=0.12, bottom=0.26)
33  #x-, y-Position, Laenge, Hoehe der Slider
34  xyA = fig.add_axes([0.1, 0.15, 0.8, 0.03])
35  xyB = fig.add_axes([0.1, 0.10, 0.8, 0.03])
36  xyC = fig.add_axes([0.1, 0.05, 0.8, 0.03])
37  xyD = fig.add_axes([0.1, 0.00, 0.8, 0.03])
38  #Slider Objekte erzeugen
39  sldA=Slider(xyA,'a', -5.0, 5.0,valinit=1,valstep=0.1)
40  sldB=Slider(xyB,'b', -10.0, 10.0,valinit=1,valstep=0.1)
41  sldC=Slider(xyC,'c',-5.0,5.0,valinit=-1,valstep=0.1)
42  sldD=Slider(xyD,'d',-10.0,10.0,valinit=5,valstep=0.1)
43  #Änderungen abfragen
44  sldA.on_changed(update)
45  sldB.on_changed(update)
46  sldC.on_changed(update)
47  sldD.on_changed(update)
48  plt.show()

Listing 6.1     Den Schnittpunkt von zwei Geraden simulieren

Analyse

Das Matplotlib-Programm besteht aus fünf Teilen: 1. aus dem Importieren der Module (Zeilen 02 bis 04), 2. aus den Funktionsdefinitionen (Zeilen 06 bis 19), 3. aus der Initialisierung der Funktionsgraphen (Zeilen 22 und 23), 4. aus der Gestaltung der Programmoberfläche (Zeilen 24 bis 42) und 5. aus der Ereignisverarbeitung (Zeilen 44 bis 47).

In Zeile 21 erzeugt die Matplotlib-Methode subplots() die beiden Objekte fig und ax.

Die Axes-Methoden set_title(), set_xlabel() und set_ylabel() erzeugen die Überschrift (Zeile 24) und die Achsenbeschriftungen (Zeilen 25 und 26). In Zeile 27 bewirkt die Axes-Methode grid(True), dass die Gitternetzlinien des Funktionsplots angezeigt werden. Die Axes-Methoden set_xlim() und set_ylim() begrenzen den Wertebereich auf der x- und y-Achse (Zeilen 28 und 29). In den Zeilen 30 und 31 legen die Axes-Methoden set_xticks() und set_yticks() die Skalierung der x- und y-Achse fest.

Das fig-Objekt dient dazu, die Zeichenfläche und die Slider-Steuerelemente auf der Programmoberfläche zu positionieren. Die Figure-Methode subplots_adjust() in Zeile 32 legt den linken und unteren Rand der Zeichenfläche fest. In den Zeilen 34 bis 37 bestimmt die Figure-Methode add_axes() die Positionen und Abmessungen der Slider-Steuerelemente.

Die Matplotlib-Methoden on_changed(update) in den Zeilen 44 bis 47 überprüfen ständig, ob eine Slider-Einstellung geändert wurde. Wenn das der Fall ist, rufen sie die selbst definierte Funktion update(val) (Zeile 12) auf und übergeben dieser Funktion die aktuelle Slider-Einstellung. Die Matplotlib-Methoden set_data(x, ...) (Zeilen 17 und 18) speichern die aktuellen Koordinatendaten der mathematischen Funktionen temporär ab, und plt.show() in Zeile 48 bewirkt, dass die mathematischen Funktionen auf dem Bildschirm angezeigt werden.

Die Figure-Methode canvas.draw_idle() in Zeile 19 hat eigentlich die Aufgabe, die Funktionsplots zu aktualisieren, wenn die Slider-Einstellungen geändert wurden. Sie muss aber nicht zwingend verwendet werden. Sie können diese Methode auskommentieren und das Programm erneut testen.

6.1.2    Der Gauß-Algorithmus

Der Gauß-Algorithmus stellt ein Eliminationsverfahren bereit, mit dem sich lineare Gleichungssysteme nach einem formal festgelegten Schema einfach lösen lassen. Das Gleichungssystem wird zunächst in Matrixschreibweise

dargestellt.

Abbildung 6.2 zeigt eine Übersicht über die allgemein üblichen Begriffsbestimmungen.

Abbildung 6.2     Allgemeines Schema der Matrixdarstellung eines linearen Gleichungssystems

Die 3×3-Matrix wird als Koeffizientenmatrix bezeichnet. Rechts daneben steht der Lösungsvektor, und auf der rechten Seite des Systems befindet sich der Inhomogenitätsvektor.

Erweiterung des Vektorbegriffs

Gewöhnlich denken Sie vielleicht bei Vektoren an Begriffe aus der Physik, z. B. an Kräfte oder Feldstärken. In der linearen Algebra wird der Vektor-Begriff jedoch allgemeiner definiert: Vektoren sind hier Strukturen in einem n-dimensionalen Vektorraum. Nach diesem Konzept können auch Funktionen und andere Objekte Vektoren sein [Scherfner: 188].

Als Nächstes wird das Gleichungssystem in die sogenannte erweiterte Koeffizientenmatrix überführt:

Folgende Umformungen können an dem Gleichungssystem vorgenommen werden, ohne dass sich der Lösungsvektor ändert:

1. Zeilen können vertauscht werden.

2. Jede Zeile kann mit einem Faktor multipliziert werden.

3. Zu einer ausgewählten Zeile kann eine andere Zeile addiert werden.

Das Ziel des Gauß-Algorithmus ist es, die Koeffizientenmatrix mithilfe der drei genannten Operationen so umzuformen, dass unterhalb der Hauptdiagonalen nur Nullen stehen. Diese Form wird auch als Stufen- oder obere Dreiecksform bezeichnet.

Hauptdiagonale

Die Hauptdiagonale einer quadratischen Matrix wird durch die Elemente von der oberen linken Ecke bis zur unteren rechten Ecke der Koeffizientenmatrix gebildet: a₁₁, a₂₂ … a_ii.

Liegt das System in der oberen Dreiecksform vor, dann lässt sich der Lösungsvektor durch Rückwärtseinsetzen (Rücksubstitution) einfach berechnen.

Grobstruktur des Gauß-Algorithmus

1. Die erste Gleichung wird bei allen folgenden Umformungen unverändert übernommen.

2. Mithilfe der ersten Gleichung wird die erste Unbekannte x₁ in allen folgenden Gleichungen eliminiert.

3. Die zweite Gleichung des letzten Systems wird unverändert übernommen. Sie ist jetzt Bezugsgleichung. Mit ihrer Hilfe wird in den folgenden Gleichungen die zweite Unbekannte x₂ eliminiert.

4. Das Verfahren wird so lange fortgesetzt, bis ein lineares Gleichungssystem entstanden ist, in dem jede folgende Gleichung mindestens eine Unbekannte weniger erhält als die Gleichung darüber. Wenn das Gleichungssystem lösbar ist, entsteht so eine obere Dreiecksform [vgl. Bossek: 98].

Der Gauß-Algorithmus soll am Beispiel eines Gleichungssystems mit drei Unbekannten besprochen werden:

Wir starten mit der erweiterten Koeffizientenmatrix:

1. Elimination von a₂₁: z₂ = z₂ – (1/3) z₁, allgemein: (a₂₁/a₁₁)

2. Elimination von a₃₁: z₃ = z₃ – (2/3) z₁, allgemein: (a₃₁/a₁₁)

3. Elimination von a₃₂: z₃ = z₃ – (–1/4) z₂, allgemein: (a₃₂/a₂₂)

Durch Rückwärtseinsetzen erhält man für x₃ = 3, für x₂ = 2 und für x₁ = 1.

Die Teiler a₁₁ und a₂₂ auf der Hauptdiagonalen werden als Pivotelemente bezeichnet.

Pivotelement

Die Diagonalelemente a_ii mit Ausnahme des letzten Diagonalelements heißen Pivotelemente. Die erste Zeile wird als Pivotzeile bezeichnet. Der aus dem Französischen stammende Begriff pivot bedeutet so viel wie »Dreh- und Angelpunkt«.

Für die Implementierung des Gauß-Algorithmus benötigt man drei ineinander verschachtelte Schleifen. Die äußere Schleife iteriert die Diagonalelemente. Die mittlere Schleife berechnet die aktuellen Werte des Inhomogenitätsvektors, und die innere Schleife berechnet die Koeffizienten der oberen Dreiecksmatrix.

Berechnung der oberen Dreiecksmatrix

Das Problem »Wie löst man ein LGS mit dem Gauß-Algorithmus?« sollten Sie in zwei Teilbereiche aufteilen, um es besser nachvollziehen zu können. Als Erstes erzeugt Listing 6.2 nur die obere Dreiecksmatrix, anschließend berechnet Listing 6.3 den Lösungsvektor durch Rücksubstitution.

01  #02_gauss1.py
02  import numpy as np
03  #Koeffizientenmatrix
04  a = np.array([[3,2,1],
05                [1,2,3],
06                [2,1,4]],dtype=float)
07  #Inhomogenitätsvektor
08  b = np.array([10,14,16],dtype=float)
09  #obere Dreiecksmatrix erzeugen
10  n=len(b)
11  for k in range(0,n):         #Diagonale
12      pivot=a[k,k]
13      for i in range(k+1,n):   #Zeilen
14          f = a[i,k]/pivot
15          b[i]=b[i]-f*b[k]
16          for j in range(k,n): #Spalten
17              a[i,j]=a[i,j]-f*a[k,j]
18          print(a)
19  print("Inhomogenitätsvektor\n",b)

Listing 6.2     Gauß-Algorithmus für Dreiecksmatrix

Ausgabe

[[3.         2.         1.        ]
 [0.         1.33333333 2.66666667]
 [2.         1.         4.        ]]
[[ 3.          2.          1.        ]
 [ 0.          1.33333333  2.66666667]
 [ 0.         -0.33333333  3.33333333]]
[[3.         2.         1.        ]
 [0.         1.33333333 2.66666667]
 [0.         0.         4.        ]]
Inhomogenitätsvektor
 [10.         10.66666667 12.        ]

Analyse

Zeilen 11 bis 17: Die erste for-Schleife iteriert über alle Zeilen und Diagonalelemente a[k,k].

Zeile 12: Die Elemente der Hauptdiagonalen a[k,k] sind die Pivotelemente. Sie werden in die Variable pivot gespeichert.

Zeilen 13 bis 15: Die zweite for-Schleife iteriert über alle Zeilen unterhalb der Diagonalelemente. In Zeile 15 wird der aktuelle Inhomogenitätsvektor b[i] berechnet.

Zeilen 16 und 17: Die dritte for-Schleife iteriert über alle Elemente der Zeile aus der zweiten Schleife. Ein Vielfaches von a[k,j] wird von der aktuellen Zeile subtrahiert. In Zeile 17 werden unterhalb der Diagonalelemente Nullen erzeugt. Der in Zeile 14 berechnete Faktor f wird so gewählt, dass diese Forderung erfüllt wird.

Das Programm hat jedoch zwei Mängel: Es überprüft nicht, ob ein Diagonalelement 0 ist. Wäre das der Fall, dann käme es in Zeile 14 zu einer Division durch 0, und das Programm würde abstürzen. Ein weiterer Mangel besteht darin, dass keine Pivotisierung durchgeführt wird, d. h., die Zeilen des Gleichungssystems werden nicht so vertauscht, dass jeweils links oben das betragsmäßig größte Element steht. In dem Testbeispiel war diese Bedingung erfüllt, daher konnte der Gauß-Algorithmus umgesetzt werden.

Pivotisierung

Auch wenn sich der Gauß-Algorithmus leicht auf Rechenmaschinen implementieren lässt, müssen Sie beachten, dass bei der naiven Anwendung dieses Verfahrens nicht mehr akzeptable Rundungsfehler entstehen können. Durch die Pivotisierung können diese Fehler reduziert werden: Das Programm wählt für die Elimination das betragsmäßig größte Element aus einer Spalte und tauscht die entsprechenden Zeilen [Arens: 522].

Rücksubstitution

Listing 6.3 berechnet die obere Dreiecksmatrix des Testbeispiels und ermittelt anschließend den Lösungsvektor durch Rücksubstitution mit der selbst definierten Python-Funktion gauss(a,b):

01  #03_gauss2.py
02  import numpy as np
03  a = np.array([[3,2,1],
04                [1,2,3],
05                [2,1,4]],dtype=float)
06  b = np.array([10,14,16],dtype=float)
07  #Gauss-Algorithmus
08  def gauss(a,b):
09      #obere Dreiecksmatrix
10      n=len(b) #Anzahl der Gleichungen
11      for k in range(0,n):         #Diagonale
12          pivot=a[k,k]
13          for i in range(k+1,n):   #Zeilen
14              f = a[i,k]/pivot
15              b[i]=b[i]-f*b[k]
16              for j in range(k,n): #Spalten
17                  a[i,j]=a[i,j]-f*a[k,j] 
18  #Rücksubstitution
19      x = np.empty(n)
20      for i in range(n-1, -1, -1):
21          for j in range(i+1,n):
22              b[i]=b[i]-a[i,j]*x[j]
23          x[i] = b[i]/a[i,i]
24      return x
25  #Ausgabe 
26  print("Koeffizientenmatrix\n", A)
27  print("Inhomogenitätsvektor\n",b)
28  print("Lösungsvektor\n",gauss(A,b))

Listing 6.3     Gauß-Algorithmus mit Rücksubstitution

Ausgabe

Koeffizientenmatrix
 [[3. 2. 1.]
 [1. 2. 3.]
 [2. 1. 4.]]
Inhomogenitätsvektor
 [10. 14. 16.]
Lösungsvektor
 [1. 2. 3.]

Analyse

Der Algorithmus für die Rücksubstitution arbeitet wie folgt: In der 1. Schleife (Zeile 20) werden die Zeilen des Gleichungssystems mit dem Startwert n-1 bis -1 mit der Schrittweite -1 von unten nach oben durchlaufen.

In der 2. Schleife (Zeile 21) werden die bekannten Variablen auf die rechte Seite gebracht. Zum Schluss wird durch den Koeffizienten a[i,i] des aktuellen x-Werts geteilt. In Zeile 24 wird der Lösungsvektor x zurückgegeben.

6.1.3    Der Gauß-Jordan-Algorithmus

Der Gauß-Jordan-Algorithmus erzeugt aus der Koeffizientenmatrix eine Diagonalmatrix. Die Gleichungen werden so lange umgeformt, bis nur noch die Elemente der Hauptdiagonalen vorkommen und diese den Wert 1 haben. Alle Elemente über und unter der Hauptdiagonalen der Koeffizientenmatrix verschwinden.

Dieser zusätzliche Aufwand lohnt sich eigentlich nicht. Sie sollten sich den Gauß-Jordan-Algorithmus aber anschauen, weil er wieder bei der Berechnung der Inversen einer Koeffizientenmatrix vorkommt. Darum geht es in Abschnitt 6.1.5.

Die Grobstruktur des Gauß-Jordan-Algorithmus

Der Gauß-Jordan-Algorithmus lässt sich kurz wie folgt umreißen:

1. Schritt: die 1. Zeile durch a₁₁ teilen

2. Schritt: alle Elemente unterhalb des Elements a₁₁ eliminieren

3. Schritt: die 2. Zeile durch a₂₂ teilen

4. Schritt: alle Elemente unterhalb und oberhalb des Elements a₂₂ eliminieren

5. Schritt: die 3. Zeile durch a₃₃ teilen

6. Schritt: alle Elemente unterhalb und oberhalb des Elements a₃₃ eliminieren

usw.

Wie dieser Algorithmus praktisch anzuwenden ist, demonstriert ein Eingangsbeispiel. Der Ausgangspunkt ist wieder die erweiterte Koeffizientenmatrix. Führen Sie dann die einzelnen Schritte durch:

1 a) z₁ = z₁/a₁₁, (a₁₁ = 3)

1 b) Elimination von a₂₁: z₂ = z₂ – a₂₁z₁, (a₂₁ = 1)

1 c) Elimination von a₃₁: z₃ = z₃ – a₃₁z₁, (a₃₁ = 2)

2 a) z₂ = z₂/a_22, (a₂₂ = 4/3)

2 b) Elimination von a₁₂: z₁ = z₁ – a₁₂z₂, (a₁₂ = 2/3)

2 c) Elimination von a₃₂: z₃ = z₃ – a₃₂z₂, (a₃₂ = –1/3)

3 a) z₃ = z₃/a₃₃, (a₃₃ = 4)

3 b) Elimination von a₁₃: z₁ = z₁ – a₁₃z₂, (a₁₃ = –1)

3 c) Elimination von a₂₃: z₂ = z₂ – a₂₃z₃, (a₂₃ = 2)

Die Lösungen können jetzt direkt abgelesen werden.

Listing 6.4 löst ein lineares Gleichungssystem mit dem Gauß-Jordan-Algorithmus. Wenn der Kommentar in Zeile 20 entfernt wird, dann gibt das Programm den Fortschritt der einzelnen Umformungen aus.

01  #04_jordan.py
02  import numpy as np
03  A = np.array([[3,2,1],
04                [1,2,3],
05                [2,1,4]],dtype=float)
06  b=np.array([10,14,16],dtype=float)
07  #Gauss-Jordan-Algorithmus
08  def jordan(a,b):
09      n=len(b)           #Anzahl der Zeilen
10      for k in range(n): #Diagonale 
11          pivot=a[k,k]
12          for j in range(k,n): #Spalten
13              a[k,j]=a[k,j]/pivot
14          b[k]=b[k]/pivot 
15          for i in range(n):   #Zeilen
16              if i==k:continue
17              f=a[i,k]
18              for j in range(k,n):
19                  a[i,j]=a[i,j] - f*a[k,j]
20              #print(a)
21              b[i]=b[i]-f*b[k]
22      return b,a
23  x,A_t=jordan(A,b)
24  #Ausgabe
25  print("Lösungsvektor\n ",x)
26  print("Transformierte Koeffizientenmatrix\n",A_t)

Listing 6.4     Der Gauß-Jordan-Algorithmus

Ausgabe

[[1.         0.66666667 0.33333333]
 [0.         1.33333333 2.66666667]
 [2.         1.         4.        ]]
[[ 1.          0.66666667  0.33333333]
 [ 0.          1.33333333  2.66666667]
 [ 0.         -0.33333333  3.33333333]]
[[ 1.          0.         -1.        ]
 [ 0.          1.          2.        ]
 [ 0.         -0.33333333  3.33333333]]
[[ 1.  0. -1.]
 [ 0.  1.  2.]
 [ 0.  0.  4.]]
[[1. 0. 0.]
 [0. 1. 2.]
 [0. 0. 1.]]
[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]
Lösungsvektor
  [1. 2. 3.]
Transformierte Koeffizientenmatrix
 [[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]

Analyse

Die Laufvariable k der äußeren Schleife (Zeile 10) indiziert die Hauptdiagonale der Koeffizientenmatrix A.

In Zeile 11 werden die Diagonalelemente a[k,k] in die Hilfsvariable pivot gespeichert.

In Zeile 13 werden die Elemente a[k,j] spaltenweise durch die Variable pivot dividiert.

In Zeile 14 werden die Zeilenelemente des Inhomogenitätsvektors b[k] durch die Variable pivot dividiert.

In Zeile 15 indiziert die Laufvariable i der for-Schleife den Zeilenindex i. Elemente der Hauptdiagonalen dürfen nicht berücksichtigt werden (Zeile 16). In Zeile 19 wird die transformierte Koeffizientenmatrix berechnet. In Zeile 21 wird der Lösungsvektor b berechnet.

In Zeile 22 gibt die return-Anweisung den Lösungsvektor b und die transformierte Koeffizientenmatrix zurück. Die transformierte Koeffizientenmatrix wurde nur zur Kontrolle mit zurückgegeben.

Vergleich der Laufzeitkomplexitäten

Der Gauß-Algorithmus benötigt für die Berechnung der oberen Dreiecksmatrix drei ineinander verschachtelte Schleifen. Hierfür beträgt der Aufwand O(n³). Für die Rücksubstitution benötigt er zwei ineinander verschachtelte Schleifen, der Aufwand beträgt hier O(n²). Insgesamt beträgt die Laufzeitkomplexität also O(n³) + O(n²).

Der Gauß-Jordan-Algorithmus hat zwar die gleiche formale Struktur, er benötigt für die Elimination der Elemente oberhalb der Hauptdiagonalen jedoch zusätzliche Rechenoperationen.

Beide Algorithmen sind deshalb nicht für die Lösung großer linearer Gleichungssysteme geeignet.

6.1.4    Lösung mit Determinanten: Die cramersche Regel

Die cramersche Regel beschreibt, wie sich Gleichungssysteme mithilfe von Determinanten lösen lassen. Was unter einer Determinante D zu verstehen ist, zeigt wieder ein einfaches Beispiel. Betrachten wir folgendes Gleichungssystem:

Dieses Gleichungssystem kann nach den Unbekannten x₁ und x₂ aufgelöst werden:

Für die Zähler ist auch folgende Notation üblich:

und

Den Nenner des Bruchs bezeichnet man als Systemdeterminante.

In Determinantenschreibweise gilt für die Lösung von x₁:

und x₂:

Die ausführliche Herleitung finden Sie in [Papula2: 24].

Determinante

Den Eintrag zwischen den beiden senkrechten Strichen bezeichnet man als Determinante. Neben den zweireihigen Determinanten gibt es noch mehrreihige Determinanten. Die beiden senkrechten Striche sind Operatoren. Sie verlangen, dass die Determinante einer quadratischen Matrix berechnet werden soll. Das Ergebnis ist eine Zahl.

Ein lineares Gleichungssystem ist nur dann lösbar, wenn die Systemdeterminante D ungleich 0 ist. Eine andere Bestimmung besagt, dass die Koeffizientenmatrix des Gleichungssystems regulär sein muss.

Reguläre Matrix

Eine Matrix A heißt regulär, wenn ihre Determinante einen von 0 verschiedenen Wert hat. Andernfalls ist sie singulär. Diese Definition gilt nur für quadratische Matrizen. Nur wenn die Koeffizientenmatrix eines Gleichungssystems regulär ist, ist es lösbar.

Aus der allgemeinen Lösung eines Gleichungssystems mit zwei Unbekannten kann abgelesen werden, wie eine zweireihige Determinante berechnet werden muss:

Eine dreireihige Determinante lässt sich noch relativ einfach mit der Regel von Sarrus berechnen:

Der Rechenaufwand steigt aber unverhältnismäßig stark an, wenn die Determinante mehr als drei Reihen hat. Deshalb wird dringend davon abgeraten, Gleichungssysteme mit der cramerschen Regel zu berechnen:

»Es wäre wahnsinnig, auf diese Weise Gleichungen lösen zu wollen.« [Strang: 273]

Ich möchte dennoch nicht darauf verzichten, Ihnen zu zeigen, wie Sie mit der cramerschen Regel ein lineares Gleichungssystem mit drei Unbekannten lösen können.

Wenn die obere Dreiecksmatrix gegeben ist, dann müssen nur die Elemente der Hauptdiagonalen a_ii miteinander multipliziert werden:

Zuerst zeige ich Ihnen in Listing 6.5, wie Sie mit dem Gauß-Algorithmus eine Determinante berechnen können:

01  #05_determinate.py
02  import numpy as np
03  from numpy.linalg import det
04  #Koeffizientenmatrix
05  a = np.array([[3,2,1],
06                [1,2,3],
07                [2,1,4]],dtype=float)
08  #Gauss-Algorithmus
09  def determinante(a):
10      dm=1
11      n=a.shape[0]
12      for k in range(0,n): 
13          for i in range(k+1,n): 
14              f = a[i,k]/a[k,k]
15              for j in range(k,n): 
16                  a[i,j]=a[i,j]-f*a[k,j]
17          dm=dm*a[k,k]
18      return dm
19  #Ausgabe
20  print("Determinante")
21  print("Gauss:",determinante(a))
22  print("NumPy:",det(a))

Listing 6.5     Eine Determinante mit dem Gauß-Algorithmus berechnen

Ausgabe

Determinante
Gauss: 16.0
NumPy: 16.000000000000007

Analyse

Zeile 03 importiert die Funktion det aus dem Untermodul linalg des Moduls numpy.

In den Zeilen 12 bis 16 berechnet die selbst definierte Python-Funktion determinante(a) (Zeile 09) zunächst aus der Koeffizientenmatrix a die obere Dreiecksmatrix. Die Determinante dm wird in Zeile 17 iterativ aus den Produkten der Hauptdiagonalelemente a[k,k] berechnet.

Anwendung der cramerschen Regel

Das nächste Python-Programm berechnet den Lösungsvektor des Testbeispiels mit der NumPy-Funktion det() aus dem Untermodul linalg.

Nach der Regel von Cramer werden die drei Unbekannten des Gleichungssystems mithilfe von Determinanten berechnet. Für die Unbekannte x₁ gilt:

Für die Unbekannte x₂ gilt:

Für die Unbekannte x₃ gilt:

Diese Rechenschritte lassen sich direkt in einen Python-Quelltext übertragen. Listing 6.6 zeigt die Umsetzung:

01  #06_cramer.py
02  import numpy as np
03  from numpy.linalg import det
04  #Koeffizientenmatrix
05  a = np.array([[3,2,1],
06                [1,2,3],
07                [2,1,4]],dtype=float)
08  #Inhomogenitätsvektor
09  b = np.array([10.,14.,16.],dtype=float)
10  #für die 1. Lösung
11  A1= np.array([[b[0],a[0,1],a[0,2]],
12                [b[1],a[1,1],a[1,2]],
13                [b[2],a[2,1],a[2,2]]],dtype=float)
14  #für die 2. Lösung
15  A2= np.array([[a[0,0],b[0],a[0,2]],
16                [a[1,0],b[1],a[1,2]],
17                [a[2,0],b[2],a[2,2]]],dtype=float)
18  #für die 3. Lösung
19  A3= np.array([[a[0,0],a[0,1],b[0]],
20                [a[1,0],a[1,1],b[1]],
21                [a[2,0],a[2,1],b[2]]],dtype=float)
22  #Systemdeterminante
23  D=det(a)
24  D1=det(A1)
25  D2=det(A2)
26  D3=det(A3)
27  #Lösungsvektoren
28  x1=D1/D
29  x2=D2/D
30  x3=D3/D
31  #Ausgabe
32  print("Koeffizientenmatrix\n",a)
33  print("Inhomogenitätsvektor\n",b)
34  print("D  =",D)
35  print("D1 =",D1)
36  print("D2 =",D2)
37  print("D3 =",D3)
38  print("x1 =",x1)
39  print("x2 =",x2)
40  print("x3 =",x3)

Listing 6.6     Lösung eines linearen Gleichungssystems mit der cramerschen Regel

Ausgabe

Koeffizientenmatrix
 [[3. 2. 1.]
 [1. 2. 3.]
 [2. 1. 4.]]
Inhomogenitätsvektor
 [10. 14. 16.]
D  = 16.000000000000007
D1 = 15.999999999999998
D2 = 32.00000000000003
D3 = 48.000000000000014
x1 = 0.9999999999999994
x2 = 2.000000000000001
x3 = 2.9999999999999996

Analyse

Schon die aufwendige Implementierung zeigt, dass es sich eigentlich nicht lohnt, lineare Gleichungssysteme mit der cramerschen Regel zu lösen. Auffällig sind im Vergleich zum Gauß-Algorithmus auch die Rundungsfehler.

Die cramersche Regel sollten Sie vorzugweise nur dann anwenden, wenn Gleichungssysteme mit bis zu drei Unbekannten allgemein gelöst werden müssen.

Vergleich des Gauß-Algorithmus mit der cramerschen Regel

Wie absurd es ist, lineare Gleichungssysteme mit der cramerschen Regel lösen zu wollen, demonstriert auch Tabelle 6.1. Sie vergleicht den Rechenaufwand bei der cramerschen Regel mit dem beim Gauß-Algorithmus.

Zahl der Unbekannten	Cramersche Regel	Gauß-Algorithmus
10	1,08 sec	≈ 10–6 sec
20	95.000 Jahre	≈ 10–5 sec

Tabelle 6.1     Rechenzeit (3 ns pro Punktoperation) [Engeln-Müllges: 117]

6.1.5    Lösung mit der Inversen einer Matrix

In Matrixschreibweise kann ein lineares Gleichungssystem kompakt mit

formuliert werden.

Durch Bildung der inversen Koeffizientenmatrix A^–1 wird diese Gleichung nach dem unbekannten Lösungsvektor

aufgelöst. Da die Inverse einer Matrix auch mit dem Gauß-Jordan-Algorithmus berechnet werden kann, ist die Lösung eines linearen Gleichungssystems mit dieser Methode eigentlich nicht notwendig [Knorrenschild: 43]. Es soll aber trotzdem gezeigt werden, wie mit NumPy die Multiplikation einer Inversen mit einem Spaltenvektor durchgeführt werden kann, denn auch diese Problemstellung kommt in Kapitel 13, »Ausgleichsrechnungen«, wieder vor.

Die Inverse der Matrix A wird wie folgt berechnet:

SymPy ist in der Lage, die Inverse einer Matrix symbolisch zu berechnen (siehe Listing 6.7):

01  #07_sympyInverse.py
02  from sympy import symbols, Matrix
03  a,b,c,d=symbols('a,b,c,d')
04  #Koeffizientenmatrix
05  A=Matrix([[a,b],
06            [c,d]])
07  #Ausgabe
08  print(A)
09  print("\ninverse Matrix\n",A.inv())

Listing 6.7     Inverse einer zweireihigen Matrix, mit SymPy berechnet

Ausgabe

Matrix([[a, b], [c, d]])
inverse Matrix
Matrix([[d/(a*d - b*c),-b/(a*d - b*c)], 
       [-c/(a*d - b*c), a/(a*d - b*c)]])

Analyse

In Zeile 03 werden die symbolischen Variablen a, b, c und d vereinbart. In Zeile 05 erzeugt die SymPy-Methode Matrix() die Matrix A. In Zeile 09 wird die Inverse mit A.inv() berechnet und ausgegeben.

Berechnung der Inversen mit dem Gauß-Jordan-Algorithmus

Die Inverse einer Matrix kann auch mit dem Gauß-Jordan-Algorithmus wie folgt berechnet werden [Weltner2: 140]:

1. Schritt: Erweiterung der Matrix A mit der Einheitsmatrix E

2. Schritt: Durchführung des gauß-jordanschen Eliminationsverfahrens, bis auf der linken Matrixhälfte die Einheitsmatrix erscheint. Die rechte Hälfte der Matrix ist dann die Inverse.

Für das Testbeispiel gilt der Ansatz:

Mit dem Gauß-Jordan-Algorithmus erhält man nach mehreren Umformungen:

Für die Inverse gilt also:

Mit dem Inhomogenitätsvektor b = (10, 14, 16)^T aus dem Testbeispiel erhalten Sie mit der Inversen das Gleichungssystem

mit dem Lösungsvektor L = (1, 2, 3)^T.

NumPy berechnet die Inverse einer Matrix numerisch mit der Funktion inv(A). Die Matrizenmultiplikation erfolgt mit dem Infix-Operator @.

Listing 6.8 löst das lineare Gleichungssystem des Testbeispiels mit der Funktion inv() aus dem Untermodul linalg von NumPy:

01  #08_numpyInverse.py
02  import numpy as np
03  from numpy.linalg import inv
04  #Koeffizientenmatrix
05  A = np.array([[3,2,1],
06                [1,2,3],
07                [2,1,4]],dtype=float)
08  #Inhomogenitätsvektor
09  b = np.array([10.,14.,16.],dtype=float)
10  #Matrizenmultiplikation
11  x=inv(A)@b
12  #Ausgabe
13  print("Koeffizientenmatrix\n", A)
14  print("Inhomogenitätsvektor\n",b)
15  print("Inverse der Koeffizientenmatrix\n",inv(A))
16  print("Lösungsvektor\n",x)
17  print("x1 =",x[0])
18  print("x2 =",x[1])
19  print("x2 =",x[2])

Listing 6.8     Lösung mit der inversen Koeffizientenmatrix

Ausgabe

Koeffizientenmatrix
 [[3. 2. 1.]
 [1. 2. 3.]
 [2. 1. 4.]]
Inhomogenitätsvektor
 [10. 14. 16.]
Inverse der Koeffizientenmatrix
 [[ 0.3125 -0.4375  0.25  ]
 [ 0.125   0.625  -0.5   ]
 [-0.1875  0.0625  0.25  ]]
Lösungsvektor
 [1. 2. 3.]
x1 = 1.0000000000000009
x2 = 1.9999999999999991
x2 = 2.9999999999999996

Analyse

In Zeile 03 wird aus dem NumPy-Untermodul linalg die Funktion inv importiert. In Zeile 11 berechnet diese Funktion die Inverse inv(A) der Koeffizientenmatrix A. Die Matrizenmultiplikation der Inversen von A mit dem Inhomogenitätsvektor b wird mit dem Infix-Operator @ durchgeführt. Der Lösungsvektor wird in die Objektvariable x gespeichert.

Alternativen zu dem Infix-Operator finden Sie mit den NumPy-Funktion dot(inv(A, b) oder matmul(inv(A,b)).

6.1.6    Lösung mit der NumPy-Funktion solve()

Die kompakteste Lösung eines linearen Gleichungssystems besteht in der NumPy-Funktion solve(A, b) aus dem Untermodul linalg. Die Funktion erwartet also nur zwei Argumente: die Koeffizientenmatrix A und den Inhomogenitätsvektor b.

Listing 6.9 zeigt die Umsetzung für das Testbeispiel:

01  #09_numpy_solve.py
02  import numpy as np
03  from numpy.linalg import solve
04  #Koeffizientenmatrix
05  A = np.array([[3,2,1],
06                [1,2,3],
07                [2,1,4]],dtype=float)
08  #Inhomogenitätsvektor
09  b = np.array([10,14,16],dtype=float)
10  #Lösungsvektor
11  L=solve(A,b)
12  #Ausgabe
13  print("Koeffizientenmatrix\n", A)
14  print("Inhomogenitätsvektor\n",b)
15  print("Lösungsvektor\n",L)
16  print("x1 =",L[0])
17  print("x2 =",L[1])
18  print("x3 =",L[2])

Listing 6.9     Lösung mit der NumPy-Funktion »solve«

Ausgabe

Koeffizientenmatrix
 [[3. 2. 1.]
 [1. 2. 3.]
 [2. 1. 4.]]
Inhomogenitätsvektor
 [10. 14. 16.]
Lösungsvektor
 [1. 2. 3.]
x1 = 0.9999999999999997
x2 = 2.0000000000000004
x3 = 3.0

Analyse

In Zeile 03 wird die Funktion solve aus dem Untermodul linalg des Moduls numpy importiert.

In Zeile 11 löst diese Funktion das lineare Gleichungssystem. Sie erwartet die Übergabe der Koeffizientenmatrix A (erster Parameter) und des Spaltenvektors b (zweiter Parameter). Der Lösungsvektor wird in die Objektvariable L gespeichert.

In Zeile 15 wird der Lösungsvektor als ndarray ausgegeben. In den Zeilen 16 bis 18 werden die einzelnen Komponenten x1, x2 und x3 des Lösungsvektors ausgegeben.

6.1.7    Lösung mit SymPy-Methoden

Mit dem Python-Modul SymPy lassen sich lineare Gleichungssysteme ebenfalls lösen. SymPy-Programme benötigen hierfür nur wenige Programmzeilen. Der besondere Vorteil der SymPy-Methoden besteht darin, dass Gleichungssysteme sowohl numerisch als auch symbolisch gelöst werden können. SymPy stellt die vier Methoden

• solve([g1,g2,g3],[x1,x2,x3])

• A.inv()*b

• linsolve([A,b],[x1,x2,x3])

• A.gauss_jordan_solve(b)

für die Lösung linearer Gleichungssysteme zur Verfügung.

Die Methode solve()

Das erste Beispiel aus Listing 6.10 zeigt, wie Sie mit der SymPy-Methode solve() ein lineares Gleichungssystem lösen können:

01  #10_sympy_lgs1.py
02  from sympy import symbols, solve
03  x1,x2,x3=symbols('x1,x2,x3')
04  g1= 3*x1+2*x2+  x3-10
05  g2=   x1+2*x2+3*x3-14
06  g3= 2*x1+  x2+4*x3-16
07  L=solve([g1,g2,g3],(x1,x2,x3),dict=True)
08  print("Lösungsvektor\n",L)
09  print("Schlüssel\n",L[0].keys())
10  print("Werte\n",L[0].values())
11  print("x1 =",L[0][x1])
12  print("x2 =",L[0][x2])
13  print("x3 =",L[0][x3])

Listing 6.10     Ein lineares Gleichungssystem mit der SymPy-Methode »solve« lösen

Ausgabe

Lösungsvektor
 [{x1: 1, x2: 2, x3: 3}]
Schlüssel
 dict_keys([x1, x2, x3])
Werte
 dict_values([1, 2, 3])
x1 = 1
x2 = 2
x3 = 3

Analyse

In Zeile 03 werden die symbolischen Variablen x1, x2 und x3 für die Unbekannten vereinbart.

In den Zeilen 04 bis 06 werden die einzelnen Gleichungen des Gleichungssystems in die Objektvariablen g1, g2 und g3 gespeichert.

In Zeile 07 löst die Methode solve() das Gleichungssystem. Die Objektvariablen und die symbolischen Variablen müssen durch eckige oder runde Klammern voneinander getrennt werden. Der Lösungsvektor wird als Dictionary in die Objektvariable L gespeichert.

Zeile 08 gibt den Lösungsvektor als list (Datenstruktur Liste) aus.

Zeile 09 gibt die Schlüssel dict_keys des Lösungsvektors L aus.

Zeile 10 gibt die Werte dict_values des Lösungsvektors L aus.

Wie auf die einzelnen Komponenten des Lösungsvektors zugegriffen werden kann, zeigen die Zeilen 11 bis 13.

Die Methoden inv(), gauss_jordan_solve(b) und linsolve()

Listing 6.11 zeigt, wie Sie lineare Gleichungssysteme mit der Inversen inv() und den SymPy-Methoden gauss_jordan_solve() und linsolve() lösen können:

01  #11_sympy_lgs2.py
02  from sympy import symbols,Matrix,linsolve
03  x1,x2,x3=symbols('x1,x2,x3')
04  A=Matrix([[3,2,1],
05            [1,2,3],
06            [2,1,4]])
07  b=Matrix([10,14,16])
08  inverse=A.inv() #A**-1 
09  L1=inverse*b
10  L2=A.gauss_jordan_solve(b)
11  L3=linsolve([A,b],(x1,x2,x3))
12  #Ausgaben
13  print("Inverse von A\n",inverse)
14  print("Lösung mit Inverse\n",L1)
15  print("x1 = %3.3f\nx2 = %3.3f\nx3 = %3.3f"%\
16        (L1[0],L1[1],L1[2]))
17  print("Lösung mit gauss_jordan_solve\n",L2)
18  print("x1 = %3.3f\nx2 = %3.3f\nx3 = %3.3f"%\
19        (L2[0][0],L2[0][1],L2[0][2]))
20  print("Lösung mit linsolve\n",L3)
21  print("x1 = %3.3f\nx2 = %3.3f\nx3 = %3.3f"%\
22        (L3.args[0][0],L3.args[0][1],L3.args[0][2]))

Listing 6.11     SymPy-Methoden für die Lösung linearer Gleichungssysteme

Ausgabe

Inverse von A
 Matrix([[5/16, -7/16, 1/4], 
         [1/8, 5/8, -1/2], 
         [-3/16, 1/16, 1/4]])
Lösung mit Inverse
 Matrix([[1], [2], [3]])
x1 = 1.000
x2 = 2.000
x3 = 3.000
Lösung mit gauss_jordan_solve
 (Matrix([
[1],
[2],
[3]]), Matrix(0, 1, []))
x1 = 1.000
x2 = 2.000
x3 = 3.000
Lösung mit linsolve
 FiniteSet((1, 2, 3))
x1 = 1.000
x2 = 2.000
x3 = 3.000

Analyse

Alle drei Methoden speichern den Lösungsvektor in die Objektvariablen L1, L2 und L3. Wenn das Programm die numerischen Werte der Lösungen ausgeben oder in nachfolgenden Berechnungen weiter damit arbeiten soll, dann müssen Sie beachten, dass jede einzelne Methode unterschiedliche DatenstruktureWn zurückgibt. Die Wahl der Methode bestimmt also die Datenstruktur der Objektvariablen. Die entsprechende Datenstruktur können Sie mit type() ermitteln.

Für die Berechnung der Inversen bietet SymPy zwei Möglichkeiten: die Lösung mit der Methode A.inv()*b oder mit der Potenznotation A**-1*b (Zeile 08). Im Lösungsvektor L1 wird eine eindimensionale Matrix gespeichert. Entsprechend muss mit L1[i] auf die einzelnen Lösungen zugegriffen werden (Zeile 16). Eine Typabfrage mit type(L1) liefert die Datenstruktur <class 'sympy.matrices.dense.MutableDenseMatrix'>.

In Zeile 10 löst die Methode A.gauss_jordan_solve(b) das LGS mit dem Gauß-Jordan-Verfahren. Diese Methode greift über den Punktoperator . auf die Koeffizientenmatrix A zu. Der Inhomogenitätsvektor b wird als Parameter übergeben. Im Lösungsvektor L2 wird eine zweidimensionale Matrix gespeichert. Entsprechend muss daher mit L2[i][j] auf die einzelnen Lösungen zugegriffen werden (Zeile 19). Eine Typabfrage mit type(L2) liefert die Datenstruktur <class 'tuple'>.

In Zeile 11 löst die Methode linsolve([A,b],(x1,x2,x3)) das LGS. Die Koeffizientenmatrix A und der Inhomogenitätsvektor b müssen ebenso wie die Unbekannten x1, x2 und x3 durch Klammern zusammengefasst werden. Sie können in beiden Fällen eckige oder runde Klammern verwenden. Im Lösungsvektor L3 wird ein FiniteSet gespeichert. Bei dieser Datenstruktur können mit L3.args[i][j] die einzelnen numerischen Lösungen ermittelt werden (Zeile 22). Eine Typabfrage mit type(L3) liefert die Datenstruktur <class 'sympy.sets.sets.FiniteSet'>.

Die Methode solve_linear_system()

Mit Listing 6.12 können Sie lineare Gleichungssysteme mit der Sympy-Methode solve_linear_system() lösen. Diese Methode erwartet die Übergabe der erweiterten Koeffizientenmatrix.

01  #12_sympy_lgs3.py
02  from sympy import *
03  x1, x2, x3 = symbols('x1 x2 x3')
04  b = Matrix([10,14,16])
05  A = Matrix([[3, 2, 1],
06              [1, 2, 3],
07              [2, 1, 4]])
08  Ab=A.col_insert(3,b)
09  L=solve_linear_system(Ab, x1, x2, x3)
10  print("Type von solve_linear_system()\n",type(L))
11  print("Lösungsvektor\n",L)
12  print(L.keys())
13  print(L.values())
14  print("x1 =",L[x1])
15  print("x2 =",L[x2])
16  print("x3 =",L[x3])

Listing 6.12     Lösung eines linearen Gleichungssystems mit der SymPy-Methode »solve_linear_system()«

Ausgabe

Type von solve_linear_system()
 <class 'dict'>
Lösungsvektor
 {x1: 1, x2: 2, x3: 3}
dict_keys([x1, x2, x3])
dict_values([1, 2, 3])
x1 = 1
x2 = 2
x3 = 3

Analyse

In Zeile 04 wird der Inhomogenitätsvektor b definiert. In den Zeilen 05 bis 07 steht die Koeffizientenmatrix A.

In Zeile 08 erzeugt die SymPy-Methode A.col_insert(3,b) die erweiterte Koeffizientenmatrix Ab. In Zeile 09 löst die Methode solve_linear_system() das lineare Gleichungssystem. Als erster Parameter muss die erweiterte Koeffizientenmatrix Ab übergeben werden. Danach erfolgt die Übergabe der symbolischen Variablen für die Unbekannten.

Der Lösungsvektor L wird als Dictionary ausgegeben. Deshalb können Sie sich die numerischen Lösungen mit L[key] ausgeben lassen, wobei Sie für key die Symbole x1, x2 oder x3 einsetzen.