12.10 Lineares Fitten

Ziel:

Es soll das Fitten von Polynomen durch eine gegebene Datenmenge erlernt werden.

Voraussetzung:

Speichern Sie die Datei polydat.dat in Ihrem MATLAB-Directory ab (Skript: uebungsdaten, Übung 10). Damit bekommen Sie auch eine Reihe von Datenfiles für Kegelschnitte, die alle mit kschnitt beginnen. Die Befehle polyval und polyfit erleichtern den Umgang mit Polynomen. Das Lösen von überbestimmten Gleichungssystemen kann mit dem $\backslash$-Operator erfolgen. Das Kapitel 10 beinhaltet eine detailierte Beschreibung der Probleme.

Skriptnamen:

fitlin.m, skegel.m, kegelfit.m

12.10.1 Fitten von Polynomen

1. Die Datei polydat.dat enthält eine Messkurve. Die erste Spalte enthält die $x$-Werte, die zweite die $y$-Werte und die dritte den Fehler $\Delta y$. Lesen Sie diese Daten ein.
2. Bestimmen Sie die besten Polynome nullter, erster, zweiter und dritter Ordnung, die durch diese Daten passen.
3. Zeichnen Sie die Datenpunkte mit Fehlerbalken (errorbar) und die gefitteten Polynome (Beschriftung, Legende, ...).

12.10.2 Fitten von Polynomen mit teilweise vorgegebenen Koeffizienten

1. Schreiben Sie ein Programm, das durch die Daten polydat.dat ein Polynom 3. Ordnung legt, bei dem jedoch zwei Koeffizienten bekannt sind:
   • Koeffizient der 3. Ordnung (also der Koeffizient vor $x^3$): $a_3=1$ (Achtung der Koeffizient zur höchsten vorkommenden Potenz von $x$ steht in der MATLAB-Polynomdarstellung immer an erster Stelle im Koeffizientenvektor)
   • Koeffizient der 1. Ordnung (also der Koeffizient vor $x$): $a_1=-3$

   Dafür müssen Sie nun selbst ein lineares Gleichungssystem erstellen und dieses mit Hilfe des $\backslash$-Operators lösen. Dies liefert Ihnen dann die gewünschten Koeffizienten. Überlegen Sie sich gut, wie die rechte Seite dieses Gleichungssystems auszusehen hat.

2. Zeichnen Sie die Datenpunkte mit Fehlerbalken und das gefittete Polynom (Beschriftung, Legende, ...).

12.10.3 Erzeugen von verrauschten Daten zu Testzwecken

1. Erzeugen Sie für ein beliebiges Polynom verrauschte Testdaten. Wählen Sie dazu einen Vektor mit $x$-Werten und berechnen den entsprechenden Vektor mit $y$-Werten. Verrauschen Sie die Daten, indem Sie einen normalverteilten Zufallsvektor mit der Varianz $\sigma$ addieren (sig*randn(size(x))).
2. Stellen Sie das Polynom und die verrauschten Daten mit Fehlerbalken ( $\Delta y=\sigma$) grafisch dar.
3. Fitten Sie ein Polynom entsprechender Ordnung (mit oder ohne teilweise vorgegebene Koeffizienten) und stellen Sie es in derselben Figur grafisch dar.

12.10.4 Fitten von Kegelschnitten

Wie schon in Übung 12.9.4 besprochen, kann man Kegelschnitte durch verallgemeinerte quadratische Formen darstellen

$$z(x,y) = s_1 x^2 + 2 s_2 x y + s_3 y^2 + s_4 x + s_5 y + s_6 = 0 \; .$$ (12.36)

Liegen nun $n$ Datenpunkte für $x$ und $y$ in den Spaltenvektren $x_d$ und $y_d$ vor, kann man einen der Koeffizienten $s_i$ mit einem vorgegebenen Wert belegen (z.B.: $s_6=1$) und den entsprechenden Term auf die rechte Seite des Gleichungssystems bringen

$$s_1 x_d^2 + 2 s_2 x_d y_d + s_3 y_d^2 + s_4 x_d + s_5 y_d = - s_6 \; ,$$ (12.37)

wobei dieses Gleichungssystem nun nicht mehr homogen ist. Ist $n=5$ bekommt man als Lösung des Gleichungssystems einen Kegelschnitt auf dem alle Datenpunkte liegen. Ist $n>5$ hat man es mit einem überbestimmten linearen Gleichungssystem zu tun, das von MATLAB im ''Least Squares'' Verfahren gelöst wird. Damit wird jener Kegelschnitt gefunden für den die Summe der Abstandsquadrate

$$S = \sum_{d=1}^{n} z(x_d,y_d)^2 \;$$ (12.38)

ein Minimum ist.

Probleme mit der Lösung des Gleichungssystems bekommt man, wenn ein ''ungeeigneter'' Term auf die rechte Seite verschoben wird. Ungeeignet sind Terme, deren Koeffizienten eigentlich Null ergeben würden, wie z.B. bei der Parabel $y=x^2$ der oben gewählte Trem $s_6$. Daher lohnt es sich, alle $6$ Terme nacheinander auf die rechte Seite zu bringen, das jeweilige Gleichungssystem zu lösen und die Summe der Fehlerquadrate zu berechnen. Damit kann man dann jene Lösung mit der geringsten Summe der Fehlerquadrate als besten Fit verwenden.

1. Erzeugen Sie jetzt ein MATLAB-Skript skegel.m, das die Namen aller zugehörigen Datenfiles einliest. Dazu gibt es die Möglichkeit diese Namen in eine MATLAB-Struktur einzulesen:

     files = dir('kschnitt*.dat');  % Struktur
     anz   = length(files);         % Anzahl
     k = 5; file = files(k).name;   % Name des 5. Files
   

   Der Benutzer soll dann eine Zahl und damit einen Filenamen auswählen, den File einlesen und die Datenvektoren xd und yd belegen.
2. Diese Daten werden der MATLAB-Funktion
   [s,err2] = kegelfit(xd,yd)
   übergeben, die dann den besten Fit für die Koeffizienten $s_i$ und die zugehörige Summe der Fehlerquadrate err2 zurückgibt.

   Folgen Sie bei der Erstellung der Funktion der oben beschriebenen Strategie alle $6$ Terme abwechselnd auf die rechte Seite zu bringen. Für die Aufteilung in Koeffizientenmatrix und Inhomogenitätsvektor eignet sich die logische Indizierung hervorragend. Damit kann man das Problem leicht in einer for-Schleife erledigen.

3. Erstellen Sie danach im Skript eine Darstellung der Fläche $z(x,y)$ für $x$- und $y$-Werte ungefähr im Bereich der Daten, zeichnen Sie ein Konturlinie bei $z=0$ (entspricht Kegelschnitt) ein und plotten die zugehörigen Datenpunkte. Diese sollten, wenn alles richtig ist, in der nahen Umgebung des Kegelschnittes liegen.