前置芝士#

会Python就行。数学水平因人而异，根据自己的需要选用相应的功能。

什么是符号计算？#

前面提到，SymPy是一个 符号计算 库，那么什么是符号计算。符号计算（Symbolic Computation），就是在计算机上用符号精确表示数学对象，而不是得到一个近似的数值解。

看一个例子：

1
>>>import math
2
>>>math.sqrt(8)
3
2.8284271247461903

Python的math库给出的计算$\sqrt 8$的结果是$2.8284271247461903$，然而实际上，这只是$\sqrt 8$的一个近似值，我们知道，它的精确值是不能用有限的小数表示的。

而在SymPy中：

1
>>> from sympy import *
2
>>> sqrt(8)
3
2*sqrt(2)

SymPy给我们的结果是$\sqrt 8 = 2 \sqrt 2$，这才是一个精确的化简的结果。在SymPy中，非完全平方数的平方根默认会被保留。

定义符号与计算#

在SymPy中，符号需要在使用前用Symbol()或symbols()来定义：

1
a = Symbol('a')
2
x, y = symbols('x y')
3
# symbols接受一系列用空格或逗号隔开的符号，返回对应的变量

最好将变量名设置得与符号名一致，方便记忆和使用。除非符号中包含了Python不允许的符号，或者符号名太长，想使用一个短一点的变量名来表示。

SymPy支持直接使用+、-、*、/、**来操作这些符号。

1
>>> x, y = symbols("x y")
2
>>> expr = x + 2*y
3
>>> expr
4
x + 2*y
5
>>> expr + 1 - x
6
2*y + 1
7
>>> expr / 2
8
x/2 + y

输入x * expr你会发现SymPy并没有为我们展开表达式为x**2 + 2*x*y，而是保留了x*(x + 2*y)的形式。因为因式分解是SymPy默认的化简操作，你也可以利用expand来展开，用factor来分解。

1
>>> expand(x * expr)
2
x**2 + 2*x*y
3
>>> factor(2*x**2 - 7*x - 4)
4
(x - 4)*(2*x + 1)

SymPy计算出了展开$x(x + 2y)$的结果$x^2 + 2xy$，和因式分解$2x^2 - 7x - 4$的结果$(x-4)(2x+1)$。

输出#

可以通过调整一些选项改变SymPy输出的格式。

用init_printing(use_unicode=True)可以允许SymPy使用Unicode输出。

计算$\int^{+\infty}_{-\infty}\sin(x^2)dx$

1
>>> init_printing(use_unicode=False)
2
>>> integrate(sin(x**2), (x, -oo, oo))
3
  ___   ____
4
\/ 2 *\/ pi
5
------------
6
     2
7
>>> init_printing(use_unicode=True)
8
>>> integrate(sin(x**2), (x, -oo, oo))
9
√2⋅√π
10
─────
11
  2

第一个输出中SymPy使用了字符画的方式画出 $\sqrt{}$，使用了pi表示$\pi$。第二个输出使用了Unicode字符“√”和“π”。

SymPy还可以使用 $\LaTeX$:

1
>>> latex(Integral(cos(x)**2, (x, 0, pi)))
2
\int\limits_{0}^{\pi} \cos^{2}{\left(x \right)}\, dx

将\int\limits_{0}^{\pi} \cos^{2}{\left(x \right)}\, dx使用任意一个 $\LaTeX$ 编译器编译，就可以得到下面这样漂亮的结果：

符号的替换#

用subs()，但注意subs()并不改变原表达式。

1
>>> expr = cos(x) + 1
2
>>> expr.subs(x, y)
3
cos(y) + 1
4
>>> expr
5
cos(x) + 1

也可以替换数字，进行在某一点求值的操作。

1
>>> expr = sin(2*x) + cos(2*x)
2
>>> expr.subs(x, pi/5)
3
-1/4 + sqrt(5)/4 + sqrt(sqrt(5)/8 + 5/8)

也可以替换多个符号

1
>>> expr = x**3 + 4*x*y - z
2
>>> expr.subs([(x, 2), (y, 4), (z, 0)])
3
40

字符串转表达式#

用sympify()（不是symplify（）

1
>>> str_expr = "x**2 + 3*x - 1/2"
2
>>> expr = sympify(str_expr)
3
>>> expr
4
x**2 + 3*x - 1/2
5
>>> expr.subs(x, 2)
6
19/2

表达式转浮点数#

用evalf()

1
>>> expr = sqrt(8)
2
>>> expr.evalf()
3
2.82842712474619

默认情况下，会给出15位的精度，但是也可以指定所需要的精度

1
expr.evalf(100)
2
2828427124746190097603377448419396157139343750753896146353359475981464956924214077700775068655283145

如果想要用subs()进行单点求值后再用evalf()进行数值估计，用expr.evalf(subs={<Symbol>: <Value>})是比expr.subs(<Symbol>, <Value>).evalf()更推荐的方法。因为这样会使计算更高效且稳定。

1
>>> expr = cos(2*x)
2
>>> expr.evalf(subs={x: 2.4})
3
0.0874989834394464

假如指定了chop=True，SymPy会自动移除小于期望精度的舍入误差

1
>>> one = cos(1)**2 + sin(1)**2
2
>>> (one - 1).evalf()
3
-0.e-124
4
>>> (one - 1).evalf(chop=True)
5
0

多点求值#

使用SymPy对表达式做多点的求值（比如成千上万点后）是比较慢的。如果有这种需求，可以对表达式lambdify()后交给其他库处理（比如NumPy或SciPy）

1
>>> import numpy
2
>>> a = numpy.arange(10)
3
>>> expr = sin(x)
4
>>> f = lambdify(x, expr, "numpy")
5
>>> f(a)
6
array([ 0.        ,  0.84147098,  0.90929743,  0.14112001, -0.7568025 ,
7
       -0.95892427, -0.2794155 ,  0.6569866 ,  0.98935825,  0.41211849])

lambdify()将一个SymPy表达式转换成一个计算表达式在某点的值的函数，以便于使用其他库进行计算。

坑一：符号与变量#

符号（Symbol）指的是SymPy提供的，用symbols()定义的符号，而变量是Python提供的。

1
>>> x = symbols('x')
2
>>> expr = x + 1
3
>>> x = 2
4
>>> print(expr)
5
x + 1

以上这个程序，并没有输出3，而是输出了x + 1，这是因为x = 2只是将变量x的值，从符号$x$，变成了数字2，而符号$x$并没有发生改变，因此expr的值也就没有改变。

想要计算expr在x=2条件下的值，可以使用subs:

1
>>> x = symbols('x')
2
>>> expr = x + 1
3
>>> expr.subs(x, 2)
4
3

坑二：等于号#

假如你用==检测两个式子是否相等，你将会得到错误的结果：

1
>>> (x + 1)**2 == x**2 + 2*x + 1
2
False

SymPy并没有扩展Python的语法，=依然表示赋值，==表示变量的相等，因此在SymPy中，主要用这两种方式表示两个式子相等：

1. 这是官方文档中推荐的方法：假如你想要验证$a = b$是否成立，可以检测$a - b = 0$是否成立，只要将两个式子相减，然后使用simplify()函数化简，如果结果是0，则必定相等，如果不是0，大多数情况下不相等，但也存在极少数情况，式子确实为0，但是SymPy无法将其化简，对于常见的数学式子，可以放心地使用这个方法检验其是否相等。

   1
   >>> a = (x + 1)**2
   2
   >>> b = x**2 + 2*x + 1
   3
   >>> simplify(a - b)
   4
   0
   5
   >>> c = x**2 - 2*x + 1
   6
   >>> simplify(a - c)
   7
   4*x
   

2. 另一种方法是使用a.equals(b)，但是这个方法并不是使用式子的化简得到的结论，而是通过随机地给式子赋多个值，比较二者是否相等。

   1
   >>> a = cos(x)**2 - sin(x)**2
   2
   >>> b = cos(2*x)
   3
   >>> a.equals(b)
   4
   True
   

坑三：^符号#

SymPy使用和Python一样的符号约定，用**表示指数，^表示异或（Xor），所以不应该尝试使用^表示指数。

1
>>> Xor(x, y)
2
x ^ y

坑四：/符号#

两个SymPy对象相除，或者一个SymPy对象相除的时候，返回值都是SymPy对象，这没有问题。但是当分子分母都是Python的int时，要小心了。两个SymPy对象相除会返回一个分数，但是两个Python的int相除会返回一个浮点数（Python3）。

1
>>> x / y
2
x/y
3
>>> x / 19
4
x/19
5
>>> 3 / x
6
3/x
7
>>> 3 / 19
8
0.15789473684210525

0.15789473684210525不是我们想要的结果，这个值不是精确的，当处理两个整数相除的时候，可以使用Rational()来避免这种结果：

1
Rational(3, 19)
2
3/19

如果和其他符号连接起来的时候这种错误就更难以察觉：

1
>>> x + 1/3
2
x + 0.333333333333333
3
>>> x + Rational(1, 3)
4
x + 1/3

初始化#

init_printing()会为你启用当前环境下可用的最佳的输出效果。

1
from sympy import init_printing
2
init_printing()

如果你使用交互式命令行，init_session() 函数将自动导入SymPy中的所有内容，创建一些常用符号，设置绘图，并运行 init_printing()。

1
>>> from sympy import init_session
2
>>> init_session()
3
Python console for SymPy 0.7.3 (Python 2.7.5-64-bit) (ground types: gmpy)
4

5
These commands were executed:
6
>>> from __future__ import division
7
>>> from sympy import *
8
>>> x, y, z, t = symbols('x y z t')
9
>>> k, m, n = symbols('k m n', integer=True)
10
>>> f, g, h = symbols('f g h', cls=Function)
11
>>> init_printing() # doctest: +SKIP
12

13
Documentation can be found at https://www.sympy.org/

用Jupyter美化输出#

但是终端的输出效果实在是有点丑

这是设置了init_printing(use_unicode=False)，不允许使用Unicode字符的输出

1
>>> Integral(sqrt(1/x), x)
2
  /
3
 |
4
 |     ___
5
 |    / 1
6
 |   /  -  dx
7
 | \/   x
8
 |
9
/

这是设置了init_printing(use_unicode=True)，允许使用Unicode字符的输出

1
>>> Integral(sqrt(1/x), x)
2
⌠
3
⎮     ___
4
⎮    ╱ 1
5
⎮   ╱  ─  dx
6
⎮ ╲╱   x
7
⌡

这是Jupyter notebook里的输出，使用了MathML渲染的$\LaTeX$：

NOTE

确保你安装好了Jupyter Notebook，在Jupyter中才能达成上面的效果

输出了很漂亮的公式，以后都可以这样使用。

一些常用的化简函数#

1. expand()展开，如

   $$(x + 2)*(x - 3) = x^2 - x - 6$$

2. factor()因式分解，如

   $$x^{2} z + 4 x y z + 4 y^{2} z = z \left(x + 2 y\right)^{2}$$

   如果希望得到所有因子的列表，用factor_list()

3. collect(expr, x)将x作为主元整理expr。如

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

4. cancel()对分子分母自动消去共同因子，如

   $$\frac{x^{2} + 2 x + 1}{x^{2} + x} = \frac{x + 1}{x}$$

5. apart()会执行部分因式分解（Partial fraction decomposition）

6. trigsimp()是三角函数版本的simplify()会智能地化简三角函数表达式。

7. expand_trig()会使用和差角、二倍角等公式展开表达式，如

   $$\sin{\left(x + y \right)} + \tan{\left(2 x \right)} = \sin{\left(x \right)} \cos{\left(y \right)} + \sin{\left(y \right)} \cos{\left(x \right)} + \frac{2 \tan{\left(x \right)}}{1 - \tan^{2}{\left(x \right)}}$$

8. expr.rewrite(function)尝试用function表示expr，如

   

SymPy中的化简函数还有很多，这里只是一些常用的。

求导#

diff(expr, x)可以求expr关于x的导数，也可以写成expr.diff(x)。

diff(expr, x，x, x)可以对x求三阶导，也可以写成diff(expr, x，3)或expr.diff(x, 3)。

Derivative()可以创建一个导数但不计算它

对未计算的导数使用.doit()可以计算它。

积分#

integrate(exp(-x), (x, 0, oo))将会计算

其中oo是两个“o”，用来表示无穷大。

integrate(exp(-x**2 - y**2), (x, -oo, oo), (y, -oo, oo))将会计算

和上面求导一样，Integral会创建一个积分但是不计算它，.doit()可以计算未计算的积分。

极限#

limit(sin(x)/x, x, 0)将会计算

计算单侧极限，可以向第三个参数传入'+'或'-'。limit(1/x, x, 0, '-')将会计算

同上，Integral会创建一个未计算的极限，.doit()可以计算它。

计算某点处的级数展开#

用expr.series(x, a, n)会给出表达式在 x = a 处的 n 阶展开。

不想要那个表示余项的O，就用.removeO()

求解方程#

在前面已经提到过，SymPy中表达相等应该用Eq()而不是==或=

用solveset(equation, x)可以求equation关于x的解集，如果等式右边是0，也可以不写成solveset(Eq(expr, 0), x)，直接写成solveset(expr, x)，例如：

1
>>> solveset(Eq(x**2, 1), x)
2
{-1, 1}
3
>>> solveset(x**2 - 1, x)
4
{-1, 1}

如果无解，将返回 $\varnothing$ ，如果无法求解，会返回一个条件集合。

1
>>> solveset(exp(x), x)     # 无解
2
∅
3
>>> solveset(cos(x) - x, x)  # 无法求解
4
{x │ x ∊ ℂ ∧ (-x + cos(x) = 0)}

求解线性方程组#

使用linsolve()

1
>>> linsolve([x + y + z - 1, x + y + 2*z - 3 ], (x, y, z)) # 等式列表形式
2
{(-y - 1, y, 2)}
3
>>> linsolve(Matrix(([1, 1, 1, 1], [1, 1, 2, 3])), (x, y, z)) # 增广矩阵形式
4
{(-y - 1, y, 2)}
5
>>> M = Matrix(((1, 1, 1, 1), (1, 1, 2, 3))) # A*x = b 形式
6
system = A, b = M[:, :-1], M[:, -1]
7
linsolve(system, x, y, z)
8
{(-y - 1, y, 2)}

求解非线性方程组#

使用nonlinsolve()

获得多项式的重根数#

solveset()返回的都是不重复的根，要获得多项式的重根数，可以使用roots()

1
>>> solveset(x**3 - 6*x**2 + 9*x, x)
2
{0, 3}
3
>>> roots(x**3 - 6*x**2 + 9*x, x)
4
{0: 1, 3: 2}

求解微分方程#

使用dsolve()

基本操作#

+ * **，分别代表加法，乘法和乘方，矩阵求逆就是M**-1

M.T返回矩阵的转置。

eye(n)、ones(n, m)、zeros(n, m)分别可以用来创建单位矩阵、全1矩阵和零矩阵。

diag()可以创建对角矩阵

1
>>> diag(-1, ones(2, 2), Matrix([5, 7, 5]))
2
⎡-1  0  0  0⎤
3
⎢           ⎥
4
⎢0   1  1  0⎥
5
⎢           ⎥
6
⎢0   1  1  0⎥
7
⎢           ⎥
8
⎢0   0  0  5⎥
9
⎢           ⎥
10
⎢0   0  0  7⎥
11
⎢           ⎥
12
⎣0   0  0  5⎦

其他操作#

• 行列式：.det()

• 转化为简化行阶梯矩阵：.rref()，返回一个元组，第一个元素是简化行阶梯矩阵，第二个元素是一个存放了阶梯头所在列的列标的元组。

  

• 子空间：.nullspace() 零空间、 .columnspace()列空间、 rowspace()行空间。

• .eigenvals()特征值、.eigenvects()特征向量、.charpoly(lamda)特征多项式、.diagonalize()对角化。

一个实例#

矩阵$A=\left[ \begin{array}{llll} 1 & 2 & 3 & 4 \\ 2 & -1 & -2 & a^2 \\ -1 & -7 & -11 & a \end{array}\right]$是一个线性方程组的增广矩阵，问$a$取多少的时候，该方程组有解？

输入以下代码，即可解出，当 $a = -3$ 或 $a = 4$ 时方程组有解。

1
a, x1, x2, x3 = symbols('a x1 x2 x3')
2
A = Matrix([
3
    [1, 2, 3],
4
    [2, -1, -2],
5
    [-1, -7, -11]
6
])
7
b = Matrix([4, a**2, a])
8
x = Matrix([x1, x2, x3])
9
solve(A*x - b)