SymPy库常用函数
简介
SymPy是一个符号计算的Python库。它的目标是成为一个全功能的计算机代数系统,同时保持代码简 洁、易于理解和扩展。它完全由Python写成,不依赖于外部库。SymPy支持符号计算、高精度计算、模式匹配、绘图、解方程、微积分、组合数学、离散 数学、几何学、概率与统计、物理学等方面的功能。(来自维基百科的描述)
Sympy安装方法
安装命令:pip install sympy
基本数值类型
实数,有理数和整数
SymPy有三个内建的数值类型:实数,有理数和整数。有理数类用两个整数来表示一个有理数。分子与分母,所以Rational(1,2)代表1/2,Rational(5,2)代表5/2,等等。
from sympy import *
a = Rational(1,2)
a
1/2
a*2
1
Rational(2)**50/Rational(10)**50
1/88817841970012523233890533447265625
当利用Python的整数计算时要注意一下,Python只会截取除法的整数部分:
1/2
0
1.0/2
0.5
然而你可以:
from future import division
1/2 #doctest: +SKIP
0.5
正确的除法在python3k和isympy中这样做,是标准的。
特殊的常数
我们也可以有一些特殊的常数,像e和pi,它们会被当作符号去对待。(1+pi不会求得值,反而它会保持为1+pi),例如:
pi**2
pi**2
pi.evalf()
3.14159265358979
(pi+exp(1)).evalf()
5.85987448204884
求表达式的浮点数-evalf()函数
正如你看到的,evalf()函数可以用求出表达式的浮点数。
有一个无穷大的类型,被成为oo:
oo > 99999
True
oo + 1
oo
If the substitution will be followed by numerical evaluation, it is better to pass the substitution to evalf as
(1/x).evalf(subs={x: 3.0}, n=21)
0.333333333333333333333
rather than
(1/x).subs({x: 3.0}).evalf(21)
0.333333333333333314830
Sympy基本使用
定义变量-Symbols函数
对比与其他的计算机代数系统,在SymPy中要明确声明符号变量:
x = symbols(‘x’)
x + 1
x + 1
x,y,z=symbols(‘x y z’)
crazy = symbols(‘unrelated’)
crazy + 1
unrelated + 1
x = symbols(‘x’)
expr = x + 1
x = 2
print(expr)
x + 1
Changing x to 2 had no effect on expr. This is because x = 2 changes the Python variable x to 2, but has no effect on the SymPy Symbol x, which was what we used in creating expr.
变量替换subs函数
x = symbols(‘x’)
expr = x + 1
expr.subs(x, 2)
3
from sympy import pi, exp, limit, oo
from sympy.abc import x, y
(1 + x*y).subs(x, pi)
pi*y + 1
(1 + x*y).subs({x:pi, y:2})
1 + 2*pi
(1 + x*y).subs([(x, pi), (y, 2)])
1 + 2*pi
reps = [(y, x**2), (x, 2)]
(x + y).subs(reps)
6
(x + y).subs(reversed(reps))
x**2 + 2
(x**2 + x**4).subs(x**2, y)
y**2 + y
(x**2 + x**4).xreplace({x**2: y})
x**4 + y
(x/y).subs([(x, 0), (y, 0)])
0
(x/y).subs([(x, 0), (y, 0)], simultaneous=True)
nan
((x + y)/y).subs({x + y: y, y: x + y})
1
((x + y)/y).subs({x + y: y, y: x + y}, simultaneous=True)
y/(x + y)
limit(x**3 - 3*x, x, oo)
oo
调用方式:[subs(*args, **kwargs)]
代数
局部的代数式展开,使用apart(expr, x):
In [1]: 1/( (x+2)*(x+1) )
Out[1]:
1
───────────────
(2 + x)*(1 + x)
In [2]: apart(1/( (x+2)*(x+1) ), x)
Out[2]:
1 1
───── - ─────
1 + x 2 + x
In [3]: (x+1)/(x-1)
Out[3]:
-(1 + x)
────────
1 - x
In [4]: apart((x+1)/(x-1), x)
Out[4]:
2
1 - ─────
1 - x
代数式的合并
(相当于展开的逆运算),使用together(expr, x):
In [7]: together(1/x + 1/y + 1/z)
Out[7]:
x*y + x*z + y*z
───────────────
x*y*z
In [8]: together(apart((x+1)/(x-1), x), x)
Out[8]:
-1 - x
──────
1 - x
In [9]: together(apart(1/( (x+2)*(x+1) ), x), x)
Out[9]:
1
───────────────
(2 + x)*(1 + x)
微积分
极限
在sympy中极限容易求出,它们遵循极限语法 limit(function, variable, point) ,所以计算x->0时f(x)的极限,即limit(f, x, 0):
from sympy import *
x=Symbol(“x”)
limit(sin(x)/x, x, 0)
1
limit(x, x, oo)
oo
limit(1/x, x, oo)
0
limit(x**x, x, 0)
1
有一些特殊的极限的例子,可以阅读文件test_demidovich.py
微分
可以对任意SymPy表达式微分。diff(func, var)。例如:
from sympy import *
x = Symbol(‘x’)
diff(sin(x), x)
cos(x)
diff(sin(2*x), x)
2*cos(2*x)
diff(tan(x), x)
1 + tan(x)**2
可以通过以下验证:
limit((tan(x+y)-tan(x))/y, y, 0)
1 + tan(x)**2
计算高阶微分 diff(func, var, n) :
diff(sin(2*x), x, 1)
2*cos(2*x)
diff(sin(2*x), x, 2)
-4*sin(2*x)
diff(sin(2*x), x, 3)
-8*cos(2*x)
级数展开
函数 series(var, point, order):
from sympy import *
x = Symbol(‘x’)
cos(x).series(x, 0, 10)
1 - x**2/2 + x**4/24 - x**6/720 + x**8/40320 + O(x**10)
(1/cos(x)).series(x, 0, 10)
1 + x**2/2 + 5*x**4/24 + 61*x**6/720 + 277*x**8/8064 + O(x**10)
积分
SymPy支持不定积分,超越函数与特殊函数的定积分。SymPy有力的扩展Risch-Norman 算法和模型匹配算法。
from sympy import *
x, y = symbols(‘xy’)
初等函数:
integrate(6*x**5, x)
x**6
integrate(sin(x), x)
-cos(x)
integrate(log(x), x)
-x + x*log(x)
integrate(2*x + sinh(x), x)
cosh(x) + x**2
特殊函数:
integrate(exp(-x**2)*erf(x), x)
pi**(1/2)*erf(x)**2/4
定积分:
integrate(x**3, (x, -1, 1))
0
integrate(sin(x), (x, 0, pi/2))
1
integrate(cos(x), (x, -pi/2, pi/2))
2
一些广义积分也可以被支持:
integrate(exp(-x), (x, 0, oo))
1
integrate(log(x), (x, 0, 1))
-1
复数
from sympy import Symbol, exp, I
x = Symbol(“x”)
exp(I*x).expand()
exp(I*x)
exp(I*x).expand(complex=True)
I*exp(-im(x))*sin(re(x)) + cos(re(x))*exp(-im(x))
x = Symbol(“x”, real=True)
exp(I*x).expand(complex=True)
I*sin(x) + cos(x)
函数
三角函数::
In [1]: sin(x+y).expand(trig=True)
Out[1]: cos(x)*sin(y) + cos(y)*sin(x)
In [2]: cos(x+y).expand(trig=True)
Out[2]: cos(x)*cos(y) - sin(x)*sin(y)
In [3]: sin(I*x)
Out[3]: I*sinh(x)
In [4]: sinh(I*x)
Out[4]: I*sin(x)
In [5]: asinh(I)
Out[5]:
π*I
───
2
In [6]: asinh(I*x)
Out[6]: I*asin(x)
In [15]: sin(x).series(x, 0, 10)
Out[15]:
3 5 7 9
x x x x
x - ── + ─── - ──── + ────── + O(x**10)
6 120 5040 362880
In [16]: sinh(x).series(x, 0, 10)
Out[16]:
3 5 7 9
x x x x
x + ── + ─── + ──── + ────── + O(x**10)
6 120 5040 362880
In [17]: asin(x).series(x, 0, 10)
Out[17]:
3 5 7 9
x 3*x 5*x 35*x
x + ── + ──── + ──── + ───── + O(x**10)
6 40 112 1152
In [18]: asinh(x).series(x, 0, 10)
Out[18]:
3 5 7 9
x 3*x 5*x 35*x
x - ── + ──── - ──── + ───── + O(x**10)
6 40 112 1152
球谐函数:
In [1]: from sympy.abc import theta, phi
In [2]: Ylm(1, 0, theta, phi)
Out[2]:
————
╲╱ 3 *cos(θ)
────────────
——
2*╲╱ π
In [3]: Ylm(1, 1, theta, phi)
Out[3]:
—— I*φ
-╲╱ 6 │sin(θ)│ℯ
────────────────────
——
4*╲╱ π
In [4]: Ylm(2, 1, theta, phi)
Out[4]:
——— I*φ
-╲╱ 30 │sin(θ)│*cos(θ)ℯ
────────────────────────────
——
4*╲╱ π
阶乘和伽玛函数:
In [1]: x = Symbol(“x”)
In [2]: y = Symbol(“y”, integer=True)
In [3]: factorial(x)
Out[3]: Γ(1 + x)
In [4]: factorial(y)
Out[4]: y!
In [5]: factorial(x).series(x, 0, 3)
Out[5]:
2 2 2 2
x *EulerGamma π *x
1 - x*EulerGamma + ────────────── + ───── + O(x**3)
2 12
Zeta函数:
In [18]: zeta(4, x)
Out[18]: ζ(4, x)
In [19]: zeta(4, 1)
Out[19]:
4
π
──
90
In [20]: zeta(4, 2)
Out[20]:
4
π
-1 + ──
90
In [21]: zeta(4, 3)
Out[21]:
4
17 π
- ── + ──
16 90
多项式
In [1]: chebyshevt(2, x)
Out[1]:
2
-1 + 2*x
In [2]: chebyshevt(4, x)
Out[2]:
2 4
1 - 8*x + 8*x
In [3]: legendre(2, x)
Out[3]:
2
3*x
-1/2 + ────
2
In [4]: legendre(8, x)
Out[4]:
2 4 6 8
35 315*x 3465*x 3003*x 6435*x
─── - ────── + ─────── - ─────── + ───────
128 32 64 32 128
In [5]: assoc_legendre(2, 1, x)
Out[5]:
—————
╱ 2
-3*x*╲╱ 1 - x
In [6]: assoc_legendre(2, 2, x)
Out[6]:
2
3 - 3*x
In [7]: hermite(3, x)
Out[7]:
3
-12*x + 8*x
微分方程
在isympy中:
In [4]: f(x).diff(x, x) + f(x) #注意在使用输入该命令之前,一定要声明f=Function(‘f’)
Out[4]:
2
d
─────(f(x)) + f(x)
dx dx
In [5]: dsolve(f(x).diff(x, x) + f(x), f(x))
Out[5]: C₁*sin(x) + C₂*cos(x)
代数方程
在isympy中:
In [7]: solve(x**4 - 1, x)
Out[7]: [i, 1, -1, -i]
In [8]: solve([x + 5*y - 2, -3*x + 6*y - 15], [x, y])
Out[8]: {y: 1, x: -3}
线性代数
矩阵
矩阵由矩阵类创立建:
from sympy import Matrix
Matrix([[1,0], [0,1]])
[1, 0]
[0, 1]
不只是数值矩阵,亦可为代数矩阵,即矩阵中存在符号:
x = Symbol(‘x’)
y = Symbol(‘y’)
A = Matrix([[1,x], [y,1]])
A
[1, x]
[y, 1]
A**2
[1 + x*y, 2*x]
[ 2*y, 1 + x*y]
关于矩阵更多的例子,请看线性代数教程。
系数匹配
使用 .match()方法,引用Wild类,来执行表达式的匹配。该方法会返回一个字典。
from sympy import *
x = Symbol(‘x’)
p = Wild(‘p’)
(5*x**2).match(p*x**2)
{p_: 5}
q = Wild(‘q’)
(x**2).match(p*x**q)
{p_: 1, q_: 2}
如果匹配不成功,则返回None:
print (x+1).match(p**x)
None
可以使用Wild类的‘exclude’参数(排除参数),排除不需要和无意义的匹配结果,来保证结论中的显示是唯一的:
x = Symbol(‘x’)
p = Wild(‘p’, exclude=[1,x])
print (x+1).match(x+p) # 1 is excluded
None
print (x+1).match(p+1) # x is excluded
None
print (x+1).match(x+2+p) # -1 is not excluded
{p_: -1}
打印输出
标准
str(expression)返回如下:
from sympy import Integral
from sympy.abc import x
print x**2
x**2
print 1/x
1/x
print Integral(x**2, x)
Integral(x**2, x)
Pretty Printing
用pprint函数可以输出不错的ascii艺术:
from sympy import Integral, pprint
from sympy.abc import x
pprint(x**2) #doctest: +NORMALIZE_WHITESPACE
2
x
pprint(1/x)
1
x
pprint(Integral(x**2, x))
/
|
| 2
| x dx
|
/
[Pretty PrintingWiki]
提示:在python解释器中,为使pretty printing为默认输出,使用:
$ python
Python 2.5.2 (r252:60911, Jun 25 2008, 17:58:32)
[GCC 4.3.1] on linux2
Type “help”, “copyright”, “credits” or “license” for more information.
from sympy import *
import sys
sys.displayhook = pprint
var(“x”)
x
x**3/3
3
x
3
Integral(x**2, x) #doctest: +NORMALIZE_WHITESPACE
/
|
| 2
| x dx
|
/
Python printing
from sympy.printing.python import python
from sympy import Integral
from sympy.abc import x
print python(x**2)
x = Symbol(‘x’)
e = x**2
print python(1/x)
x = Symbol(‘x’)
e = 1/x
print python(Integral(x**2, x))
x = Symbol(‘x’)
e = Integral(x**2, x)
LaTeX printing
from sympy import Integral, latex
from sympy.abc import x
latex(x**2)
x2x2
latex(1/x)
1x1x
latex(Integral(x**2, x))
∫x2dx∫x2dx
MathML
from sympy.printing.mathml import mathml
from sympy import Integral, latex
from sympy.abc import x
print mathml(x**2)
x2
print mathml(1/x)
x-1
Pyglet
from sympy import Integral, preview
from sympy.abc import x
preview(Integral(x**2, x)) #doctest:+SKIP
注解
Isympy默认调用pprint,所以这就是为什么看到pretty printing为默认的。
有一个打印的有效模块,sympy.printing。用这个模块实现其他的打印:
pretty(expr), pretty_print(expr), pprint(expr): 分别返回或者输出,,表达式的漂亮描述。这是相同
latex(expr), print_latex(expr):分别返回或者输出,LaTex描写的表达式
mathml(expr), print_mathml(expr):分别返回或者输出,MathML描写的表达式
print_gtk(expr): 表达式打印到Gtkmathview , 这是一个GTK小配件显示MathML代码。Gtkmathview程序是必须的。
【例】
#例1:解下列二元一次方程
# 2x-y=3
# 3x+y=7
import sympy
#将变量x,y符号化
x = sympy.Symbol('x')
y = sympy.Symbol('y')
#print(type(x))查看x类型
#以上代码也可写作
#x, y = sympy.symbols('x y')
#表达式写法
#加号 +
#减号 -
#除号 /
#乘号 *
#指数 **
#对数 sympy.log()
#e的指数次幂 sympy.exp()
#解方程用solve函数,解一元方程用法
#solve(x*2-4,x),solve返回类型为dict
i=sympy.solve([2 * x - y - 3, 3 * x + y - 7],[x, y])
print (i)
#print(type(i))
#运行结果为{y: 1, x: 2}
#例2:limit函数求极限
#[(n+3)/(n+2)]^n ,n趋紧无穷大,求极值
#print ( sympy.limit(((x+3)/(x+2))**x, x, sympy.oo) )
#以上代码也可写作
n = sympy.Symbol('n')
s = ((n+3)/(n+2))**n
#print(type(s))
print ( sympy.limit(s, n, sympy.oo) )
#符号介绍:
#sympy.oo 无穷大(标识方式是两个小写字母o连接在一起)
#sympy.E e
#sympy.pi 圆周率
#例3:不定积分
x = sympy.Symbol('x')
print( sympy.integrate(sympy.cos(x),x) )
#例4: 定积分
x, t = sympy.symbols('x t')
m = sympy.integrate( sympy.sin(t)/(sympy.pi-t) , (t,0,x) )
#t表示积分变量,0表示积分下限,x表示积分上限
print(m)
#例5: 微分函数
x = sympy.Symbol('x')
print( sympy.diff(sympy.sin(x),x) )
print( sympy.diff(sympy.sin(x),x,2) )#2阶导数
#例6:解微分方程 y'=2xy
f = sympy.Function('f')
x = sympy.Symbol('x')
print ( sympy.dsolve(sympy.diff(f(x),x) - 2*f(x)*x,f(x)) )
#例7: 矩阵化简
# (a11 a12 a13) (x1)
# (x1,x2,x3)*|a12 a22 a23|*|x2|
# (a13 a23 a33) (x3)
x1,x2,x3 = sympy.symbols('x1 x2 x3')
a11,a12,a13,a22,a23,a33 = sympy.symbols('a11 a12 a13 a22 a23 a33')
m = sympy.Matrix([[x1,x2,x3]])
n = sympy.Matrix([[a11,a12,a13],[a12,a22,a23],[a13,a23,a33]])
v = sympy.Matrix([[x1],[x2],[x3]])
f = m * n * v
print (f)
print ( f[0].subs({x1:1, x2:1, x3:1}) ) #令x1,x2,x3都为1
Sympy特点—————数学符号计算
import math
math.sqrt(8)
2.8284271247461903
我们看看Python中结果
math.sqrt(8).math.sqrt(8)
8.000000000000002
本以为会得到8.0,但没想到得到8.000000000000002。
一、为什么会这样?
如果我们平常计算的任务常常有类似于上面的例子这样的表达式,那么直接用python计算其结果只是真实值的逼近。如果这样的计算很大很多,误差会逐渐积累,这是我们不能忍受的,所以这时候就需要Python能处理这种数学符号计算。
二、什么是数学符号计算?
数学符号计算能处理表征数字的符号计算。这意味着数学对象被精确地表示,而不是近似地表示,而具有未被计算的变量的数学表达式被留在符号形式中。
sympy库简介
Sympy是Python的一个数学符号计算库。它目的在于成为一个富有特色的计算机代数系统。它保证自身的代码尽可能的简单,且易于理解,容易扩展。Sympy完全由Python写成,不需要额外的库。
sympy的表达式与我们平常的手写的数学表达式略微有所区别,下面是sympy的方程表示符号
- 加号 +
- 减号 -
- 除号 /
- 乘号 *
- 等号 Eq()
- 指数 **
- 对数 log()
- e的指数次幂 exp()
上面的例子我们用Python实现一下。
import sympy
sympy.sqrt(8)
2*sqrt(2)
用sympy计算
sympy.sqrt(8)*sympy.sqrt(8)
8
三、简单学一下sympy中的几个实例
- 定义数学符号(类似于数学中的变量)
- 展开与折叠
- 简化表达式
- 解方程
- 赋值计算
- log计算
- 导数
- 积分
- 求极限
3.1 定义数学符号
让我们定义一个符号表达式代表数学表达式 x+2yx+2y。首先我们要注意到python中的变量必须赋值才能使用,所以无法表达该数学表达式。所以这里一定要引入特殊的符号,这里有两种方法
- 方法一
from sympy import symbols
x,y = symbols('x y')
expr = x + 2*y
expr
x + 2*y
- 方法二
from sympy.abc import x,y
expr2 = x + 2*y
expr2
x + 2*y
**当数学表达式中的变量不是x,y这种单一字符,而是result这种多个字符长度的变量时,只能用方法一。
3.2 展开与折叠
from sympy import expand,factor
from sympy.abc import x,y
expr = x**2+x*y+3*x
expr
x**2 + x*y + 3*x
- 折叠
factor(expr)
x**2 + x*y + 3*x
- 展开
expr2 = x*(x+y+3)
expand(expr2)
x**2 + x*y + 3*x
3.3 简化表达式
有时候我们需要简化表达式
- 普通的化简
from sympy import simplify
from sympy.abc import x
simplify((x**3 + x**2 - x - 1)/(x**2 + 2*x + 1))
x - 1
- 三角化简trigsimp
from sympy import trigsimp,sin,cos
from sympy.abc import x,y
y = sin(x)/cos(x)
trigsimp(y)
tan(x)
- 指数化简
from sympy import powsimp
from sympy.abc import x,a,b
y = x**a * x**b
y
x**a*x**b
#指数化简
powsimp(y)
x**(a + b)
3.4 解方程
注意在python中=是赋值的意思,==虽然表示等于,但是会有很大的问题。在sympy中,我们使用Eq(x,y)表示x=y
from sympy.abc import x,y
from sympy import solve,linsolve,Eq
#对一个方程求解,使用solve
solve(Eq(2*x-1,3), x)
[2]
使用linsolve([方程1,方程2,...],(变量1,变量2,...))
#对多个方程求解,使用linsolve。方程的解为x=-1,y=3
linsolve([x+2*y-5,2*x+y-1], (x,y))
{(-1, 3)}
3.5 赋值计算
from sympy.abc import x,y
from sympy import sin,cos
y = sin(x)+cos(x)
y
sin(x) + cos(x)
y.subs(x, x**2)
sin(x**2) + cos(x**2)
这里的赋值,不仅可以实现变量的替换,还可以赋与数字,进行计算。
y.subs(x, 0)
1
3.6 log运算
from sympy import log,expand_log
from sympy.abc import x,y,e
#expand_log为展开log,但需要将force=True,展开才能发生
expand_log(log(x**3), force=True)
3*log(x)
#expand_log为展开log,但需要将force=True,展开才能发生
expand_log(log(x**3))
log(x**3)
expand_log(log(e**x), force=True)
x*log(e)
3.7 导数
from sympy import diff,sin,cos
from sympy.abc import x,y,z,f
#对sin(x)求导
diff(sin(x))
cos(x)
diff(cos(x))
-sin(x)
偏导
#求偏导
f = 3*x**2*y*z
diff(f, x,y)
6*x*z
# 对f关于x_1求1次导 (一阶偏导)
diff(f, x_1, 1)
# 先对x_1求偏导,再对x_2求偏导(二阶偏导)
diff(f, x_1,x_2)
# 先对x_1求偏导,再对x_2求两次偏导(二阶偏导)
diff(f, x_1,x_2,2)
3.8 积分
from sympy.abc import pi,x
from sympy import integrate,sin
integrate(sin(x), (x,0,pi))
-cos(pi) + 1
3.9 极限
from sympy.abc import x
from sympy import limit
limit(1/x, x, 0, '+')
oo
3.10 展开式
高数中有泰勒展开式,拉格朗日展开式。
e^x=1+x+x^2/2!+x^3/3!+x^4/4!+...+x^n/n!+o(x^n)
比如当n=3时,
e^x=1+x+x^2/2+o(x^3)
这里实现的方法是:sympy表达式.series(变量, 0, n)
from sympy import exp,symbols
x = symbols('x')
expr = exp(x)
expr.series(x, 0, 3)
1 + x + x**2/2 + O(x**3)