作者 : Fabian Pedregosa
目的
为什么是SymPy? SymPy是符号数学的Python库。它的目的是成为Mathematica或Maple等系统的替代品,同时让代码尽可能简单并且可扩展。SymPy完全是用Python写的,并不需要外部的库。
Sympy文档及库安装见http://www.sympy.org/
章节内容
SymPy定义了三种数字类型:实数、有理数和整数。
有理数类将有理数表征为两个整数对: 分子和分母,因此Rational(1,2)代表1/2, Rational(5,2)代表5/2等等:
from sympy import *
a = Rational(1,2)
a
1/2
a*2
1
SymPy在底层使用mpmath, 这使它可以用任意精度的算术进行计算。这样,一些特殊的常数,比如e, pi, oo (无限), 可以被作为符号处理并且可以以任意精度来评估:
pi**2
pi**2
pi.evalf()
3.14159265358979
(pi + exp(1)).evalf()
5.85987448204884
如你所见,将表达式评估为浮点数。
也有一个类代表数学的无限, 称为 oo:
oo > 99999
True
oo + 1
oo
from sympy import *
x = Symbol('x')
y = Symbol('y')
然后你可以计算他们:
x + y + x - y
2*x
(x + y)**2
(x + y)**2
符号可以使用一些Python操作符操作: +, -, *, ** (算术), &, |, ~ , >>, << (布尔逻辑).
打印 这里我们使用下列设置打印
sympy.init_printing(use_unicode=False, wrap_line=True)
expand((x + y)**3)
x**3 + 3*x**2*y + 3*x*y**2 + y**3
3*x*y**2 + 3*y*x**2 + x**3 + y**3
x**3 + 3*x**2*y + 3*x*y**2 + y**3
可以通过关键词的形式使用更多的选项:
expand(x + y, complex=True)
re(x) + re(y) + I*im(x) + I*im(y)
I*im(x) + I*im(y) + re(x) + re(y)
re(x) + re(y) + I*im(x) + I*im(y)
expand(cos(x + y), trig=True)
-sin(x)*sin(y) + cos(x)*cos(y)
cos(x)*cos(y) - sin(x)*sin(y)
-sin(x)*sin(y) + cos(x)*cos(y)
如果可以将表达式转化为更简单的形式,可以使用化简:
simplify((x + x*y) / x)
y + 1
limit(sin(x)/x, x, 0)
1
你也可以计算一下在无限时候的极限:
limit(x, x, oo)
oo
limit(1/x, x, oo)
0
limit(x**x, x, 0)
1
你可以使用diff(func, var)
微分任何SymPy表达式。例如:
diff(sin(x), x)
cos(x)
diff(sin(2*x), x)
2*cos(2*x)
diff(tan(x), x)
tan(x)**2 + 1
你可以用下列方法检查是否正确:
limit((tan(x+y) - tan(x))/y, y, 0)
tan(x)**2 + 1
可以用diff(func, var, n)
方法来计算更高的导数:
diff(sin(2*x), x, 1)
2*cos(2*x)
diff(sin(2*x), x, 2)
-4*sin(2*x)
diff(sin(2*x), x, 3)
-8*cos(2*x)
SymPy也知道如何计算一个表达式在一个点的Taylor序列。使用series(expr, var)
:
series(cos(x), x)
1 - x**2/2 + x**4/24 + O(x**6)
series(1/cos(x), x)
1 + x**2/2 + 5*x**4/24 + O(x**6)
练习
计算$\lim_{x\rightarrow 0} \sin(x)/x$
计算log(x)
对于x的导数。
SymPy支持超验基础和特殊函数的无限和有限积分,通过integrate()
功能, 使用了强大的扩展的Risch-Norman算法和启发式和模式匹配。你可以积分基本函数:
integrate(6*x**5, x)
x**6
integrate(sin(x), x)
-cos(x)
integrate(log(x), x)
x*log(x) - x
integrate(2*x + sinh(x), x)
x**2 + cosh(x)
也可以很简单的处理特殊函数:
integrate(exp(-x**2)*erf(x), x)
sqrt(pi)*erf(x)**2/4
也可以计算一下有限积分:
integrate(x**3, (x, -1, 1))
0
integrate(sin(x), (x, 0, pi/2))
1
integrate(cos(x), (x, -pi/2, pi/2))
2
不标准积分也支持:
integrate(exp(-x), (x, 0, oo))
1
integrate(exp(-x**2), (x, -oo, oo))
sqrt(pi)
solve(x**4 - 1, x)
[-1, 1, -I, I]
如你所见,第一个参数是假设等于0的表达式。它可以解一个很大的多项式方程,也可以有能力求解多个方程,可以将各自的多个变量作为元组以第二个参数给出:
solve([x + 5*y - 2, -3*x + 6*y - 15], [x, y])
{x: -3, y: 1}
也直接求解超越方程(有限的):
solve(exp(x) + 1, x)
[I*pi]
多项式方程的另一个应用是factor
。factor
将多项式因式分解为可化简的项,并且可以计算不同域的因式:
f = x**4 - 3*x**2 + 1
factor(f)
(x**2 - x - 1)*(x**2 + x - 1)
factor(f, modulus=5)
(x - 2)**2*(x + 2)**2
SymPy也可以解布尔方程,即,判断一个布尔表达式是否满足。对于这个情况,我们可以使用satisfiable
函数:
satisfiable(x & y)
{x: True, y: True}
这告诉我们(x & y)
是真,当x和y都是True的时候。如果一个表达式不是True,即它的任何参数值都无法使表达式为真,那么它将返回False:
satisfiable(x & ~x)
False
from sympy import Matrix
Matrix([[1,0], [0,1]])
Matrix([ [1, 0], [0, 1]])
与NumPy数组不同,你也可以在里面放入符号:
x = Symbol('x')
y = Symbol('y')
A = Matrix([[1,x], [y,1]])
A
Matrix([ [1, x], [y, 1]])
A**2
Matrix([ [x*y + 1, 2*x], [ 2*y, x*y + 1]])
SymPy可以解 (一些) 常规微分。要求解一个微分方程,使用dsolve
。首先,通过传递cls=Function来创建一个未定义的符号函数:
f, g = symbols('f g', cls=Function)
f 和 g是未定义函数。我们可以调用f(x), 并且它可以代表未知的函数:
f(x)
f(x)
f(x).diff(x, x) + f(x)
f(x) + Derivative(f(x), x, x)
dsolve(f(x).diff(x, x) + f(x), f(x))
f(x) == C1*sin(x) + C2*cos(x)
关键词参数可以向这个函数传递,以便帮助确认是否找到最适合的解决系统。例如,你知道它是独立的方程,你可以使用关键词hint=’separable’来强制dsolve
来将它作为独立方程来求解:
dsolve(sin(x)*cos(f(x)) + cos(x)*sin(f(x))*f(x).diff(x), f(x), hint='separable')
[f(x) == -asin(sqrt(C1/cos(x)**2 + 1)) + pi, f(x) == asin(sqrt(C1/cos(x)**2 + 1)) + pi, f(x) == -asin(sqrt(C1/cos(x)**2 + 1)), f(x) == asin(sqrt(C1/cos(x)**2 + 1))]
练习
求解Bernoulli微分方程
$x \frac{d f(x)}{x} + f(x) - f(x)^2=0$
使用hint=’Bernoulli’求解相同的公式。可以观察到什么?