12.11 Nichtlineares Fitten

Ziel:

Das Ziel der Übung ist das Erlernen der Kurvenanpassung mit Hilfe nichtlinearer Modelfunktionen. Außerdem soll die Auswertung solcher Daten (Nullstellen, Minima, Fläche) erlernt werden. Abzugeben sind die beiden Skriptfiles expfit.m und sinfit.m.

Voraussetzung:

Mitschrift der Vorlesungsstunde über die Verwendung von Inline-Funktionen und über Funktionen zum Minimieren, Nullstellensuchen und Integrieren. Laden Sie die Files expfun1.dat und sinfun1.dat in Ihr MATLABDirectory.

1. Exponentialfunktion: Beim radioaktiven Zerfall ergibt sich bei einer Zerfallskette von Isotop 1 zu Isotop 2 (Halbwertszeit $\tau_1$) und dem weiteren Zerfall von Isotop 2 (Halbwertszeit $\tau_2$) folgender Typ von Gleichung für die Anzahl der Teilchensorte 2

   $$y(t) = a_1 \left[1 - \exp(-a_2 t)\right] \exp(-a_3 t) \; .$$ (12.39)

   
   
   Im File expfun1.dat finden Sie in jeder Zeile Messwerte für $t$ und $y(t)$. Fitten Sie dazu die obige Modellfunktion und stellen Sie sie zusammen mit den Datenwerten dar. Verwenden Sie bei der Darstellung eine größere Endzeit, Sie sollten dabei einen schönen exponentiellen Abfall der Kurve sehen.

   Für die Definition obiger Gleichung in MATLAB verwendet man den MATLAB-Befehl inline. Damit können einfache Funktionen direkt in einer Zeile eingegeben werden, ohne dass man ein Skriptfile schreiben muss. Zur nichtlinearen Anpassung kann man die Funktion nlinfit verwenden. Im Gegensatz zum linearen Fitten, wo sich aus dem Lösen des linearen Gleichungssystems immer eine exakte Lösung ergibt, braucht man hier Anfangswerte für die Parameter $a$, die nicht zu weit von der Lösung entfernt sein sollen. In diesem Fall liegen sie in der Nähe von 0.8, 4, 0.5. Bedenken Sie auch, dass bei nlinfit die Inline-Funktion für $y$ als $y(a,t)$ (zuerst der Paramtervektor $a$, dann der Zeitvektor $t$) aufgerufen wird.

   Der Zusammenhang zwischen den Parametern $a$ und den physikalischen Größen beim Zerfall lautet:

   $$\tau_1 = \ln 2 /a_3 \;\; , \;\; \tau_2 = \ln 2 /(a_2+a_3) \;\; , \;\; y_{01} = a_1 a_2 / a_3 \; ,$$ (12.40)

   
   
   wobei $y_{01}$ die Anzahl der Teilchen 1 zum Zeitpunkt $t=0$ ist.

   Außerdem soll das Maximum der Funktion $y_{max}=y(t_{max})$ gefunden und im Plot mit einem Marker dargestellt werden. Dazu verwendet man die Funktion fminsearch, wobei man als Startwert wieder einen Zeitwert in der Nähe des Maximums angeben muss. Vom Programm wird nun aber bei der Suche nach dem Maximum jeweils $t$ bei fixen Werten von $a$ verändert. Daher muss man eine modifizierte Inline-Funktion mit vertauschten Parametern schreiben. Ein zweiter Punkt zum Aufpassen liegt darin, dass wir hier ein Maximum suchen, fminsearch aber, wie schon der Name sagt, ein Minimum finden will. Daher muss diese Inline-Funktion $-y(t,a)$ enthalten.

   Sie können Ihr numerisches Ergebnis mit dem analytischen vergleichen (ein bisschen Analysis schadet nicht)

   $$t_{max} = -\ln(a_3/(a_2+a_3))/a_2 \;\; , \;\; y_{max} = \frac{a_1 a_2 \exp(a_3 \ln(a_3/(a_2+a_3))/a_2)}{a_2+a_3} \; .$$ (12.41)

   
   

   Auch das Integral für die Fläche unter der gesamten Kurve

   $$\int_0^{\infty}y(t)dt = \frac{a_1 a_2}{(a_2+a_3) a_3} \; ,$$ (12.42)

   
   
   soll mit dem numerischen Resultat verglichen werden. Dafür kann man die Routine quadl verwenden. Der Aufruf der Funktion geht auch hier mit $y(t,a)$. Aufpassen muss man bei der oberen Grenze, da man bei numerischen Rechnungen nicht $\infty$ einsetzen kann. Man darf daher die obere Grenze weder zu klein noch zu groß wählen.

2. Sinusfunktion Fitten Sie die Daten in sinfun1.dat mit Hilfe der nichtlinearen Modelfunktion

   $$y(a,x) = a_1 x \sin^2 (a_2 x) \; ,$$ (12.43)

   
   
   und stellen Sie sie zusammen mit den Daten im Intervall $[0,4\pi]$ dar. Die Startwerte liegen im Bereich von 1.8, 1.1.

   Finden Sie die ersten vier Minima außer jenem bei Null (die Minima entsprechen hier auch den Nullstellen) und die ersten vier Maxima und markieren Sie sie im Plot (Tipp: Zwei neue Inline-Funktionen $y(x,a)$ und $-y(x,a)$).

   Berechnen Sie das Integral der Funktion zwischen 0 und dem vierten Minimum. Der Wert sollte ungefähr $60$ ergeben.

3. FREIWILLIG
   MATLAB kann nicht nur numerisch sondern auch symbolisch rechnen. Die symbolischen Rechnungen zum ersten Teil der Übung kann man sich hier anschauen und auch ausprobieren.
   % Definition der Konstanten als symbolische Groessen
   
   a1 = sym('a1','positive');
   
   a2 = sym('a2','positive');
   
   a3 = sym('a3','positive');
   
   % Definition der Zeit (hier xx) und der Grenzen
   
   xx = sym('xx','positive');
   
   syms gl gu
   
   % Funktion
   
   yy = a1*(1-exp(-a2*xx)).*exp(-a3*xx)
   
   dyy = diff(yy,xx) % Ableitung nach xx
   
   % Maximum (Nullstelle der ersten Ableitung)
   
   maxx = solve(dyy,xx) % Loese dyy=0 nach xx auf
   
   maxy = factor(subs(yy,xx,maxx)) % Setze fuer xx ein
   
   % Integral
   
   iyy = int(yy,xx,gl,gu)
   
   % Unter Grenze 0; Obere Grenze Limes gegen unendlich
   
   diyy = limit(subs(iyy,gl,0),gu,inf)
   
   Hier wird zuerst mit symbolischen Größen gearbeitet und erst durch subs wird z.B. für die unter Grenze 0 zugewiesen. Die Ergebnisse kann man mit den obigen Formeln vergleichen.

   Die symbolischen Rechnungen für den zweiten Teil der Übung finden Sie hier.

   % Definition der symbolischen Groessen
   
   a1 = sym('a1','positive');
   
   a2 = sym('a2','positive');
   
   xx = sym('xx','positive');
   
   syms gl gu
   
   % Funktion und Ableitung
   
   ee = a1*xx*sin(a2*xx)2
   
   dee = diff(ee,xx)
   
   % Unbestimmtes und bestimmtes Integral
   
   iiee = int(ee,xx)
   
   iee  = int(ee,xx,gl,gu)
   
   % Einsetzen der Werte, wobei davon ausgegangen wird, dass die
   
   % Parameter a im Vektor af und dass die Minima im Vektor xmin
   
   % stehen.
   
   fl = subs(iee,{a1,a2,gl,gu},{af(1),af(2),0,xmin(4)});
   
   Mit solve findet man hier nur das triviale Minimum bei 0, dies hilft einem also hier nicht weiter. Nach Anschauen der Funktion erkennt man aber, dass die Minima natürlich bei $x_{min}=k\pi/a_2$ liegen müssen, wobei $k$ eine ganze Zahl ist. Für die Maxima gibt es keine analytische Lösung.