SymPy
- это пакет для символьных вычислений на питоне, подобный системе Mathematica. Он работает с выражениями, содержащими символы.
from sympy import *
init_printing()
Основными кирпичиками, из которых строятся выражения, являются символы. Символ имеет имя, которое используется при печати выражений. Объекты класса Symbol
нужно создавать и присваивать переменным питона, чтобы их можно было использовать. В принципе, имя символа и имя переменной, которой мы присваиваем этот символ - две независимые вещи, и можно написать abc=Symbol('xyz')
. Но тогда при вводе программы Вы будете использовать abc
, а при печати результатов SymPy
будет использовать xyz
, что приведёт к ненужной путанице. Поэтому лучше, чтобы имя символа совпадало с именем переменной питона, которой он присваивается.
В языках, специально предназначенных для символьных вычислений, таких, как Mathematica, если Вы используете переменную, которой ничего не было присвоено, то она автоматически воспринимается как символ с тем же именем. Питон не был изначально предназначен для символьных вычислений. Если Вы используете переменную, которой ничего не было присвоено, Вы получите сообщение об ошибке. Объекты типа Symbol
нужно создавать явно.
x=Symbol('x')
a=x**2-1
a
Можно определить несколько символов одновременно. Строка разбивается на имена по пробелам.
y,z=symbols('y z')
Подставим вместо $x$ выражение $y+1$.
a.subs(x,y+1)
SymPy
не раскрывает скобки автоматически. Для этого используется функция expand
.
a=(x+y-z)**6
a
a=expand(a)
a
Степень многочлена $a$ по $x$.
degree(a,x)
Соберём вместе члены с определёнными степенями $x$.
collect(a,x)
Многочлен с целыми коэффициентами можно записать в виде произведения таких многочленов (причём каждый сомножитель уже невозможно расфакторизовать дальше, оставаясь в рамках многочленов с целыми коэффициентами). Существуют эффективные алгоритмы для решения этой задачи.
a=factor(a)
a
SymPy
не сокращает отношения многочленов на их наибольший общий делитель автоматически. Для этого используется функция cancel
.
a=(x**3-y**3)/(x**2-y**2)
a
cancel(a)
SymPy
не приводит суммы рациональных выражений к общему знаменателю автоматически. Для этого используется функция together
.
a=y/(x-y)+x/(x+y)
a
together(a)
Функция simplify
пытается переписать выражение в наиболее простом виде. Это понятие не имеет чёткого определения (в разных ситуациях наиболее простыми могут считаться разные формы выражения), и не существует алгоритма такого упрощения. Функция symplify
работает эвристически, и невозможно заранее предугадать, какие упрощения она попытается сделать. Поэтому её удобно использовать в интерактивных сессиях, чтобы посмотреть, удастся ли ей записать выражение в каком-нибудь разумном виде, но нежелательно использовать в программах. В них лучше применять более специализированные функции, которые выполняют одно определённое преобразование выражения.
simplify(a)
Разложение на элементарные дроби по отношению к $x$.
apart(a,x)
Подставим конкретные численные значения вместо переменных $x$ и $y$.
a=a.subs({x:1,y:2})
a
А сколько это будет численно?
a.n()
SymPy
автоматически применяет упрощения элементарных функция (которые справедливы во всех случаях).
sin(-x)
cos(pi/4),tan(5*pi/6)
SymPy
может работать с числами с плавающей точкой, имеющими сколь угодно большую точность. Вот $\pi$ с 100 значащими цифрами.
pi.n(100)
E
- это основание натуральных логарифмов.
log(1),log(E)
exp(log(x)),log(exp(x))
А почему не $x$? Попробуйте подставить $x=2\pi i$.
sqrt(0)
sqrt(x)**4,sqrt(x**4)
А почему не $x^2$? Попробуйте подставить $x=i$.
Символы могут иметь некоторые свойства. Например, они могут быть положительными. Тогда SymPy
может сильнее упростить квадратные корни.
p,q=symbols('p q',positive=True)
sqrt(p**2)
sqrt(12*x**2*y),sqrt(12*p**2*y)
Пусть символ $n$ будет целым (I
- это мнимая единица).
n=Symbol('n',integer=True)
exp(2*pi*I*n)
Метод rewrite
пытается переписать выражение в терминах заданной функции.
cos(x).rewrite(exp),exp(I*x).rewrite(cos)
asin(x).rewrite(log)
Функция trigsimp
пытается переписать тригонометрическое выражение в наиболее простом виде. В программах лучше использовать более специализированные функции.
trigsimp(2*sin(x)**2+3*cos(x)**2)
Функция expand_trig
разлагает синусы и косинусы сумм и кратных углов.
expand_trig(sin(x-y)),expand_trig(sin(2*x))
Чаще нужно обратное преобразование - произведений и степеней синусов и косинусов в выражения, линейные по этим функциям. Например, пусть мы работаем с отрезком ряда Фурье.
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
Мы хотим возвести его в квадрат и опять получить отрезок ряда Фурье.
a=(a**2).rewrite(exp).expand().rewrite(cos).expand()
a
a.collect([cos(x),cos(2*x),cos(3*x),sin(x),sin(2*x),sin(3*x)])
Функция expand_log
преобразует логарифмы произведений и степеней в суммы логарифмов (только для положительных величин); logcombine
производит обратное преобразование.
a=expand_log(log(p*q**2))
a
logcombine(a)
Функция expand_power_exp
переписывает степени, показатели которых - суммы, через произведения степеней.
expand_power_exp(x**(p+q))
Функция expand_power_base
переписывает степени, основания которых - произведения, через произведения степеней.
expand_power_base((x*y)**n)
Функция powsimp
выполняет обратные преобразования.
powsimp(exp(x)*exp(2*y)),powsimp(x**n*y**n)
Можно вводить функции пользователя. Они могут иметь произвольное число аргументов.
f=Function('f')
f(x)+f(x,y)
Внутреннее представление выражения - это дерево. Функция srepr
возвращает строку, представляющую его.
srepr(x+1)
srepr(x-1)
srepr(x-y)
srepr(2*x*y/3)
srepr(x/y)
Вместо бинарных операций +
, *
, **
и т.д. можно использовать функции Add
, Mul
, Pow
и т.д.
Mul(x,Pow(y,-1))
srepr(f(x,y))
Атрибут func
- это функция верхнего уровня в выражении, а args
- список её аргументов.
a=2*x*y**2
a.func
a.args
for i in a.args:
print(i)
Функция subs
заменяет переменную на выражение.
a.subs(y,2)
Она можнет заменить несколько переменных. Для этого ей передаётся список кортежей или словарь.
a.subs([(x,pi),(y,2)])
a.subs({x:pi,y:2})
Она может заменить не переменную, а подвыражение - функцию с аргументами.
a=f(x)+f(y)
a.subs(f(y),1)
(2*x*y*z).subs(x*y,z)
(x+x**2+x**3+x**4).subs(x**2,y)
Подстановки производятся последовательно. В данном случае сначала $x$ заменился на $y$, получилось $y^3+y^2$; потом в этом результате $y$ заменилось на $x$.
a=x**2+y**3
a.subs([(x,y),(y,x)])
Если написать эти подстановки в другом порядке, результат будет другим.
a.subs([(y,x),(x,y)])
Но можно передать функции subs
ключевой параметр simultaneous=True
, тогда подстановки будут производиться одновременно. Таким образом можно, например, поменять местами $x$ и $y$.
a.subs([(x,y),(y,x)],simultaneous=True)
Можно заменить функцию на другую функцию.
g=Function('g')
a=f(x)+f(y)
a.subs(f,g)
Метод replace
ищет подвыражения, соответствующие образцу (содержащему произвольные переменные), и заменяет их на заданное выражение (оно может содержать те же произвольные переменные).
a=Wild('a')
(f(x)+f(x+y)).replace(f(a),a**2)
(f(x,x)+f(x,y)).replace(f(a,a),a**2)
a=x**2+y**2
a.replace(x,x+1)
Соответствовать образцу должно целое подвыражение, это не может быть часть сомножителей в произведении или меньшая степеть в большей.
a=2*x*y*z
a.replace(x*y,z)
(x+x**2+x**3+x**4).replace(x**2,y)
a,b,c,d,e,f=symbols('a b c d e f')
Уравнение записывается как функция Eq
с двумя параметрами. Функция solve
возврящает список решений.
solve(Eq(a*x,b),x)
Впрочем, можно передать функции solve
просто выражение. Подразумевается уравнение, что это выражение равно 0.
solve(a*x+b,x)
Квадратное уравнение имеет 2 решения.
solve(a*x**2+b*x+c,x)
Система линейных уравнений.
solve([a*x+b*y-e,c*x+d*y-f],[x,y])
Функция roots
возвращает корни многочлена с их множественностями.
roots(x**3-3*x+2,x)
Функция solve_poly_system
решает систему полиномиальных уравнений, строя их базис Грёбнера.
p1=x**2+y**2-1
p2=4*x*y-1
solve_poly_system([p1,p2],x,y)
exp(x).series(x,0,5)
Ряд может начинаться с отрицательной степени.
cot(x).series(x,n=5)
И даже идти по полуцелым степеням.
sqrt(x*(1-x)).series(x,n=5)
log(gamma(1+x)).series(x,n=6).rewrite(zeta)
Подготовим 3 ряда.
sinx=series(sin(x),x,0,8)
sinx
cosx=series(cos(x),x,n=8)
cosx
tanx=series(tan(x),x,n=8)
tanx
Произведения и частные рядов не вычисляются автоматически, к ним надо применить функцию series
.
series(tanx*cosx,n=8)
series(sinx/cosx,n=8)
А этот ряд должен быть равен 1. Но поскольку sinx
и cosx
известны лишь с ограниченной точностью, мы получаем 1 с той же точностью.
series(sinx**2+cosx**2,n=8)
Здесь первые члены сократились, и ответ можно получить лишь с меньшей точностью.
series((1-cosx)/x**2,n=6)
Ряды можно дифференцировать и интегрировать.
diff(sinx,x)
integrate(cosx,x)
Можно подставить ряд (если он начинается с малого члена) вместо переменной разложения в другой ряд. Вот ряды для $\sin(\tan(x))$ и $\tan(\sin(x))$.
st=series(sinx.subs(x,tanx),n=8)
st
ts=series(tanx.subs(x,sinx),n=8)
ts
series(ts-st,n=8)
В ряд нельзя подставлять численное значение переменной разложения (а значит, нельзя и строить график). Для этого нужно сначала убрать $\mathcal{O}$ член, превратив отрезок ряда в многочлен.
sinx.removeO()
a=x*sin(x+y)
diff(a,x)
diff(a,y)
Вторая производная по $x$ и первая по $y$.
diff(a,x,2,y)
Можно дифференцировать выражения, содержащие неопределённые функции.
a=x*f(x**2)
b=diff(a,x)
b
Что это за зверь такой получился?
print(b)
Функция Derivative
представляет невычисленную производную. Её можно вычислить методом doit
.
a=Derivative(sin(x),x)
Eq(a,a.doit())
integrate(1/(x*(x**2-2)**2),x)
integrate(1/(exp(x)+1),x)
integrate(log(x),x)
integrate(x*sin(x),x)
integrate(x*exp(-x**2),x)
a=integrate(x**x,x)
a
Получился невычисленный интеграл.
print(a)
a=Integral(sin(x),x)
Eq(a,a.doit())
Определённые интегралы.
integrate(sin(x),(x,0,pi))
oo
- это $\infty$.
integrate(exp(-x**2),(x,0,oo))
integrate(log(x)/(1-x),(x,0,1))
summation(1/n**2,(n,1,oo))
summation((-1)**n/n**2,(n,1,oo))
summation(1/n**4,(n,1,oo))
Невычисленная сумма обозначается Sum
.
a=Sum(x**n/factorial(n),(n,0,oo))
Eq(a,a.doit())
limit((tan(sin(x))-sin(tan(x)))/x**7,x,0)
Ну это простой предел, считается разложением числителя и знаменателя в ряды. А вот если в $x=0$ существенно особая точка, дело сложнее. Посчитаем односторонние пределы.
limit((tan(sin(x))-sin(tan(x)))/(x**7+exp(-1/x)),x,0,'+')
limit((tan(sin(x))-sin(tan(x)))/(x**7+exp(-1/x)),x,0,'-')
t=Symbol('t')
x=Function('x')
p=Function('p')
Первого порядка.
dsolve(diff(x(t),t)+x(t),x(t))
Второго порядка.
dsolve(diff(x(t),t,2)+x(t),x(t))
Система уравнений первого порядка.
dsolve((diff(x(t),t)-p(t),diff(p(t),t)+x(t)))
a,b,c,d,e,f=symbols('a b c d e f')
Матрицу можно построить из списка списков.
M=Matrix([[a,b,c],[d,e,f]])
M
M.shape
Матрица-строка.
Matrix([[1,2,3]])
Матрица-столбец.
Matrix([1,2,3])
Можно построить матрицу из функции.
def g(i,j):
return Rational(1,i+j+1)
Matrix(3,3,g)
Или из неопределённой функции.
g=Function('g')
M=Matrix(3,3,g)
M
M[1,2]
M[1,2]=0
M
M[2,:]
M[:,1]
M[0:2,1:3]
Единичная матрица.
eye(3)
Матрица из нулей.
zeros(3)
zeros(2,3)
Диагональная матрица.
diag(1,2,3)
M=Matrix([[a,1],[0,a]])
diag(1,M,2)
Операции с матрицами.
A=Matrix([[a,b],[c,d]])
B=Matrix([[1,2],[3,4]])
A+B
A*B,B*A
A*B-B*A
simplify(A**(-1))
det(A)
x=Symbol('x',real=True)
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
det(M)
Значит, у этой матрицы есть нулевое подпространство (она обращает векторы из этого подпространства в 0). Базис этого подпространства.
v=M.nullspace()
len(v)
Оно одномерно.
v=simplify(v[0])
v
Проверим.
simplify(M*v)
Собственные значения и их кратности.
M.eigenvals()
Если нужны не только собственные значения, но и собственные векторы, то нужно использовать метод eigenvects
. Он возвращает список кортежей. В каждом из них нулевой элемент - собственное значение, первый - его кратность, и последний - список собственных векторов, образующих базис (их столько, какова кратность).
v=M.eigenvects()
len(v)
for i in range(len(v)):
v[i][2][0]=simplify(v[i][2][0])
v
Проверим.
for i in range(len(v)):
z=M*v[i][2][0]-v[i][0]*v[i][2][0]
pprint(simplify(z))
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
Метод M.jordan_form()
возвращает пару матриц, матрицу преобразования $P$ и собственно жорданову форму $J$: $M = P J P^{-1}$.
P,J=M.jordan_form()
J
P=simplify(P)
P
Проверим.
Z=P*J*P**(-1)-M
simplify(Z)
SymPy
использует matplotlib
. Однако он распределяет точки по $x$ адаптивно, а не равномерно.
%matplotlib inline
Одна функция.
plot(sin(x)/x,(x,-10,10))
Несколько функций.
plot(x,x**2,x**3,(x,0,2))
Другие функции надо импортировать из пакета sympy.plotting
.
from sympy.plotting import (plot_parametric,plot_implicit,
plot3d,plot3d_parametric_line,
plot3d_parametric_surface)
Параметрический график - фигура Лиссажу.
t=Symbol('t')
plot_parametric(sin(2*t),cos(3*t),(t,0,2*pi),
title='Lissajous',xlabel='x',ylabel='y')
Неявный график - окружность.
plot_implicit(x**2+y**2-1,(x,-1,1),(y,-1,1))
Поверхность. Если она строится не inline
, а в отдельном окне, то её можно вертеть мышкой.
plot3d(x*y,(x,-2,2),(y,-2,2))
Несколько поверхностей.
plot3d(x**2+y**2,x*y,(x,-2,2),(y,-2,2))
Параметрическая пространственная линия - спираль.
a=0.1
plot3d_parametric_line(cos(t),sin(t),a*t,(t,0,4*pi))
Параметрическая поверхность - тор.
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))