SymPy 符号计算基本教程

Longson

好奇心泛滥

SymPy 是一个由 Python 语言编写的符号计算库。我将在本文中简要地介绍如何利用 SymPy 进行符号计算。在介绍 SymPy 之前，我们首先要明确何谓符号计算？计算机代数系统又是什么？

什么是符号计算？

处理数学对象的计算称为符号计算。在符号计算中，数学对象是精确表示的，而不是近似的，未计算的数学表达式会以符号形式保留。与符号计算相对应的是数值计算，下面将以两个例子来展示二者之间的区别。

数值计算示例

下面是一个计算 \sin\pi 数值解的例子：

import math
math.pi
print(math.sin(math.pi))

\displaystyle{1.2246467991473532e-16}

符号计算示例

下面是一个计算 \sin\pi 解析解的例子：

from sympy import *
sin(pi)

\displaystyle 0

对比 \sin\pi 的数值和符号计算结果可以发现，数值计算结果无法精确地表示出 \sin\pi = 0 ，只能用一个很小的浮点数 (1.22\times 10^{-16}) 表示，而符号计算结果则得出 \sin\pi = 0 。

明确了数值计算和符号计算之间的区别后，让我们再来认识什么是计算机代数系统。

什么是计算机代数系统？

计算机代数系统（Computer Algebra System，缩写作：CAS）是进行符号运算的软件。在计算机代数系统中运算的对象是数学表达式，通常表达式有如下几类：

多元多项式
标准函数（三角函数、指数函数等）
特殊函数（ \Gamma 函数、Bessel 函数等）
多种函数组成的复合函数
表达式的导数、积分、和与积等
级数
矩阵

以下列出了几种典型的符号计算：

表达式化简
表达式求值
表达式的变形：展开、积、幂、部分分式表示、将三角函数转换为指数函数等
一元或多元微分
带条件的化简
部分或完整的因式分解
求解线性或非线性方程
求解微分方程或差分方程
求极限
求函数的定积分、不定积分
泰勒展开、洛朗展开等
无穷级数展开
级数求和
矩阵运算
数学公式的 \TeX 或 \LaTeX 显示

通常符号计算软件也具备一定的数值运算能力，例如可以进行如下运算：

求函数确切值
求高精度值，如 \pi
线性代数的数值运算

此外符号计算软件也具有描绘二维、三维函数图像的功能。

实际上，目前存在众多的计算机代数系统，下面列出了几种：

Maple
MuPAD
Maxima
Mathcad
Mathematica
MATLAB Symbolic Math Toolbox
SageMath

为什么选择 SymPy ？

那么，是什么让 SymPy 从这众多软件中脱颖而出，让我们选择它呢？我觉得有如下 4 个原因：

SymPy 是自由软件，免费开源，在遵守许可证条款的前提下，用户可以自由地根据自己的需求修改其源代码。与之形成对比的是，Maple、MuPad、Mathcad、MATLAB、Mathematica 等都是商业软件，价格昂贵；
SymPy 使用 Python 编写而成，与使用自己发明的语言的计算机代数系统相比（如 Maxima 由 LISP 编写)，SymPy 具有很强的通用性。SymPy 完全用 Python 编写，完全在 Python 中执行。这样，只要您熟悉 Python，那么 SymPy 将会很容易入门；
与另一个使用 Python 的符号计算软件——SageMath 相比，SymPy 的优点是安装包体积小；
SymPy 的一个重要特性是，它可以作为库集成到其他软件中，为别的软件所用。SageMath 便是将 SymPy 作为其子模块，然后再加入其他模块，构建出一个功能齐全的数学软件。

准备知识

在学习如何使用 SymPy 进行符号计算之前，请确保您满足如下几个条件：

学习过微积分
学习过线性代数
熟悉 Python 基本语法
了解 Python 面向对象编程方法
会使用 JupyterLab Notebook 交互式开发环境
了解 \LaTeX 是什么东西

如何使用 SymPy ？

前面的第 1 个符号计算示例展示了如何利用 SymPy 精确地计算三角函数，实际上，它的功能远不仅于此。作为一个强大的符号计算库，它几乎能够计算所有带符号变量的表达式。下面从本节开始将介绍如何使用 SymPy。

导入 SymPy 库

在使用 SymPy 之前需要先将其导入，有两种方式：

直接导入：

import sympy

2. 利用 from 语句导入：

from sympy import *

两种方式都导入了 SymPy 库中的所有函数、对象、变量等。区别是调用方式不同。比如在调用 sqrt ( \sqrt{\;\;} )函数时，前者应写成 sympy.sqrt(2) ，后者则直接写成 sqrt(2) 。为了力求简洁，我们使用第 2 种方式导入 SymPy 。

注意：为了防止命名空间冲突，PEP 标准推荐使用第一种方式导入库。但是，通常一个符号运算 Python 源文件是单独使用的，稍加注意就可以避免命名空间冲突的问题。

新建符号

在使用符号之前，先要利用 symbols 函数定义符号，语句是：

# 新建符号 x, y
x, y = symbols('x y')

还有一个更简洁的方法是，利用 SymPy 的 abc 子模块导入所有拉丁、希腊字母：

# 利用 SymPy 的 abc 子模块新建符号 x, y
from sympy.abc import x, y

注意：希腊字母 \lambda (lambda) 是 Python 保留关键字，当用户需要使用这个字母时，请写成 lamda (不写中间的 'b')。

新建符号变量时可以指定其定义域，比如指定 x\in \mathbf{R}, x > 0 :

x = symbols('x', positive = True)

这样在求解过程中 x 必须满足这个前提条件。

可以利用 symbols 函数依次新建类似 x_{1},\,x_{2},\,x_{3},\,x_{4} 的多个变量：

vars = symbols('x_1:5')
(x_1, x_2, x_3, x_4)
vars[0]

\displaystyle x_{1}

下面是一个符号计算的完整例子：

from sympy import *
x, y, z = symbols('x y z')
y = expand((x + 1)**2) # expand() 是展开函数
y

\displaystyle x^{2} + 2 x + 1

z = Rational(1, 2) # 构造分数 1/2
z

\displaystyle \frac{1}{2}

符号计算基本操作

在本节中，我将介绍几个符号计算的基本操作。

替换

采用符号变量的 subs 方法进行替换操作，例如：

x = symbols('x')
expr = cos(x) + 1
expr.subs(x, 0)

\displaystyle 2

将字符串转换为 SymPy 表达式

利用 sympify 函数可以将字符串表达式转换为 SymPy 表达式。

注意： sympify 是符号化，与另一个函数 simplify （化简）拼写相近，不要混淆。

str_expr = 'x**2 + 2*x + 1'
expr = sympify(str_expr)
expr

\displaystyle x^{2} + 2 x + 1

转换为指定精度的数值解

可以使用符号变量的 evalf 方法将其转换为指定精度的数值解，例如：

pi.evalf(3) # pi 保留 3 位有效数字

\displaystyle 3.14

利用 `lambdify` 函数将 SymPy 表达式转换为 NumPy 可使用的函数

如果进行简单的计算，使用 subs 和 evalf 是可行的，但要获得更高的精度，则需要使用更加有效的方法。例如，要保留小数点后 1000 位，则使用 SymPy 的速度会很慢。这时，您就需要使用 NumPy 库。

lambdify 函数的功能就是可以将 SymPy 表达式转换为 NumPy 可以使用的函数，然后用户可以利用 NumPy 计算获得更高的精度。

import numpy
a = numpy.pi / 3
x = symbols('x')
expr = sin(x)
f = lambdify(x, expr, 'numpy')
f(a)
0.8660254037844386
expr.subs(x, pi/3)

\displaystyle \frac{\sqrt{3}}{2}

使用 `simplify` (化简)

在符号计算中，最常用的操作就是利用 simplify 函数对表达式化简。默认情况下， simplify 函数将自行寻找它认为的最简单的表达形式，呈现给用户。

simplify(sin(x)**2 + cos(x)**2)

\displaystyle 1

alpha_mu = symbols('alpha_mu')
simplify(2*sin(alpha_mu)*cos(alpha_mu))

\displaystyle \sin{\left(2 \alpha_{\mu} \right)}

由于 simplify 函数执行过程是启发式的，它需要寻找它认为的最简形式，所以有时它的响应会比较慢。所以，当您知道化简形式是什么类型时，不要使用 simplify 函数，而应该使用专门的函数，如 factor （后续将会介绍）。

多项式和有理函数化简

下面介绍几个用于多项式或有理函数化简的函数。

`expand` (展开)

将多项式展开，使用 expand 函数。例如：

x_1 = symbols('x_1')
expand((x_1 + 1)**2)

\displaystyle x_{1}^{2} + 2 x_{1} + 1

`factor` (因式分解)

用 factor 函数可以对多项式进行因式分解，例如：

factor(x**3 - x**2 + x - 1)

\displaystyle \left(x - 1\right) \left(x^{2} + 1\right)

实际上，多项式的展开和因式分解是互逆过程，因此 factor 和 expand 也是相对的。

`collect` (合并同类项)

利用 collect 合并同类项，例如：

expr = x*y + x - 3 + 2*x**2 - z*x**2 + x**3
collect(expr, x)

\displaystyle x^{3} + x^{2} \left(2 - z\right) + x \left(y + 1\right) - 3

`cancel` (有理分式化简)

消去分子分母的公因式使用 cancel 函数，例如：

cancel((x**2 + 2*x + 1)/(x**2 + x))

\displaystyle \frac{x + 1}{x}

`apart` (部分分式展开)

使用 apart 函数可以将分式展开，例如：

expr = (4*x**3 + 21*x**2 + 10*x + 12)/(x**4 + 5*x**3 + 5*x**2 + 4*x)
expr

\displaystyle \frac{4 x^{3} + 21 x^{2} + 10 x + 12}{x^{4} + 5 x^{3} + 5 x^{2} + 4 x}

apart(expr)

\displaystyle \frac{2 x - 1}{x^{2} + x + 1} - \frac{1}{x + 4} + \frac{3}{x}

微积分符号计算

在本节中，将介绍使用 SymPy 进行微积分的基本操作。

一元函数求导函数

求导函数使用 diff 函数，例如：

# 求一阶导数
diff(cos(x), x)

\displaystyle - \sin{\left(x \right)}

# 求 3 阶导数
diff(x**4, x, 3)

\displaystyle 24 x

我们也可以用符号变量的 diff 方法求微分，例如：

expr = cos(x)
expr.diff(x, 2)

\displaystyle - \cos{\left(x \right)}

多元函数求偏导函数

可以用 diff 函数求多元函数的偏导数，例如：

f(x, y, z) = \text{e}^{xyz}\\ \frac{\partial \, f(x, y, z)}{\partial \, x}= yz \,\text{e}^{xyz}

expr = exp(x*y*z)
diff(expr, x)

\displaystyle y z e^{x y z}

`integrate` (积分)

使用 integrate 函数求积分，例如：

# 求不定积分
integrate(cos(x), x)

\displaystyle \sin{\left(x \right)}

求 \text{e}^{-x} 的定积分：

\int_{0}^{\infty}\text{e}^{-x} \,\text{d}x

注意：在 SymPy 中，我们用 'oo' 表示 \infty 。

integrate(exp(-x), (x, 0, oo))

\displaystyle 1

求函数 f(x,y) 在 [-\infty,\,\infty] 的二重积分：

f(x,y) = \text{e}^{-x^2 - y^2}\\ \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\,f(x,y)\,\text{d}x \text{d}y

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

\displaystyle \pi

`limit` (求极限)

使用 limit 函数求极限，例如：

limit(sin(x)/x, x, 0)

\displaystyle 1

当 x\to 0^{+} 时，求 \dfrac{1}{x} 的极限：

limit(1/x, x, 0, '+')

\displaystyle \infty

`series` (级数展开)

使用符号变量的 series 方法可以对函数 f(x) 在 x=x_0 处进行 n 阶展开。例如，对函数 f(x)=\sin(x) 在 x=0 处进行 4 阶展开：

expr = sin(x)
expr.series(x, 0, 4)

\displaystyle x - \frac{x^{3}}{6} + O\left(x^{4}\right)

解方程

使用 solveset 求解方程。

求解一元二次方程

求解方程 x^2 - x =0 ，首先要构造方程，使用 Eq 函数构造等式：

Eq(x**2 - x, 0)

注意：在 SymPy 中，我们用 Eq (左边表达式, 右边表达式) 表示左边表达式与右边表达式相等。

solveset(Eq(x**2 - x, 0), x, domain = S.Reals)

\displaystyle \{0, 1\}

求解微分方程

使用 dsolve 函数求解微分方程。首先需要建立符号函数变量：

f = symbols('f', cls = Function)

然后求解微分方程：

f''(x) - 2f'(x) + f(x) = \sin(x)

diffeq = Eq(f(x).diff(x, 2) - 2*f(x).diff(x) + f(x), sin(x))
diffeq

\displaystyle f{\left(x \right)} - 2 \frac{d}{d x} f{\left(x \right)} + \frac{d^{2}}{d x^{2}} f{\left(x \right)} = \sin{\left(x \right)}

dsolve(diffeq, f(x))

\displaystyle f{\left(x \right)} = \left(C_{1} + C_{2} x\right) e^{x} + \frac{\cos{\left(x \right)}}{2}

矩阵运算

我们在进行矩阵运算之前，需要用 Matrix 构造矩阵，例如：

# 构造矩阵
Matrix([[1, -1], [3, 4], [0, 2]])

\displaystyle \left[\begin{matrix}1 & -1\\3 & 4\\0 & 2\end{matrix}\right]

# 构造列向量
Matrix([1, 2, 3])

\displaystyle \left[\begin{matrix}1\\2\\3\end{matrix}\right]

# 构造行向量
Matrix([[1], [2], [3]]).T

\displaystyle \left[\begin{matrix}1 & 2 & 3\end{matrix}\right]

矩阵转置用矩阵变量的 T 方法。

# 构造单位矩阵
eye(4)

\displaystyle \left[\begin{matrix}1 & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & 1\end{matrix}\right]

# 构造零矩阵
zeros(4)

\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\end{matrix}\right]

# 构造壹矩阵
ones(4)

\displaystyle \left[\begin{matrix}1 & 1 & 1 & 1\\1 & 1 & 1 & 1\\1 & 1 & 1 & 1\\1 & 1 & 1 & 1\end{matrix}\right]

# 构造对角矩阵
diag(1, 2, 3, 4)

\displaystyle \left[\begin{matrix}1 & 0 & 0 & 0\\0 & 2 & 0 & 0\\0 & 0 & 3 & 0\\0 & 0 & 0 & 4\end{matrix}\right]

矩阵转置

矩阵转置用矩阵变量的 T 方法。例如：

a = Matrix([[1, -1], [3, 4], [0, 2]])
a

\displaystyle \left[\begin{matrix}1 & -1\\3 & 4\\0 & 2\end{matrix}\right]

# 求矩阵 a 的转置
a.T

\displaystyle \left[\begin{matrix}1 & 3 & 0\\-1 & 4 & 2\end{matrix}\right]

求矩阵的幂

求矩阵 \boldsymbol{M} 的 2 次幂：

# 求矩阵 M 的 2 次幂
M = Matrix([[1, 3], [-2, 3]])
M**2

\displaystyle \left[\begin{matrix}-5 & 12\\-8 & 3\end{matrix}\right]

特殊地，矩阵的 -1 次幂就是矩阵的逆。

# 求矩阵 M 的逆
M**-1

\displaystyle \left[\begin{matrix}\frac{1}{3} & - \frac{1}{3}\\\frac{2}{9} & \frac{1}{9}\end{matrix}\right]

求矩阵的行列式

用矩阵变量的 det 方法可以求其行列式：

M = Matrix([[1, 0, 1], [2, -1, 3], [4, 3, 2]])
M

\displaystyle \left[\begin{matrix}1 & 0 & 1\\2 & -1 & 3\\4 & 3 & 2\end{matrix}\right]

M.det()

\displaystyle -1

求矩阵的特征值和特征多项式

用矩阵变量的 eigenvals 和 charpoly 方法求其特征值和特征多项式。

M = Matrix([[3, -2,  4, -2], [5,  3, -3, -2], [5, -2,  2, -2], [5, -2, -3,  3]])
M

\displaystyle \left[\begin{matrix}3 & -2 & 4 & -2\\5 & 3 & -3 & -2\\5 & -2 & 2 & -2\\5 & -2 & -3 & 3\end{matrix}\right]

M.eigenvals()
{3: 1, -2: 1, 5: 2}
lamda = symbols('lamda')
p = M.charpoly(lamda)
factor(p)

\displaystyle \left(\lambda - 5\right)^{2} \left(\lambda - 3\right) \left(\lambda + 2\right)

Laplace 变换

可以利用 laplace_transform 函数进行 Laplace 变换，例如：

#  Laplace (拉普拉斯)变换
from sympy.abc import t, s
expr = sin(t)
laplace_transform(expr, t, s)

\displaystyle \texttt{(1/(s**2 + 1), 0, True)}

利用 inverse_laplace_transform 函数进行逆 Laplace 变换：

expr = 1/(s - 1)
inverse_laplace_transform(expr, s, t)

\displaystyle e^{t} \theta\left(t\right)

利用 SymPy 画函数图像

使用 plot 函数绘制二维函数图像，例如：

from sympy.plotting import plot
from sympy.abc import x
plot(x**2, (x, -2, 2))

<sympy.plotting.plot.Plot at 0x20d094def40>

导入 SymPy 的 plot_implicit 函数绘制隐函数图像：

from sympy import plot_implicit

什么是符号计算 ？

数值计算示例

\displaystyle{1.2246467991473532e-16}

符号计算示例

什么是计算机代数系统 ？

为什么选择 SymPy ？

准备知识

如何使用 SymPy ？

导入 SymPy 库

新建符号

符号计算基本操作

替换

将字符串转换为 SymPy 表达式

转换为指定精度的数值解

利用 lambdify 函数将 SymPy 表达式转换为 NumPy 可使用的函数

使用 simplify (化简)

多项式和有理函数化简

expand (展开)

factor (因式分解)

collect (合并同类项)

cancel (有理分式化简)

apart (部分分式展开)

微积分符号计算

一元函数求导函数

多元函数求偏导函数

integrate (积分)

limit (求极限)

series (级数展开)

解方程

求解一元二次方程

求解微分方程

矩阵运算

矩阵转置

求矩阵的幂

求矩阵的行列式

求矩阵的特征值和特征多项式

Laplace 变换

利用 SymPy 画函数图像

什么是符号计算？

什么是计算机代数系统？

利用 `lambdify` 函数将 SymPy 表达式转换为 NumPy 可使用的函数

使用 `simplify` (化简)

`expand` (展开)

`factor` (因式分解)

`collect` (合并同类项)

`cancel` (有理分式化简)

`apart` (部分分式展开)

`integrate` (积分)

`limit` (求极限)

`series` (级数展开)