Es bietet sich an zuerst alle Pakete zu importieren, die zum Ausführen des Notebooks benötigt werden. Stellt man während der Arbeit am Notebook fest, dass zusätzliche Pakete benötigt werden, sollte die Liste hier ergänzt werden.
In [1]:
# matplotlib und numpy importieren
%pylab nbagg
# Arbeiten mit Daten
import pandas as pd
# Zum Fitten
import lmfit
# Fehlerrechnung
import uncertainties as uct
from uncertainties.umath import *
Die Fitroutine minimiert das $\chi^2$, welches aus den Quadraten der Residuen (=Abweichung der Modellvorhersage von den Messwerten) berechnet wird. Die folgende Funktion bestimmt für ein beliebiges Modell die dazugehörige Residuenfunktion. Es ist nicht wichtig, diesen Quellcode zu verstehen. Im Zweifel reicht es, diese Zelle einfach auszuführen, damit die Funktion definiert ist.
In [2]:
def residual(userfcn):
"""
Gibt für eine Modellfunktion die entsprechende Residuenfunktion
zurück, die für gewisse Parameter und x/y-Wertepaare die Residuen
(modell(x) - y) berechnet. Wird ein Fehler angegeben, dann werden
die Residuen noch mit dem Kehrwert des Fehlers gewichtet.
Parameters
----------
userfcn : callable
The model function for which the residuals should be calculated.
This should take two arguments params and x, where params
denotes the parameters (lmfit.Parameters.parameter instance) to
be used in the model and x denotes the values for which the model
function should be evaluated.
Returns
-------
residualfunction : callable
This function calculates the residuals . It takes the model parameters and
x/y-values as arguments and optionally the uncertainties of the y-values
as a keyword argument 'eps='.
Note
----
In Summary the userfcn should be a callable of the form
model(params, x)
which returns a ndarray with the y-values corresponding to the x-values and the
model with the given parameters. The residual function which is returned is of
the form
residual(params, x, y, eps=yerr)
and returns a ndarray with the corresponding residuals for the parameters and the
given x/y-value pairs.
"""
def out(*args, **kwargs):
p, x, y = args
assert isinstance(p, lmfit.parameter.Parameters)
if "eps" in kwargs:
eps = kwargs.pop("eps")
else:
eps = None
modell = userfcn(p.valuesdict(), x)
if eps is None:
residual = abs(modell-y)
else:
residual = abs((modell-y)/eps)
return residual
return out
Zusätzlich ändern wir einige globale Optionen zur Formatierung der Plots. Es ist sinnvoll, das auch am Anfang des Notebooks zu tun, damit diese Änderungen für alle Plots wirksam sind.
In [3]:
plt.rcParams['font.size'] = 10
plt.rcParams['axes.grid'] = True
plt.rcParams['figure.figsize'] = (4,3)
plt.rcParams['figure.dpi'] = 144
Nun zur eigentlichen Auswertung. Wir beginnen damit, die Daten von Aufgabe 1 zu importieren.
Wir importieren erst einmal die Daten. Dazu nutzen wir das Paket pandas (http://pandas.pydata.org/), welches exzellent ist, um mit Daten zu arbeiten. Wir haben es als pd importiert. Das heißt, immer, wenn wir pd schreiben, weiß der Python-Interpreter, dass eigentlich pandas gemeint ist. Wir verwenden den Befehl read_csv für den Datenimport.
In [4]:
data = pd.read_csv('dat/a1-100ohm.csv')
Jetzt sind die Daten importiert und wir schauen sie uns an.
In [5]:
# Betrachte die ersten 4 (=Standardwert) Zeilen
data.head()
Out[5]:
Das sind nur Einstellungen von dem Oszilloskop, die brauchen wir gar nicht! wir betrachten noch ein paar mehr der ersten Zeilen, um herauszufinden, wo die Daten anfangen.
In [6]:
data.head(20)
Out[6]:
Offensichtlich fangen die Daten in Zeile 16 an. Wir können also getrost die ersten 16 Zeilen beim Import überspringen (die erste Zeile der Datei wird standardmäßig für die Spaltennamen verwendet).
In [7]:
data = pd.read_csv('dat/a1-100ohm.csv', skiprows=16)
Jetzt betrachten wir erneut die ersten Zeilen der importierten Daten.
In [8]:
data.head()
Out[8]:
Wir brauchen auch nur die vierte und fünfte Spalte, also importieren wir nur diese. Achtung! Python beginnt immer bei 0 mit Zählen, also hat die vierte Spalte den Index 3 und die fünfte Spalte den Index 4. Außerdem betrachten wir direkt die ersten Zeilen der importierten Daten.
In [9]:
data = pd.read_csv('dat/a1-100ohm.csv', skiprows=16, usecols=[3, 4])
data.head()
Out[9]:
Jetzt ergänzen wir noch sinnvolle Spaltenüberschriften in unserem Import-Befehl und betrachten die Daten.
In [10]:
data = pd.read_csv('dat/a1-100ohm.csv', skiprows=16, usecols=[3, 4], names=["t", "U"])
data.head()
Out[10]:
Das sieht gut aus, damit können wir arbeiten. Jetzt plotten wir die Daten. Es ist an dieser Stelle sinnvoll, die Zeitdaten direkt nach dem Import in Nanosekunden umzurechnen.
In [11]:
data = pd.read_csv('dat/a1-100ohm.csv', skiprows=16, usecols=[3, 4], names=["t", "U"])
# Zeit von Sekunden in Nanosekunden umrechnen
data.t *= 1e9
# Daten plotten, um einen Überblick zu gewinnen
data.plot(x="t", y="U")
Out[11]:
Diesen Datensatz wollen wir Fitten. Zum Fitten benötigt man immer ein Modell, Startparameter und Messwerte, für die der Fit durchgeführt werden soll. Wir definieren alles in der nächsten Zelle.
In [12]:
# Fitmodell definieren
def aufladekurve_kondensator(params, t):
try :
p = params # needed for fitting with lmfit
except:
p = params.valuesdict() # needed later for plotting a fitresult
return p["U0"] * (1 - np.exp(-(t - p["t0"]) / p["tau"])) + p["offset"]
modell = aufladekurve_kondensator
# Startparameter festlegen
params = lmfit.Parameters()
params.add("U0", value=11.6, vary=True)
params.add("tau", 500)
params.add("t0", -50)
params.add("offset", -12)
# Messwerte für das Fitten auswählen
# Die Daten für t<0 werden für den Fit hier ignoriert, da zu sehen ist,
# dass der Anstieg der Spannung nicht durch unser Modell beschrieben wird.
# Das ist hier auch nicht schlimm, da noch genügend Messdaten vorliegen.
# Es ist auch sinnvoll, da die Rechteckspannung an den Rändern nicht perfekt
# ist.
fitdata = data[(data.t > 0)]
# Fitdaten auswählen
X, Y = fitdata.t, fitdata.U
# Fehler der y-Werte definieren
Yerr = 0.03 * fitdata.U
Jetzt erstellen wir einen $\chi^2$-Minimierer und führen den Fit durch. Definiert man das Modell immer in der Variable modell, die Parameter immer in params und die Messwerte in (X, Y, Yerr). Dann kann die nächste Zelle für jeden Fit 1zu1 übernommen werden.
In [13]:
# chi2 minimizer erstellen
minimizer = lmfit.Minimizer(userfcn=residual(modell),
params=params,
fcn_args=(X,Y),
fcn_kws=dict(eps=Yerr),
nan_policy="omit"
)
# fit auch tatsächlich durchführen
result = minimizer.minimize()
# ergebnis ausgeben
print(lmfit.fit_report(result))
Jetzt plotten wir das Resultat.
In [14]:
# plotten
fig, ax = plt.subplots()
data.plot(ax=ax, x="t", y="U", label="Messwerte")
ax.plot(data.t, modell(result.params, data.t), label="Fit", c="r")
ax.legend(loc="lower right", title='Kondensator ($R=100\,\Omega$)')
ax.set_ylim(-20,20)
ax.set_xlim(-500,None)
ax.set_xlabel('Zeit $t$ [ns]')
ax.set_ylabel('Spannung $U$ [V]')
plt.tight_layout()
plt.savefig('a1-kondensator-100ohm.pdf')
plt.show()
Jetzt fitten wir die Kurve für die Spule. Die Dokumentation wird abgekürzt.
In [15]:
# Daten importieren
data = pd.read_csv('data/a1-aufladekurve-spule-100ohm.csv', skiprows=16, usecols=[3,4], names=["t", "U"])
# Zeit von Sekunden in millisekunden umrechnen
data.t = data.t * 1e6
# Fitmodell definieren
def aufladekurve_spule(params, t):
try :
p = params # needed for fitting with lmfit
except:
p = params.valuesdict() # needed later for plotting a fitresult
return p["U0"] * (np.exp(-(t - p["t0"]) / p["tau"])) + p["offset"]
modell = aufladekurve_spule
# Startparameter festlegen
params = lmfit.Parameters()
params.add("U0", 22.5, vary='False')
params.add("tau", 100)
params.add("t0", -50)
params.add("offset", -12)
# Fitten
fitdata = data[(data.t > 10)]
X, Y = fitdata.t, fitdata.U
Yerr = 0.03 * fitdata.U
minimizer = lmfit.Minimizer(userfcn=residual(modell),
params=params,
fcn_args=(X,Y),
fcn_kws=dict(eps=Yerr),
nan_policy="omit"
)
result = minimizer.minimize()
print(lmfit.fit_report(result))
# Plot erstellen
fig, ax = plt.subplots()
data.plot(ax=ax, x="t", y="U", label="Messwerte")
ax.plot(data.t, aufladekurve_spule(result.params, data.t), label="Fit", c="r")
ax.legend(loc="upper right", title='Kondensator ($R=100\,\Omega$)')
ax.set_xlim(-50,200)
ax.set_ylim(-5,25)
ax.set_xlabel('Zeit $t$ [$\mu$s]')
ax.set_ylabel('Spannung $U$ [V]')
plt.tight_layout()
plt.savefig('a1-spule-100ohm.pdf')
plt.show()
Dazu benutzen wir das Paket uncertainties, das uns die Fehlerrechnung abnimmt.
In [16]:
# Definiere Zeitkonstanten und Widerstand mit Unsicherheiten aus den Fits
tauC= uct.ufloat(242.68,1.27)
tauL = uct.ufloat(50.67, 0.25)
R = uct.ufloat(100, 100*0.05)
# Berechne Kapazität und Induktivität sowie Resonanzfrequenz
C = tauC/R * 1e-9
L = tauL * R * 1e-6
omega = 1 / (sqrt(L*C))
nu = omega / (2*pi)
# Gebe Ergebnis aus
for i,j in zip(['C = ', 'L = ', 'omega = ', 'nu = '], [C, L, omega, nu]):
print('{}{}'.format(i,j))
Wir plotten die Messung.
In [17]:
# Daten plotten
data = pd.read_csv('dat/a2-100ohm.csv', skiprows=15, usecols=[3,4], names=["t", "U"])
data.t *= 1e6
fig, ax = plt.subplots(dpi=144)
data.plot(ax=ax, x='t', legend=False)
ax.set_xlabel('Zeit $t$ [$\mu$s]')
ax.set_ylabel('Spannung $U$ [V]')
ax.set_ylim(-.3,.3)
plt.tight_layout()
plt.savefig('a2-100ohm.pdf')
plt.show()
In [18]:
# abgelesene Werte zweier Punkte gleicher Phase nach N Perioden
t1 = uct.ufloat(12.09, .5)
t2 = uct.ufloat(194.6, 1)
N = 8
# periodendauer
T = (t2 - t1) / N
# frequenz
f = 1/T
# kreisfrequenz
omega = 2*pi*f
# ausgabe der werte
for i,j in zip(['T', 'f', 'omega'], [T,f,omega]):
print('{} = '.format(i), '{:f}'.format(j))
In [19]:
# Werte aus Graphen ablesen
U1 = uct.ufloat(0.25, 0.01)
U2 = uct.ufloat(0.05, 0.01)
N = 8
# log Dekrement und Güte berechen
delta = 1 / N * log(U1 / U2)
q = pi / delta
tau = q / omega
# Ausgabe der Werte
for i,j in zip(['delta', 'q', 'tau'], [delta, q, tau]):
print('{} = '.format(i), '{:f}'.format(j))
In [20]:
# Plot erstellen
data1 = pd.read_csv('dat/a3-a-kleinefrequenz-100ohm.csv', skiprows=16, usecols=[3,4], names=["t", "U"])
data1.t *= 1e6
data2 = pd.read_csv('dat/a3-a-resonanzfrequenz-100ohm.csv', skiprows=16, usecols=[3,4], names=["t", "U"])
data2.t *= 1e6
fig, (ax1, ax2) = plt.subplots(2)
data1.plot(ax=ax1, x='t', title='$\\omega \\ll \\omega_0$', legend=False)
data2.plot(ax=ax2, x='t', title='$\\omega = \\omega_0$', legend=False)
for ax in (ax1, ax2):
ax.set_xlabel('Zeit $t$ [$\mu$s]')
ax.set_ylabel('Spannung $U$ [V]')
ax1.set_yticks([-2,-1,0,1,2])
ax2.set_yticks([-40,-20,0,20,40])
plt.tight_layout()
plt.savefig('a3-a-100ohm.pdf')
plt.show()
In [21]:
# Werte ablesen
Uklein = uct.ufloat(1.87, 0.1)
Ures = uct.ufloat(32.5, 1.2)
# Güte bestimmen
q = Ures / Uklein
# Ausgabe
print("q = ", q)
In [22]:
# Plot
data = pd.read_csv('dat/a3-b-resonanz.csv')
data.columns = ['t', 'U']
fig, ax = plt.subplots()
data.plot(ax=ax, x='t', legend=False)
ax.set_xlabel('Zeit $t$ [$\mu$s]')
ax.set_ylabel('Spannung $U$ [V]')
plt.tight_layout()
plt.savefig('a3-b-resonanz.pdf')
plt.show()
In [23]:
# Ablesen der Werte
Umax = Ures
U1 = uct.ufloat(4, 2)
U2 = uct.ufloat(30,2)
N = 14
t1 = uct.ufloat(44e-6,1e-6)
t2 = uct.ufloat(469e-6, 1e-6)
Nperioden = 19
T = (t2 - t1) / Nperioden
omega = 2 * pi / T
# Berechnung der Zeitkonstanten
delta_anschwing = 1 / N * log((Umax - U1) / (Umax - U2))
q = pi / delta_anschwing
tau = pi / (delta_anschwing * omega)
# Ausgabe
for i,j, in zip(['delta_anschwing', 'q', 'T', 'omega', 'tau'], [delta_anschwing, q, T, omega, tau]):
print('{} = {}'.format(i,j))
In [24]:
data1 = pd.read_csv('dat/a3-b-schwebung+5kHz-100ohm.csv', skiprows=16, usecols=[3,4], names=['t', 'U'])
data2 = pd.read_csv('dat/a3-b-schwebung-5kHz-100ohm.csv', skiprows=16, usecols=[3,4], names=['t', 'U'])
data3 = pd.read_csv('dat/a3-b-schwebung+10kHz-100ohm.csv', skiprows=16, usecols=[3,4], names=['t', 'U'])
data4 = pd.read_csv('dat/a3-b-schwebung-10kHz-100ohm.csv', skiprows=16, usecols=[3,4], names=['t', 'U'])
for data in [data1, data2, data3, data4]:
data.t *= 1e6
fig, ((ax1, ax3), (ax2, ax4)) = plt.subplots(2,2, figsize=(8,6), dpi=100)
for data, ax, title in zip(
[data1, data2, data3, data4],
[ax1, ax2, ax3, ax4],
['$+5$ kHz', '$-5$ kHz', '$+10$ kHz', '$-10$ kHz']
):
data.plot(ax=ax, x='t', legend=False, title=title)
for ax in (ax1, ax2, ax3, ax4):
ax.set_xlim(0,400)
ax.set_ylim(-12, 12)
for ax in (ax1, ax2):
ax.set_ylabel('Spannung $U$ [V]')
for ax in (ax2, ax4):
ax.set_xlabel('Zeit $t$ [$\mu$s]')
for ax in (ax1, ax3):
ax.set_xlabel('')
plt.tight_layout()
plt.savefig('a3-b-schwebungen.pdf')
plt.show()
Hier wurde die Datei 'a4-data-100ohm' selbst erstellt, deshalb sieht der Import-Befehl anders aus.
In [25]:
# Daten einlesen
data = pd.read_csv('dat/a4-100ohm.csv', comment='#')
# div's in Spannungen umrechnen
data.U0 = data.U0 * data.skal
data.U1 = data.U1 * data.skal
# alpha berechnen
data['alpha'] = arcsin(np.sqrt(data.U1 / data.U0)) / pi
for i in arange(14,28):
data.ix[i,'alpha'] = 1 - data.ix[i,'alpha']
data
# fehler für alpha und für f berechnen
deltaU = 0.2 * data.skal
data['alpha_err'] = np.sqrt(
deltaU**2 / data.U0**2 / (1 - (data.U1**2 / data.U0**2)) \
+ deltaU**2 * data.U1**2 / data.U0**4 / (1 - (data.U1**2 / data.U0**2))
)
data['f_err'] = 0.1
# Plotten
fig, ax = plt.subplots()
data.plot(ax=ax, kind='scatter', x='f',y='alpha', yerr='alpha_err', legend=False, edgecolor='none', s=10)
ax.set_xlabel('Frequenz $f$ [kHz]')
ax.set_ylabel('Phasenverschiebung $\\alpha$')
ax.set_yticks([0,0.25,0.5,0.75,1])
ax.set_yticklabels(['$0$', '$\\frac{\\pi}{4}$', '$\\frac{\\pi}{2}$', '$\\frac{3\\pi}{4}$', '$\\pi$'])
ax.set_ylim(0,1)
plt.tight_layout()
plt.savefig('a4-100ohm.pdf')
plt.show()
In [26]:
data.U0.max()/sqrt(2)
fres = uct.ufloat(44.70, 0.01)
# fmin und fmax aus graph ablesen
fmin = uct.ufloat(42.70, 0.01)
fmax = uct.ufloat(46.65, 0.01)
q = fres / (fmax - fmin)
tau = q / fres
for i, j in zip(['fres', 'fmin', 'fmax', 'q', 'tau'], [fres, fmin, fmax, q, tau]):
print('{} = {}'.format(i,j))
In [27]:
data = pd.read_csv('dat/a5-100ohm.csv', skiprows=16, usecols=[3,4,10], names=['t', 'CH1', 'CH2'])
# Bestimmung der Frequenzskalierung
# ---------------------------------
# Diese Frequenzen habt ihr (theoretisch) mit einer Rechteckspannung am Wobblegenerator bestimmt
fmin = uct.ufloat(10, 2)
fmax = uct.ufloat(83, 2)
# Den Wert hättet ihr (theoretisch) beim Anlegen der Rechteckspannung ablesen sollen.
# Er war aber bei jeder Praktikumsgruppe identisch mit dem der Dreieckspannung
Urechteck = uct.ufloat(40, 0.5)
# Hier muss natürlich die zweifache Amplitude der Dreieckspannung abgelesen werden,
# da die Frequenzen ja auch vom Minimum bis zum Maximum der Dreieckspannung durch-
# laufen werden
Udreieck = uct.ufloat(40, 0.5)
# Der Umrechnungsfaktor zwischen Rechteck und Dreieckspannung, da die Grenzfrequenzen
# (eigentlich) bei Rechteckspannung bestimmt werden, aber die Messung bei Dreieckspannung
# durchgeführt wird
gamma = Urechteck / Udreieck
# Die Resonanzspannung von CH1 ist eigentlich nur wichtig, um die Höhe bei dem
# 1/sqrt(2)-fachen einzuzeichnen
Ures = uct.ufloat(12.33, 0.5)
# Die Anzahl der Kästchen auf der x-Achse
smin = uct.ufloat(0.069, 0.001)
smax = uct.ufloat(0.717, 0.001)
sres = uct.ufloat(0.393, 0.001)
s = smax - smin
# Skalierung berechnen
skal = (fmax * gamma - fmin) / s
# Bestimmung von delta omega
Ured = Ures / sqrt(2)
# sklein und sgross können natürlich erst mithilfe des unteren Plots abgelesen werden, indem
# die axhline bei y=Ures/sqrt(2) hinzugefügt wird
sklein = uct.ufloat(0.382, 0.001)
sgross = uct.ufloat(0.405, 0.001)
# Güte bestimmen
delta_omega = (sgross - sklein) * skal
omega = fmin + (sres - smin) * skal
q = omega / delta_omega
for i,j in zip(['gamma', 's', 'skal', 'delta_omega', 'omega', 'q'], [gamma, s, skal, delta_omega, omega, q]):
print('{} = {}'.format(i,j))
In [28]:
# Plot
fig, ax = plt.subplots()
data.plot(ax=ax, x='t', y="CH2", color="g", legend=False)
ax2 = ax.twinx()
data.plot(ax=ax2, x='t', y="CH1", color="b", legend=False)
ax2.set_ylabel("Amplitude CH1 [div]")
ax.set_xlabel('Frequenzachse [div]')
ax.set_ylabel('Amplitude CH2 [div]')
ax.set_xlim(0.,0.8)
ax.set_ylim(-21,21)
ax.axhline(y=Ures.nominal_value / sqrt(2), xmin=0.2, xmax=0.8, c='r', linewidth=1.)
ax.axvline(x=smin.nominal_value, c='g', linewidth=1., alpha=0.5)
ax.axvline(x=smax.nominal_value, c='g', linewidth=1., alpha=0.5)
ax.axvline(x=sklein.nominal_value, c='y', linewidth=1., alpha=0.5)
ax.axvline(x=sgross.nominal_value, c='y', linewidth=1., alpha=0.5)
ax.axvline(x=sres.nominal_value, c='m', linewidth=1., alpha=0.5)
ax.set_title('Resonanzkurve, 100$\Omega$')
ax.annotate('$\\frac{U_{\\rm{res}}}{\\sqrt{2}}$',
xy=(0.655, Ures.nominal_value / sqrt(2)),
xycoords='data',
ha='left', va='center',
color="r")
handles, labels = ax.get_legend_handles_labels()
hndls2, lbls2 = ax2.get_legend_handles_labels()
handles.append(hndls2[0])
labels.append(lbls2[0])
ax.legend(reversed(handles), reversed(labels), loc="upper left", fontsize="small")
plt.grid()
plt.tight_layout()
plt.savefig('a5-100ohm.pdf')
plt.show()
Vergrößerter Plot um die Resonanz
In [29]:
# Plot
fig, ax = plt.subplots()
data.plot(ax=ax, x='t', y="CH2", color="g", legend=False)
ax2 = ax.twinx()
data.plot(ax=ax2, x='t', y="CH1", color="b", legend=False)
ax2.set_ylabel("Amplitude CH1 [div]")
ax.set_xlabel('Frequenzachse [div]')
ax.set_ylabel('Amplitude CH2 [div]')
ax.set_xlim(0.3,0.5)
ax.set_ylim(-21,21)
ax.axhline(y=Ures.nominal_value / sqrt(2), xmin=0.2, xmax=0.8, c='r', linewidth=1.)
ax.axvline(x=sklein.nominal_value, c='y', linewidth=1., alpha=0.5)
ax.axvline(x=sgross.nominal_value, c='y', linewidth=1., alpha=0.5)
ax.axvline(x=sres.nominal_value, c='m', linewidth=1., alpha=0.5)
ax.set_title('Resonanzkurve (vergrößert), 100$\Omega$')
ax.annotate('$\\frac{U_{\\rm{res}}}{\\sqrt{2}}$',
xy=(0.655, Ures.nominal_value / sqrt(2)),
xycoords='data',
ha='left', va='center',
color="r")
handles, labels = ax.get_legend_handles_labels()
hndls2, lbls2 = ax2.get_legend_handles_labels()
handles.append(hndls2[0])
labels.append(lbls2[0])
ax.legend(reversed(handles), reversed(labels), loc="upper left", fontsize="small")
plt.grid()
plt.tight_layout()
plt.savefig('a5-100ohm.pdf')
plt.show()
In [ ]: