Lineare Gleichungssysteme

Einführendes Beispiel - Astrometrie

Um die Positionen von Sternen und anderen Himmelskörpern zu bestimmen, benötigt man sogenannte Plattenkonstanten. Heute sind das natürlich keine Fotoplatten mehr, sondern elektronische Aufnahmechips. Es heißt aber immer noch Plattenkonstanten. Worum geht es?

Um die Koordinaten eines Kometen an Hand bekannter Sterne zu bestimmen, muss man die Plattenkonstanten kennen. Foto: Michael Jäger

Man muss von zweierlei ausgehen: 1) Das Bildfeld ist nicht parallel zum Meridian ausgerichtet 2) Die Brennweite der Optik ist nicht überall konstant. Was man herausbekommen möchte sind Koordinaten X und Y, die parallel zum Meridian sind und nur von der nennbrennweite abhängen. Die Eigenheiten der Optik sind aber überhaupt nicht formelmäßig bekannt, so dass man mit einer Näherung arbeiten muss. Eine solche Näherung bezeichnet man als "Ansatz", weil man einfach annimmt (ansetzt), dass die wirklichen Gegebenheiten so sind. Als Ansatz für die Koordinate X wird benutzt:

X = a · x + b · y + c

Für Y macht man einen ähnlichen Ansatz. Für uns genügt jetzt diese eine Koordinate.
Die Werte x und y sind die auf der Platte gemessenen Koordinaten sind. Eigentlich sollte X nur von x abhängen. Wenn aber das Bildfeld etwas verdreht ist, dann hängt es eben doch ein wenig von y ab. Dieser Ansatz hat drei unbekannte Koeffizienten, a, b und c. Um diese zu ermitteln, kann man drei bekannte Sterne messen, deren X-Koordinate man aus einer Tabelle entnimmt. Für jeden Stern ist X also bekannt, und x und y kann man messen. Man erhält dann drei Gleichungen

X₁ = a · x₁ + b · y₁ + c
X₂ = a · x₂ + b · y₂ + c
X₃ = a · x₃ + b · y₃ + c

Hiermit kann man a, b, und c berechnen. Das einfache Einsetzverfahren erfordert, dass man z.B. die erste Gleichung nach a umstellt, das alles in die zweite Gleichung einsetzt (damit gibt es dort kein a mehr). Die zweite Gleichung stellt man nach b um. In die dritte Gleichung setzt man zuerst die umgestellte erste Gleichung ein, damit entfällt das a. Dann setzt man für jedes b die zweite Gleichung ein, damit entfällt auch das b und wir haben nur noch c. Jetzt stellen wir nach c um und können das ausrechnen. Dieses c setzen wir in die umgestellte zweite Gleichung ein, rechnen b aus. b und c setzen wir in die umgestellte erste Gleichung ein und rechnen a aus.
Genauso umständlich, wie es klingt, ist es auch. Bei zwei Unbekannten kann man das noch machen, bei drei ist es fraglich, über drei Unbekannte muss etwas anderes her.

Schematisierung - Speicherung in einer Matrix

Eine Matrix ist ein rechteckiges Rechentableau, weiter nichts. In unserem Fall enthält die Matrix die gemessenen Koordinaten der Referenzsterne und als Ergebnisspalte die wirklichen X-Werte. Die zu berechnenden Faktoren a, b, un c erden nicht mitgeschrieben, ich stelle sie aber in die Überschrift:

a       b      c
___________
x₁y₁ 1    X₁
x₂y₂1    X₂
x₃y₃1    X₃

Die erste Zeile der Matrix besagt also a · x₁ + b · y₁ + c · 1 = X₁, also die erste Gleichung.

Mit dieser Matrix können wir nun Rechenoperationen ausführen. Diejenigen, welche wir benötigen, beziehen sich immer auf ganze Zeilen:

Man kann zwei Zeilen vertauschen. Das bedeutet nichts anderes, als hätten wir die Gleichungen in einer anderen Reihenfolge hingeschrieben
Man darf eine Zeile mit einem Faktor multiplizieren. Dabei würden wir die linke und rechte Seite der Gleichung mit demselben Faktor malnehmen, die Gleichung bliebe immer noch wahr.
Man darf eine Zeile zu einer anderen addieren. Das ist vielleicht nicht so intuitiv verständlich. Aber wenn jede linke Seite gleich der rechten Seite ist, dann sind die Summen zweier linken Seiten auch gleich der Summe zweier rechten Seiten.

Der Gauss-Algorithmus zur Lösung eines linearen Gleichungssystems

Mit Hilfe dieser erlaubten Operationen wird nun folgendes durchgeführt:

Man benutzt die erste Zeile dazu, um in der ersten Spalte überall Nullen zu erzeugen, mit Ausnahme dieser ersten Zeile selbst. Dazu muss man diese erste Zeile mit x₂ / x₁ multiplizieren und von der zweiten Zeile abziehen. Multiplizieren und subtrahieren geschieht in einem Schritt, die Zeile 1 brauchen wir nicht zu verändern. Ähnliches passiert mit der dritten Zeile.
Man benutzt die zweite Zeile, um in der zweiten Spalte alle Elemente zu Null zu bringen, wieder mit Ausnahme dieser zweiten Zeile
Dasselbe geschieht für die letzte Zeile.

Durch Anwendung erlaubter Rechenoperationen haben wir dann eine Matrix, die so aussieht:

a          b          c
__________________
m_{1,1    0          0}           m_3,1
0m_{2,2     0}          m_4,2
00m_3,2 m_4,3

Die Elemente m_1,1 usw. sind die nach der Rechnung in der Matrix stehenden Element, also irgendwelche Zahlen, welche die Rechnung ergab. In das "Gleichungsdeutsch" übersetzt heißt die erste Zeile a · m_1,1= m_3,1. Damit haben wir a, aus der zweiten Zeile b und aus der dritten c. Im Studium habe ich so was per Hand berechnet, bis hin zu Gleichungssystemen vierter Ordnung. Weniger fehlerträchtig ist das mit dem Computer. Ich gebe deshalb ein kleines Python-Programm an, welches auch ein Beispiel enthält:

def gauss(n, m): 
    for i in range (n): # Pivotelement in der Diagonalen
        for j in range (n): # jede Zeile mit diesem Pivot-E.
            f=m[j][i]/m[i][i]
            for k in range(i,n+1): # alle Elemente der Zeile korrigieren
                if (i!=j):
                    m[j][k]=m[j][k]-f*m[i][k]

# Hauptprogramm
n=3 # Grad des LGS
m = [[0 for x in range(n+1)] for y in range(n)] 
# Matrix befüllen
m[0][0]=1; m[0][1]=1;   m[0][2]=1; m[0][3]=6
m[1][0]=1; m[1][1]=1.5; m[1][2]=2; m[1][3]=10
m[2][0]=2; m[2][1]=1;   m[2][2]=2; m[2][3]=10
gauss(n, m)
for i in range (n):
    print(i, m[i][n]/m[i][i])

In diesem Programm enthält die Beispiel-Matrix nacheinander die folgenden Werte:

1 1   1  6       1  1   1  6         1 0   -1 -2       1 0   0 1
1 1.5 2 10       0  0.5 1  4         0 0.5  1  4       0 0.5 0 1
2 1   2 10       0 -1   0 -2         0 0    2  6       0 0   2 6

Das Ergebnis is a=1, b=2 und c=3.

Verbesserungen dieses Verfahrens

Das zum Nullsetzen benutzte Diagonalelement bezeichnet man auch als Pivot-Element und die betreffende Zeile Pivot-Zeile. Das rührt vom französischen pivot = Angelpunkt her. Es kann nun vorkommen, dass ein Pivotelement gerade den Wert 0 erhält. n diesem Fall scheitert der Algorithmus. Das Ergebnis wird unnötig ungenau, wenn ein Pivotelement einen sehr kleinen Wert enthält. Es ist besser, als Pivot-Zeile eine andere Zeile zu benutzen, das ist ja zulässig. Ich gebe hier eine Modifikation an, welche alle noch zulässigen Zeilen untersucht, das betragsgrößte Pivotelement auswählt und durch Zeilen-vertauschen an die richtige Stelle versetzt.

def gauss(n, m): 
    for i in range (n): # Pivot-Zeile.
        # bestes Pivotelement suchen
        max=0; k=0
        for j in range(i,n): # benutzte Zeilen auslassen = ab i
            if abs(m[j][i])>max:
                   max=abs(m[j][i]); k=j
        # Zeilen vertauschen, so dass Pivot-E = Diagonal-E.
        for j in range(n+1):
            h=m[i][j]
            m[i][j]=m[k][j]
            m[k][j]=h
        # weiter wie in gewöhnlichen Gauss-Algorithmus
        for j in range (n): # jede Zeile mit diesem Pivot-E.
            f=m[j][i]/m[i][i]
            for k in range(i,n+1): # alle Elemente der Zeile korr.
                if (i!=j):
                    m[j][k]=m[j][k]-f*m[i][k]

Das Hauptprogramm ist gleich. Die Matrix nimmt jetzt nacheinander die folgenden Werte an:

1 1   1  6       2 1   2 10         2 0  1    5         2 0  0    2
1 1.5 2 10       0 1   1  5         0 1  1    5         0 1  0    2
2 1   2 10       0 0.5 0  1         0 0 -0.5 -1.5       0 0 -0.4 -1.5

Mit dieser Art Pivotwahl kann man durchaus einigermaßen gutartige Gleichungssysteme mit 100 Unbekannten lösen.

(Uwe Pilz, Oktober 2019)