from pylab import *
%matplotlib inline
from scipy import integrate
from matplotlib.patches import ConnectionPatch
# nastavitve za izris grafov (http://matplotlib.org/1.3.1/users/customizing.html)
rc('text', usetex=True)
rc('font', size=12, family='serif', serif=['Computer Modern'])
rc('xtick', labelsize='small')
rc('ytick', labelsize='small')
rc('legend', fontsize='medium')
rc('figure', figsize=(5, 3))
rc('lines', linewidth=2.0)
rc('axes', color_cycle=['k'])
Diferencialne operacije¶
Odvod¶
Količino, ki je definirana z odvodom merjene količine \(y\), lahko določimo le približno, saj v tabeli ne moremo izpeljati limitnega procesa \(h\rightarrow0\), ki ga terja matematična definicija. Zadovoljimo se z diferenčnim približkom, v katerem je \(h\) kar korak naše tabele merjenih vrednosti: \(y_i^\prime=\frac{y_{i+1}-y_i}{x_{i+1}-x_i}\). Vrednosti odvoda \(y_i^\prime\) je v tem primeru treba pripisati vrednostim \(x\) na sredi intervala, torej smo pripravili novo tabelo parov \(y_i^\prime\), \(\tilde{x}_i=\frac{x_i+x_{i+1}}{2}\).
Kadar bi radi obdržali količini \(y\) in \(y^\prime\) v isti tabeli, torej pripisani istim \(x\), lahko uporabimo naslednji postopek: skozi tri sosednje točke \((x_{i-1},y_{i-1})\), \((x_i, y_i)\) in \((x_{i+1}, y_{i-1})\) povlečemo parabolo
\[y(x)=y_{i-1}\frac{(x-x_{i})(x-x_{i+1})}{(x_{i-1}-x_{i})(x_{i-1}-x_{i+1})}+ y_{i}\frac{(x-x_{i-1})(x-x_{i+1})}{(x_{i}-x_{i-1})(x_{i}-x_{i+1})}+ y_{i+1}\frac{(x-x_{i-1})(x-x_{i})}{(x_{i+1}-x_{i-1})(x_{i+1}-x_{i})}\]
in izračunano njen odvod pri \(x=x_i\):
\[y_i^{\prime}=y_{i-1}\frac{x_i-x_{i+1}}{(x_{i-1}-x_{i})(x_{i-1}-x_{i+1})}+ y_{i}\frac{-x_{i-1}+2x_i-x_{i+1}}{(x_{i}-x_{i-1})(x_{i}-x_{i+1})}+ y_{i+1}\frac{x_i-x_{i-1}}{(x_{i+1}-x_{i-1})(x_{i+1}-x_{i})}.\]
Pri ekvidistantni tabeli, \(x_{i+1}-x_i=h\) za vsak \(i\), se izraz poenostavi:
\[y_i^{\prime}=\frac{y_{i+1}-y_{i-1}}{2h}.\]
V prvi in zadnji točki, \(x_1\) in \(x_N\), zgoraj opisani postopek ne deluje. Tam moramo odvod parabole izračunati v eni od robnih točk. Navedimo izraza za ekvidistantno tabelo:
\[y_1^{\prime}=-\frac{3y_1-4y_2+y_3}{2h},\]
\[y_N^{\prime}=\frac{3y_N-4y_{N-1}+y_{N-2}}{2h}.\]
Odvajanje podatkov prinaša še eno neprijetnost, izgubo natančnosti. Pri količkaj gladki odvisnosti \(y(x)\) sta si sosednji vrednosti v tabeli blizu tudi po velikosti in je njuna razlika v števcu odvoda majhna. Njena natančnost (absolutna napaka) je istega reda velikosti kot napaka izmerkov – v resnici je za faktor \(\sqrt{2}\) večja – zato pa je relativna napaka lahko dosti večja. Povedano drugače, napaka se začne na določeni decimalki merjenega rezultata. Ko dve taki vrednosti odštejemo, lahko izgubimo eno ali več vodilnih decimalk, tako da ima napaka v rezultatu odštevanja večji relativni delež. Zato so grafi odvodov vedno bolj zašumljeni (zobati) kot grafi izvirnih količin. Če se zanesemo na gladkost \(y\), lahko zobatost odvoda popravimo z zglajenimi formulami, ki upoštevajo širši interval sosednjih vrednosti. Ob povedanem je očitno, da potrebujemo za drugi odvod količine zelo natančne podatke, da se po dvakratnem odštevanju rezultat ne izgubi v šumu meritve.
Integral¶
Integracija, nasprotno, je stabilna operacija, ki ne ojačuje šuma meritve. Integralsko količino \(Y=\int_{x_1}^{x_N}\mathrm{d}x\,y(x)\) lahko aproksimiramo s trapezno integracijo:
\[Y_1=0,\]
\[Y_{i+1}=Y_{i}+\frac{1}{2}\left(x_{i+1}-x_{i}\right)\left(y_i+y_{i+1}\right).\]
Primer¶
V datoteki kamen.dat so podatki o poti \(s\) (2. stolpec, v m), hitrosti \(v\) (3. stolpec, v m/s) in pospešku \(a\) (4. stolpec, v m/s\(^2\)) kamna v odvisnosti od časa \(t\) (1. stolpec, v s). Pot, hitrost in pospešek so podani na tri signifikantna mesta.
podatki=loadtxt('../podatki/kamen.dat')
t=podatki[:, 0]
s=podatki[:, 1]
v=podatki[:, 2]
a=podatki[:, 3]
print podatki
Izračunajmo hitrost in pospešek kamna z odvajanjem poti po času. Na slikah (b) in (c) vidimo, da pri tem izgubimo natančnost, kar je še posebej opazno pri pospešku.
def odvod(f, h):
f1=convolve(f, [1, 0, -1], mode='same') / (2 * h)
f1[0] = -(3 * f[0] - 4 * f[1] + f[2]) / (2 * h)
f1[-1] = (f[-3] - 4 * f[-2] + 3 * f[-1]) / (2 * h)
return f1
h=t[1] - t[0]
v1=odvod(s, h)
a1=odvod(v1, h)
Izračunajmo še hitrost in pot z integriranjem pospeška kamna. Na slikah (d) in (e) vidimo, da pri integraciji izgube natančnosti ni.
v2=integrate.cumtrapz(a, t, initial=0.0)
s2=integrate.cumtrapz(v2, t, initial=0.0)
fig, ((axa, axd), (axb, axe), (axc, axf)) = subplots(3, 2, sharex='col', sharey='row', figsize=(10, 9))
axa.plot(t, s);
axa.set_ylabel(r'Pot (m)')
axa.grid(alpha=0.5)
axa.annotate('(a)', xy=(0.03, 0.92), xycoords='axes fraction')
axb.plot(t, v, t, v1, 'ro', markersize=5.0);
axb.set_ylabel(r'Hitrost (m/s)')
axb.set_ylim([0, 11])
axb.legend((r'to\v cno', r'izra\v cunano'), loc=7, frameon=True)
axb.grid(alpha=0.5)
axb.annotate('(b)', xy=(0.03, 0.92), xycoords='axes fraction')
axc.plot(t, a, t, a1, 'ro', markersize=5.0);
axc.set_ylabel(r'Pospe\v sek (m/s$^2$)')
axc.set_xlabel(r'\v Cas (s)')
axc.set_ylim([-3, 11.5])
axc.grid(alpha=0.5)
axc.annotate('(c)', xy=(0.03, 0.92), xycoords='axes fraction')
axd.plot(t, s, t, s2, 'bo', markersize=5.0);
axd.grid(alpha=0.5)
axd.annotate('(d)', xy=(0.03, 0.92), xycoords='axes fraction')
axe.plot(t, v, t, v2, 'bo', markersize=5.0);
axe.legend((r'to\v cno', r'izra\v cunano'), loc=7, frameon=True)
axe.grid(alpha=0.5)
axe.annotate('(e)', xy=(0.03, 0.92), xycoords='axes fraction')
axf.plot(t, a);
axf.set_xlabel(r'\v Cas (s)')
axf.grid(alpha=0.5)
axf.annotate('(f)', xy=(0.03, 0.92), xycoords='axes fraction')
axb.add_artist(ConnectionPatch(xyA=(0.5, 0.1), xyB=(0.5, 0.9), zorder=100, coordsA="axes fraction", coordsB="axes fraction", axesA=axa, axesB=axb, arrowstyle="->", linewidth=5, color="red"))
axc.add_artist(ConnectionPatch(xyA=(0.5, 0.1), xyB=(0.5, 0.9), zorder=100, coordsA="axes fraction", coordsB="axes fraction", axesA=axb, axesB=axc, arrowstyle="->", linewidth=5, color="red"))
axe.add_artist(ConnectionPatch(xyA=(0.5, 0.9), xyB=(0.5, 0.1), zorder=100, coordsA="axes fraction", coordsB="axes fraction", axesA=axe, axesB=axd, arrowstyle="->", linewidth=5, color="blue"))
axf.add_artist(ConnectionPatch(xyA=(0.5, 0.9), xyB=(0.5, 0.1), zorder=100, coordsA="axes fraction", coordsB="axes fraction", axesA=axf, axesB=axe, arrowstyle="->", linewidth=5, color="blue"))
axa.annotate('odvod', xy=(0.52, 0.05), xycoords="axes fraction", color="red", fontsize="large")
axb.annotate('odvod', xy=(0.52, 0.05), xycoords="axes fraction", color="red", fontsize="large")
axe.annotate('integral', xy=(0.30, 0.9), xycoords="axes fraction", color="blue", fontsize="large")
axf.annotate('integral', xy=(0.30, 0.9), xycoords="axes fraction", color="blue", fontsize="large")
subplots_adjust(hspace=0.05, wspace=0.075)
Naloge¶
-
Nariši graf diferencialne upornosti \(\mathrm{d}U/\mathrm{d}I\) za tokovno odvisnost v datoteki Korozija.dat. To je meritev karakteristike \(I\)-\(U\) za kovinsko elektrodo v določeni korozivni raztopini. V prvem stolpcu je \(U\) (mV), v drugem pa \(I\) (A).
-
Datoteka body_acc.txt vsebuje podatke o kartezičnih komponentah vektorja pospeška premikajoče se osebe (komponente \(x\), \(y\) in \(z\) pospeška v enotah \(g=9.81\,\)m/s\(^2\) so v prvem, drugem in tretjem stolpcu datoteke). Akcelerometer v mobilnem telefonu, ki ga je oseba nosila na pasu, je podatke o pospešku zapisoval s frekvenco \(50\,\)Hz. Nariši graf položaja osebe v odvisnosti od časa.
- V datoteki prevajanje.dat je podan temperaturni profil v \(0.1\,\)m debeli steni v odvisnosti od časa. V prvem stolpcu je podana oddaljenost (v metrih) merilne točke od levega roba stene, v ostalih stolpcih pa so temperature (v stopinjah Celzija) na teh mestih ob časih od \(0\,\)s do \(180\,\)s s korakom \(1\,\)s. (a) Kako se s časom spreminjata toplotna tokova na obeh robovih stene? (b) Ob času \(t=1\,\)min primerjaj krajevno odvisnost odvoda temperature po času s krajevno odvisnostjo drugega odvoda temperature po kraju.