1
0
Fork 0
mirror of https://gitlab.rlp.net/pgp/pgp1-python-einfuehrung synced 2024-10-12 13:24:22 +00:00
pgp1-python-einfuehrung/Erweiterete_Musterloesung_Fitten_mit_der_Schiefenebene.ipynb

1103 lines
38 KiB
Text
Raw Normal View History

{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 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",
"# Zeitfehler in s umrechnen:\n",
"dtime = [i/1000 for i in dtime]\n",
"\n",
"def fallhoehe(t, g):\n",
" return 0.5 * g * t**2"
]
},
{
"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",
2020-08-26 06:42:16 +00:00
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}