SymPy 的简单介绍

SymPy 是一个符号计算的 Python 库，完全由 Python 写成，为许多数值分析，符号计算提供了重要的工具。SymPy 的第一个版本于 2007 年开源，并且经历了十几个版本的迭代，在 2019 年已经基于修正的 BSD 许可证开源了 1.4 版本。SymPy 的开源地址和官方网站分别是：

GitHub 链接：https://github.com/sympy/sympy
SymPy 官方网站：https://www.sympy.org/en/index.html

从 SymPy 的 1.4 版本文档中，可以看出，SymPy 可以支持很多初等数学，高等数学，甚至研究生数学的符号计算。在初等数学和高等数学中，SymPy 可以支持的内容包括但不限于：

基础计算（Basic Operations）；
公式简化（Simplification）；
微积分（Calculus）；
解方程（Solver）；
矩阵（Matrices）；
几何（geometry）；
级数（Series）；

在更多的数学领域中，SymPy 可以支持的内容包括但不限于：

范畴论（Category Theory）；
微分几何（Differential Geometry）；
常微分方程（ODE）；
偏微分方程（PDE）；
傅立叶变换（Fourier Transform）；
集合论（Set Theory）；
逻辑计算（Logic Theory）。

SymPy 的工具库介绍

SymPy 的基础计算

在数学中，基础的计算包括实数和复数的加减乘除，那么就需要在程序中描述出实数与复数。著名的欧拉公式

$e^{i\pi}+1 = 0$ 。

正好用到了数学中最常见的五个实数。在 SymPy 里面， $e, i, \pi, \infty$ 是用以下符号来表示的：其中 sympy.exp() 表示以 $e$ 为底的函数。

sympy.exp(1), sympy.I, sympy.pi, sympy.oo

而想要计算欧拉公式的话，只需要输入下面的公式即可：

>>> sympy.exp(sympy.I * sympy.pi) + 1
0

如果需要看 $e, \pi$ 的小数值，可以使用 evalf() 函数，其中 evalf() 函数里面的值表示有效数字的位数。例如下面就是精确到 10 位有效数字。当然，也可以不输入。

>>> sympy.E.evalf(10)
2.718281828
>>> sympy.E.evalf()
2.71828182845905
>>> sympy.pi.evalf(10)
3.141592654
>>> sympy.pi.evalf()
3.14159265358979

除此之外，如果需要查看某个实数的有效数字，也是类似操作的：

>>> expr = sympy.sqrt(8)
>>> expr.evalf()
2.82842712474619

而对于实数的加减乘除，则可以如下操作：

>>> x, y= sympy.symbols("x y")
>>> x + y
x + y
>>> x - y
x - y
>>> x * y
x*y
>>> x / y
x/y

而对于复数的加减乘除，则是类似的操作，令两个复数分别是 $z_{1} = x_{1} + i y_{1}$ 和 $z_{2} = x_{2} + i y_{2}$ 。

>>> x1, y1, x2, y2 = sympy.symbols("x1 y1 x2 y2")
>>> z1 = x1 + y1 * sympy.I
x1 + I*y1
>>>  z2 = x2 + y2 * sympy.I
x2 + I*y2
>>> z1 + z2
x1 + x2 + I*y1 + I*y2
>>> z1 - z2
x1 - x2 + I*y1 - I*y2
>>> z1 * z2
(x1 + I*y1)*(x2 + I*y2)
>>> z1 / z2
(x1 + I*y1)/(x2 + I*y2)

对于多项式而言，有的时候我们希望将其展开，有的时候则需要将其合并，最终将其简化成最简单的形式。

>>> sympy.expand((x+1)**2)
x**2 + 2*x + 1
>>> sympy.expand((x+1)**5)
x**5 + 5*x**4 + 10*x**3 + 10*x**2 + 5*x + 1
>>> sympy.factor(x**3+1)
(x + 1)*(x**2 - x + 1)
>>> sympy.factor(x**2+3*x+2)
(x + 1)*(x + 2)
>>> sympy.simplify(x**2 + x + 1 - x)
x**2 + 1
>>> sympy.simplify(sympy.sin(x)**2 + sympy.cos(x)**2)
1

在多变量的场景下，SymPy 也可以对其中的某个变量合并同类项，同时还可以计算某个变量的某个次数所对应的系数是多少，例如：

>>> expr = x*y + x - 3 + 2*x**2 - x**2 + x**3 + y**2 + x**2*y**2
>>> sympy.collect(expr,x)
x**3 + x**2*(y**2 + 1) + x*(y + 1) + y**2 - 3
>>> sympy.collect(expr,y)
x**3 + x**2 + x*y + x + y**2*(x**2 + 1) - 3
>>> expr.coeff(x, 2)
y**2 + 1
>>> expr.coeff(y, 1)
x

有理函数形如 $f(x) = p(x)/q(x)$ ，其中 $p(x)$ 和 $q(x)$ 都是多项式。一般情况下，我们希望对有理函数进行简化，合并或者分解的数学计算。

在需要合并的情形下，如果想把有理函数处理成标准格式 $p(x)/q(x)$ 并且去除公因子，那么可以使用 cancel 函数。另一个类似的就是 together 函数，但是不同之处在于 cancel 会消除公因子，together 不会消除公因子。例如：

$expr = \frac{x^{2}+3x+2}{x^{2}+x}$ ，

>>> expr = (x**2 + 3*x + 2)/(x**2 + x)
>>> sympy.cancel(expr)
(x + 2)/x
>>> sympy.together(expr)
(x**2 + 3*x + 2)/(x*(x + 1))

除了合并和消除公因子之外，有的时候还希望对分子和分母进行因式分解，例如：

expr = (x**2 + 3*x + 2)/(x**2 + x)
>>> sympy.factor(expr)
(x + 2)/x
>>> expr = (x**3 + 3*x**2 + 2*x)/(x**5+x)
>>> sympy.factor(expr)
(x + 1)*(x + 2)/(x**4 + 1)
>>> expr = x**2 + (2*x+1)/(x**3+1)
>>> sympy.factor(expr)
(x**5 + x**2 + 2*x + 1)/((x + 1)*(x**2 - x + 1))

合并的反面就是部分分式展开（Partial Fraction Decomposition），它是把有理函数分解成多个次数较低的有理函数和的形式。这里需要用 apart 函数：

>>> expr = (x**4 + 3*x**2 + 2*x)/(x**2+x)
>>> sympy.apart(expr)
x**2 - x + 4 - 2/(x + 1)
>>> expr = (x**5 + 1)/(x**3+1)
>>> sympy.apart(expr)
x**2 - (x - 1)/(x**2 - x + 1)

在 SymPy 里面，同样支持各种各样的三角函数，例如：三角函数的简化函数 trigsimp，三角函数的展开 expand_trig，

>>> expr = sympy.sin(x)**2 + sympy.cos(x)**2
>>> sympy.trigsimp(expr)
1
>>> sympy.expand_trig(sympy.sin(x+y))
sin(x)*cos(y) + sin(y)*cos(x)
>>> sympy.expand_trig(sympy.cos(x+y))
-sin(x)*sin(y) + cos(x)*cos(y)
>>> sympy.trigsimp(sympy.sin(x)*sympy.cos(y) + sympy.sin(y)*sympy.cos(x))
sin(x + y)
>>> sympy.trigsimp(-sympy.sin(x)*sympy.sin(y) + sympy.cos(x)*sympy.cos(y))
cos(x + y)

同样的，在乘幂上面，同样有简化函数 powsimp，效果与之前提到的 simplify 一样。除此之外，还可以根据底数来做合并，即分别使用 expand_power_exp 函数与 expand_power_base 函数。

>>> sympy.powsimp(x**z*y**z*x**z)
x**(2*z)*y**z
>>> sympy.simplify(x**z*y**z*x**z)
x**(2*z)*y**z
>>> sympy.expand_power_exp(x**(y + z))
x**y*x**z
>>> sympy.expand_power_base(x**(y + z))
x**(y + z)

作为指数的反函数对数，sympy 也是有着类似的展开合并函数，expand_log，logcombine 承担了这样的角色。

$\ln(xy) = \ln(x) + \ln(y)$ ，

$\ln(x/y) = \ln(x) - \ln(y)$ 。

>>> sympy.expand_log(sympy.log(x*y), force=True)
log(x) + log(y)
>>> sympy.expand_log(sympy.log(x/y), force=True)
log(x) - log(y)

SymPy 的微积分工具

下面，我们会从一个最基本的函数 $f(x) = 1/x$ 出发，来介绍 SymPy 的各种函数使用方法。如果想进行变量替换，例如把 $x$ 变成 $y$ ，那么可以使用 substitution 方法。除此之外，有的时候也希望能够得到函数 $f$ 在某个点的取值，例如 $f(1) = 1/1 = 1$ ，那么可以把参数换成 1 即可得到函数的取值。例如，

>>> import sympy
>>> x = sympy.Symbol("x")
>>> f = 1 / x
1/x
>>> y = sympy.Symbol("y")
>>> f = f.subs(x,y)
1/y
>>> f = f.subs(y,1)
1

在微积分里面，最常见的概念就是极限，SymPy 里面的极限函数就是 limit。使用方法如下：

>>> f = 1/x
>>> sympy.limit(f,x,0)
oo
>>> sympy.limit(f,x,2)
1/2
>>> sympy.limit(f,x,sympy.oo)
0
>>> g = x * sympy.log(x)
>>> sympy.limit(g,x,0)
0

对于函数 $f(x) = 1/x$ 而言，它的导数计算函数是 diff， $n$ 阶导数也可以用这个函数算。

>>> f = 1/x
>>> sympy.diff(f,x)
-1/x**2
>>> sympy.diff(f,x,2)
2/x**3
>>> sympy.diff(f,x,3)
-6/x**4
>>> sympy.diff(f,x,4)
24/x**5

提到 $n$ 阶导数，就必须要提一下 Taylor Series 了。对于常见函数的 Taylor Series，SymPy 也是有非常简便的方法，那就是 series 函数。其参数包括 expr, x, x0, n, dir，分别对应着表达式，函数的自变量，Taylor Series 的中心点，n 表示阶数，dir 表示方向，包括”+-“，”-“，”+”，分别表示 $x\rightarrow x0$ , $x\rightarrow x0^{-}$ , $x\rightarrow x0^{+}$ 。

sympy.series.series.series(expr, x=None, x0=0, n=6, dir='+') >>> g = sympy.cos(x) >>> sympy.series(g, x) 1 - x**2/2 + x**4/24 + O(x**6) >>> sympy.series(g, x, x0=1, n=10) cos(1) - (x - 1)*sin(1) - (x - 1)**2*cos(1)/2 + (x - 1)**3*sin(1)/6 + (x - 1)**4*cos(1)/24 - (x - 1)**5*sin(1)/120 - (x - 1)**6*cos(1)/720 + (x - 1)**7*sin(1)/5040 + (x - 1)**8*cos(1)/40320 - (x - 1)**9*sin(1)/362880 + O((x - 1)**10, (x, 1))

积分的计算函数是 integrate，包括定积分与不定积分：

$\int\frac{1}{x}dx = \ln(x)+C$ ，

$\int_{1}^{2}\frac{1}{x}dx = \ln(2)$ 。

>>> f = 1/x
>>> sympy.integrate(f,x)
log(x)
>>> sympy.integrate(f, (x,1,2))
log(2)

对于广义积分而言，就需要用到 $\infty$ 这个概念了，但是在 SymPy 里面的写法还是一样的。

$\int_{-\infty}^{0}e^{-x^{2}}dx=\frac{\sqrt{\pi}}{2}$ ，

$\int_{0}^{+\infty}e^{-x}dx = 1$ ，

$\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}e^{-x^{2}-y^{2}}dxdy = \pi$ 。

>>> g = sympy.exp(-x**2)
>>> sympy.integrate(g, (x,-sympy.oo,0))
sqrt(pi)/2
>>> g = sympy.exp(-x)
>>> sympy.integrate(g, (x, 0, sympy.oo))
1
>>> h = sympy.exp(-x**2 - y**2)
>>> sympy.integrate(h, (x,-sympy.oo, sympy.oo), (y, -sympy.oo, sympy.oo))
pi

SymPy 的方程工具

在初等数学中，通常都存在一元一次方程，一元二次方程等，并且在不同的域上有着不同的解。SymPy 里面的相应函数就是 solveset，根据定义域的不同，可以获得完全不同的解。

$\{x\in\mathbb{R}: x^{3}-1=0\}$ ，

$\{x\in\mathbb{C}:x^{3}-1=0\}$ ，

$\{x\in\mathbb{R}:e^{x}-x=0\}$ ，

$\{x\in\mathbb{R}:e^{x}-1=0\}$ ，

$\{x\in\mathbb{C}:e^{x}-1=0\}$ 。

>>> sympy.solveset(sympy.Eq(x**3,1), x, domain=sympy.S.Reals)
{1}
>>> sympy.solveset(sympy.Eq(x**3,1), x, domain=sympy.S.Complexes)
{1, -1/2 - sqrt(3)*I/2, -1/2 + sqrt(3)*I/2}
>>> sympy.solveset(sympy.Eq(x**3 - 1,0), x, domain=sympy.S.Reals)
{1}
>>> sympy.solveset(sympy.Eq(x**3 - 1,0), x, domain=sympy.S.Complexes)
{1, -1/2 - sqrt(3)*I/2, -1/2 + sqrt(3)*I/2}
>>> sympy.solveset(sympy.exp(x),x)
EmptySet()
>>> sympy.solveset(sympy.exp(x)-1,x,domain=sympy.S.Reals)
{0}
>>> sympy.solveset(sympy.exp(x)-1,x,domain=sympy.S.Complexes)
ImageSet(Lambda(_n, 2*_n*I*pi), Integers)

在这里，Lambda 表示的是数学公式，第一个是自变量，第二个是函数，最后是自变量的定义域。

在线性代数中，最常见的还是多元一次方程组，那么解法是一样的：

$\begin{cases}x+y-10=0 \\ x-y-2=0\end{cases}$

>>> sympy.solve([x+y-10, x-y-2], [x,y])
{x: 6, y: 4}

对于三角函数，也是类似的写法：

$\begin{cases} \sin(x-y)=0 \\ \cos(x+y)=0 \end{cases}$ ，

>>> sympy.solve([sympy.sin(x-y), sympy.cos(x+y)], [x,y])
[(-pi/4, 3*pi/4), (pi/4, pi/4), (3*pi/4, 3*pi/4), (5*pi/4, pi/4)]

SymPy 的矩阵工具

在矩阵论中，最常见的就是单位矩阵了，而单位矩阵只与一个参数有关，那就是矩阵的大小。下面就是 3*3，3*2，2*3 大小的矩阵。

>>> sympy.eye(3)
Matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]])
>>> sympy.eye(3,2)
Matrix([
[1, 0],
[0, 1],
[0, 0]])
>>> sympy.eye(2,3)
Matrix([
[1, 0, 0],
[0, 1, 0]])

另外还有全零和全一矩阵，意思就是矩阵中的所有值全部是零和一。

>>> sympy.ones(2,3)
Matrix([
[1, 1, 1],
[1, 1, 1]])
>>> sympy.zeros(3,2)
Matrix([
[0, 0],
[0, 0],
[0, 0]])

而对角矩阵也可以使用 diag 轻松获得：

>>> sympy.diag(1,1,2)
Matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 2]])

而矩阵的加法，减法，乘法，逆运算，转置，行列式，SymPy 都是可以支持的：

$A = \left(\begin{array}{cc} 1 & 1 \\ 0 & 2 \end{array}\right)$ ，

$B = \left(\begin{array}{cc} 1 & 0 \\ 1 & 1 \end{array}\right)$ ，

>>> A = sympy.Matrix([[1,1],[0,2]])
>>> B = sympy.Matrix([[1,0],[1,1]])
>>> A
Matrix([
[1, 1],
[0, 2]])
>>> B
Matrix([
[1, 0],
[1, 1]])
>>> A + B
Matrix([
[2, 1],
[1, 3]])
>>> A - B
Matrix([
[ 0, 1],
[-1, 1]])
>>> A * B
Matrix([
[2, 1],
[2, 2]])
>>> A * B**-1
Matrix([
[ 0, 1],
[-2, 2]])
>>> B**-1
Matrix([
[ 1, 0],
[-1, 1]])
>>> A.T
Matrix([
[1, 0],
[1, 2]])
>>> A.det()
2

在某些情况下，需要对矩阵进行加上一行或者加上一列的操作，在这里就可以使用 row_insert 或者 col_insert 函数：第一个参数表示插入的位置，第二个参数就是相应的行向量或者列向量。而在删除上就很简单了，直接使用 row_del 或者 col_del 即可。

>>> A
Matrix([
[1, 1],
[0, 2]])
>>> A.row_insert(1, sympy.Matrix([[10,10]]))
Matrix([
[ 1, 1],
[10, 10],
[ 0, 2]])
>>> A.col_insert(0, sympy.Matrix([3,3]))
Matrix([
[3, 1, 1],
[3, 0, 2]])
>>> A.row_del(0)
>>> A
Matrix([[0, 2]])
>>> A.col_del(1)
>>> A
Matrix([[0]])

在对角化方面，同样可以使用 eigenvals()，eigenvecs(), diagonalize() 函数：

>>> A
Matrix([
[1, 1],
[0, 2]])
>>> A.eigenvals()
{2: 1, 1: 1}
>>> A.eigenvects()
[(1, 1, [Matrix([
[1],
[0]])]), (2, 1, [Matrix([
[1],
[1]])])]
>>> A.diagonalize()
(Matrix([
[1, 1],
[0, 1]]), Matrix([
[1, 0],
[0, 2]]))

在 eigenvals() 返回的结果中，第一个表示特征值，第二个表示该特征值的重数。在特征向量 eigenvecs() 中，第一个表示特征值，第二个表示特征值的重数，第三个表示特征向量。在对角化 diagonalize() 中，第一个矩阵表示 $P$ ，第二个矩阵表示 $D$ ， $A = P*D*P^{-1}$ 。

在矩阵中，最常见的还是多元一次方程的解。如果要求 $Ax =b$ 的解，可以有以下方案：

>>> A = sympy.Matrix([[1,1],[0,2]])
>>> A
Matrix([
[1, 1],
[0, 2]])
>>> b = sympy.Matrix([3,5])
>>> b
Matrix([
[3],
[5]])
>>> A**-1*b
Matrix([
[1/2],
[5/2]])
>>> sympy.linsolve((A,b))
{(1/2, 5/2)}
>>> sympy.linsolve((A,b),[x,y])
{(1/2, 5/2)}

SymPy 的集合论工具

集合论可以说是数学的基础，在任何数学的方向上都能够看到集合论的身影。在 SymPy 里面，有一个类叫做 sympy.sets.sets.set。在集合论里面，常见的就是边界，补集，包含，并集，交集等常见的操作。但是感觉 SymPy 中的集合论操作主要集中在实数域或者复数域上。

对于闭区间 $I=[0,1]$ 和开区间 $J = (0,1)$ 而言，在 SymPy 中使用以下方法来表示：

I = sympy.Interval(0,1)
J = sympy.Interval.open(0,1)
K = sympy.Interval(0.5,2)

其开始和结束的点可以分别使用 start 和 end 来表示：

>>> I.start
0
>>> I.end
1

其长度用 measure 来表示：

>>> I.measure
1

其闭包用 closure 来表示：

>>> I.closure
Interval(0, 1)

其内点用 interior 来表示：

>>> I.interior
Interval.open(0, 1)

判断其边界条件可以使用 left_open 或者 right_open 来做：

>>> I.left_open
False
>>> I.right_open
False

对于两个集合之间的操作，可以参考以下方法：

I = sympy.Interval(0,1)
K = sympy.Interval(0.5,2)
>>> I.intersect(K)
Interval(0.500000000000000, 1)
>>> I.union(K)
Interval(0, 2)
>>> I-K
Interval.Ropen(0, 0.500000000000000)
>>> K-I
Interval.Lopen(1, 2)
>>> I.symmetric_difference(K)
Union(Interval.Ropen(0, 0.500000000000000), Interval.Lopen(1, 2))

实数集 $\mathbb{R}$ 在 SymPy 中用 sympy.S.Reals 来表示，自然数使用 sympy.S.Naturals，非负整数用 sympy.S.Naturals0，整数用 sympy.S.Integers 来表示。补集的计算可以用减号，也可以使用 complement 函数。

>>> sympy.S.Reals
Reals
>>> sympy.S.Reals-I
Union(Interval.open(-oo, 0), Interval.open(1, oo))
>>> I.complement(sympy.S.Reals)
Union(Interval.open(-oo, 0), Interval.open(1, oo))
>>> sympy.S.Reals.complement(I)
EmptySet()
>>> I.complement(K)
Interval.Lopen(1, 2)
>>> I.complement(sympy.S.Reals)
Union(Interval.open(-oo, 0), Interval.open(1, oo))

SymPy 的逻辑工具

在逻辑运算中，我们可以使用 $A, B, C$ 来代表元素。&, |, ~, >> 分别表示 AND，OR，NOT，imply。而逻辑运算同样可以使用 sympy.simplify_logic 简化。

A, B, C = sympy.symbols("A B C")
>>> sympy.simplify_logic(A | (A & B))
A
>>> sympy.simplify_logic((A>>B) & (B>>A))
(A & B) | (~A & ~B)
>>> A>>B
Implies(A, B)

SymPy 的级数工具

SymPy 的级数工具有一部分放在具体数学（Concrete Mathematics）章节了。有的时候，我们希望计算某个级数是发散的，还是收敛的，就可以使用 is_convergence 函数。考虑最常见的级数：

$\sum_{n=1}^{\infty}\frac{1}{n} = +\infty$ ，

$\sum_{n=1}^{\infty}\frac{1}{n^{2}} = \frac{\pi^{2}}{6}$ 。

>>> n = sympy.Symbol("n", integer=True)
>>> sympy.Sum(1/n, (n,1,sympy.oo)).is_convergent()
False
>>> sympy.Sum(1/n**2, (n,1,sympy.oo)).is_convergent()
True

如果想计算出收敛级数的值，加上 doit() 函数即可；如果想计算有效数字，加上 evalf() 函数即可。

>>> sympy.Sum(1/n**2, (n,1,sympy.oo)).evalf()
1.64493406684823
>>> sympy.Sum(1/n**2, (n,1,sympy.oo)).doit()
pi**2/6
>>> sympy.Sum(1/n**3, (n,1,sympy.oo)).evalf()
1.20205690315959
>>> sympy.Sum(1/n**3, (n,1,sympy.oo)).doit()
zeta(3)

除了加法之外，SymPy 也支持连乘，其符号是 sympy.Product，考虑

$\prod_{n=1}^{+\infty}\frac{n}{n+1}$ ，

$\prod_{n=1}^{+\infty}\cos\left(\frac{\pi}{n}\right)$ 。

>>> sympy.Product(n/(n+1), (n,1,sympy.oo)).is_convergent()
False
>>> sympy.Product(sympy.cos(sympy.pi/n), (n, 1, sympy.oo)).is_convergent()
True

SymPy 的 ODE 工具

在常微分方程（Ordinary Differential Equation）中，最常见的就是解方程，而解方程主要是靠 dsolve 函数。例如想求解以下的常微分方程：

$df/dx + f(x) = 0$ ,

$d^{2}f/dx^{2} + f(x) = 0$ ，

$d^{3}f/dx^{3} + f(x) = 0$ ，

可以使用 dsolve 函数：

>>> f = sympy.Function('f')
>>> sympy.dsolve(sympy.Derivative(f(x),x) + f(x), f(x))
Eq(f(x), C1*exp(-x))
>>> sympy.dsolve(sympy.Derivative(f(x),x,2) + f(x), f(x))
Eq(f(x), C1*sin(x) + C2*cos(x))
>>> sympy.dsolve(sympy.Derivative(f(x),x,3) + f(x), f(x))
Eq(f(x), C3*exp(-x) + (C1*sin(sqrt(3)*x/2) + C2*cos(sqrt(3)*x/2))*sqrt(exp(x)))

而常微分方程对于不同的方程类型也有着不同的解法，可以使用 classify_ode 来判断常微分方程的类型：

>>> sympy.classify_ode(sympy.Derivative(f(x),x) + f(x), f(x))
('separable', '1st_exact', '1st_linear', 'almost_linear', '1st_power_series', 'lie_group', 'nth_linear_constant_coeff_homogeneous', 'separable_Integral', '1st_exact_Integral', '1st_linear_Integral', 'almost_linear_Integral')
>>> sympy.classify_ode(sympy.Derivative(f(x),x,2) + f(x), f(x))
('nth_linear_constant_coeff_homogeneous', '2nd_power_series_ordinary')
>>> sympy.classify_ode(sympy.Derivative(f(x),x,3) + f(x), f(x))
('nth_linear_constant_coeff_homogeneous',)

SymPy 的 PDE 工具

在偏微分方程（Partitial Differential Equation）中，同样可以直接求解和判断偏微分方程的类型，分别使用函数 pdsolve() 和 classify_pde()。假设 $f = f(x,y)$ 是一个二元函数，分别满足以下偏微分方程：

$\partial f/\partial x + \partial f/\partial y =0$ ，

$\partial f/\partial x + \partial f/\partial y + f = 0$ ，

$\partial f/\partial x + \partial f/\partial y + f + 10 = 0$ 。

>>> f = sympy.Function("f")(x,y)
>>> sympy.pdsolve(sympy.Derivative(f,x)+sympy.Derivative(f,y),f)
Eq(f(x, y), F(x - y))
>>> sympy.pdsolve(f.diff(x)+f.diff(y)+f,f)
Eq(f(x, y), F(x - y)*exp(-x/2 - y/2))
>>> sympy.pdsolve(f.diff(x)+f.diff(y)+f+10,f)
Eq(f(x, y), F(x - y)*exp(-x/2 - y/2) - 10)

查看类型就用 classify_pde() 函数：

>>> sympy.classify_pde(f.diff(x)+f.diff(y)+f)
('1st_linear_constant_coeff_homogeneous',)
>>> sympy.classify_pde(f.diff(x)+f.diff(y)+f+10,f)
('1st_linear_constant_coeff', '1st_linear_constant_coeff_Integral')
>>> sympy.classify_pde(f.diff(x)+f.diff(y)+f+10,f)
('1st_linear_constant_coeff', '1st_linear_constant_coeff_Integral')

不过目前的 PDE 解法貌似只支持一阶偏导数，二阶或者以上的偏导数就不支持了。

SymPy 的数论工具

在数论中，素数就是一个最基本的概念之一。而素数的批量计算，比较快的方法就是筛法（sieve method）。在 sympy 中，同样有 sympy.sieve 这个工具，用于计算素数。如果想输出前100个素数，那么

>>> sympy.sieve._reset()
>>> sympy.sieve.extend_to_no(100)
>>> sympy.sieve._list
array('l', [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631])

如果想输出一个区间内的所有素数，可以使用 primerange(a,b) 函数：

>>> [i for i in sympy.sieve.primerange(10,100)]
[11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

search() 函数是为了计算某个数附近是第几个素数：

>>> sympy.sieve.search(10)
(4, 5)
>>> sympy.sieve.search(11)
(5, 5)

如果只想获得第 n 个素数，则使用函数 sympy.ntheory.generate.prime(n) 即可。如果是希望计算 x 后面的下一个素数，使用 sympy.ntheory.generate.nextprime(x) 即可。判断 x 是否是素数，可以使用 sympy.ntheory.generate.isprime(x)。

>>> sympy.ntheory.generate.prime(10)
29
>>> sympy.ntheory.generate.nextprime(10)
11
>>> sympy.ntheory.generate.nextprime(11)
13
>>> sympy.ntheory.generate.isprime(11)
True
>>> sympy.ntheory.generate.isprime(12)
False

除此之外，SymPy 的数论方法还有很多，需要读者根据 SymPy 的官方文档自行探索。

SymPy 的范畴论工具

SymPy 还支持范畴论（Category Theory）的一些计算方法，在这里简要地列举一下。

>>> A = sympy.categories.Object("A")
>>> B = sympy.categories.Object("B")
>>> f = sympy.categories.NamedMorphism(A,B,"f")
>>> f.domain
Object("A")
>>> f.codomain
Object("B")

由于范畴论是数学的“黑话”，因此其余方法留给范畴论的科研人员自行翻阅。

总结：

整体来看，SymPy 是一个非常卓越的 Python 开源符号计算库。在符号计算领域，不仅支持常见的微积分，线性代数，几何运算，还支持集合论，微分方程，数论等诸多数学方向。后续笔者将会持续跟进并研究这一卓越的开源工具库。

参考文献：

Meurer A, Smith C P, Paprocki M, et al. SymPy: symbolic computing in Python[J]. PeerJ Computer Science, 2017, 3: e103.
GitHub：https://github.com/sympy/sympy
SymPy：https://www.sympy.org/en/index.html
Sympy 维基百科：https://en.wikipedia.org/wiki/SymPy
GreatX’s Blog：数值 Python：符号计算：https://vlight.me/2018/04/01/Numerical-Python-Symbolic-Computing/
SymPy 符号计算-让Python帮我们推公式：https://zhuanlan.zhihu.com/p/83822118
Python 科学计算利器—SymPy库：https://www.jianshu.com/p/339c91ae9f41

ZHANG RONG

Monthly Archives: December 2019

用 Python 来研究数学 — SymPy 符号工具包介绍