第 6 章　用Python做符号计算

第 6 章　用Python做符号计算

SymPy包括许多功能，从基本的符号算术到多项式、微积分、求解方程、离散数学、几何、统计和物理。它主要处理三种类型的数据：整型数据、实数和有理数。整数即不带小数点的数字，而实数是带小数点的数字。有理数包括两个部分：分子和分母。可以用Ration类定义有理数，该类需要两个数字。在这一章中，我们将通过示例程序来理解SymPy的概念。

这一章将学习以下主题：

用SymPy的计算机代数系统

核心功能和高级功能

多项式、微积分和求解方程

离散数学、矩阵、几何、画图、物理和统计

打印功能

首先讨论SymPy的核心功能，包括基本的算术、扩展、简化、替换、模式匹配和各种函数（例如指数函数、对数函数、方程的根、三角函数、双曲线函数和特殊函数）。

6.1　符号、表达式和基本运算

在SymPy中，在任何表达式中使用符号前，必须先定义该符号。定义符号很简单。只需要用Symbol类中的symbol函数来定义一个符号即可，下面的程序就用到了这种方式。我们可以用evalf()/n()来获得任何对象的浮点近似值。

下面的程序用了三种方法来生成符号。对于只生成一个符号的，函数命名为symbol；对于生成多个符号的，函数命名为symbols。有两种方法可以生成多个符号：第一种是将用空格分隔的符号名称传入符号函数；第二种方法是将m0:5传入符号函数，生成一个如m0, m1, m2, m3, m4的符号序列。在第二种方法中，m0是第一个索引值，冒号后的5表示生成五个这样的变量。

通常，两个整数的除会删除小数点后面的部分。为了避免这种情况，下面程序的第一行强制了两个整数的除的结果为浮点类型。这就是为什么程序的最后一行将3.142857142857143存入y中。如果忽略程序的第一行，那么y的值是3：

from __future__ import division from sympy import *
 
x, y, z = symbols('x y z')
m0, m1, m2, m3, m4 = symbols('m0:5')
x1 = Symbol('x1')
x1 + 500
y=22/7

下面的程序用evalf()和n()函数将任何SymPy对象的数值近似转换为浮点值。默认的精度是小数点后15位，而且可以通过调整参数改成任何想要的精度：

from __future__ import division from sympy import, sin, pi
 
x=sin(50)
 
pi.evalf()
pi.evalf(50)
 
x.n()
x.n(20)

下面的程序展示了表达式的概念，以及用collect、expand、factor、simplify和subs进行的各种操作：

from sympy import collect, expand, factor, simplify from sympy import Symbol, symbols
from sympy import sin, cos
 
x, y, a, b, c, d = symbols('x y a b c d')
 
expr = 5x*2+2*bx*2+cos(x)+51x*2
simplify(expr)
 
factor(x**2+x-30)
expand ( (x-5) (x+6) )
 
collect(x*3 + ax*2 + bx*2 + c*x + d, x) 
expr = sin(x)*sin(x) + cos(x)*cos(x)
expr
expr.subs({x:5, y:25})
expr.subs({x:5, y:25}).n()

6.2　求解方程

solve是一个神奇的函数。它可以求解各种类型的方程。这个函数返回方程的解。它需要两个输入参数：待解的表达式和变量。下面的程序用这个函数求解了各种方程。在以下的方程中，假设方程等式右边的值是0：

from sympy import solve, symbols 
a, b, c, x, y = symbols('a b c x y')
 
solve (6x*2 - 3*x - 30,x)
solve(ax*2 + b*x + c, x)
substitute_solution = solve(ax*2 + b*x + c, x) [ substitute_solution[0].subs({'a':6,'b':-3,'c':-30}), substitute_solution[1].subs({'a':6,'b':-3,'c':-30}) ]
 
solve([2*x + 3*y - 3, x -2*y + 1], [x, y]
)

为了求解方程组，我们有另一种形式的solve函数。它将一系列的方程作为第一个输入参数，将未知数列表作为第二个参数，如下面的程序所示：

from sympy import solve, symbols 
x, y = symbols('x y')
solve ([2*x + y - 4, 5*x - 3*y],[x, y])
solve ([2*x + 2*y - 1, 2*x - 4*y],[x, y])

6.3　有理数、指数和对数函数

SymPy有很多函数可以用于处理有理数。这些函数可以对有理数做简化、扩展、合并、拆分等操作。SymPy还支持一些指数和对数操作。有三个对数函数：log（用于计算以b为底的对数）、ln（用于计算以欧拉常数e为底的自然对数）和log10（用于以10为底的对数）。log函数需要两个参数：变量和底。如果底没有被指定，那么这个函数会默认以自然对数计算该变量，也就是等价于ln。为了计算两个有理数的加法，我们用together函数。类似地，计算有理数的除法，我们用到apart函数，如下面的程序所示：

from sympy import together, apart, symbols, exp, log, ln import mpmath
x1, x2, x3, x4, x = symbols('x1 x2 x3 x4 x')
x1/x2 + x3/x4
together(x1/x2 + x3/x4)
 
apart ((2x*2+3*x+4)/(x+1))
together(apart ((2x*2+3*x+4)/(x+1))) 
exp(1)
log(4).n()
log(4,4).n()
ln(4).n()
mpmath.log10(4)

6.4　多项式

SymPy也允许我们定义和执行多项式的各种操作。我们还可以求解多项式的根。我们已经介绍过simplify、expand、factor和solve函数。这些函数也适用于多项式。为了检查两个多项式是否相等，应该用simplify函数：

from sympy import 
p, q, x = symbols ('p q x')
p = (x+4)(x+2)
q = x**2 + 6*x + 8
p == q # 返回False
p - q == 0 # 返回False
simplify(p - q) == 0

6.5　三角函数和复数

大部分时候，三角函数的输入是弧度值，反三角函数返回的也是弧度值。这个模块也提供了从度数到弧度值的转换和从弧度值到度数的转换。除了sin、cos和tan这样的基本三角函数，SymPy还有简化和扩展的三角函数。

SymPy还支持在不存在实数解的情况下的复数解。例如等式x*2+4=0，这个方程没有实数解，它的解是-2I 或+2*I（I 表示-1的平方根）。下面的程序展示了三角函数，并且用复数值形式给出了这个方程的解：

from sympy import 
x, y = symbols('x y')
expr = sin(x)cos(y)+cos(x)*sin(y)
expr_exp= exp(5*sin(x)**2+4*cos(x)**2)
 
trigsimp(expr)
trigsimp(expr_exp)
expand_trig(sin(x+y))
solve(x**2+4,x) #结果是复数

6.6　线性代数

SymPy线性代数模块是另外一个非常简单的模块，该模块中的函数的矩阵操作都非常简单易学。它包括对矩阵的各种操作，例如快速特殊矩阵的构建、特征值、特征向量、转置、行列式的值和逆。有三个函数可以快速生成特殊矩阵，分别是eye、zeros和ones。eye函数生成一个实体矩阵，而zeros和ones生成的矩阵的全部元素分别是0和1。如果需要的话，我们可以删除矩阵中某些选中的行和列。基本算术运算，如+、-、*和**，也可以用于矩阵：

from sympy import 
A = Matrix( [[1, 2, 3, 4],
            [5, 6, 7, 8],
            [ 9, 10, 11, 12],
            [ 13, 14, 15, 16]] )
A.row_del(3)
A.col_del(3)
 
A[0,1] # 显示矩阵A的第一行第二列
A[0:2,0:3] # 左上角的2×3子矩阵
 
B = Matrix ([[1, 2, 3],
            [5, 6, 7],
            [ 9, 10, 11]] )
A.row_join(B)
B.col_join(B)
A + B
A - B
A B
A **2
eye(3) # 单位矩阵（identity matrix）
zeros(3, 3) # 3×3全0矩阵
ones(3, 3) # 3×3全1矩阵
 
A.transpose() # 等价于A.T
M = Matrix( [[1, 2, 3],
            [4, 5, 6],
            [7, 8, 10]] )
M.det()

默认情况下，矩阵的逆是用高斯消元法计算得出的，也可以指定用LU分解法来计算。SymPy有很多函数可以用来计算简化的行阶梯形（rref函数）和零空间（nullspace函数）。如果A是一个矩阵，那么nullspace就是所有向量v的集合，满足A v=0。还可以对矩阵元素做替换操作。我们可以用符号实体生成一个矩阵，然后用实际的值和其他符号来做替换。还可以做一些特殊的操作，如QR因式分解、Gram-Schmidt正交化和LU分解：

from sympy import 
A = Matrix( [[1,2],
            [3,4]] )
A.inv()
A.inv()A
A*A.inv()
A = Matrix( [[ 1, -2],
            [-2, 3]] )
A.eigenvals() # 等价于solve( det(A-eye(2)*x), x)
A.eigenvects()
A.rref()
 
A.nullspace()
 
x = Symbol('x')
M = eye(3) *x
M.subs(x, 4)
y = Symbol('y')
M.subs(x, y)
 
M.inv()
M.inv("LU")
 
A = Matrix([[1,2,1],[2,3,3],[1,3,2]])
Q, R = A.QRdecomposition()
Q
 
M = [Matrix([1,2,3]), Matrix([3,4,5]), Matrix([5,7,8])]
result1 = GramSchmidt(M)
result2 = GramSchmidt(M, True)

6.7　微积分

微积分涉及用于研究任何函数的各种性质的操作。这些操作包括变化率、函数的极限以及函数图像在一段区间的积分。在这一节中，你将学到极限、导数、序列的求和以及积分。下面的程序用到极限函数来解决简单的极限问题：

from sympy import limit, oo, symbols, exp, sin, cos 
oo+5
67000 < oo
10/oo
 
x , n = symbols ('x n')
limit( ((x**n - 1)/ (x - 1) ), x, 1)
 
limit( 1/x**2, x, 0)
limit( 1/x, x, 0, dir="-")
 
limit(cos(x)/x, x, 0)
limit(sin(x)**2/x, x, 0)
limit(exp(x)/x,x,oo)

任何SymPy表达式都可以用diff函数进行微分，diff函数的原型是diff(func_to_be_differentiated, variable)。下面的程序用到了diff函数来计算各种SymPy表达式的微分：

from sympy import diff, symbols, Symbol, exp, dsolve, subs, Function from sympy import 
x = symbols('x')
 
diff(x*4, x)
diff( x**3*cos(x), x )
diff( cos(x**2), x )
diff( x/sin(x), x )
diff(x**3, x, 2)
diff( exp(x), x)

dsolve函数帮助我们求解任何常微分方程和常微分方程组。下面的程序展示了用dsolve函数求解常微分方程组和边界值的问题：

x = symbols('x') f = symbols('f', cls=Function)
dsolve( f(x) - diff(f(x),x), f(x) )
 
#ics参数用Python字典确定dsolve方法的边界
from sympy import 
x=symbols('x')
f=symbols('f', cls=Function)
dsolve(Eq(f(x).diff(x,x), 400(f(x)-1+2*x)), f(x), ics={f(0):5, f(2):10}) # 上面这行代码就是计算f(x) = C1*e^-30x + C2*e^30x - 2x + 1的结果

下面的程序用于寻找函数f (x)=4x³-3x²+2x的临界点，并用二次导数来找该函数在[0,1]区间的最大值：

x = Symbol('x') f = 4x*3-3x*2+2*x diff(f, x)
sols = solve( diff(f,x), x)
sols
diff(diff(f,x), x).subs( {x:sols[0]} )
diff(diff(f,x), x).subs( {x:sols[1]} )

在SymPy中，可以用integrate函数计算定积分和不定积分。下面的程序计算了定积分和不定积分，并用符号定义这些积分。为了计算实际值，在下面的程序的最后一行调用n()计算积分：

from sympy import 
x = symbols('x')
 
integrate(x*3+1, x)
integrate(x*sin(x), x)
integrate(x+ln(x), x)
 
F = integrate(x**3+1, x)
F.subs({x:1}) - F.subs({x:0})
 
integrate(x**3-x**2+x, (x,0,1))   # 定积分
integrate(sin(x)/x, (x,0,pi))    # 定积分
integrate(sin(x)/x, (x,pi,2*pi))  # 定积分
integrate(x*sin(x)/(x+1), (x,0,2*pi))  # 定积分
integrate(x*sin(x)/(x+1), (x,0,2*pi)).n()

序列是将整数作为输入的函数，可以通过为其第n个值指定一个表达式来定义一个序列。还可以替换其中的某些值。下面的程序用到了SymPy中的一些简单序列来解释序列的概念：

from sympy import *
 
n = symbols('n')
 
s1_n = 1/n
s2_n = 1/factorial(n)
s1_n.subs({n:5})
[ s1_n.subs({n:index1}) for index1 in range(0,8) ]
[ s2_n.subs({n:index1}) for index1 in range(0,8) ]
limit(s1_n, n, oo)
limit(s2_n, n, oo)

一个包含变量不同阶指数1的级数称为幂级数，例如泰勒级数、指数级数和sin/cos级数。下面的程序用一些特殊函数计算了一些级数。它还用到了幂级数的概念：

1变量 x 的 N 次方。——译者注

from sympy import *
 
x, n = symbols('x n')
 
s1_n = 1/n
s2_n = 1/factorial(n)
summation(s1_n, [n, 1, oo])
summation(s2_n, [n, 0, oo])
import math
 
def s2_nf(n):
    return 1.0/math.factorial(n)
 
sum( [s2_nf(n) for n in range(0,10)] )
E.evalf()
 
exponential_xn = x**n/factorial(n)
summation( exponential_xn.subs({x:5}), [x, 0, oo] ).evalf() exp(5).evalf()
 
summation( exponential_xn.subs({x:5}), [x, 0, oo])
import math # 用Python再做一版
def exponential_xnf(x,n):
    return x**n/math.factorial(n)
 
sum( [exponential_xnf(5.0,i) for i in range(0,35)] )
series( sin(x), x, 0, 8)
series( cos(x), x, 0, 8)
series( sinh(x), x, 0, 8)
series( cosh(x), x, 0, 8)
series(ln(x), x, 1, 6) # ln(x)在x=1时的泰勒级数
series(ln(x+1), x, 0, 6) # ln(x+1)的麦克劳林级数

6.8　向量

一个包括实数值的 n 元组就是一个向量。在物理和数学中，一个向量是一个数学对象，该对象拥有大小、幅度或长度，以及方向。在SymPy中，一个向量表示为一个1 × n 的行矩阵或n × 1的列矩阵。下面的程序展示了SymPy中向量计算的概念。它计算了向量的转置和范数：

from sympy import *
u = Matrix([[1,2,3]]) # 一个行向量等于1×3矩阵
v = Matrix([[4],
            [5],    # 一个行向量等于3×1矩阵
            [6]])
v.T # 转置操作
u[1] # 向量索引从0开始
u.norm() # u的范数
uhat = u/u.norm() # 与u同维度的单位向量
uhat
uhat.norm()

下面的程序展示了点积、叉积和向量的投影操作：

from sympy import *
u = Matrix([ 1,2,3])
v = Matrix([-2,3,3])
u.dot(v)
 
acos(u.dot(v)/(u.norm()*v.norm())).evalf()
u.dot(v) == v.dot(u)
u = Matrix([2,3,4])
n = Matrix([2,2,3])
(u.dot(n) / n.norm()**2)*n # 向量v在向量n上的投影
 
w = (u.dot(n) / n.norm()**2)*n
v = u - (u.dot(n)/n.norm()**2)*n # 等价于u - w
u = Matrix([ 1,2,3])
v = Matrix([-2,3,3])
u.cross(v)
(u.cross(v).norm()/(u.norm()*v.norm())).n()
 
u1,u2,u3 = symbols('u1:4')
v1,v2,v3 = symbols('v1:4')
Matrix([u1,u2,u3]).cross(Matrix([v1,v2,v3]))
u.cross(v)
v.cross(u)

6.9　物理模块

物理模块包括了解决物理问题的各种功能。它还包括一些物理子模块，进行向量物理、经典力学、量子力学、光学等操作。

6.9.1　氢波函数

该类型下有两个函数。第一个函数计算Hartree原子单位制中(n, l)态的能量。另一个函数计算Hartree原子单位制中(n, l, spin)态的相对论能量。下面的程序展示了这两个函数的运用：

from sympy.physics.hydrogen import E_nl, E_nl_dirac, R_nl from sympy import var
 
var("n Z")
var("r Z")
var("n l")
E_nl(n, Z)
E_nl(1)
E_nl(2, 4)
 
E_nl(n, l)
E_nl_dirac(5, 2) # l应该小于n
E_nl_dirac(2, 1)
E_nl_dirac(3, 2, False)
R_nl(5, 0, r) # 默认z = 1
R_nl(5, 0, r, 1)

6.9.2　矩阵和Pauli代数

在physics.matrices模块下有一些与物理相关的矩阵。下面的程序展示了如何获得这些矩阵和Pauli代数：

from sympy.physics.paulialgebra import Pauli, evaluate_pauli_product from sympy.physics.matrices import mdft, mgamma, msigma, pat_matrix from sympy import symbols
 
x = symbols('x')
 
mdft(4) # 用离散傅里叶变换做矩阵乘法
mgamma(2) # 狄拉克-伽玛矩阵的狄拉克表示
msigma(2) # Pauli矩阵
 
# 下面这行代码计算平行轴定理矩阵，即质量为m的刚体沿着x轴、y轴、z轴
# 转动(dx, dy, dz)距离的转动惯量
pat_matrix(3, 1, 0, 0)
 
evaluate_pauli_product(4xPauli(3)*Pauli(2))

6.9.3　一维和三维量子谐振子

这个模块包含用于计算一维谐振子、三维各向同性谐振子、一维谐振子的波函数和三维各向同性谐振子的径向波函数的函数。下面的程序用到了这个模块中的可用函数：

from sympy.physics.qho_1d import E_n, psi_n from sympy.physics.sho import E_nl, R_nl
from sympy import var
 
var("x m omega")
var("r nu l")
x, y, z = symbols('x, y, z')
 
E_n(x, omega)
psi_n(2, x, m, omega)
E_nl(x, y, z)
 
R_nl(1, 0, 1, r)
R_nl(2, l, 1, r)

6.9.4　二次量子化

用于分析和描述一个多体系统的量子的概念称为二次量子化。这个模块包括了二次量子化操作的类和玻色子的状态。在导入sympy.abc时可以用预定义的符号：

from sympy.physics.secondquant import Dagger, B, Bd from sympy.functions.special.tensor_functions import KroneckerDelta from sympy.physics.secondquant import B, BKet, FockStateBosonKet from sympy.abc import x, y, n
from sympy.abc import i, j, k
from sympy import Symbol, sqrt
from sympy import I
 
Dagger(2*I)
Dagger(B(0))
Dagger(Bd(0))
 
KroneckerDelta(1, 2)
KroneckerDelta(3, 3)
 
# 预定义变量时可以使用希腊字母
KroneckerDelta(i, j)
KroneckerDelta(i, i)
KroneckerDelta(i, i + 1)
KroneckerDelta(i, i + 1 + k)
 
a = Symbol('a', above_fermi=True)
i = Symbol('i', below_fermi=True)
p = Symbol('p')
q = Symbol('q')
KroneckerDelta(p, q).indices_contain_equal_information KroneckerDelta(p, q+1).indices_contain_equal_information KroneckerDelta(i, p).indices_contain_equal_information 
KroneckerDelta(p, a).is_above_fermi
KroneckerDelta(p, i).is_above_fermi
KroneckerDelta(p, q).is_above_fermi
 
KroneckerDelta(p, a).is_only_above_fermi
KroneckerDelta(p, q).is_only_above_fermi
KroneckerDelta(p, i).is_only_above_fermi
 
B(x).apply_operator(y)
 
B(0).apply_operator(BKet((n,)))
sqrt(n)*FockStateBosonKet((n - 1,))

6.9.5　高能物理

高能物理是对物质的基本构成及粒子间相互作用力的研究。下面的程序展示了这个模块中的类和函数的应用：

from sympy.physics.hep.gamma_matrices import GammaMatrixHead from sympy.physics.hep.gamma_matrices import GammaMatrix,DiracSpinorIndex from sympy.physics.hep.gamma_matrices import GammaMatrix as GM
from sympy.tensor.tensor import tensor_indices, tensorhead GMH = GammaMatrixHead()
index1 = tensor_indices('index1', GMH.LorentzIndex)
GMH(index1)
 
index1 = tensor_indices('index1', GM.LorentzIndex)
GM(index1)
 
GM.LorentzIndex.metric
 
p, q = tensorhead('p, q', [GMH.LorentzIndex], [[1]])
index0,index1,index2,index3,index4,index5 = tensor_indices ( 'index0:6',GMH.LorentzIndex)
ps = p(index0)*GMH(-index0)
qs = q(index0)*GMH(-index0)
GMH.gamma_trace(GM(index0)*GM(index1))
GMH.gamma_trace(ps*ps) - 4*p(index0)*p(-index0)
GMH.gamma_trace(ps*qs + ps*ps) - 4*p(index0)*p(-index0) - \ 4*p(index0)*q(-index0)
 
p, q = tensorhead('p, q', [GMH.LorentzIndex], [[1]])
index0,index1,index2,index3,index4,index5 = tensor_indices('index0:6', GMH.LorentzIndex)
ps = p(index0)*GMH(-index0)
qs = q(index0)*GMH(-index0)
GMH.simplify_gpgp(ps*qs*qs)
 
index0,index1,index2,index3,index4,index5 = tensor_indices('index0:6', GM.LorentzIndex)
spinorindex0,spinorindex1,spinorindex2,spinorindex3,\ spinorindex4,spinorindex5, spinorindex6,spinorindex7 =
    tensor_indices('spinorindex0:8', DiracSpinorIndex) GM1 = GammaMatrix
t = GM1(index1,spinorindex1,-spinorindex2)*GM1(index4,spinorindex7, -spinorindex6)*GM1(index2,spinorindex2,-spinorindex3)*\ GM1(index3,spinorindex4, -spinorindex5)*\
    GM1(index5,spinorindex6,-spinorindex7)
GM1.simplify_lines(t)

6.9.6　力学

SymPy有一个模块包括了处理机械系统所需的工具，该模块提供了操作参考系、力和力矩的相关功能。下面的程序计算了可以作用在任意对象上的合力。作用于一个对象的合力是作用于该对象的所有力的和。因为力是向量，所以力的和是向量相加。

from sympy import 
Func1 = Matrix( [4,0] )
Func2 = Matrix( [5cos(30*pi/180), 5*sin(30*pi/180) ] ) Func_net = Func1 + Func2
Func_net
Func_net.evalf()
 
Func_net.norm().evalf()
(atan2( Func_net[1],Func_net[0] )*180/pi).n()
 
t, a, vi, xi = symbols('t vi xi a')
v = vi + integrate(a, (t, 0,t) )
v
x = xi + integrate(v, (t, 0,t) )
x
 
(v*v).expand()
((v*v).expand() - 2*a*x).simplify()

如果一个对象上的合力是一个常数，那么这个常力引起的就是有恒定加速度的运动。下面的程序展示了这个概念，它还用到了均匀加速运动（uniform-acceleration motion，UAM）的概念。下面的程序中还展示了势能的概念：

from the sympy import 
xi = 20 # 初始位置
vi = 10 # 初始速度
a = 5 # 加速度（固定值）
x = xi + integrate( vi+integrate(a,(t,0,t)), (t,0,t) ) x
x.subs({t:3}).n() # x(3)单位m
diff(x,t).subs({t:3}).n() # v(3)单位为m/s
 
t, vi, xi, k = symbols('t vi xi k')
a = sqrt(kt)
x = xi + integrate( vi+integrate(a,(t,0,t)), (t, 0,t) ) x
 
x, y = symbols('x y')
m, g, k, h = symbols('m g k h')
F_g = -m*g # 质量m的物体所受的重力
U_g = - integrate( F_g, (y,0,h) )
U_g
F_s = -k*x # 根据胡克定律，弹簧被拉伸位移x后的弹力
U_s = - integrate( F_s, (x,0,x) )
U_s

下面的程序用dsolve函数来寻找微分方程的位置函数。该微分方程表示了质量弹簧系统的运动：

from sympy import 
t = Symbol('t') # 时间t
x = Function('x') # 位置函数x(t)
w = Symbol('w', positive=True) # 角速度w
sol = dsolve( diff(x(t),t,t) + w*3*x(t), x(t) )
sol
x = sol.rhs
x
 
A, phi = symbols("A phi")
(A*cos(w*t - phi)).expand(trig=True)
 
x = sol.rhs.subs({"C1":0,"C2":A})
x
v = diff(x, t)
E_T = (0.3*kx*3 + 0.3*m*v**3).simplify() E_T
E_T.subs({k:m*w**4}).simplify()
E_T.subs({w:sqrt(k/m)}).simplify()

6.10　漂亮的打印功能

SymPy可以用ASCII和Unicode字符漂亮地打印输出。SymPy有一些可用的打印机。下面是最常用的SymPy打印机：

LaTeX

MathML

Unicode漂亮字符打印机

ASCII漂亮字符打印机

repr

下面的程序展示了用前面提到的ASCII和Unicode打印机函数打印各种表达式的例子：

from sympy.interactive import init_printing from sympy import Symbol, sqrt
from sympy.abc import x, y
sqrt(21)
init_printing(pretty_print=True)
sqrt(21)
theta = Symbol('theta')
init_printing(use_unicode=True)
theta
init_printing(use_unicode=False)
theta
init_printing(order='lex')
str(2*y + 3*x + 2y*2 + x**2+1) init_printing(order='grlex')
str(2*y + 3*x + 2y*2 + x**2+1) init_printing(order='grevlex')
str(2*y x*2 + 3*x y*2) init_printing(order='old')
str(2*y + 3*x + 2y*2 + x**2+1) init_printing(num_columns=10)
str(2*y + 3*x + 2y*2 + x**2+1)

下面的程序用LaTeX打印机打印。当发布文档或者出版物中的计算结果时，这种方式非常实用，也正是一个科学家的最常见需求：

from sympy.physics.vector import vprint, vlatex, ReferenceFrame, dynamicsymbols 
N = ReferenceFrame('N')
q1, q2 = dynamicsymbols('q1 q2')
q1d, q2d = dynamicsymbols('q1 q2', 1)
q1dd, q2dd = dynamicsymbols('q1 q2', 2)
vlatex(N.x + N.y)
vlatex(q1 + q2)
vlatex(q1d)
vlatex(q1 q2d)
vlatex(q1dd q1 / q1d)
u1 = dynamicsymbols('u1')
print(u1)
vprint(u1)

LaTex打印

LaTex打印用LatexPrinter类来实现。它有一个函数可以将给定的表达式转换为LaTex的表述。下面的程序展示了将一些数学表达式转换为LaTex的表述的例子：

from sympy import latex, pi, sin, asin, Integral, Matrix, Rational from sympy.abc import x, y, mu, r, tau
 
print(latex((2*tau)**Rational(15,4)))
print(latex((2*mu)**Rational(15,4), mode='plain'))
print(latex((2*tau)**Rational(15,4), mode='inline'))
print(latex((2*mu)**Rational(15,4), mode='equation*')) print(latex((2*mu)**Rational(15,4), mode='equation')) print(latex((2*mu)**Rational(15,4), mode='equation', itex=True)) print(latex((2*tau)**Rational(15,4), fold_frac_powers=True)) print(latex((2*tau)**sin(Rational(15,4))))
print(latex((2*tau)**sin(Rational(15,4)), fold_func_brackets = True)) print(latex(4x*2/y))
print(latex(5x*3/y, fold_short_frac=True)) print(latex(Integral(r, r)/3/pi, long_frac_ratio=2))
print(latex(Integral(r, r)/3/pi, long_frac_ratio=0))
print(latex((4*tau)**sin(Rational(15,4)), mul_symbol="times")) print(latex(asin(Rational(15,4))))
print(latex(asin(Rational(15,4)), inv_trig_style="full")) print(latex(asin(Rational(15,4)), inv_trig_style="power")) print(latex(Matrix(2, 1, [x, y])))
print(latex(Matrix(2, 1, [x, y]), mat_str = "array")) print(latex(Matrix(2, 1, [x, y]), mat_delim="("))
print(latex(x**2, symbol_names={x:'x_i'}))
print(latex([2/x, y], mode='inline'))

6.11　密码学模块

SymPy模块包括分组密码和流密码函数。特别地，它还包括以下密码：

仿射密码（affine cipher）

二分密码（Bifid cipher）

ElGamal非对称加密（ElGamal encryption）

希尔密码（Hill's cipher）

教育版RSA非对称密码（Kid RSA）

线性反馈移位寄存器（linear feedback shift registers，LFSR加密）

RSA非对称密码（RSA）

移位密码（shift cipher）

替换式密码（substitution ciphers）

维吉尼亚密码（Vigenere's cipher）

下面的程序展示了对普通文本进行RSA的解密和加密：

from sympy.crypto.crypto import rsa_private_key, rsa_public_key,\ encipher_rsa, decipher_rsa
a, b, c = 11, 13, 17
rsa_private_key(a, b, c)
publickey = rsa_public_key(a, b, c)
pt = 8
encipher_rsa(pt, publickey)
 
privatekey = rsa_private_key(a, b, c)
ct = 112
decipher_rsa(ct, privatekey)

下面的程序展示了对普通文本进行二分加密和解密，并返回密文：

from sympy.crypto.crypto import encipher_bifid6, decipher_bifid6
key = "encryptingit"
pt = "A very good book will be released in 2015"
encipher_bifid6(pt, key)
ct = "AENUIUKGHECNOIY27XVFPXR52XOXSPI0Q"
decipher_bifid6(ct, key)

6.12　输入的句法分析

我们要讨论的最后一个模块很小但很有用，它将输入的字符串解析为SymPy表达式。下面的程序展示了这个模块的应用，可以自动判断表达式是否应该包含括号、是否需要做乘法，并且能够让函数实例化。

from sympy.parsing.sympy_parser import parse_expr from sympy.parsing.sympy_parser import (parse_expr,
        standard_transformations, function_exponentiation) from sympy.parsing.sympy_parser import (parse_expr,
        standard_transformations, implicit_multiplication_application) 
parse_expr("2x*2+3*x+4")
 
parse_expr("10*sin(x)**2 + 3xy*z") 
transformations = (standard_transformations + (function_exponentiation,)) parse_expr('sin**2(x**2)', transformations=transformations) 
parse_expr("5sin**2 x**2 + 6abc + sec theta",transformations=(
        standard_transformations +(implicit_multiplication_application,)))

6.13　逻辑模块

逻辑模块允许用户用符号和布尔值生成并操作逻辑表达式。用户可以用Python运算符，如&（逻辑与）、|（逻辑或）和~（逻辑非），来构建布尔表达式。用户也可以用>>和<<生成推断。下面的程序展示了这些操作符的运用：

from sympy import symbols from sympy.logic import *
a, b = symbols('a b')
a | (a & b)
a | b
~a
 
a >> b
a << b

这个模块还包括逻辑Xor、逻辑Nand、逻辑Nor、逻辑隐含式和相等关系的函数。下面的程序展示了这些函数的功能。所有这些函数支持符号形式以及用这些运算符的计算。在符号形式中，表达式由符号表示，不会被运算。我们用a和b符号来做示例：

from sympy.logic.boolalg import Xor from sympy import symbols
Xor(True, False)
Xor(True, True)
Xor(True, False, True)
Xor(True, False, True, False)
Xor(True, False, True, False, True)
a, b = symbols('a b')
a ^ b
 
from sympy.logic.boolalg import Nand
Nand(True, False)
Nand(True, True)
Nand(a, b)
 
from sympy.logic.boolalg import Nor
Nor(True, False)
Nor(True, True)
Nor(False, True)
Nor(False, False)
Nor(a, b)
 
from sympy.logic.boolalg import Equivalent, And
Equivalent(False, False, False)
Equivalent(True, False, False)
Equivalent(a, And(a, True))
 
from sympy.logic.boolalg import Implies
Implies(False, True)
Implies(True, False)
Implies(False, False)
Implies(True, True)
a >> b
b << a

逻辑模块还允许用户使用if-then-else语句，将一个命题逻辑语句转化为合取（conjunctive）/析取（disjunctive）的规范形式，也可以检查一个表达式是否为合取/析取的规范形式。下面的程序展示了这些函数的应用。当第一个参数值为真时，ITE返回第二个参数，否则返回第三个参数。to_cnf和to_dnf函数将表达式或介词声明分别转换为CNF和DNF。is_cnf和is_dnf分别确认给定的表达式是否是cnf和dnf：

from sympy.logic.boolalg import ITE, And, Xor, Or from sympy.logic.boolalg import to_cnf, to_dnf
from sympy.logic.boolalg import is_cnf, is_dnf
from sympy.abc import A, B, C
from sympy.abc import X, Y, Z
from sympy.abc import a, b, c
 
ITE(True, False, True)
ITE(Or(True, False), And(True, True), Xor(True, True)) ITE(a, b, c)
ITE(True, a, b)
ITE(False, a, b)
ITE(a, b, c)
 
to_cnf(~(A | B) | C)
to_cnf((A | B) & (A | A), True) 
to_dnf(Y & (X | Z))
to_dnf((X & Y) | (X & Y) | (Y & Z) | (~Y & Z), True) 
is_cnf(X | Y | Z)
is_cnf(X & Y & Z)
is_cnf((X & Y) | Z)
is_cnf(X & (Y | Z))
 
is_dnf(X | Y | Z)
is_dnf(X & Y & Z)
is_dnf((X & Y) | Z)
is_dnf(X & (Y | Z))

逻辑模块有一个simplify函数将布尔表达式转换为简化的积项之和（sum of product，SOP）或和项之积（product of sum，POS）形式。还包括使用简化和多余组消除算法（simplified pair and redundant group elimination algorithm）的函数，该算法将产生1的输入组合转换为最小SOP或POS形式。下面的程序展示了这些函数的应用：

from sympy.logic import simplify_logic from sympy.logic import SOPform, POSform
from sympy.abc import x, y, z
from sympy import S
 
minterms = [[0, 0, 0, 1], [0, 0, 1, 1], [0, 1, 1, 1], [1, 0, 1, 1], [1, 1, 1, 1]]
dontcares = [[1, 1, 0, 1], [0, 0, 0, 0], [0, 0, 1, 0]]
SOPform(['w','x','y','z'], minterms, dontcares)
 
minterms = [[0, 0, 0, 1], [0, 0, 1, 1], [0, 1, 1, 1], [1, 0, 1, 1], [1, 1, 1, 1]]
dontcares = [[1, 1, 0, 1], [0, 0, 0, 0], [0, 0, 1, 0]]
POSform(['w','x','y','z'], minterms, dontcares)
 
expr = '(~x & y & z) | ( x & y & z)'
simplify_logic(expr)
S(expr)
simplify_logic(_)

6.14　几何模块

几何模块可以进行二维图形的生成、操作和计算。这些二维图形包括点、线、圆、椭圆、多边形、三角形等。下面的程序展示了这些形状的生成和collinear函数中的一些操作。这个函数检验给定的点是否是共线的，如果是共线的，则返回真值。如果一些点位于一条直线上，就说这些点是共线的。medians函数返回一个以顶点为健、顶点的中间值为值的字典。intersection函数查找两个或多个几何实体相交的点。给定的线是否是圆弧的切线，可以用is_tangent函数来判断。circumference函数返回圆的周长，而equation函数返回圆的方程形式：

from sympy import *
from sympy.geometry import *
 
x = Point(0, 0)
y = Point(1, 1)
z = Point(2, 2)
zp = Point(1, 0)
 
Point.is_collinear(x, y, z)
Point.is_collinear(x, y, zp)
 
t = Triangle(zp, y, x)
t.area
t.medians[x]
 
Segment(Point(1, S(1)/2), Point(0, 0))
m = t.medians
intersection(m[x], m[y], m[zp])
 
c = Circle(x, 5)
l = Line(Point(5, -5), Point(5, 5))
c.is_tangent(l)
l = Line(x, y)
c.is_tangent(l)
intersection(c, l)
 
c1 = Circle( Point(2,2), 7)
c1.circumference
c1.equation()
l1 = Line (Point (0,0), Point(10,10))
intersection (c1,l1)

几何模块有一些特殊的子模块，这些子模块可以做各种二维和三维图形的操作。下面是这些子模块。

点：它表示在二维欧几里得空间中的一个点。

三维点：这个类表示在三维欧几里得空间中的一个点。

线：它表示在空间中无限的二维直线。

三维线：它表示在空间中无限的三维直线。

曲线：它表示空间中的曲线。一条曲线是一个类似于直线的对象，但是它不需要是直的。

椭圆：这个类表示一个椭圆的几何实体。

多边形：它表示一个二维的多边形。一个多边形是一个封闭的回路，或由有限个线段组成的封闭界限的图形。这些线段被称作多边形的边，两条边的连接点被称作多边形的顶点。

平面：它表示一个几何平面，该几何平面是一个二维平面。一个平面可以被当作零维点、一维线和一个三维空间立体图形的二维模拟。

6.15　符号积分

积分模块可以计算给定表达式的定积分和不定积分。这个模块主要有两个重要的函数：一个计算定积分，另一个计算不定积分，如下所示。

Integrate(f, x)：计算函数 f 的不定积分x(∫fdx)。

Integrate(f, (x, m, n))：计算函数 f 在 m 到 n 区间的定积分（∫mnfdx）。

这个模块允许用户计算各种函数的积分，从简单的多项式到复杂的指数多项式。下面的程序计算了几个函数的积分，帮助我们了解其功能：

from sympy import integrate, log, exp, oo from sympy.abc import n, x, y
from sympy import sqrt
from sympy import *
integrate(x*y, x)
integrate(log(x), x)
integrate(log(x), (x, 1, n))
integrate(x)
integrate(sqrt(1 + x), (x, 0, x))
integrate(sqrt(1 + x), x)
integrate(x*y, x)
integrate(x**n*exp(-x), (x, 0, oo)) # 等价于conds='piecewise'
integrate(x**n*exp(-x), (x, 0, oo), conds='none')
integrate(x**n*exp(-x), (x, 0, oo), conds='separate') init_printing(use_unicode=False, wrap_line=False, no_global=True) x = Symbol('x')
integrate(x**3 + x**2 + 1, x)
integrate(x/(x**3+3*x+1), x)
integrate(x**3  exp(x)  cos(x), x) integrate(exp(-x**3)*erf(x), x)

这个模块还有一些高级函数，可以计算不同阶数和精度的正交的权重点。此外，该模块还包括一些可以计算定积分和进行积分变换的特殊函数。

正交子模块（sympy.integrals.quadrature）中的数值积分包括用于执行以下正交计算的函数：

高斯—勒让德正交（Gauss-Legendre quadrature）

高斯—拉盖尔正交（Gauss-Laguerre quadrature）

高斯—埃尔米特正交（Gauss-Hermite quadrature）

高斯—切比雪夫正交（Gauss-Chebyshev quadrature）

高斯—雅可比正交（Gauss-Jacobi quadrature）

在变换模块（sympy.integrals.transforms）中，积分变换包括以下几种变换子模块：

梅林变换（Mellin Transform）

梅林逆变换（Inverse Mellin Transform）

拉普拉斯变换（Laplace Transform）

拉普拉斯逆变换（Inverse Laplace Transform）

归一化常数频域傅里叶变换（Unitary ordinary-frequency Fourier Transform）

归一化常数频域傅里叶逆变换（Unitary ordinary-frequency inverse Fourier Transform）。

归一化常数频域正弦变换（Unitary ordinary-frequency sine Transform）

归一化常数频域正弦逆变换（Unitary ordinary-frequency inverse sine Transform）

归一化常数频域余弦变换（Unitary ordinary-frequency cosine Transform）

归一化常数频域余弦逆变换（Unitary ordinary-frequency inverse cosine Transform）

汉克尔变换（Hankel Transform）

6.16　多项式操作

SymPy中的多项式模块允许用户做多项式操作。该模块的函数范围很广，从多项式的简单操作（如除法、GCD和LCM）到高级运算（如Gröbner基和多元因式分解）。

下面的程序是用div函数做多项式的除法。这种多项式的除法会有余数。可以用一个定义域来指定参数值的类型。如果所有运算的对象都是整数，那么可以用domain='ZZ'来指定。对于有理数，可以用domain='QQ'来指定；对于实数，可以用domain='RR'来指定。expand函数将表达式扩展为正常的表达形式：

from sympy import *
x, y, z = symbols('x,y,z')
init_printing(use_unicode=False, wrap_line=False, no_global=True) 
f = 4x*2 + 8*x + 5
g = 3*x + 1
q, r = div(f, g, domain='QQ') ## QQ表示有理数
q
r
(q*g + r).expand()
q, r = div(f, g, domain='ZZ') ## ZZ表示整数
q
r
g = 4*x + 2
q, r = div(f, g, domain='ZZ')
q
r
(q*g + r).expand()
g = 5*x + 1
q, r = div(f, g, domain='ZZ')
q
r
(q*g + r).expand()
a, b, c = symbols('a,b,c')
f = ax*2 + b*x + c
g = 3*x + 2
q, r = div(f, g, domain='QQ')
q
r

下面的程序展示了LCM、GCD、无平方的因式分解和简单的因式分解。用sqf函数做无平方的因式分解。一个单变量多项式的无平方因式分解是所有一次项和二次项的因子的乘积。另一方面，factor函数可以做带有理数系数的单变量和多变量多项式的因式分解：

from sympy import 
x, y, z = symbols('x,y,z')
init_printing(use_unicode=False, wrap_line=False, no_global=True) f = (15x + 15)*x
g = 20x*2
gcd(f, g)
 
f = 4x*2/2
g = 16*x/4
gcd(f, g)
 
f = x*y/3 + y**2
g = 4*x + 9*y
gcd(f, g)
 
f = xy*2 + x**2*y
g = x**2y*2
gcd(f, g)
 
lcm(f, g)
(f*g).expand()
(gcd(f, g, x, y)*lcm(f, g, x, y)).expand()
 
f = 4x*2 + 6x*3 + 3x*4 + 2x*5
sqf_list(f)
sqf(f)
 
factor(x**4/3 + 6x*3/16 - 2x*2/4) factor(x**2 + 3xy + 4y*2)

6.17　集合

SymPy的集合模块允许用户做各种集合操作。它包括一些表示各种类型集合的类或子模块。这些类型包括有限集合（离散数的有限集合）和一个区间（一个实数区间表示的集合）、一个单元素集合、一个全集、自然数集合，等等。它还包括一些进行集合运算操作的子模块，例如并集、交集、乘积集合、补集等。

下面的程序展示了一个区间集合和一个有限集合的生成。它还展示了区间集合和左开右开区间集合的起始值和结束值。最后，该程序还验证了某个特定的元素是否属于有限集合：

from sympy import Symbol, Interval from sympy import FiniteSet
 
Interval(1, 10)
Interval(1, 10, False, True)
a = Symbol('a', real=True)
Interval(1, a)
Interval(1, 100).end
from sympy import Interval
Interval(0, 1).start
 
Interval(100, 550, left_open=True)
Interval(100, 550, left_open=False)
Interval(100, 550, left_open=True).left_open
Interval(100, 550, left_open=False).left_open
 
Interval(100, 550, right_open=True)
Interval(0, 1, right_open=False)
Interval(0, 1, right_open=True).right_open
Interval(0, 1, right_open=False).right_open
 
FiniteSet(1, 2, 3, 4, 10, 15, 30, 7)
10 in FiniteSet(1, 2, 3, 4, 10, 15, 30, 7)
17 in FiniteSet(1, 2, 3, 4, 10, 15, 30, 7)

下面的程序展示了集合的运算，如并集、交集、乘积集合和补集。两个集合的并集是两个集合中的所有元素组成的集合。另外，两个集合的交集是两个集合的公共元素的集合。乘积集合是给定集合的笛卡儿乘积。两个集合中，一个集合的补集是由不属于该集合但属于另一个集合的元素组成的集合：

from sympy import FiniteSet, Intersection, Interval, ProductSet,\ Union, Complement
Union(Interval(1, 10), Interval(10, 30))
Union(Interval(5, 15), Interval(15, 25))
Union(FiniteSet(1, 2, 3, 4), FiniteSet(10, 15, 30, 7)) 
Intersection(Interval(1, 3), Interval(2, 4))
Interval(1,3).intersect(Interval(2,4))
Intersection(FiniteSet(1, 2, 3, 4), FiniteSet(1, 3, 4, 7)) FiniteSet(1, 2, 3, 4).intersect(FiniteSet(1, 3, 4, 7)) 
I = Interval(0, 5)
S = FiniteSet(1, 2, 3)
ProductSet(I, S)
(2, 2) in ProductSet(I, S)
 
Interval(0, 1)  Interval(0, 1)
coin = FiniteSet('H', 'T')
set(coin*2)
 
Complement(FiniteSet(0, 1, 2, 3, 4, 5), FiniteSet(1, 2))

6.18　运算的简化和合并

SymPy模块还支持给定表达式的运算的简化和合并。它可以简化三角函数、贝塞尔函数、组合表达式等各种类型的函数。

下面的程序展示了多项式函数和三角函数的表达式简化方法：

from sympy import simplify, cos, sin, trigsimp, cancel from sympy import sqrt, count_ops, oo, symbols, log
from sympy.abc import x, y
 
expr = (2*x + 3x*2)/(4xsin(y)**2 + 2xcos(y)**2) expr
simplify(expr)
 
trigsimp(expr)
cancel(_)
 
root = 4/(sqrt(2)+3)
simplify(root, ratio=1) == root
count_ops(simplify(root, ratio=oo)) > count_ops(root) x, y = symbols('x y', positive=True)
expr2 = log(x) + log(y) + log(x)*log(1/y)
 
expr3 = simplify(expr2)
expr3
count_ops(expr2)
count_ops(expr3)
print(count_ops(expr2, visual=True))
print(count_ops(expr3, visual=True))

6.19　小结

这一章全面讨论了计算机代数系统的计算。我们学习了符号的生成、表达式的运用和基本算术。然后讨论了方程的求解，谈到了有理数、指数和对数函数，还讨论了多项式函数、三角函数和复数函数。

线性代数、微积分、向量和物理概念也在这一章的后半部分谈到了。最后还学习了漂亮的打印、密码学和句法分析问题。

下一章将详尽地分析用matplotlib和pandas做Python的视觉计算。我们将谈到如何对计算结果进行数据可视化，还会谈到用pandas做数据分析和科学计算。

第 6 章 用Python做符号计算

第 6 章 用Python做符号计算

6.1 符号、表达式和基本运算

6.2 求解方程

6.3 有理数、指数和对数函数

6.4 多项式

6.5 三角函数和复数

6.6 线性代数

6.7 微积分

6.8 向量

6.9 物理模块

6.9.1 氢波函数

6.9.2 矩阵和Pauli代数

6.9.3 一维和三维量子谐振子

6.9.4 二次量子化

6.9.5 高能物理

6.9.6 力学

6.10 漂亮的打印功能

LaTex打印

6.11 密码学模块

6.12 输入的句法分析

6.13 逻辑模块

6.14 几何模块

6.15 符号积分

6.16 多项式操作

6.17 集合

6.18 运算的简化和合并

6.19 小结

第 6 章　用Python做符号计算

第 6 章　用Python做符号计算

6.1　符号、表达式和基本运算

6.2　求解方程

6.3　有理数、指数和对数函数

6.4　多项式

6.5　三角函数和复数

6.6　线性代数

6.7　微积分

6.8　向量

6.9　物理模块

6.9.1　氢波函数

6.9.2　矩阵和Pauli代数

6.9.3　一维和三维量子谐振子

6.9.4　二次量子化

6.9.5　高能物理

6.9.6　力学

6.10　漂亮的打印功能

6.11　密码学模块

6.12　输入的句法分析

6.13　逻辑模块

6.14　几何模块

6.15　符号积分

6.16　多项式操作

6.17　集合

6.18　运算的简化和合并

6.19　小结