SymPy入门

什么是SymPy？SymPy是一个用于符号数学的Python库。 它旨在成为一个功能齐全的计算机代数系统（CAS），同时保持代码尽可能简单，以便易于理解和易于扩展。为什么要写这篇文章呢，Python本身已经是一个强大的计算器了，内置函数就可以进行三角函数和复数的运算等等，之前也提到了cmath库，它的全部用法也就是那么多，作为对复数运算的一个简单扩充。但是，在高等数学中，常见的极限、求导和积分运算却难以在Python中表达，因此SymPy这一符号库应运而生。谈谈我个人对“易于理解和易于扩展”的理解：SymPy与Python配合的很好，SymPy使用Python的语法来构建表达式，换言之，熟悉Python编程的可以很快地适应并掌握SymPy；同时，SymPy支持LaTeX，MathML和Dot等格式输出，有很好的扩展性，我在初次使用SymPy时是被它的简洁优雅而打动。下面的代码绝大多数来自官方文档，引用在此作为对常用属性和方法的一个总结，不求全，但求精。

abc模块

该模块的作用就是方便调用符号变量，提高效率。提供英文字母大小写和希腊字母，以下是abc模块中部分代码：

a, b, c, d, e, f, g, h, i, j = symbols('a, b, c, d, e, f, g, h, i, j')
k, l, m, n, o, p, q, r, s, t = symbols('k, l, m, n, o, p, q, r, s, t')
u, v, w, x, y, z = symbols('u, v, w, x, y, z')

A, B, C, D, E, F, G, H, I, J = symbols('A, B, C, D, E, F, G, H, I, J')
K, L, M, N, O, P, Q, R, S, T = symbols('K, L, M, N, O, P, Q, R, S, T')
U, V, W, X, Y, Z = symbols('U, V, W, X, Y, Z')

alpha, beta, gamma, delta = symbols('alpha, beta, gamma, delta')
epsilon, zeta, eta, theta = symbols('epsilon, zeta, eta, theta')
iota, kappa, lamda, mu = symbols('iota, kappa, lamda, mu')
nu, xi, omicron, pi = symbols('nu, xi, omicron, pi')
rho, sigma, tau, upsilon = symbols('rho, sigma, tau, upsilon')
phi, chi, psi, omega = symbols('phi, chi, psi, omega')

这样做的好处就是，可以很快创建符号变量，需要哪个字母(或希腊字母)直接调用就好，例如：

>>> from sympy.abc import x, y

创建符号变量

符号变量和变量有什么不同呢？打个比方，符号变量就如同数学函数f(x)中的自变量x，而变量就如同一个容器，里面存放着数值。注意，符号变量不能被赋值，但是可以计算在某一点的值，而变量是可以被赋值的，这是很大的不同。在SymPy中，存在大写开头的Symbol()和小写开头的symbols()，前者是一个类，创建Symbol对象。后者是一个操作函数，功能和结果一样，只是实现方法上是有区别的。(ps：参数传递单双引号皆可。)

>>> from sympy import *
# 单一符号化变量
>>> x = Symbol('x')
>>> x = symbols('x')
# 多重符号化变量
>>> x, y, z = symbols('x y z')
# 另一种写法，字符串可用空格或者逗号分隔
>>> x, y, z = symbols('x,y,z')

迭代创建变量

>>> symbols('x:5')
(x0, x1, x2, x3, x4)
>>> symbols('x5:x10')
(x5, x6, x7, x8, x9)
>>> symbols('x5(:2)')
(x50, x51)
>>> symbols('x5:10,y:5')
(x5, x6, x7, x8, x9, y0, y1, y2, y3, y4)
>>> symbols(('x5:10', 'y:5'))
((x5, x6, x7, x8, x9), (y0, y1, y2, y3, y4))

>>> symbols('x:z')
(x, y, z)
>>> symbols('x:c')
()
>>> symbols('x(:c)')
(xa, xb, xc)
>>> symbols(':c')
(a, b, c)
>>> symbols('a:d','x:z')
(a, b, c, d, x, y, z)
>>> symbols(('a:d','x:z'))
((a, b, c, d), (x, y, z))

创建指定特点的符号变量

>>> a = symbols('a', integer=True)
>>> x, y, z = symbols('x, y, z', real=True)
# 查看符号的属性
>>> [attr for attr in dir(x) if attr.startswith("is_") and attr.lower() == attr]
['is_algebraic',
 'is_algebraic_expr',
 'is_antihermitian',
 'is_commutative',
 'is_comparable',
 'is_complex',
 'is_composite',
 'is_constant',
 'is_even',
 'is_finite',
 'is_hermitian',
 'is_hypergeometric',
 'is_imaginary',
 'is_infinite',
 'is_integer',
 'is_irrational',
 'is_negative',
 'is_noninteger',
 'is_nonnegative',
 'is_nonpositive',
 'is_nonzero',
 'is_number',
 'is_odd',
 'is_polar',
 'is_polynomial',
 'is_positive',
 'is_prime',
 'is_rational',
 'is_rational_function',
 'is_real',
 'is_symbol',
 'is_transcendental',
 'is_zero']

基础运算

变量替换(Substitution)

>>> expr = cos(x) + 1
>>> expr.subs(x, y)
cos(y) + 1

变量替换的两个作用：

1. 求出某一点的值

   >>> expr.subs(x, 0)
   2

2. 用某个表达式去替换另一个表达式

   >>> expr = x**y
   >>> expr = expr.subs(y, x**y)
   >>> expr
   x**(x**y)

subs方法返回新的表达式，但是不改变SymPy对象。(SymPy表达式是不可变的，所有函数不改变他们，所有函数返回新的表达式。)

>>> expr = cos(x)
>>> expr.subs(x, 0)
1
>>> expr
cos(x)
>>> x
x

转换字符串为SymPy表达式

>>> str_expr = "x**2 + 3*x - 1/2"
>>> expr = sympify(str_expr)
>>> expr
x**2 + 3*x - 1/2
>>> expr.subs(x, 2)
19/2

evalf函数(evaluate floating point number)

>>> expr = sqrt(8)
>>> expr.evalf()
2.82842712474619
# 指定位数
>>> pi.evalf(20)
3.1415926535897932385
>>> expr = cos(2*x)
# a dictionary of Symbol: point pairs
>>> expr.evalf(subs={x: 2.4})
0.0874989834394464
# 浮点数过小可置chop为True
>>> one = cos(1)**2 + sin(1)**2
>>> (one - 1).evalf()
-0.e-124
>>> (one - 1).evalf(chop=True)
0

打印输出

临时交互测试SymPy的时候，可用以下命令，其中init_printing()函数会自动使用当前环境中最好的打印方式。

>>> from sympy import init_session
>>> init_session()
IPython console for SymPy 1.2 (Python 3.6.6-64-bit) (ground types: python)

These commands were executed:
from __future__ import division
from sympy import *
x, y, z, t = symbols('x y z t')
k, m, n = symbols('k m n', integer=True)
f, g, h = symbols('f g h', cls=Function)
init_printing()

Documentation can be found at http://docs.sympy.org/1.2/

在SymPy中，有以下常用的打印方式：

• str

  >>> from sympy.abc import x
  >>> print(Integral(sqrt(1/x), x))
  Integral(sqrt(1/x), x)

• srepr

• ASCII pretty printer

  >>> pprint(Integral(sqrt(1/x), x), use_unicode=False)
    /          
   |           
   |     ___   
   |    / 1    
   |   /  -  dx
   | \/   x    
   |           
  /

• Unicode pretty printer

  >>> pprint(Integral(sqrt(1/x), x), use_unicode=True)
  ⌠           
  ⎮     ___   
  ⎮    ╱ 1    
  ⎮   ╱  ─  dx
  ⎮ ╲╱   x    
  ⌡

• LaTeX

  >>> print(latex(Integral(sqrt(1/x), x)))
  \int \sqrt{\frac{1}{x}}\, dx

• MathML

• Dot

化简

simplify

函数的作用是简化数学表达式，该函数是一个通用的函数，计算会很耗时；若使用特定的函数则会提高运算速度，接下来的一些函数都是一些特定的化简方法。但是要注意它只适用主要的简化操作，因为简化的定义不是非常明确。

>>> simplify(sin(x)**2 + cos(x)**2)
1
>>> simplify(x**2 + 2*x + 1)
 2
x  + 2⋅x + 1

expand

SymPy中最常用的化简函数之一，用于展开多项表达式。可能展开与化简是对立的，因为它使表达式变得更加复杂，但有些情况下的确会变短。

>>> expand((x+1)**2)
x**2 + 2*x + 1
>>> expand((x + 2)*(x - 3))
x**2 - x - 6
>>> expand((x + 1)*(x - 2) - (x - 1)*x)
-2

factor

对于多项式，factor()是与expand()相反的。

>>> factor(x**3 - x**2 + x - 1)
(x - 1)*(x**2 + 1)
>>> factor(x**2*z + 4*x*y*z + 4*y**2*z)
z*(x + 2*y)**2

trigsimp(三角化简)

>>> trigsimp(sin(x)**2 + cos(x)**2)
1
>>> trigsimp(sin(x)**4 - 2*cos(x)**2*sin(x)**2 + cos(x)**4)
cos(4*x)/2 + 1/2
>>> trigsimp(sin(x)*tan(x)/sec(x))
sin(x)**2

微积分

求导—diff

>>> diff(sin(x), x)
cos(x)
# 可以对自变量多次求导
>>> diff(x**4, x, x, x)
24*x
# 数字跟在变量后面指明求导次数
>>> diff(x**4, x, 3)
24*x
# 求偏导
>>> expr = exp(x*y*z)
>>> diff(expr, x, y, z)
(x**2*y**2*z**2 + 3*x*y*z + 1)*exp(x*y*z)
# diff作为方法传递参数
>>> expr.diff(x, y, z)
(x**2*y**2*z**2 + 3*x*y*z + 1)*exp(x*y*z)
# 使用Derivative类，创建一个先不计算的对象
>>> deriv = Derivative(expr, x, y, z)
>>> deriv
Derivative(exp(x*y*z), x, y, z)
# 类的计算调用方法doit()
>>> deriv.doit()
(x**2*y**2*z**2 + 3*x*y*z + 1)*exp(x*y*z)

积分—integrate

积分分为定积分和不定积分。对于不定积分就是求导的逆运算，因此在表达式后传入变量即可。对于定积分，参数的格式是(integration_variable, lower_limit, upper_limit)。FYI，无穷大在SymPy中的表示为oo。例如：

>>> integrate(cos(x), x)
sin(x)
>>> integrate(exp(-x), (x, 0, oo))
1

对于双重积分，依次传入参数，例如：

>>> integrate(exp(-x**2 - y**2), (x, -oo, oo), (y, -oo, oo))
pi

如果是不可积的积分，SymPy会返回未计算的积分对象。

>>> expr = integrate(x**x, x)
>>> print(expr)
Integral(x**x, x)

对应于求导的类Derivative，积分的类是Integral。既保留了待积分式，又可以进行积分运算。

>>> expr = Integral(log(x)**2, x)
>>> expr
Integral(log(x)**2, x)
>>> expr.doit()
x*log(x)**2 - 2*x*log(x) + 2*x

极限—limit

参数格式：limit(f(x), x, x0)

>>> limit(sin(x)/x, x, 0)
1
>>> expr = sin(x)/x
# 注意此处subs的局限性
>>> expr.subs(x, 0)
nan
>>> limit(expr, x, 0)
1

同样地，极限的类是Limit。

>>> expr = Limit((cos(x) - 1)/x, x, 0)
>>> expr
Limit((cos(x) - 1)/x, x, 0)
>>> expr.doit()
0

计算单侧极限，添加第三个参数，例如求：

>>> limit(1/x, x, 0, '+')
oo
>>> limit(1/x, x, 0, '-')
-oo

级数展开

SymPy可计算函数在某一点的渐进级数展开。参数格式：f(x).series(x, x0, n)，其中x0和n是可省略的，省略时x0=0且n=6。

>>> expr = exp(sin(x))
>>> expr.series(x, 0, 4)
1 + x + x**2/2 + O(x**4)
>>> expr.series(x, 0, 4).removeO()
x**2/2 + x + 1
>>> exp(x - 6).series(x, x0=6)
-5 + (x - 6)**2/2 + (x - 6)**3/6 + (x - 6)**4/24 + (x - 6)**5/120 + x + O((x - 6)**6, (x, 6))

SymPy给开发者提供了很好的程序接口，除了以上提到的函数方法，还包括矩阵运算，解方程，解微分方程等其他诸多的内容，主要内容列举于此以供参考。