SciPy je nadgradnja NumPy paketa, i sadrži veliki broj numeričkih algoritama za cijeli niz područja. Ovdje su pobrojana neka nama zanimljivija:
Za zadnja dva područja postoje i napredniji paketi. Za slike smo već npr. koristili scikit-image), a za statistiku ćemo korititi pandas paket.
SciPy paket učitavamo pomoću scipy modula.
In [1]:
from scipy import *
Narvno, možemo učitati i samo podpaket koji nas zanima, u ovom slučaju za linearnu algebru.
In [2]:
import scipy.linalg as la
In [3]:
# jn, yn: Besselove funkcije prvog i drugog reda s realnim stupnjem
# jn_zeros, yn_zeros: računaju pripadne nultočke
from scipy.special import jn, yn, jn_zeros, yn_zeros
In [4]:
n = 0 # stupanj
x = 0.0
print ("J_{}({}) = {:f}".format(n, x, jn(n, x)))
x = 1.0
print ("Y_{}({}) = {:f}".format(n, x, yn(n, x)))
In [2]:
from pylab import *
%matplotlib inline
In [6]:
x = linspace(0, 10, 100)
fig, ax = subplots()
for n in range(4):
ax.plot(x, jn(n, x), label=r"$J_%d(x)$" % n)
ax.legend();
In [7]:
n = 0 # stupanj
m = 4 # broj nultočaka za izračunati
jn_zeros(n, m)
Out[7]:
In [8]:
from scipy.integrate import quad, dblquad, tplquad
quad funkcija se koriste za numeričku integracije (quad ... jer se na engleskom taj proces zove kvadratura). dblquad služi za dvostruke, a tplquad za trostruke integrale.
Jednostavan primjer, računamo $$\begin{equation*} \int_0^1 x\, \mathrm{d}x \end{equation*}$$
In [9]:
def f(x):
return x
In [10]:
x_donje = 0
x_gornje = 1
rez, abserr = quad(f, x_donje, x_gornje)
print ("Rezultat = {}, apsolutna greška = {}".format(rez,abserr))
Ove funkcije imaju puno opcionalnih argumenata. Ako želimo funkciji koju integriramo proslijediti dodatne parametre (vidi ovdje), možemo korisiti varijablu args.
In [11]:
def integrand(x, n):
"""
Besselova funkcija prvog tipa stupnja n.
"""
return jn(n, x)
x_d = 0
x_g = 10
rez, abserr = quad(integrand, x_d, x_g, args=(3,))
print (rez, abserr)
Korištenje anonimnih funkcija:
In [12]:
rez, abserr = quad(lambda x: exp(-x ** 2), -Inf, Inf)
print ("numerički = {}, {}".format(rez, abserr))
egzaktno = sqrt(pi)
print ("egzaktno = {}".format(egzaktno))
U više dimenzija:
\begin{equation*} \int_a^b \int_{g(x)}^{h(x)} f(x,y)\,\mathrm{d}y\mathrm{d}x \end{equation*}
In [13]:
def integrand(x, y):
return exp(-x**2-y**2)
x_d = 0
x_g = 10
y_d = 0
y_g = 10
# ovdje je a = x_d, b = x_g, g(x) = y_d, h(x) = y_g
# g(x) i h(x) trebaju biti funkcije!
rez, abserr = dblquad(integrand, x_d, x_g, lambda x : y_d, lambda x: y_g)
print (rez, abserr)
In [3]:
from scipy.integrate import odeint, ode
Sustav ODJ zapisujemo kao:
$y' = f(y, t)$
gdje je
$y = [y_1(t), y_2(t), ..., y_n(t)]$,
Još trebamo i početne uvjete $y(0)$.
Ovo je sintaksa:
y_t = odeint(f, y_0, t)
t je niz vremena za koje želimo riješiti ODJ y_t je niz s jednim retkom za svaki trenutak iz t, a stupci daje rješenje y_i(t) u tom trenutkuOpis problema: http://en.wikipedia.org/wiki/Double_pendulum
In [15]:
from IPython.display import Image
Image(url='http://upload.wikimedia.org/wikipedia/commons/c/c9/Double-compound-pendulum-dimensioned.svg')
Out[15]:
Ovo su jednadže s wiki stranice:
${\dot \theta_1} = \frac{6}{m\ell^2} \frac{ 2 p_{\theta_1} - 3 \cos(\theta_1-\theta_2) p_{\theta_2}}{16 - 9 \cos^2(\theta_1-\theta_2)}$
${\dot \theta_2} = \frac{6}{m\ell^2} \frac{ 8 p_{\theta_2} - 3 \cos(\theta_1-\theta_2) p_{\theta_1}}{16 - 9 \cos^2(\theta_1-\theta_2)}.$
${\dot p_{\theta_1}} = -\frac{1}{2} m \ell^2 \left [ {\dot \theta_1} {\dot \theta_2} \sin (\theta_1-\theta_2) + 3 \frac{g}{\ell} \sin \theta_1 \right ]$
${\dot p_{\theta_2}} = -\frac{1}{2} m \ell^2 \left [ -{\dot \theta_1} {\dot \theta_2} \sin (\theta_1-\theta_2) + \frac{g}{\ell} \sin \theta_2 \right]$
Definiramo:
$x = [\theta_1, \theta_2, p_{\theta_1}, p_{\theta_2}]$
In [4]:
g = 9.82
L = 0.5
m = 0.1
def dx(x, t):
"""
Desna strana ODJ
"""
x1, x2, x3, x4 = x[0], x[1], x[2], x[3]
dx1 = 6.0/(m*L**2) * (2 * x3 - 3 * cos(x1-x2) * x4)/(16 - 9 * cos(x1-x2)**2)
dx2 = 6.0/(m*L**2) * (8 * x4 - 3 * cos(x1-x2) * x3)/(16 - 9 * cos(x1-x2)**2)
dx3 = -0.5 * m * L**2 * ( dx1 * dx2 * sin(x1-x2) + 3 * (g/L) * sin(x1))
dx4 = -0.5 * m * L**2 * (-dx1 * dx2 * sin(x1-x2) + (g/L) * sin(x2))
return [dx1, dx2, dx3, dx4]
In [5]:
# početni uvjet
x0 = [pi/4, pi/2, 0, 0]
In [6]:
# niz vremena
t = linspace(0, 10, 250)
In [7]:
# rješenje ODJ
x = odeint(dx, x0, t)
In [20]:
# nacrtajmo rješenje
# crtamo kuteve
fig, axes = subplots(1,2, figsize=(12,4))
axes[0].plot(t, x[:, 0], 'r', label="theta1")
axes[0].plot(t, x[:, 1], 'b', label="theta2")
x1 = + L * sin(x[:, 0])
y1 = - L * cos(x[:, 0])
x2 = x1 + L * sin(x[:, 1])
y2 = y1 - L * cos(x[:, 1])
axes[1].plot(x1, y1, 'r', label="njihalo1")
axes[1].plot(x2, y2, 'b', label="njihalo2")
axes[1].set_ylim([-1, 0])
axes[1].set_xlim([1, -1]);
Jednostavna animacija, kasnije ćemo vidjeti kako možemo napravit bolju animaciju.
In [8]:
from IPython.display import display,clear_output
import time
In [9]:
fig, ax = subplots(figsize=(4,4))
for t_idx, tt in enumerate(t[:200]):
x1 = + L * sin(x[t_idx, 0])
y1 = - L * cos(x[t_idx, 0])
x2 = x1 + L * sin(x[t_idx, 1])
y2 = y1 - L * cos(x[t_idx, 1])
ax.cla()
ax.plot([0, x1], [0, y1], 'r.-')
ax.plot([x1, x2], [y1, y2], 'b.-')
ax.set_ylim([-1.5, 0.5])
ax.set_xlim([1, -1])
display(fig)
clear_output(wait=True)
time.sleep(0.03)
Opis problema možete pročitati ovdje: http://en.wikipedia.org/wiki/Damping
Jednadžba je
\begin{equation} \frac{\mathrm{d}^2x}{\mathrm{d}t^2} + 2\zeta\omega_0\frac{\mathrm{d}x}{\mathrm{d}t} + \omega^2_0 x = 0 \end{equation}Definiramo $p = \frac{\mathrm{d}x}{\mathrm{d}t}$:
\begin{equation} \frac{\mathrm{d}p}{\mathrm{d}t} = - 2\zeta\omega_0 p - \omega^2_0 x \end{equation}\begin{equation} \frac{\mathrm{d}x}{\mathrm{d}t} = p \end{equation}
In [23]:
def dy(y, t, zeta, w0):
"""
Desna strana ODJ za harmonički oscilator
"""
x, p = y[0], y[1]
dx = p
dp = -2 * zeta * w0 * p - w0**2 * x
return [dx, dp]
In [24]:
# početno stanje:
y0 = [1.0, 0.0]
In [25]:
# vremena, frekvencija
t = linspace(0, 10, 1000)
w0 = 2*pi*1.0
In [26]:
# rješavamo ODJ za tri vrste prigušenja
y1 = odeint(dy, y0, t, args=(0.0, w0)) # negušeno
y2 = odeint(dy, y0, t, args=(0.2, w0)) # podgušeno
y3 = odeint(dy, y0, t, args=(1.0, w0)) # kritičko gušenje
y4 = odeint(dy, y0, t, args=(5.0, w0)) # pregušeno
In [27]:
fig, ax = subplots()
ax.plot(t, y1[:,0], 'k', label="negušeno", linewidth=0.25)
ax.plot(t, y2[:,0], 'r', label="podgušeno")
ax.plot(t, y3[:,0], 'b', label="kritičko gušenje")
ax.plot(t, y4[:,0], 'g', label="pregušeno")
ax.legend();
In [28]:
from scipy.fftpack import *
Primjenimo Fourierovu transformaciju na prethodni primjer harmoničkog oscilatora.
In [29]:
N = len(t)
dt = t[1]-t[0]
# y2 je rješenje podgušenog harmoničkog oscilatora
F = fft(y2[:,0])
# izračunajmo frekvencije
w = fftfreq(N, dt)
In [30]:
fig, ax = subplots(figsize=(9,3))
ax.plot(w, abs(F));
Kako je signal realan, spektar je simetričan. Stoga nam je dosta nacrtati pozitivne frekvencije.
In [31]:
indeksi = where(w > 0)
w_pos = w[indeksi]
F_pos = F[indeksi]
In [32]:
fig, ax = subplots(figsize=(9,3))
ax.plot(w_pos, abs(F_pos))
ax.set_xlim(0, 5);
Detaljna dokumentacija: http://docs.scipy.org/doc/scipy/reference/linalg.html
Nećemo prolaziti kroz sve funkcije.
In [35]:
A = array([[1,2,-1], [4,5,6], [7,8,9]])
b = array([1,2,3])
In [36]:
x = solve(A, b)
x
Out[36]:
In [37]:
# provjera
dot(A, x) - b
Out[37]:
$A X = B$
In [38]:
A = rand(3,3)
B = rand(3,3)
In [39]:
X = solve(A, B)
X
Out[39]:
In [40]:
# provjera
norm(dot(A, X) - B)
Out[40]:
In [41]:
evals = eigvals(A)
evals
Out[41]:
In [42]:
evals, evecs = eig(A)
In [43]:
evals
Out[43]:
In [44]:
evecs
Out[44]:
Svojstveni vektori su stupci u evecs:
In [45]:
n = 1
norm(dot(A, evecs[:,n]) - evals[n] * evecs[:,n])
Out[45]:
To nije sve, postoje i specijalizirane funkcije, kao npr. eigh za hermitske matrice
In [46]:
# inverz
inv(A)
Out[46]:
In [47]:
# determinanta
det(A)
Out[47]:
In [48]:
# razne norme
norm(A, ord=2), norm(A, ord=Inf)
Out[48]:
Više informacija na http://en.wikipedia.org/wiki/Sparse_matrix
Postoji više formata rijetkih matrica, mi nećemo ulaziti u detalje.
In [49]:
from scipy.sparse import *
In [50]:
# gusta matrica
M = array([[1,0,0,0], [0,3,0,0], [0,1,1,0], [1,0,0,1]]); M
Out[50]:
In [51]:
# pretvorimo je u rijetku matricu
A = csr_matrix(M); A
Out[51]:
In [52]:
# vratimo natrag
A.todense()
Out[52]:
Pametniji način kreiranja rijetke matrice.
In [53]:
A = lil_matrix((4,4)) # prazna 4x4 rijetka matrica
A[0,0] = 1
A[1,1] = 3
A[2,2] = A[2,1] = 1
A[3,3] = A[3,0] = 1
A
Out[53]:
In [54]:
A.todense()
Out[54]:
In [55]:
# konvertiranje
A = csr_matrix(A); A
Out[55]:
In [56]:
A = csc_matrix(A); A
Out[56]:
In [57]:
A.todense()
Out[57]:
In [58]:
(A * A).todense()
Out[58]:
In [59]:
(A @ A).todense()
Out[59]:
In [61]:
dot(A,A)
Out[61]:
In [62]:
v = array([1,2,3,4])[:,newaxis]; v
Out[62]:
Vektor v smo mogli konstruirati i drugačije (vidjeli smo primjere u predavanju o NumPy-ju), no uvijek trebamo doći do dvodimenzionalnog niza. Za razliku od MATLAB-a u kojemu su svi nizovi 2D, u NumPy-ju 1D niz nije isto što i matrica $n\times 1$ ili $1\times n$.
Npr. jedna mogućnost je
v = array([[1,2,3,4]]).T
In [63]:
# rijetka matrica puta vektor
A * v
Out[63]:
In [64]:
A.todense() * v
Out[64]:
Više na http://scipy-lectures.github.com/advanced/mathematical_optimization/index.html
Modul je optimize:
In [65]:
from scipy import optimize
In [66]:
def f(x):
return 4*x**3 + (x-2)**2 + x**4
In [67]:
fig, ax = subplots()
x = linspace(-5, 3, 100)
ax.plot(x, f(x));
In [68]:
x_min = optimize.fmin_bfgs(f, -2)
x_min
Out[68]:
In [69]:
optimize.fmin_bfgs(f, 0.5)
Out[69]:
In [70]:
optimize.brent(f)
Out[70]:
In [71]:
optimize.fminbound(f, -4, 2)
Out[71]:
In [72]:
omega_c = 3.0
def f(omega):
return tan(2*pi*omega) - omega_c/omega
In [73]:
import numpy as np
np.seterr(divide='ignore')
fig, ax = subplots(figsize=(10,4))
x = linspace(0, 3, 1000)
y = f(x)
maska = where(abs(y) > 50)
x[maska] = y[maska] = NaN # da se riješimo asimptote
ax.plot(x, y)
ax.plot([0, 3], [0, 0], 'k')
ax.set_ylim(-5,5);
In [74]:
optimize.fsolve(f, 0.1)
Out[74]:
In [75]:
optimize.fsolve(f, 0.6)
Out[75]:
In [76]:
optimize.fsolve(f, 1.1)
Out[76]:
In [77]:
from scipy.interpolate import *
In [78]:
n = arange(0, 10)
x = linspace(0, 9, 100)
y_meas = sin(n) + 0.1 * randn(len(n)) # ubacujemo malo šuma
y_real = sin(x)
linear_interpolation = interp1d(n, y_meas)
y_interp1 = linear_interpolation(x)
cubic_interpolation = interp1d(n, y_meas, kind='cubic')
y_interp2 = cubic_interpolation(x)
In [79]:
fig, ax = subplots(figsize=(10,4))
ax.plot(n, y_meas, 'bs', label='podaci sa šumom')
ax.plot(x, y_real, 'k', lw=2, label='originalna funkcija')
ax.plot(x, y_interp1, 'r', label='linearna interpolacija')
ax.plot(x, y_interp2, 'g', label='kubična interpolacija')
ax.legend(loc=3);
Više na http://docs.scipy.org/doc/scipy/reference/stats.html.
Mi ćemo kasnije raditi s moćnijim paketom pandas.
In [80]:
from scipy import stats
In [81]:
# slučajna varijabla s Poissionovom distribucijom
X = stats.poisson(3.5)
In [82]:
n = arange(0,15)
fig, axes = subplots(2,1, sharex=True)
# kumulativna distribucija (CDF)
axes[0].step(n, X.cdf(n))
# histogram 1000 slučajnih realizacija od X
axes[1].hist(X.rvs(size=1000));
In [83]:
# normalna distribucija
Y = stats.norm()
In [84]:
x = linspace(-5,5,100)
fig, axes = subplots(3,1, sharex=True)
# PDF
axes[0].plot(x, Y.pdf(x))
# CDF
axes[1].plot(x, Y.cdf(x));
# histogram
axes[2].hist(Y.rvs(size=1000), bins=50);
Osnovna statistika:
In [85]:
X.mean(), X.std(), X.var()
Out[85]:
In [86]:
Y.mean(), Y.std(), Y.var()
Out[86]:
In [87]:
from verzije import *
from IPython.display import HTML
HTML(print_sysinfo()+info_packages('numpy,scipy,matplotlib'))
Out[87]:
lil_matrix, konvertirajte je u CSR format i riješite $A x = b$ za neki $b$.Slika moonlanding.png je puna šuma. Zadaća je očistiti sliku koristeći Fourierovu transformaciju.
pylab.imread().scipy.fftpack i nacrtajte spektar slike.