3.2 Sympy：Python中的符号数学
最后更新于：2022-04-01 11:21:59
# 3.2 Sympy：Python中的符号数学
> **作者** : Fabian Pedregosa
**目的**
* 从任意的精度评估表达式。
* 在符号表达式上进行代数运算。
* 用符号表达式进行基本的微积分任务 (极限、微分法和积分法)。
* 求解多项式和超越方程。
* 求解一些微分方程。
为什么是SymPy? SymPy是符号数学的Python库。它的目的是成为Mathematica或Maple等系统的替代品，同时让代码尽可能简单并且可扩展。SymPy完全是用Python写的，并不需要外部的库。
Sympy文档及库安装见[http://www.sympy.org/](http://www.sympy.org/)
**章节内容**
* SymPy第一步
* 使用SymPy作为计算器
* 练习
* 符号
* 代数运算
* 展开
* 化简
* 微积分
* 极限
* 微分法
* 序列扩展
* 积分法
* 练习
* 方程求解
* 练习
* 线性代数
* 矩阵
* 微分方程
## 3.2.1 SymPy第一步
### 3.2.1.1 使用SymPy作为计算器
SymPy定义了三种数字类型：实数、有理数和整数。
有理数类将有理数表征为两个整数对: 分子和分母，因此Rational(1,2)代表1/2, Rational(5,2)代表5/2等等:
In [2]:
```
from sympy import *
a = Rational(1,2)
```
In [2]:
```
a
```
Out[2]:
```
1/2
```
In [3]:
```
a*2
```
Out[3]:
```
1
```
SymPy在底层使用mpmath, 这使它可以用任意精度的算术进行计算。这样，一些特殊的常数，比如e, pi, oo (无限), 可以被作为符号处理并且可以以任意精度来评估:
In [4]:
```
pi**2
```
Out[4]:
```
pi**2
```
In [5]:
```
pi.evalf()
```
Out[5]:
```
3.14159265358979
```
In [6]:
```
(pi + exp(1)).evalf()
```
Out[6]:
```
5.85987448204884
```
如你所见，将表达式评估为浮点数。
也有一个类代表数学的无限, 称为 oo:
In [7]:
```
oo > 99999
```
Out[7]:
```
True
```
In [8]:
```
oo + 1
```
Out[8]:
```
oo
```
### 3.2.1.2 练习
* 计算 $\sqrt{2}$ 小数点后一百位。
* 用有理数算术计算1/2 + 1/3 in rational arithmetic.
### 3.2.1.3 符号
与其他计算机代数系统不同，在SymPy你需要显性声明符号变量:
In [4]:
```
from sympy import *
x = Symbol('x')
y = Symbol('y')
```
然后你可以计算他们:
In [10]:
```
x + y + x - y
```
Out[10]:
```
2*x
```
In [11]:
```
(x + y)**2
```
Out[11]:
```
(x + y)**2
```
符号可以使用一些Python操作符操作: +, -, *, ** (算术), &, |, ~ , >>, << (布尔逻辑).
**打印** 这里我们使用下列设置打印
In [ ]:
```
sympy.init_printing(use_unicode=False, wrap_line=True)
```
## 3.2.2 代数运算
SymPy可以进行强大的代数运算。我们将看一下最常使用的：展开和化简。
### 3.2.2.1 展开
使用这个模块展开代数表达式。它将试着密集的乘方和相乘:
In [13]:
```
expand((x + y)**3)
```
Out[13]:
```
x**3 + 3*x**2*y + 3*x*y**2 + y**3
```
In [14]:
```
3*x*y**2 + 3*y*x**2 + x**3 + y**3
```
Out[14]:
```
x**3 + 3*x**2*y + 3*x*y**2 + y**3
```
可以通过关键词的形式使用更多的选项:
In [15]:
```
expand(x + y, complex=True)
```
Out[15]:
```
re(x) + re(y) + I*im(x) + I*im(y)
```
In [16]:
```
I*im(x) + I*im(y) + re(x) + re(y)
```
Out[16]:
```
re(x) + re(y) + I*im(x) + I*im(y)
```
In [17]:
```
expand(cos(x + y), trig=True)
```
Out[17]:
```
-sin(x)*sin(y) + cos(x)*cos(y)
```
In [18]:
```
cos(x)*cos(y) - sin(x)*sin(y)
```
Out[18]:
```
-sin(x)*sin(y) + cos(x)*cos(y)
```
## 3.2.2.2 化简
如果可以将表达式转化为更简单的形式，可以使用化简:
In [19]:
```
simplify((x + x*y) / x)
```
Out[19]:
```
y + 1
```
化简是一个模糊的术语，更准确的词应该是：powsimp (指数化简)、 trigsimp (三角表达式)、logcombine、radsimp一起。
**练习**
* 计算$(x+y)^6$的展开。
* 化简三角表达式$ \sin(x) / \cos(x)$
## 3.2.3 微积分
### 3.2.3.1 极限
在SymPy中使用极限很简单，允许语法limit(function, variable, point), 因此要计算f(x)类似$x \rightarrow 0$, 你应该使用limit(f, x, 0):
In [5]:
```
limit(sin(x)/x, x, 0)
```
Out[5]:
```
1
```
你也可以计算一下在无限时候的极限:
In [6]:
```
limit(x, x, oo)
```
Out[6]:
```
oo
```
In [7]:
```
limit(1/x, x, oo)
```
Out[7]:
```
0
```
In [8]:
```
limit(x**x, x, 0)
```
Out[8]:
```
1
```
### 3.2.3.2 微分法
你可以使用`diff(func, var)`微分任何SymPy表达式。例如:
In [9]:
```
diff(sin(x), x)
```
Out[9]:
```
cos(x)
```
In [10]:
```
diff(sin(2*x), x)
```
Out[10]:
```
2*cos(2*x)
```
In [11]:
```
diff(tan(x), x)
```
Out[11]:
```
tan(x)**2 + 1
```
你可以用下列方法检查是否正确:
In [12]:
```
limit((tan(x+y) - tan(x))/y, y, 0)
```
Out[12]:
```
tan(x)**2 + 1
```
可以用`diff(func, var, n)`方法来计算更高的导数:
In [13]:
```
diff(sin(2*x), x, 1)
```
Out[13]:
```
2*cos(2*x)
```
In [14]:
```
diff(sin(2*x), x, 2)
```
Out[14]:
```
-4*sin(2*x)
```
In [15]:
```
diff(sin(2*x), x, 3)
```
Out[15]:
```
-8*cos(2*x)
```
### 3.2.3.3 序列展开
SymPy也知道如何计算一个表达式在一个点的Taylor序列。使用`series(expr, var)`:
In [16]:
```
series(cos(x), x)
```
Out[16]:
```
1 - x**2/2 + x**4/24 + O(x**6)
```
In [17]:
```
series(1/cos(x), x)
```
Out[17]:
```
1 + x**2/2 + 5*x**4/24 + O(x**6)
```
**练习**
计算$\lim_{x\rightarrow 0} \sin(x)/x$
计算`log(x)`对于x的导数。
### 3.2.3.4 积分法
SymPy支持超验基础和特殊函数的无限和有限积分，通过`integrate()` 功能, 使用了强大的扩展的Risch-Norman算法和启发式和模式匹配。你可以积分基本函数:
In [18]:
```
integrate(6*x**5, x)
```
Out[18]:
```
x**6
```
In [19]:
```
integrate(sin(x), x)
```
Out[19]:
```
-cos(x)
```
In [20]:
```
integrate(log(x), x)
```
Out[20]:
```
x*log(x) - x
```
In [21]:
```
integrate(2*x + sinh(x), x)
```
Out[21]:
```
x**2 + cosh(x)
```
也可以很简单的处理特殊函数:
In [22]:
```
integrate(exp(-x**2)*erf(x), x)
```
Out[22]:
```
sqrt(pi)*erf(x)**2/4
```
也可以计算一下有限积分:
In [23]:
```
integrate(x**3, (x, -1, 1))
```
Out[23]:
```
0
```
In [24]:
```
integrate(sin(x), (x, 0, pi/2))
```
Out[24]:
```
1
```
In [25]:
```
integrate(cos(x), (x, -pi/2, pi/2))
```
Out[25]:
```
2
```
不标准积分也支持:
In [26]:
```
integrate(exp(-x), (x, 0, oo))
```
Out[26]:
```
1
```
In [27]:
```
integrate(exp(-x**2), (x, -oo, oo))
```
Out[27]:
```
sqrt(pi)
```
#### 3.2.3.5 练习
### 3.2.4 方程求解
SymPy可以求解线性代数方程，一个或多个变量:
In [28]:
```
solve(x**4 - 1, x)
```
Out[28]:
```
[-1, 1, -I, I]
```
如你所见，第一个参数是假设等于0的表达式。它可以解一个很大的多项式方程，也可以有能力求解多个方程，可以将各自的多个变量作为元组以第二个参数给出:
In [29]:
```
solve([x + 5*y - 2, -3*x + 6*y - 15], [x, y])
```
Out[29]:
```
{x: -3, y: 1}
```
也直接求解超越方程（有限的）:
In [30]:
```
solve(exp(x) + 1, x)
```
Out[30]:
```
[I*pi]
```
多项式方程的另一个应用是`factor`。`factor`将多项式因式分解为可化简的项，并且可以计算不同域的因式:
In [31]:
```
f = x**4 - 3*x**2 + 1
factor(f)
```
Out[31]:
```
(x**2 - x - 1)*(x**2 + x - 1)
```
In [32]:
```
factor(f, modulus=5)
```
Out[32]:
```
(x - 2)**2*(x + 2)**2
```
SymPy也可以解布尔方程，即，判断一个布尔表达式是否满足。对于这个情况，我们可以使用`satisfiable`函数:
In [33]:
```
satisfiable(x & y)
```
Out[33]:
```
{x: True, y: True}
```
这告诉我们`(x & y)`是真，当x和y都是True的时候。如果一个表达式不是True，即它的任何参数值都无法使表达式为真，那么它将返回False:
In [34]:
```
satisfiable(x & ~x)
```
Out[34]:
```
False
```
### 3.2.4.1 练习
* 求解系统方程$x + y = 2$, $2\cdot x + y = 0$
* 是否存在布尔值，使$(~x | y) & (~y | x)$为真?
### 3.2.5 线性代数
#### 3.2.5.1 矩阵
矩阵通过Matrix类的一个实例来创建:
In [35]:
```
from sympy import Matrix
Matrix([[1,0], [0,1]])
```
Out[35]:
```
Matrix([
[1, 0],
[0, 1]])
```
与NumPy数组不同，你也可以在里面放入符号:
In [36]:
```
x = Symbol('x')
y = Symbol('y')
A = Matrix([[1,x], [y,1]])
A
```
Out[36]:
```
Matrix([
[1, x],
[y, 1]])
```
In [37]:
```
A**2
```
Out[37]:
```
Matrix([
[x*y + 1, 2*x],
[ 2*y, x*y + 1]])
```
### 3.2.5.2 微分方程
SymPy可以解 (一些) 常规微分。要求解一个微分方程，使用`dsolve`。首先，通过传递cls=Function来创建一个未定义的符号函数:
In [38]:
```
f, g = symbols('f g', cls=Function)
```
f 和 g是未定义函数。我们可以调用f(x), 并且它可以代表未知的函数:
In [39]:
```
f(x)
```
Out[39]:
```
f(x)
```
In [40]:
```
f(x).diff(x, x) + f(x)
```
Out[40]:
```
f(x) + Derivative(f(x), x, x)
```
In [41]:
```
dsolve(f(x).diff(x, x) + f(x), f(x))
```
Out[41]:
```
f(x) == C1*sin(x) + C2*cos(x)
```
关键词参数可以向这个函数传递，以便帮助确认是否找到最适合的解决系统。例如，你知道它是独立的方程，你可以使用关键词hint=’separable’来强制`dsolve`来将它作为独立方程来求解:
In [42]:
```
dsolve(sin(x)*cos(f(x)) + cos(x)*sin(f(x))*f(x).diff(x), f(x), hint='separable')
```
Out[42]:
```
[f(x) == -asin(sqrt(C1/cos(x)**2 + 1)) + pi,
f(x) == asin(sqrt(C1/cos(x)**2 + 1)) + pi,
f(x) == -asin(sqrt(C1/cos(x)**2 + 1)),
f(x) == asin(sqrt(C1/cos(x)**2 + 1))]
```
**练习**
* 求解Bernoulli微分方程
$x \frac{d f(x)}{x} + f(x) - f(x)^2=0$
* 使用hint=’Bernoulli’求解相同的公式。可以观察到什么?