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.
from sympy import *
Da bi dobili lijepi $\LaTeX$ izlaz:
from sympy import init_printing
init_printing()
Koristit ćemo i interaktivne widgete, pa ih ovdje učitavamo
from IPython.display import display
from ipywidgets import interact, fixed, interact_manual
import ipywidgets as widgets
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:
x = Symbol('x')
# ili x,y,z = symbols('x,y,z')
# ili from sympy.abc import x,y,z
# ili var(x:z)
(pi + x)**2
a, b, c = symbols("stranica_a, stranica_b, stranica_c")
type(a)
sympy.core.symbol.Symbol
a
a, b, c = symbols("alpha, beta, gamma")
a**2+b**2+c**2
symbols("x:5")
Možemo navoditi i dodatne pretpostavke:
x = Symbol('x', real=True)
x.is_imaginary
False
x = Symbol('x', positive=True)
x > 0
Možemo kreirati i apstraktne funkcije:
f = Function('f')
f(0)
g = Function('g')(x)
g.diff(x), g.diff(a)
Imaginarna jedinica se označava s I
.
1+1*I
I**2
(x * I + 1)**2
Postoje tri numerička tipa: Real
, Rational
, Integer
:
r1 = Rational(4,5)
r2 = Rational(5,4)
r1
r1+r2
r1/r2
denom(r1)
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.
pi.evalf(n=50)
y = (x + pi)**2
N(y, 5)
Ukoliko želimo zamijeniti varijablu s konkretnim brojem, to možemo učiniti koristeći funkciju subs
:
y.subs(x, 1.5)
N(y.subs(x, 1.5))
No subs
možemo korisiti i općenitije:
y.subs(x, a+pi)
Sympy i Numpy se mogu simultano koristiti:
import numpy
x_vec = numpy.arange(0, 10, 0.1)
y_vec = numpy.array([N(((x + pi)**2).subs(x, xx)) for xx in x_vec])
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:
# prvi argument je lista varijabli funkcije f, u ovom slučaju funckcija je x -> f(x)
f = lambdify([x], (x + pi)**2, 'numpy')
y_vec = f(x_vec)
Razlika u brzini izvođenja:
%%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
%%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:
string = '1/(x-1) + 1/(x+1) + x + 1'
izraz = sympify(string)
izraz
Jedan interaktivan primjer:
x = Symbol('x')
def factorit(n):
return display(Eq(x ** n - 1, factor(x ** n - 1)))
Eq
kreira matematičke jednakosti, tj. jednadžbe.
factorit(18)
interact(factorit,n=(2,20));
interact(factorit,n=(1,20,2));
interact(factorit,n=widgets.widget_int.IntSlider(min=2,max=20,step=1,value=2));
together(izraz)
cancel(together(izraz))
(x+1)*(x+2)*(x+3)
expand((x+1)*(x+2)*(x+3))
expand
prima dodatne argumente. Npr. trig=True
:
sin(a+b)
expand(sin(a+b), trig=True)
simplify(sin(a)**2 + cos(a)**2)
simplify(cos(x)/sin(x))
f1 = 1/((a+1)*(a+2))
apart(f1)
f2 = 1/(a+2) + 1/(a+3)
together(f2)
y
diff(y**2, x)
Više derivacije:
diff(y**2, x, x)
diff(y**2, x, 2)
from sympy.abc import x,y,z
# ili npr. symbols ('x:z')
f = sin(x*y) + cos(y*z)
Želimo izračunati $$\frac{\partial^3f}{\partial x \partial y^2}$$
diff(f, x, 1, y, 2)
def deriv(f):
display(diff(f,x))
interact_manual(deriv, f='x');
f
integrate(f, x)
Definitni integrali:
integrate(f, (x, -1, 1))
Nepravi integrali:
integrate(exp(-x**2), (x, -oo, oo))
n = Symbol("n")
Sum(1/n**2, (n, 1, 10))
Sum(1/n**2, (n,1, 10)).evalf()
Sum(1/n**2, (n, 1, oo)).evalf()
Product(n, (n, 1, 10))
limit(sin(x)/x, x, 0)
f
diff(f, x)
h = Symbol("h")
limit((f.subs(x, x+h) - f)/h, h, 0)
limit(1/x, x, 0, dir="+")
limit(1/x, x, 0, dir="-")
series(exp(x), x)
Rastav oko $x=1$:
series(exp(x), x, 1)
series(exp(x), x, 1, 10)
tan(x).series(x,pi/2)
s1 = cos(x).series(x, 0, 5)
s1
s2 = sin(x).series(x, 0, 2)
s2
expand(s1 * s2)
S metodom removeO
se možemo riješiti $\mathcal{O}$ dijela:
expand(s1.removeO() * s2.removeO())
Ali oprezno s time:
(cos(x)*sin(x)).series(x, 0, 6)
Reziduumi:
residue(2/sin(x), x, 0)
m11, m12, m21, m22 = symbols("m11, m12, m21, m22")
b1, b2 = symbols("b1, b2")
A = Matrix([[m11, m12],[m21, m22]])
A
b = Matrix([[b1], [b2]])
b
A**2
A * b
def funkcija(A,f):
return display(getattr(A,f)())
interact(funkcija,A = fixed(A), f=['det','inv','adjoint','charpoly']);
solve(x**2 - 1, x)
solve(x**4 - x**2 - 1, x)
eq = Eq(x**3 + 2*x**2 + 4*x + 8, 0)
eq
solve(eq, x)
Sustavi jednadžbi:
solve([x + y - 1, x - y - 1], [x,y])
solve([x + y - a, x - y - c], [x,y])
Više o interaktivnim widgetima možete naučiti preko primjera koji se nalaze ovdje.
from verzije import *
from IPython.display import HTML
HTML(print_sysinfo()+info_packages('sympy,matplotlib,IPython,numpy, ipywidgets'))
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 |
def evaluiraj(izrazi, x, x0):
"""
Za svaki izraz iz izrazi funkcija ispisuje vrijednost izraza za x = x0
"""
$\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.$
Izraz (x**2 + 3*x + 1)/(x**3 + 2*x**2 + x)
pojednostavite do izraza 1/(x**2 + 2*x + 1) + 1/x
Izraz (a**b)**c)
pojednostavite do izraza a**(b*c)
Napišite funkciju koja rješava (egzaktno) kvadratnu jednadžbu.
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).
Izračunajte vektorski produkt vektora $(a,b,b)$ i $(b,a,a)$
Riješite Cauchyjev problem
i pojednostavite dobijeno rješenje. Ovdje su $a$, $b$ i $I$ nepoznate konstante.
i nacrtajte partikularna rješenja koristeći sympy modul za crtanje.