Sympy

Sympy je Python biblioteka za simboličku matematiku. Prednost Sympy-ja je što je potpuno napisan u Pythonu (što je katkad i mana). Mi ćemo u nastavku kolegiju obraditi i puno moćniji Sage, koji je CAS u klasi Mathematice i Maplea. No Sage nije biblioteka u Pythonu, već CAS koji koristi Python kao programski jezik.

Korištenje Sympy-ja počinje kao i kod ostalih biblioteka, s importiranjem.



In [45]:

    
from sympy import *

Da bi dobili lijepi $\LaTeX$ izlaz:



In [46]:

    
from sympy import init_printing
init_printing()

Koristit ćemo i interaktivne widgete, pa ih ovdje učitavamo



In [47]:

    
from IPython.display import display
from ipywidgets import interact, fixed, interact_manual
import ipywidgets as widgets

Simboličke varijable

Kako je Sympy samo Python paket, trebamo deklarirati koje simbole ćemo koristiti kao simboličke vatrijable. To možemo napraviti na više načina:



In [48]:

    
x = Symbol('x')
# ili x,y,z = symbols('x,y,z')
# ili from sympy.abc import x,y,z
# ili var(x:z)



In [5]:

    
(pi + x)**2









    Out[5]:





$$\left(x + \pi\right)^{2}$$



In [6]:

    
a, b, c = symbols("stranica_a, stranica_b, stranica_c")



In [7]:

    
type(a)









    Out[7]:





sympy.core.symbol.Symbol



In [8]:

    
a









    Out[8]:





$$stranica_{a}$$



In [9]:

    
a, b, c = symbols("alpha, beta, gamma")
a**2+b**2+c**2









    Out[9]:





$$\alpha^{2} + \beta^{2} + \gamma^{2}$$



In [10]:

    
symbols("x:5")









    Out[10]:





$$\left ( x_{0}, \quad x_{1}, \quad x_{2}, \quad x_{3}, \quad x_{4}\right )$$

Možemo navoditi i dodatne pretpostavke:



In [11]:

    
x = Symbol('x', real=True)



In [12]:

    
x.is_imaginary









    Out[12]:





False



In [13]:

    
x = Symbol('x', positive=True)



In [14]:

    
x > 0









    Out[14]:





$$\mathrm{True}$$

Možemo kreirati i apstraktne funkcije:



In [15]:

    
f = Function('f')
f(0)









    Out[15]:





$$f{\left (0 \right )}$$



In [16]:

    
g = Function('g')(x)
g.diff(x), g.diff(a)









    Out[16]:





$$\left ( \frac{d}{d x} g{\left (x \right )}, \quad 0\right )$$

Kompleksni brojevi

Imaginarna jedinica se označava s I.



In [17]:

    
1+1*I









    Out[17]:





$$1 + i$$



In [18]:

    
I**2









    Out[18]:





$$-1$$



In [19]:

    
(x * I + 1)**2









    Out[19]:





$$\left(i x + 1\right)^{2}$$

Razlomci

Postoje tri numerička tipa: Real, Rational, Integer:



In [20]:

    
r1 = Rational(4,5)
r2 = Rational(5,4)



In [21]:

    
r1









    Out[21]:





$$\frac{4}{5}$$



In [22]:

    
r1+r2









    Out[22]:





$$\frac{41}{20}$$



In [23]:

    
r1/r2









    Out[23]:





$$\frac{16}{25}$$



In [24]:

    
denom(r1)









    Out[24]:





$$5$$

Numerička evaluacija

SymPy može računati u proizvoljnoj točnosti te ima predefinirane matematičke konstante kao: pi, e te oo za beskonačnost.

Funkcija evalf ili metoda N s ulaznom varijablom n računaju izraz na n decimala.



In [25]:

    
pi.evalf(n=50)









    Out[25]:





$$3.1415926535897932384626433832795028841971693993751$$



In [49]:

    
y = (x + pi)**2



In [27]:

    
N(y, 5)









    Out[27]:





$$\left(x + 3.1416\right)^{2}$$

Ukoliko želimo zamijeniti varijablu s konkretnim brojem, to možemo učiniti koristeći funkciju subs:



In [28]:

    
y.subs(x, 1.5)









    Out[28]:





$$\left(1.5 + \pi\right)^{2}$$



In [29]:

    
N(y.subs(x, 1.5))









    Out[29]:





$$21.5443823618587$$

No subs možemo korisiti i općenitije:



In [30]:

    
y.subs(x, a+pi)









    Out[30]:





$$\left(\alpha + 2 \pi\right)^{2}$$

Sympy i Numpy se mogu simultano koristiti:



In [31]:

    
import numpy



In [32]:

    
x_vec = numpy.arange(0, 10, 0.1)



In [33]:

    
y_vec = numpy.array([N(((x + pi)**2).subs(x, xx)) for xx in x_vec])



In [34]:

    
from matplotlib.pyplot import subplots
%matplotlib inline
fig, ax = subplots()
ax.plot(x_vec, y_vec);

Efikasniji kod se postiže funkcijom lambdify koja kompajlira Sympy izraz u funkciju:



In [35]:

    
# prvi argument je lista varijabli funkcije f, u ovom slučaju funckcija je x -> f(x)
f = lambdify([x], (x + pi)**2, 'numpy')



In [36]:

    
y_vec = f(x_vec)

Razlika u brzini izvođenja:



In [37]:

    
%%timeit

y_vec = numpy.array([N(((x + pi)**2).subs(x, xx)) for xx in x_vec])









    



10 loops, best of 3: 22.6 ms per loop



In [38]:

    
%%timeit

y_vec = f(x_vec)









    



The slowest run took 19.02 times longer than the fastest. This could mean that an intermediate result is being cached.
1000000 loops, best of 3: 1.73 µs per loop

Ovdje smo mogli koristiti i theano ili uFuncify.

Pretvaranje stringa u Sympy izraz:



In [39]:

    
string = '1/(x-1) + 1/(x+1) + x + 1'
izraz = sympify(string)
izraz









    Out[39]:





$$x + 1 + \frac{1}{x + 1} + \frac{1}{x - 1}$$

Jedan interaktivan primjer:



In [50]:

    
x = Symbol('x')
def factorit(n):
    return display(Eq(x ** n - 1, factor(x ** n - 1)))

Eq kreira matematičke jednakosti, tj. jednadžbe.



In [41]:

    
factorit(18)









    





$$x^{18} - 1 = \left(x - 1\right) \left(x + 1\right) \left(x^{2} - x + 1\right) \left(x^{2} + x + 1\right) \left(x^{6} - x^{3} + 1\right) \left(x^{6} + x^{3} + 1\right)$$



In [51]:

    
interact(factorit,n=(2,20));









    





$$x^{17} - 1 = \left(x - 1\right) \left(x^{16} + x^{15} + x^{14} + x^{13} + x^{12} + x^{11} + x^{10} + x^{9} + x^{8} + x^{7} + x^{6} + x^{5} + x^{4} + x^{3} + x^{2} + x + 1\right)$$



In [52]:

    
interact(factorit,n=(1,20,2));









    





$$x^{11} - 1 = \left(x - 1\right) \left(x^{10} + x^{9} + x^{8} + x^{7} + x^{6} + x^{5} + x^{4} + x^{3} + x^{2} + x + 1\right)$$



In [53]:

    
interact(factorit,n=widgets.widget_int.IntSlider(min=2,max=20,step=1,value=2));









    





$$x^{9} - 1 = \left(x - 1\right) \left(x^{2} + x + 1\right) \left(x^{6} + x^{3} + 1\right)$$

Algebarske manipulacije



In [45]:

    
together(izraz)









    Out[45]:





$$\frac{1}{\left(x - 1\right) \left(x + 1\right)} \left(x \left(x - 1\right) \left(x + 1\right) + 2 x + \left(x - 1\right) \left(x + 1\right)\right)$$



In [46]:

    
cancel(together(izraz))









    Out[46]:





$$\frac{1}{x^{2} - 1} \left(x^{3} + x^{2} + x - 1\right)$$



In [47]:

    
(x+1)*(x+2)*(x+3)









    Out[47]:





$$\left(x + 1\right) \left(x + 2\right) \left(x + 3\right)$$



In [48]:

    
expand((x+1)*(x+2)*(x+3))









    Out[48]:





$$x^{3} + 6 x^{2} + 11 x + 6$$

expand prima dodatne argumente. Npr. trig=True:



In [49]:

    
sin(a+b)









    Out[49]:





$$\sin{\left (\alpha + \beta \right )}$$



In [50]:

    
expand(sin(a+b), trig=True)









    Out[50]:





$$\sin{\left (\alpha \right )} \cos{\left (\beta \right )} + \sin{\left (\beta \right )} \cos{\left (\alpha \right )}$$



In [51]:

    
simplify(sin(a)**2 + cos(a)**2)









    Out[51]:





$$1$$



In [52]:

    
simplify(cos(x)/sin(x))









    Out[52]:





$$\frac{1}{\tan{\left (x \right )}}$$



In [53]:

    
f1 = 1/((a+1)*(a+2))



In [54]:

    
apart(f1)









    Out[54]:





$$- \frac{1}{\alpha + 2} + \frac{1}{\alpha + 1}$$



In [55]:

    
f2 = 1/(a+2) + 1/(a+3)



In [56]:

    
together(f2)









    Out[56]:





$$\frac{2 \alpha + 5}{\left(\alpha + 2\right) \left(\alpha + 3\right)}$$

Analiza

Deriviranje



In [6]:

    
y









    Out[6]:





$$\left(x + \pi\right)^{2}$$



In [7]:

    
diff(y**2, x)









    Out[7]:





$$4 \left(x + \pi\right)^{3}$$

Više derivacije:



In [8]:

    
diff(y**2, x, x)









    Out[8]:





$$12 \left(x + \pi\right)^{2}$$



In [9]:

    
diff(y**2, x, 2)









    Out[9]:





$$12 \left(x + \pi\right)^{2}$$



In [54]:

    
from sympy.abc import x,y,z
# ili npr. symbols ('x:z')



In [55]:

    
f = sin(x*y) + cos(y*z)

Želimo izračunati $$\frac{\partial^3f}{\partial x \partial y^2}$$



In [39]:

    
diff(f, x, 1, y, 2)









    Out[39]:





$$- x \left(x y \cos{\left (x y \right )} + 2 \sin{\left (x y \right )}\right)$$



In [56]:

    
def deriv(f):
    display(diff(f,x))
interact_manual(deriv, f='x');









    





$$- \sin{\left (x \right )}$$

Integracija



In [65]:

    
f









    Out[65]:





$$\sin{\left (x y \right )} + \cos{\left (y z \right )}$$



In [66]:

    
integrate(f, x)









    Out[66]:





$$x \cos{\left (y z \right )} + \begin{cases} 0 & \text{for}\: y = 0 \\- \frac{1}{y} \cos{\left (x y \right )} & \text{otherwise} \end{cases}$$

Definitni integrali:



In [67]:

    
integrate(f, (x, -1, 1))









    Out[67]:





$$2 \cos{\left (y z \right )}$$

Nepravi integrali:



In [68]:

    
integrate(exp(-x**2), (x, -oo, oo))









    Out[68]:





$$\sqrt{\pi}$$

Sume i produkti



In [69]:

    
n = Symbol("n")



In [70]:

    
Sum(1/n**2, (n, 1, 10))









    Out[70]:





$$\sum_{n=1}^{10} \frac{1}{n^{2}}$$



In [71]:

    
Sum(1/n**2, (n,1, 10)).evalf()









    Out[71]:





$$1.54976773116654$$



In [72]:

    
Sum(1/n**2, (n, 1, oo)).evalf()









    Out[72]:





$$1.64493406684823$$



In [73]:

    
Product(n, (n, 1, 10))









    Out[73]:





$$\prod_{n=1}^{10} n$$

Limesi



In [74]:

    
limit(sin(x)/x, x, 0)









    Out[74]:





$$1$$



In [75]:

    
f









    Out[75]:





$$\sin{\left (x y \right )} + \cos{\left (y z \right )}$$



In [76]:

    
diff(f, x)









    Out[76]:





$$y \cos{\left (x y \right )}$$

$$ \frac{\partial f(x,y)}{\partial x} = \lim_{h\to 0}\frac{f(x+h,y)-f(x,y)}{h}$$



In [77]:

    
h = Symbol("h")



In [78]:

    
limit((f.subs(x, x+h) - f)/h, h, 0)









    Out[78]:





$$y \cos{\left (x y \right )}$$



In [79]:

    
limit(1/x, x, 0, dir="+")









    Out[79]:





$$\infty$$



In [80]:

    
limit(1/x, x, 0, dir="-")









    Out[80]:





$$-\infty$$

(Taylorovi) redovi



In [81]:

    
series(exp(x), x)









    Out[81]:





$$1 + x + \frac{x^{2}}{2} + \frac{x^{3}}{6} + \frac{x^{4}}{24} + \frac{x^{5}}{120} + \mathcal{O}\left(x^{6}\right)$$

Rastav oko $x=1$:



In [82]:

    
series(exp(x), x, 1)









    Out[82]:





$$e + e \left(x - 1\right) + \frac{e}{2} \left(x - 1\right)^{2} + \frac{e}{6} \left(x - 1\right)^{3} + \frac{e}{24} \left(x - 1\right)^{4} + \frac{e}{120} \left(x - 1\right)^{5} + \mathcal{O}\left(\left(x - 1\right)^{6}; x\rightarrow1\right)$$



In [83]:

    
series(exp(x), x, 1, 10)









    Out[83]:





$$e + e \left(x - 1\right) + \frac{e}{2} \left(x - 1\right)^{2} + \frac{e}{6} \left(x - 1\right)^{3} + \frac{e}{24} \left(x - 1\right)^{4} + \frac{e}{120} \left(x - 1\right)^{5} + \frac{e}{720} \left(x - 1\right)^{6} + \frac{e}{5040} \left(x - 1\right)^{7} + \frac{e}{40320} \left(x - 1\right)^{8} + \frac{e}{362880} \left(x - 1\right)^{9} + \mathcal{O}\left(\left(x - 1\right)^{10}; x\rightarrow1\right)$$



In [84]:

    
tan(x).series(x,pi/2)









    Out[84]:





$$- \frac{1}{x - \frac{\pi}{2}} - \frac{\pi}{6} + \frac{1}{45} \left(x - \frac{\pi}{2}\right)^{3} + \frac{2}{945} \left(x - \frac{\pi}{2}\right)^{5} + \frac{x}{3} + \mathcal{O}\left(\left(x - \frac{\pi}{2}\right)^{6}; x\rightarrow\frac{\pi}{2}\right)$$



In [85]:

    
s1 = cos(x).series(x, 0, 5)
s1









    Out[85]:





$$1 - \frac{x^{2}}{2} + \frac{x^{4}}{24} + \mathcal{O}\left(x^{5}\right)$$



In [86]:

    
s2 = sin(x).series(x, 0, 2)
s2









    Out[86]:





$$x + \mathcal{O}\left(x^{2}\right)$$



In [87]:

    
expand(s1 * s2)









    Out[87]:





$$x + \mathcal{O}\left(x^{2}\right)$$

S metodom removeO se možemo riješiti $\mathcal{O}$ dijela:



In [88]:

    
expand(s1.removeO() * s2.removeO())









    Out[88]:





$$\frac{x^{5}}{24} - \frac{x^{3}}{2} + x$$

Ali oprezno s time:



In [89]:

    
(cos(x)*sin(x)).series(x, 0, 6)









    Out[89]:





$$x - \frac{2 x^{3}}{3} + \frac{2 x^{5}}{15} + \mathcal{O}\left(x^{6}\right)$$

Reziduumi:



In [90]:

    
residue(2/sin(x), x, 0)









    Out[90]:





$$2$$

Linearna algebra

Matrice



In [57]:

    
m11, m12, m21, m22 = symbols("m11, m12, m21, m22")
b1, b2 = symbols("b1, b2")



In [58]:

    
A = Matrix([[m11, m12],[m21, m22]])
A









    Out[58]:





$$\left[\begin{matrix}m_{11} & m_{12}\\m_{21} & m_{22}\end{matrix}\right]$$



In [31]:

    
b = Matrix([[b1], [b2]])
b









    Out[31]:





$$\left[\begin{matrix}b_{1}\\b_{2}\end{matrix}\right]$$



In [32]:

    
A**2









    Out[32]:





$$\left[\begin{matrix}m_{11}^{2} + m_{12} m_{21} & m_{11} m_{12} + m_{12} m_{22}\\m_{11} m_{21} + m_{21} m_{22} & m_{12} m_{21} + m_{22}^{2}\end{matrix}\right]$$



In [14]:

    
A * b









    Out[14]:





$$\left[\begin{matrix}b_{1} m_{11} + b_{2} m_{12}\\b_{1} m_{21} + b_{2} m_{22}\end{matrix}\right]$$



In [59]:

    
def funkcija(A,f):
    return display(getattr(A,f)())
interact(funkcija,A = fixed(A), f=['det','inv','adjoint','charpoly']);









    





$$m_{11} m_{22} - m_{12} m_{21}$$

Rješavanje jednadžbi



In [99]:

    
solve(x**2 - 1, x)









    Out[99]:





$$\left [ -1, \quad 1\right ]$$



In [100]:

    
solve(x**4 - x**2 - 1, x)









    Out[100]:





$$\left [ - i \sqrt{- \frac{1}{2} + \frac{\sqrt{5}}{2}}, \quad i \sqrt{- \frac{1}{2} + \frac{\sqrt{5}}{2}}, \quad - \sqrt{\frac{1}{2} + \frac{\sqrt{5}}{2}}, \quad \sqrt{\frac{1}{2} + \frac{\sqrt{5}}{2}}\right ]$$



In [101]:

    
eq = Eq(x**3 + 2*x**2 + 4*x + 8, 0)
eq









    Out[101]:





$$x^{3} + 2 x^{2} + 4 x + 8 = 0$$



In [102]:

    
solve(eq, x)









    Out[102]:





$$\left [ -2, \quad - 2 i, \quad 2 i\right ]$$

Sustavi jednadžbi:



In [103]:

    
solve([x + y - 1, x - y - 1], [x,y])









    Out[103]:





$$\left \{ x : 1, \quad y : 0\right \}$$



In [104]:

    
solve([x + y - a, x - y - c], [x,y])









    Out[104]:





$$\left \{ x : \frac{\alpha}{2} + \frac{\gamma}{2}, \quad y : \frac{\alpha}{2} - \frac{\gamma}{2}\right \}$$

Više o interaktivnim widgetima možete naučiti preko primjera koji se nalaze ovdje.



In [105]:

    
from verzije import *
from IPython.display import HTML
HTML(print_sysinfo()+info_packages('sympy,matplotlib,IPython,numpy, ipywidgets'))









    Out[105]:




Python verzija 3.5.3
kompajler GCC 4.8.2 20140120 (Red Hat 4.8.2-15)
sustav Linux
broj CPU-a 8
interpreter 64bit
sympy verzija 1.0
matplotlib verzija 2.0.0
IPython verzija 5.3.0
numpy verzija 1.11.3
 ipywidgets verzija 6.0.0

Zadaci za vježbanje

1) Napišite funkciju koja prima listu izraza, varijablu i točku, a ispisuje vrijednost svakog izraza s liste u toj točki. Funkcija treba izgledati ovako

def evaluiraj(izrazi, x, x0):
"""
Za svaki izraz iz izrazi funkcija ispisuje vrijednost izraza za x = x0
"""

2) Izračunajte $\frac{d}{dx}\sin(x)e^x,\, \frac{\partial}{\partial x}\sin(xy)e^x,\, \frac{\partial^2}{\partial x\partial y}\sin(xy)e^x.$

3) Izraz (x**2 + 3*x + 1)/(x**3 + 2*x**2 + x) pojednostavite do izraza 1/(x**2 + 2*x + 1) + 1/x

4) Izraz (a**b)**c) pojednostavite do izraza a**(b*c)

5) Napišite funkciju koja rješava (egzaktno) kvadratnu jednadžbu.

6) Napišite funkciju koja za ulazne parametre prima funkciju (danu simboličkim izrazom) te točku (tj. broj) a crta danu funkciju i njenu tangentu u danoj točki. Učinite funkciju interaktivnom na način da se može birati funkcija (upisivanjem u polje) te točka (klizačem).

7) Izračunajte vektorski produkt vektora $(a,b,b)$ i $(b,a,a)$

8) Riješite Cauchyjev problem $$\begin{align} u'(t)&=-au(t)+b,\\ u(0)&=I \end{align}$$ i pojednostavite dobijeno rješenje. Ovdje su $a$, $b$ i $I$ nepoznate konstante.

9) Riješite diferencijalnu jednadžbu $$ y(t) + \frac{d^2}{dt^2}y(t) = 0$$ i nacrtajte partikularna rješenja koristeći sympy modul za crtanje.

Python verzija	3.5.3
kompajler	GCC 4.8.2 20140120 (Red Hat 4.8.2-15)
sustav	Linux
broj CPU-a	8
interpreter	64bit

sympy verzija	1.0
matplotlib verzija	2.0.0
IPython verzija	5.3.0
numpy verzija	1.11.3
ipywidgets verzija	6.0.0