Moto: panda v koruně pevného stromu
Panda za nás prováděla hodně věcí automaticky -- typické pro Python. Řešení problémů je obvykle představováno kombinací několika základních kamenů. Příkladem budiž kreslení grafů.
Matplotlib.org je grafický modul (sada nástrojů pro kreslení grafů).
Zpřístupní se pomocí kouzelného slůvka
In [1]:
import matplotlib.pyplot as plt # plt je vseobecne uzivana zkratka, grafy si kreslime primo do notebookove stranky:
%matplotlib inline
Vykreslení dat do grafu zajistí
In [2]:
plt.plot(1,2) # bod v x = 1, y = 2
Out[2]:
Fakt je tam... jen není vidět.
In [3]:
plt.plot(1,2,'*') # muj oblibeny symbol je pentagram
Out[3]:
Žádáme-li víc bodů, musíme je uzavřít do hranatých závorek.
In [4]:
plt.plot([1,2,3],[11,12,13],'+') # pokud nedefinujeme kreslici symbol, spoji se carou
Out[4]:
In [5]:
plt.plot?
Zjevně bude potřeba upravit měřítka na osách:
In [6]:
plt.xlim(0.9,3.1)
plt.ylim(10,15)
plt.plot([1,2,3],[11,12,13],'*')
Out[6]:
Poznámky:
Chceme-li vykreslit mnoho bodů, je nutné je zapsat do seznamu. Seznam je něco jako pojmenovaný sloupec v tabulkovém procesoru.
Vytváří se výčtem nebo vygenerováním nebo ...:
In [7]:
x = [1,2,3] # vycet
y = [] # prazdny seznam
Vlastnosti seznamu:
In [8]:
x.append(4) # pridani prvku s hodnotou 4 nakonec
print(x)
In [9]:
x.pop() # odebrani posledniho
print(x)
In [10]:
print(x[2]) # vypsani s indexem 2 -- pozice 3
In [11]:
print(len(x)) # zjisteni delky
In [12]:
x = [1,2,3] # prvni sloupecek
y = []
y.append(1 + x[0]**2)
y.append(1 + x[1]**2)
y.append(1 + x[2]**2)
# kapanek neprehledne
In [13]:
plt.xlim(0.9,3.1)
plt.ylim(0,12)
plt.plot(x,y,'*') # presna obdoba
Out[13]:
Pythoní obdoba udělej něco pro několik prvků je
$$y_n = 1 + x_n^2, \quad \mathrm{pro} \ n = 1, N$$
In [14]:
y = [1 + x[n]**2 for n in range(0,len(x))]
print(y)
In [15]:
plt.xlim(0.9,3.1)
plt.ylim(0,12)
plt.plot(x,y,'*')
Out[15]:
Funkci $y = f(x)$ vykreslíme tak, že vygenerujeme sadu (vektor) bodů $\pmb{x} = [x_1, x_2, \ldots ,x_N]$ a vypočteme $y_n = f(x_n)$ pro $n=1,2,\ldots, N$.
Jenže jak vygenerovat vektor?
In [16]:
import numpy as np
numpy definuje operace s vektory, maticemi, umí je násobit, hledat vlastní čísla,...
Definice vektoru pomocí listu (seznamu):
In [17]:
x = np.array([1,2,3,4,5]) # vektor x
$\pmb{y} = f(\pmb{x}) = 1 + \pmb{x}^2$
In [18]:
y = 1 + x**2 # vektor y je vypocitany z vektoru x, prvek po prvku....
In [19]:
plt.plot(x,y,'*')
Out[19]:
In [20]:
x = np.zeros(5) # petiprvkovy vektor plny nul
print(x)
In [21]:
x = np.linspace(1,5,50) # vektor s prvky 1 az 5, 50 prvku
print(x)
print(x[0],x[1],x[50-1],(5-1)/49)
A právě toto je trik, který potřebujeme na vygenerování $x_n$.
In [22]:
y = 1 + x**2
In [23]:
plt.plot(x,y)
Out[23]:
(Už vidíme, proč má funkce
Násobení vektorů ve mně vyvolává pnutí:
In [24]:
x = np.array([0,1,0]) # jednotkovy vektor podel y osy
y = np.array([1,0,0]) # jednotkovy vektor podel x osy
s1 = np.inner(x,y) # inner() je skalarni soucin
s2 = np.cross(x,y) # cross() je vektorovy soucin
print("Skalární součin ",x,"*",y,"=",s1)
print("Vektorový součin ",x,"x",y,"=",s2)
Zkusíme několik speciálních funkcí (takových, které se nedají zapsat konečným počtem operací).
In [25]:
import scipy.special as spec
In [26]:
x = np.linspace(0,10)
y = np.sin(x)
z = np.exp(-x)
plt.plot(x,y)
plt.plot(x,z)
Out[26]:
In [27]:
x = np.linspace(0,30) # bude to zubate
bessel = spec.j1(x) # funkci je tolik, ze si kazdy muze zvolit svou vlastni...
plt.plot(x,bessel)
Out[27]:
In [28]:
x = np.linspace(0,30,1000) # zvolime jemnejsi deleni intervalu
bessel = spec.j1(x)
plt.plot(x,bessel)
Out[28]:
def jmeno(argumenty):
prikaz(y)
return vysledek
In [29]:
# priklad
def fun(x):
print("Argument:",x)
return x + 1
x = 3
print("Fun with fun:",fun(x))
Nyní máme vše připraveno pro kreslení složených funkcí. Například aproximace pily:
$$ f(x) \approx \sum_{k=1}^N \frac{(-1)^k}{k} \sin(2\pi k x)$$ Pila na wiki
In [30]:
def fourier(x,n):
t = 2*np.pi*x # skalar
s = np.array([np.sin(t*k) for k in range(0,n) if k >0])
a = np.array([(-1)**k/k for k in range(0,n) if k >0])
return sum(a*s)
In [31]:
x = np.linspace(0,1.5,1000)
# tady nelze pouzit y = fourier(x), protoze mame jen skalarni argument
y = [fourier(x[i],n=5) for i in range(0,1000)] # volbou n menime presnost aproximace
plt.plot(x,y)
Out[31]:
wiki.
Rozlišením je konvenčně myšlena poloha prvního minima intensity $$I(x) = I_0 \left( \frac{J_1(x)}{x} \right)^2,$$ kde $$x=\frac{2\pi}{\lambda} R \sin \theta$$.
In [32]:
x = np.linspace(0,20,1000)
y = (spec.j1(x)/x)**2
plt.plot(x,y)
Out[32]:
Hledání minima může být snadno nahrazeno řešením rovnice $f(x) = J_1(x) / x = 0$
In [33]:
y = spec.j1(x)/x
In [34]:
plt.grid() # pridani mrizky do obrazku
plt.plot(x,y)
Out[34]:
In [35]:
plt.xlim(3,5) # zmensime interval
plt.grid()
plt.plot(x,y)
Out[35]:
In [36]:
plt.xlim(3.5,4.0) # znova zmensime interval
plt.grid()
plt.plot(x,y)
Out[36]:
In [37]:
plt.xlim(3.8,3.9) # a pokracujem ...
plt.ylim(-0.1,0.1)
plt.grid()
plt.plot(x,y)
Out[37]:
In [38]:
plt.xlim(3.82,3.85) # a porad pokracujem ...
plt.ylim(-0.01,0.01)
plt.grid()
plt.plot(x,y)
Out[38]:
In [39]:
plt.xlim(3.83,3.835) # porad nekoncime...
plt.ylim(-1e-3,1e-3)
plt.grid()
plt.plot(x,y)
Out[39]:
In [40]:
plt.xlim(3.831,3.832) # tak naposledy....
plt.ylim(-1e-4,1e-4)
plt.grid()
plt.plot(x,y)
Out[40]:
Výsledek řešení rovnice $f(x) = J_1(x)/x = 0$ je pro $x=3.8317 \pm 0.0001$.
In [41]:
R = 0.1 # [m]
lam = 550e-9 # [m]
rad = 180 / np.pi # radian je asi 57.3 stupne
In [42]:
theta = 3.8317 / (2*np.pi / lam * (R / 2) ) * rad * 3600
In [43]:
print("Rozlišení dalekohledu o průměru",100*R,"cm ve viditelném světle je",theta,"arcsec.") # vomit
In [44]:
print("Rozlišení dalekohledu o průměru {0:.0f} cm ve viditelném světle je {1:.1f} arcsec.".format(100*R,theta))
Numerické řešení rovnice $f(x) = 0$ je založeno na počátečním odhadu krajů intervalu $a,b$ pro které musí platit $f(a) \, f(b) < 0$. Ke správnému řešení se propracujeme pomocí půlení intervalu a testování předchozí podmínky. Tedy pro bod v polovině $$x = \frac{a + b}{2}$$ a víme, že kořen je v tom intervalu, pokud $f(a) f(x) < 0$ nebo $f(b) f(x) < 0$.
In [45]:
def f(x):
return spec.j1(x)/x
a = 3
b = 5
print(f(a)*f(b))
x = (a + b) / 2
print(x,f(x)*f(a),f(x)*f(b))
b = x
x = (a + b) / 2
print(x,f(x)*f(a),f(x)*f(b))
a = x
x = (a + b) / 2
print(x,f(x)*f(a),f(x)*f(b))
a = x
x = (a + b) / 2
print(x,f(x)*f(a),f(x)*f(b))
b = x
x = (a + b) / 2
print(x,f(x)*f(a),f(x)*f(b))
a = x
x = (a + b) / 2
print(x,f(x)*f(a),f(x)*f(b))
b = x
x = (a + b) / 2
print(x,f(x)*f(a),f(x)*f(b))
a = x
x = (a + b) / 2
print(x,f(x)*f(a),f(x)*f(b))
b = x
x = (a + b) / 2
print(x,f(x)*f(a),f(x)*f(b))
b = x
x = (a + b) / 2
print(x,f(x)*f(a),f(x)*f(b))
a = x
x = (a + b) / 2
print(x,f(x)*f(a),f(x)*f(b))
a = x
x = (a + b) / 2
print(x,f(x)*f(a),f(x)*f(b))
a = x
x = (a + b) / 2
print(x,f(x)*f(a),f(x)*f(b))
b = x
x = (a + b) / 2
print(x,f(x)*f(a),f(x)*f(b))
a = x
x = (a + b) / 2
print(x,f(x)*f(a),f(x)*f(b))
#for k in range(0,17):
# x = (a + b)/2
# if
Potřebujeme nějak učinně opakovat výpočty a rozhodovat se na základě hodnot funkce.
Cyklické opakování je možné zařídit prostřednictvím příkazu
for k in [vektor]:
prikaz(y)
In [46]:
# priklad cyklu
for i in range(0,3):
x = 1 + i**2
print("#",i," x=",x)
Občas je potřeba se rozhodnout v závislosti na hodnotě nějaké podmínky. To je možné prostřednictvím příkazu:
if podminka:
prikaz(y), pokud je splnena
else:
prikaz(y), pokud neni splnena
In [47]:
# priklad podminky
a = 1
b = 2
if (a > b):
print("Podmínka splněna:", a > b)
else:
print("Podmínka nesplněna:",a > b)
In [48]:
a = 3
b = 5
for k in range(0,20):
x = (a + b) / 2
if (f(x)*f(a) < 0):
b = x
else:
a = x
print("#",k," x=",x)
In [49]:
def pulky(f,a,b):
for k in range(0,20):
x = (a + b) / 2
if (f(x)*f(a) < 0):
b = x
else:
a = x
return x
In [50]:
print("Funkce na pulky {0:.4f}.".format(pulky(f,3,5)))
In [51]:
import scipy.optimize
x = scipy.optimize.brentq(f,3,5)
print("Brent's method for x=",x)
In [ ]: