CoCalc Shared Files06. Sympy / Sympy.slides.htmlOpen in CoCalc with one click!
Author: Ivica Nakić
Sympy slides

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 [34]:
from sympy import *

Da bi dobili lijepi LaTeX\LaTeX izlaz:

In [35]:
from sympy import init_printing
init_printing()

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

In [42]:
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 [19]:
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]:
(x+π)2\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]:
stranicaastranica_{a}
In [9]:
a, b, c = symbols("alpha, beta, gamma")
a**2+b**2+c**2
Out[9]:
α2+β2+γ2\alpha^{2} + \beta^{2} + \gamma^{2}
In [10]:
symbols("x:5")
Out[10]:
(x0,x1,x2,x3,x4)\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]:
True\mathrm{True}

Možemo kreirati i apstraktne funkcije:

In [15]:
f = Function('f')
f(0)
Out[15]:
f(0)f{\left (0 \right )}
In [16]:
g = Function('g')(x)
g.diff(x), g.diff(a)
Out[16]:
(ddxg(x),0)\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+i1 + i
In [18]:
I**2
Out[18]:
1-1
In [19]:
(x * I + 1)**2
Out[19]:
(ix+1)2\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]:
45\frac{4}{5}
In [22]:
r1+r2
Out[22]:
4120\frac{41}{20}
In [23]:
r1/r2
Out[23]:
1625\frac{16}{25}
In [24]:
denom(r1)
Out[24]:
55

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.14159265358979323846264338327950288419716939937513.1415926535897932384626433832795028841971693993751
In [20]:
y = (x + pi)**2
In [27]:
N(y, 5)
Out[27]:
(x+3.1416)2\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]:
(1.5+π)2\left(1.5 + \pi\right)^{2}
In [29]:
N(y.subs(x, 1.5))
Out[29]:
21.544382361858721.5443823618587

No subs možemo korisiti i općenitije:

In [30]:
y.subs(x, a+pi)
Out[30]:
(α+2π)2\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+1x+1+1x1x + 1 + \frac{1}{x + 1} + \frac{1}{x - 1}

Jedan interaktivan primjer:

In [21]:
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)
x181=(x1)(x+1)(x2x+1)(x2+x+1)(x6x3+1)(x6+x3+1)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 [22]:
interact(factorit,n=(2,20));
x151=(x1)(x2+x+1)(x4+x3+x2+x+1)(x8x7+x5x4+x3x+1)x^{15} - 1 = \left(x - 1\right) \left(x^{2} + x + 1\right) \left(x^{4} + x^{3} + x^{2} + x + 1\right) \left(x^{8} - x^{7} + x^{5} - x^{4} + x^{3} - x + 1\right)
In [23]:
interact(factorit,n=(1,20,2));
x151=(x1)(x2+x+1)(x4+x3+x2+x+1)(x8x7+x5x4+x3x+1)x^{15} - 1 = \left(x - 1\right) \left(x^{2} + x + 1\right) \left(x^{4} + x^{3} + x^{2} + x + 1\right) \left(x^{8} - x^{7} + x^{5} - x^{4} + x^{3} - x + 1\right)
In [24]:
interact(factorit,n=widgets.widget_int.IntSlider(min=2,max=20,step=1,value=2));
x121=(x1)(x+1)(x2+1)(x2x+1)(x2+x+1)(x4x2+1)x^{12} - 1 = \left(x - 1\right) \left(x + 1\right) \left(x^{2} + 1\right) \left(x^{2} - x + 1\right) \left(x^{2} + x + 1\right) \left(x^{4} - x^{2} + 1\right)

Algebarske manipulacije

In [45]:
together(izraz)
Out[45]:
1(x1)(x+1)(x(x1)(x+1)+2x+(x1)(x+1))\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]:
1x21(x3+x2+x1)\frac{1}{x^{2} - 1} \left(x^{3} + x^{2} + x - 1\right)
In [47]:
(x+1)*(x+2)*(x+3)
Out[47]:
(x+1)(x+2)(x+3)\left(x + 1\right) \left(x + 2\right) \left(x + 3\right)
In [48]:
expand((x+1)*(x+2)*(x+3))
Out[48]:
x3+6x2+11x+6x^{3} + 6 x^{2} + 11 x + 6

expand prima dodatne argumente. Npr. trig=True:

In [49]:
sin(a+b)
Out[49]:
sin(α+β)\sin{\left (\alpha + \beta \right )}
In [50]:
expand(sin(a+b), trig=True)
Out[50]:
sin(α)cos(β)+sin(β)cos(α)\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]:
11
In [52]:
simplify(cos(x)/sin(x))
Out[52]:
1tan(x)\frac{1}{\tan{\left (x \right )}}
In [53]:
f1 = 1/((a+1)*(a+2))
In [54]:
apart(f1)
Out[54]:
1α+2+1α+1- \frac{1}{\alpha + 2} + \frac{1}{\alpha + 1}
In [55]:
f2 = 1/(a+2) + 1/(a+3)
In [56]:
together(f2)
Out[56]:
2α+5(α+2)(α+3)\frac{2 \alpha + 5}{\left(\alpha + 2\right) \left(\alpha + 3\right)}

Analiza

Deriviranje

In [6]:
y
Out[6]:
(x+π)2\left(x + \pi\right)^{2}
In [7]:
diff(y**2, x)
Out[7]:
4(x+π)34 \left(x + \pi\right)^{3}

Više derivacije:

In [8]:
diff(y**2, x, x)
Out[8]:
12(x+π)212 \left(x + \pi\right)^{2}
In [9]:
diff(y**2, x, 2)
Out[9]:
12(x+π)212 \left(x + \pi\right)^{2}
In [37]:
from sympy.abc import x,y,z
# ili npr. symbols ('x:z')
In [38]:
f = sin(x*y) + cos(y*z)

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

In [39]:
diff(f, x, 1, y, 2)
Out[39]:
x(xycos(xy)+2sin(xy))- x \left(x y \cos{\left (x y \right )} + 2 \sin{\left (x y \right )}\right)
In [44]:
def deriv(f):
    display(diff(f,x))
interact_manual(deriv, f='x');
11

Integracija

In [65]:
f
Out[65]:
sin(xy)+cos(yz)\sin{\left (x y \right )} + \cos{\left (y z \right )}
In [66]:
integrate(f, x)
Out[66]:
xcos(yz)+{0fory=01ycos(xy)otherwisex \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]:
2cos(yz)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]:
n=1101n2\sum_{n=1}^{10} \frac{1}{n^{2}}
In [71]:
Sum(1/n**2, (n,1, 10)).evalf()
Out[71]:
1.549767731166541.54976773116654
In [72]:
Sum(1/n**2, (n, 1, oo)).evalf()
Out[72]:
1.644934066848231.64493406684823
In [73]:
Product(n, (n, 1, 10))
Out[73]:
n=110n\prod_{n=1}^{10} n

Limesi

In [74]:
limit(sin(x)/x, x, 0)
Out[74]:
11
In [75]:
f
Out[75]:
sin(xy)+cos(yz)\sin{\left (x y \right )} + \cos{\left (y z \right )}
In [76]:
diff(f, x)
Out[76]:
ycos(xy)y \cos{\left (x y \right )}
f(x,y)x=limh0f(x+h,y)f(x,y)h \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]:
ycos(xy)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+x22+x36+x424+x5120+O(x6)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=1x=1:

In [82]:
series(exp(x), x, 1)
Out[82]:
e+e(x1)+e2(x1)2+e6(x1)3+e24(x1)4+e120(x1)5+O((x1)6;x1)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(x1)+e2(x1)2+e6(x1)3+e24(x1)4+e120(x1)5+e720(x1)6+e5040(x1)7+e40320(x1)8+e362880(x1)9+O((x1)10;x1)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]:
1xπ2π6+145(xπ2)3+2945(xπ2)5+x3+O((xπ2)6;xπ2)- \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]:
1x22+x424+O(x5)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+O(x2)x + \mathcal{O}\left(x^{2}\right)
In [87]:
expand(s1 * s2)
Out[87]:
x+O(x2)x + \mathcal{O}\left(x^{2}\right)

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

In [88]:
expand(s1.removeO() * s2.removeO())
Out[88]:
x524x32+x\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]:
x2x33+2x515+O(x6)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]:
22

Linearna algebra

Matrice

In [29]:
m11, m12, m21, m22 = symbols("m11, m12, m21, m22")
b1, b2 = symbols("b1, b2")
In [30]:
A = Matrix([[m11, m12],[m21, m22]])
A
Out[30]:
[m11m12m21m22]\left[\begin{matrix}m_{11} & m_{12}\\m_{21} & m_{22}\end{matrix}\right]
In [31]:
b = Matrix([[b1], [b2]])
b
Out[31]:
[b1b2]\left[\begin{matrix}b_{1}\\b_{2}\end{matrix}\right]
In [32]:
A**2
Out[32]:
[m112+m12m21m11m12+m12m22m11m21+m21m22m12m21+m222]\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]:
[b1m11+b2m12b1m21+b2m22]\left[\begin{matrix}b_{1} m_{11} + b_{2} m_{12}\\b_{1} m_{21} + b_{2} m_{22}\end{matrix}\right]
In [33]:
def funkcija(A,f):
    return display(getattr(A,f)())
interact(funkcija,A = fixed(A), f=['det','inv','adjoint','charpoly']);
m11m22m12m21m_{11} m_{22} - m_{12} m_{21}

Rješavanje jednadžbi

In [99]:
solve(x**2 - 1, x)
Out[99]:
[1,1]\left [ -1, \quad 1\right ]
In [100]:
solve(x**4 - x**2 - 1, x)
Out[100]:
[i12+52,i12+52,12+52,12+52]\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]:
x3+2x2+4x+8=0x^{3} + 2 x^{2} + 4 x + 8 = 0
In [102]:
solve(eq, x)
Out[102]:
[2,2i,2i]\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]:
{x:1,y:0}\left \{ x : 1, \quad y : 0\right \}
In [104]:
solve([x + y - a, x - y - c], [x,y])
Out[104]:
{x:α2+γ2,y:α2γ2}\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 verzija3.5.3
kompajlerGCC 4.8.2 20140120 (Red Hat 4.8.2-15)
sustavLinux
broj CPU-a8
interpreter64bit
sympy verzija1.0
matplotlib verzija2.0.0
IPython verzija5.3.0
numpy verzija1.11.3
ipywidgets verzija6.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 ddxsin(x)ex,xsin(xy)ex,2xysin(xy)ex.\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)(a,b,b) i (b,a,a)(b,a,a)

8) Riješite Cauchyjev problem i pojednostavite dobijeno rješenje. Ovdje su aa, bb i II nepoznate konstante.

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