Ein Löser für 2-2-Bandmatrizen

Worum geht's?

Die Berechnung von Splinekurven erfordert die Lösung von Gleichungssystemen, wobei die Anzahl Gleichungen nahe der Anzahl Stützstellen liegt. Für kubische Splines sind die Matrizen der entsprechenden Gleichungssysteme nur entlang der Diagonalen sowie beiderseitss je einer Nebendiagonalen besetzt. Für diese Gleichungssysteme gibt es bekannte einfache Lösungen.

Für Splinefunktionen höherer Ordnung oder andere Randwertprobleme kann es nützlich sein, Gleichungssysteme mit je zwei Nebendiagonalen lösen zu können. Eine solche Lösung wird im Folgenden hergeleitet und gezeigt.

LR-Zerlegung

Die Lösung des Gleichungssystems beruht auf der LU-Zerlegung und der Tatsache - der Beweis sei hier unterlassen - dass die LU-Zerlegung eine Bandmatrix wieder eine Bandmatrix derselben Breite ist.

Gehen wir von der Ausgangsmatrix A mit 2-2-Bandstruktur aus und betrachten die beiden ersten Schritte der LR-Zerlegung:

1. Schritt

Der erste Schritt macht aus der Matrix A = A⁽¹⁾ die Matrix A⁽²⁾

Welche Elemente werden wie verändert?

• Die 1. Zeile bleibt unverändert
• Die 1. Spalte ab Zeile 2 entsteht durch Division der Elemente a₂₁...a_n1 durch a₁₁, wobei wegen der Bandstruktur nur a₂₁ und a₃₁ von Belang sind, die anderen sind Null
• Die Untermatrix a⁽²⁾₂₂ bis a⁽²⁾₃₃ entsteht durch das Subtraktion des dyadischen Vektorprodukts der Vektoren [a⁽²⁾₂₁, a⁽²⁾₃₁]^T mit [a⁽¹⁾₁₂, a⁽¹⁾₁₃] von [[a⁽¹⁾₂₂, a⁽¹⁾₂₃], [a⁽¹⁾₃₂, a⁽¹⁾₃₃]]
  a⁽²⁾₂₂ = a⁽¹⁾₂₂ – (a⁽¹⁾₁₂a⁽¹⁾₂₁)/a⁽¹⁾₁₁ usw.
• Spalte 4, in der Grafik orange eingefärbt, entsteht ebenfalls aus Spalte 4 durch Subtraktion eines Produkts aus dem Vektor der Spalte 1 mit Zeile 1. Diese enthält jedoch Null, so bleiben Spalte 4 und alle folgenden unverändert.
• Analoges gilt für die Zeile 4, gelb eingefärbt, und weiter bis Zeile n

Alles in Allem werden also nur Elemente in der Untermatrix A⁽¹⁾_2..3,1..3 verändert.

2. Schritt

Im 2. Schritt wird die Matrix A⁽²⁾ in die Matrix A⁽³⁾ überführt, und wieder werden nur die Elemente in einem 2x3-Block verändert.

Welche Elemente werden wir verändert?

• Zeile 2 bleibt unverändert
• Spalte 2, Elemente a⁽³⁾₃₂ und a⁽³⁾₄₂ entstehen aus a⁽²⁾₃₂ und a⁽²⁾₄₂ durch Division mit a⁽²⁾₂₂ aus dem 1. Schritt. Diese beiden Manipulationen müssen in einer Formel zusammengefasst werden: a⁽³⁾₃₂ = (a⁽¹⁾₃₂ – (a⁽²⁾₁₂a⁽²⁾₃₁)/a⁽²⁾₁₁)/a⁽²⁾₂₂) usw.
• Die Untermatrix a⁽³⁾₃₃ bis a⁽³⁾₄₄ entsteht durch Subtraktion des dyadischen Vektorprodukts der Vektoren [a⁽³⁾₃₂, a⁽³⁾₄₂]^T mit [a⁽²⁾₂₃, a⁽²⁾₂₄] von [[a⁽²⁾₃₃, a⁽²⁾₃₄], [a⁽²⁾₄₃, a⁽²⁾₄₄]]
  Davon ist speziell das Diagonalelement a⁽³⁾₃₃ zu erwähnen, welches aus a⁽²⁾₃₃ hervorgeht. Element a⁽²⁾₃₃ kann nicht zwischengespeichert werden, und deshalb wird der Berechnungsschritt von a⁽¹⁾₃₃ nach a⁽²⁾₃₃ auch in die Formel aufgenommen. Für die Elemente der 4. Zeile und 4. Spalte können die Eingangswerte aus der Originalmatrix verwendet werden

Für die weiteren Berechnungsschritte k liegen die Eingangswerte der Berechnungen entweder in A^(k-1), d.h. im Bereich des Rechenschemas, oder in der Originalmatrix A vor.

Beispiel

Gleichungssystem Ax - b = 0

Matrix A  und Bedingungsvektor b

Das implementierte Rechenschema kann 20 Gleichungen verarbeiten, aber auch weniger. Es werden keine Pivotelemente gesucht. Im Eingangsbereich von A sind die Spalten a₂ und a₁ die unteren Nebendiagonalen, Spalten c₁ und c₂ die oberen.

LR-Zerlegung und Vektor y

Die LR-Zerlegung von A in einem Format analog zur speichersparenden Zusammenlegung in der vollen Matrizenform: die L-Matrix, mit M bezeichnet, hat eine Diagonale mit lauter Einsen, die nicht gespeichert werden, sondern im Rechenschema eingesetzt sind.

Lösungsvektor x

Die Lösung des Gleichungssystems.