Am Ende von Kapitel 3 stand ein Gleichungssystem an der Tafel — das reduzierte 3×3-System der Federkette. Wir haben es von Hand gelöst, durch geschicktes Einsetzen, und drei Zahlen geerntet: \(u_2 = 0{,}5\), \(u_3 = 1{,}0\), \(u_4 = 1{,}5\ \mathrm{mm}\). Das war überschaubar, weil es nur drei Unbekannte waren.
Aber die Botschaft der Federkette war ja gerade: Der Weg bleibt gleich, egal wie lang die Kette wird. Das Netz, mit dem wir am Ende dieses Buches den Träger im Brandfall durchrechnen, hat rund 400 Unbekannte. Ein Handy-Crashtest bei einem Hersteller hat leicht zehn Millionen. Niemand setzt zehn Millionen Gleichungen von Hand ineinander ein. Die Frage dieses Kapitels lautet darum:
Wie bringt man einer Maschine bei, ein Gleichungssystem zu lösen — und warum dauert das bei zehn Millionen Unbekannten nicht Jahre?
Wir bauen dazu drei Löser, die den Rest des Buches antreiben werden. Alle drei bekommen zur Prüfung dasselbe reduzierte Federketten-System vorgelegt, und alle drei müssen am Ende dieselben drei Zahlen liefern wie unsere Handrechnung: \(0{,}5\), \(1{,}0\), \(1{,}5\ \mathrm{mm}\). Die Handrechnung aus Kapitel 3 ist von jetzt an unsere Probe.
Lernziele
Nach diesem Kapitel kannst du …
… ein 3×3-System mit Gauß-Elimination von Hand lösen — leerräumen, dann rückwärts einsetzen,
… erklären, warum die Federketten-Matrix fast nur Nullen enthält (dünn besetzt, Bandmatrix) und warum das Rechenzeit spart,
… den Thomas-Algorithmus als „Gauß, der nur das Band anfasst” anwenden,
… Gauß-Seidel als Entspannungsverfahren erklären und beobachten, wie sich die Lösung Runde um Runde einpendelt,
… abschätzen, warum volle Elimination bei großen Systemen unbezahlbar wird — mit konkreten Zahlen, ohne Fachbegriffe.
WarnungNaheliegende Vermutung
Vermutung:„Doppelt so viele Unbekannte — doppelt so lange Rechenzeit. Zehnmal so viele, zehnmal so lange.”
Warum sie naheliegt: Das ist die Erfahrung aus dem Alltag: doppelt so viele Kartoffeln schälen dauert doppelt so lange. Rechenarbeit müsste sich doch genauso verhalten.
Was stattdessen stimmt: Bei der vollen Gauß-Elimination verachtfacht sich die Arbeit, wenn sich die Zahl der Unbekannten nur verdoppelt — nicht mal zwei, sondern mal zwei hoch drei. Von 10 auf 20 Unbekannte ist es kaum zu merken; von einer Million auf zwei Millionen wird aus einem Tag mehr als eine Woche. Genau deshalb reicht es nicht, irgend einen Löser zu haben. Der Ausweg ist, die vielen Nullen der FEM-Matrix auszunutzen: Beim Thomas-Algorithmus ist doppelt so groß dann wirklich nur doppelt so teuer. Diesen Unterschied — Fluch und Rettung — rechnen wir in diesem Kapitel mit echten Zahlen aus.
4.1 Was „lösen” überhaupt heißt
Bevor wir eine Maschine losschicken, halten wir fest, was sie eigentlich suchen soll. Ein Gleichungssystem zu lösen heißt: Zahlen finden, die alle Gleichungen gleichzeitig erfüllen. Nicht die eine oder die andere — alle auf einmal.
Am kleinsten Fall sieht man es am besten. Nimm zwei Gleichungen mit zwei Unbekannten \(x\) und \(y\):
\[
\begin{aligned}
x + y &= 3,\\
x - y &= 1.
\end{aligned}
\]
Jede dieser Gleichungen beschreibt eine Gerade in einem \(x\)-\(y\)-Diagramm — eine ganze Straße von Punkten, die sie erfüllen. Die erste erfüllen alle Punkte mit \(x + y = 3\), die zweite alle mit \(x - y = 1\). Gesucht ist der eine Punkt, der auf beiden Geraden liegt: ihr Schnittpunkt.
Code
import matplotlib.pyplot as pltxs = [-1.0, 4.0]gerade1 = [3.0- x for x in xs] # x + y = 3 -> y = 3 - xgerade2 = [x -1.0for x in xs] # x - y = 1 -> y = x - 1fig, ax = plt.subplots(figsize=(6.2, 5.0))ax.plot(xs, gerade1, color="tab:blue", lw=2.0, label="x + y = 3")ax.plot(xs, gerade2, color="tab:orange", lw=2.0, label="x − y = 1")ax.plot(2.0, 1.0, "o", ms=13, color="tab:green", markeredgecolor="black", zorder=5)ax.annotate("Lösung (2, 1)", xy=(2.0, 1.0), xytext=(2.5, 2.1), fontsize=11, color="tab:green", arrowprops=dict(arrowstyle="-|>", color="tab:green", lw=1.6))ax.axhline(0, color="0.7", lw=0.8)ax.axvline(0, color="0.7", lw=0.8)ax.set_xlabel("x")ax.set_ylabel("y")ax.set_xlim(-1.0, 4.0)ax.set_ylim(-2.0, 4.0)ax.set_aspect("equal")ax.legend(loc="upper left")ax.set_title("Zwei Geraden, ein Schnittpunkt")plt.tight_layout()plt.show()
Abbildung 4.1: Lösen heißt schneiden: Jede der beiden Gleichungen ist eine Gerade. Die erste (blau) sammelt alle Punkte mit x + y = 3, die zweite (orange) alle mit x − y = 1. Nur ein einziger Punkt liegt auf beiden — der Schnittpunkt bei x = 2, y = 1. Das ist die Lösung des Systems.
Was man hier sieht: Zwei Geraden, ein Schnittpunkt. Am Punkt \((2, 1)\) gilt \(2 + 1 = 3\)und\(2 - 1 = 1\) — beide Gleichungen stimmen zugleich. Setz die Zahlen ein und probiere es nach. Genau das ist eine Lösung: ein Satz von Werten, an dem alle Gleichungen aufgehen.
Bei drei Unbekannten — wie im reduzierten Federketten-System — wären es drei „Flächen” im Raum, die sich in einem Punkt treffen; bei vieren kann man es nicht mehr zeichnen. Deshalb rechnen wir statt zu zeichnen, und zwar so stur und schematisch, dass es auch eine Maschine tun kann. Das erste Verfahren dafür ist über zweihundert Jahre alt und trägt den Namen des Mannes, der es in eine feste Form goss: Carl Friedrich Gauß.
4.2 Gauß-Elimination: erst leerräumen, dann ernten
Die Idee der Gauß-Elimination ist bestechend einfach. Ein System mit drei Unbekannten ist schwer, weil in jeder Zeile alle drei zugleich vorkommen. Wäre in der letzten Zeile nur noch eine Unbekannte, könnten wir sie sofort ausrechnen. Stünden in der vorletzten nur noch zwei, davon eine schon bekannt, ginge auch die. Und so weiter nach oben.
Das Verfahren stellt genau das her. Es hat zwei Phasen:
Leerräumen (Vorwärtselimination): Wir machen mit geschickten Zeilen-Kombinationen alle Zahlen unter der Diagonale zu null. Übrig bleibt eine Treppe — eine Dreiecksform, in der die letzte Zeile nur noch eine Unbekannte enthält.
Ernten (Rückwärtseinsetzen): Von unten nach oben rechnen wir eine Unbekannte nach der anderen aus, jede aus den schon bekannten darüber.
Führen wir es an unserem reduzierten Federketten-System vor. Es lautet, Zeile für Zeile als Gleichung geschrieben (Einheiten N und mm lassen wir beim Rechnen weg):
Bevor du weiterrechnest: Wir kennen die Antwort schon aus Kapitel 3 (\(0{,}5\), \(1{,}0\), \(1{,}5\)). Diesmal geht es um den Weg. Der erste Schritt wird die \(-100\) am Anfang von Zeile (II) beseitigen, indem er ein Vielfaches von Zeile (I) dazuzählt. Welches Vielfache? Überlege: Was muss man mit den \(200\) aus Zeile (I) tun, damit gerade \(+100\) herauskommt, das die \(-100\) ausgleicht? Leg dich fest.
4.2.1 Phase 1, Schritt 1: die erste Spalte leerräumen
Wir wollen in Zeile (II) das \(-100\,u_2\) loswerden. In Zeile (I) steht an derselben Stelle \(200\,u_2\). Addieren wir die Hälfte von Zeile (I) zu Zeile (II), so wird aus \(-100 + \tfrac{1}{2}\cdot 200 = -100 + 100 = 0\) — genau null. Die halbe Zeile (I) ist unser Werkzeug (das Vielfache aus dem Vorhersage-Punkt: der Faktor \(\tfrac{1}{2}\)).
Rechnen wir Zeile (II) \(+\ \tfrac{1}{2}\cdot\) Zeile (I) Glied für Glied aus:
Die erste Spalte ist unter der Diagonale leer. \(u_2\) kommt nur noch in der ersten Zeile vor.
4.2.2 Phase 1, Schritt 2: die zweite Spalte leerräumen
Jetzt dasselbe eine Etage tiefer: In Zeile (III) soll das \(-100\,u_3\) verschwinden. Darüber, in Zeile (II’), steht jetzt \(150\,u_3\). Diesmal brauchen wir nicht die Hälfte, sondern \(\tfrac{100}{150} = \tfrac{2}{3}\) von Zeile (II’), damit \(-100 + \tfrac{2}{3}\cdot 150 = -100 + 100 = 0\) wird.
Abbildung 4.2: Die Vorwärtselimination als Schrittbild. Links das volle System, in der Mitte nach dem ersten Schritt (die erste Spalte unter der Diagonale ist geräumt), rechts nach dem zweiten (auch die zweite Spalte). Die grün schraffierte Treppe ist die Dreiecksform: unter ihr steht nur noch null. Die rot umrandeten Felder sind die, die im jeweiligen Schritt gerade zu null gemacht wurden.
Was man hier sieht: Wie das System von links nach rechts „leergeräumt” wird. In jedem Schritt fällt ein Feld unter der Diagonale auf null (rot umrandet), bis nur noch die grüne Dreieckstreppe Zahlen trägt. Diese Form ist das Ziel der ersten Phase — an ihr kann man von unten nach oben ernten.
4.2.3 Phase 2: die Ernte (Rückwärtseinsetzen)
Jetzt kommt der schöne Teil. Wir lesen die Dreiecksform von unten nach oben.
Aus Zeile (III’) — der einzelnen Unbekannten — folgt sofort
Geerntet: \(u_2 = 0{,}5\), \(u_3 = 1{,}0\), \(u_4 = 1{,}5\ \mathrm{mm}\) — Zahl für Zahl dieselbe Lösung wie die Handrechnung aus Kapitel 3. Die Probe stimmt. Wir haben sie diesmal nicht durch pfiffiges Einsetzen gefunden, sondern durch ein stures Schema, das man Zeile für Zeile abarbeitet, ohne nachzudenken — und genau darum kann es eine Maschine.
4.2.4 Dasselbe in Python
Die Prosa hat erklärt, was passiert; jetzt macht der Code die Schritte anfassbar. Das folgende Programm ist die Vorwärtselimination in reinem Python — dieselbe Doppelschleife, die wir eben auf Papier ausgeführt haben, mit denselben Namen: für jede Spalte den Pivot nehmen, darunter jede Zeile um ein Vielfaches der Pivotzeile verringern. Nach jedem Schritt druckt es die Matrix, damit du das Leerräumen mitverfolgen kannst (die Nullen als Punkte, wie beim Matrix-Drucker aus Kapitel 3).
Code
def zeige_elimination(matrix, rechte_seite):"""Fuehrt die Vorwaertselimination aus und druckt jeden Schritt. matrix: die Systemmatrix als verschachtelte Liste. rechte_seite: Liste der rechten Seiten. Rueckgabe: nichts; die Funktion druckt nur mit. Reine Vorfuehrung: Spalte fuer Spalte wird unter der Diagonale leergeraeumt, nach jeder Spalte wird die Matrix gezeigt. """ n =len(matrix)# Auf Kopien arbeiten, damit das Original erhalten bleibt. werk = []for zeile in matrix: werk.append(list(zeile)) seite =list(rechte_seite)for spalte inrange(n): pivot = werk[spalte][spalte]for zeile inrange(spalte +1, n): faktor = werk[zeile][spalte] / pivotfor k inrange(spalte, n): werk[zeile][k] = werk[zeile][k] - faktor * werk[spalte][k] seite[zeile] = seite[zeile] - faktor * seite[spalte]print(f"Nach dem Leerraeumen von Spalte {spalte +1}:") drucke(werk, seite)def drucke(matrix, seite):"""Druckt Matrix und rechte Seite; Nullen werden zu Punkten."""for zeile inrange(len(matrix)):for wert in matrix[zeile]:ifabs(wert) <0.0001:print(f"{'.':>8}", end="")else:print(f"{wert:8.2f}", end="")print(f" | {seite[zeile]:7.2f}")print()system = [[200.0, -100.0, 0.0], [-100.0, 200.0, -100.0], [0.0, -100.0, 100.0]]kraefte = [0.0, 0.0, 50.0]zeige_elimination(system, kraefte)
Nach dem Leerraeumen von Spalte 1:
200.00 -100.00 . | 0.00
. 150.00 -100.00 | 0.00
. -100.00 100.00 | 50.00
Nach dem Leerraeumen von Spalte 2:
200.00 -100.00 . | 0.00
. 150.00 -100.00 | 0.00
. . 33.33 | 50.00
Nach dem Leerraeumen von Spalte 3:
200.00 -100.00 . | 0.00
. 150.00 -100.00 | 0.00
. . 33.33 | 50.00
Interpretation der Ausgabe: Genau unsere zwei Papier-Schritte. Nach Spalte 1 ist die erste Spalte unter der Diagonale leer und in der Mitte steht die \(150\); nach Spalte 2 auch die zweite, und unten rechts steht die \(33{,}33\). Der Code hat nichts anderes getan als wir — nur schneller und ohne sich zu verrechnen. Das Rückwärtseinsetzen und die drei Verfahren dieses Kapitels stecken gebündelt im wiederverwendbaren Modul programme/gemeinsam/loeser.py, das die späteren Kapitel importieren.
4.3 Warum die FEM-Matrix fast nur aus Nullen besteht
Die Gauß-Elimination funktioniert für jedes lösbare System. Aber sie ist verschwenderisch, wenn das System viele Nullen enthält — und die FEM-Matrizen enthalten fast nur Nullen. Erinnern wir uns an den Matrix-Drucker aus Kapitel 3: Bei fünf Federn waren von den 36 Feldern der 6×6-Matrix nur 16 belegt, der Rest war leer. Drucken wir zum Vergleich die Matrix einer längeren Kette und daneben, wie eine „volle” Matrix gleicher Größe aussähe.
Code
import matplotlib.pyplot as pltdef zeichne_muster(ax, belegt, titel): n =len(belegt)for r inrange(n):for c inrange(n): ax.add_patch(plt.Rectangle((c, n -1- r), 1, 1, facecolor="tab:blue"if belegt[r][c] else"white", edgecolor="0.85", lw=0.6)) ax.set_xlim(0, n) ax.set_ylim(0, n) ax.set_aspect("equal") ax.axis("off") ax.set_title(titel, fontsize=10)n =10voll = [[Truefor _ inrange(n)] for _ inrange(n)]tri = [[abs(r - c) <=1for c inrange(n)] for r inrange(n)]# 2D-Netzmatrix: Tridiagonal plus zwei Nebenbaender im Abstand einer# „Netzbreite" (hier 3) -- schematisch fuer den Ausblick.breite =3band2d = [[(abs(r - c) <=1) or (abs(r - c) == breite) for c inrange(n)]for r inrange(n)]fig, achsen = plt.subplots(1, 3, figsize=(10.5, 3.9))zeichne_muster(achsen[0], voll, "voll besetzt\n(100 Zahlen)")zeichne_muster(achsen[1], tri, "Federkette / 1D-Träger\n(28 Zahlen)")zeichne_muster(achsen[2], band2d, "2D-Netz (Ausblick)\n(dünn, breiteres Band)")plt.tight_layout()plt.show()
Abbildung 4.3: Drei Besetzungsmuster im Vergleich (blau = hier steht eine Zahl, weiß = null). Links eine voll besetzte 10×10-Matrix: jedes Feld trägt eine Zahl. In der Mitte die Federkette (tridiagonal): nur die Diagonale und ihre beiden Nachbarn — ein schmales Band. Rechts, als Ausblick, eine 2D-Netzmatrix: auch dünn besetzt, aber mit einem breiteren Band und zwei zusätzlichen Nebenbändern. Je dünner das Band, desto billiger die Lösung.
Was man hier sieht: Dieselbe Größe, drei Welten. Die volle Matrix (links) hat in jedem der 100 Felder eine Zahl. Die Federkette (Mitte) hat nur 28 — die Diagonale und je eine Nebendiagonale; alles andere ist leer. Man nennt eine solche Matrix dünn besetzt: Der allergrößte Teil ihrer Felder ist null. Und weil die wenigen Zahlen sich als schmales Band um die Diagonale drängen, heißt sie eine Bandmatrix. Beides hat denselben physikalischen Grund, den wir in Kapitel 3 gesehen haben: Jeder Knoten hängt nur mit seinen unmittelbaren Nachbarn zusammen, mit keinem Knoten weiter weg. Rechts deutet sich an, wie es in 2D aussieht (ab Kapitel 10): immer noch dünn, aber mit einem etwas breiteren Band.
Diese Leere ist bares Geld. Ein Löser, der die Nullen erst gar nicht ansieht, spart fast die gesamte Arbeit. Genau das tut der nächste Löser.
4.4 Der Thomas-Algorithmus: Gauß, der nur das Band anfasst
Schau noch einmal auf unsere Gauß-Elimination von eben. Beim Leerräumen der ersten Spalte mussten wir nur Zeile (II) anfassen — Zeile (III) hatte an der ersten Stelle ohnehin schon eine Null, dort war nichts zu eliminieren. Und beim Leerräumen der zweiten Spalte gab es überhaupt nur eine Zeile darunter. Die vielen Nullen haben uns Arbeit geschenkt, ohne dass wir es geplant hatten.
Der Thomas-Algorithmus macht aus diesem Geschenk ein System. Er ist nichts anderes als die Gauß-Elimination, aber er weiß von vornherein, dass die Matrix tridiagonal ist — dass also außerhalb des schmalen Bandes nur Nullen stehen. Darum fasst er in jeder Zeile nur drei Zahlen an: die auf der Diagonale und ihre beiden Nachbarn. Alles andere überspringt er, weil es sowieso null ist und null bleibt.
Man beschreibt ein tridiagonales System darum gar nicht mehr als volles Quadrat, sondern nur durch seine drei Diagonalen — drei Listen:
die untere Nebendiagonale (die Zahlen links der Diagonale),
die Hauptdiagonale,
die obere Nebendiagonale (die Zahlen rechts der Diagonale).
Für unser reduziertes Federketten-System sind das:
Und jetzt kommt das Schöne: Wenn der Thomas-Algorithmus dieses Band abarbeitet, tut er genau dieselben Rechnungen wie unsere Gauß-Handrechnung vorhin — dieselben Faktoren \(\tfrac{1}{2}\) und \(\tfrac{2}{3}\), dieselben Zwischenzahlen \(150\) und \(\tfrac{100}{3}\). Er läuft einmal von oben nach unten durch das Band (die Vorwärtswelle, das Leerräumen) und einmal von unten nach oben zurück (die Rückwärtswelle, das Ernten). Nur berührt er nie ein Feld außerhalb des Bandes. Bei drei Unbekannten fällt der Unterschied nicht auf; bei tausend spart er fast die gesamte Arbeit.
HinweisGeschichte: der Namensgeber des Thomas-Algorithmus
Der Algorithmus trägt den Namen von Llewellyn Hilleth Thomas (1903–1992), einem britischen Physiker. Bekannt wurde er zunächst durch einen ganz anderen Beitrag — die „Thomas-Präzession” in der Atomphysik der 1920er Jahre. Das Verfahren zum schnellen Lösen tridiagonaler Systeme beschrieb er 1949 in einem Bericht des Watson Scientific Computing Laboratory an der Columbia University (Thomas 1949). Dass ausgerechnet solche Systeme so oft auftauchen — bei Wärmeleitung, bei Schwingungen, überall dort, wo nur Nachbarn miteinander reden —, macht seinen kleinen, schnellen Algorithmus bis heute zu einem der meistbenutzten der Numerik.
Den Thomas-Löser haben wir als wiederverwendbaren Baustein in programme/gemeinsam/loeser.py abgelegt. Sein Kern ist kurz genug, um ihn ganz zu zeigen — die Vorwärtswelle, die das Band leerräumt, und die Rückwärtswelle, die erntet:
Code
def thomas(unter, diagonale, ober, rechte_seite):"""Loest ein tridiagonales System, ohne die Nullen anzufassen. unter: Nebendiagonale links (unter[0] wird nicht benutzt). diagonale: die Hauptdiagonale (n Zahlen). ober: Nebendiagonale rechts (ober[n-1] wird nicht benutzt). rechte_seite: die n rechten Seiten. Rueckgabe: Liste der n Loesungswerte. """ n =len(diagonale) ober_neu =list(ober) seite_neu =list(rechte_seite)# Vorwaertswelle: von oben nach unten das Band leerraeumen. ober_neu[0] = ober[0] / diagonale[0] seite_neu[0] = rechte_seite[0] / diagonale[0]for i inrange(1, n): nenner = diagonale[i] - unter[i] * ober_neu[i -1]if i < n -1: ober_neu[i] = ober[i] / nenner seite_neu[i] = (rechte_seite[i] - unter[i] * seite_neu[i -1]) / nenner# Rueckwaertswelle: von unten nach oben ernten. loesung = [0.0] * n loesung[n -1] = seite_neu[n -1]for i inrange(n -2, -1, -1): loesung[i] = seite_neu[i] - ober_neu[i] * loesung[i +1]return loesungunter = [0.0, -100.0, -100.0] # unter[0] gibt es nichtdiagonale = [200.0, 200.0, 100.0]ober = [-100.0, -100.0, 0.0] # ober[2] gibt es nichtrechte_seite = [0.0, 0.0, 50.0]loesung = thomas(unter, diagonale, ober, rechte_seite)gerundet = []for wert in loesung: gerundet.append(round(wert, 3))print("Thomas-Lösung:", gerundet, "mm")
Thomas-Lösung: [0.5, 1.0, 1.5] mm
Interpretation der Ausgabe:\(0{,}5\), \(1{,}0\), \(1{,}5\) — wieder dieselbe Probe. Der Thomas-Löser hat nie ein Feld außerhalb des Bandes berührt und liefert doch aufs Tausendstel dasselbe wie die volle Gauß-Elimination.
4.4.1 Der Thomas-Lauf als Welle durch den Träger
Weil der Thomas-Algorithmus das Band strikt der Reihe nach abarbeitet, kann man ihm zusehen wie einer Welle, die durch den Träger läuft. Die folgende Animation zeigt es an einer längeren Kette: Erst wandert die Vorwärtswelle von links nach rechts und räumt Knoten für Knoten leer (das Leerräumen), dann kehrt die Rückwärtswelle von rechts nach links zurück und trägt Knoten für Knoten den fertigen Verschiebungswert ein (das Ernten).
Code
import matplotlib.pyplot as pltfrom matplotlib.animation import FuncAnimationfrom IPython.display import HTMLanzahl =9# neun Knoten (acht Federn)loesung = [0.5* i for i inrange(anzahl)] # 0, 0.5, 1.0, ... 4.0 mmpositionen =list(range(anzahl))fig, ax = plt.subplots(figsize=(9.0, 2.6))def zeichne(frame): ax.clear()# frame 0..anzahl-1: Vorwaertswelle; danach Rueckwaertswelle. vor = frame < anzahlif vor: aktiv = frameelse: aktiv = anzahl -1- (frame - anzahl)for i inrange(anzahl):if vor: fertig = i < aktiv farbe ="tab:green"if i == aktiv else ("0.6"if fertig else"0.85") hoehe =0.0else: geerntet = i >= aktiv farbe ="tab:orange"if i == aktiv else ("tab:blue"if geerntet else"0.6") hoehe = loesung[i] if geerntet else0.0 ax.plot(positionen[i], hoehe, "o", ms=15, color=farbe, markeredgecolor="black", zorder=4)if (not vor) and i >= aktiv: ax.text(positionen[i], hoehe +0.25, "%.1f"% loesung[i], ha="center", fontsize=8) ax.add_patch(plt.Rectangle((-0.6, -0.4), 0.35, 0.8, hatch="///", facecolor="lightsteelblue", edgecolor="black")) titel = ("Vorwärtswelle: Knoten %d leerräumen"% (aktiv +1) if vorelse"Rückwärtswelle: Verschiebung an Knoten %d ernten"% (aktiv +1)) ax.set_title(titel, fontsize=11) ax.set_xlim(-0.8, anzahl -0.2) ax.set_ylim(-0.6, 4.6) ax.set_xlabel("Knoten (Verschiebung in mm, überhöht)") ax.set_yticks([])return []ani = FuncAnimation(fig, zeichne, frames=2* anzahl, interval=450, blit=False)plt.close(fig)HTML(ani.to_jshtml())
Abbildung 4.4: Der Thomas-Lauf an einer Kette aus acht Federn. Grün eingefärbt der Knoten, den die Vorwärtswelle gerade leerräumt (links nach rechts); orange der Knoten, an dem die Rückwärtswelle gerade den fertigen Verschiebungswert einträgt (rechts nach links). Der Algorithmus berührt jeden Knoten genau zweimal — einmal hin, einmal zurück.
Was man sieht: Der Algorithmus berührt jeden Knoten genau zweimal — einmal auf dem Hinweg (leerräumen), einmal auf dem Rückweg (ernten). Zwei Durchläufe durch die Kette, mehr nicht. Verdoppelt man die Kettenlänge, verdoppelt sich schlicht die Länge der beiden Wege: doppelt so viele Knoten, doppelt so viel Arbeit. Genau das Versprechen aus dem Vermutungs-Kasten.
4.5 Gauß-Seidel: die Lösung pendelt sich ein
Die beiden ersten Verfahren rechnen die Lösung in einem Zug aus — sie räumen leer und ernten, fertig. Das dritte Verfahren geht ganz anders vor, und es lohnt sich, weil es später die großen 2D-Netze bändigt und weil es ein Bild liefert, das den Kern der ganzen Physik trifft.
Die Idee heißt Entspannung (englisch relaxation). Stell dir vor, jeder Knoten wüsste seine eigene Gleichung — aber nur seine eigene. Knoten 3 etwa weiß: „Meine Verschiebung sollte in der Mitte zwischen meinen beiden Nachbarn liegen.” Er schaut also auf seine Nachbarn und rückt sich dorthin, wohin seine Gleichung ihn haben will. Dann macht es Knoten 4 genauso, dann Knoten 2, und dann wieder von vorn. Kein Knoten sieht das ganze System; jeder korrigiert nur sich selbst, immer wieder. Und erstaunlicherweise pendelt sich die ganze Kette dabei in die richtige Lösung ein.
Machen wir das konkret. Wir stellen jede Zeile des reduzierten Systems nach ihrer eigenen Unbekannten um — nach der Zahl auf der Diagonale:
Lies die mittlere Zeile als Satz: „\(u_3\) ist der Mittelwert seiner beiden Nachbarn.” Das ist die Entspannungsregel in Reinform. Jetzt raten wir einen Startwert — sagen wir, wir wissen nichts und setzen alle Verschiebungen auf 0 — und wenden die drei Regeln immer wieder an. Dabei benutzen wir jeden neuen Wert sofort weiter, sobald wir ihn haben (das ist der Unterschied zwischen Gauß-Seidel und seinem etwas langsameren Bruder, dem Jacobi-Verfahren).
WichtigVorhersage-Punkt
Bevor du die Tabelle liest: Wir starten bei \(0/0/0\), die richtige Lösung ist \(0{,}5/1{,}0/1{,}5\). In der ersten Runde ändert sich nur \(u_4\) (weil nur seine Gleichung eine Zahl von außen, die \(0{,}5\), mitbringt) — die beiden anderen bleiben zunächst bei null. Was glaubst du: Sind wir nach 5 Runden schon nah dran, oder braucht es eher 20? Und rückt die Lösung von unten heran (immer zu klein) oder pendelt sie über das Ziel hinaus? Leg dich fest.
Rechnen wir die ersten Runden von Hand mit (jede Zahl aus der Regel darüber, immer mit den frischesten Nachbarwerten):
Tabelle 4.1: Die ersten fünf Gauß-Seidel-Runden am reduzierten Federketten-System. Alle Werte kriechen von unten auf ihre Zielwerte 0,5 / 1,0 / 1,5 zu; die „größte Änderung” wird von Runde zu Runde kleiner.
Runde
\(u_2\)
\(u_3\)
\(u_4\)
größte Änderung
Start
0
0
0
—
1
0
0
0,500
0,500
2
0
0,250
0,750
0,250
3
0,125
0,438
0,938
0,188
4
0,219
0,578
1,078
0,141
5
0,289
0,684
1,184
0,105
Man sieht es in Tabelle 4.1: Die Werte nähern sich von unten und pendeln nicht über das Ziel hinaus — jede Runde ein Stück näher, aber die Schritte werden kleiner. Nach fünf Runden sind wir noch nicht da (\(u_2\) steht erst bei \(0{,}29\) statt \(0{,}5\)). Das ist typisch: Gauß-Seidel braucht viele kleine Runden, wo Gauß und Thomas in einem Zug fertig sind. Sein Vorteil liegt woanders — bei den großen 2D-Netzen, wo eine volle Elimination zu teuer wird.
4.5.1 Woher weiß man, wann man aufhören darf?
Anders als Gauß hat Gauß-Seidel kein natürliches Ende — man könnte ewig weiterrechnen und käme der Lösung immer noch ein Härchen näher. Man braucht also ein Abbruchkriterium: eine Regel, die sagt „gut genug ist gut genug”. Die einfachste und ehrlichste ist die letzte Spalte von Tabelle 4.1: Wenn sich in einer ganzen Runde kein Wert mehr um mehr als eine kleine Schwelle ändert — etwa um weniger als \(0{,}001\ \mathrm{mm}\) —, dann bewegt sich nichts mehr nennenswert, und wir hören auf. Man verlangt also nicht die perfekte Lösung (die man ohne Gauß gar nicht kennt), sondern eine, die sich nicht mehr merklich verbessert.
4.5.2 Das Einpendeln als Animation
Die Tabelle zeigt Zahlen; das Bild zeigt die Bewegung. Die folgende Animation lässt die Federkette Runde um Runde entspannen. Sie startet als flache Linie bei null (alle Verschiebungen falsch) und schwingt sich in die gleichmäßige Treppe \(0/0{,}5/1{,}0/1{,}5\) ein — genau die Lösung, die wir in Kapitel 3 von Hand gefunden haben.
Abbildung 4.5: Gauß-Seidel am Werk: Das Verschiebungsprofil der Kette startet flach bei null (grau) und pendelt sich Runde um Runde (blau, immer dunkler) in die endgültige Treppe 0 / 0,5 / 1,0 / 1,5 mm ein (grün gestrichelt das Ziel). Knoten 1 bleibt fest im Auflager. Jede Runde rückt die Kurve ein Stück näher, die Schritte werden kleiner.
Was man sieht: Aus der falschen flachen Linie wächst Runde um Runde die richtige Treppe. Das ist mehr als eine Rechenspielerei — es ist ein Bild für etwas Physikalisches. Später, wenn dieselbe Sorte Gleichungssystem nicht Verschiebungen, sondern Temperaturen trägt (ab Kapitel 7), beschreibt genau dieses Einpendeln, wie sich der Träger im Brandfall auf sein endgültiges Temperaturprofil einschwingt. Die Entspannungsregel „jeder Knoten ist der Mittelwert seiner Nachbarn” ist dann kein Rechentrick mehr, sondern die Physik selbst. Doch das ist die Geschichte von Teil III — hier bleibt es bei der Federkette.
4.5.3 Wie schnell wird der Fehler kleiner?
Man möchte wissen, wie viele Runden Gauß-Seidel braucht. Tragen wir dazu den Fehler — die größte Änderung pro Runde aus Tabelle 4.1 — über der Rundenzahl auf. Weil der Fehler in jeder Runde um denselben Bruchteil schrumpft (grob um ein Viertel), fällt er anfangs schnell und dann immer zäher; eine gewöhnliche Achse würde das schlecht zeigen. Darum benutzen wir für den Fehler eine besondere Achsenteilung.
TippWas ist eine logarithmische Achse?
Auf einer gewöhnlichen Achse haben gleiche Abstände gleiche Differenzen: von 1 bis 2 ist es genauso weit wie von 100 bis 101. Auf einer logarithmischen Achse haben gleiche Abstände gleiche Faktoren: von 1 bis 10 ist es genauso weit wie von 10 bis 100 oder von 0,01 bis 0,1 — jeder gleiche Schritt bedeutet „mal zehn”. Das ist genau richtig für etwas, das in jeder Runde denselben Bruchteil wegnimmt: Auf der Log-Achse wird aus dem zäh abfallenden Bogen eine gerade Linie, und ihre Steilheit sagt unmittelbar, wie schnell das Verfahren ist. Fällt die Gerade steil, halbiert sich der Fehler in wenigen Runden; fällt sie flach, braucht es viele.
Code
import matplotlib.pyplot as plt# Fehler = groesste Aenderung je Runde, aus dem Verfahren selbst.u2, u3, u4 =0.0, 0.0, 0.0runden = []fehler = []for runde inrange(1, 30): alt2, alt3, alt4 = u2, u3, u4 u2 =0.5* u3 u3 =0.5* (u2 + u4) u4 =0.5+ u3 aenderung =max(abs(u2 - alt2), abs(u3 - alt3), abs(u4 - alt4)) runden.append(runde) fehler.append(aenderung)fig, ax = plt.subplots(figsize=(7.4, 4.2))ax.semilogy(runden, fehler, "o-", color="tab:blue", lw=1.8, ms=6)ax.axhline(0.001, color="tab:red", ls="--", lw=1.4)ax.text(0.5, 0.0012, "Abbruchschwelle 0,001", color="tab:red", fontsize=9)ax.set_xlabel("Runde")ax.set_ylabel("größte Änderung pro Runde (mm, log)")ax.set_title("Konvergenz von Gauß-Seidel: fast eine Gerade")ax.grid(True, which="both", ls=":", alpha=0.5)plt.tight_layout()plt.show()
Abbildung 4.6: Der Gauß-Seidel-Fehler (größte Änderung pro Runde) über der Rundenzahl, mit logarithmisch geteilter senkrechter Achse. Weil jede Runde denselben Bruchteil des Fehlers wegnimmt, liegen die Punkte fast auf einer Geraden — die Signatur eines Verfahrens mit gleichmäßiger Konvergenz. Die gestrichelte Linie markiert die Abbruchschwelle 0,001; sie wird nach 22 Runden unterschritten.
Was man hier sieht: Auf der logarithmischen Achse liegen die Punkte nahezu auf einer Geraden, die stetig fällt. Das ist die Handschrift eines Verfahrens, das den Fehler Runde um Runde um denselben Bruchteil verkleinert. Wo die Punkte die rote Schwelle \(0{,}001\) unterschreiten — nach 22 Runden —, gilt die Lösung als gut genug. Für dieses winzige System sind 22 Runden viel (Thomas war in einem Durchlauf fertig); der Lohn kommt erst bei großen 2D-Netzen.
4.6 Was kostet das alles?
Jetzt können wir die Frage aus dem Vermutungs-Kasten mit Zahlen beantworten. „Kosten” meint hier die Anzahl der einfachen Rechenschritte (eine Multiplikation und eine Subtraktion), die ein Löser braucht — denn die Rechenzeit hängt fast nur daran.
Für die volle Gauß-Elimination muss man beim Leerräumen für jede der \(n\) Spalten alle darunter liegenden Zeilen anfassen, und in jeder Zeile alle Einträge. Das sind grob \(n \times n \times n\), also \(n^3\) Schritte — genauer etwa ein Drittel davon, \(n^3/3\). Für den Thomas-Algorithmus dagegen ist in jeder Zeile nur das schmale Band zu bearbeiten, immer gleich viel Arbeit pro Zeile; das ergibt grob \(8\,n\) Schritte — proportional zu \(n\), nicht zu \(n^3\). Setzen wir Zahlen ein:
Tabelle 4.2: Rechenschritte für volle Gauß-Elimination gegenüber dem Thomas-Löser. Bei jeder Verzehnfachung der Unbekannten wächst die volle Elimination um das Tausendfache (zehn hoch drei), der Thomas-Löser nur um das Zehnfache.
Unbekannte \(n\)
volle Elimination (\(\approx n^3/3\))
Thomas (\(\approx 8n\))
10
≈ 330
80
100
≈ 330 000
800
1 000
≈ 330 Millionen
8 000
1 000 000 (Crashtest)
≈ 3 · 10¹⁷
8 Millionen
Tabelle 4.2 macht den Vermutungs-Kasten greifbar. Übersetzen wir die beiden Extremzeilen in Zeit — für einen Rechner, der eine Milliarde einfache Schritte pro Sekunde schafft (eine vorsichtige Annahme für einen heutigen Computer):
Bei 1 000 Unbekannten: volle Elimination rund \(0{,}3\) Sekunden, Thomas ein Wimpernschlag von acht Millionstel Sekunden. Beides schnell — hier merkt man den Unterschied kaum.
Bei einer Million Unbekannten (Größenordnung eines echten Crashtests): volle Elimination rund \(3 \cdot 10^{17}\) Schritte — das wären über zehn Jahre Rechnen. Der Thomas-Löser bräuchte \(8\) Millionen Schritte, also acht Tausendstelsekunden.
Zehn Jahre gegen einen Wimpernschlag — das ist der ganze Unterschied zwischen „unmöglich” und „geht sofort”, und er kommt allein daher, dass der eine Löser die vielen Nullen ausnutzt und der andere nicht. Deshalb ist die Bandstruktur aus Kapitel 3 kein hübsches Detail, sondern die Bedingung dafür, dass große Simulationen überhaupt existieren.
HinweisDas Kleingedruckte: Warum wir nie umsortieren müssen
In Lehrbüchern zur Gauß-Elimination liest man oft von Pivotisierung — dem Vertauschen von Zeilen, damit auf der Diagonale keine Null (oder keine zu kleine Zahl) steht, durch die man dann nicht teilen könnte. Wir haben das im ganzen Kapitel nie gebraucht, und das ist kein Glück, sondern hat einen Grund: Die FEM-Matrizen dieses Buches sind diagonaldominant — auf der Diagonale steht immer die betragsgrößte Zahl der Zeile (bei der Federkette die \(200\) gegen zwei \(-100\) daneben). Das kommt direkt aus der Physik: Ein Knoten reagiert auf seine eigene Bewegung stärker als auf die jedes einzelnen Nachbarn. Solche Systeme kann man stur von oben nach unten eliminieren, ohne je umzusortieren — und Gauß-Seidel ist bei ihnen sogar garantiert gutmütig, es pendelt sich immer ein. Wann genau ein Verfahren konvergiert, ist eine eigene Wissenschaft, die wir hier nicht betreten; für unsere gutartigen Matrizen dürfen wir uns auf sie verlassen.
4.7 Die interaktive Einheit: der gläserne Löser
Jetzt bist du dran. Das folgende Programm ist ein gläserner Löser — er rechnet nicht nur, sondern zeigt bei jedem Schritt sein Inneres. Der erste Teil führt die Gauß-Elimination vor und druckt die Matrix nach jedem geräumten Spalten-Schritt (die Nullen als Punkte). Du kannst die Systemgröße zwischen 3 und 6 einstellen — das Programm baut dann eine Federketten-Matrix dieser Größe und räumt sie vor deinen Augen leer.
Der Code ist genau die Papier-Elimination von oben, in Python gegossen: dieselbe Doppelschleife, dieselben Namen (pivot, faktor, leerräumen).
WichtigVorhersage-Punkt
Bevor du ausführst: Stell GROESSE = 4 ein. Wie viele Eliminationsschritte (geräumte Spalten) wird es geben, bis die Dreiecksform steht? Und im zweiten Teil, dem Gauß-Seidel: Wie viele Runden braucht es wohl, bis sich kein Wert mehr um mehr als \(0{,}001\) ändert — mehr oder weniger als die 22 aus dem kleinen Beispiel? Leg dich fest, dann probiere es.
Abbildung 4.7: Vorgerenderte Fassung des gläsernen Lösers für vier freie Knoten: die gelöste Federkette als gleichmäßig ansteigende Treppe. Im Browser ersetzt dein eigenes Ergebnis dieses Bild, sobald du die Zelle darüber ausführst — stell die Systemgröße um und sieh zu, wie sich die Treppe verlängert.
Was man sieht: Für vier freie Knoten steht die Treppe \(0{,}5\), \(1{,}0\), \(1{,}5\), \(2{,}0\ \mathrm{mm}\) — vier gleiche Federn, jede um denselben halben Millimeter gedehnt. In der Textausgabe darüber siehst du, wie die Matrix Spalte für Spalte leergeräumt wird, bis die Dreiecksform steht. Probiere andere Größen und beobachte, wie viele Schritte das Leerräumen jeweils braucht — bei \(n\) freien Knoten sind es \(n-1\) geräumte Spalten.
4.8 Das Kapitel-Programm
Das vollständige, eigenständig lauffähige Skript zu diesem Kapitel liegt in programme/kap04/kap04_federkette_maschinell.py. Es ist die im Buch lange versprochene Verheiratung der beiden Bausteine: Der Matrix-Drucker aus Kapitel 3 (baue_systemmatrix) baut die Systemmatrix, der Löser aus diesem Kapitel (thomas) knackt sie. Zur Kontrolle löst es dasselbe System noch ein zweites Mal mit der vollen gauss_elimination — beide müssen aufs Tausendstel übereinstimmen.
Als Aufgabe nimmt es eine Kette aus acht Federn — mehr, als man bequem von Hand rechnet. Es assembliert die 9×9-Systemmatrix, druckt ihre Bandstruktur, baut die Randbedingung ein (Knoten 1 fest), liest die drei Diagonalen ab und löst. Heraus kommt die erwartete Treppe in \(0{,}5\)-mm-Schritten bis zu \(4{,}0\ \mathrm{mm}\) am freien Ende. So schließt sich der Bogen von Kapitel 3 (assemblieren) zu Kapitel 4 (lösen): Die FEM-Maschinerie steht jetzt vollständig — für die Federkette. Ab Kapitel 7 trägt dieselbe Maschinerie Wärme.
TippMerkkasten
Ein Gleichungssystem lösen heißt: Zahlen finden, die alle Gleichungen zugleich erfüllen — im 2×2-Fall der Schnittpunkt zweier Geraden.
Gauß-Elimination: erst unter der Diagonale leerräumen (Dreiecksform), dann von unten nach oben ernten (Rückwärtseinsetzen). Funktioniert immer, kostet aber rund \(n^3/3\) Schritte.
FEM-Matrizen sind dünn besetzt (Bandmatrizen): fast nur Nullen, weil jeder Knoten nur mit seinen Nachbarn zusammenhängt.
Der Thomas-Algorithmus ist die Gauß-Elimination, die nur das Band anfasst — eine Vorwärts- und eine Rückwärtswelle, Kosten nur \(\approx 8n\).
Gauß-Seidel entspannt: Jeder Knoten rückt sich wiederholt in den Mittelwert seiner Nachbarn; die Lösung pendelt sich ein. Abbruch, wenn sich nichts mehr merklich ändert.
Die vielen Nullen auszunutzen ist kein Luxus: Bei einer Million Unbekannten entscheidet es über zehn Jahre gegen einen Wimpernschlag.
Roter Faden
Zurück: Das 4×4-System der Federkette aus Kapitel 3 haben wir dort von Hand gelöst und das reduzierte 3×3-System als Probe behalten. Dieses Kapitel hat es maschinell gelöst — dreimal, mit Gauß, Thomas und Gauß-Seidel, und jedes Mal kamen die vertrauten \(0{,}5/1{,}0/1{,}5\ \mathrm{mm}\) heraus. Die Bandstruktur, die Kapitel 3 als hübsches Muster einführte, wurde hier zur Rechenersparnis.
Vor: Der Thomas-Löser wird ab Kapitel 7 den 1D-Träger lösen — dann tragen die Gleichungen Temperaturen statt Verschiebungen, aber die Matrix ist dieselbe tridiagonale Bandmatrix. Das Einpendeln von Gauß-Seidel wird in Kapitel 9 zur physikalischen Anschauung, wie sich der Träger zeitlich auf sein Temperaturprofil einschwingt, und in Kapitel 10 bändigt Gauß-Seidel die größeren 2D-Netze. Und die Kostenrechnung dieses Kapitels ist der Grund, warum das interaktive Finale in Kapitel 15 seine Netze deckelt: Reines Python im Browser verträgt nur so viele Unbekannte, wie ein Bandlöser in Sekunden schafft.
Übungen
Ü 4.1 (Verstehen). Löse das folgende 3×3-System von Hand mit Gauß-Elimination — erst leerräumen, dann ernten:
Schreibe jeden Schritt hin und prüfe dein Ergebnis anschließend mit dem gläsernen Löser, indem du das System dort einträgst.
HinweisMusterlösung zu Ü 4.1
Leerräumen, Spalte 1: Zeile 2 hat vorne eine \(1\), Zeile 1 eine \(2\). Ziehe die Hälfte von Zeile 1 von Zeile 2 ab (Faktor \(\tfrac{1}{2}\)): Zeile 2 wird \(0\ x + (3 - \tfrac{1}{2})\,y + z = 7 - 2\), also \(2{,}5\,y + z = 5\). Zeile 3 hat vorne schon eine \(0\).
Leerräumen, Spalte 2: Unter der \(2{,}5\) in Zeile 2 steht in Zeile 3 eine \(1\). Ziehe \(\tfrac{1}{2{,}5} = 0{,}4\) von Zeile 2 von Zeile 3 ab: Zeile 3 wird $(1 - 0{,}4)? $ … konkret \(y\)-Glied \(1 - 0{,}4\cdot 2{,}5 = 0\), \(z\)-Glied \(2 - 0{,}4\cdot 1 = 1{,}6\), rechte Seite \(5 - 0{,}4\cdot 5 = 3\). Also \(1{,}6\,z = 3\).
Ernten:\(z = 3/1{,}6 = 1{,}875\); dann \(2{,}5\,y = 5 - z = 3{,}125\), also \(y = 1{,}25\); dann \(2x = 4 - y = 2{,}75\), also \(x = 1{,}375\). Probe durch Einsetzen in die erste Zeile: \(2\cdot 1{,}375 + 1{,}25 = 4\). ✓
Ü 4.2 (Verändern). Gauß-Seidel braucht am reduzierten Federketten-System vom Nullstart aus 22 Runden. Füttere es stattdessen mit einem besseren Startwert — einem linear ansteigenden Profil, das die richtige Form schon erahnt (etwa \(0{,}4\), \(0{,}8\), \(1{,}2\)). Sage zuerst vorher: Spart das viele Runden oder wenige? Rechne dann nach; das Skript loesungen/kap04_ue2.py vergleicht beide Startwerte.
HinweisMusterlösung zu Ü 4.2
Der lineare Startwert \(0{,}4/0{,}8/1{,}2\) braucht nur 16 Runden statt 22 — sechs gespart. Der Grund: Die wahre Lösung ist eine gleichmäßige Treppe, also linear ansteigend. Der Nullstart hat von dieser Form gar nichts und muss sie erst aufbauen; das lineare Profil hat die Form schon und muss nur noch die Höhe nachjustieren. Beide landen bei denselben \(0{,}5/1{,}0/1{,}5\ \mathrm{mm}\) — ein guter Startwert ändert das Tempo, nie das Ziel. (Das ist der Grund, warum man in der Praxis eine grobe Näherung als Startwert verwendet, wo man eine hat — etwa die Lösung des vorigen Zeitschritts.)
Ü 4.3 (Übertragen). Verheirate den Matrix-Drucker aus Kapitel 3 mit dem Thomas-Löser aus diesem Kapitel und löse eine Kette aus acht gleichen Federn maschinell: assemblieren, Randbedingung einbauen, die drei Diagonalen ablesen, lösen. Welche Verschiebung hat das freie Ende? Das Skript loesungen/kap04_ue3.py baut es Schritt für Schritt (und das Kapitel-Programm programme/kap04/ zeigt die ausführliche Fassung).
HinweisMusterlösung zu Ü 4.3
Acht gleiche Federn, jede unter der durchgereichten Kraft von \(50\ \mathrm{N}\) um \(0{,}5\ \mathrm{mm}\) gedehnt: Die Knoten stehen bei \(0{,}5\), \(1{,}0\), \(1{,}5\), …, \(4{,}0\ \mathrm{mm}\). Das freie Ende (Knoten 9) kommt genau \(8 \cdot 0{,}5 = 4{,}0\ \mathrm{mm}\) weit. Der Thomas-Löser liefert das, ohne je ein Feld außerhalb des schmalen Bandes anzufassen — und die Gauß-Elimination als Kontrolle im Kapitel-Programm bestätigt es aufs Tausendstel.
Das Kleingedruckte
Ein paar ehrliche Feinheiten abseits des Hauptpfads. Erstens: Wir haben Gauß-Seidel nur an einem gutmütigen System gezeigt und behauptet, es pendle sich „immer” ein. Das gilt für die diagonaldominanten Matrizen dieses Buches; für beliebige Systeme kann Gauß-Seidel auch davonlaufen statt zu konvergieren. Wann genau es funktioniert, ist eine eigene Theorie, die wir bewusst nicht betreten — unsere FEM-Matrizen sind stets auf der sicheren Seite. Zweitens: Die Kostenformeln \(n^3/3\) und \(8n\) sind gerundete Faustzahlen; die genaue Schrittzahl hängt von Feinheiten der Umsetzung ab. Für die Größenordnung — und darum geht es hier — sind sie belastbar. Drittens: Es gibt weit mehr Löser als diese drei (Cholesky, konjugierte Gradienten, Mehrgitter …), und schnellere für riesige Netze. Wir nehmen die drei, weil man sie ganz durchschaut und weil sie den Rest des Buches tragen. Wer später ein großes Profi-Programm benutzt, wird andere Namen lesen — aber das Grundprinzip, die vielen Nullen auszunutzen, steckt in ihnen allen.
Thomas, Llewellyn H. 1949. Elliptic Problems in Linear Difference Equations over a Network. Watson Scientific Computing Laboratory, Columbia University.
