sympy_python代数符号运算

1 SymPy介绍

SymPy 是一个符号计算的 Python 库，完全由 Python 写成，为许多数值分析，符号计算提供了重要的工具。它是一种计算机代数系统computer algebra system (CAS), 既可以用作单独的应用也可以作为其它的应用中的包。SymPy 的开源地址和官方网站分别是：

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

SymPy 包含大量的可用功能, 涵盖了基本符号代数, 微积分, 代数学, 离散数学, 量子物理等学科. SymPy 可将结果格式化输出为 LaTeX, MathML等样式。

SymPy 可以支持的内容包括但不限于：

• 基础计算（Basic Operations）
• 公式简化（Simplification）
• 微积分（Calculus）
• 解方程（Solver）
• 矩阵（Matrices）
• 几何（geometry）
• 级数（Series）
• 范畴论（Category Theory）
• 微分几何（Differential Geometry）
• 常微分方程（ODE）
• 偏微分方程（PDE）
• 傅立叶变换（Fourier Transform）
• 集合论（Set Theory）
• 逻辑计算（Logic Theory）
• 统计（Statistics）
• 组合数学（Combinatorics）
• 作图（Plotting）

1.1 安装 SymPy

直接使用 pip 安装即可：

pip install sympy

2 基础计算

2.1 常数符号表示

在数学中，基础的计算包括实数和复数的加减乘除，那么就需要在程序中描述出实数与复数。著名的欧拉公式 $e^{i\pi }+1=0$，正好用到了数学中最常见的五个实数。在 SymPy 里面，$e,i,\pi ,\infty$ 用以下符号来表示的：

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

其中sympy.exp()表示以 $e$ 为底的函数。

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

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

2.2 evalf() 函数查看有效位数

如果需要看 $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

2.3 数字的加减乘除

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

>>> 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)

2.3 多项式处理

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

>>> 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)=\frac{p(x)}{q(x)}$，其中 $p(x),q(x)$ 都是多项式。一般情况下，我们希望对有理函数进行简化，合并或者分解的数学计算。在需要合并的情形下，如果想把有理函数处理成标准格式 $\frac{p(x)}{q(x)}$可以使用 cancel 函数。另一个类似的就是 together 函数，但是不同之处在于 cancel 会消除公因子，together 不会消除公因子。例如: $\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)

2.4 三角函数变换trigsimp，expand_trig

在 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)

2.5 幂函数化简

同样的，在乘幂上面，同样有简化函数 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)

2.6 对数化简expand_log，logcombine

作为指数的反函数对数，sympy 也是有着类似的展开合并函数，expand_log，logcombine 承担了这样的角色。
$$\begin{gathered} \ln (xy)=\ln (x)+\ln (y) \\ \ln (x/y)=\ln (x)-\ln (y) \end{gathered}$$

>>> 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)

3 微分工具

3.1 变量替换substitution

对于函数 $f(x)=\frac{1}{x}$，如果想进行变量替换，例如把 $x$ 替换成 $y$，那么可以使用 substitution 方法；如果需要计算 $f(x)$ 在 $x=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

3.2 极限 limit

在微积分里面，最常见的概念就是极限，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

3.3 导数 diff

对于函数 $f(x)=\frac{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

3.4 泰勒展开series

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

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))

3.5 不定积分integrate

积分的计算函数是 integrate，包括定积分与不定积分：
$$\begin{gathered} \int \frac{1}{x}dx=\ln (x)+C \\ {\int }_1^2\frac{1}{x}dx=\ln (2) \end{gathered}$$

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

对于广义积分，
$$\begin{gathered} {\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 \end{gathered}$$

>>> 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

4 方程工具

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

$$\begin{gathered} 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 \end{gathered}$$

>>> 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{array}{l}x+y-10=0,\\ x-y-2=0\end{array}$$

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

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

$$\{\begin{array}{l}sin(x-y)=0,\\ cos(x+y)=0\end{array}$$

>>> 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)]

5 矩阵工具

在矩阵论中，最常见的就是单位矩阵了，而单位矩阵只与一个参数有关，那就是矩阵的大小。下面就是 33，32，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 都是可以支持的：

$$\begin{gathered} A=\left(\begin{array}{c}11\\ 02\end{array}\right), \\ B=\left(\begin{array}{c}10\\ 11\end{array}\right) \end{gathered}$$

>>> 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\ast D\ast 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)}

6 集合论工具

在 SymPy 里面，有一个类叫做 sympy.sets.sets.set。在集合论里面，常见的就是边界，补集，包含，并集，交集等常见的操作。

对于闭区间 $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))

7 逻辑工具

在逻辑运算中，我们可以使用 $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)

8 级数工具

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

$$\begin{gathered} \sum _{n=1}^{\infty }\frac{1}{n}=+\infty \\ \sum _{n=1}^{\infty }\frac{1}{n^2}=\frac{{\pi }^2}{6} \end{gathered}$$

>>> 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，考虑

$$\begin{gathered} \prod _{n=1}^{\infty }\frac{n}{n+1}, \\ \prod _{n=1}^{\infty }\cos (\frac{\pi }{n}) \end{gathered}$$

>>> 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

9 ODE 工具

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

$$\begin{gathered} \frac{df}{dx}+f(x)=0, \\ \frac{d^{2}f}{d{x}^2}+f(x)=0, \\ \frac{d^{3}f}{d{x}^3}+f(x)=0 \end{gathered}$$

可以使用 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',)

10 PDE 工具

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

$$\begin{gathered} \frac{\partial f}{\partial x}+\frac{\partial f}{\partial y}=0, \\ \frac{\partial f}{\partial x}+\frac{\partial f}{\partial y}+f=0, \\ \frac{\partial f}{\partial x}+\frac{\partial f}{\partial y}+f+10=0, \end{gathered}$$

>>> 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 解法貌似只支持一阶偏导数，二阶或者以上的偏导数就不支持了。

11 数论工具

在数论中，素数就是一个最基本的概念之一。而素数的批量计算，比较快的方法就是筛法（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 的官方文档自行探索。

12 范畴论工具

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")

13 绘图

SymPy 使用 Matplotlib 库作为后端来渲染2维和3维数学函数图像. 要确保在当前安装的Python 中Matplotlib 库是可用的。

支持的作图命令定义在模块 sympy.plotting中. 如下的函数位于 plotting 模块：

• plot, 绘制二维曲线
• plot3d，绘制三维线型图
• plot_parametric，二维参数绘图
• plot3d_parametric，三维参数绘图

13.1 plot()

函数 plot() 返回一个 Plot 类的实例. 一个 plot 命令可以包含一个或多个 SymPy 表达式. 尽管可以在后端使用 Matplotlib, 其它的后端(诸如 texplot, pyglet 或者 Google charts API) 也可以使用.

plot(expr, range, kwargs)

其中 expr 是任何有效的 symPy 表达式. 如果没有提及, 默认使用的作图使用的范围是 (-10, 10).

下例中对于在范围(-10,10)中的每个值, 做出了 x^2 的值

>>> from sympy.plotting import plot 
>>> from sympy import * 
>>> x=Symbol('x') 
>>> plot(x**2, line_color='red')

为了在相同范围内画出多个图像, 在范围之前给出多个表达式.

>>> plot( sin(x),cos(x), (x, -pi, pi))

可以对每个表达式指定专门的范围

plot((expr1, range1), (expr2, range2))

>>> plot( (sin(x),(x, -2*pi, 2*pi)),(cos(x), (x, -pi, pi)))

下面的关键词参数可以在plot() 函数中指定：

• line_color − 指定绘图线的颜色.
• title − 作为标题显示的字符串
• xlabel − 作为x轴标签显示的字符串
• ylabel − 作为y轴标签显示的字符串

13.2 plot3d()

plot3d(expr, xrange, yrange, kwargs)

>>> from sympy.plotting import plot3d 
>>> x,y=symbols('x y') 
>>> plot3d(x*y, (x, -10,10), (y, -10,10))

和2维作图类似, 三维图形命令plot3d也可以对多个函数在不同范围内分别作图.

>>> plot3d(x*y, x/y, (x, -5, 5), (y, -5, 5))

13.3 plot3d_parametric_line()

>>> from sympy.plotting import plot3d_parametric_line 
>>> plot3d_parametric_line(cos(x), sin(x), x, (x, -5, 5))

13.3 plot3d_parametric_surface()

为了做出参数曲面使用函数 plot3d_parametric_surface(). 其语法是:

plot3d_parametric_surface(xexpr, yexpr, zexpr, rangex, rangey, kwargs)

>>> from sympy.plotting import plot3d_parametric_surface 
>>> plot3d_parametric_surface(cos(x+y), sin(x-y), x-y, (x, -5, 5), (y, -5, 5))