13    Ausgleichsrechnungen

In diesem Kapitel lernen Sie, wie Sie aus Messdaten Parameter für das lineare und nichtlineare Ausgleichsproblem mit linearen Gleichungssystemen und dem Vektoransatz lösen können.

In der Praxis der Natur- und Ingenieurwissenschaften kommen Aufgabenstellungen vor, bei denen aus vorliegenden Messdaten Parameter berechnet werden müssen, die den Verlauf einer bekannten mathematischen Funktion bestimmen. Beispielsweise weiß der Astronom, dass sich Planeten auf elliptischen Bahnen bewegen. Will er die beiden Halbachsen der elliptischen Umlaufbahn eines Planeten bestimmen, so muss er zu unterschiedlichen Zeitpunkten die zugehörigen Planetenpositionen messen, diese Daten in x,y-Koordinaten umrechnen und aus diesen Wertepaaren die gesuchten Parameter berechnen. Da Messdaten prinzipiell fehlerbehaftet sind, muss er mit der Methode der Ausgleichsrechnung die Parameter so bestimmen, dass der Funktionsgraph möglichst genau an die Wertepaare der Messdaten angepasst wird.

Ein anderes Beispiel für eine nützliche Anwendung ist die Prognose eines Bremsweges. Für kleine Geschwindigkeiten sind die Parameter der mathematischen Bremsfunktion bekannt. Die Frage lautet nun: Welcher Bremsweg ist bei höheren Geschwindigkeiten zu erwarten?

Neben diesen linearen Ausgleichsproblemen gibt es auch nichtlineare Ausgleichsprobleme. Dazu zählt z. B. das exponentielle Wachstum der Population einer Bakterienkultur oder die Entwicklung einer Epidemie. Das exponentielle Wachstum kann anhand der Wachstumsrate vorausgesagt werden, wenn sich die Wachstumsbedingungen nicht ändern. Aus den Messdaten des bisherigen Infektionsgeschehens kann mit der Ausgleichsrechnung die Wachstumsrate bestimmt werden, um damit den zukünftigen Verlauf der Epidemie vorauszusagen.

13.1    Lineare Ausgleichsprobleme

Nehmen Sie an, dass Messdaten in Form von Wertepaaren (x, y) vorliegen. Mit der Methode der kleinsten Fehlerquadrate wird aus diesen Messdaten ein lineares Gleichungssystem mit n Unbekannten aufgestellt. Bei den Unbekannten handelt es sich um die Parameter einer bekannten Funktion. Bei der Lösung eines linearen Ausgleichsproblems müssen also lineare Gleichungssysteme aufgestellt und gelöst werden.

Neben dieser Methode gibt es noch die Möglichkeit, direkt aus den Messdaten ein überbestimmtes Gleichungssystem aufzustellen, das mit Matrizenoperationen, wie Transponieren, Bildung der Inversen und Multiplikation, gelöst werden kann.

13.1.1    Lösung mit linearen Gleichungssystemen

Am Beispiel einer linearen Funktion und eines Polynoms 2. Grades soll gezeigt werden, wie für die vorgegebenen Messdaten ein lineares Gleichungssystem aufgestellt werden muss und wie dieses Gleichungssystem mit der NumPy-Funktion solve() gelöst werden kann. Ich kann jedoch bereits vorwegnehmen, dass diese Methode sehr aufwendig wird, wenn mehr als drei Parameter berechnet werden müssen. Der Vektoransatz, bei dem überbestimmte Gleichungssysteme gelöst werden müssen, ist dagegen wesentlich einfacher zu handhaben, verlangt aber einen größeren maschinellen Rechenaufwand.

Ausgleichsgerade

Bei der Bestimmung einer linearen Ausgleichsgeraden geht es darum, für eine Messreihe eine Gerade der Form

zu finden. Diese Gerade soll zu allen Messpunkten jeweils den kleinsten möglichen Abstand haben (siehe Abbildung 13.1).

Betrachten Sie die Messreihe aus Tabelle 13.1.

xi	1	2	3	4	5
yi	6	6,8	10	10,5	10,2

Tabelle 13.1     Wertetabelle für einen linearen Zusammenhang

Ein erster Blick auf die Messwerte lässt die Vermutung zu, dass hier ein annähernd linearer Zusammenhang zwischen x und y besteht.

Die Fehlerquadratsumme

liefert den Ansatz zur Bestimmung der Parameter a und b. Sie wird auch als Fehlerfunktional [Knorrenschild: 94] bezeichnet.

Um den kleinsten Wert der Fehlerquadratsumme zu bestimmen, müssen die partiellen Ableitungen und gleich 0 gesetzt werden:

Daraus erhalten Sie ein lineares Gleichungssystem mit den beiden unbekannten Parametern a und b:

Dieses Gleichungssystem wird als das System der Normalgleichungen bezeichnet. Die 2×2-Matrix bezeichnet man auch als Koeffizientenmatrix. Der Spaltenvektor (a, b)^T heißt Lösungsvektor, und der Spaltenvektor auf der rechten Seite des Gleichungssystems wird als Inhomogenitätsvektor bezeichnet.

Wenn Sie die Werte aus Tabelle 13.1 in die Gleichung einsetzen, erhalten Sie das lineare Gleichungssystem

mit den Lösungen a = 1,21 und b = 5,07.

Listing 13.1 berechnet aus den Messdaten die Parameter a und b der Ausgleichsgerade. Die Messdaten werden dann in ein NumPy-Array np.array([...]) gespeichert. Die NumPy-Methode np.sum() berechnet die Summen. Das lineare Gleichungssystem wird anschließend mit der NumPy-Methode solve(A,b) gelöst.

01  #01_plot_gerade.py
02  import numpy as np
03  from numpy.linalg import solve
04  import matplotlib.pyplot as plt
05  xi=np.array([1,2,3,4,5])
06  yi=np.array([6,6.8,10,10.5,10.2])
07  n=len(xi)
08  #Koeffizientenmatrix
09  A=np.array([[np.sum(xi**2) ,np.sum(xi)],
10              [np.sum(xi),n]])
11  #Inhomogenitätsvektor
12  y=np.array([np.sum(xi*yi),np.sum(yi)])
13  #Lösungsvektor
14  L=solve(A,y)
15  a,b=L[0],L[1]
16  E=0 #Fehlerfunktional
17  for i in range(n):
18      E=E + (yi[i] - (a*xi[i] + b))**2
19  #Ausgaben
20  print("Fehlerfunktional\n E =",E)
21  print("Koeffizientenmatrix\n",A)
22  print("Inhomogenitätsvektor\n",y)
23  print("Lösungsvektor\n",L)
24  print("Steigung        a=",a)
25  print("Achsenabschnitt b=",b)
26  #Funktionsplot
27  xa=np.linspace(0,np.max(xi),500)
28  ya=a*xa+b
29  fig, ax = plt.subplots()
30  ax.set(xlabel="x",ylabel="y",title=r"$y=ax+b$")
31  ax.plot(xi,yi,'r+')
32  ax.plot(xa,ya,'b')
33  ax.set_ylim(0,15)
34  plt.show()

Listing 13.1     Linearer Ausgleich

Ausgabe

Fehlerfunktional
 E = 3.439000000000002
Koeffizientenmatrix
[[55 15]
 [15  5]]
Inhomogenitätsvektor
[142.6  43.5]
Lösungsvektor
[1.21 5.07]
Steigung        a= 1.209999999999998
Achsenabschnitt b= 5.070000000000007

Abbildung 13.1     Funktionsplot für den linearen Ausgleich

Analyse

In den Zeilen 02 bis 04 werden die erforderlichen Module importiert. Aus dem Untermodul linalg wird nur die Methode solve() für die Lösung von linearen Gleichungen importiert.

In den Zeilen 05 und 06 werden die Messdaten in die Objektvariablen xi und yi gespeichert.

Zeile 07 ermittelt die Anzahl n der Messdaten.

Die Koeffizientenmatrix A steht in den Zeilen 09 bis 10. Die NumPy-Methode sum(...) berechnet die Summen aus den Messwerten der unabhängigen Variablen xi.

In Zeile 12 werden die Messwerte des Inhomogenitätsvektors in die Objektvariable y gespeichert.

In Zeile 14 löst die Methode solve(A,y) das lineare Gleichungssystem und speichert den Lösungsvektor in die Objektvariable L.

In Zeile 15 werden die Einträge des Lösungsvektors L[0] und L[1] in die Variablen a und b gespeichert.

In den Zeilen 17 und 18 berechnet der Summenalgorithmus das Fehlerfunktional E.

Die numerischen Ausgaben erfolgen in den Zeilen 20 bis 25.

Die Zeilen 31 bis 32 sind für den Funktionsplot der Messpunkte (Zeile 31) und der Ausgleichsgeraden (Zeile 32) zuständig.

Polynom 2. Grades

Für das Polynom 2. Grades

sollen die Parameter bestimmt werden. In dem Funktionsterm kommt das Quadrat der Variablen x vor.

Für ein Polynom 2. Grades lassen sich mit dem Ansatz der Methode der kleinsten Quadrate folgende Normalgleichungen aufstellen [Papula3: 716]:

Für die Messdaten aus Tabelle 13.2 sollen die Parameter a, b und c bestimmt werden.

xi	1	2	3	4	5
yi	5,1	2,1	1	1,8	5,2

Tabelle 13.2     Messwerte für ein Polynom 2. Grades

Listing 13.2 berechnet diese Parameter mit der NumPy-Methode solve(...) aus den Normalgleichungen:

01  #02_plot_parabel.py
02  import numpy as np
03  from numpy.linalg import solve
04  import matplotlib.pyplot as plt
05  xi=np.array([1,2,3,4,5])
06  yi=np.array([5.1,2.1,1,1.8,5.2])
07  n=len(xi)
08  #Koeffizientenmatix
09  A=np.array([[np.sum(xi**4),np.sum(xi**3),np.sum(xi**2),],
10              [np.sum(xi**3),np.sum(xi**2),np.sum(xi)],
11              [np.sum(xi**2),np.sum(xi),n]])
12  #Inhomogenitätsvektor
13  y=np.array([np.sum(xi**2*yi),np.sum(xi*yi),np.sum(yi)])
14  #Lösungsvektor
15  L=solve(A,y)
16  a,b,c=L[0],L[1],L[2]
17  #Ausgabe der numerischen Werte
18  print(" a\t b\tc\n",L)
19  print("%2.3fx^2 %2.3fx %2.3f" %(a,b,c))
20  #Funktionsplot
21  xa=np.linspace(0,np.max(xi+1),500)
22  ya=a*xa**2 + b*xa + c
23  fig, ax = plt.subplots()
24  ax.set(xlabel="x",ylabel="y",title=r"$y=ax^2 + bx +c$")
25  ax.plot(xi,yi,'r+')
26  ax.plot(xa,ya,'b')
27  ax.set_xlim(0,6)
28  ax.set_ylim(0,12)
29  plt.show()

Listing 13.2     Ausgleichsproblem für ein Polynom 2. Grades

Ausgabe

   a     b     c
[ 1.05 -6.31 10.42]
1.050x^2 -6.310x 10.420

Abbildung 13.2     Polynom 2. Grades

Analyse

Das Programm ist im Prinzip genauso aufgebaut wie Listing 13.1. An der Koeffizientenmatrix in Zeile 09 bis 11 zeigt sich, dass der Python-Interpreter einen vergleichsweise höheren Aufwand für die Berechnung der Koeffizientenmatrix und für die Lösung des Gleichungssystems (Zeile 15) zu bewältigen hat.

13.1.2    Lösung mit dem Vektoransatz

Sollen mehr als drei Parameter mit einem linearen Gleichungssystem berechnet werden, dann steigt der Aufwand unverhältnismäßig stark an. Mit dem Matrix-Vektor-Ansatz kann der Aufwand für das Aufstellen der Gleichungssysteme erheblich reduziert werden. Das System der Normalgleichungen kann mit der Matrix-Vektor-Schreibweise direkt aus den Wertetabellen als Koeffizientenmatrix und Inhomogenitätsvektor notiert werden. Allerdings erhält man so ein überbestimmtes lineares Gleichungssystem, das mit einem höheren Rechenaufwand gelöst werden muss.

Für das überbestimmte Gleichungssystem gilt der Ansatz:

Mit dem Lösungsvektor

erhält man das überbestimmte Gleichungssystem in Matrix-Schreibweise:

Dieses Gleichungssystem wird auch als Fehlergleichungssystem [Knorrenschild: 96] bezeichnet. Das Fehlergleichungssystem enthält n Gleichungen und zwei Unbekannte. Es ist also überbestimmt und muss deshalb auf beiden Seiten der Gleichung mit der transponierten Matrix A^T multipliziert werden:

Die Lösungsmenge λ erhält man durch Bildung der Inversen:

Normalgleichung

Man nennt A^TA λ = A^Ty die Normalgleichung von A λ = y. Das lineare Ausgleichsproblem hat eine eindeutige Lösung λ, die durch A^TA λ = A^Ty gegeben ist.

Algorithmus für den Lösungsweg

1. Schritt: Aufstellen der Koeffizientenmatrix

2. Schritt: Transponieren der Koeffizientenmatrix

3. Schritt: Multiplikation der transponierten Koeffizientenmatrix mit der Koeffizientenmatrix A^TA. Es entsteht eine quadratische Koeffizientenmatrix.

4. Schritt: Berechnen des Inhomogenitätsvektors A^Ty

5. Schritt: Berechnen der inversen Matrix (A^TA)^–1

6. Schritt: Berechnen der Lösung λ= (A^TA)^–1 A^Ty

Mit dem Modul NumPy können Sie die erforderlichen Matrizenoperationen besonders einfach ohne großen Programmieraufwand durchführen. Die Matrizenmultiplikation wird am einfachsten mit dem At-Operator @ durchgeführt. Transponiert wird eine Matrix mit A.T. Mit NumPy können Sie die Matrizenmultiplikation zusammen mit der Invertierung inv() in einer Anweisung unterbringen: L=inv(A.T@A)@A.T@y.

Listing 13.3 berechnet mit dem Vektoransatz die Parameter a und b einer linearen Funktion für die Messdaten aus Tabelle 13.1.

01  #03_vektor_gerade.py
02  import numpy as np
03  from numpy.linalg import inv
04  x=np.array([1,2,3,4,5])
05  y=np.array([6,6.8,10,10.5,10.2])
06  n=len(x)
07  A=np.array([x,np.ones(n)]).T
08  L=inv(A.T@A)@A.T@y
09  print("1. Koeffizientenmatrix des überbestimmten Gleichungssystems\n",A)
10  print("2. transponierte Koeffizientenmatrix\n",A.T)
11  print("3. Koeffizientenmatrix der Normalgleichungen\n",A.T@A)
12  print("4. normalisierter Inhomogenitätsvektor\n",A.T@y)
13  print("5. invertierte Koeffizientenmatrix der Normalgleichung\n", inv(A.T@A))
14  print("6. Lösungsvektor\n",L)
15  a,b=L[0],L[1]
16  print("Steigung        a =",a)
17  print("Achsenabschnitt b =",b)

Listing 13.3     Lösen eines überbestimmten linearen Gleichungssystems

Ausgabe

1. Koeffizientenmatrix des überbestimmten Gleichungssystems
 [[1. 1.]
 [2. 1.]
 [3. 1.]
 [4. 1.]
 [5. 1.]]
2. transponierte Koeffizientenmatrix
 [[1. 2. 3. 4. 5.]
 [1. 1. 1. 1. 1.]]
3. Koeffizientenmatrix der Normalgleichungen
 [[55. 15.]
 [15.  5.]]
4. normalisierter Inhomogenitätsvektor
 [142.6  43.5]
5. invertierte Koeffizientenmatrix der Normalgleichung
 [[ 0.1 -0.3]
 [-0.3  1.1]]
6. Lösungsvektor
 [1.21 5.07]
Steigung        a = 1.21
Achsenabschnitt b = 5.070000000000006

Analyse

Zeile 03 importiert die Funktion inv aus dem Unterpaket linalg von NumPy. Diese Funktion invertiert Matrizen.

In Zeile 07 wird die Koeffizientenmatrix mit den Einträgen der x-Werte und dem Einheitsvektor np.ones(n) aufgestellt. Sie muss mit dem Operator T transponiert werden, weil das NumPy-array() nur Zeilenvektoren erzeugt.

Für die Berechnung des Lösungsvektors L ist nur Zeile 08 notwendig. Alle anderen Ausgaben dienen nur dazu, die Abfolge der einzelnen Rechenoperationen zu verdeutlichen.

Der Vektoransatz berechnet den gleichen Wert wie Listing 13.1.

13.1.3    Lösung mit dem erweiterten Vektoransatz

Für lineare Ausgleichsprobleme, bei denen mehr als zwei Parameter berechnet werden müssen, kann allgemein der Ansatz

gewählt werden [Knorrenschild: 96].

Oder kürzer notiert:

Da mehr Wertepaare als Unbekannte vorliegen, muss die n × m-Matrix in eine quadratische Matrix mit A^TA umgewandelt werden. Die Anzahl der Zeilen des Inhomogenitätsvektors y wird mit A^Ty der quadratischen Matrix angepasst. Es gilt also:

Beispiel: Für das Polynom 2. Grades

sollen die Parameter a, b und c mit dem erweiterten Vektoransatz berechnet werden. Die Daten für den Inhomogenitätsvektor y werden aus Tabelle 13.2 entnommen.

Die Koeffizientenmatrix A beschreibt den funktionalen Zusammenhang:

In der 1. Spalte stehen die Quadrate der x-Werte, in der 2. Spalte stehen die x-Werte aus der Wertetabelle, und in der 3. Spalte steht der Einheitsvektor e.

Für die transponierte Koeffizientenmatrix gilt:

Die Multiplikation der transponierten Koeffizientenmatrix A^T mit der Koeffizientenmatrix A ergibt die Koeffizientenmatrix der Normalgleichungen:

Der Inhomogenitätsvektor wird wie folgt berechnet:

Die Koeffizientenmatrix der Normalgleichungen muss noch invertiert werden:

Den Lösungsvektor λ erhält man durch Multiplikation der Inversen der Koeffizientenmatrix mit dem Inhomogenitätsvektor:

Die einzelnen Lösungsschritte wurden zur besseren Nachvollziehbarkeit hier ausführlich dargestellt. NumPy berechnet den Lösungsvektor des linearen Gleichungssystems in einer einzigen Programmzeile mit der Anweisung L=inv(A.T@A)@A.T@y. Die einzelnen Operationen werden von links nach rechts ausgeführt: Als Erstes wird das Matrizenprodukt A.T@A invertiert, dann wird dieses Produkt mit der transponierten Koeffizientenmatrix A.T multipliziert und zum Abschluss noch mit dem Inhomogenitätsvektor y multipliziert. Es werden also insgesamt drei Operationen benötigt: die Invertierung mit der Funktion inv(), die Transponierung mit A.T und die Matrizenmultiplikation mit dem Operator @. Die Reihenfolge der Multiplikationen darf nicht vertauscht werden.

Listing 13.4 gibt alle sechs Zwischenergebnisse einzeln aus, damit Sie den Lösungsweg besser nachvollziehen können. Der Lösungsvektor wird in die Objektvariable L gespeichert.

01  #04_vektor_parabel.py
02  import numpy as np
03  from numpy.linalg import inv
04  x=np.array([1,2,3,4,5])
05  y=np.array([5.1,2.1,1,1.8,5.2])
06  n=len(x)
07  A=np.array([x**2,x,np.ones(n)]).T
08  L=inv(A.T@A)@A.T@y #Lösungsvektor
09  print("1. Koeffizientenmatrix des überbestimmten Gleichungssystems\n",A)
10  print("2. transponierte Matrix\n",A.T)
11  print("3. Koeffizientenmatrix der Normalgleichungen\n",A.T@A)
12  print("4. normalisierter Inhomogenitätsvektor\n",A.T@y)
13  print("5. invertierte Koeffizientenmatrix der Normalgleichung\n", inv(A.T@A))
14  print("6. Lösungsvektor\n",L)
15  a,b,c=L[0],L[1],L[2]
16  print("   a  \t   b\t c")
17  print("%2.3fx^2 %2.3fx %2.3f" %(a,b,c))

Listing 13.4     Ausgleichsproblem für Polynom 2. Grades, mit dem Vektoransatz gelöst

Ausgabe

1. Koeffizientenmatrix des überbestimmten Gleichungssystems
 [[ 1.  1.  1.]
 [ 4.  2.  1.]
 [ 9.  3.  1.]
 [16.  4.  1.]
 [25.  5.  1.]]
2. transponierte Koeffizientenmatrix
 [[ 1.  4.  9. 16. 25.]
 [ 1.  2.  3.  4.  5.]
 [ 1.  1.  1.  1.  1.]]
3. Koeffizientenmatrix der Normalgleichungen
 [[979. 225.  55.]
 [225.  55.  15.]
 [ 55.  15.   5.]]
4. normalisierter Inhomogenitätsvektor
 [181.3  45.5  15.2]
5. invertierte Koeffizientenmatrix der Normalgleichung
 [[ 0.07142857 -0.42857143  0.5       ]
 [-0.42857143  2.67142857 -3.3       ]
 [ 0.5        -3.3         4.6       ]]
6. Lösungsvektor
 [ 1.05 -6.31 10.42]
   a       b   c
1.050x^2 -6.310x 10.420

Analyse

Das Programm hat prinzipiell den gleichen Aufbau wie Listing 13.3.

Die berechneten Parameter a, b und c stimmen mit den Ergebnissen aus Listing 13.2 überein.

13.1.4    Anwendungsbeispiel: Bremsweg

Schauen Sie sich ein weiteres Beispiel an, das den Prognoseaspekt und die praktische Relevanz der Ausgleichsrechnung betont und das Verständnis für das Verfahren des Vektoransatzes noch weiter festigt.

Für den Verlauf eines Bremsweges sind die Messdaten aus Tabelle 13.3 ermittelt worden:

vi in km/h	10	20	30	40	50
si in m	1,2	3,8	9,2	17	24,9

Tabelle 13.3     Bremswege in Abhängigkeit von den Geschwindigkeiten

Die Daten stammen aus [Meister: 218]. Es wurden absichtlich die gleichen Werte übernommen, damit eine Vergleichbarkeit der Ergebnisse gewährleistet ist.

Es wird angenommen, dass der Bremsweg quadratisch mit der Geschwindigkeit steigt. Ein Polynom 2. Grades

kann also als Ansatzfunktion gewählt werden. Es sind die Parameter a und b zu bestimmen.

Listing 13.5 berechnet die Parameter a und b sowie die zu prognostizierenden Bremswege s für die Geschwindigkeiten 80 km/h, 100 km/h und 150 km/h (siehe Tabelle 13.4).

v in km/h	80	100	150
s in m	63,9	99,44	222,56

Tabelle 13.4     Prognosewerte für höhere Geschwindigkeiten

Außerdem wird der Bremsweg s(v) als Funktionsplot dargestellt.

01  #05_plot_bremsweg.py
02  import numpy as np
03  from numpy.linalg import inv
04  import matplotlib.pyplot as plt
05  v=np.array([10,20,30,40,50])
06  s=np.array([1.2,3.8,9.2,17,24.9])
07  n=len(v) #np.size(v)
08  A=np.array([v**2,v]).T
09  L=inv(A.T@A)@A.T@s
10  print("Lösungsvektor L\n",L)
11  a,b =L[0],L[1]
12  #Ausgleichsfunktion für Bremsweg
13  def sa(a,b,v):
14      return a*v**2 +b*v
15  #Prognose
16  v1,v2,v3=80,100,150
17  print("s = %3.4fv^2 + %3.4fv"%(a,b))
18  print("Bremswege für:")
19  print("v = %3.2f km/h s = %3.2f m"%(v1,sa(a,b,v1)))
20  print("v = %3.2f km/h s = %3.2f m"%(v2,sa(a,b,v2)))
21  print("v = %3.2f km/h s = %3.2f m"%(v3,sa(a,b,v3)))
22  #Funktionsplot
23  fig, ax = plt.subplots()
24  va=np.linspace(0,np.max(v),500)
25  ax.set_title(r"$s\left( v\right) = av^{2}$ + bv")
26  ax.set(xlabel="Geschwindigkeit v",ylabel="Bremsweg s")
27  ax.plot(v,s,'r+')
28  ax.plot(va,sa(a,b,va),'b')
29  plt.show()

Listing 13.5     Funktionsparameter für den Bremsweg ermitteln

Ausgabe

Lösungsvektor L
[0.00978571 0.01585714]
s = 0.0098v^2 + 0.0159v
Bremswege für:
v = 80.00 km/h s = 63.90 m
v = 100.00 km/h s = 99.44 m
v = 150.00 km/h s = 222.56 m

Abbildung 13.3     Bremsweg als Funktion der Geschwindigkeit

Analyse

Entscheidend sind die Einträge in Zeile 08: np.array([v**2, v]).T. Das erste Element repräsentiert das Quadrat der Geschwindigkeit. Das zweite Element repräsentiert die Geschwindigkeit selbst. Die Reihenfolge der Einträge legt fest, in welcher Reihenfolge die Parameter a und b berechnet werden.

Die Bremswege für die Geschwindigkeiten von 80 km/h, 100 km/h und 150 km/h werden in den Zeilen 19 bis 21 berechnet und ausgegeben. Als Faustformel sollte sich der Autofahrer merken, dass sich der Bremsweg vervierfacht, wenn sich die Geschwindigkeit verdoppelt.

13.1.5    Anwendungsbeispiel: Planetenbahn

Das zweite Anwendungsbeispiel zeigt, wie die Halbachsen einer Ellipse mit dem Vektoransatz bestimmt werden können. Es stammt aus [Richter: 122].

Der Italiener Giuseppe Piazzi (1746–1826) aus Palermo beobachtete im Jahr 1801 die Positionen des Kleinplaneten Ceres innerhalb eines Zeitraums von 40 Tagen. Danach verlor er die Position des Planeten wegen dessen Nähe zur Sonne [Meister: 215]. Carl Friedrich Gauß (1777–1855) übernahm Piazzis Daten und berechnete mit der Methode der kleinsten Quadrate die Achsen a und b der elliptischen Umlaufbahn von Ceres.

Die Daten für die x,y-Koordinaten sind in Tabelle 13.5 aufgelistet.

xi	0,2	–0,7	–0,3	1,2	2,0	–2,4
yi	1,8	1,6	–1,8	–1,6	1,1	0,5

Tabelle 13.5     Koordinatendaten des Kleinplaneten Ceres

Eine Ellipse wird durch die Gleichung

beschrieben.

Mit den Daten aus Tabelle 13.5 kann mit der Substitution c = 1/a² und d = 1/b² folgendes überbestimmtes Gleichungssystem aufgestellt werden. Allgemein gilt für das lineare Gleichungssystem in Matrizenschreibweise:

In der 1. Spalte der Koeffizientenmatrix A werden die Quadrate der x-Werte und in der 2. Spalte die Quadrate der y-Werte aus Tabelle 13.5 eingesetzt:

Dieses Gleichungssystem ist überbestimmt und muss deshalb noch mit der Matrizenmultiplikation A^TA aus der transponierten Koeffizientenmatrix mit sich selbst in eine quadratische Form gebracht werden. Das NumPy-Programm aus Listing 13.6 führt diese aufwendige Operation durch.

Für den Inhomogenitätsvektor berechnet das Programm:

Das Programm löst das lineare Gleichungssystem durch Matrizenmultiplikation:

Aus dem Kehrwert des Lösungsvektors wird die Wurzel gezogen:

Die große Halbachse hat einen Wert a = 2,513 AE, und der Wert der kleinen Halbachse beträgt b = 1,783 AE (astronomische Einheit).

Listing 13.6 berechnet die Halbachsen der elliptischen Umlaufbahn des Kleinplaneten Ceres mit der NumPy-Datenstruktur array[...]. Die Matrizenmultiplikation wird mit dem Operator @ durchgeführt.

01  #06_plot_planetenbahn.py
02  import numpy as np
03  from numpy.linalg import inv
04  import matplotlib.pyplot as plt
05  #Koordinatendaten
06  x=np.array([0.2,-0.7,-0.3,1.2,2.0,-2.4])
07  y=np.array([1.8,1.6,-1.8,-1.6,1.1,0.5])
08  n=np.size(x)
09  #Koeffizientenmatrix
10  A=np.array([x**2,y**2]).T
11  #Inhomogenitätsvektor
12  e=np.ones(n)
13  #Lösung des Gleichungssystems
14  L=inv(A.T@A)@A.T@e
15  #Halbachsen berechnen
16  a,b=1./np.sqrt(L[0]),1./np.sqrt(L[1])
17  #Ausgabe
18  print("Koeffizientenmatrix des überbestimmten Gleichungssystems\n",A)
19  print("Koeffizientenmatrix der Normalgleichung\n",A.T@A)
20  print("Inhomogenitätsvektor der Normalgleichung\n",A.T@e)
21  print("Lösungsvektor\n",L)
22  print("Halbachsen der Ellipse")
23  print("a =",a)
24  print("b =",b)
25  t=np.linspace(0,2*np.pi,500)
26  xa=a*np.sin(t)
27  ya=b*np.cos(t)
28  fig, ax = plt.subplots()
29  ax.set_title(r"$\frac{x^{2}}{a^{2}} +\frac{y^{2}}{b^{2}} =1$")
30  ax.set(xlabel="x",ylabel="y")
31  ax.plot(x,y,'ro',lw=2)
32  ax.plot(xa,ya,'b--')
33  plt.show()

Listing 13.6     Berechnung der Halbachsen der elliptischen Umlaufbahn für den Kleinplaneten Ceres

Ausgabe

Koeffizientenmatrix des überbestimmten Gleichungssystems
 [[0.04 3.24]
 [0.49 2.56]
 [0.09 3.24]
 [1.44 2.56]
 [4.   1.21]
 [5.76 0.25]]
Koeffizientenmatrix der Normalgleichung
 [[51.501 11.642]
 [11.642 35.629]]
Inhomogenitätsvektor der Normalgleichung
 [11.82 13.06]
Lösungsvektor
 [0.15834493 0.31481513]
Halbachsen der Ellipse
a = 2.5130314219781282
b = 1.7822646805095228

Abbildung 13.4     Planetenbahn des Zwergplaneten Ceres

Analyse

Die Einträge für die Koeffizientenmatrix A in Zeile 10 enthalten die Quadrate der x,y-Koordinaten.

In Zeile 12 erzeugt die NumPy-Methode ones(n) den Inhomogenitätsvektor e als Einheitsmatrix mit sechs Einträgen. In der Spalte des Vektors stehen also sechs Einsen.

In Zeile 14 wir der Lösungsvektor L für die Ersatzvariablen c und d berechnet.

In Zeile 16 erfolgt die Umrechnung in die Halbachsen a und b.