1
0
Fork 0
mirror of https://gitlab.rlp.net/pgp/pgp1-python-einfuehrung synced 2024-11-23 13:51:35 +00:00
pgp1-python-einfuehrung/Erweiterete_Musterloesung_Fitten_mit_der_Schiefenebene.ipynb

1377 lines
47 KiB
Text
Raw Normal View History

{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Musterlösung:\n",
"\n",
"Die Musterlösugen der einzelnen Aufgaben befindet sich direkt in den entsprechenden Jupyter-Notebooks. In diesem Notebook möchte ich nochmal auf einige beliebte Fehler so wie die Auswertungsaufgabe zum Thema Fitten eingehen. \n",
"\n",
"## Fehlermeldungen in Python:\n",
"Einigen von euch sind über ein paar Fehlermeldungen gestolpert und wussten nicht genau wie diese zu beheben sind. Im Grunde sind Fehlermeldungen in Python sehr einfach zu verstehen und es bedarf lediglich ein wenig Übung. Schauen wir uns zunächst einmal eine typische Fehlermeldung in Python an:\n",
"\n",
"<img src=\"images/ExampleTraceback.png\" alt=\"Tab-Taste\" width=\"1000\"/>\n",
"\n",
"Diese sieht erst einmal schrecklich kompliziert aus, die wichtigsten Informationen sind jedoch farbig hervorgehoben. Gehen wir diese doch einmal Schritt für Schritt durch.\n",
"\n",
"1. Art der Fehlermeldung. Meistens ist dies ein erster guter Indikator wie der Fehler zustanden gekommen ist.\n",
"2. Beschreibung des Fehlers, hier steht meistens etwas ausführlicher was genau das Problem ist. Wir werde uns im Nachfolgenden noch ein paar Beispiele angucken.\n",
"3. Ort in eurem Hauptprogramm an dem der Fehler aufgetreten ist. \n",
"4. Exakte Position an welcher der Fehler aufgetreten ist. Dies kann gleich sein mit Punkt 3. sofern eure Funktion nicht innerhalb einer weiteren Funktion (hier `curve_fit`) verwendet wird. \n",
"\n",
"Der schwarz umrandete Teil kann bei sehr komplexen Funktionen sehr länglich ausfallen. In der Regel können wir dies jedoch ignorieren und uns nur auf den Anfang und das Ende der Fehlermeldung konzentrieren.\n",
"\n",
"### Beispiele:\n",
"\n",
"Im Folgenden möchte euch ein paar beliebte Fehler vorstellen. Führt einfach die entsprechenden Zellen aus um die Fehlermeldung euch anzugucken.\n",
"\n",
"Fangen wir einfach an:\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"nicht_definiert"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Diese Fehlermeldung ist relative selbsterklärend. Der Fehlertyp ist ein `NameError` und der dazugehörige Fehlertext sagt uns, dass die Variable `nicht_definiert` noch nicht definiert wurde. \n",
"\n",
"Außerdem verrät uns Python, dass `nicht_definiert` in Zeile 1 unserer Zelle ausgeführt wird. Um die Zellennummerierung angezeigt zubekommen müsst ihr die Zelle anwählen und mit Esc in den Editiermodus wechseln. Nun könnt ihr mit dem Buchstaben \"L\" die Zeilennummerierung aktivieren. Dies kann bei längeren Code-Blöcken sehr hilfreich sein. \n",
"\n",
"Ein ähnlicher Fehler tritt auf sofern ihr versucht eine Funktion aus einem Package zu laden welche nicht existiert. Zum Beispiel, weil ihr euch verschrieben habt:\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"plt.plott([1,2,3,4]) # vertippt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Da es sich hierbei jedoch um eine Funktion innerhalb eins Packages handelt spricht man von einem `AttributeError` (Warum es sich um ein Attribut handelt lernt ihr noch ausführlicher in anderen Programmierveranstaltungen). \n",
"\n",
"Richtig wäre:\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.plot([1,2,3,4]) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Ein weiterer beliebter Fehler ist der folgende:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = 1+1\n",
"b = 2\n",
"[[1,a,b,4] # Beispiel verschachtelte Liste"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Bei diesem Fehler handelt es sich um einen Syntaxfehler, der besagt, dass Python das Ende eures Codes erreicht hat (end of file EOF), jedoch nicht alle Code-Blöcke abgeschlossen sind. In diesem konkreten Beispiel liegt es daran, dass wir vergessen haben die Klammer der zweiten Liste zu schließen. Korrigiert sieht die Zelle wie folgt aus:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = 1+1\n",
"b = 2\n",
"[[1,a,b,4]] # Beispiel verschachtelte Liste"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Dieser Fehler kann auch vorkommen sofern ihr Code-Blöcke in for-Schleifen vergesst."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for i in range(4):"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Manche von euch haben versehentlich einen Zeilenumbruch innerhalb eines string verwendet. Hier bekommt ihr eine ähnliche Fehlermeldung."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(\"String mit \n",
" Zeilenumbruch\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So wäre es richtig:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(\"String mit \" \n",
" \"Zeilenumbruch\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Ein anderer Fehler, welcher mal leicht durch das überdefinieren von Variablen passieren kann, ist der folgende:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def quadrat(x):\n",
" return x**2\n",
"\n",
"quadrat = [1**2, 2**2, 3**2]\n",
"\n",
"quadrat(3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"oder"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"quadrat = 2\n",
"quadrat(3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Bei dieser Art von Fehler handelt es sich um einen TypeError. Ein TypeError besagt, dass ein Objekt (hier `quadrat`) nicht mit dem vom Python erwartenden Typ übereinstimmt. In unserem obigen Beispiel erwarten wir ein Objekt vom Typ `callable` was im Allgemeinen eine Funktion represntiert z.B.: unsere Funktion:\n",
"\n",
"```python\n",
"def quadrat(x):\n",
" return x**2\n",
"```\n",
"\n",
"Leider haben wir diese jedoch versehentlich einmal mit einer Liste und einmal mit einem Intiger überdefiniert. Python weiß hierdurch leider nicht, was `[1**2, 2**2, 3**2](3)` bzw. `3(3)` bedeuten soll. Versucht daher eindeutige Variablennamen zu benutzen. \n",
"\n",
"Als letztes hier noch ein Fehler, welcher euch beim Plotten und Fitten unterlaufen kann:\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"xWerte = [1, 2, 3, 4]\n",
"yWerte = [1, 2, 4]\n",
"\n",
"plt.plot(xWerte, yWerte)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Dieser Fehler besagt euch, dass eure Liste für die x-Werte und eure Liste für die y-Werte nicht die gleiche länge haben. Dies kann leicht beim Abtippen von Messdaten passieren.\n",
"\n",
"Ich hoffe, dass diese Beispiele euch helfen werden in Zukunft besser mit Fehlermeldungen umzugehen. Hier noch drei allgemeine Tipps zum Programmieren und Fehlerbeheben in Jupyter-Notebooks:\n",
"\n",
"1. Wenn ihr euren Code schreibt führt nach jeder neuen Zeile/Funktion euren Code aus und überprüft ob das Resultat mit euren Erwartungen übereinstimmt. \n",
"2. Solltet ihr eine Fehlermeldung bekommen und ihr seht nicht genau wie dieser Fehler zu beheben ist, kopiert euch Zeile für Zeile (bzw. Funktion) den Code in eine neue Zelle und führte diesen aus. Meistens findet man hierdurch sehr schnell wo das Problem liegt. \n",
"3. Da ihr in einem Notebook Zellen in willkürlicher Reihenfolge ausführen könnt, kann es leicht passieren, dass ihr Variablen überdefiniert. Solltet ihr also mal komplett verloren sein kann es hilfreich sein den Kernal neu zu starten, um alle Variablen zu löschen.\n",
" * Geht hierzu in der Menüleiste auf Kernal.\n",
" * Klickt anschließend auf Restart.\n",
" * Anschließend müsst ihr erneut alle Zellen (in der richtigen Reihenfolge) ausführen.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Lösungen zum Thema Fitten:\n",
"\n",
"In diesem Abschnitt wollen wir uns den Lösungen zum Thema Fitten widmen. Neben den Musterlösungen zu den Aufgaben, befindet sich am Ende nochmal ein weiteres Beispiel zu illustration.\n",
"\n",
"## Bestimmen der Fallbeschleunigung des Planeten X:\n",
"\n",
"\n",
"**Versuchsbeschreibung:**\n",
"Stellen Sie sich den folgenden Versuch vor: Jahr 2132, die Firma SpaceYpsilon hat Sie auf eine Außenmission auf den Planeten ?? geschickt. Hier sollen Sie zusammen mit ihrem Versuchspartner/in die Fallbeschleunigung g?? des Planeten bestimmen. Als Versuch lassen sie eine Kugel aus unterschiedlichen Fallhöhen innerhalb einer evakuierten Glasröhre fallen. Sie lassen die Kugel insgesamt aus 10 unterschiedlichen\n",
"Höhen Fallen.\n",
"\n",
"Basierend auf der Versuchsbeschreibung wissen wir, dass es sich bei dem Versuch um einen Freien Fall handelt, welcher als eine gleichförmig beschleunigte Bewegung beschrieben werden kann. D.h. es liegt der folgende Zusammenhang zwischen den gemessenen Höhen und Fallzeiten vor:\n",
"\n",
"$$h(t, g) = 1/2 \\cdot g \\cdot t^2$$ \n",
"\n",
"### Aufgabenstellung:\n",
"\n",
"Bestimmen Sie mit Hilfe ihrer Vorbereitungsaufgabe 1 und der entsprechenden Funktion die\n",
"Fallbeschleunigung g?? mittels eines Chi²-Fits. Diskutieren Sie anschließend mittels der Güte\n",
"Ihres Fits ob ihre Fitfunktion die gemessenen Daten gut widerspiegelt. Auf welchen Planeten\n",
"in unserem Sonnensystem befinden Sie sich?\n",
"Testen Sie anschließend ob nicht ein linearere Fit besser geeignet wäre. Begründen Sie Ihre\n",
"Antwort."
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:33:22.218156Z",
"start_time": "2020-08-26T06:33:19.532136Z"
}
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"from scipy.optimize import curve_fit\n",
"height = [1, 1.2, 1.4, 1.6, 2, 2.2, 2.4, 2.6, 2.8] # in m\n",
"dheight = [0.01]*len(height) # in m\n",
"time = [0.74, 0.8, 0.87, 0.94, 1.03, 1.1, 1.15, 1.17, 1.24] # in s\n",
"dtime = [12, 11, 9, 8, 11, 12, 13, 80, 10] # in ms\n",
"\n",
"# Achtung: Zeitfehler in s umrechnen, da in ms angegeben.\n",
"dtime = [i/1000 for i in dtime]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Beim fitten von Messdaten versuchen wir immer eine eine Funktion $y$ \n",
"\n",
"$$y(x, p1, p2, p3) = ... $$\n",
"\n",
"mit den Parameter $p_i$ an unsere Messdaten anzupassen. Genau nach diesem Schema müsst ihr auch eure Funktion in Python definieren:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def fallhoehe(t, g):\n",
" return 0.5 * g * t**2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Ein paar von euch haben leider die Funktion als $h(g, t)$ oder g(t, h) definiert statt $h(t, g)$. Dies funktioniert leider nicht. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Anschließend können wir uns die Daten erstmal angucken:"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:33:22.513340Z",
"start_time": "2020-08-26T06:33:22.220084Z"
}
},
2020-08-26 06:42:16 +00:00
"outputs": [],
"source": [
"# Plotten der Messdaten:\n",
"plt.figure(dpi=100)\n",
"plt.errorbar(time, \n",
" height, \n",
" xerr=dtime, \n",
" yerr=dheight, \n",
" ls='', \n",
" marker='.')\n",
"plt.xlabel('Zeit [s]')\n",
"plt.ylabel('Fallhöhe [m]')\n",
"plt.grid()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Die Messdaten sehen bereits leicht parabelförmig aus. Als nächstes wollen wir unser Model `fallhoehe` an unsere Daten fitten:"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:33:24.562437Z",
"start_time": "2020-08-26T06:33:24.550504Z"
}
},
"outputs": [],
"source": [
"para, pcov = curve_fit(fallhoehe, \n",
" time,\n",
" height,\n",
" sigma=dheight,\n",
" absolute_sigma=True\n",
" )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Der Fit hat funktioniert gucken wir uns also doch mal die Fitgüte und den Wert für g an:"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:33:25.833057Z",
"start_time": "2020-08-26T06:33:25.823083Z"
}
},
2020-08-26 06:42:16 +00:00
"outputs": [],
"source": [
"chi = sum([(fallhoehe(t, para[0]) - h)**2/dh**2 for t, h, dh in zip(time, height, dheight)])\n",
"\n",
"print(f'''\n",
"Das die Fitgüte Chi²/ndof lautet {chi:.2f}/{len(height)-1}\n",
"und der Wert für g ist {para[0]:.2f} +/- {pcov[0,0]:.2f} m/s\n",
"''')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Wie ihr seht ist unser $\\chi^2$/ndof (wobei ndof die Anzahl an Freiheitsgrade sind) nicht schlecht. Darüber hinaus scheint der Fehler des Wertes kleiner zu sein als die Anazahl an signifikanten Stellen die uns zur Verfügung stehen. \n",
"\n",
"Als zweites sollen wir das ganze nochmal wiederholen um zu testen ob ein lineare Fit nicht besser an die Daten passt."
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:33:27.412181Z",
"start_time": "2020-08-26T06:33:27.398218Z"
}
},
2020-08-26 06:42:16 +00:00
"outputs": [],
"source": [
"def fallhoehe2(t, v, h0):\n",
" return t * v + h0\n",
"\n",
"para2, pcov2 = curve_fit(fallhoehe2, \n",
" time,\n",
" height,\n",
" sigma=dheight,\n",
" absolute_sigma=True\n",
" )\n",
"\n",
"chi = sum([(fallhoehe2(t, para2[0], para2[1]) - h)**2/dh**2 for t, h, dh in zip(time, height, dheight)])\n",
"\n",
"print(f'''\n",
"Das die Fitgüte Chi²/ndof lautet {chi:.2f}/{len(height)-1}\n",
"und der Wert für g ist {para2[0]:.2f} +/- {pcov2[0,0]:.2f} m/s\n",
"''')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Wie ihr sehen könnt ist hier das $\\chi^2$ wesentlich schlechter und somit ist die Hypothese, dass es sich um ein lineares Model handeln könnte abgelehnt. Zu letzt können wir noch unsere beiden Ergebnisse gemeinsam mit den Daten plotten:"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:33:29.443916Z",
"start_time": "2020-08-26T06:33:29.222481Z"
}
},
2020-08-26 06:42:16 +00:00
"outputs": [],
"source": [
"# Plotten der Messdaten:\n",
"plt.figure(dpi=100)\n",
"plt.errorbar(time, \n",
" height, \n",
" xerr=dtime, \n",
" yerr=dheight, \n",
" ls='', \n",
" marker='.')\n",
"times2 = [i/100 for i in range(70, 130)]\n",
"plt.plot(times2, \n",
" [fallhoehe(t, para[0]) for t in times2],\n",
" label='Gleichförmig Beschleunigte Bewegung')\n",
"plt.plot(times2, \n",
" [fallhoehe2(t, para2[0], para2[1]) for t in times2],\n",
" label='Bewegung konstanter Geschwindigkeit')\n",
"plt.legend()\n",
"plt.xlabel('Zeit [s]')\n",
"plt.ylabel('Fallhöhe [m]')\n",
"plt.grid()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Zusatz:\n",
"\n",
"Wie wir bereits am Versuchstag selbst festgestellt haben, können bei einem einfachen $\\chi^2$ Fit lediglich die Fehler des Funktionswertes berücksichtigt werden. Da in unseren obigen Messdaten der Zeitfehler dominiert wollen wir nun noch einmal angucken, was denn passiert sofern wir die beiden Achsen tauschen. D.h. dieses mal wollen wir eine Funktion t(h, g) an unsere Messdaten anpassen:\n",
"\n",
"$$t(h, g) = \\sqrt{2 \\cdot h/g}$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Gucken wir uns zunächst wieder die Messdaten an:"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:33:37.852532Z",
"start_time": "2020-08-26T06:33:37.655061Z"
}
},
2020-08-26 06:42:16 +00:00
"outputs": [],
"source": [
"# Plotten der Messdaten:\n",
"plt.figure(dpi=100)\n",
"plt.errorbar(height, \n",
" time, \n",
" xerr=dheight, \n",
" yerr=dtime, \n",
" ls='', \n",
" marker='.')\n",
"plt.xlabel('Fallhöhe [m]')\n",
"plt.ylabel('Fallzeit [s]')\n",
"plt.grid()\n",
"plt.show()\n",
"\n",
"def fallzeit(h, g):\n",
" return (2 * h/g)**0.5"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:33:39.764472Z",
"start_time": "2020-08-26T06:33:39.757481Z"
}
},
2020-08-26 06:42:16 +00:00
"outputs": [],
"source": [
"parat, pcovt = curve_fit(fallzeit, \n",
" height,\n",
" time,\n",
" sigma=dtime,\n",
" absolute_sigma=True\n",
" )\n",
"\n",
"chit = sum([(fallzeit(h, para[0]) - t)**2/dt**2 for t, h, dt in zip(time, height, dtime)])\n",
"\n",
"print(f'''\n",
"Das die Fitgüte Chi²/ndof lautet {chit:.2f}/{len(height)-1}\n",
"und der Wert für g ist {parat[0]:.2f} +/- {pcovt[0,0]:.2f} m/s\n",
"''')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Wie ihr sehen könnt sind die Werte für $g$ fast identisch mit den Werten des vorherigen Fits, jedoch sieht das $\\chi^2$ dieses mal aufgrund der größeren Fehler besser aus. Dies spricht dafür, dass die Fehler von der Fallhöhe unterschätzt worden sind."
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:33:48.115891Z",
"start_time": "2020-08-26T06:33:47.902422Z"
}
},
2020-08-26 06:42:16 +00:00
"outputs": [],
"source": [
"plt.figure(dpi=100)\n",
"plt.errorbar(height, \n",
" time, \n",
" xerr=dheight, \n",
" yerr=dtime, \n",
" ls='', \n",
" marker='.',\n",
" label='Messdaten')\n",
"x = [i/10 for i in range(10, 30)]\n",
"plt.plot(x, \n",
" [fallzeit(i, parat[0]) for i in x],\n",
" label='t(h, g)')\n",
"plt.legend()\n",
"plt.xlabel('Fallhöhe [m]')\n",
"plt.ylabel('Fallzeit [s]')\n",
"plt.grid()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Bestimmen der Erdbeschleunigung mittels curve_fit und einer schiefen Ebene:\n",
"\n",
"Der Zusammenhang der zurückgelegten Strecke $s$ einer Vollkugel auf einer schiefen Ebene ist gegeben durch:\n",
"\n",
"$$ s(t) = \\frac{5}{14} \\cdot g \\cdot \\sin(\\alpha) \\cdot t^2$$\n",
"\n",
"Während des Versuchs wurden jedoch die Startfallhöhe $h$ und die benötigte Zeit $t_\\text{i}$ gemessen (wobei $i$ der Index für eure drei Messversuche ist.). Dies bedeutet, dass wir unsere obige Formel in Abhängigkeit dieser beiden Parameter ausdrücken müssen um $g$ mittels eines Fits bestimmen zu können. Das können wir erreichen sofern wir den folgenden Zusammenhang verwenden\n",
"\n",
"$$ \\sin(\\alpha) = \\frac{h}{l} $$\n",
"\n",
"wobei $l$ die Länge unserer schiefen Ebene ist. Setzen wir dies in unsere obige Formel ein und lösen die Gleichung nach $h$ auf so erhalten wir\n",
"\n",
"$$ h = \\frac{14 \\cdot l^2}{5} \\cdot \\frac{1}{g} \\cdot \\frac{1}{t_\\text{i}^2} $$\n",
"\n",
"wobei wir hier noch verwendet haben das die maximal zurückgelegte Strecke der Kugel nach einer Zeit $t_\\text{i}$ der Gesammtlänge der schiefen Ebene entspricht. Diese Formel für $h$ können wir nun auf unterschiedliche Arten und Weisen in Abhängigkeit der Zeit setzen wobei die bereits gezeigte Variante die erste ist:\n",
"\n",
"**Variante 2. Parabel:**\n",
"$$h(x=\\frac{1}{t_\\text{i}}) = \\frac{14 \\cdot l^2}{5} \\cdot \\frac{1}{g} \\cdot x^2$$\n",
"\n",
"**Variante 3. Ursprungsgerade:**\n",
"$$h(x=\\frac{1}{t_\\text{i}^2}) = \\frac{14 \\cdot l^2}{5} \\cdot \\frac{1}{g} \\cdot x$$\n",
"\n",
"\n",
2020-08-26 06:42:16 +00:00
"Da die erste Variante von uns am wenigsten Arbeit verlangt werden wir im folgenden diese Verwenden. \n",
"\n",
"Eure Messwerte sollten so oder so ähnlich aussehen. Es wurden für verschiedene Fallhöhen jeweils dreimal die Zeit gemessen:"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:35:01.100071Z",
"start_time": "2020-08-26T06:35:01.090098Z"
}
},
"outputs": [],
"source": [
"# Eingelesene Messwerte:\n",
"h = [0.095, 0.112, 0.134, 0.148, 0.17, 0.188, 0.21, 0.235, 0.25, 0.276] # [m]\n",
"t1 = [2.65, 2.4, 2.17, 2.06, 1.91, 1.8, 1.68, 1.6, 1.52, 1.46] # [s]\n",
"t2 = [2.71, 2.36, 2.19, 2.06, 1.9, 1.78, 1.69, 1.69, 1.53, 1.44] # [s]\n",
"t3 = [2.66, 2.36, 2.19, 2.06, 1.9, 1.8, 1.68, 1.59, 1.52, 1.44] # [s]\n",
"delta_t = [0.1]*len(h) # [s]\n",
"delta_h = [0.005]*len(h) # [m]\n",
"\n",
"l = 1.507 #[m]\n",
"delta_l = 0.005 #[m] "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Nun müssen wir uns unsere Formel definieren um unsere Messdaten fitten zu können. Hierbei ist es wichtig das eine Fitfunktion $\\lambda$ von dieser Form ist $\\lambda(x, \\Theta)$ woebei $\\Theta$ unsere Parameter sind. In unserem Fall entspricht dies einer Funktion $h(t, g)$ bzw. in den anderen beiden Varianten $h(x, g)$ wobei $x=1/t$ und $x=1/t^2$ in den Varianten 2 und 3 ist. "
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:35:12.883418Z",
"start_time": "2020-08-26T06:35:12.878431Z"
}
},
"outputs": [],
"source": [
"def fallhoehe(t, g):\n",
" l = 1.507\n",
" return 14/5 *l**2 * 1/g * 1/t**2 "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Als nächstes sollten wir zu jeder Fallhöhe den Mittelwert unserer drei Zeitmessungen bilden und den entsprechenden Fehler berechnen. Hierdurch erhalten wir ein genaueres Ergebnis für unsere gesuchte Zahl von $g$.\n",
"\n",
"D.h. wir müssen uns zunächst eine Funktion definieren, welche den Mittelwert für unsere Messwerte berechnet\n",
"\n",
"$$ \\bar{x} = \\frac{1}{n} \\sum_i^n x_i $$\n",
"\n",
"sowie dessen Standardabweichung \n",
"\n",
"$$\\sigma_\\text{n-1} = \\sqrt{\\frac{1}{n-1} \\sum_\\text{i}^\\text{n} (\\bar{x} - x_\\text{i})^2}$$"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:35:17.670513Z",
"start_time": "2020-08-26T06:35:17.662534Z"
}
},
"outputs": [],
"source": [
"# Mittelwert:\n",
"def mittelwert(l1, l2, l3):\n",
" '''\n",
" Funktion welche den Mittelwert für drei Messspalten \n",
" berechnet.\n",
" '''\n",
" result = []\n",
" for i,j,k in zip(l1, l2, l3):\n",
" result.append((i + j + k)/3)\n",
" \n",
" return result\n",
"\n",
"# Standardabweichung\n",
"def standardabweichung(l1, l2, l3):\n",
" '''\n",
" Funktion welche die Standardabweichung für drei Messspalten \n",
" berechnet. \n",
" '''\n",
" mean = mittelwert(l1, l2, l3) # <-- hier rufen wir unsere Funktion des \n",
" # Mittelwertes in einer weitere Funktion\n",
" # auf.\n",
" result = []\n",
" for m, i,j,k in zip(mean, l1, l2, l3):\n",
" result.append( (1/2 *(m-i)**2 + (m-j)**2 + (m-k)**2 )**(1/2))\n",
" \n",
" return result\n",
" \n",
"t_mean = mittelwert(t1, t2, t3)\n",
"t_std = standardabweichung(t1, t2, t3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Nun können wir unsere Messwerte ersteinmal plotten:"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:35:20.098903Z",
"start_time": "2020-08-26T06:35:20.094915Z"
}
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:35:20.926337Z",
"start_time": "2020-08-26T06:35:20.456594Z"
}
},
2020-08-26 06:42:16 +00:00
"outputs": [],
"source": [
"plt.errorbar(t_mean, \n",
" h, \n",
" xerr=t_std,\n",
" yerr=delta_h,\n",
" ls='',\n",
" marker='.',\n",
" label='Rollzeit Vollkugel')\n",
"plt.grid()\n",
"plt.xlabel('Gemessene Rollzeit $t$ [s]')\n",
"plt.ylabel('Starthöhe $h$ [m]')\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Man kann anhand des Plots bereits schön den nicht linearen Zusammenhang zwischen der Starthöhe $h$ und der Rollzeit $t$ erkennen. Als nächstes wollen wir nun aus diesen Daten unsere Erdbeschleunigung $g$ bestimmen. Hierzu verwenden wir wieder unsere Funktion `curve_fit`. Im folgenden möchte ich jedoch nochmal anhand dieses Datensatzes illustrieren, was curve_fit genau macht. Hierzu gucken wir uns die nachfolgenden Plots an (ignoriert zunächst einmal den Code)."
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:35:34.756775Z",
"start_time": "2020-08-26T06:35:33.818248Z"
}
},
2020-08-26 06:42:16 +00:00
"outputs": [],
"source": [
"for g in [7.1, 14.3, 9.81, 10.5,]:\n",
" plt.errorbar(t_mean, \n",
" h, \n",
" xerr=t_std,\n",
" yerr=delta_h,\n",
" ls='',\n",
" marker='.',\n",
" label='Rollzeit Vollkugel')\n",
"\n",
" time = [i/10 for i in range(1,30)]\n",
"\n",
" plt.plot(time, \n",
" [fallhoehe(t, g) for t in time],\n",
" label=f'h(t, g={g})')\n",
"\n",
" plt.xlim(1.4, 2.7)\n",
" plt.ylim(0, 0.4)\n",
" plt.grid()\n",
" plt.xlabel('Gemessene Rollzeit $t$ [s]')\n",
" plt.ylabel('Starthöhe $h$ [m]')\n",
" plt.legend()\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Wie ihr sehen könnt, können wir durch einfaches ausprobieren feststellen welcher Wert von $g$ am besten zu unseren Messwerten passt. Und genau das macht curve_fit für euch. Curve_fit probiert solange (nach der Methode der kleinsten Quadrate), verschiedene Werte von $g$ aus bis es den am besten passenden Wert gefunden hat. Probieren wir dies nochmal aus:"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:36:34.611760Z",
"start_time": "2020-08-26T06:36:34.606773Z"
}
},
"outputs": [],
"source": [
"from scipy.optimize import curve_fit"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:36:34.901660Z",
"start_time": "2020-08-26T06:36:34.894687Z"
}
},
"outputs": [],
"source": [
"parameter, covariance_matrix = curve_fit(fallhoehe,\n",
" t_mean,\n",
" h,\n",
" sigma=delta_h,\n",
" absolute_sigma=True\n",
" )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Hierbei schreibt curve_fit das Ergebnis der besten Werte in die Variable `parameter` und den deren Fehler in eine so genannte Kovarianzmatrix. Sprich das Ergebnis sieht für eine Funktion mit drei Parametern `def f(x, p1, p2, p3):` allgemein so aus:\n",
"\n",
"```\n",
"paramter = [p1, p2, p3]\n",
"covariance = [[cov_1,1, cov_1,2, cov_1,3], \n",
" [cov_2,1, cov_2,2, cov_2,3],\n",
" [cov_3,1, cov_3,2, cov_3,3]]\n",
"```\n",
"wobei `cov_i,i` der Varianz sprich dem $\\sigma^2$ des Parameters i entspricht. Da wir in unserem Beispiel lediglich einen Parameter haben sieht das Ergebnis wie folgt aus:"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:36:47.493940Z",
"start_time": "2020-08-26T06:36:47.487956Z"
}
},
2020-08-26 06:42:16 +00:00
"outputs": [],
"source": [
"print(f'g ist {parameter[0]:.2f} m/s^2') # <-- einfache liste\n",
"print(f'Delta g ist {(covariance_matrix[0][0])**(1/2):.2f} m/s^2') # <-- doppel liste"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Nun sollten wir noch das $\\chi^2$ und die Anzahl an Freiheitsgraden berechnen um zu gucken wie gut unser Fit funktioniert hat:"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:38:27.965486Z",
"start_time": "2020-08-26T06:38:27.958504Z"
}
},
"outputs": [],
"source": [
"def chiquadrat(xwerte, ywerte, dywerte, fun, g):\n",
" chi = 0\n",
" for x,y,dy in zip(xwerte, ywerte, dywerte):\n",
2020-08-26 06:42:16 +00:00
" chi += (fun(x, g) - y)**2/dy**2 # Der Operator += addiert zu dem vorherigen Wert von chi unseren\n",
" # neuen Wert drauf und aktualisiert gleich danach chi.\n",
" # Sprich es ist die kurzschreibweise für chi = chi + ...\n",
" return chi"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:38:33.333572Z",
"start_time": "2020-08-26T06:38:33.327134Z"
}
},
2020-08-26 06:42:16 +00:00
"outputs": [],
"source": [
"chi = chiquadrat(t_mean, h, delta_h, fallhoehe, parameter[0])\n",
"ndof = len(h) - 1 # Anzahl Messwerte - Anzahl der Fitparamter\n",
"\n",
"print(f' Das chi-quadrat und die Anzhal der Freiheitsgrade sind: {chi:.0f}/{ndof}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Wie ihr seht ist das $\\chi^2/$ndof > 1, was bedeutet, dass unsere Fitfunktion die Daten nicht ganz wiederspiegelt. Um dies nachvollziehen zu können gucken wir uns erst einmal den finalen Plot an:"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:38:38.267629Z",
"start_time": "2020-08-26T06:38:37.982393Z"
}
},
2020-08-26 06:42:16 +00:00
"outputs": [],
"source": [
"plt.figure(figsize=(5.6, 3.8), dpi=150) # <-- Größe eines A4-Blatts ausnutzen\n",
"# Plot der Messdaten\n",
"plt.errorbar(t_mean, \n",
" h, \n",
" xerr=t_std,\n",
" yerr=delta_h,\n",
" ls='',\n",
" marker='.',\n",
" label='Rollzeit Vollkugel')\n",
"\n",
"# Fitergebnis:\n",
"time = [i/10 for i in range(1,30)]\n",
"plt.plot(time, \n",
" [fallhoehe(t, parameter[0]) for t in time],\n",
" label=f'Fitparamter:\\ng: ({parameter[0]:.2f}+/-{(covariance_matrix[0][0])**(1/2):.2f}) m/s$^2$\\n'\n",
" f'$\\chi^2/$ndof: {chi:.0f}/{ndof}')\n",
"\n",
"plt.xlim(1.4, 2.7)\n",
"plt.ylim(0, 0.4)\n",
"plt.grid()\n",
"plt.xlabel('Gemessene Rollzeit $t$ [s]', fontsize=10)\n",
"plt.ylabel('Starthöhe $h$ [m]', fontsize=10)\n",
"plt.legend(fontsize=10)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Wie ihr seht weicht unsere Funktion lediglich leicht von dem Großteil unserer Messdaten ab. Lediglich für kleine und große Fallhöhen ist die Abweichung stärker und unsere Funktion beschreibt die Messdaten nicht genau genug. Dies erklärt den leicht erhöhten Wert für unser $\\chi^2$/ndof. Woran könnte dies liegen? "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Bestimmen von mehren Parameter mittels curve_fit:\n",
"\n",
"Und weil es so schön ist hier noch ein letztes Beispiel zum Thema fitten. Einige von euch haben sich gefragt warum ich überhaupt fitten muss wenn ich bei einfachen Funktionen wie zum Beispiel\n",
"\n",
"$$ s(t) = 1/2 \\cdot g \\cdot t^2$$\n",
"\n",
"die Funktion einfach nach $g$ auflösen und den Mittelwert und die Standardabweichung von $g$ berechnen kann. Bei Funktionen welche lediglich nur von einem Parameter abhängen geht das relative einfach aber wie sieht es im folgenden Beispiel aus:\n",
"\n",
"$$ T(t, T_0, \\tau, t_0) = \\tau \\cdot \\cos\\bigg(2 \\cdot \\pi \\cdot \\bigg(\\frac{t-t0}{365 d}\\bigg)\\bigg) + T_0 $$\n",
"\n",
"Die Funktion $T(t, T_0, \\tau, t_0)$ soll die jährlichen Temperaturschwankungen an einem bestimmten Ort auf der Erde wiederspiegeln. Hierbei ist $T_0$ die Durchschnittstemperatur, $\\tau$ der Temperatur unterschied und $t0$ eine Verschiebung des Cosinus entsprechend des Tages an dem die maximale Temperatur innerhalb eines Jahres gemessen wurde. Gucken wir uns zunächst die Messdaten an: "
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:38:46.482971Z",
"start_time": "2020-08-26T06:38:46.478982Z"
}
},
"outputs": [],
"source": [
"import numpy as np # trigonometrische Funktionen findet ihr ebenfalls in numpy \n",
"import matplotlib.pyplot as plt\n",
"from scipy.optimize import curve_fit"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:38:46.870189Z",
"start_time": "2020-08-26T06:38:46.855230Z"
}
},
"outputs": [],
"source": [
"# Gemessene Werte:\n",
"tage = [0.28, 10.36, 20.4, 30.23, 40.22, 50.11,\n",
" 60.25, 70.22, 80.25, 90.03, 100.24, 110.21,\n",
" 120.22, 130.25, 140.14, 150.09, 160.33, 170.31,\n",
" 180.27, 190.28, 200.25, 210.33, 220.18, 230.15,\n",
" 240.19, 250.37, 260.39, 270.35, 280.56, 290.23, \n",
" 300.31, 310.17, 320.2, 330.11, 340.28, 350.48, 360.26] # d\n",
"\n",
"gemessene_temperatur = [15.17, 15.31, 14.46, 16.2, 15.49,\n",
" 16.18, 17.18, 16.17, 17.43, 18.24,\n",
" 18.96, 19.69, 20.19, 21.33, 22.27, \n",
" 23.14, 23.6, 23.37, 23.39, 25.27, \n",
" 25.2, 24.63, 23.22, 23.95, 23.53, \n",
" 22.9, 22.59, 21.84, 20.77, 20.12, \n",
" 18.79, 18.29, 17.87, 16.86, 16.48, \n",
" 15.41, 14.2] # °C\n",
"\n",
"fehler_temperatur = [0.52, 0.54, 0.54, 0.54, 0.51, 0.48, \n",
" 0.44, 0.39, 0.33, 0.28, 0.24, 0.22, \n",
" 0.24, 0.28, 0.33, 0.39, 0.44, 0.48, \n",
" 0.52, 0.54, 0.55, 0.55, 0.53, 0.5, \n",
" 0.46, 0.41, 0.35, 0.29, 0.25, 0.22, \n",
" 0.23, 0.26, 0.31, 0.37, 0.42, 0.47, \n",
" 0.5]\n",
"\n",
"def temp(t, T0, tau, t0):\n",
" \"\"\"\n",
" Jahrestemperaturverlauf für Ort x.\n",
" \n",
" Args:\n",
" t: Zeit in Tagen\n",
" T0: Mittlere Temperatur in °C\n",
" tau: Temperaturschwankungsamplitude in °C\n",
" t0: Zeitpunkt des heißesten Tages in Tagen\n",
" \"\"\"\n",
" return tau * np.cos(2*np.pi*(t - t0)/365) + T0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Zunächst können wir uns ja mal die Messdaten angucken:"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:38:48.953130Z",
"start_time": "2020-08-26T06:38:48.766630Z"
}
},
2020-08-26 06:42:16 +00:00
"outputs": [],
"source": [
"# Plot der Messdaten:\n",
"plt.figure(dpi=100)\n",
"plt.errorbar(tage, \n",
" gemessene_temperatur,\n",
" fehler_temperatur, \n",
" ls='',\n",
" marker='.')\n",
"plt.xlabel('Tage des Jahres')\n",
"plt.ylabel('Temperatur [°C]')\n",
"plt.grid()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Als nächstes führen wir, wie in den anderen Beispielen auch, den Fit mittelts curve_fit durch. "
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:38:50.796960Z",
"start_time": "2020-08-26T06:38:50.789977Z"
}
},
"outputs": [],
"source": [
"# Fitten der Messdaten\n",
"para, pcov = curve_fit(temp,\n",
" tage,\n",
" gemessene_temperatur,\n",
" sigma=fehler_temperatur,\n",
" absolute_sigma=True\n",
" )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Im vergleich zu vorher haben wir diesesmal mehre Parameter sprich para ist jetzt eine Liste von Werten und die Kovarianzmatrix pcov eine verschachtelte Liste:"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:38:55.005310Z",
"start_time": "2020-08-26T06:38:54.998325Z"
}
},
2020-08-26 06:42:16 +00:00
"outputs": [],
"source": [
"para"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:38:55.422567Z",
"start_time": "2020-08-26T06:38:55.414589Z"
}
},
2020-08-26 06:42:16 +00:00
"outputs": [],
"source": [
"pcov"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Denkt daran, dass die Fehler euer Parameter der Wurzel der Hauptdiagonalen der Kovarianzmatrix entsprechen. Gucken wir uns doch nochmal die durch den Fit berechneten Werte etwas genauer an:"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:38:58.444536Z",
"start_time": "2020-08-26T06:38:58.438539Z"
}
},
2020-08-26 06:42:16 +00:00
"outputs": [],
"source": [
"# Printausgabe der Parameter:\n",
"for ind, (pname, einheit) in enumerate(zip(('T0', 'tau', 't0'), ('°C', '°C', 'd'))):\n",
" wert = para[ind]\n",
" fehler = pcov[ind, ind]**0.5 # dies entspricht der Hauptdiagnolen mit den Indizes 1,1 2,2 etc.\n",
" print(f'Der Wert für ist {pname} ({wert:.2f} +/- {fehler:.2f}) {einheit}')\n",
" \n",
"# Zusatzinfo:\n",
"# In unserer for-Schleife haben wir einen weiteren nützlichen Befehl eingebaut:\n",
"# enumerate gibt euch einen Index enstsprechend dem aktuellen Schritt in euer Schleife \n",
"# Probiert doch mal das folgende aus:\n",
"# for ind, buchstabe in enumerate(['A', 'B', 'C'], start=0):\n",
"# print(ind, buchstabe)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Nun sollten wir auch das $\\chi^2$ berechnen um ein Gefühl für die Fitgüte zu bekommen:"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:39:02.007699Z",
"start_time": "2020-08-26T06:39:01.998689Z"
}
},
2020-08-26 06:42:16 +00:00
"outputs": [],
"source": [
"# Berechnen des Chi**2:\n",
"chi_liste = []\n",
"for t, T, dT in zip(tage, gemessene_temperatur, fehler_temperatur):\n",
" chi_liste.append((temp(t, para[0], para[1], para[2]) - T)**2/dT**2)\n",
"chi = sum(chi_liste) \n",
"print(f'Das Chi-Quadrat beträgt {chi:.2f} mit {len(gemessene_temperatur) - 3} Freiheitsgraden.')\n",
"\n",
"\n",
"# Zusatzinfo: \n",
"# Ihr könnt das ganze auch wieder etwas kompakter als list comprehension schreiben.\n",
"# Außerdem könnt ihr die Fitparameter anstatt als para[0], para[1], para[2] mit Hilfe von\n",
"# *para an die Funktion geben. Der Stern ordnet die Werte in euerer Liste der Reihe nach den \n",
"# Argumenten euer Funktion zu:\n",
"# chi = sum([(temp(t, *para) - T)**2/dT**2 for t, T, dT in zip(tage, gemessene_temperatur, fehler_temperatur)])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Wie ihr seht scheinen die Fitparameter den Funktionsverlauf ganz gut zu beschreiben, gucken wir uns also mal das Resultat zusammen mit den Messwerten an."
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:39:04.861580Z",
"start_time": "2020-08-26T06:39:04.639174Z"
}
},
2020-08-26 06:42:16 +00:00
"outputs": [],
"source": [
"plt.figure(dpi=100)\n",
"plt.errorbar(tage, \n",
" gemessene_temperatur,\n",
" fehler_temperatur, \n",
" ls='',\n",
" marker='.',\n",
" label='T an Ort x')\n",
"tage2 = [t/10 for t in range(3650)]\n",
"plt.plot(tage2, temp(tage2, *para), label='Fitfunktion')\n",
"plt.legend(loc=2)\n",
"plt.xlabel('Tage des Jahres')\n",
"plt.ylabel('Temperatur [°C]')\n",
"plt.grid()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Hier noch ein kleiner Zusatz: Ist euch etwas aufgefallen? Der Fitparameter $T_0$ hat einen Wert von ca. 20 Tagen obwohl das Maximum eher bei 200 Tagen liegt. Dennoch beschreibt der Fit den Verlauf der Messdaten sehr gut. Dies liegt daran, das der Cosinus eine periodische Funktion ist. Bei der Methode der kleinsten Quadrate werden die Fitparameter so lange, nach einem gewissen Schema variiert, bis das Chi-Qudarat minimal ist. Bei einer periodischen Funktion gibt es mehre dieser Minima. \n",
"\n",
"Dies kann auch bei anderen komplexeren Funktionen der Fall sein. Da die meisten Funktionen jedoch nicht periodisch sind handelt es sich in der Regel bei diesen zusätzlichen Minima nur um lokale Minima. Um die besten Fitparameter zu finden wollen wir jedoch das globale Minimum finden. Um dies zu erreichen können wir curve_fit ein wenig helfen und zum Beispiel noch zusätzlich Startwerte für unsere Fitparameter mitgeben: "
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:39:40.646282Z",
"start_time": "2020-08-26T06:39:40.638303Z"
}
},
"outputs": [],
"source": [
"# Die Startwerte sind entsprechend der Reihenfolge der Parameter der Fitfunktion anzugeben.\n",
"# In unserem Fall ist dies T0, tau, t0. Sollten wir keine Idee für einen Startwert haben \n",
"# können wir einfach den entsprechenden Wert auf einen passenden Wert setzen hier z.B, 1.\n",
"startwerte = [1, 1, 200]\n",
"\n",
"# Fitten der Messdaten\n",
"para, pcov = curve_fit(temp,\n",
" tage,\n",
" gemessene_temperatur,\n",
" sigma=fehler_temperatur,\n",
" absolute_sigma=True,\n",
" p0=startwerte # <-- Übergeben der Startwerte an die Funktion:\n",
" )"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:39:41.286672Z",
"start_time": "2020-08-26T06:39:41.278694Z"
}
},
2020-08-26 06:42:16 +00:00
"outputs": [],
"source": [
"# Erneutes printen der Fitwerte und des Chi**2\n",
"for ind, (pname, einheit) in enumerate(zip(('T0', 'tau', 't0'), ('°C', '°C', 'd'))):\n",
" wert = para[ind]\n",
" fehler = pcov[ind, ind]**0.5 # dies entspricht der Hauptdiagnolen mit den Indizes 1,1 2,2 etc.\n",
" print(f'Der Wert für ist {pname} ({wert:.2f} +/- {fehler:.2f}) {einheit}')\n",
" \n",
"chi = sum([(temp(t, *para) - T)**2/dT**2 for t, T, dT in zip(tage, gemessene_temperatur, fehler_temperatur)])\n",
"chi = sum(chi_liste) \n",
"print(f'Das Chi-Quadrat beträgt {chi:.2f} mit {len(gemessene_temperatur) - 3} Freiheitsgraden.')"
]
},
{
"cell_type": "code",
2020-08-26 06:42:16 +00:00
"execution_count": null,
"metadata": {
"ExecuteTime": {
2020-08-26 06:42:16 +00:00
"end_time": "2020-08-26T06:39:42.279956Z",
"start_time": "2020-08-26T06:39:42.070512Z"
}
},
2020-08-26 06:42:16 +00:00
"outputs": [],
"source": [
"# Erneutes plotten der Funktion:\n",
"plt.figure(dpi=100)\n",
"plt.errorbar(tage, \n",
" gemessene_temperatur,\n",
" fehler_temperatur, \n",
" ls='',\n",
" marker='.',\n",
" label='T an Ort x')\n",
"tage2 = [t/10 for t in range(3650)]\n",
"plt.plot(tage2, temp(tage2, *para), label='Fitfunktion')\n",
"plt.legend(loc=2)\n",
"plt.xlabel('Tage des Jahres')\n",
"plt.ylabel('Temperatur [°C]')\n",
"plt.grid()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Wie ihr sehen könnt konnten wir mit Hilfe der Startwerte den Fit so beeinflussen, dass curve_fit dieses mal das \"richtige\" Minimum finden konnte. Daher empfehle ich euch bei komplexeren Problem sich immer erst die Messdaten an zugucken und ein paar Startwerte für den Fit zu raten/schätzen. "
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
2020-08-25 14:48:41 +00:00
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}