Ing Mathematik: Numerisches Lösen von linearen Gleichungssystemen

Einleitung

Lineare Gleichungssysteme kommen in der Technik häufig vor. In der Mathematik lernt man in Einführungsvorlesungen z. B. das Gauß-Verfahren als händisches Verfahren zur Lösung linearer Gleichungssysteme kennen, dann aber i. A. für Matrizen mit nur wenigen Zeilen/Spalten. Nun treten in der Praxis (z.B. beim Finite-Elemente-Verfahren) aber Gleichungssysteme mit Tausenden Zeilen auf. Diese händisch zu lösen ist unmöglich. Auch mit einem "computerisierten" Gauß-Verfahren stößt man an die Grenzen von Speicherplatz und Rechenzeit. Aus diesem Grunde wurden zur numerischen Lösung linearer Gleichungssysteme eine Vielzahl von verfeinerten und spezialisierten Verfahren entwickelt. Dies ist nach wie vor ein aktives Feld von mathematischer Forschung. Das nachfolgende Bild gibt einen ersten Eindruck von der Vielfalt an Verfahren (es sind aber nicht alle verfügbaren Verfahren dargestellt).

Als Ingenieur muss man jetzt normalerweise nicht jedes einzelne Verfahren bis ins Detail kennen (dafür gibt es spezialisierte Mathematiker). Zumal die bereits vorgestellten Computersprachen (Octave, Python etc.) schon Funktionen zur Lösung linearer Gleichungssysteme bereitstellen. Mit diesen wird man möglicherweise schon das Auslangen finden. Wenn nicht, dann benötigt man tiefergehende Kenntnisse.

Lineare Gleichungssysteme bei Wikipedia

Direkte vs. iterative Verfahren

Gegeben sei ein lineares Gleichungssystem $Ax=b$ . $A$ heiße Koeffizientenmatrix, $x$ ist der Lösungsvektor und $b$ die bekannte rechte Seite des Gleichungssystems. n sei die Anzahl der Unbekannten des Gleichungssystems.

Ein direktes Verfahren liefert nach einer endlichen Zahl von Rechenschritten den Lösungsvektor. Ein iteratives Verfahren verbessert schrittweise eine Anfangsnäherung. Direkte Verfahren werden bis zu ca. n = 10000 eingesetzt. Iterative Verfahren eignen sich insbesondere für Gleichungssysteme mit n $\gg$ 10000. (Quelle: [1])

Lineare Gleichungssysteme mit Python und SciPy lösen

Das Python-Paket SciPy stellt zahlreiche Solver zur Verfügung. Nachfolgend seien diesbezügliche Beispiele dargestellt. Details zu den einzelnen Solvern folgen später.

Der "Standardsolver"

Quellcode:

import scipy.linalg

A = [[  3,   0,  -2,  0,  0],
     [  0,   6,   0,  0,  0],
     [  0,   0, -12,  1,  0],
     [  0,   0,   1,  0,  0],
     [  0,   0,   0, 10,  7]]

b = [5, -4, 2, 7, 1]

x = scipy.linalg.solve(A, b)
print(x)

Ausgabe:

[   6.33333333   -0.66666667    7.           86.         -122.71428571]

Einige Krylov-Unterraum-Solver

Quellcode:

import numpy as np
from scipy.sparse import csc_matrix
from scipy.sparse.linalg import bicg, gmres, cgs

A = csc_matrix([[  3,   0,  -2,  0,  0],
                [  0,   6,   0,  0,  0],
                [  0,   0, -12,  1,  0],
                [  0,   0,   1,  0,  0],
                [  0,   0,   0, 10,  7]])

b = np.array([5, -4, 2, 7, 1])

x1, exitCode1 = bicg(A, b)
x2, exitCode2 = gmres(A, b)
x3, exitCode3 = cgs(A, b)

print(x1)
print(x2)
print(x3)

Ausgabe:

[   6.33333333   -0.66666667    7.           86.         -122.71428571]
[   6.33333333   -0.66666667    7.           86.         -122.71428571]
[   6.33333333   -0.66666667    7.           86.         -122.71428571]

Das csc in der Funktion csc_matrix steht ausgeschrieben für Compressed Sparse Column. Details zu dieser Funktion siehe z.B. in der SciPy-Dokumentation [2]. Generell ist es empfehlenswert diese Dokumentation bei spezifischen Fragen zu bemühen.

Direkte Verfahren

Das gaußsche Eliminationsverfahren

Ohne Pivotierung

Hier sei ein Programm zu diesem Thema in der Sprache Python in der einfachsten Variante gezeigt. Generell ist zu sagen, dass man die von Python und den zugehörigen Bibliotheken (z.B. SciPy) zur Verfügung gestellten Routinen beim konkreten Einsatz bevorzugt verwenden soll (Gründe: Fehlerfreiheit und Performance). Eigene Implementierungen sind nur für Lehr-, Lern- oder Testzwecke sinnvoll. Oder auch für Verfahren und Weiterentwicklungen, für die es keine SciPy-Prozeduren gibt.

Python-Code:

import numpy as np 

# Gauss-Solver
def gauss(A, b):
  n = b.shape[0]
  x = np.zeros((n,1))

  for p in np.arange(0, n):
    for i in np.arange(p+1, n):
      c_ip = A[i,p] / A[p,p]
      b[i] -= c_ip * b[p]
      for k in np.arange(p, n):
        A[i,k] -=  c_ip * A[p,k]

  for i in np.arange(n-1, -1, -1):
    x[i] = b[i]

    if i<n:
      for k in np.arange(i+1, n):
        x[i] -= A[i,k] * x[k]
    x[i] /= A[i,i]

  return x

# Hauptprogramm
A = np.array([[ 3.,  1.,  2., 1.,  1.],
              [-1.,  2., -1., 3.,  2.],
              [ 0.,  2.,  4., 3.,  9.],
              [-4., -4.,  3., 5., -1.],
              [-5.,  3.,  1., 0., -7.]])

b = np.array([5., -4., -2., 0., -1.])

x = gauss(A, b)
print(x)

Ausgabe:

[[ 1.22537669]
 [-0.21141951]
 [ 1.08929421]
 [ 0.02410785]
 [-0.66740682]]

Im Buch Burg, Haf, Wille, Meister: Höhere Mathematik für Ingenieure, Band II: Lineare Algebra. 7. Auflage, Springer Vieweg, 2012, Seite 91 ist dieser Algorithmus in der Programmiersprache MATLAB implementiert. MATLAB ist weitgehend kompatibel mit GNU Octave. Das dort abgedruckte Programm funktioniert somit auch mit Octave (bei Interesse siehe man im genannten Buch).

Dieses Verfahren ohne Pivotierung versagt, wenn im Laufe der Berechnung ein Diagonalelement $a_{kk}=0$ auftritt. Dann kann das Gauß-Verfahren mit Pivotierung (aus dem französischen, pivot, Dreh- oder Angelpunkt) eine Lösungsmöglichkeit darstellen.

Mit Pivotierung

Bricht obiges Verfahren ab, weil ein Diagonalelement $a_{kk}=0$ ist, so kann dies mit einer Pivotierung behoben werden. Vorausgesetzt natürlich die Koeffizientenmatrix ist regulär. Dabei werden Zeilen vertauscht (oder auch Spalten). Man spricht dann von Zeilen-, Spalten- oder Totalpivotierung. Eine Pivotierung kann auch im Sinne der Minimierung von Rundungsfehlern sinnvoll sein.

Das Cholesky-Verfahren

Für positiv definite Matrizen kann auch das Verfahren von Cholesky eingesetzt werden. Positiv definit heißen Matrizen für die Folgendes gilt: $x^{T}Ax>0$ für alle $x\in \mathbb {R} ^{n},x\neq 0$ und $A$ reell und symmetrisch.

Festgestellt werden kann die positive Definitheit mit dem Kriterium von Hadamard (Burg, Haf, Wille, Meister: Seite 194f und 206ff) bzw. dem Satz von Sylvester (siehe Riemer, Seemann, Wauer, Wedig: Mathematische Methoden der Technischen Mechanik. 3. Auflage, Springer, 2019, Seite 10). Beide Sätze oder Kriterien sagen das Gleiche aus, nämlich: Sind alle Hauptminoren der Matrix positiv, so ist die Matrix positiv definit. Es gilt auch die Regel: Eine quadratisch positiv definite Matrix ist immer regulär. Zum Begriff Hauptminor wird ein Beispiel geliefert:

$A={\begin{pmatrix}1&2&3\\4&5&6\\7&8&9\end{pmatrix}}$

Die Hauptminoren sind folgende Unterdeterminanten:

$\det(A_{1})=1;\quad \det(A_{2})={\begin{vmatrix}1&2\\4&5\\\end{vmatrix}}=-3;\quad \det(A_{3})={\begin{vmatrix}1&2&3\\4&5&6\\7&8&9\end{vmatrix}}=0.$

Die Matrix $A$ in diesem Beispiel ist somit nicht positiv definit (weil nicht alle Hauptminoren positiv sind)! Im konkreten Fall ist $A$ sogar singulär, weil $\det(A)=0$ . Gleichwertig zum Kriterium von Hadamard ist folgendes: Alle Eigenwerte von $A$ sind positiv.

Aber nun wieder zurück zum Cholesky-Verfahren. Die positiv definite Matrix $A$ kann zerlegt werden in $A=LL^{T}$ . Wobei $L$ eine links-untere Dreiecksmatrix ist.

Skripten im WWW

Einfach nach den Stichwörtern "numerische Mathematik Skript" o.ä. googeln. Es finden sich jede Menge Skripten im PDF-Format zu diesem Thema im Internet.

Gedruckte Werke (auszugsweise)

Burg, Haf, Wille, Meister: Höhere Mathematik für Ingenieure, Band II: Lineare Algebra. 7. Auflage, Springer Vieweg, 2012, ISBN 978-3-8348-1853-9
Hanke-Bourgeois: Grundlagen der Numerischen Mathematik und des Wissenschaftlichen Rechnens. 3. Aufl., Vieweg+Teubner, 2009, ISBN 978-3-8348-0708-3
Knorrenschild: Numerische Mathematik, Eine beispielorientierte Einführung. 6. Aufl., Hanser, 2017, ISBN 978-3-446-45161-2
Meister: Numerik linearer Gleichungssysteme. 4. Aufl., Vieweg+Teubner, 2011, ISBN 978-3-8348-1550-7

Zurück zu Verisonsverwaltung mit git |

Hoch zu Gesamtinhaltsverzeichnis |