一、简介
sympy官网地址:https://docs.sympy.org/latest/index.html
当前最新正式版本:1.10.1
sympy是一个Python的科学计算库,用一套强大的符号计算体系完成诸如多项式求值、求极限、解方程、求积分、微分方程、级数展开、矩阵运算等等计算问题。虽然Matlab的类似科学计算能力也很强大,但是Python以其语法简单、易上手、异常丰富的三方库生态,个人认为可以更优雅地解决日常遇到的各种计算问题。
如果你遇到了一个难题,不要犹豫,来找Python,它几乎不会让你失望的。写本文的初衷也是妹子在做金融作业的时候遇到大量的计算,用普通计算器算真是工程量浩大,所以找到sympy这么个库,写代码帮忙辅助计算一下,然后在这里写博客记录下来。
安装sympy库
pip3 install sympy
二、常用的sympy内置符号
虚数单位i
In [13]: import sympy
In [14]: sympy.I
Out[14]: I
In [15]: sympy.I ** 2
Out[15]: -1
# 求-1的平方根
In [16]: sympy.sqrt(-1)
Out[16]: I
注:本文后面的示例都省略导包语句:import sympy
自然对数的底e
In [18]: sympy.E
Out[18]: E
# 求对数
In [20]: sympy.log(sympy.E)
Out[20]: 1
无穷大oo
In [26]: 1/sympy.oo
Out[26]: 0
In [27]: 1 + sympy.oo
Out[27]: oo
圆周率pi
In [60]: sympy.pi
Out[60]: pi
In [61]: sympy.sin(sympy.pi/2)
Out[61]: 1
三、用sympy进行初等运算
Python 2.x中用除号/
做两个整数的除法,实际上是整除运算,为了防止这种情况的发生,避免不必要的麻烦,下文的所有示例一开始都加上一句:from __future__ import division
,这个时候除号/
本身就变成了真实除法,而//
才是整除,比如:
# 导入division包之前
In [1]: 1/2
Out[1]: 0
In [2]: from __future__ import division
# 导入division包之后
In [3]: 1/2
Out[3]: 0.5
In [4]: 1//2
Out[4]: 0
求对数
# 自然对数
In [10]: sympy.log(sympy.E)
Out[10]: 1
In [11]: sympy.log(sympy.E ** 3)
Out[11]: 3
# 以10为底1000的对数
In [12]: sympy.log(1000,10)
Out[12]: 3
求平方根
In [13]: sympy.sqrt(4)
Out[13]: 2
In [14]: sympy.sqrt(-1)
Out[14]: I
求n次方根
# 求8的3次方根
In [15]: sympy.root(8,3)
Out[15]: 2
求k次方
In [21]: 2 ** 3
Out[21]: 8
In [22]: 16 ** (1/2)
Out[22]: 4.0
求阶乘
In [35]: sympy.factorial(4)
Out[35]: 24
求三角函数
以sin
函数为例:
In [86]: sympy.sin(sympy.pi)
Out[86]: 0
In [87]: sympy.sin(sympy.pi/2)
Out[87]: 1
四、表达式与表达式求值
sympy可以用一套符号系统来表示一个表达式,如函数、多项式等,并且可以进行求值,比如:
# 首先定义x为一个符号,表示一个变量
In [96]: x = sympy.Symbol('x')
In [97]: fx = 2*x + 1
# 可以看到fx是一个sympy.core.add.Add类型的对象,也就是一个表达式
In [98]: type(fx)
Out[98]: sympy.core.add.Add
# 用evalf函数,传入变量的值,对表达式进行求值
In [101]: fx.evalf(subs={x:2})
Out[101]: 5.00000000000000
还支持多元表达式:
In [102]: x,y = sympy.symbols('x y')
In [103]: f = 2 * x + y
# 以字典的形式传入多个变量的值
In [104]: f.evalf(subs = {x:1,y:2})
Out[104]: 4.00000000000000
# 如果只传入一个变量的值,则原本输出原来的表达式
In [105]: f.evalf(subs = {x:1})
Out[105]: 2.0*x + y
五、用sympy解方程(组)
使用sympy.solve
函数解方程,该函数通常传入两个参数,第1个参数是方程的表达式(把方程所有的项移到等号的同一边形成的式子),第2个参数是方程中的未知数。函数的返回值是一个列表,代表方程的所有根(可能为复数根)。
解最简单的方程
比如下面我们来求两个方程:
# 首先定义 `x`为一个符号,代表一个未知数
In [24]: x = sympy.Symbol('x')
# 解方程:x - 1 = 0
In [25]: sympy.solve(x - 1,x)
Out[25]: [1]
# 解方程:x ^ 2 - 1 = 0
In [26]: sympy.solve(x ** 2 - 1,x)
Out[26]: [-1, 1]
# 解方程:x ^ 2 + 1 = 0
In [27]: sympy.solve(x ** 2 + 1,x)
Out[27]: [-I, I]
把函数式赋给一个变量
有时候为了书写起来简洁,可以把一个函数式起个名字,比如:
In [30]: x = sympy.Symbol('x')
In [31]: f = x + 1
In [32]: sympy.solve(f,x)
Out[32]: [-1]
解方程组
比如要解这么个二元一次方程组:
代码如下:
# 一次性定义多个符号
In [28]: x,y = sympy.symbols('x y')
In [29]: sympy.solve([x + y - 1,x - y -3],[x,y])
Out[29]: {x: 2, y: -1}
计算求和式
计算求和式可以使用sympy.summation
函数,其函数原型为:sympy.summation(f, *symbols, **kwargs)
。
话不多少,举个栗子,比如求下面这个求和式子的值:
我们用初中的知识可以知道,这个式子的结果为:5050 * 2 = 10100
下面用代码来求:
In [37]: n = sympy.Symbol('n')
In [38]: sympy.summation(2 * n,(n,1,100))
Out[38]: 10100
可见结果是正确的。
如果sympy.summation
函数无法计算出具体的结果,那么会返回求和表达式。
解带有求和式的方程
比如求这么一个方程:
代码如下:
In [43]: x = sympy.Symbol('x')
In [44]: i = sympy.Symbol('i',integer = True)
In [46]: f = sympy.summation(x,(i,1,5)) + 10 * x - 15
In [47]: sympy.solve(f,x)
Out[47]: [1]
六、求极限
求极限用sympy.limit
函数,其函数文档如下:
Signature: sympy.limit(e, z, z0, dir='+')
Docstring:
Compute the limit of e(z) at the point z0.
z0 can be any expression, including oo and -oo.
For dir="+" (default) it calculates the limit from the right
(z->z0+) and for dir="-" the limit from the left (z->z0-). For infinite
z0 (oo or -oo), the dir argument is determined from the direction
of the infinity (i.e., dir="-" for oo).
函数文档中已经说得很清楚了,下面用代码示例来求几个极限。
如果学过微积分,就会知道微积分中有3个重要的极限:
下面就用sympy.limit
函数来分别求这3个极限:
In [53]: x = sympy.Symbol('x')
In [54]: f1 = sympy.sin(x)/x
In [55]: sympy.limit(f1,x,0)
Out[55]: 1
In [56]: f2 = (1+x)**(1/x)
In [57]: sympy.limit(f2,x,0)
Out[57]: E
In [58]: f3 = (1+1/x)**x
In [59]: sympy.limit(f3,x,sympy.oo)
Out[59]: E
可见三个极限的计算结果都完全正确。
三、求导
求导使用sympy.diff
函数,传入2个参数:函数表达式和变量名,举例如下:
In [63]: x = sympy.Symbol('x')
In [64]: f = x ** 2 + 2 * x + 1
In [65]: sympy.diff(f,x)
Out[65]: 2*x + 2
In [66]: f2 = sympy.sin(x)
In [67]: sympy.diff(f2,x)
Out[67]: cos(x)
# 多元函数求偏导
In [68]: y = sympy.Symbol('y')
In [70]: f3 = x**2 + 2*x + y**3
In [71]: sympy.diff(f3,x)
Out[71]: 2*x + 2
In [72]: sympy.diff(f3,y)
Out[72]: 3*y**2
八、求定积分
使用sympy.integrate
函数求定积分,其功能比较复杂,非常强大,下面仅仅举几个比较简单的例子。
先来求一个最简单的积分:
用牛顿-莱布尼兹公式
可以立马口算出上面这个式子的结果是1,用代码计算如下:
n [74]: x = sympy.Symbol('x')
n [75]: f = 2 * x
# 传入函数表达式和积分变量、积分下限、上限
n [76]: sympy.integrate(f,(x,0,1))
ut[76]: 1
下面来算一个复杂一点的多重积分:
其中:
我们通过口算可以求出f(x)
:
所以:
下面用代码来计算上述过程:
In [82]: t,x = sympy.symbols('t x')
In [83]: f = 2 * t
In [84]: g = sympy.integrate(f,(t,0,x))
In [85]: sympy.integrate(g,(x,0,3))
Out[85]: 9
九、求不定积分
同样也是使用sympy.integrate
函数求不定积分,下面仅仅举几个比较简单的例子。
比如求下面这个不定积分:
通过观察我们知道它的结果是:
下面用代码来计算这个不定积分的结果:
In [79]: x = sympy.Symbol('x')
In [80]: f = sympy.E ** x + 2 * x
In [81]: sympy.integrate(f,x)
Out[81]: x**2 + exp(x)
十、矩阵运算
>>> from sympy import *
>>> init_printing(use_unicode=True)
To make a matrix in SymPy, use the Matrix
object. A matrix is constructed by providing a list of row vectors that make up the matrix. For example, to construct the matrix
>>> Matrix([[1, -1], [3, 4], [0, 2]])⎡1 -1⎤
⎢ ⎥
⎢3 4 ⎥
⎢ ⎥
⎣0 2 ⎦
To make it easy to make column vectors, a list of elements is considered to be a column vector.
>>> Matrix([1, 2, 3])⎡1⎤
⎢ ⎥
⎢2⎥
⎢ ⎥
⎣3⎦
Matrices are manipulated just like any other object in SymPy or Python.
>>> M = Matrix([[1, 2, 3], [3, 2, 1]])>>> N = Matrix([0, 1, 1])>>> M*N⎡5⎤
⎢ ⎥
⎣3⎦
One important thing to note about SymPy matrices is that, unlike every other object in SymPy, they are mutable. This means that they can be modified in place, as we will see below. The downside to this is that Matrix
cannot be used in places that require immutability, such as inside other SymPy expressions or as keys to dictionaries. If you need an immutable version of Matrix
, use ImmutableMatrix
.
(一)基本运算
Here are some basic operations on Matrix
.
Shape
To get the shape of a matrix, use shape() function.
>>> from sympy import shape>>> M = Matrix([[1, 2, 3], [-2, 0, 4]])>>> M⎡1 2 3⎤
⎢ ⎥
⎣-2 0 4⎦
>>> shape(M)(2, 3)
Accessing Rows and Columns
To get an individual row or column of a matrix, use row
or col
. For example, M.row(0)
will get the first row. M.col(-1)
will get the last column.
>>> M.row(0)[1 2 3]
>>> M.col(-1)⎡3⎤
⎢ ⎥
⎣4⎦
Deleting and Inserting Rows and Columns
To delete a row or column, use row_del
or col_del
. These operations will modify the Matrix in place.
>>> M.col_del(0)>>> M⎡2 3⎤
⎢ ⎥
⎣0 4⎦
>>> M.row_del(1)>>> M[2 3]
To insert rows or columns, use row_insert
or col_insert
. These operations do not operate in place.
>>> M[2 3]
>>> M = M.row_insert(1, Matrix([[0, 4]]))>>> M⎡2 3⎤
⎢ ⎥
⎣0 4⎦
>>> M = M.col_insert(0, Matrix([1, -2]))>>> M⎡1 2 3⎤
⎢ ⎥
⎣-2 0 4⎦
Unless explicitly stated, the methods mentioned below do not operate in place. In general, a method that does not operate in place will return a new Matrix
and a method that does operate in place will return None
.
(二)基本方法
As noted above, simple operations like addition and multiplication are done just by using +
, *
, and **
. To find the inverse of a matrix, just raise it to the -1
power.
>>> M = Matrix([[1, 3], [-2, 3]])>>> N = Matrix([[0, 3], [0, 7]])>>> M + N⎡1 6 ⎤
⎢ ⎥
⎣-2 10⎦
>>> M*N⎡0 24⎤
⎢ ⎥
⎣0 15⎦
>>> 3*M⎡3 9⎤
⎢ ⎥
⎣-6 9⎦
>>> M**2
⎡-5 12⎤
⎢ ⎥
⎣-8 3 ⎦
>>> M**-1
⎡1/3 -1/3⎤
⎢ ⎥
⎣2/9 1/9 ⎦
>>> N**-1
Traceback (most recent call last):
...
NonInvertibleMatrixError: Matrix det == 0; not invertible.
To take the transpose of a Matrix, use T
.
>>> M = Matrix([[1, 2, 3], [4, 5, 6]])>>> M⎡1 2 3⎤
⎢ ⎥
⎣4 5 6⎦
>>> M.T⎡1 4⎤
⎢ ⎥
⎢2 5⎥
⎢ ⎥
⎣3 6⎦
(三)矩阵构造器
Several constructors exist for creating common matrices. To create an identity matrix, use eye
. eye(n)
will create an n×n identity matrix.
>>> eye(3)⎡1 0 0⎤
⎢ ⎥
⎢0 1 0⎥
⎢ ⎥
⎣0 0 1⎦
>>> eye(4)⎡1 0 0 0⎤
⎢ ⎥
⎢0 1 0 0⎥
⎢ ⎥
⎢0 0 1 0⎥
⎢ ⎥
⎣0 0 0 1⎦
To create a matrix of all zeros, use zeros
. zeros(n, m)
creates an n×m matrix of 0s.
>>> zeros(2, 3)⎡0 0 0⎤
⎢ ⎥
⎣0 0 0⎦
Similarly, ones
creates a matrix of ones.
>>> ones(3, 2)⎡1 1⎤
⎢ ⎥
⎢1 1⎥
⎢ ⎥
⎣1 1⎦
To create diagonal matrices, use diag
. The arguments to diag
can be either numbers or matrices. A number is interpreted as a 1×1 matrix. The matrices are stacked diagonally. The remaining elements are filled with 0s.
>>> diag(1, 2, 3)⎡1 0 0⎤
⎢ ⎥
⎢0 2 0⎥
⎢ ⎥
⎣0 0 3⎦
>>> diag(-1, ones(2, 2), Matrix([5, 7, 5]))⎡-1 0 0 0⎤
⎢ ⎥
⎢0 1 1 0⎥
⎢ ⎥
⎢0 1 1 0⎥
⎢ ⎥
⎢0 0 0 5⎥
⎢ ⎥
⎢0 0 0 7⎥
⎢ ⎥
⎣0 0 0 5⎦
(四)高级方法
行列式
To compute the determinant of a matrix, use det
.
>>> M = Matrix([[1, 0, 1], [2, -1, 3], [4, 3, 2]])>>> M⎡1 0 1⎤
⎢ ⎥
⎢2 -1 3⎥
⎢ ⎥
⎣4 3 2⎦
>>> M.det()-1
RREF
To put a matrix into reduced row echelon form, use rref
. rref
returns a tuple of two elements. The first is the reduced row echelon form, and the second is a tuple of indices of the pivot columns.
>>> M = Matrix([[1, 0, 1, 3], [2, 3, 4, 7], [-1, -3, -3, -4]])>>> M⎡1 0 1 3 ⎤
⎢ ⎥
⎢2 3 4 7 ⎥
⎢ ⎥
⎣-1 -3 -3 -4⎦
>>> M.rref()⎛⎡1 0 1 3 ⎤ ⎞
⎜⎢ ⎥ ⎟
⎜⎢0 1 2/3 1/3⎥, (0, 1)⎟
⎜⎢ ⎥ ⎟
⎝⎣0 0 0 0 ⎦ ⎠
Note
The first element of the tuple returned by rref
is of type Matrix
. The second is of type tuple
.
Nullspace
To find the nullspace of a matrix, use nullspace
. nullspace
returns a list
of column vectors that span the nullspace of the matrix.
>>> M = Matrix([[1, 2, 3, 0, 0], [4, 10, 0, 0, 1]])>>> M⎡1 2 3 0 0⎤
⎢ ⎥
⎣4 10 0 0 1⎦
>>> M.nullspace()⎡⎡-15⎤ ⎡0⎤ ⎡ 1 ⎤⎤
⎢⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎥
⎢⎢ 6 ⎥ ⎢0⎥ ⎢-1/2⎥⎥
⎢⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎥
⎢⎢ 1 ⎥, ⎢0⎥, ⎢ 0 ⎥⎥
⎢⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎥
⎢⎢ 0 ⎥ ⎢1⎥ ⎢ 0 ⎥⎥
⎢⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎥
⎣⎣ 0 ⎦ ⎣0⎦ ⎣ 1 ⎦⎦
Columnspace
To find the columnspace of a matrix, use columnspace
. columnspace
returns a list
of column vectors that span the columnspace of the matrix.
>>> M = Matrix([[1, 1, 2], [2 ,1 , 3], [3 , 1, 4]])>>> M⎡1 1 2⎤
⎢ ⎥
⎢2 1 3⎥
⎢ ⎥
⎣3 1 4⎦
>>> M.columnspace()⎡⎡1⎤ ⎡1⎤⎤
⎢⎢ ⎥ ⎢ ⎥⎥
⎢⎢2⎥, ⎢1⎥⎥
⎢⎢ ⎥ ⎢ ⎥⎥
⎣⎣3⎦ ⎣1⎦⎦
Eigenvalues, Eigenvectors, and Diagonalization
To find the eigenvalues of a matrix, use eigenvals
. eigenvals
returns a dictionary of eigenvalue: algebraic_multiplicity
pairs (similar to the output of roots).
>>> M = Matrix([[3, -2, 4, -2], [5, 3, -3, -2], [5, -2, 2, -2], [5, -2, -3, 3]])>>> M⎡3 -2 4 -2⎤
⎢ ⎥
⎢5 3 -3 -2⎥
⎢ ⎥
⎢5 -2 2 -2⎥
⎢ ⎥
⎣5 -2 -3 3 ⎦
>>> M.eigenvals(){-2: 1, 3: 1, 5: 2}
This means that M
has eigenvalues -2, 3, and 5, and that the eigenvalues -2 and 3 have algebraic multiplicity 1 and that the eigenvalue 5 has algebraic multiplicity 2.
To find the eigenvectors of a matrix, use eigenvects
. eigenvects
returns a list of tuples of the form (eigenvalue, algebraic_multiplicity, [eigenvectors])
.
>>> M.eigenvects()⎡⎛ ⎡⎡0⎤⎤⎞ ⎛ ⎡⎡1⎤⎤⎞ ⎛ ⎡⎡1⎤ ⎡0 ⎤⎤⎞⎤
⎢⎜ ⎢⎢ ⎥⎥⎟ ⎜ ⎢⎢ ⎥⎥⎟ ⎜ ⎢⎢ ⎥ ⎢ ⎥⎥⎟⎥
⎢⎜ ⎢⎢1⎥⎥⎟ ⎜ ⎢⎢1⎥⎥⎟ ⎜ ⎢⎢1⎥ ⎢-1⎥⎥⎟⎥
⎢⎜-2, 1, ⎢⎢ ⎥⎥⎟, ⎜3, 1, ⎢⎢ ⎥⎥⎟, ⎜5, 2, ⎢⎢ ⎥, ⎢ ⎥⎥⎟⎥
⎢⎜ ⎢⎢1⎥⎥⎟ ⎜ ⎢⎢1⎥⎥⎟ ⎜ ⎢⎢1⎥ ⎢0 ⎥⎥⎟⎥
⎢⎜ ⎢⎢ ⎥⎥⎟ ⎜ ⎢⎢ ⎥⎥⎟ ⎜ ⎢⎢ ⎥ ⎢ ⎥⎥⎟⎥
⎣⎝ ⎣⎣1⎦⎦⎠ ⎝ ⎣⎣1⎦⎦⎠ ⎝ ⎣⎣0⎦ ⎣1 ⎦⎦⎠⎦
This shows us that, for example, the eigenvalue 5 also has geometric multiplicity 2, because it has two eigenvectors. Because the algebraic and geometric multiplicities are the same for all the eigenvalues, M
is diagonalizable.
To diagonalize a matrix, use diagonalize
. diagonalize
returns a tuple (P,D), where D is diagonal and M=PDP−1.
>>> P, D = M.diagonalize()>>> P⎡0 1 1 0 ⎤
⎢ ⎥
⎢1 1 1 -1⎥
⎢ ⎥
⎢1 1 1 0 ⎥
⎢ ⎥
⎣1 1 0 1 ⎦
>>> D⎡-2 0 0 0⎤
⎢ ⎥
⎢0 3 0 0⎥
⎢ ⎥
⎢0 0 5 0⎥
⎢ ⎥
⎣0 0 0 5⎦
>>> P*D*P**-1
⎡3 -2 4 -2⎤
⎢ ⎥
⎢5 3 -3 -2⎥
⎢ ⎥
⎢5 -2 2 -2⎥
⎢ ⎥
⎣5 -2 -3 3 ⎦
>>> P*D*P**-1 == MTrue
lambda
is a reserved keyword in Python, so to create a Symbol called λ, while using the same names for SymPy Symbols and Python variables, use lamda
(without the b
). It will still pretty print as λ.
Note that since eigenvects
also includes the eigenvalues, you should use it instead of eigenvals
if you also want the eigenvectors. However, as computing the eigenvectors may often be costly, eigenvals
should be preferred if you only wish to find the eigenvalues.
If all you want is the characteristic polynomial, use charpoly
. This is more efficient than eigenvals
, because sometimes symbolic roots can be expensive to calculate.
>>> lamda = symbols('lamda')>>> p = M.charpoly(lamda)>>> factor(p.as_expr()) 2
(λ - 5) ⋅(λ - 3)⋅(λ + 2)
十二、其他
1、数学符合补充
#数学符合
#虚数单位i
sympy.I
#自然对数低e
sympy.E
#无穷大
sympy.oo
#圆周率
sympy.pi
#求n次方根
sympy.root(8,3)
#求对数
sympy.log(1024,2)
#求阶乘
sympy.factorial(4)
#三角函数
sympy.sin(sympy.pi)
sympy.tan(sympy.pi/4)
sympy.cos(sympy.pi/2)
2、公式展开与折叠
x=sympy.Symbol('x')
#公式展开用expand方法
f=(1+2*x)*x**2
ff=sympy.expand(f)
print(ff)
#公式折叠用factor方法
f=x**2+1+2*x
ff=sympy.factor(f)
print(ff)
3、公式分离与合并(分数的分离与合并)
x=sympy.Symbol('x')
y=sympy.Symbol('y')
#公式展开用apart方法,和expand区别不是很大,常用于分数进行分离
f=(x+2)/(x+1)
ff=sympy.apart(f)
print(ff)
#公式折叠用tegother方法
f=(1/x+1/y)
ff=sympy.together(f)
print(ff)
4、表达式化简
#simplify( )普通的化简
simplify((x**3 + x**2 - x - 1)/(x**2 + 2*x + 1))
#trigsimp( )三角化简
trigsimp(sin(x)/cos(x))
#powsimp( )指数化简
powsimp(x**a*x**b)
十三、更多...
本文列举的仅仅是常用基本类型的SymPy库中的支持功能,另外还有物理学计算支持,等等。有兴趣的朋友,请自行到官网参考。
十四、引用资料
https://www.jianshu.com/p/339c91ae9f41
https://docs.sympy.org/latest/index.html
http://www.cppcns.com/jiaoben/python/414232.html