Exponentialfunktion: Beim radioaktiven Zerfall ergibt sich bei
einer Zerfallskette von Isotop 1 zu Isotop 2 (Halbwertszeit
) und
dem weiteren Zerfall von Isotop 2 (Halbwertszeit
) folgender Typ
von Gleichung für die Anzahl der Teilchensorte 2
![$\displaystyle y(t) = a_1 \left[1 - \exp(-a_2 t)\right] \exp(-a_3 t) \; .$](img698.gif) |
(12.39) |
Im File expfun1.dat
finden Sie in jeder Zeile
Messwerte für
und
. 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
, 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
als
(zuerst der Paramtervektor
,
dann der Zeitvektor
) aufgerufen wird.
Der Zusammenhang zwischen den Parametern
und den physikalischen Größen
beim Zerfall lautet:
![$\displaystyle \tau_1 = \ln 2 /a_3 \;\; , \;\; \tau_2 = \ln 2 /(a_2+a_3) \;\; , \;\; y_{01} = a_1 a_2 / a_3 \; ,$](img701.gif) |
(12.40) |
wobei
die Anzahl der Teilchen 1 zum Zeitpunkt
ist.
Außerdem soll das Maximum der Funktion
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
bei fixen Werten von
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
enthalten.
Sie können Ihr numerisches Ergebnis mit dem analytischen vergleichen (ein
bisschen Analysis schadet nicht)
![$\displaystyle 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} \; .$](img706.gif) |
(12.41) |
Auch das Integral für die Fläche unter der gesamten Kurve
![$\displaystyle \int_0^{\infty}y(t)dt = \frac{a_1 a_2}{(a_2+a_3) a_3} \; ,$](img707.gif) |
(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
. Aufpassen muss man bei der oberen Grenze, da man bei
numerischen Rechnungen nicht
einsetzen kann. Man darf daher die
obere Grenze weder zu klein noch zu groß wählen.
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
liegen müssen, wobei
eine ganze Zahl ist. Für die Maxima gibt es keine analytische Lösung.