питон

SymPy

SymPy - это пакет для символьных вычислений на питоне, подобный системе Mathematica. Он работает с выражениями, содержащими символы.

In [1]:
from sympy import *
init_printing()

Основными кирпичиками, из которых строятся выражения, являются символы. Символ имеет имя, которое используется при печати выражений. Объекты класса Symbol нужно создавать и присваивать переменным питона, чтобы их можно было использовать. В принципе, имя символа и имя переменной, которой мы присваиваем этот символ - две независимые вещи, и можно написать abc=Symbol('xyz'). Но тогда при вводе программы Вы будете использовать abc, а при печати результатов SymPy будет использовать xyz, что приведёт к ненужной путанице. Поэтому лучше, чтобы имя символа совпадало с именем переменной питона, которой он присваивается.

В языках, специально предназначенных для символьных вычислений, таких, как Mathematica, если Вы используете переменную, которой ничего не было присвоено, то она автоматически воспринимается как символ с тем же именем. Питон не был изначально предназначен для символьных вычислений. Если Вы используете переменную, которой ничего не было присвоено, Вы получите сообщение об ошибке. Объекты типа Symbol нужно создавать явно.

In [2]:
x=Symbol('x')
In [3]:
a=x**2-1
a
Out[3]:
$$x^{2} - 1$$

Можно определить несколько символов одновременно. Строка разбивается на имена по пробелам.

In [4]:
y,z=symbols('y z')

Подставим вместо $x$ выражение $y+1$.

In [5]:
a.subs(x,y+1)
Out[5]:
$$\left(y + 1\right)^{2} - 1$$

Многочлены и рациональные функции

SymPy не раскрывает скобки автоматически. Для этого используется функция expand.

In [6]:
a=(x+y-z)**6
a
Out[6]:
$$\left(x + y - z\right)^{6}$$
In [7]:
a=expand(a)
a
Out[7]:
$$x^{6} + 6 x^{5} y - 6 x^{5} z + 15 x^{4} y^{2} - 30 x^{4} y z + 15 x^{4} z^{2} + 20 x^{3} y^{3} - 60 x^{3} y^{2} z + 60 x^{3} y z^{2} - 20 x^{3} z^{3} + 15 x^{2} y^{4} - 60 x^{2} y^{3} z + 90 x^{2} y^{2} z^{2} - 60 x^{2} y z^{3} + 15 x^{2} z^{4} + 6 x y^{5} - 30 x y^{4} z + 60 x y^{3} z^{2} - 60 x y^{2} z^{3} + 30 x y z^{4} - 6 x z^{5} + y^{6} - 6 y^{5} z + 15 y^{4} z^{2} - 20 y^{3} z^{3} + 15 y^{2} z^{4} - 6 y z^{5} + z^{6}$$

Степень многочлена $a$ по $x$.

In [8]:
degree(a,x)
Out[8]:
$$6$$

Соберём вместе члены с определёнными степенями $x$.

In [9]:
collect(a,x)
Out[9]:
$$x^{6} + x^{5} \left(6 y - 6 z\right) + x^{4} \left(15 y^{2} - 30 y z + 15 z^{2}\right) + x^{3} \left(20 y^{3} - 60 y^{2} z + 60 y z^{2} - 20 z^{3}\right) + x^{2} \left(15 y^{4} - 60 y^{3} z + 90 y^{2} z^{2} - 60 y z^{3} + 15 z^{4}\right) + x \left(6 y^{5} - 30 y^{4} z + 60 y^{3} z^{2} - 60 y^{2} z^{3} + 30 y z^{4} - 6 z^{5}\right) + y^{6} - 6 y^{5} z + 15 y^{4} z^{2} - 20 y^{3} z^{3} + 15 y^{2} z^{4} - 6 y z^{5} + z^{6}$$

Многочлен с целыми коэффициентами можно записать в виде произведения таких многочленов (причём каждый сомножитель уже невозможно расфакторизовать дальше, оставаясь в рамках многочленов с целыми коэффициентами). Существуют эффективные алгоритмы для решения этой задачи.

In [10]:
a=factor(a)
a
Out[10]:
$$\left(x + y - z\right)^{6}$$

SymPy не сокращает отношения многочленов на их наибольший общий делитель автоматически. Для этого используется функция cancel.

In [11]:
a=(x**3-y**3)/(x**2-y**2)
a
Out[11]:
$$\frac{x^{3} - y^{3}}{x^{2} - y^{2}}$$
In [12]:
cancel(a)
Out[12]:
$$\frac{x^{2} + x y + y^{2}}{x + y}$$

SymPy не приводит суммы рациональных выражений к общему знаменателю автоматически. Для этого используется функция together.

In [13]:
a=y/(x-y)+x/(x+y)
a
Out[13]:
$$\frac{x}{x + y} + \frac{y}{x - y}$$
In [14]:
together(a)
Out[14]:
$$\frac{x \left(x - y\right) + y \left(x + y\right)}{\left(x - y\right) \left(x + y\right)}$$

Функция simplify пытается переписать выражение в наиболее простом виде. Это понятие не имеет чёткого определения (в разных ситуациях наиболее простыми могут считаться разные формы выражения), и не существует алгоритма такого упрощения. Функция symplify работает эвристически, и невозможно заранее предугадать, какие упрощения она попытается сделать. Поэтому её удобно использовать в интерактивных сессиях, чтобы посмотреть, удастся ли ей записать выражение в каком-нибудь разумном виде, но нежелательно использовать в программах. В них лучше применять более специализированные функции, которые выполняют одно определённое преобразование выражения.

In [15]:
simplify(a)
Out[15]:
$$\frac{x^{2} + y^{2}}{x^{2} - y^{2}}$$

Разложение на элементарные дроби по отношению к $x$.

In [16]:
apart(a,x)
Out[16]:
$$- \frac{y}{x + y} + \frac{y}{x - y} + 1$$

Подставим конкретные численные значения вместо переменных $x$ и $y$.

In [17]:
a=a.subs({x:1,y:2})
a
Out[17]:
$$- \frac{5}{3}$$

А сколько это будет численно?

In [18]:
a.n()
Out[18]:
$$-1.66666666666667$$

Элементарные функции

SymPy автоматически применяет упрощения элементарных функция (которые справедливы во всех случаях).

In [19]:
sin(-x)
Out[19]:
$$- \sin{\left (x \right )}$$
In [20]:
cos(pi/4),tan(5*pi/6)
Out[20]:
$$\left ( \frac{\sqrt{2}}{2}, \quad - \frac{\sqrt{3}}{3}\right )$$

SymPy может работать с числами с плавающей точкой, имеющими сколь угодно большую точность. Вот $\pi$ с 100 значащими цифрами.

In [21]:
pi.n(100)
Out[21]:
$$3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068$$

E - это основание натуральных логарифмов.

In [22]:
log(1),log(E)
Out[22]:
$$\left ( 0, \quad 1\right )$$
In [23]:
exp(log(x)),log(exp(x))
Out[23]:
$$\left ( x, \quad \log{\left (e^{x} \right )}\right )$$

А почему не $x$? Попробуйте подставить $x=2\pi i$.

In [24]:
sqrt(0)
Out[24]:
$$0$$
In [25]:
sqrt(x)**4,sqrt(x**4)
Out[25]:
$$\left ( x^{2}, \quad \sqrt{x^{4}}\right )$$

А почему не $x^2$? Попробуйте подставить $x=i$.

Символы могут иметь некоторые свойства. Например, они могут быть положительными. Тогда SymPy может сильнее упростить квадратные корни.

In [26]:
p,q=symbols('p q',positive=True)
sqrt(p**2)
Out[26]:
$$p$$
In [27]:
sqrt(12*x**2*y),sqrt(12*p**2*y)
Out[27]:
$$\left ( 2 \sqrt{3} \sqrt{x^{2} y}, \quad 2 \sqrt{3} p \sqrt{y}\right )$$

Пусть символ $n$ будет целым (I - это мнимая единица).

In [28]:
n=Symbol('n',integer=True)
exp(2*pi*I*n)
Out[28]:
$$1$$

Метод rewrite пытается переписать выражение в терминах заданной функции.

In [29]:
cos(x).rewrite(exp),exp(I*x).rewrite(cos)
Out[29]:
$$\left ( \frac{e^{i x}}{2} + \frac{1}{2} e^{- i x}, \quad i \sin{\left (x \right )} + \cos{\left (x \right )}\right )$$
In [30]:
asin(x).rewrite(log)
Out[30]:
$$- i \log{\left (i x + \sqrt{- x^{2} + 1} \right )}$$

Функция trigsimp пытается переписать тригонометрическое выражение в наиболее простом виде. В программах лучше использовать более специализированные функции.

In [31]:
trigsimp(2*sin(x)**2+3*cos(x)**2)
Out[31]:
$$\cos^{2}{\left (x \right )} + 2$$

Функция expand_trig разлагает синусы и косинусы сумм и кратных углов.

In [32]:
expand_trig(sin(x-y)),expand_trig(sin(2*x))
Out[32]:
$$\left ( \sin{\left (x \right )} \cos{\left (y \right )} - \sin{\left (y \right )} \cos{\left (x \right )}, \quad 2 \sin{\left (x \right )} \cos{\left (x \right )}\right )$$

Чаще нужно обратное преобразование - произведений и степеней синусов и косинусов в выражения, линейные по этим функциям. Например, пусть мы работаем с отрезком ряда Фурье.

In [33]:
a1,a2,b1,b2=symbols('a1 a2 b1 b2')
a=a1*cos(x)+a2*cos(2*x)+b1*sin(x)+b2*sin(2*x)
a
Out[33]:
$$a_{1} \cos{\left (x \right )} + a_{2} \cos{\left (2 x \right )} + b_{1} \sin{\left (x \right )} + b_{2} \sin{\left (2 x \right )}$$

Мы хотим возвести его в квадрат и опять получить отрезок ряда Фурье.

In [34]:
a=(a**2).rewrite(exp).expand().rewrite(cos).expand()
a
Out[34]:
$$\frac{a_{1}^{2}}{2} \cos{\left (2 x \right )} + \frac{a_{1}^{2}}{2} + a_{1} a_{2} \cos{\left (x \right )} + a_{1} a_{2} \cos{\left (3 x \right )} + a_{1} b_{1} \sin{\left (2 x \right )} + a_{1} b_{2} \sin{\left (x \right )} + a_{1} b_{2} \sin{\left (3 x \right )} + \frac{a_{2}^{2}}{2} \cos{\left (4 x \right )} + \frac{a_{2}^{2}}{2} - a_{2} b_{1} \sin{\left (x \right )} + a_{2} b_{1} \sin{\left (3 x \right )} + a_{2} b_{2} \sin{\left (4 x \right )} - \frac{b_{1}^{2}}{2} \cos{\left (2 x \right )} + \frac{b_{1}^{2}}{2} + b_{1} b_{2} \cos{\left (x \right )} - b_{1} b_{2} \cos{\left (3 x \right )} - \frac{b_{2}^{2}}{2} \cos{\left (4 x \right )} + \frac{b_{2}^{2}}{2}$$
In [35]:
a.collect([cos(x),cos(2*x),cos(3*x),sin(x),sin(2*x),sin(3*x)])
Out[35]:
$$\frac{a_{1}^{2}}{2} + a_{1} b_{1} \sin{\left (2 x \right )} + \frac{a_{2}^{2}}{2} \cos{\left (4 x \right )} + \frac{a_{2}^{2}}{2} + a_{2} b_{2} \sin{\left (4 x \right )} + \frac{b_{1}^{2}}{2} - \frac{b_{2}^{2}}{2} \cos{\left (4 x \right )} + \frac{b_{2}^{2}}{2} + \left(\frac{a_{1}^{2}}{2} - \frac{b_{1}^{2}}{2}\right) \cos{\left (2 x \right )} + \left(a_{1} a_{2} - b_{1} b_{2}\right) \cos{\left (3 x \right )} + \left(a_{1} a_{2} + b_{1} b_{2}\right) \cos{\left (x \right )} + \left(a_{1} b_{2} - a_{2} b_{1}\right) \sin{\left (x \right )} + \left(a_{1} b_{2} + a_{2} b_{1}\right) \sin{\left (3 x \right )}$$

Функция expand_log преобразует логарифмы произведений и степеней в суммы логарифмов (только для положительных величин); logcombine производит обратное преобразование.

In [36]:
a=expand_log(log(p*q**2))
a
Out[36]:
$$\log{\left (p \right )} + 2 \log{\left (q \right )}$$
In [37]:
logcombine(a)
Out[37]:
$$\log{\left (p q^{2} \right )}$$

Функция expand_power_exp переписывает степени, показатели которых - суммы, через произведения степеней.

In [38]:
expand_power_exp(x**(p+q))
Out[38]:
$$x^{p} x^{q}$$

Функция expand_power_base переписывает степени, основания которых - произведения, через произведения степеней.

In [39]:
expand_power_base((x*y)**n)
Out[39]:
$$x^{n} y^{n}$$

Функция powsimp выполняет обратные преобразования.

In [40]:
powsimp(exp(x)*exp(2*y)),powsimp(x**n*y**n)
Out[40]:
$$\left ( e^{x + 2 y}, \quad \left(x y\right)^{n}\right )$$

Можно вводить функции пользователя. Они могут иметь произвольное число аргументов.

In [41]:
f=Function('f')
f(x)+f(x,y)
Out[41]:
$$f{\left (x \right )} + f{\left (x,y \right )}$$

Структура выражений

Внутреннее представление выражения - это дерево. Функция srepr возвращает строку, представляющую его.

In [42]:
srepr(x+1)
Out[42]:
"Add(Symbol('x'), Integer(1))"
In [43]:
srepr(x-1)
Out[43]:
"Add(Symbol('x'), Integer(-1))"
In [44]:
srepr(x-y)
Out[44]:
"Add(Symbol('x'), Mul(Integer(-1), Symbol('y')))"
In [45]:
srepr(2*x*y/3)
Out[45]:
"Mul(Rational(2, 3), Symbol('x'), Symbol('y'))"
In [46]:
srepr(x/y)
Out[46]:
"Mul(Symbol('x'), Pow(Symbol('y'), Integer(-1)))"

Вместо бинарных операций +, *, ** и т.д. можно использовать функции Add, Mul, Pow и т.д.

In [47]:
Mul(x,Pow(y,-1))
Out[47]:
$$\frac{x}{y}$$
In [48]:
srepr(f(x,y))
Out[48]:
"Function('f')(Symbol('x'), Symbol('y'))"

Атрибут func - это функция верхнего уровня в выражении, а args - список её аргументов.

In [49]:
a=2*x*y**2
a.func
Out[49]:
sympy.core.mul.Mul
In [50]:
a.args
Out[50]:
$$\left ( 2, \quad x, \quad y^{2}\right )$$
In [51]:
for i in a.args:
    print(i)
2
x
y**2

Функция subs заменяет переменную на выражение.

In [52]:
a.subs(y,2)
Out[52]:
$$8 x$$

Она можнет заменить несколько переменных. Для этого ей передаётся список кортежей или словарь.

In [53]:
a.subs([(x,pi),(y,2)])
Out[53]:
$$8 \pi$$
In [54]:
a.subs({x:pi,y:2})
Out[54]:
$$8 \pi$$

Она может заменить не переменную, а подвыражение - функцию с аргументами.

In [55]:
a=f(x)+f(y)
a.subs(f(y),1)
Out[55]:
$$f{\left (x \right )} + 1$$
In [56]:
(2*x*y*z).subs(x*y,z)
Out[56]:
$$2 z^{2}$$
In [57]:
(x+x**2+x**3+x**4).subs(x**2,y)
Out[57]:
$$x^{3} + x + y^{2} + y$$

Подстановки производятся последовательно. В данном случае сначала $x$ заменился на $y$, получилось $y^3+y^2$; потом в этом результате $y$ заменилось на $x$.

In [58]:
a=x**2+y**3
a.subs([(x,y),(y,x)])
Out[58]:
$$x^{3} + x^{2}$$

Если написать эти подстановки в другом порядке, результат будет другим.

In [59]:
a.subs([(y,x),(x,y)])
Out[59]:
$$y^{3} + y^{2}$$

Но можно передать функции subs ключевой параметр simultaneous=True, тогда подстановки будут производиться одновременно. Таким образом можно, например, поменять местами $x$ и $y$.

In [60]:
a.subs([(x,y),(y,x)],simultaneous=True)
Out[60]:
$$x^{3} + y^{2}$$

Можно заменить функцию на другую функцию.

In [61]:
g=Function('g')
a=f(x)+f(y)
a.subs(f,g)
Out[61]:
$$g{\left (x \right )} + g{\left (y \right )}$$

Метод replace ищет подвыражения, соответствующие образцу (содержащему произвольные переменные), и заменяет их на заданное выражение (оно может содержать те же произвольные переменные).

In [62]:
a=Wild('a')
(f(x)+f(x+y)).replace(f(a),a**2)
Out[62]:
$$x^{2} + \left(x + y\right)^{2}$$
In [63]:
(f(x,x)+f(x,y)).replace(f(a,a),a**2)
Out[63]:
$$x^{2} + f{\left (x,y \right )}$$
In [64]:
a=x**2+y**2
a.replace(x,x+1)
Out[64]:
$$y^{2} + \left(x + 1\right)^{2}$$

Соответствовать образцу должно целое подвыражение, это не может быть часть сомножителей в произведении или меньшая степеть в большей.

In [65]:
a=2*x*y*z
a.replace(x*y,z)
Out[65]:
$$2 x y z$$
In [66]:
(x+x**2+x**3+x**4).replace(x**2,y)
Out[66]:
$$x^{4} + x^{3} + x + y$$

Решение уравнений

In [67]:
a,b,c,d,e,f=symbols('a b c d e f')

Уравнение записывается как функция Eq с двумя параметрами. Функция solve возврящает список решений.

In [68]:
solve(Eq(a*x,b),x)
Out[68]:
$$\left [ \frac{b}{a}\right ]$$

Впрочем, можно передать функции solve просто выражение. Подразумевается уравнение, что это выражение равно 0.

In [69]:
solve(a*x+b,x)
Out[69]:
$$\left [ - \frac{b}{a}\right ]$$

Квадратное уравнение имеет 2 решения.

In [70]:
solve(a*x**2+b*x+c,x)
Out[70]:
$$\left [ \frac{1}{2 a} \left(- b + \sqrt{- 4 a c + b^{2}}\right), \quad - \frac{1}{2 a} \left(b + \sqrt{- 4 a c + b^{2}}\right)\right ]$$

Система линейных уравнений.

In [71]:
solve([a*x+b*y-e,c*x+d*y-f],[x,y])
Out[71]:
$$\left \{ x : \frac{- b f + d e}{a d - b c}, \quad y : \frac{a f - c e}{a d - b c}\right \}$$

Функция roots возвращает корни многочлена с их множественностями.

In [72]:
roots(x**3-3*x+2,x)
Out[72]:
$$\left \{ -2 : 1, \quad 1 : 2\right \}$$

Функция solve_poly_system решает систему полиномиальных уравнений, строя их базис Грёбнера.

In [73]:
p1=x**2+y**2-1
p2=4*x*y-1
solve_poly_system([p1,p2],x,y)
Out[73]:
$$\left [ \left ( 4 \left(-1 - \sqrt{- \frac{\sqrt{3}}{4} + \frac{1}{2}}\right) \sqrt{- \frac{\sqrt{3}}{4} + \frac{1}{2}} \left(- \sqrt{- \frac{\sqrt{3}}{4} + \frac{1}{2}} + 1\right), \quad - \sqrt{- \frac{\sqrt{3}}{4} + \frac{1}{2}}\right ), \quad \left ( - 4 \left(-1 + \sqrt{- \frac{\sqrt{3}}{4} + \frac{1}{2}}\right) \sqrt{- \frac{\sqrt{3}}{4} + \frac{1}{2}} \left(\sqrt{- \frac{\sqrt{3}}{4} + \frac{1}{2}} + 1\right), \quad \sqrt{- \frac{\sqrt{3}}{4} + \frac{1}{2}}\right ), \quad \left ( 4 \left(-1 - \sqrt{\frac{\sqrt{3}}{4} + \frac{1}{2}}\right) \sqrt{\frac{\sqrt{3}}{4} + \frac{1}{2}} \left(- \sqrt{\frac{\sqrt{3}}{4} + \frac{1}{2}} + 1\right), \quad - \sqrt{\frac{\sqrt{3}}{4} + \frac{1}{2}}\right ), \quad \left ( - 4 \left(-1 + \sqrt{\frac{\sqrt{3}}{4} + \frac{1}{2}}\right) \sqrt{\frac{\sqrt{3}}{4} + \frac{1}{2}} \left(\sqrt{\frac{\sqrt{3}}{4} + \frac{1}{2}} + 1\right), \quad \sqrt{\frac{\sqrt{3}}{4} + \frac{1}{2}}\right )\right ]$$

Ряды

In [74]:
exp(x).series(x,0,5)
Out[74]:
$$1 + x + \frac{x^{2}}{2} + \frac{x^{3}}{6} + \frac{x^{4}}{24} + \mathcal{O}\left(x^{5}\right)$$

Ряд может начинаться с отрицательной степени.

In [75]:
cot(x).series(x,n=5)
Out[75]:
$$\frac{1}{x} - \frac{x}{3} - \frac{x^{3}}{45} + \mathcal{O}\left(x^{5}\right)$$

И даже идти по полуцелым степеням.

In [76]:
sqrt(x*(1-x)).series(x,n=5)
Out[76]:
$$\sqrt{x} - \frac{x^{\frac{3}{2}}}{2} - \frac{x^{\frac{5}{2}}}{8} - \frac{x^{\frac{7}{2}}}{16} - \frac{5 x^{\frac{9}{2}}}{128} + \mathcal{O}\left(x^{5}\right)$$
In [77]:
log(gamma(1+x)).series(x,n=6).rewrite(zeta)
Out[77]:
$$- \gamma x + \frac{\pi^{2} x^{2}}{12} - \frac{x^{3} \zeta\left(3\right)}{3} + \frac{\pi^{4} x^{4}}{360} - \frac{x^{5} \zeta\left(5\right)}{5} + \mathcal{O}\left(x^{6}\right)$$

Подготовим 3 ряда.

In [78]:
sinx=series(sin(x),x,0,8)
sinx
Out[78]:
$$x - \frac{x^{3}}{6} + \frac{x^{5}}{120} - \frac{x^{7}}{5040} + \mathcal{O}\left(x^{8}\right)$$
In [79]:
cosx=series(cos(x),x,n=8)
cosx
Out[79]:
$$1 - \frac{x^{2}}{2} + \frac{x^{4}}{24} - \frac{x^{6}}{720} + \mathcal{O}\left(x^{8}\right)$$
In [80]:
tanx=series(tan(x),x,n=8)
tanx
Out[80]:
$$x + \frac{x^{3}}{3} + \frac{2 x^{5}}{15} + \frac{17 x^{7}}{315} + \mathcal{O}\left(x^{8}\right)$$

Произведения и частные рядов не вычисляются автоматически, к ним надо применить функцию series.

In [81]:
series(tanx*cosx,n=8)
Out[81]:
$$x - \frac{x^{3}}{6} + \frac{x^{5}}{120} - \frac{x^{7}}{5040} + \mathcal{O}\left(x^{8}\right)$$
In [82]:
series(sinx/cosx,n=8)
Out[82]:
$$x + \frac{x^{3}}{3} + \frac{2 x^{5}}{15} + \frac{17 x^{7}}{315} + \mathcal{O}\left(x^{8}\right)$$

А этот ряд должен быть равен 1. Но поскольку sinx и cosx известны лишь с ограниченной точностью, мы получаем 1 с той же точностью.

In [83]:
series(sinx**2+cosx**2,n=8)
Out[83]:
$$1 + \mathcal{O}\left(x^{8}\right)$$

Здесь первые члены сократились, и ответ можно получить лишь с меньшей точностью.

In [84]:
series((1-cosx)/x**2,n=6)
Out[84]:
$$\frac{1}{2} - \frac{x^{2}}{24} + \frac{x^{4}}{720} + \mathcal{O}\left(x^{6}\right)$$

Ряды можно дифференцировать и интегрировать.

In [85]:
diff(sinx,x)
Out[85]:
$$1 - \frac{x^{2}}{2} + \frac{x^{4}}{24} - \frac{x^{6}}{720} + \mathcal{O}\left(x^{7}\right)$$
In [86]:
integrate(cosx,x)
Out[86]:
$$x - \frac{x^{3}}{6} + \frac{x^{5}}{120} - \frac{x^{7}}{5040} + \mathcal{O}\left(x^{9}\right)$$

Можно подставить ряд (если он начинается с малого члена) вместо переменной разложения в другой ряд. Вот ряды для $\sin(\tan(x))$ и $\tan(\sin(x))$.

In [87]:
st=series(sinx.subs(x,tanx),n=8)
st
Out[87]:
$$x + \frac{x^{3}}{6} - \frac{x^{5}}{40} - \frac{55 x^{7}}{1008} + \mathcal{O}\left(x^{8}\right)$$
In [88]:
ts=series(tanx.subs(x,sinx),n=8)
ts
Out[88]:
$$x + \frac{x^{3}}{6} - \frac{x^{5}}{40} - \frac{107 x^{7}}{5040} + \mathcal{O}\left(x^{8}\right)$$
In [89]:
series(ts-st,n=8)
Out[89]:
$$\frac{x^{7}}{30} + \mathcal{O}\left(x^{8}\right)$$

В ряд нельзя подставлять численное значение переменной разложения (а значит, нельзя и строить график). Для этого нужно сначала убрать $\mathcal{O}$ член, превратив отрезок ряда в многочлен.

In [90]:
sinx.removeO()
Out[90]:
$$- \frac{x^{7}}{5040} + \frac{x^{5}}{120} - \frac{x^{3}}{6} + x$$

Производные

In [91]:
a=x*sin(x+y)
diff(a,x)
Out[91]:
$$x \cos{\left (x + y \right )} + \sin{\left (x + y \right )}$$
In [92]:
diff(a,y)
Out[92]:
$$x \cos{\left (x + y \right )}$$

Вторая производная по $x$ и первая по $y$.

In [93]:
diff(a,x,2,y)
Out[93]:
$$- x \cos{\left (x + y \right )} + 2 \sin{\left (x + y \right )}$$

Можно дифференцировать выражения, содержащие неопределённые функции.

In [94]:
a=x*f(x**2)
b=diff(a,x)
b
Out[94]:
$$2 x^{2} \left. \frac{d}{d \xi_{1}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x^{2} }} + f{\left (x^{2} \right )}$$

Что это за зверь такой получился?

In [95]:
print(b)
2*x**2*Subs(Derivative(f(_xi_1), _xi_1), (_xi_1,), (x**2,)) + f(x**2)

Функция Derivative представляет невычисленную производную. Её можно вычислить методом doit.

In [96]:
a=Derivative(sin(x),x)
Eq(a,a.doit())
Out[96]:
$$\frac{d}{d x} \sin{\left (x \right )} = \cos{\left (x \right )}$$

Интегралы

In [97]:
integrate(1/(x*(x**2-2)**2),x)
Out[97]:
$$\frac{1}{4} \log{\left (x \right )} - \frac{1}{8} \log{\left (x^{2} - 2 \right )} - \frac{1}{4 x^{2} - 8}$$
In [98]:
integrate(1/(exp(x)+1),x)
Out[98]:
$$x - \log{\left (e^{x} + 1 \right )}$$
In [99]:
integrate(log(x),x)
Out[99]:
$$x \log{\left (x \right )} - x$$
In [100]:
integrate(x*sin(x),x)
Out[100]:
$$- x \cos{\left (x \right )} + \sin{\left (x \right )}$$
In [101]:
integrate(x*exp(-x**2),x)
Out[101]:
$$- \frac{e^{- x^{2}}}{2}$$
In [102]:
a=integrate(x**x,x)
a
Out[102]:
$$\int x^{x}\, dx$$

Получился невычисленный интеграл.

In [103]:
print(a)
Integral(x**x, x)
In [104]:
a=Integral(sin(x),x)
Eq(a,a.doit())
Out[104]:
$$\int \sin{\left (x \right )}\, dx = - \cos{\left (x \right )}$$

Определённые интегралы.

In [105]:
integrate(sin(x),(x,0,pi))
Out[105]:
$$2$$

oo - это $\infty$.

In [106]:
integrate(exp(-x**2),(x,0,oo))
Out[106]:
$$\frac{\sqrt{\pi}}{2}$$
In [107]:
integrate(log(x)/(1-x),(x,0,1))
Out[107]:
$$- \frac{\pi^{2}}{6}$$

Суммирование рядов

In [108]:
summation(1/n**2,(n,1,oo))
Out[108]:
$$\frac{\pi^{2}}{6}$$
In [109]:
summation((-1)**n/n**2,(n,1,oo))
Out[109]:
$$- \frac{\pi^{2}}{12}$$
In [110]:
summation(1/n**4,(n,1,oo))
Out[110]:
$$\frac{\pi^{4}}{90}$$

Невычисленная сумма обозначается Sum.

In [111]:
a=Sum(x**n/factorial(n),(n,0,oo))
Eq(a,a.doit())
Out[111]:
$$\sum_{n=0}^{\infty} \frac{x^{n}}{n!} = e^{x}$$

Пределы

In [112]:
limit((tan(sin(x))-sin(tan(x)))/x**7,x,0)
Out[112]:
$$\frac{1}{30}$$

Ну это простой предел, считается разложением числителя и знаменателя в ряды. А вот если в $x=0$ существенно особая точка, дело сложнее. Посчитаем односторонние пределы.

In [113]:
limit((tan(sin(x))-sin(tan(x)))/(x**7+exp(-1/x)),x,0,'+')
Out[113]:
$$\frac{1}{30}$$
In [114]:
limit((tan(sin(x))-sin(tan(x)))/(x**7+exp(-1/x)),x,0,'-')
Out[114]:
$$0$$

Дифференциальные уравнения

In [115]:
t=Symbol('t')
x=Function('x')
p=Function('p')

Первого порядка.

In [116]:
dsolve(diff(x(t),t)+x(t),x(t))
Out[116]:
$$x{\left (t \right )} = C_{1} e^{- t}$$

Второго порядка.

In [117]:
dsolve(diff(x(t),t,2)+x(t),x(t))
Out[117]:
$$x{\left (t \right )} = C_{1} \sin{\left (t \right )} + C_{2} \cos{\left (t \right )}$$

Система уравнений первого порядка.

In [118]:
dsolve((diff(x(t),t)-p(t),diff(p(t),t)+x(t)))
Out[118]:
$$\left [ x{\left (t \right )} = C_{1} \sin{\left (t \right )} + C_{2} \cos{\left (t \right )}, \quad p{\left (t \right )} = C_{1} \cos{\left (t \right )} - C_{2} \sin{\left (t \right )}\right ]$$

Линейная алгебра

In [119]:
a,b,c,d,e,f=symbols('a b c d e f')

Матрицу можно построить из списка списков.

In [120]:
M=Matrix([[a,b,c],[d,e,f]])
M
Out[120]:
$$\left[\begin{matrix}a & b & c\\d & e & f\end{matrix}\right]$$
In [121]:
M.shape
Out[121]:
$$\left ( 2, \quad 3\right )$$

Матрица-строка.

In [122]:
Matrix([[1,2,3]])
Out[122]:
$$\left[\begin{matrix}1 & 2 & 3\end{matrix}\right]$$

Матрица-столбец.

In [123]:
Matrix([1,2,3])
Out[123]:
$$\left[\begin{matrix}1\\2\\3\end{matrix}\right]$$

Можно построить матрицу из функции.

In [124]:
def g(i,j):
    return Rational(1,i+j+1)
Matrix(3,3,g)
Out[124]:
$$\left[\begin{matrix}1 & \frac{1}{2} & \frac{1}{3}\\\frac{1}{2} & \frac{1}{3} & \frac{1}{4}\\\frac{1}{3} & \frac{1}{4} & \frac{1}{5}\end{matrix}\right]$$

Или из неопределённой функции.

In [125]:
g=Function('g')
M=Matrix(3,3,g)
M
Out[125]:
$$\left[\begin{matrix}g{\left (0,0 \right )} & g{\left (0,1 \right )} & g{\left (0,2 \right )}\\g{\left (1,0 \right )} & g{\left (1,1 \right )} & g{\left (1,2 \right )}\\g{\left (2,0 \right )} & g{\left (2,1 \right )} & g{\left (2,2 \right )}\end{matrix}\right]$$
In [126]:
M[1,2]
Out[126]:
$$g{\left (1,2 \right )}$$
In [127]:
M[1,2]=0
M
Out[127]:
$$\left[\begin{matrix}g{\left (0,0 \right )} & g{\left (0,1 \right )} & g{\left (0,2 \right )}\\g{\left (1,0 \right )} & g{\left (1,1 \right )} & 0\\g{\left (2,0 \right )} & g{\left (2,1 \right )} & g{\left (2,2 \right )}\end{matrix}\right]$$
In [128]:
M[2,:]
Out[128]:
$$\left[\begin{matrix}g{\left (2,0 \right )} & g{\left (2,1 \right )} & g{\left (2,2 \right )}\end{matrix}\right]$$
In [129]:
M[:,1]
Out[129]:
$$\left[\begin{matrix}g{\left (0,1 \right )}\\g{\left (1,1 \right )}\\g{\left (2,1 \right )}\end{matrix}\right]$$
In [130]:
M[0:2,1:3]
Out[130]:
$$\left[\begin{matrix}g{\left (0,1 \right )} & g{\left (0,2 \right )}\\g{\left (1,1 \right )} & 0\end{matrix}\right]$$

Единичная матрица.

In [131]:
eye(3)
Out[131]:
$$\left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]$$

Матрица из нулей.

In [132]:
zeros(3)
Out[132]:
$$\left[\begin{matrix}0 & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right]$$
In [133]:
zeros(2,3)
Out[133]:
$$\left[\begin{matrix}0 & 0 & 0\\0 & 0 & 0\end{matrix}\right]$$

Диагональная матрица.

In [134]:
diag(1,2,3)
Out[134]:
$$\left[\begin{matrix}1 & 0 & 0\\0 & 2 & 0\\0 & 0 & 3\end{matrix}\right]$$
In [135]:
M=Matrix([[a,1],[0,a]])
diag(1,M,2)
Out[135]:
$$\left[\begin{matrix}1 & 0 & 0 & 0\\0 & a & 1 & 0\\0 & 0 & a & 0\\0 & 0 & 0 & 2\end{matrix}\right]$$

Операции с матрицами.

In [136]:
A=Matrix([[a,b],[c,d]])
B=Matrix([[1,2],[3,4]])
A+B
Out[136]:
$$\left[\begin{matrix}a + 1 & b + 2\\c + 3 & d + 4\end{matrix}\right]$$
In [137]:
A*B,B*A
Out[137]:
$$\left ( \left[\begin{matrix}a + 3 b & 2 a + 4 b\\c + 3 d & 2 c + 4 d\end{matrix}\right], \quad \left[\begin{matrix}a + 2 c & b + 2 d\\3 a + 4 c & 3 b + 4 d\end{matrix}\right]\right )$$
In [138]:
A*B-B*A
Out[138]:
$$\left[\begin{matrix}3 b - 2 c & 2 a + 3 b - 2 d\\- 3 a - 3 c + 3 d & - 3 b + 2 c\end{matrix}\right]$$
In [139]:
simplify(A**(-1))
Out[139]:
$$\left[\begin{matrix}\frac{d}{a d - b c} & - \frac{b}{a d - b c}\\- \frac{c}{a d - b c} & \frac{a}{a d - b c}\end{matrix}\right]$$
In [140]:
det(A)
Out[140]:
$$a d - b c$$

Собственные значения и векторы

In [141]:
x=Symbol('x',real=True)
In [142]:
M=Matrix([[(1-x)**3*(3+x),4*x*(1-x**2),-2*(1-x**2)*(3-x)],
          [4*x*(1-x**2),-(1+x)**3*(3-x),2*(1-x**2)*(3+x)],
          [-2*(1-x**2)*(3-x),2*(1-x**2)*(3+x),16*x]])
M
Out[142]:
$$\left[\begin{matrix}\left(- x + 1\right)^{3} \left(x + 3\right) & 4 x \left(- x^{2} + 1\right) & \left(- x + 3\right) \left(2 x^{2} - 2\right)\\4 x \left(- x^{2} + 1\right) & - \left(- x + 3\right) \left(x + 1\right)^{3} & \left(x + 3\right) \left(- 2 x^{2} + 2\right)\\\left(- x + 3\right) \left(2 x^{2} - 2\right) & \left(x + 3\right) \left(- 2 x^{2} + 2\right) & 16 x\end{matrix}\right]$$
In [143]:
det(M)
Out[143]:
$$0$$

Значит, у этой матрицы есть нулевое подпространство (она обращает векторы из этого подпространства в 0). Базис этого подпространства.

In [144]:
v=M.nullspace()
len(v)
Out[144]:
$$1$$

Оно одномерно.

In [145]:
v=simplify(v[0])
v
Out[145]:
$$\left[\begin{matrix}- \frac{2}{x - 1}\\\frac{2}{x + 1}\\1\end{matrix}\right]$$

Проверим.

In [146]:
simplify(M*v)
Out[146]:
$$\left[\begin{matrix}0\\0\\0\end{matrix}\right]$$

Собственные значения и их кратности.

In [147]:
M.eigenvals()
Out[147]:
$$\left \{ 0 : 1, \quad - \left(x^{2} + 3\right)^{2} : 1, \quad \left(x^{2} + 3\right)^{2} : 1\right \}$$

Если нужны не только собственные значения, но и собственные векторы, то нужно использовать метод eigenvects. Он возвращает список кортежей. В каждом из них нулевой элемент - собственное значение, первый - его кратность, и последний - список собственных векторов, образующих базис (их столько, какова кратность).

In [148]:
v=M.eigenvects()
len(v)
Out[148]:
$$3$$
In [149]:
for i in range(len(v)):
    v[i][2][0]=simplify(v[i][2][0])
v
Out[149]:
$$\left [ \left ( 0, \quad 1, \quad \left [ \left[\begin{matrix}- \frac{2}{x - 1}\\\frac{2}{x + 1}\\1\end{matrix}\right]\right ]\right ), \quad \left ( - \left(x^{2} + 3\right)^{2}, \quad 1, \quad \left [ \left[\begin{matrix}\frac{x}{2} + \frac{1}{2}\\\frac{x + 1}{x - 1}\\1\end{matrix}\right]\right ]\right ), \quad \left ( \left(x^{2} + 3\right)^{2}, \quad 1, \quad \left [ \left[\begin{matrix}\frac{x - 1}{x + 1}\\- \frac{x}{2} + \frac{1}{2}\\1\end{matrix}\right]\right ]\right )\right ]$$

Проверим.

In [150]:
for i in range(len(v)):
    z=M*v[i][2][0]-v[i][0]*v[i][2][0]
    pprint(simplify(z))
⎡0⎤
⎢ ⎥
⎢0⎥
⎢ ⎥
⎣0⎦
⎡0⎤
⎢ ⎥
⎢0⎥
⎢ ⎥
⎣0⎦
⎡0⎤
⎢ ⎥
⎢0⎥
⎢ ⎥
⎣0⎦

Жорданова нормальная форма

In [151]:
M=Matrix([[Rational(13,9),-Rational(2,9),Rational(1,3),Rational(4,9),Rational(2,3)],
          [-Rational(2,9),Rational(10,9),Rational(2,15),-Rational(2,9),-Rational(11,15)],
          [Rational(1,5),-Rational(2,5),Rational(41,25),-Rational(2,5),Rational(12,25)],
          [Rational(4,9),-Rational(2,9),Rational(14,15),Rational(13,9),-Rational(2,15)],
          [-Rational(4,15),Rational(8,15),Rational(12,25),Rational(8,15),Rational(34,25)]])
M
Out[151]:
$$\left[\begin{matrix}\frac{13}{9} & - \frac{2}{9} & \frac{1}{3} & \frac{4}{9} & \frac{2}{3}\\- \frac{2}{9} & \frac{10}{9} & \frac{2}{15} & - \frac{2}{9} & - \frac{11}{15}\\\frac{1}{5} & - \frac{2}{5} & \frac{41}{25} & - \frac{2}{5} & \frac{12}{25}\\\frac{4}{9} & - \frac{2}{9} & \frac{14}{15} & \frac{13}{9} & - \frac{2}{15}\\- \frac{4}{15} & \frac{8}{15} & \frac{12}{25} & \frac{8}{15} & \frac{34}{25}\end{matrix}\right]$$

Метод M.jordan_form() возвращает пару матриц, матрицу преобразования $P$ и собственно жорданову форму $J$: $M = P J P^{-1}$.

In [152]:
P,J=M.jordan_form()
J
Out[152]:
$$\left[\begin{matrix}1 & 0 & 0 & 0 & 0\\0 & 2 & 1 & 0 & 0\\0 & 0 & 2 & 0 & 0\\0 & 0 & 0 & 1 - i & 0\\0 & 0 & 0 & 0 & 1 + i\end{matrix}\right]$$
In [153]:
P=simplify(P)
P
Out[153]:
$$\left[\begin{matrix}-2 & \frac{10}{9} & 0 & \frac{5 i}{12} & - \frac{5 i}{12}\\-2 & - \frac{5}{9} & 0 & - \frac{5 i}{6} & \frac{5 i}{6}\\0 & 0 & \frac{4}{3} & - \frac{3}{4} & - \frac{3}{4}\\1 & \frac{10}{9} & 0 & - \frac{5 i}{6} & \frac{5 i}{6}\\0 & 0 & 1 & 1 & 1\end{matrix}\right]$$

Проверим.

In [154]:
Z=P*J*P**(-1)-M
simplify(Z)
Out[154]:
$$\left[\begin{matrix}0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\end{matrix}\right]$$

Графики

SymPy использует matplotlib. Однако он распределяет точки по $x$ адаптивно, а не равномерно.

In [155]:
%matplotlib inline

Одна функция.

In [156]:
plot(sin(x)/x,(x,-10,10))
Out[156]:
<sympy.plotting.plot.Plot at 0xaf83218c>

Несколько функций.

In [157]:
plot(x,x**2,x**3,(x,0,2))
Out[157]:
<sympy.plotting.plot.Plot at 0xb14ab02c>

Другие функции надо импортировать из пакета sympy.plotting.

In [158]:
from sympy.plotting import (plot_parametric,plot_implicit,
                            plot3d,plot3d_parametric_line,
                            plot3d_parametric_surface)

Параметрический график - фигура Лиссажу.

In [159]:
t=Symbol('t')
plot_parametric(sin(2*t),cos(3*t),(t,0,2*pi),
                title='Lissajous',xlabel='x',ylabel='y')
Out[159]:
<sympy.plotting.plot.Plot at 0xaf6a188c>

Неявный график - окружность.

In [160]:
plot_implicit(x**2+y**2-1,(x,-1,1),(y,-1,1))
Out[160]:
<sympy.plotting.plot.Plot at 0xaf71c9cc>

Поверхность. Если она строится не inline, а в отдельном окне, то её можно вертеть мышкой.

In [161]:
plot3d(x*y,(x,-2,2),(y,-2,2))
Out[161]:
<sympy.plotting.plot.Plot at 0xaf6bea8c>

Несколько поверхностей.

In [162]:
plot3d(x**2+y**2,x*y,(x,-2,2),(y,-2,2))
Out[162]:
<sympy.plotting.plot.Plot at 0xaf64c36c>

Параметрическая пространственная линия - спираль.

In [163]:
a=0.1
plot3d_parametric_line(cos(t),sin(t),a*t,(t,0,4*pi))
Out[163]:
<sympy.plotting.plot.Plot at 0xaf714c0c>

Параметрическая поверхность - тор.

In [164]:
u,v=symbols('u v')
a=0.3
plot3d_parametric_surface((1+a*cos(u))*cos(v),
                          (1+a*cos(u))*sin(v),a*sin(u),
                          (u,0,2*pi),(v,0,2*pi))
Out[164]:
<sympy.plotting.plot.Plot at 0xaf646cac>