from __future__ import unicode_literals #Ez teszi lehetővé az ékezetes ábrafeliratokat!!!
#Ha nem nem működik az alábbi pár feladat akkor az következő import sorok kikommentálásával próbálkozzatok!
#import numpy #numerikus számolásokhoz használt csomag
#import matplotlib #ábrázoláshoz használt csomag
Egyszerű adatbeolvasás, adatok ábrázolása és mozgó átlag számolása¶
A Nap forró plazmája a differenciális forgás és elektromágneses effektusok összjátékának eredményeképp nagyon bonyolult viselkedést mutat. Többek közt 22 évente megváltoztatja az "észak-déli" polaritás irányát, és ezzel együtt a kiáramlásokat jellemző napfoltok száma is fluktuál. A napfoltok számát sok év óta rendszeresen gyűjtik, pl. innen letölthetőek. A letöltött monthssn.dat fájlt a "+ New" menü segítségével tölthetjük fel a szerveren lévő könyvtárunkba. A fájl első néhány sorát így nézhetjük meg:
with open("monthssn.dat") as myfile: #a myfile nevű változó "hozzá kötése" a filehoz
#Figyelem ez nem azt jelenti hogy a myfile változó tartalmazni fogja a filet
#inkább úgy kell erre gondolni mint egy címkére
for i in range(1,10): # A file első 10 sorának megjelenítése
print myfile.readline() # a readline parancs segítségével a file egy sorát egy string-be tudjuk tenni
# ezt a string-et később fel tudjuk dolgozni
A második oszlop tartalmazz a tizedes tört formájában kifejezett, években mért időt, a harmadik pedig az átlagos havi napfoltszámot. Az adatokat beolvashatjuk egy két indexes tömbbe, és a megfelelő oszlopokat ábrázolhatjuk. Figyelj a 0-val kezdődő indexelésre!
data = loadtxt("monthssn.dat") # Ez egy egyszerű lehetőség csak számokat tartalmazó adatok beolvasására.
# A következő példában összetettebb formázott fájl olvasást használunk.
x = data[:,1] # az x változó tartalmazza a file-ből beolvasott adatok második oszlopát (jelen esetben ez az idő évben)
# az első oszlop
y = data[:,2] # y pedig a harmadik oszlopot (jelen esetben a napfoltok száma)
#figsize(6.0, 4.0) # az ábrák alapértelmezett mérete
#figsize(16,8) # ez az adatsor elég hosszú így célszerűbb az ábra méretét megnövelni
plot(x,y) # Az ábra gyártás leg alapvetőbb parancsa.
Itt a matplotlib könyvtár pyplot csomagjának plot parancsát használtuk és beolvasott adatokat ábrázoltunk. Természetesen mi is generálhatunk adatokat a programmal.
(Ha nem jelenne meg az ábra akkor a notebook elején lévő import parancsokat futassuk le a comment-elt sorokból!)
Mozgó átlag¶
Sokszor előfordul hogy egy mért mennyiséget sokkal sűrűbben mintavételezünk mint azt az általunk megfigyelni kívánt jelenség megkövetelné. Jó példa erre az előző feladatban látott napfoltokat leíró adat sor is. Jól láthatóan "tüskés" az adatsor. Ilyen adatsorok kisímitására egy egyserű módszer a mozgó átlag számítása. Egy \(x_i\) \(i=1\dots N\) adatsor mozgó átlagát az alábbiak szerint definiáljuk: \[\bar{x}_i=\frac{1}{2w}\sum_{j=-w}^w x_j, \] ahol \(w\) egy egész szám az átlagolási ablak fele.
Írjunk egy egyszerű függvényt a mozgó átlag kiszámítására!
#Egyszerű mozgóátlag függvény
#input:
# adat: adatsor, átlagolni kívánt tetszőleges hosszúságú adatsor
# ablak: az átlagolási ablak, egész szám
#output:
# átlagolt len(adat)/ablak hosszú adatsor
def mozgoatlag(adat,ablak):
"Mozgó átlag: Az adat változóban tárolt adatsor mozgóátlaga, ablak számú adatpont alapján"
atlagolt=[] #inicializálás
for i in range(len(adat)/ablak): #figyelem itt két egész számót osztunk egymással!
#ez határozza meg az eredmény hosszát!
atlag=sum(adat[range(ablak*i,ablak*(1+i))])*1.0/ablak #az adott ablakban szerplő pontok átlagolása a sum parancs segítségével
atlagolt.append(atlag) #az átlag hozzáfűzése a végeredményhez
return atlagolt
Ábrázoljuk a napfoltokat ábrázoló adatsort a mozgó átlaggal együtt!
figsize(16,8) # Ez is legyen nagy ábra!
#Két görbét úgy tudunk egyszerűen ugyan azon az ábrán megjeleníteni ha a plot parancsot kétszer hívjuk meg!
#Későbbiekben ki fogunk térni arra is hogy hogy kell több ábrát egymás mellé tenni.
plot(x,y, linewidth=0.2 ) #Az adatsor, vékony vonallal (0.2-es vastagság)
plot(mozgoatlag(x,10),mozgoatlag(y,10),'-',linewidth=2.0) # a mozgó átlag vastag vonallal
# Figyeljük meg hogy a mozgó átlag x értékeit és az y értékeit is a "mozgoatlag" függvény hívásával határoztuk meg!
title('Napfoltok',fontsize=25);#Az ábra címe
# tengelyfeliratozások
xlabel('t[év]',fontsize=20); # x tengely
ylabel('Napfolotk száma',fontsize=20); # y tengely
Adat feldolgozás sorról sorra, illesztés, numerikus derivált¶
File olvasás sorról sorra, adatok egyenkénti kinyerése¶
Töltsük be a "baumgartner_data" neű filet, mely tartalmazza Felix Baumgartner ugrásának magasság idő diagrammját. Az adatfile beolvasása után nyerjük ki a pálya adatpontjait és ábrázoljuk a magasság idő/diagrammot illetve a sebesség/idő diagrammot! A sebesség/idő diagrammon jelöljük meg hogy mikor volt maximális a sebesség!
datafile = open('baumgartner_data') #file megnyitása
lines = datafile.readlines() #a file összes sorának beolvasása egy stringekből álló list-be
#a readline parancs csak egy sort olvas be,
# míg itt a readlines (TÖBBESSZÁM!!) parancs az össze sort beolvassa
datafile.close() #mivel már nem kell az adatfile be is zárjuk
lines #mi volt a file-ban ? itt használhatnánk print-et is
Ez az adatfile egy kicsit bonyolultabb mint az elő példában látott. Tartalmaz három sor fejlécet, illetve egy sor láblécet, a számok vesszővel vannak elválasztva egymástól és exponenciális formában vannak megadva. Általában mielőtt egy adatfile-ból számunkra értékes információt ki tudunk nyerni fel kell dolgoznunk azt. Ez a fenti példában azt jelenti hogy a lines listából csak azokat a sorokkal szeretnénk dolgozni amiben a számok vannak illetve hogy a számok között lévő vesszőtől meg szertnénk szabadulni! Az alábbi példában használni fogjuk a reguláris kifejezések csomagot https://docs.python.org/2/library/re.html, aminek a segítségével egyszerűen kereshetünk string változókban.
datalines=((lines[4:])[:-1]) #az adatokat tartalmazó sor kiválasztása, az 5. sortól az utolsó előtti-ig !
# lines[4:] csak az 5. sor után
# ((lines[4:])[:-1]) az 5.sor után lévő rész kivéve az utolsó sor
import re # regularis kifejezések csomag
t=zeros(len(datalines)) # inicializálás, a számokat tartalmazó tömbök mérete akkora mint a számokat tartalmazó sorok száma
h=zeros(len(datalines))
for i in range(len(datalines)):
# az első oszlop az idő a második a megtett út
t[i], h[i] = map(double, re.findall(r'[+-]?\d+.\d+...',datalines[i])) # az adott sorban találd meg azokat az alstringeket
# melyek olyan alakúak hogy:
# lehet hogy + vagy - al kezdődnek de nem feltétlenül
# egy vagy több egész szám után egy tetszőleges karakter
# (tizedes pont)
# egy vagy több egész szám után három tetszőleges karakter
# (e betű, + vagy -, exponens)
# ebből a mintából találd meg az összes-et (jelen esetben csak 2)
# és generálj az alstring-ekből valós számokat (double)
Szükségünk lesz a magasság idő numerikus deriváltjára is, ezt az alábbi egyserű függvény szolgáltatja.
# numerikus derivált függvény
def nderiv(y,x):
"Első szomszéd differenciál f"
n = len(y) # adatpontok száma
d = zeros(n,'d') # változó inicializálás
# mivel a legegyszerűbb numerikus differenciálás nem szimmetrikus a végpontokat
# kicsit másképp kezeljük mint a tömb belsejében lévő pontokat
for i in range(1,n-1):
d[i] = (y[i+1]-y[i])/(x[i+1]-x[i]) #egy általános pont deriváltja
d[0] = (y[1]-y[0])/(x[1]-x[0]) #az első pont deriváltja
d[n-1] = (y[n-1]-y[n-2])/(x[n-1]-x[n-2]) # az utolsó pont deriváltja
return d
v=nderiv(h,t); # a sebesség a megtett út idő szerinti deriváltja
[t[argmin(v)],min(v)] # mivel Baumgartner esik ezért a sebessége negatív! tehát a v(t) görbe minimumát keressük
Ábrázoljuk a magasság idő illetve a sebesség idő grafikont! Jelöljük be hogy mikor volt maximális a sebesség és hogy ebben a pillanatban milyen magasan volt Bamugartner!
fig=figure(); #mivel most két ábránk is lesz ezért szükség lesz egy mind a kettőt magában foglaló "figure" objektumra,
#amire fig-ként fogunk hivatkozni.
figsize(20,4) #Széles ábrát akarunk!
# Az első ábra
# az ax1 címkével hivatkozunk az első ábrára
ax1=fig.add_subplot(1,2,1); # az alábra helyének létrehozása
# az add_subplot(Ny,Nx,i) parancs jelentése:
# Hozz létre egy Ny sorból és Nx oszlopból álló ábra táblázatot (ha még nem tetted meg...)
# és abba tedd bele az alábbi ábrát az i-edik helyre!
# Mi két egymás mellett lévő ábrát szeretnénk amelyikből most jön az első
ax1.plot(t,h,'o-'); # a h(t) függvény ábrázolása kerek pontokkal és őket összekötő egyenessel
ax1.plot(t[argmin(v)],h[argmin(v)],'ro',markersize=10); # egyetlen pont a maximális sebesség helyén
# az argmin() függvény megondja hogy hanyadik eleme a legkisebb egy vektornak
title('h(t) Baumgartner ugrás',fontsize=20); # keret és tengely feliratok
xlabel('t[s]',fontsize=20);
ylabel('h[m]',fontsize=20);
grid(True); # rácsozás
# Az második ábra
# az ax2 címkével hivatkozunk az második ábrára
ax2=fig.add_subplot(1,2,2);
ax2.plot(t,v,'o-');
ax2.plot(t[argmin(v)],min(v),'ro',markersize=10);
# nyíl húzása egy pontba és feliratozása
# a nyilat az xytext vektorból az xy pontba húzza
# a koordináták az adathoz vannak rögzítve és nem az ábra dimenzióihoz (xycoords='data',textcoords='data')
ax2.annotate(r'$v_{max}$', xy=(t[argmin(v)],min(v)), xycoords='data',
xytext=(200, -250), textcoords='data',
arrowprops=dict(facecolor='grey', shrink=0.05), # szürke legyen a nyíl
horizontalalignment='center', verticalalignment='center',fontsize=20) # a felirat tulajdonságai
title('v(t) Baumgartner ugrás',fontsize=20); # keret és tengely feliratok
xlabel('t[s]',fontsize=20);
ylabel('v[m/s]',fontsize=20);
grid(True); # rácsozás
fig.tight_layout() # hogy ne lógjanak egymásba az ábrák
Illesztés¶
Sokszor előfordul hogy egy kísérlet eredményeit valamilyen elméleti becsléssel szeretnénk összevetni, illetve hogy az elméletben szereplő paramétereket egy kísérlethez szeretnénk igazítani. Az egyik legelterjedtebb eljárás az úgynevezett legkisebb négyzetek módszere. Tegyük fel hogy a kísérleti adataink \((x_i,y_i)\) számpárokként álnak rendelkezésünkre és van \(N\) db mintapontunk, azaz \(i=1\dots N\) ! Erre az adatsorra szeretnénk illeszteni egy \(y=f(x,a,b,c,d,\dots)\) függvényt! A feladat az hogy olyan \(a,b,c,d,\dots\) paramétereket találjunk amik a legjobban közelítik az adatpontjainkat. Ezt a feladatot az \[S(a,b,c,d,\dots)=\sum_{i=1}^N (y_i-f(x_i,a,b,c,d,\dots))^2\] függvény minimalizálásával oldhatjuk meg. Ez a legkisebb négyzetek módszer lényege. Ilyen jellegű illesztési feladatokra is van pyhon csomag, az alábbiakban ezzel fogunk megismerkedni egy példán keresztül.
Az előző adatsor eleje jó közelítéssel szabad esést ír le a vége pedig egyenes vonalú egyenletes mozgást! Határozzuk meg a nehézségi gyorsulás értékét a pálya elején szereplő pontokból, illetve határozzuk meg Felix Baumgartner légellenállását a pálya végén szereplő pontokból! Kezdjük a feladatot a felhasznált fizikai jelenségek összefoglalásával!
Szabadesés az ugrás elején
Mivel a szabad esés egyenletesen gyorsuló mozgásnak felel meg ezért a magasság idő függvényt a \[h(t)=h_0+v_0t-\frac{g}{2}t^2,\] alakban kreshetjük! A nehézségi gyorsulás tehát meghatározható a pálya elejére illesztett parabola együtthatóiból hiszen a négyzetes tag együtthatójának a kétszerese épen \(g\)!
Egyenletes sebesség az ugrás végén
Mivel nincs gyorsulás ezért ezen a ponton a közeg ellenállás és a gravitációs erő kiegyenlíti egymást! Azaz \[mg=\alpha v^2 \] amiből a közeg ellenállás \(\alpha=\frac{mg}{v^2}\).
from scipy.optimize import curve_fit # görbe illesztéshez használt csomag betöltése
#Az illesztés során használt egyszerű polinomok definiálása
def axpb(x,a,b): return a*x+b
def ax2pbxpc(x,a,b,c): return a*x**2+b*x+c
def constfun(x,c): return 0*x+c
#Fittelés
params,conv = curve_fit(ax2pbxpc,t[0:5],h[0:5]) #a h(t) függvény elejére illesztünk parabolát g számolásához
ah,bh,ch=params # a params változó 3 értéket tartalmaz!! g=2*ah!
#Figyeljük meg hogy a curv_fit rutin a bemenő függvény első argumentumát tekinti az illesztendő függvény változójának
#a többi paraméter pedig illesztési paraméter
params,conv = curve_fit(axpb,t[-6:],h[-6:]) #a h(t) függvény végére illesztünk egyenest \alpha számolásához
dh,eh=params #\alpha=m*g/dh**2=m*2*ah/dh**2
Figyeljük meg: hogy ha egy rutinnak több output változója van akkor azokat a
ou1,ou2,ou3, ... outN = fuggveny(in1,in2,...,inN)
szintaxist kell alkalmazni!
#Az ábra készítése
tp=linspace(0,300,1000);
figsize(16,6)
#a label kapcsoló opció segítségével határozhatjuk meg az adott görbe, adatso r ábraleírását
plot(t,h,'o-',label="Mérési adatok")
plot(tp,ax2pbxpc(tp,ah,bh,ch),label="g="+str(ah*2)[0:8]+' $\mathrm{m}/\mathrm{s}^2$')
plot(tp,axpb(tp,dh,eh),label=r'$\alpha/m$='+str(2*ah/dh**2)[0:8]+' $1/\mathrm{m}$')
legend(loc='upper right',fontsize=30) # az ábramagyarázat pozicionálása
ylim((0,40000)); # az y tengely határai
tick_params(axis='both', which='major', labelsize=20) # a tengelyeken szereplő számok méretezése és formázása
# a szokásos tengely és ábra cimkék
title('Baumgartner ugrása',fontsize=40)
xlabel('t[s]',fontsize=30);
ylabel('h[m]',fontsize=30);
grid(True);
