Python 符号计算简介

贡献者： shizy0808; addis

1. 什么是符号计算？

　　符号计算又称计算机代数，通俗地说就是用计算机推导数学公式，如对表达式进行因式分解、化简、微分、积分、解代数方程、求解常微分方程等。在 SciPy 数值微分与积分部分我们已经介绍了如何利用 python 实现相关数值计算，这里面我们将进一步介绍符号计算在 python 中的实现。那么数值计算与符号计算有什么区别与联系呢？个人认为：首先在数值计算过程中，所有出现的 变量 或者 参数 在使用之前必须给定具体取值，并且计算结果大多数是近似的；相反的是，在符号计算过程中，变量可以预先不给定取值，计算结果是准确的，解析的。不太准确的表述为：符号计算就是对表达式进行的操作；数值计算是对数据进行的操作。

2. SymPy 库

　　在 python 中，专门进行符号计算的库是 sympy (symbol python 的简写)。利用这个库可以进行符号表达式的加减乘除等四则运算、符号化简、求导、积分、极限、解方程（组）、解微分方程（组）等等。下面我们将进行逐一介绍。

符号变量的定义

　　我们首先导入 sympy 库：

import sympy as sm

为了代码简洁，后面默认在代码之前均引入了这个库。例如

　　定义单个变量（类型为 sympy.core.symbol.Symbol）

x0 = sm.symbols('x0')

同时定义多个变量

x0 = sm.symbols('x0')
x1 = sm.symbols('x1')

这样定义理论上没有问题，但是显得繁琐。上述代码可以改写为

x0, x1 = sm.symbols('x0, x1')

对于多个变量的定义方法类似。需要说明的是 sm.symbols('x0, x1') 中的 , 用空格隔开也是可以的。即 sm.symbols('x0 x1')。

　　当定义多个连续变量时，也可以这样

x, y, z = sm.symbols('x:z')
x4, x5, x6, x7 = sm.symbols('x4:8')

注意不包括 x8.

符号常量与函数

　　 sympy 库中预置了许多常量：圆周率 $\pi$,自然常数 $e$,无穷大 $\infty$,虚数单位 $ \mathrm{i} $ 等等。它们分别为

sm.pi # 类型为 sympy.core.numbers.Pi
sm.E  # 类型为 sympy.core.numbers.Exp1
sm.I  # 类型为 sympy.core.numbers.ImaginaryUnit
sm.oo # 类型为 sympy.core.numbers.Infinity

里面也有许多常用函数，比如 $ \sin\left(x\right) $, $ \arcsin\left(x\right) $, $ \sinh\left(x\right) $, $ \exp\left(x\right) $, $ \log\left(x\right) $, $\sqrt{x}$。对应 python 代码为：

sm.sin(), sm.asin(), sm.sinh(), sm.exp(), sm.log(), sm.sqrt()
# 返回的类型由返回的表达式决定， 例如
# type(sm.sqrt(2)) = sympy.core.power.Pow
# type(sm.sqrt(4)) = sympy.core.numbers.Integer
# type(sm.sqrt(8)) = sympy.core.mul.Mul
# type(sm.sin(sm.pi)) = sympy.core.numbers.Zero
# type(sm.sin(sm.pi/2)) = sympy.core.numbers.One
# type(sm.sin(sm.pi/4)) = sympy.core.mul.Mul

上面只是列举了部分函数与常量，sympy 库只能还有非常多，感兴趣的读者可以参考官方文档另外需要注意的是，sympy 中的符号不能与 numpy, scipy 中的混用，否则会报错，例如

>>> import numpy as np
>>> x = sm.symbols('x')
>>> np.exp(x) # 尝试用 numpy 来计算 exp(x)，注意此处 x 是符号变量.

此处返回错误

Traceback (most recent call last):
File "<stdin>", line 1, in <module>
AttributeError: 'Symbol' object has no attribute 'exp'

正确的用法为

>>> sm.exp(x) # Use SymPy's version instead.
exp(x)

我们在看这个例子

>>> x = sm.symbols('x')
>>> (2/3) * sm.sin(x) # 2/3 返回一个浮点数而不是有理数.
0.666666666666667*sin(x)
>>> sm.Rational(2, 3) * sm.sin(x) # 使用有理数 2/3

为了巩固上面符号计算的一些细节我们来做一个练习。

例 1　

　　用 sympy 库的函数写一个函数返回如下表达式：

\begin{align} \frac{2}{5} \mathrm{e} ^{x^2-y} \cosh\left(x+y\right) +\frac{1}{2} \log\left(xy+1\right) ~. \end{align}

　　对应代码如下

def cal():
   x, y  = sm.symbols('x,y')
   express = sm.Rational(2,5)*sm.exp(x**2-y)* \
   sm.cosh(x+y)+Rational(1,2)*sm.log(x*y+1)
   return express

符号加与符号乘运算

　　利用 sympy 计算 $\sum_{i=1}^4{x+ \mathrm{i} y}$ 与 $\prod_{i=1}^5{x+ \mathrm{i} y}$

>>> x, y, i = sm.symbols('x y i') # 定义三个符号变量
>>> sm.summation(x + i*y, (i, 1, 4)) # 求和
4*x + 10*y
>>> sm.product(x + i*y, (i, 0, 5)) # 求积
x*(x + y)*(x + 2*y)*(x + 3*y)*(x + 4*y)*(x + 5*y)

符号表达式化简

表1：常用的化简函数

`sm.cancel()`	分子分母化简
`sm.expand()`	展开表达式
`sm.factor()`	因式分解表达式
`sm.radsimp()`	分母有理化
`sm.simplify()`	化简表达式
`sm.trigsimp()`	仅仅化简表达式的三角函数部分

　　例如，化简表达式

>>> x = sm.symbols('x')
>>> expr = (x**2 + 2*x + 1) / ((x+1)*((sm.sin(x)/sm.cos(x))**2 + 1))
>>> print(expr)
(x**2 + 2*x + 1)/((x + 1)*(sin(x)**2/cos(x)**2 + 1))
>>> sm.simplify(expr) # 化简表达式
(x + 1)*cos(x)**2

　　展开表达式

>>> expr = sm.product(x + i*y, (i, 0, 3))
>>> print(expr)
x*(x + y)*(x + 2*y)*(x + 3*y)
>>> expr_long = sm.expand(expr) # 展开表达式
>>> print(expr_long)
x**4 + 6*x**3*y + 11*x**2*y**2 + 6*x*y**3

　　约分表达式

>>> expr_long = expr_long / (x + 3*y)
>>> expr_short = sm.cancel(expr_long) # Cancel out the denominator.
x**3 + 3*x**2*y + 2*x*y**2

　　因式分解

>>> sm.factor(expr_short) 
x*(x + y)*(x + 2*y)

　　三角函数化简

>>> sm.trigsimp(2*sm.sin(x)*sm.cos(x))
sin(2*x)

需要说明的是以上各种表达式化简不会改变原有表达式。

替换与计算表达式

　　 sympy 库中提供了表达式替换函数 subs 与计算函数 evalf，它的用法为

# 定义两个符号变量
>>> x,y = sm.symbols('x y')
>>> expr = sm.expand((x + y)**3)
>>> print(expr)
x**3 + 3*x**2*y + 3*x*y**2 + y**3
# 用2*x替换原来表达式中的y
>>> expr.subs(y, 2*x)
27*x**3
# 同时替换多个变量
>>> new_expr = expr.subs({x:sm.pi, y:1})
>>> print(new_expr)
1 + 3*pi + 3*pi**2 + pi**3
# 将符号计算转为数值计算
>>> new_expr.evalf() 
71.0398678443373
>>> expr.evalf(subs={x:1, y:2})
27.0000000000000

另外，sympy 库还提供了函数 lambdify 函数，此函数相当于匿名函数，它的用法是

sm.lambdify(变量或变量列表, 符号表达式)

例如

# 定义单个符号变量
>>> f = sm.lambdify(x, sm.sin(x)**2)
>>> print(f(0), f(np.pi/2), f(np.pi), sep=' ')
0.0 1.0 1.4997597826618576e-32
# 定义多个符号变量
>>> f = sm.lambdify((x,y), sm.sin(x)**2 + sm.cos(y)**2)
>>> print(f(0,1), f(1,0), f(np.pi, np.pi), sep=' ')
0.2919265817264289 1.708073418273571 1.0

仔细的观察我们发现，虽然我们定义的是符号变量与符号表达式，但是该函数返回值是浮点型数据。默认情况下，lambdify 函数调用了 math 模块，也就是说 sm.sin() 被转化为 math.sin().

符号方程(组)求解

　　 sympy 库中提供了 solve 函数进行方程的求解，它的用法为

sm.solve(表达式，求解变量)

它的返回值是列表形式。例如

# 定义两个符号变量x y
>>> x,y = sm.symbols('x y')
# 求解方程 x^2 - 2x + 1 = 0，在第一个参数中不需要写都等于0
>>> sm.solve(x**2 - 2*x + 1, x)
[1] 
#求解x^2 - 1 = 0
>>> sm.solve(x**2 - 1, x)
[-1, 1]
>>> sm.solve(x/(y-x) + (x-y)/y, x)
[y*(-sqrt(5) + 3)/2, y*(sqrt(5) + 3)/2]
#求解方程组
>>> sm.solve([x + y - 1,x - y -3],[x,y])
 {x: 2, y: -1}

　　利用 solve_linear_system 求解线性方程组，例如

\begin{align} x+y+z=5 ~.\\ 2 x+4 y+3 z=2~. \\ 5 x+10 y+2 z=4~. \end{align}

　　对应的代码为：

>>> x, y, z = sm.symbols('x y z')
# 首先写出系数增广矩阵
>>> M = sm.Matrix([ [1, 1, 1, 5],
[2, 4, 3, 2],
[5, 10, 2, 4] ])
# 调用函数
>>> sm.solve_linear_system(M, x, y, z)
{x: 98/11, y: -45/11, z: 2/11}
# 它的返回值为字典类型。

符号求导

　　 sympy 提供了两种方法计算符号表达式的导数，分别为 sm.Derivative() 与 sm.diff(),用法为：

>>> x, y = sm.symbols('x y')
>>> f = sm.sin(y) * sm.cos(x)**2
# 对符号变量x求导
>>> df = sm.Derivative(f, x)
>>> print(df)
Derivative(sin(y)*cos(x)**2, x)
# Derivative与doit结合使用才能实际求导，
# Derivative只返回一个形式导数
>>> df.doit() 
-2*sin(x)*sin(y)*cos(x)
# diff方法直接求导
>>> sm.diff(f, x)
-2*sin(x)*sin(y)*cos(x)
# 分别对x,y,x进行多次求导
>>> sm.diff(f, x, y, x)
2*(sin(x)**2 - cos(x)**2)*cos(y) # N

符号积分

　　 sympy 提供了两种方法计算符号表达式的积分，分别为 Integral() 与 sm.integrate(),这两个方法与求导的两种方法一一对应。用法为：

# 计算不定积分
>>> sm.integrate(sm.sec(x), x)
-log(sin(x) - 1)/2 + log(sin(x) + 1)/2
# 计算定积分
>>> sm.integrate(sm.cos(x)**2, (x,0,sm.pi/2))
pi/4
#计算多重积分
>>> sm.integrate(y**2 * x**2, (x,0,2), (y,-1,1))
16/9

致读者：小时百科一直以来坚持所有内容免费，这导致我们处于严重的亏损状态。长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。因此，我们请求广大读者热心打赏 ，使网站得以健康发展。如果看到这条信息的每位读者能慷慨打赏 10 元，我们一个星期内就能脱离亏损，并保证在接下来的一整年里向所有读者继续免费提供优质内容。但遗憾的是只有不到 1% 的读者愿意捐款，他们的付出帮助了 99% 的读者免费获取知识，我们在此表示感谢。