Numerische Methoden zur Lösung linearer, inhomogener Gleichungssysteme

Das grundsätzliche Problem

Gegeben sind die $n$ Gleichungen

$$\begin{array}{c} a_{11}x_{1}+a_{12}x_{2}+ \cdots +a_{1n}x_{n... ... a_{n1}x_{1}+a_{n2}x_{2}+ \cdots +a_{nn}x_{n}=b_{n} \end{array}$$ (2.1)

wobei die $a_{ij}$ und die $b_{i}$ reellwertige Größen sind. Außerdem soll gelten:

$$\mid b_{1} \mid + \mid b_{2} \mid + \cdots + \mid b_{n} \mid \ne 0$$

Unter diesen Voraussetzungen stellt (2.1) ein reellwertiges, lineares, inhomogenes Gleichungssystem $n$-ter Ordnung dar. Die Größen $x_{1},\ldots,x_{n}$, welche die gegebenen Gleichungen simultan erfüllen, nennt man die Lösungen des Systems.

(2.1) wird gewöhnlich in Form der Matrixgleichung

$$A \cdot {\bf x} = {\bf b}$$ (2.2)

angeschrieben.
Die Lösung derartiger Systeme stellt ein zentrales Problem der numerischen Mathematik dar, weil eine ganze Reihe von numerischen Verfahren, wie z.B. die Interpolationsverfahren, die Least-Squares-Verfahren, die Differenzenverfahren usw. auf die Lösung inhomogener linearer Gleichungssysteme zurückgeführt werden können.

Vom theoretischen Standpunkt aus bereitet die Lösung von (2.1) keinerlei Schwierigkeiten, solange die Determinante der Koeffizientenmatrix $A$ nicht verschwindet, solange also das Problem nicht singulär ist:

$$\mbox{det} (A) \ne 0$$

In diesem Fall kann das Problem mit Hilfe der Cramer'schen Regel angegangen werden. In der Praxis zeigt sich jedoch, daß die Anwendung dieser Regel für $n \geq 4$ recht kompliziert und aus verschiedenen Gründen für die Anwendung am Computer ungeeignet ist.

Wie auch in anderen Bereichen der numerischen Mathematik gibt es auch hier zwei Gruppen von Methoden:

• Direkte Methoden:
  

  Diese enthalten keine Verfahrensfehler und führen deshalb - abgesehen von Rundungsfehler-Einflüssen - stets zur exakten Lösung. Allerdings können Rundungsfehler gerade bei direkten Verfahren, die häufig sehr rechenintensive Algorithmen aufweisen, massiv in Erscheinung treten.

  Als Beispiel dazu wird im folgenden das Eliminationsverfahren von Gauss in der Formulierung von Doolittle und Crout (LU decomposition) behandelt.

• Iterative Methoden
  

  zeichnen sich oft durch einen besonders einfachen Rundungsfehler-stabilen und Speicherplatz-sparenden Algorithmus aus, sie haben jedoch den Nachteil, daß sie nicht in jedem Fall konvergieren, also nicht auf alle an sich lösbaren linearen Gleichungssysteme angewendet werden können.

  Außerdem erhält man auch im Falle der Konvergenz nie die exakte Problemlösung, sondern nur eine mit einem Abbruchfehler behaftete Näherungslösung.

  Als Beispiel für ein Iterationsverfahren zur Lösung linearer Gleichungssysteme steht in diesem Kapitel das Gauss-Seidel-Verfahren.

Ziel der direkten Methoden: Überführung der Koeffizientenmatrix in eine obere Dreiecksmatrix.

Bei allen direkten Methoden geht es letzten Endes darum, das gegebene System

$$A \cdot {\bf x} = {\bf b}$$

in ein äquivalentes (d.h. denselben Lösungsvektor aufweisendes) System

$$U \cdot {\bf x} = {\bf y}$$ (2.3)

überzuführen, wobei $U$ eine sogenannte obere Dreiecksmatrix mit der Eigenschaft

$$U=\left[ u_{ij} \right] \qquad \mbox{mit} \qquad u_{ij}=0 \qquad \mbox{f\uml ur} \qquad i>j$$

darstellt. Der Grund dafür ist einfach: Systeme der Art (2.3) sind durch Rück-Substitution (s. Glg.en (2.12)und (2.13)) leicht aufzulösen.

Das Eliminationsverfahren von Gauss in der Formulierung von Doolittle und Crout (LU-Aufspaltung).

Die Gauss'sche Elimination besteht im wesentlichen aus zwei Schritten, nämlich aus der Überführung von $A{\bf x}={\bf b}$ in das äquivalente System
$U{\bf x}={\bf y}$ und der Lösung dieses Systems mittels einer Rück-Substitution.

Das Verfahren beruht auf einem Satz aus der linearen Algebra:
Ein lineares Gleichungssystem bleibt unverändert, wenn zu einer seiner Gleichungen eine Linearkombination der restlichen Zeilen addiert wird.

Doolittle und Crout haben nun gezeigt, daß der Gauss'sche Algorithmus wie folgt formuliert werden kann: Eine reelle Matrix $A$ kann stets als Produkt zweier reeller Matrizen $L$ und $U$ dargestellt werden, also

$$A=L \cdot U \qquad$$ (2.4)

wobei $U$ eine obere Dreieckmatrix ist:

$$U = \left( \begin{array}{cccc} u_{11} & u_{12} & ....... ... ... . & & & . \\ 0 & 0 & ....... & u_{nn} \end{array} \right)$$

Die Matrix $L$ stellt eine untere Dreiecksmatrix der Form

$$L=\left( \begin{array}{ccccc} 1 & 0 & 0 & ....... & 0 ... ... m_{n1} & m_{n2} & m_{n3} & ....... & 1 \end{array} \right)$$

dar.

Man nennt die Zerlegung (2.4) einer Matrix $A$ die $LU$- Zerlegung ($LU$ decomposition) von $A$.

Ohne Ableitung: ein einfacher Formelapparat für die LU-Zerlegung.

Die Umformung erfolgt spaltenweise, also $j=1,2,\ldots,n$ :

$$u_{ij} = a_{ij} - \sum_{k=1}^{i-1} m_{ik} u_{kj} \qquad i=1,\ldots,j-1$$ (2.5)

$$\gamma_{ij}=a_{ij}-\sum_{k=1}^{j-1} m_{ik} u_{kj} \qquad i=j,\ldots,n$$ (2.6)

$$u_{jj}=\gamma_{jj}$$ (2.7)

$$m_{ij}=\frac{\gamma_{ij}}{\gamma_{jj}} \qquad i=j+1,\ldots,n$$ (2.8)

Da die Hauptdiagonalelemente von $L$ stets Einser sind, braucht man diese Werte (die $m_{jj}$) nicht extra abzuspeichern, sondern man kann alle relevanten Werte in einer $n$ x $n$-Matrix (der 'LU-Matrix') unterbringen:

$$\mbox{LU-Matrix}: \left( \begin{array}{ccccc} u_{11} & u_{12... ...m_{n1} & m_{n2} & m_{n3} & ....... & u_{nn} \end{array}\right)$$ (2.9)

Die Aufspaltung der Matrix $A$ in die Matrizen $L$ und $U$ kann nun zur Berechnung des Lösungsvektors herangezogen werden:

$$A \cdot {\bf x} \equiv L \cdot U \cdot {\bf x} = L \cdot (U \cdot {\bf x}) = {\bf b}$$

Mit $U\cdot {\bf x} \equiv {\bf y}$ erhält man das System $L \cdot {\bf y} = {\bf b}$, das wegen der besonderen Form der Matrix $L$ (untere Dreiecksmatrix) mittels einer Vorwärts-Substitution gelöst werden kann:

$$y_{1}=b_{1}$$ (2.10)

$$y_{i}=b_{i} - \sum_{j=1}^{i-1} m_{ij}y_{j} \qquad i=2,3,\ldots,n$$ (2.11)

Da $U$ eine obere Dreiecksmatrix ist, ist bei Kenntnis des Hilfsvektors ${\bf y}$ das System $U \cdot {\bf x}={\bf y}$ durch Rückwärts-Substitution ebenfalls leicht zu lösen:

$$x_{n}=\frac{y_{n}}{u_{nn}}$$ (2.12)

$$x_{i}=\frac{1}{u_{ii}} \left[ y_{i}-\sum_{j=i+1}^{n} u_{ij}x_{j} \right] \qquad i=n-1,n-2,\ldots,1$$ (2.13)

Bei der konkreten Anwendung des Doolittle-Crout-Verfahrens geht man am besten so vor, daß man zwei Programme verwendet, von denen eines die LU-Zerlegung der Koeffizientenmatrix durchführt, während das zweite den Lösungsvektor unter Verwendung der Formeln (2.10-2.13) berechnet.

Demonstration einer speicherplatz-sparenden LU-Decomposition

an Hand einer allgemeinen 3x3-Matrix.

Gegeben:

$$\left(\begin{array}{ccc}a_{11}&a_{12}&a_{13} a_{21}&a_{22}&a_{23}\\ a_{31}&a_{32}&a_{33}\end{array}\right)$$

j = 1: erste Spalte

(2.6)	i=1	$\gamma_{11}=a_{11}$
	i=2	$\gamma_{21}=a_{21}$
	i=3	$\gamma_{31}=a_{31}$
.
.
(2.7)		$u_{11}=\gamma_{11}$
.
.
(2.8)	i=2	$m_{21}=\gamma_{21}/ \gamma_{11}$
	i=3	$m_{31}=\gamma_{31}/\gamma_{11}$

Speichert man nun die eben berechneten Koeffizienten $u_{11}$, $m_{21}$ und $m_{31}$ anstelle der entsprechenden und nie wieder gebrauchten Koeffizienten der Matrix A ab, so ergibt sich die Matrix

$$\left(\begin{array}{ccc} u_{11} & a_{12} & a_{13} \\ m_{21} & a_{22} & a_{23} m_{31} & a_{32} & a_{33} \end{array}\right)$$

j=2: zweite Spalte

(2.5)	i=1	$u_{12}=a_{12}$
.
.
(2.6)	i=2	$\gamma_{22}=a_{22}- m_{21} u_{12}$
	i=3	$\gamma_{32}=a_{32}-m_{31} u_{12}$
.
.
(2.7)		$u_{22}=\gamma_{22}$
.
.
(2.8)	i=3	$m_{32}= \gamma_{32}/\gamma_{22}$

$$\left(\begin{array}{ccc} u_{11} & u_{12} & a_{13} \\ m_{21} & u_{22} & a_{23} m_{31} & m_{32} & a_{33} \end{array}\right)$$

j=3: dritte Spalte

(2.5)	i=1	$u_{13}=a_{13}$
	i=2	$u_{23}=a_{23}- m_{21} u_{13}$
.
.
(2.6)i=3		$\gamma_{33}=a_{33}-m_{31}u_{13}-m_{32}u_{23}$
.
.
(2.7)		$u_{33}=\gamma_{33}$

$$\left(\begin{array}{ccc} u_{11} & u_{12} & u_{13} \\ m_{21}... ...{33} \end{array}\right) = \mbox{LU-Matrix\quad [s. Glg.(2.9)]}$$

Entscheidend ist also, daß durch die sukzessive Überspeicherung der gegebenen Koeffizientenmatrix $A$ durch die Spalte für Spalte entstehende LU-Matrix zu keiner Zeit ein Informationsmanko entsteht!

Da andererseits auch die neuen $u_{ij}$- bzw. $m_{ij}$-Werte stets verschiedene Speicherplätze einnehmen, kann im Programm auf eine unterschiedliche Bezeichnung der Elemente $a_{ij}$, $u_{ij}$, $\gamma_{ij}$ und $m_{ij}$ überhaupt verzichtet werden, und alle diese Elemente werden im Struktogramm 3 (LUDCMP) auf einem Feld A(,) abgespeichert!
Man spart damit sehr viel Speicherplatz, allerdings auf Kosten der Lesbarkeit des Programms.

Rundungsfehler-Optimierung durch partielle Pivotisierung.

Angenommen, es existiert ein Programm, das die Auswertung eines inhomogenen Gleichungssystems nach der LU-Zerlegung durchführt. Die mittels dieses Programmes erhaltenen numerischen Ergebnisse enthalten keine Verfahrensfehler (direktes Verfahren!); jede Abweichung von den 'wahren' Ergebnissen muß daher auf Rundungsfehler zurückgehen.

Dazu ein konkretes Beispiel: Gegeben sei das folgende System von 3 Gleichungen:

$$x_{1} + 5923181 x_{2} + 1608 x_{3} = 5924790$$

$$5923181 x_{1} + 337116 x_{2} - 7 x_{3} = 6260290$$ (2.14)

$$6114 x_{1} + 2 x_{2} + 9101372 x_{3} = 9107488$$

Dieses System hat die exakte Lösung $x_{1}=x_{2}=x_{3}=1$.

Die numerische Auswertung (FORTAN, einfache Genauigkeit, d.h. 4 Byte pro reeller Zahl) liefert die folgende LU-Matrix und den folgenden Lösungsvektor:

	0.1000000E+01	0.5923181E+07	0.1608000E+04
LU-Matrix:	0.5923181E+07	-.3508407E+14	-.9524475E+10
	0.6114000E+04	0.1032216E-02	0.9101371E+07
		0.9398398E+00
x:		0.1000000E+01
		0.9995983E+00

Wie man sieht, ist der Rundungsfehler-Einfluß (besonders beim $x_{1}$) beträchtlich!

Wenn man nun die Reihenfolge der Gleichungen in (2.14) abändert (was theoretisch natürlich keinen Einfluß auf den Lösungsvektor hat), und zwar so, daß die ursprüngliche Reihenfolge $1/2/3$ übergeht in die Reihenfolge $2/3/1$, so ergibt dasselbe Programm:

	0.5923181E+07	0.3371160E+06	-.7000000E+01
LU-Matrix:	0.1032216E-02	-.3459764E+03	0.9101372E+07
	0.1688282E-06	-.1712019E+05	0.1558172E+12
		0.9999961E+00
x:		0.1000068E+01
		0.1000000E+01

mit bereits wesentlich kleineren Rundungsfehlern. Mit der Zeilenfolge $2/1/3$ erhält man den korrekten Lösungsvektor in Maschinengenauigkeit:

	0.5923181E+07	0.3371160E+06	-.7000000E+01
LU-Matrix:	0.1688282E-06	0.5923181E+07	0.1608000E+04
	0.1032216E-02	-.5841058E-04	0.9101372E+07
		0.1000000E+01
x:		0.1000000E+01
		0.1000000E+01

Das bedeutet also, daß die richtige Zeilenfolge eines linearen Gleichungssystems einen entscheidenden Einfluß auf die Größe der Rundungsfehler und somit auf die Qualität der numerischen Lösung hat!

Was bedeutet nun 'richtige' Reihenfolge? Vergleicht man die $m_{ij}$-Werte in den LU-Matrizen, so erkennt man sofort, daß die Rundungsfehler dann am geringsten sind, wenn die Beträge der $m_{ij}$-Werte möglichst klein sind.
Um die $m_{ij}$ so klein als möglich zu halten, müssen die in der Gleichung (2.8) auftretenden $\gamma_{jj}$-Werte dem Betrage nach möglichst groß sein. Dies kann einfach dadurch erreicht werden, daß man aus allen $\gamma_{ij}$ nach Glg. (2.6) das betragsgrößte Element bestimmt und jene Zeile, in welcher dieses Element vorkommt, mit der $j$-ten Zeile vertauscht. Erst danach werden die Gleichungen (2.7) und (2.8) ausgewertet.
Eine derartige Strategie nennt man eine LU-Zerlegung mit teilweiser Pivotisierung. In der Regel hat man - vor allem bei Systemen höherer Ordnung - ohne Pivotisierung keine Chance auf korrekte Ergebnisse. Es gibt aber auch (s. [2], S. 56) Gleichungssysteme, bei welchen die LU-Decomposition auch ohne Pivotisierung numerisch stabil ist, nämlich dann, wenn die Koeffizientenmatrix symmetrisch, tridiagonal, zyklisch tridiagonal, diagonal-dominant bzw. allgemein positiv definit ist.

Es ist klar, daß mit der Strategie der teilweisen Pivotisierung auch das Nullwerden eines Diagonalelementes $\gamma_{jj}$ und die damit auftretenden Probleme in Glg.(2.8) vermieden werden können. Stellt sich nämlich heraus, daß während der Rechnung alle $\gamma_{ij},\quad i=j,\ldots,n$ simultan verschwinden, ist die Koeffizientenmatrix singulär und das Gleichungssystem mit dem gegebenen Verfahren nicht lösbar!

Kondition eines Gleichungssystems.

Die mit Hilfe direkter Methoden ermittelte Lösung eines linearen Gleichungssystems ist meist nicht die exakte Lösung, weil

• im Verlauf der Rechnung Rundungsfehler auftreten, deren Akkumulation zur Verfälschung der Ergebnisse führen kann.
• Ungenauigkeiten in den Ausgangsgrößen bestehen können, die Ungenauigkeiten in den Lösungen hervorrufen.

Wenn kleine Änderungen in den Ausgangsdaten große Änderungen in der Lösung hervorrufen, spricht man von einem schlecht konditionierten System (s. Abb. 2.1).

Abbildung 2.1: Gut und schlecht konditioniertes 2x2-System (aus: [22], S. 228).

Die sog. Konditionszahl der Koeffizientenmatrix kann dabei als Maß für die Güte der erhaltenen Näherungslösung für den Lösungsvektor verwendet werden.

In der mathematischen Literatur wird eine Reihe von Konditionszahlen diskutiert. Im folgenden Programm LUDCMP ist die Berechnung der sog. Hadamarschen Konditionszahl der Matrix $A$ inkludiert, welche in der Form [2]

$$K_{{\rm H}}(A)=\frac{\vert det(A)\vert}{\alpha_{1}\alpha_{2... ...pha_{i} = \sqrt{a_{i1}^{2} + a_{i2}^{2} + \ldots +a_{in}^{2}}$$

gegeben ist. Für die Praxis gelten die Erfahrungswerte

• $K_{{\rm H}}(A) < 0.01$        schlechte Kondition,
• $K_{{\rm H}}(A) > 0.1$        gute Kondition,
• $0.01 \leq K_{{\rm H}}(A) \leq 0.1$        keine Aussage.

Das Programm LUDCMP

Quelle: [9], S.38f; [10], S. 46f (mit geringen Änderungen).

Das Programm LUDCMP (LU DeCoMPosition) führt eine LU-Aufspaltung einer reellen quadratischen Matrix durch.

INPUT Parameter:

A( , ):

Koeffizientenmatrix des linearen Systems.

N:

Ordnung des Systems = Anzahl der Zeilen und Spalten von A.

OUTPUT Parameter:

A( , ):

'LU Matrix' der ursprünglichen Matrix A (s. Anmerkung 1).

INDX( ):

Indexvektor, der die Zeilenvertauschungen auf Grund der partiellen Pivotisierung speichert.

D:

Determinante der Koeffizientenmatrix (s. Abschnitt 2.3.6).

KHAD:

Hadamardsche Konditionszahl (s. Abschnitt 2.3.2).

=15cm Struktogramm 3LUDCMP(A,N,INDX,D,KHAD)KHAD:=1.0I=1(1)NSUM:=0.0J=1(1)NSUM:=SUM + A(I,J)**2KHAD:=KHAD*SQRT(SUM)D:=1.0J=1(1)NI=1(1)J-1SUM:=A(I,J)K=1(1)I-1SUM:=SUM-A(I,K)*A(K,J)A(I,J):=SUMAAMAX:=0.0I=J(1)NSUM:=A(I,J)K=1(1)J-1SUM:=SUM-A(I,K)*A(K,J)A(I,J):=SUM     DUM:=$\vert$ SUM$\vert$DUM $\ge$ AAMAXIMAX:=I     AAMAX:=DUM......J $\ne$ IMAXK=1(1)NDUM:=A(IMAX,K)
A(IMAX,K):=A(J,K)
A(J,K):=DUMD:=-D......INDX(J):=IMAXA(J,J) = 0.0A(J,J):=TINY......

...DUM:=1.0/A(J,J)I=J+1(1)NA(I,J):=A(I,J)*DUMD:=D*A(J,J)KHAD:=$\vert$ D $\vert$/KHAD(return)

Anmerkungen zum Programm LUDCMP:

1. Um Speicherplatz zu sparen, werden die Koeffizienten der LU Matrix sukzessive auf den entsprechenden Speicherplätzen der gegebenen Matrix A abgespeichert. Dies hat allerdings den Nachteil, daß die gegebene Matrix zerstört wird und daher gegebenenfalls (z.B. wenn A für spätere Rechnungen wieder gebraucht wird) vor Aufruf von
   LUDCMP gesichert werden muß.
2. Der Output-Parameter D hat Bedeutung bei der Berechnung der Determinante von A.
3. Die Konstante TINY ist in [9] mit $10^{-20}$ definiert. Wie aus dem Struktogramm 3 hervorgeht, verhindert sie ein crash-down des Programmes, wenn ein Divisor in Glg.(2.8) exakt Null wird. Zu diesem Thema siehe Kap.2.5.
4. In [9] werden die Matrixkoeffizienten für die Pivotsuche zeilenweise skaliert. Da optimale Skalierungen von linearen Systemen ein überraschend schwieriges Problem darstellen (s. z.B. [7], S.140ff und [12], S.41ff), wurde diese Skalierung aus dem hier präsentierten Programm herausgenommen.

Das Unterprogramm LUBKSB

Quelle: [9], S. 39; [10], S. 47 (mit Änderungen).

Das Programm LUBKSB (LU BacKSuBstitution) berechnet den Lösungsvektor x des Systems $LU{\bf x}={\bf b}$ unter Verwendung der Gleichungen (2.10 - 2.13).

INPUT-Parameter:

A( , ):

LU-Matrix der Koeffizientenmatrix, wie sie z.B. vom Programm LUDCMP berechnet wird.

N:

Ordnung des Systems = Zeilen und Spalten von A.

INDX( ):

Indexvektor, der Informationen über die von LUDCMP vorgenommenen Zeilenvertauschungen enthält.

B( ):

inhomogener Vektor des Systems.

OUTPUT Parameter:

X( ):

Lösungsvektor.

Interner Vektor:

Y( ):

.

Struktogramm 4LUBKSB(A,N,INDX,B,X)I=1(1)NBB(I):=B(I)I=1(1)NLL:=INDX(I)
SUM:=BB(LL)
BB(LL):=BB(I)J=1(1)I-1SUM:=SUM-A(I,J)*Y(J)Y(I):=SUMI=N(-1)1SUM:=Y(I)J=I+1(1)NSUM:=SUM-A(I,J)*X(J)X(I):=SUM/A(I,I)(return)

Die Verwendungsmöglichkeiten für die Programme LUDCMP und LUBKSB.

Lösung eines inhomogenen Gleichungssystems $A{\bf x}={\bf b}$:

                ||---------------------------||
                ||  LUDCMP(A,N,INDX,D,KHAD)  ||
                ||---------------------------||
                ||   LUBKSB(A,N,INDX,B,X)    ||
                ||---------------------------||

Ein wesentlicher Vorteil des hier vorgestellten Verfahrens besteht darin, daß bei Auftreten verschiedener inhomogener Vektoren ${\bf b}_{1}$, ${\bf b}_{2}$, ..., aber gleichbleibender Koeffizientenmatrix $A$ die LU-Zerlegung nur einmal durchgeführt werden muß:

                 ||-------------------------||
                 || LUDCMP(A,N,INDX,D,KHAD) ||
                 ||-------------------------||
                              |
                 ||-------------------------||
                 || LUBKSB(A,N,INDX,B1,X1)  ||
                 || LUBKSB(A,N,INDX,B2,X2)  ||
                 ||           .             ||
                 ||           .             ||
                 ||-------------------------||

Invertierung einer Matrix:

Mit Hilfe der Programme LUDCMP und LUBKSB kann eine gegebene Matrix $A$ Spalte für Spalte invertiert werden:

$$A \cdot X = E \qquad E = \mbox{Einheitsmatrix} \qquad X \equiv A^{-1}$$

Diese Matrixgleichung kann man nun in eine Reihe von inhomogenen linearen Gleichungssystemen der Art

$$A \left( \begin{array}{c} x_{11} . . . x_{n1} \end... ... \left( \begin{array}{c} 0 0 . . 1 \end{array}\right)$$

zerlegen, d.h. man hat $n$ Systeme mit gleicher Koeffizientenmatrix $A$ und mit wechselnden inhomogenen Vektoren (= Spaltenvektoren der Einheitsmatrix) zu berechnen:

StruktogrammInvertierung einer MatrixA,N,INDX,D,KHADLUDCMPJ=1(1)NI=1(1)NVEKTOR(I):=0.0VEKTOR(J):=1.0A,N,INDX,VEKTOR,XLUBKSBI=1(1)NAINV(I,J):=X(I)

Berechnung der Determinante einer Matrix:

Für die Determinante der Matrix A gilt

$$det(A) = det(L\cdot U) = det(L) det(U) .$$

Weil beide Matrizen $L$ und $U$ Dreiecksmatrizen sind, sind deren Determinanten gleich den jeweiligen Produkten der Diagonalelemente. Weil aber alle Diagonalelemente der $L$-Matrix den Wert 1 haben, gilt

$$det(L) = 1 ,$$

und es ergibt sich

$$det(A) = det(U) = \sum_{i=1}^n u_{ii} .$$

Die Berechnung von $det(U)$ erfolgt während der LU-Dekomposition im Programm LUDCMP; dabei ist zu berücksichtigen, daß bei jeder Zeilenvertauschung sich das Vorzeichen der Determinante ändert.

Testbeispiele für die Programme LUDCMP und LUBKSB.

Gegeben ist die Matrix mit vier Zeilen und Spalten

$$A=\left( \begin{array}{cccc} 1.1161 & 0.1254 & 0.1397 & 0.... ...71 \\ 0.2368 & 0.2471 & 0.2568 & 1.2671 \end{array} \right)$$

Die Konditionszahl lautet $K_{{\rm H}}$ = 0.752, das System ist also gut konditioniert!

• Wie lautet der Lösungsvektor ${\bf x}$ des linearen inhomogenen Gleichungssystems
  
  $$A \cdot {\bf x} = (-1.8367, 1.1944, 3.2368, -0.7232)^{T} \quad ?$$
  
  

  Das numerische Ergebnis ist bis auf Maschinengenauigkeit exakt:
  

  $${\bf x}=(-.2000000E+01, 0.1000000E+01, 0.3000000E+01, -.1000000E+01)^{T}$$
  
  

• Wie lautet die inverse Matrix von $A$?

  Das numerische Ergebnis ist fast bis auf Maschinengenauigkeit exakt (bis auf maximal zwei Einheiten in der siebenten Stelle hinterm Komma) und lautet:

  Invertierte Matrix:
  
     .9379443E+00   -.6843720E-01   -.7960770E-01   -.8592076E-01
    -.8852433E-01    .9059826E+00   -.9919081E-01   -.1055899E+00
    -.1113511E+00   -.1169667E+00    .8784252E+00   -.1270733E+00
    -.1354557E+00   -.1401826E+00   -.1438075E+00    .8516058E+00
  
  Test-Matrix:
  
     .9999999E+00    .1126516E-07    .2518815E-08   -.1120211E-08
    -.4888310E-08    .1000000E+01   -.6621389E-08   -.3829147E-08
    -.2647183E-08    .5092004E-08    .9999999E+00    .1665558E-07
     .7356128E-08   -.1727934E-07   -.1442864E-07    .1000000E+01
  

  Die Test-Matrix stellt einfach das Multiplikations-Ergebnis: gegebene Matrix mal (numerisch) invertierte Matrix dar. Zu erwarten wäre im Idealfall die Einheitsmatrix.

• Welchen Wert hat die Determinante der Matrix $A$?

  Das numerische Ergebnis (auf Maschinengenauigkeit exakt) lautet
  

  $$\mbox{Determinante} = 0.1758306E+01$$
  
  

Verbesserung eines Lösungsvektors durch Nachiteration.

Es ist klar, daß die Genauigkeit der numerischen Lösung eines linearen Gleichungssystems durch die Maschinengenauigkeit begrenzt ist. Leider ist es in vielen Fällen so, daß diese Genauigkeit wegen auftretender Rundungsfehler-Kumulationen nicht erreicht wird.

In solchen Fällen gibt es die Möglichkeit einer iterativen Verbesserung der numerischen Lösung, deren (sehr einfache) Theorie nun kurz erläutert wird:

Die numerische Lösung des Systems

$$A \cdot {\bf x} = {\bf b}$$

sei der durch Rundungsfehler verfälschte Vektor $\bar{\bf x}= {\bf x} + \delta{\bf x}$. Setzt man nun $\bar{\bf x}$ in das Gleichungssystem ein, so erhält man den verfälschten inhomogenen Vektor ${\bf b}+\delta{\bf b}$:

$$A\cdot \bar{\bf x} = A({\bf x}+\delta{\bf x}) = {\bf b}+\delta{\bf b}$$

bzw.

$$\underbrace{(A{\bf x}-{\bf b})}_{=0} + A\cdot \delta{\bf x} = \delta{\bf b} .$$

Somit ergibt sich der Fehlervektor $\delta{\bf x}$ als Lösung des Systems

$$A \cdot {\bf\delta x} = {\bf\delta b} \quad ,$$

wobei $A$ die unveränderte Koeffizientenmatrix und ${\bf\delta b}$ den sogenannten Residuenvektor $A \bar{\bf x}-{\bf b}$ darstellt.

Das Struktogramm 5 zeigt die Lösung eines linearen Gleichungssystems mit anschließender iterativer Verbesserung des Lösungsvektors:

$${\bf x} = \bar {\bf x} - \delta {\bf x}\quad .$$

Dazu noch einige Anmerkungen:

1. Da beim ersten Aufruf von LUDCMP die Matrix $A$ zerstört wird, muß diese Größe rechtzeitig gesichert werden.
2. Wegen der großen Gefahr einer subtractive cancellation sollte die Berechnung des Residuenvektors - wenn irgend möglich - mit doppelter Genauigkeit erfolgen.
3. Manche Programm-Bibliotheken (z.B. NAG) enthalten Programme zur Lösung inhomogener Gleichungssysteme, die eine iterative Verbesserung bereits inkludiert haben.

Struktogramm 5Iterative Verbesserung des LösungsvektorsI=1(1)NJ=1(1)NASP(I,J):=A(I,J)A,N,INDX,D,KHADLUDCMPA,N,INDX,B,XLUBKSBI=1(1)NSUMDP:=-B(I)J=1(1)NSUMDP:=SUMDP + DBLE(ASP(I,J))*DBLE(X(J))RES(I):=SUMDPA,N,INDX,RES,DELXLUBKSBI=1(1)NX(I):=X(I)-DELX(I)(Bedingung, unter welcher die Iteration
verlassen wird)(Ausgabe und return)

Rundungsfehler-Probleme mit schlecht konditionierten bzw. singulären Systemen.

Die Methode der partiellen Pivotisierung bewirkt eine entscheidende Reduktion der beim Gauss'schen Verfahren auftretenden Rundungsfehler. Dennoch kann bei schlecht konditionierten Gleichungssystemen der numerisch ermittelte Lösungsvektor deutlich fehlerbehaftet sein, wie das folgende Beispiel zeigt:

0.10000E+01 0.50000E+00 0.33333E+00 0.25000E+00 0.20000E+00   0.228333E+01

0.50000E+00 0.33333E+00 0.25000E+00 0.20000E+00 0.16667E+00   0.145000E+01

0.33333E+00 0.25000E+00 0.20000E+00 0.16667E+00 0.14286E+00   0.109286E+01

0.25000E+00 0.20000E+00 0.16667E+00 0.14286E+00 0.12500E+00   0.884530E+00

0.20000E+00 0.16667E+00 0.14286E+00 0.12500E+00 0.11111E+00   0.745640E+00

Die Koeffizientenmatrix stellt eine Hilbert-Matrix 5. Ordnung mit den auf 5 signifikante Stellen gerundeten Komponenten

$$a_{ij}=\frac{1}{i+j-1}$$

dar, und der inhomogene Vektor ist so gewählt, daß die exakte Lösung

$$x_{1}=x_{2}=x_{3}=x_{4}=x_{5}=1$$

lautet. Die Konditionszahl lautet $K_{{\rm H}}$ = 0.55 $\cdot 10^{-10}$; es liegt also ein sehr schlecht konditioniertes System vor, und Schwierigkeiten sind zu erwarten!

Das exakte Ergebnis wird im folgenden dem numerischen Ergebnis $\bar{\bf x}$ gegenübergestellt (FORTRAN, single precision). Weiters ist zur Kontrolle der Vektor $A \cdot \bar{\bf x}$ ausgegeben:

${\bf x}$ (exakt)	$\bar{\bf x}$ (numerisch)	$A \cdot \bar{\bf x}$
1.	0.9999459E+00	0.2283330E+01
1.	0.1000955E+01	0.1450000E+01
1.	0.9960009E+00	0.1092860E+01
1.	0.1005933E+01	0.8845300E+00
1.	0.9971328E+00	0.7456400E+00

Das Dilemma ist deutlich zu sehen: Obwohl einige Komponenten des numerisch erhaltenen Lösungsvektors beträchtlich vom exakten Wert abweichen, ist das lineare Gleichungssystem mit Maschinengenauigkeit erfüllt! Auch eine Nachiteration bringt keine Verbesserung der Ergebnisse.

Dieses Beispiel enthüllt noch ein weiteres ernstes Problem bei der Verwendung numerischer Methoden zur Lösung linearer Systeme: Es gibt trotz zahlreicher theoretischer Untersuchungen (s. z.B. [12], S.109ff) keine einfache, also in der Praxis einsetzbare Abschätzung des im Lösungsvektor enthaltenen Rundungsfehlers.

Eine weitere unangenehme Rundungsfehler-Konsequenz besteht darin, daß es bei der numerischen Auswertung praktisch unmöglich ist, klar zwischen singulären und fast-singulären Koeffizientenmatrizen zu unterscheiden. Im Programm wird eine Matrix dann als singulär erkannt, wenn alle Elemente einer Pivot-Spalte exakt Null sind. Da aber diese Elemente ihrerseits Ergebnisse arithmetischer Operationen sind, ist wegen der dabei auftretenden Rundungsfehler auch im Falle der Singularität der Koeffizientenmatrix nicht mit exakten Nullen zu rechnen (das folgende aus [12], S.63):

Die Lage wird dadurch besonders kompliziert, daß die scharfe mathematische Unterscheidung von singulären und nicht-singulären Matrizen nur in der Idealwelt des Mathematikers existiert. Sobald wir mit Matrizen in gerundeter Arithmetik arbeiten, wird die Unterscheidung notwendigerweise ungerechtfertigt. So können gewisse nicht-singuläre Matrizen als Ergebnis von durch die Rundungsfehler eingeführten Störungen singulär werden. Noch wahrscheinlicher wird eine wirklich singuläre vorgegebene Matrix durch den Rundungsfehler in eine benachbarte nicht-singuläre Matrix abgeändert.

Aus diesem Grund verzichten viele Programme auf eine Singularitäts-Diagnose und überlassen diese dem Benutzer.
Im Programm LUDCMP wird zwar abgefragt, ob im Laufe der Rechnung irgendein Diagonalelement $\gamma_{jj}$ exakt(!) Null ist, aber nur, um in diesem (unwahrscheinlichen) Fall eine Division durch Null zu verhindern!

Direkte Lösungsverfahren für Systeme mit speziellen Koeffizientenmatrizen.

Nachdem direkte Verfahren zur Lösung linearer Gleichungssysteme relativ rechen- und speicherplatz-intensiv sind, ist es in der Praxis von sehr großer Bedeutung, spezielle Koeffizientenmatrizen auszunutzen, um vereinfachte Algorithmen zu formulieren.

Lösung von Gleichungssystemen mit tridiagonaler Koeffizientenmatrix.

Eine Reihe wichtiger numerischer Verfahren (z.B. die Spline-Interpolation) führen zu linearen Gleichungssystemen der Form

$$\left( \begin{array}{cccccc} b_{1} & c_{1} & 0 & 0 & \cdots ... ... r_{2} r_{3} . . . r_{n} \end{array} \right)$$ (2.15)

wobei die tridiagonale Koeffizientenmatrix eindeutig durch die drei Vektoren ${\bf b}$ (Hauptdiagonale) sowie ${\bf a}$ und ${\bf c}$ (untere bzw. obere Nebendiagonale) gegeben ist. Wendet man in diesem Fall die LU-Zerlegung (ohne Pivotisierung, s.u.) an, so erhält man die Struktur

$$\left( \begin{array}{ccccc} 1 & 0 & 0 & \cdots & 0 \\ m_... ...} r_{1} r_{2} . . . r_{n} \end{array} \right)$$

wobei sich bei Anwendung der Glg.en (2.5-2.13) die folgenden einfachen Formeln ergeben:

$$u_{1} = b_{1}$$

$$y_{1} = r_{1}$$

$$m_{j} = a_{j}/u_{j-1}$$

$$u_{j} = b_{j} - m_{j} \cdot c_{j-1} \qquad j=2,\ldots,n$$ (2.16)

$$y_{j} = r_{j} - m_{j} \cdot y_{j-1}$$

Die Komponenten des Lösungsvektors ergeben sich wieder mittels Rück-Substitution:

$$x_{n} = y_{n}/u_{n}$$

$$x_{j} = (y_{j} - c_{j} \cdot x_{j+1})/u_{j} \qquad j=n-1,\ldots,1$$ (2.17)

Diese Gleichungen sind die Basis für das folgende Programm.

Das Programm TRID.

Quelle: [9], S. 43 (mit Änderungen).

Das Programm TRID berechnet den Lösungsvektor eines linearen, inhomogenen Gleichungssystems mit tridiagonaler Koeffizientenmatrix.

INPUT-Parameter:

A( ),B( ),C( ):

Vektoren der tridiagonalen Matrix.

R( ):

inhomogener Vektor.

N:

Ordnung des Systems.

OUTPUT-Parameter:

X( ):

Lösungsvektor.

Interne Vektoren:

U( ), Y( ):

.

Struktogramm 6TRID(A,B,C,R,N,X)Y(1):=R(1)
U(1):=B(1)U(1) = 0.0(Abbruch mit Fehlermitteilung!)......J=2(1)NM:=A(J)/U(J-1)
U(J):=B(J)-M*C(J-1)U(J) = 0.0(Abbruch mit Fehlermitteilung!)......Y(J):=R(J)-M*Y(J-1)X(N):=Y(N)/U(N)J=N-1(-1)1X(J):=(Y(J)-C(J)*X(J+1))/U(J)(return)

Anmerkungen zu TRID:

1. Der Hauptvorteil von TRID liegt natürlich in einer beträchtlichen Speicherplatz-Ersparnis, da ja als Input-Informationen nur die vier eindimensionalen Felder A, B, C und R eingegeben werden müssen. Im Programm kommen dazu noch die ebenfalls eindimensionalen internen Felder U und Y sowie der Lösungsvektor X. Bei einem Problem n-ter Ordnung bedeutet das einen Speicherplatzbedarf von $7n$ gegenüber $n^{2}+n$ beim normalen LU-Verfahren.
2. Ein weiterer Vorteil gegenüber dem normalen LU-Verfahren ist die beträchtlich verminderte Rechenzeit.
3. Allerdings bedeutet der Verzicht auf die Pivotisierung, daß es auch bei nicht-singulären Problemen zu einem Abbruch des Verfahrens kommen kann! Prinzipiell würde der Einbau einer Pivotisierung ins Programm TRID möglich sein, es würde sich jedoch dabei die Bandbreite der Matrix erhöhen (und zwar im ungünstigsten Fall verdoppeln).

   Es kann jedoch gezeigt werden, daß bei allen sogenannten diagonal- dominanten tridiagonalen Matrizen, d.h. bei solchen, die die Eigenschaft $\mid b_{j} \mid > \mid a_{j} \mid + \mid c_{j} \mid$ für alle $j=1,\ldots,n$ haben, keine derartigen Probleme auftreten. Diese Feststellung ist von Bedeutung, weil viele in der Praxis auftretenden tridiagonalen Matrizen diese Eigenschaft der Diagonal-Dominanz haben!

4. Viele Programmbibliotheken enthalten Programme für die Lösung linearer Gleichungssysteme mit tridiagonalen Koeffizientenmatrizen: siehe z.B. [2], S.304 und [9], S.42.

Weitere Sonderformen der Koeffizientenmatrix.

In der Literatur findet man noch eine große Anzahl von Algorithmen bzw. Programmen für die numerische Behandlung von linearen Gleichungssystemen, deren Koeffizientenmatrizen besondere Strukturen haben. Insbesondere soll auf [2] und auf [9] hingewiesen werden, wo man Informationen zu den Themen: symmetrische Matrizen (Cholesky-Verfahren), zyklisch-tridiagonale Matrizen, fünfdiagonale Matrizen, allgemeine Bandmatrizen, Blockmatrizen usw. erhält.

Das Gauss-Seidel-Verfahren.

Allgemeines.

Das Gauss-Seidel-Verfahren (sowie seine Variationen, die in der Literatur unter den Namen Einzelschrittverfahren, Gesamtschrittverfahren, Relaxationsverfahren, SOR-Verfahren usw. zu finden sind) ist ein iteratives Verfahren zur Lösung von linearen inhomogenen Gleichungssystemen.

Wie schon zu Beginn dieses Kapitels erwähnt wurde, weisen iterative Verfahren gegenüber den direkten Verfahren einige Vorteile (Einfachheit des Algorithmus usw.), aber auch Nachteile (eventuell Konvergenzprobleme) auf.

Der wichtigste Vorteil der iterativen Verfahren liegt aber wohl in der Tatsache, daß die gegebene Koeffizientenmatrix während des Iterationsprozesses nicht verändert wird. Dieser Vorteil wirkt sich - wie im folgenden gezeigt wird - ganz besonders positiv aus, wenn man lineare Gleichungssysteme von hoher Ordnung auszuwerten hat, die schwach besetzte Koeffizientenmatrizen haben.

Eine Matrix heißt schwach besetzt, wenn nur wenige ihrer Koeffizienten von Null verschieden sind.

Da sich einige sehr wichtige Probleme der Numerischen Mathematik (z.B. die Lösung von Randwertproblemen mittels der Differenzenmethode) auf lineare Gleichungssysteme mit schwach besetzten Matrizen zurückführen lassen, soll im folgenden vor allem auf derartige Probleme eingegangen werden.

Die Grundprinzipien des Gauss-Seidel-Verfahrens.

Ausgangspunkt ist wieder das lineare System (2.1):

$$\begin{array}{c} a_{11}x_{1}+a_{12}x_{2}+\cdots+a_{1n}x_{n... ...a_{n1}x_{1}+a_{n2}x_{2}+\cdots+a_{nn}x_{n}=f_{n} \end{array}$$

Unter der Voraussetzung, daß alle Diagonalelemente der Koeffizientenmatrix ungleich Null sind, kann nun jede Gleichung des Systems nach der zum jeweiligen Diagonalelement gehörenden x-Komponente aufgelöst werden:

$$x_{1}=-\frac{1}{a_{11}} \left[ a_{12}x_{2}+a_{13}x_{3}+\cdots+a_{1n}x_{n} -f_{1} \right]$$

$$x_{2}=-\frac{1}{a_{22}} \left[ a_{21}x_{1}+a_{23}x_{3}+\cdots+a_{2n}x_{n} -f_{2} \right]$$

bzw. allgemein

$$x_{i}=-\frac{1}{a_{ii}} \left[ \sum_{j=1(j\ne i)}^{n} a_{ij}x_{j} - f_{i} \right] \quad .$$ (2.18)

Das so umgeformte Gleichungssystem kann als Matrixgleichung

$${\bf x} = C\cdot {\bf x} + {\bf b}$$

mit

$$C = \left\{ c_{ij}\right\}\quad\mbox{mit}\quad c_{ij}=\left\... ...y}\right. \qquad \mbox{sowie}\qquad b_{i}=\frac{f_{i}}{a_{ii}}$$ (2.19)

geschrieben werden. Mit der trivialen Umformung von Glg. (2.18)

$$x_{i}=x_{i}-\left( x_{i}+\frac{1}{a_{ii}} \left[\sum_{j=1(j\ne i)}^{n} a_{ij} x_{j}-f_{i} \right] \right)$$

erhält man sofort die Iterationsvorschrift für das Gauss-Seidel-Verfahren. Setzt man nämlich in die rechte Seite der obigen Gleichung irgendwelche Näherungswerte für die $x_{i}$ ein, so ergeben sich - Konvergenz vorausgesetzt - verbesserte Näherungswerte für die Komponenten des Lösungsvektors:

$$x_{i}^{(t+1)}=x_{i}^{(t)}-\Delta x_{i}^{(t)} \qquad (i=1,\ldots,n)$$ (2.20)

mit

$$\Delta x_{i}^{t}=x_{i}^{(t)} + \frac{1}{a_{ii}} \left[\sum_{j=1(j\ne i)}^{n} a_{ij} x_{j}^{(t)} - f_{i} \right] \quad .$$ (2.21)

Auf diese Weise kann im Falle einer Konvergenz eine Folge von Näherungsvektoren

$${\bf x}^{(0)} \quad ; \quad {\bf x}^{(1)} \quad ; \quad {\bf x}^{(2)} \quad ; \quad \ldots \quad ; \quad {\bf x}^{(t)}$$

für den exakten Lösungsvektor ${\bf x}$ berechnet werden, wobei man diesen zwar niemals erhält, sich aber - abgesehen von Rundungsfehlern - beliebig nahe an ihn 'heranarbeiten' kann:

$$\lim_{t \to \infty} {\bf x}^{(t)} \rightarrow {\bf x}$$

Der Vektor ${\bf x}^{(0)}$ heißt der Startvektor der Iteration.

Typisch für alle Iterationsverfahren sind die sogenannten Abbruchbedingungen, die festlegen, wann eine Iteration zu beenden ist.

In jedem Iterationsprogramm sollten zumindest die beiden folgenden Abbruchbedingungen enthalten sein:

• Die Iteration ist zu beenden, wenn die Ergebnisse eine vorgegebene Genauigkeit erreicht haben.
• Die Iteration ist in jedem Fall zu beenden, wenn eine vorgegebene maximale Zahl von Iterationsschritten erreicht ist ('Notausgang' im Falle einer Divergenz bzw. einer sehr langsamen Konvergenz).

Das Gauss-Seidel-Verfahren für Bandmatrizen.

Eine Bandmatrix ist eine Matrix, die nur entlang einiger Diagonalen Komponenten enthält, die von Null verschieden sind (s. Abbildung 2.2 (oben)).

Es ist sofort einzusehen, daß es vom Standpunkt der Speicherplatz-Ökonomie ein Nonsens ist, eine solche Matrix als ($n$ x $n$)-Matrix abzuspeichern. Die gesamte Information kann nämlich viel platzsparender abgespeichert werden, wenn man die Diagonalen als Vektoren speichert und über die Lage der einzelnen Diagonalen innerhalb der Matrix buchführt.
Wie dies z.B. geschehen kann, ist in der Abbildung 2.2 (unten) schematisch dargestellt:

Abbildung 2.2: Zur Abspeicherung einer Bandmatrix.

$z$ sei die Anzahl der in der Matrix vorkommenden Diagonalen. Die Lage jeder Diagonale wird nun durch eine ganze Zahl angegeben, die angibt, wo sich diese Diagonale relativ zur Hauptdiagonale befindet. Damit kann die Abspeicherung der gesamten Information in einem zweidimensionalen Feld $d_{ik}$ erfolgen, wobei der Index $i$ die jeweilige Matrixzeile bedeutet. Der zweite Index ($k$) gibt an, um welche der $z$ Diagonalen es sich handelt. Zusätzlich wird ein INTEGER-Feld $s(k)$ mit den relativen Positionen definiert:

$$s(k) = \left\{ \begin{array}{ccl} 0 & \mbox{Hauptdiagonale... ...\mbox{-te obere(untere) Nebendiagonale} \end{array} \right.$$

Da beim Gauss-Seidel-Verfahren die Elemente der Hauptdiagonale eine besonders wichtige Rolle spielen, soll der Index $k$, der die Hauptdiagonale bezeichnet, als $k_{o}$ hervorgehoben werden. Es gilt also

$$s(k=k_{o})=0$$

Die Zuordnung der Matrixelemente $a_{ij}$ ($i$-te Zeile, $j$-te Spalte) zu den Elementen $d_{ik}$ ist somit gegeben durch

$$a_{ij}=d_{ik} \quad \mbox{mit} \quad j=s(k)+i \quad ,$$ (2.22)

wobei $i=1,\ldots,n$ und $k=1,\ldots,z$.

Auf Grund all dieser Überlegungen ist es kein Problem, die Iterationsvorschrift (2.21) im Hinblick auf (2.22) zu variieren:

$$x_{i}^{(t+1)} = x_{i}^{(t)} - \Delta x_{i}^{(t)}$$ (2.23)

mit

$$\Delta x_{i}^{(t)} = x_{i}^{(t)} + \frac{1}{d_{i,k_{o}}} \le... ...{k=1(k \ne k_{o})}^{z} d_{ik} x_{s(k)+i}^{(t)} - f_{i} \right]$$ (2.24)

Anmerkung: In der Summe in Glg.(2.24)werden nur jene Summanden zugelassen, für welche gilt:

$$0 < s(k)+i \leq n$$

Konvergenzkriterien und Fehlerverhalten.

Es erhebt sich nun die wichtige Frage, ob man es einem gegebenen linearen Gleichungssystem schon a priori 'ansieht', daß es Schwierigkeiten bei der Konvergenz geben wird.

Nun gibt es in der Literatur in der Tat eine ganze Reihe von Konvergenz- Kriterien, d.h. von Bedingungen, die im Falle ihrer Erfüllung die Konvergenz der Gauss-Seidel-Iteration sicherstellen. Solche Kriterien sind zwar hinreichend, aber i.a. nicht notwendig, d.h. auch manche Systeme, die diese Kriterien nicht erfüllen, können mittels des Gauss-Seidel- Verfahrens gelöst werden.

Aus diesem Grund möchte ich hier auf die Angabe von mathematisch formulierten Konvergenz-Kriterien verzichten und mich auf eine in der Praxis bewährte Faustregel beschränken:

Alle linearen Gleichungssysteme, in deren Koeffizientenmatrizen die Hauptdiagonalelemente die anderen Matrixelemente dominieren, haben gute Chancen auf eine Gauss-Seidel-Konvergenz.

Ein konvergentes Verhalten der Iteration bedeutet natürlich, daß der auftretende Verfahrensfehler $\epsilon_{V}$ mit fortschreitender Iteration immer kleiner wird, daß also gilt:

$$\lim_{t \to \infty} \epsilon_{V}^{(t)} \rightarrow 0$$

In der Praxis wird solange iteriert, bis der Korrekturvektor $\Delta {\bf x}$ [Glg. (2.24)] genügend klein geworden ist. Man verwendet in diesem Fall als Abbruchkriterien die absolute Fehlerabfrage

$$\mbox{max} \mid \Delta x_{i} \mid < \mbox{GEN}$$

bzw. die relative Fehlerabfrage

$$\mbox{max} \mid \frac{\Delta x_{i}}{x_{i}} \mid < \mbox{GEN}$$

wobei GEN eine vom Programmbenutzer festgesetzte Fehler-Toleranz ist.

Achten Sie darauf, die Größe GEN nicht zu klein anzusetzen, da im Falle von schlecht konditionierten Problemen die entsprechende Rechengenauigkeit durch auftretende Rundungsfehler nicht erreicht werden kann!

Nun gibt es jedoch die Aussage, daß iterative Methoden i.a. wesentlich geringere Rundungsfehler aufweisen als die direkten Methoden. Es ist insbesondere eine Eigenschaft von Iterationsmethoden, daß jeweils nur die Rundungsfehler des letzten Iterationsschrittes einen Effekt haben: es gibt also keine $\epsilon_{R}$-Kumulation über alle Iterationen hinweg! Auf Grund dieser Tatsache ist die Verwendung des Gauss-Seidel-Verfahrens in manchen Fällen auch dann einem direkten Verfahren vorzuziehen, wenn es höhere Rechenzeiten erfordert.

Das Unterprogramm GAUSEI.

Das Programm GAUSEI (GAUss-SEIdel) berechnet den Lösungsvektor eines inhomogenen linearen Gleichungssystems mit einer allgemeinen Bandmatrix als Koeffizientenmatrix.

INPUT-Parameter:

N:

Ordnung des Systems.

NDIAG:

Anzahl der in der Bandmatrix vorkommenden Diagonalen.

S( ):

INTEGER-Feld mit den relativen Positionen der Diagonalen.

DIAG( , ):

Feld mit den Matrix-Koeffizienten: der erste Index spezifiziert die Matrixzeile, der zweite die Diagonale.

F( ):

inhomogener Vektor des Systems.

TMAX:

maximale Anzahl der Iterationsschritte.

W:

Relaxationsparameter (s. Abschnitt 2.7.6).

IREL:

IREL $\ne$ 1: absolute Fehlerabfrage,
IREL = 1: relative Fehlerabfrage.

GEN:

absoluter oder relativer Fehler, der im Laufe der Iteration zu unterschreiten ist.

OUTPUT-Parameter:

LOES( ):

Lösungsvektor.

T:

Anzahl der von GAUSEI durchgeführten Iterationsschritte.

FEHLER:

logische Variable zur Fehlerdiagnostik:

FEHLER hat nach der Exekution von GAUSEI den Wert 'falsch', wenn die geforderte Genauigkeit erreicht wurde, und den Wert 'wahr', wenn

• nicht alle Werte der Hauptdiagonalen von Null verschieden sind.
• innerhalb von TMAX Iterationsschritten die geforderte Genauigkeit nicht erreicht wurde.

wichtige interne Variable:

K0:

Index der Hauptdiagonalen.

DX:

Iterations-Korrekturwert gemäß Glg.(2.24).

ISCH:

Steuervariable bei der Genauigkeitsüberprüfung.

Programmstruktur von GAUSEI:

1. Check, welche der übergebenen Diagonalen die Hauptdiagonale ist. Der entsprechende Index $k$ wird als K0 abgespeichert. Ist keine der übergebenen Diagonalen eine Hauptdiagonale, wird das Programm abgebrochen.
2. Check, ob die Hauptdiagonale Nullen enthält. Wenn ja, wird das Programm abgebrochen. Gleichzeitig wird der Startvektor der Iteration (= Nullvektor) definiert.
3. Durchführung der Gauss-Seidel-Iteration mit Fehlerabfrage.

=15cm Struktogramm 8GAUSEI(N,NDIAG,S,DIAG,F,TMAX,W,IREL,GEN,
LOES,T,FEHLER)FEHLER:= .false.
K0:=0K=1(1)NDIAGS(K) = 0K0:=K......K0 = 0FEHLER:= .true.
(return 'Keine Hauptdiagonale')......I=1(1)NLOES(I):=0.0DIAG(I,K0) = 0.0FEHLER:= .true.
(return 'Hauptdiagonale enthaelt Null')......T:=0ISCH:=0I=1(1)NS1:=0.0K=1(1)NDIAGJ:=S(K)+IK $\ne$ K0 .and. J $>$ 0 .and. J $\le$ NS1:=S1+DIAG(I,K)*LOES(J)......DX:=W*((S1-F(I))/DIAG(I,K0)+LOES(I)LOES(I):=LOES(I)-DXIREL = 1ERR:=$\vert$ DX/LOES(I) $\vert$ERR:= $\vert$ DX $\vert$ERROR $>$ GENISCH:=1......T:=T+1T $>$ TMAX .or. ISCH=0ISCH=0print: 'Keine Konvergenz'
FEHLER:= .true.......(return)

Eine Variation des Gauss-Seidel-Verfahrens.

Die Iterationsvorschrift (2.23)

$${\bf x}^{(t+1)} = {\bf x}^{(t)} - \Delta {\bf x}^{(t)}$$

kann durch die formale Einführung eines Faktors $\omega$ im Sinne von

$${\bf x}^{(t+1)} = {\bf x}^{(t)} - \omega \cdot \Delta{\bf x}^{(t)}$$ (2.25)

erweitert werden. Im Falle von $\omega=1$ liegt also die normale Gauss-Seidel-Iteration vor, bei $\omega>1$ spricht man von einer Überrelaxationsmethode, bei $\omega<1$ von einer Unterrelaxationsmethode.

Durch die richtige Wahl des Relaxationsparameters $\omega$ kann in manchen Fällen die Konvergenzgeschwindigkeit außerordentlich erhöht werden. Ein gutes Beispiel für diesen Effekt findet man in [7], S.192ff. Iteriert man das einfache lineare System

$$\begin{array}{c} x+2y=3 x-4y=-3 \end{array}$$

bis zum Unterschreiten des absoluten Fehlers GEN=10$^{-8}$ (C, double precision), so erhält man die folgende Abhängigkeit der erforderlichen Zahl von Iterationsschritten $t$ von $\omega$:

$\omega$	$t$
$0.65$	$20$
$0.70$	$18$
$0.75$	$15$
$0.8$	$14$
$0.85$	$12$
$0.9$	$12$
$0.95$	$21$
$1.0$	$31$
$1.05$	$48$

Es wäre eine enorme Effizienzsteigerung der Gauss-Seidel-Methode, wenn man den idealen Relaxationsparameter $\omega_{{\rm ideal}}$ a priori bestimmen könnte. Zu diesem Thema folgen nun einige Anmerkungen:

• Die Relaxationsmethode ist nur konvergent für $0\leq\omega \leq 2$.
• 'Unter gewissen mathematischen Restriktionen', die i.a. von Matrizen erfüllt werden, welche bei der Anwendung von Differenzenverfahren (s. Kap. 9) zur Lösung von (insbesondere elliptischen) Randwertproblemen auftreten, liegt $\omega_{{\rm ideal}}$ stets im Bereich zwischen 1 und 2, d.h. es liegt immer eine Überrelaxation vor.
• In solchen Fällen gilt:
  

  $$\omega_{{\rm ideal}} = \frac{2}{1+\sqrt{1-\lambda_{1}^{2}}} ,$$ (2.26)

  
  
  wobei $\lambda_{1}$ den größten Eigenwert der Matrix $C$ [Glg. (2.19)] darstellt.
• Da jedoch die genaue Berechnung von $\lambda_{1}$ häufig zu aufwendig ist, verweise ich auf einen Vorschlag in Ref. [14], S. 63, der sich in der Praxis gut bewährt hat und sehr einfach ins Programm GAUSEI eingebaut werden kann:

  Man kann zeigen, daß der Ausdruck
  

  $$\frac{\vert\Delta {\bf x}^{(t)}\vert} {\vert\Delta {\bf x}^{(t-1)}\vert}$$ (2.27)

  
  
  für $t\to\infty$ den Grenzwert $\lambda_{1}^{2}$ hat. Wenn man z.B. im Programm GAUSEI die Iteration mit $\omega$ = 1 beginnt und die ersten $t_{0}$ Iterationen durchführt^2.1, so kann man mittels Glg. (2.27) in vielen Fällen schon eine ausreichend gute Approximation für $\lambda_{1}^{2}$ und damit mittels Glg (2.26) auch für $\omega_{{\rm ideal}}$ erhalten. Mit diesem $\omega$-Wert wird dann die Iteration fortgesetzt.

Wie gut diese Strategie funktioniert, wird in den Übungen zu dieser Vorlesung demonstriert.

Effizienz des Gauss-Seidel-Verfahrens

In den vorigen Abschnitten wurde als Vorteil der Gauss-Seidel-Methode angegeben, daß die Koeffizientenmatrix sich während der Iteration nicht ändert, was vor allem im Falle schwach besetzter Bandmatrizen günstig ist, weil man nur die besetzten Diagonalen abspeichern muß.

Dies soll an Hand eines konkreten Beispiels dokumentiert werden, und zwar für die numerische Auswertung einer sehr bedeutenden partiellen Differentialgleichung, nämlich der Laplace-Gleichung. Unter Verwendung der Differenzenmethode wird diese Differentialgleichung samt ihren Nebenbedingungen in ein System von linearen Gleichungen umgewandelt. Die Ordnung dieses Systems ergibt sich dabei aus der Zahl der Stützpunkte im Grundgebiet der zu lösenden Laplace-Gleichung. Die folgende Abbildung 2.3 (a) zeigt, wie viele Prozent der Elemente der Koeffizientenmatrix als Funktion der Stützpunktzahl tatsächlich besetzt sind. Man sieht, daß dieser Wert ab etwa 200 Stützpunkten deutlich unter 10 Prozent liegt.

Wie sich das auf die im Arbeitsspeicher des Computers nötigen Speicherplätze auswirkt, ist in Abb. 2.3 (b) dargestellt.

Abbildung 2.3: (a) Prozentueller Anteil der besetzten Matrixplätze als Funktion der Zahl der Stützpunkte. (b) Erforderliche Speicherplätze als Funktion der Zahl der Stützpunkte.

Überhuber ([22], S. 392) gibt ein Beispiel aus der industriellen Anwendung linearer Methoden (Karosserie-Festigkeitsberechnung im Automobilbau): Eine Matrix BCSSTK32 hat die Ordnung 44609, also $2\cdot 10^9$ Elemente; von diesen sind aber 'nur' 1029655 Nicht-Nullelemente, d. h. der Besetzungsgrad beträgt 0.05 %.

.

Software-Angebot

Am Ende jedes Hauptkapitels dieses Skriptums finden Sie unter dem Titel 'Software-Angebot' Hinweise auf Programme, die Sie sich entweder kostenlos übers Internet beschaffen können (public domain (p.d.) Programme), oder die kommerziell vertrieben werden.

Quellen: (1)    Bücher von Ch. Überhuber [21], [22].
(2)    Erfahrungen des Autors dieses Skriptums.

Diese Anmerkungen können auf Grund des riesigen Software-Angebotes nur erste Anregungen für Sie sein, dieses Angebot zu nutzen.

Zitat aus [21], S. 322:
Von herausragender Bedeutung für die rasche, einfache und effiziente Beschaffung von frei erhältlicher (public domain)-Software aus dem Numerik-Bereich ist der Internet-Dienst NETLIB (Dongara, Grosse, 'Distribution of Mathematical Software via Electronic Mail', Comm. ACM 30 (1987), 403).

Auf Programme aus dem NETLIB kann über verschiedene Internet-Dienste zugegriffen werden (email, FTP, Internet); am einfachsten ist ein 'Downladen' aus dem www-System. Ich werde mich im folgenden auf diese Möglichkeit beschränken.
Zuvor noch ein weiteres Zitat ([22], S. 183) mit speziellem Bezug auf das Kap. 2 dieses Skriptums:
Für die numerische Lösung linearer Gleichungssysteme gibt es, mehr als in jedem anderen Gebiet der Numerischen Datenverarbeitung, eine Fülle qualitativ hochwertiger Softwareprodukte: LAPACK, TEMPLATES, ITPACK, NSPCG, SLAP, UMFPACK, Teile der IMSL- und NAG-Bibliotheken. Die eigene Entwicklung von Software zur Lösung linearer Gleichungssysteme wäre aus diesem Grund ein sehr unökonomisches Unterfangen.

• Zugriff auf NETLIB:        www.netlib.org
• Einen ersten Überblick über das Angebot:
  Klick auf Browse the Netlib repository

Die wichtigsten Libraries für den Problemkreis "Lineare Gleichungssystemeßind die folgenden:

LAPACK:

Umfangreiche Sammlung von F77-Programmen zu den Themen 'Lösung linearer Gleichungssysteme', 'Eigenwertprobleme' u.a.
Die komplette aktuelle Bibliothek (31.5.2000) inklusive aller sog. BLAS-Hilfsprogramme^2.2umfaßt ca. 5 MByte und kann durch Anklicken von
lapack.tgz downgeladen werden. Danach einfach auf Ihrem LINUX-Rechner dekomprimieren und in die einzelnen Programme aufteilen:

gunzip -c lapack.tgz | tar xvf -

Achtung: Alle Programme (ohne OBJ- und EXE-Files) in den 4 Datentypen REAL, DOUBLE PRECISION, COMPLEX und COMPLEX*16 benötigen ca. 34 MByte.

• Wenn Sie LAPACK-Programme intensiv verwenden wollen, ist es vernünftig, sich den LAPACK Users' Guide (E. Anderson et al, Society for Industrial and Applied Mathematics, 1999) zuzulegen. Dieses Manual ist auch im Internet unter www.netlib.org/lapack/lug/lapack_lug.html erhältlich.
• Wenn Sie das ganze LAPACK auf Ihren Rechner installieren wollen, ist es auch nützlich, den LAPACK Installation Guide (ps-File, 232 kByte, 25 Seiten) downzuladen.
• Wenn Sie eine bestimmte LAPACK-Routine laden wollen, gehen Sie wie folgt vor: Klicken Sie im LAPACK-Verzeichnis netlib.org/lapack/index.html die entsprechende Eintragung in den Individual Routines an.
  Angenommen, Sie wollen ein DOUBLE-PRECISION-Lösungsprogramm für lineare Gleichungssysteme in Ihren Rechner laden. In diesem Fall klicken Sie auf lib double; als erste Eintragung in der erscheinenden Programmliste finden Sie:

  file dgesv.f   dgesv.f plus dependencies
  prec double
  for  Solves a general system of linear equations AX=B.
  

  Sie könnten nun durch Anklicken von dgesv.f das Programm speichern. Allerdings hat das i.a. wenig Sinn, da die meisten LAPACK-Programme eine ganze Kaskade von weiteren Programmen aufrufen. Es ist daher besser, diese ganze 'Programmfamilie' durch Anklicken von dgesv.f plus dependencies zu laden. Im erscheinenden Dialogfenster wählen Sie als "Packaging format" Tar-gzip, und die nötigen BLAS-Files laden Sie am besten gleich dazu.
  Wenn das geklappt hat, haben Sie in Ihrem Verzeichnis den File
  netlibfiles.tgz. Die weitere Verarbeitung lautet (s.o.):

  gunzip -c netlibfiles.tgz | tar xvf -
  

  Damit haben Sie sofort, auf drei neu erstellte Verzeichnisse aufgeteilt, die gesamte Programmfamilie von dgesv.f zur Verfügung:

  blas/blas2test.f
  blas/dgemm.f
  blas/dger.f
  blas/dscal.f
  blas/dswap.f
  blas/dtrsm.f
  blas/idamax.f
  blas/xerbla.f
  lapack/double/dgesv.f
  lapack/double/dgetf2.f
  lapack/double/dgetrf.f
  lapack/double/dgetrs.f
  lapack/double/dlaswp.f
  lapack/util/ieeeck.f
  lapack/util/ilaenv.f
  

Die folgende knappe Zusammenstellung gibt Ihnen eine Übersicht über das LAPACK-Angebot ([22], S. 266). Die Begriffe simple drivers und expert drivers bedeuten:
simple drivers:    black box-Programme, die primär die numerische Lösung der jeweiligen Problemstellung liefern.
expert drivers    liefern neben den Lösungen auch noch zusätzliche Informationen (Konditionsabschätzungen, Fehlerschranken usw.).

LAPACK bietet nur direkte Verfahren in FORTRAN-77. Diese Programmiersprache ist zwar immer noch weitverbreitet im technisch-naturwissenschaftlichen Bereich; da aber neue Fortranversionen wie F90, F95 sowie andere Sprachen wie C und C++ immer wichtiger werden, gibt es inzwischen auch schon die entsprechenden LAPACK-Routinen:

LAPACK95

stellt ein Interface zu F77 dar, d. h. Sie benötigen auch die LAPACK-F77 Programme. LAPACK95 erhalten Sie über
www.netlib.org/lapack95/index.html.
Downladen von lapack95.tgz und nachfolgende Dekomprimierung erzeugt in Ihrer Directory das Unterverzeichnis LAPACK95. In diesem gibt es das Verzeichnis SRC mit allen zur Verfügung stehenden Routinen.

CLAPACK

erhalten Sie über www.netlib.org/clapack.
Es handelt sich dabei um eine vollständige C-Implementierung von LAPACK, wobei die Konvertierung von F77 in C automatisch erfolgte. Der Nachteil dieses Vorgehens ist allerdings, daß CLAPACK wertvolle Eigenschaften von C, wie z. B. dynamische Speicherallokation, nicht nutzt.
Darüber hinaus können Sie nicht einfach irgendein einzelnes CLAPACK-Programm herunterladen, sondern es ist eine komplette Installation von
clapack.tgz (Version 3.0, ca. 6 MByte) erforderlich. Dieser File enthält auch die Installationshilfen
readme.install    und    readme.maintain     sowie den Header-File clapack.h

LAPACK++

stellt eine auf C++ basierende objektorientierte Schnittstelle zu LAPACK dar. Der Sourcecode heißt lapack++.tgz und enthält in komprimierter Form 66 kByte. Genauere Information über dieses System sind in den folgenden Postcript-Files zu finden:
overview.ps        classref.ps        install.ps        user.ps

Eine weitere interessante Internet-Adresse zum Thema LAPACK++:
http://math.nist.gov/lapack++

SCALAPACK

Anwendung von verteilten und parallelen Rechnersystemen zur Lösung von Gleichungssystemen, Ausgleichsproblemen und Eigenwertproblemen mit sehr großen Koeffizientenmatrizen ($n > 5000$).
Eine kurze und prägnante Einführung in dieses Thema (Stand: 1995) s. [22], S. 283-286.

Mathematica und MATLAB:

Mathematica erlaubt sowohl symbolische als auch numerische Behandlung wichtiger Probleme der Linearen Algebra, wie z.B.:

Inverse[m]        Matrix-Inversion
Det[m]        Determinante
LinearSolve[m,b]        Lösung des inhomogenen linearen Systems.

MATLAB:^2.3 im Prinzip bietet diese Benutzersoftware die einfachste Realisierung der Lösung des linearen Gleichungssystems A.x=f :

>>% Eingabe der Koeffizientenmatrix A:
>> A=[ ... ; ... ];
>>% Eingabe des inhomogenen Vektors f:
>> f=[ ... ],
>>% Loesung des Gleichungssystems:
>> x=A\f

Eine alternative Lösungsmöglichkeit bietet Matlab durch folgende Befehle:

>>% Der inhomogene Vektor wird zur Koeffizientenmatrix als letzte Spalte
>>% hinzugefuegt:
>> Aplus=[A f];
>>% Diese erweiterte Matrix wird als Argument der Funktion  rref
>>% verwendet. Nun wird die urspruengliche Matrix mittels des Verfahrens
>>% von Gauss-Jacobi in eine Einheitsmatrix umgewandelt, und die letzte
>>% Spalte enthaelt den Loesungsvektor:
>> Aneu=rref(Aplus);
>> x=Aneu(:,end)
>>% Dieses Verfahren hat den Vorteil, dass es etwas stabiler gegenueber
>>% Rundungsfehlern ist als die obige "einfache" Loesung.

Iterative Lösung schwach-besetzter Systeme

TEMPLATES

Diese NETLIB-Library enthält eine Reihe von Programmpaketen mit iterativen Lösungsmethoden für lineare Gleichungssysteme:

Dateiname	Inhalt der Datei
sctemplates.tgz	C-Routinen, einfache Genauigkeit
dctemplates.tgz	C-Routinen, doppelte Genauigkeit
cpptemplates.tgz	C++-Routinen
sftemplates.tgz	F77-Routinen, einfache Genauigkeit
dftemplates.tgz	F77-Routinen, doppelte Genauigkeit
mltemplates.tgz	Matlab-Routinen

Eine interessante Matlab-Routine ist sor.m. Dieses Programm löst ein lineares Gleichungssystem mittels der Sukzessiven Überrelaxationsmethode (allerdings ohne Berücksichtigung einer eventuellen Bandstruktur der Koeffizientenmatrix).