赵工的个人空间


专业技术部分转网页计算转业余爱好部分


 编程语言

常用的编程语言
C#编程语言基础
C#面向对象与多线程
C#数据及文件操作
JavaScript基础
JavaScript的数据类型和变量
JavaScript的运算符和表达式
JavaScript的基本流程控制
JavaScript的函数
JavaScript对象编程
JavaScript内置对象和方法
JavaScript的浏览器对象和方法
JavaScript访问HTML DOM对象
JavaScript事件驱动编程
JavaScript与CSS样式表
Ajax与PHP
ECMAScript6的新特性
Vue.js前端开发
PHP的常量与变量
PHP的数据类型与转换
PHP的运算符和优先规则
PHP程序的流程控制语句
PHP的数组操作及函数
PHP的字符串处理与函数
PHP自定义函数
PHP的常用系统函数
PHP的图像处理函数
PHP类编程
PHP的DataTime类
PHP处理XML和JSON
PHP的正则表达式
PHP文件和目录处理
PHP表单处理
PHP处理Cookie和Session
PHP文件上传和下载
PHP加密技术
PHP的Socket编程
PHP国际化编码
MySQL数据库基础
MySQL数据库函数
MySQL数据库账户管理
MySQL数据库基本操作
MySQL数据查询
MySQL存储过程和存储函数
MySQL事务处理和触发器
PHP操作MySQL数据库
数据库抽象层PDO
Smarty模板
ThinkPHP框架
Python语言基础
Python语言结构与控制
Python的函数和模块
Python的复合数据类型
Python面向对象编程
Python的文件操作
Python的异常处理
Python的绘图模块
Python的NumPy模块
Python的SciPy模块
Python的SymPy模块
Python的数据处理
Python操作数据库
Python网络编程
Python图像处理
Python机器学习
TensorFlow深度学习
Tensorflow常用函数
TensorFlow用于卷积网络
生成对抗网络GAN


首页 > 专业技术 > 编程语言 > Python的SymPy模块
Python的SymPy模块

计算机直接对数学符号进行正确的计算称为符号计算。符号计算也称为计算机代数。SymPy是Python版的开源计算机代数系统实现,是使用纯Python代码,没有第三方库。SymPy开发是Ondrej Certik从2006年8月开始的,此后不断有开发者加入项目,规模达到几百人。现在这个程序库包括26个模块,可以满足常用的计算需求,如符号计算、积分、代数、离散数学、量子物理、画图与打印等,计算结果输出为LaTeX或其他格式。
SymPy程序库分为一个核心模块和多个高级可选模块:
·Assumptions:假设引擎
·Concrete:符号积和符号总和
·Core basic class structure:基本的,及加、乘、指数等
·Functions:基本的函数和特殊的函数
·Galgebra:几何代数
·Geometry:几何实体
·Integrals:符号积分
·Interactive:交互会话
·Logic:布尔代数和定理证明
·Matrices:线性代数和矩阵
·mpmath:快速的任意精度的数值运算
·ntheory:数论函数
·Parsing:数学的和最大化的句法分析
·Physics:物理单位和量子相关
·Plotting:用Pyglet进行二维和三维画图
·Polys:多项式代数和因式分解
·Printing:漂亮的打印和代码生成
·Series:符号极限和截断的序列
·Simplify:用其他形式改写表达式
·Solvers:代数、循环和差分
·Statistics:标准概率分布
·Utilities:测试架构和兼容性相关的内容
SymPy包括很多功能,从基本符号算术到多项式、微积分、求解方程、离散数学、几何、统计和物理,主要处理整型数据、实数和有理数三种类型数据,整数是不带小数点的数字,实数是带小数点的数字,有理数包括分子和分母,用Ration类定义有理数,该类需要两个数字。SymPy的核心功能是基本的算术、扩展、简化、替换、模式匹配和各种函数。

1. 符号、表达式和基本运算:

1)符号的定义:
在SymPy中,在任何表达式中使用符号前,必须先定义该符号,定义符号只需要用Symbol类中的symbol来定义一个符号即可。示例:
from __future__ import division
from sympy import *
x,y,z=symbols('x y x')
m0,m1,m2,m3,m4=symbols('m0:5')
x1=Symbol('x1')
print(x1+500)
y=22/7
代码中使用symbol只生成一个符号;生成多个符号使用symbols,可以用空格分隔多个符号名,也可以使用冒号:生成一个系列,冒号前是第一个值,冒号后是生成符号的数量。
2)将SymPy对象的数值转换为近似浮点值:
可以用evalf()和n()来获得任何对象的浮点近似值,默认的精度是15位有效数字,而且可以通过调整参数改为任何想要的精度。示例:
from __future__ import division
from sympy import sin,pi
x=sin(50)
print(pi.evalf())
print(pi.evalf(50)) # 50位有效数字
print(x.n())
print(x.n(20)) # 20位有效数字
3)表达式的常用操作:
表达式可以使用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=5*x**2+2*b*x**2+cos(x)+51*x**2
simplify(expr)
factor(x**2+x-30)
expand((x-5)*(x+6))
collect(x**3+a*x**2+b*x**2+c*x+d,x)
expr=sin(x)*sin(x)+cos(x)*cos(x)
print(expr)
print(expr.subs({x:5,y:25}))
print(expr.subs({x:5,y:25}).n())

2. 求解方程:

solve可以求解各种类型的方程,需要待解的表达式和变量两个输入参数。示例:
from sympy import solve,symbols
a,b,c,x,y=symbols('a b c x y')
solve(6*x**2-3*x-30,x)
solve(a*x**2+b*x+c,x)
substitute_solution=solve(a*x**2+b*x+c,x)
print([substitute_solution[0].subs({'a':6,'b':-3,'c':-30}),
 substitute_solution[1].subs({'a':6,'b':-3,'c':-30})]) # [5/2, -2]
print(solve([2*x+3*y-3,x-2*y+1],[x,y])) # {x: 3/7, y: 5/7}
还有另一种形式的solve方法,将一系列方程作为第一个输入参数,将未知数列表作为第二个参数:
from sympy import solve,symbols
x,y=symbols('x y')
print(solve([2*x+y-4,5*x-3*y],[x,y])) # {x: 12/11, y: 20/11}
print(solve([2*x+2*y-1,2*x-4*y],[x,y])) # {x: 1/3, y: 1/6}

3. 有理数、指数和对数函数:

SymPy有很多函数可以用于处理有理数,可以对有理数做简化、扩展、合并、拆分等操作,用于计算两个有理数的加法使用together函数,计算有理数的除法用到apart函数。SymPy还支持一些指数和对数操作,log用于计算以b为底的对数,ln用于计算以e为底的自然对数,log10用于计算以10为底的对数。
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
print(together(x1/x2+x3/x4)) # (x1*x4 + x2*x3)/(x2*x4)
print(apart((2*x**2+3*x+4)/(x+1))) # 2*x + 1 + 3/(x + 1)
print(together(apart((2*x**2+3*x+4)/(x+1)))) # (2*x*(x + 1) + x + 4)/(x + 1))
print(exp(1))
print(log(4).n())
print(log(4,4).n())
print(ln(4).n())
print(mpmath.log10(4))

4. 多项式:

SymPy允许定义和执行多项式的各种操作,还可以求解多项式的根。为了检查两个多项式是否相等,要使用simplify函数。
from sympy import *
p,q,x=symbols('p q x')
p=(x+4)*(x+2)
q=x**2+6*x+8
print(simplify(p-q)==0) # True

5. 三角函数和复数:

模块提供了从角度到弧度的转换,除了sin、cos、tan这些基本的三角函数,还有简化和扩展的三角函数。SymPy还支持在不存在实数解时的复数解。
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)
print(trigsimp(expr)) # sin(x + y)
print(trigsimp(expr_exp)) # exp(sin(x)**2 + 4)
print(expand_trig(sin(x+y))) # sin(x)*cos(y) + sin(y)*cos(x)
print(solve(x**2+4,x)) # [-2*I, 2*I]

6. 线性代数:

SymPy线性代数模块中的函数包括对矩阵的各种操作,例如快速特殊矩阵的构建、特征值、特征向量、转置、行列式的值和逆。有三个快速生成特殊矩阵的函数,sye、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)
print(A[0,1])
print(A[0:2,0:3])
B=Matrix([[1,2,3],[5,6,7],[9,10,11]])
A.row_join(B)
B.col_join(B)
print(A+B)
print(A-B)
print(A*B)
print(A**2)
print(eye(3))
print(zeros(3,3))
print(ones(3,3))
print(A.transpose())
M=Matrix([[1,2,3],[4,5,6],[7,8,10]])
print(A.det())
默认情况下,矩阵的逆都是用高斯消元法计算得出,也可以指定用LU分解法来计算。SymPy有很多方法可以用来计算简化的行阶梯型(rref函数)和零空间(nullspace函数)。如果A是一个矩阵,那么nullspace就是所有向量v的集合,满足Av=0,还可以对矩阵元素做替换操作。
from sympy import *
A=Matrix([[1,2],[3,4]])
A.inv()
print(A.inv()*A)
print(A*A.inv())
A=Matrix([[1,-2],[-2,3]])
print(A.eigenvals())               # solve(det(A-eye(2)*x),x)
print(A.eigenvects())
print(A.rref())
print(A.nullspace())
x=Symbol('x')
M=eye(3)*x
print(M.subs(x,4))
y=Symbol('y')
print(M.subs(x,y))
print(M.inv())
print(M.inv('LU'))
A=Matrix([[1,2,1],[2,3,3],[1,3,2]])
Q,R=A.QRdecomposition()
print(Q)
M=[Matrix([1,2,3]),Matrix([3,4,5]),Matrix([5,7,8])]
result1=GramSchmidt(M)
result2=GramSchmidt(M,True)

7. 微积分:

微积分操作包括极限、导数、序列的求和以及积分。求极限问题的示例:
from sympy import limit,oo,symbols,exp,sin,cos
print(oo+5) # oo
print(6700<oo) # True
print(10/oo) # 0
x,n=symbols('x n')
print(limit(((x**n-1)/(x-1)),x,1)) # n
print(limit(1/x**2,x,0)) # oo
print(limit(1/x,x,0,dir='-')) # -oo
print(limit(cos(x)/x,x,0)) # oo
print(limit(sin(x)**2/x,x,0)) # 0
print(limit(exp(x)/x,x,oo)) # oo
任何SymPy表达式都可以用diff函数进行微分,两个参数分别是表达式和变量。示例:
from sympy import diff,Symbol,symbols,exp,dsolve,cos,sin,Function
x=symbols('x')
print(diff(x**4,x)) # 4*x**3
print(diff(x**3*cos(x),x)) # -x**3*sin(x) + 3*x**2*cos(x)
print(diff(cos(x**2),x)) # -2*x*sin(x**2)
print(diff(x/sin(x),x)) # -x*cos(x)/sin(x)**2 + 1/sin(x)
print(diff(x**3,x,2)) # 6*x
print(diff(exp(x),x)) # exp(x)
dsolve函数可用来求解常微分方程和常微分方程组:
from sympy import diff,Symbol,symbols,exp,dsolve,cos,sin,Function
x=symbols('x')
f=symbols('f',cls=Function)
print(dsolve(f(x)-diff(f(x),x),f(x))) # Eq(f(x), C1*exp(x))
可以使用ics参数确定dsolve方法的边界:
from sympy import *
x=symbols('x')
f=symbols('f',cls=Function)
print(dsolve(Eq(f(x).diff(x,x),400*(f(x)-1+2*x)),f(x),ics={f(0):5,f(2):10}))
计算结果为:Eq(f(x), C1*exp(-20*x) + C2*exp(20*x) - 2*x + 1)
下面的代码用于寻找f(x)=4*x**3-3*x**2+2*x的临界点,并用二次导数来找该函数在[0,1]区间的最大值:
from sympy import *
x=symbols('x')
f=4*x**3-3*x**2+2*x
print(diff(f,x)) # 12*x**2 - 6*x + 2
sols=solve(diff(f,x),x)
print(sols) # [1/4 - sqrt(15)*I/12, 1/4 + sqrt(15)*I/12]
print(diff(diff(f,x),x).subs({x:sols[0]})) # -2*sqrt(15)*I
print(diff(diff(f,x),x).subs({x:sols[1]})) # 2*sqrt(15)*I
SymPy中可以用integrate函数计算定积分和不定积分。
from sympy import *
x=symbols('x')
print(integrate(x**3+1,x)) # x**4/4 + x
print(integrate(x*sin(x),x)) # -x*cos(x) + sin(x)
print(integrate(x+ln(x),x)) # x**2/2 + x*log(x) - x
F=integrate(x**3+1,x)
print(F.subs({x:1})-F.subs({x:0})) # 5/4
print(integrate(x**3-x**2+x,(x,0,1))) # 5/12
print(integrate(sin(x)/x,(x,0,pi))) # Si(pi)
print(integrate(sin(x)/x,(x,pi,2*pi))) # -Si(pi) + Si(2*pi)
print(integrate(x*sin(x)/(x+1),(x,0,2*pi))) #
print(integrate(x*sin(x)/(x+1),(x,0,2*pi)).n())
序列是将整数作为输入的函数,可以通过为其第n个值指定一个表达式来定义一个序列,还可以替换其中的某个值。
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)]
print(limit(s1_n,n,oo)) # 0
print(limit(s2_n,n,oo)) # 0
下面使用一些特殊函数计算一些级数。
from sympy import *
x,n=symbols('x n')
s1_n=1/n
s2_n=1/factorial(n)
print(summation(s1_n,[n,1,oo])) # oo
print(summation(s2_n,[n,1,oo])) # -1 + E
import math
def s2_nf(n):
    return 1.0/math.factorial(n)
def exponential_xnf(x,n):
    return x**n/math.factorial(n)
print(sum([s2_nf(n) for n in range(0,10)])) # 2.7182815255731922
print(E.evalf()) # 2.71828182845905
exponential_xn=x**n/factorial(n)
print(summation(exponential_xn.subs({x:5}),[x,0,oo]).evalf()) # inf*5.0**n/factorial(n)
print(exp(5).evalf()) # 148.413159102577
print(summation(exponential_xn.subs({x:5}),[x,0,oo])) # oo*5**n/factorial(n)
print(sum([exponential_xnf(5.0,i) for i in range(0,35)])) # 148.41315910257657
print(series(sin(x),x,0,8)) # x - x**3/6 + x**5/120 - x**7/5040 + O(x**8)
print(series(cos(x),x,0,8)) # 1 - x**2/2 + x**4/24 - x**6/720 + O(x**8)
print(series(sinh(x),x,0,8)) # x + x**3/6 + x**5/120 + x**7/5040 + O(x**8)
print(series(cosh(x),x,0,8)) # 1 + x**2/2 + x**4/24 + x**6/720 + O(x**8)
print(series(ln(x),x,1,6)) # -1 - (x - 1)**2/2 + (x - 1)**3/3 - (x - 1)**4/4 + (x - 1)**5/5 + x + O((x - 1)**6, (x, 1))
print(series(ln(x+1),x,0,6)) # x - x**2/2 + x**3/3 - x**4/4 + x**5/5 + O(x**6)

8. 向量:

一个包括实数值的n元组就是一个向量。在物理和数学中,一个向量是一个数学对象,该对象拥有大小、幅度或长度,以及方向。在SymPy中,一个向量表示为一个1xn的行矩阵或nx1的列矩阵。向量可以计算转置和范数:
from sympy import *
u=Matrix([[1,2,3]])
v=Matrix([[4],[5],[6]])
print(v.T) # Matrix([[4, 5, 6]])
print(u[1]) # 2
print(u.norm()) # sqrt(14)
uhat=u/u.norm()
print(uhat) # Matrix([[sqrt(14)/14, sqrt(14)/7, 3*sqrt(14)/14]])
print(uhat.norm()) # 1
向量可以计算点积、叉积和向量的投影:
from sympy import *
u=Matrix([[1,2,3]])
v=Matrix([[-2,3,3]])
print(u.dot(v)) # 13
print(acos(u.dot(v)/(u.norm())).evalf()) # 1.91718313311192*I
print(u.dot(v)==v.dot(u)) # True
u=Matrix([[2,3,4]])
n=Matrix([[2,2,3]])
w=(u.dot(n)/n.norm()**2)*n
v=u-((u.dot(n)/n.norm()**2)*n)
print(w) # Matrix([[44/17, 44/17, 66/17]])
print(v) # Matrix([[-10/17, 7/17, 2/17]])
u=Matrix([[1,2,3]])
v=Matrix([[-2,3,3]])
print(u.cross(v)) # Matrix([[-3, -9, 7]])
print((u.cross(v).norm()/(u.norm()*v.norm())).n()) # 0.671787690642439
u1,u2,u3=symbols('u1:4')
v1,v2,v3=symbols('v1:4')
print(Matrix([u1,u2,u3]).cross(Matrix([v1,v2,v3])))
print(u.cross(v)) # Matrix([[-3, -9, 7]])
print(v.cross(u)) # Matrix([[3, 9, -7]])

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')
print(E_nl(n,Z)) # -Z**2/(2*n**2)
print(E_nl(1)) # -1/2
print(E_nl(2,4)) # -2
print(E_nl(n,l)) # -l**2/(2*n**2)
print(E_nl_dirac(5,2)) # -0.0200000390505011
print(E_nl_dirac(2,1)) # -0.125000416028342
print(E_nl_dirac(3,2,False)) # -0.0555558020932949
print(R_nl(5,0,r))
print(R_nl(5,0,r,l))
2)矩阵和Pauli代数:
在physics.matrices模块下有一些与物理相关的矩阵。示例:
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')
print(mdft(4))                  # 用离散傅里叶变换做矩阵乘法
print(mgamma(2))                # 狄拉克-伽马矩阵的狄拉克表示
print(msigma(2))                # Pauli矩阵
print(pat_matrix(3,1,0,0))      # 计算平行轴定理矩阵,即沿各轴转动惯量
print(evaluate_pauli_product(4*x*Pauli(3)*Pauli(2)))
3)一维和三维量子谐振子:
模块包含用于计算一维谐振子、三维各向同性谐振子、一维谐振子的波函数和三维各向同性谐振子的径向波函数。示例:
from sympy.physics.qho_1d import E_n,psi_n
from sympy.physics.sho import E_nl,R_nl
from sympy import var,symbols
var('x m omega')
var('r nu l')
x,y,z=symbols('x,y,z')
print(E_n(x,omega))
print(psi_n(2,x,m,omega))
print(E_nl(x,y,z))
print(R_nl(1,0,1,r))
print(R_nl(2,1,1,r))
4)二次量子化:
用于分析和描述一个多体系统的量子的概念称为二次量子化,模块中包括了二次量子化操作的类和玻色子的状态。
from sympy.physics.secondquant import Dagger,B,Bd,BKet
from sympy.physics.secondquant import FockStateBosonKet
from sympy.functions.special.tensor_functions import KroneckerDelta
from sympy.abc import x,y,n,i,j,k
from sympy import Symbol,sqrt
from sympy import I
Dagger(2*I)
Dagger(B(0))
Dagger(Bd(0))
print(KroneckerDelta(1,2))
print(KroneckerDelta(3,3))
print(KroneckerDelta(i,j))
print(KroneckerDelta(i,i))
print(KroneckerDelta(i,i+1))
print(KroneckerDelta(i,i+1+k))
a=Symbol('a',above_fermi=True)
i=Symbol('i',above_fermi=True)
p=Symbol('p')
q=Symbol('q')
print(KroneckerDelta(p,q).indices_contain_equal_information)
print(KroneckerDelta(p,q+1).indices_contain_equal_information)
print(KroneckerDelta(i,p).indices_contain_equal_information)
print(KroneckerDelta(p,a).is_above_fermi)
print(KroneckerDelta(p,i).is_above_fermi)
print(KroneckerDelta(p,q).is_above_fermi)
print(KroneckerDelta(p,a).is_only_above_fermi)
print(KroneckerDelta(p,q).is_only_above_fermi)
print(KroneckerDelta(p,i).is_only_above_fermi)
print(B(x).apply_operator(y))
print(B(0).apply_operator(BKet((n,))))
print(sqrt(n)*FockStateBosonKet((n-1,)))
5)力学:
SymPy的力学模块包括了处理机械系统所需的工具,提供了操作参考系、力、力矩的相关功能。下面的示例计算了作用于任意对象的合力,是向量相加。
from sympy import *
Func1=Matrix([4,0])
Func2=Matrix([5*cos(30*pi/180),5*sin(30*pi/180)])
Func_net=Func1+Func2
print(Func_net) # Matrix([[4 + 5*sqrt(3)/2], [5/2]])
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))
print(v) # Matrix([[4 + 5*sqrt(3)/2], [5/2]])
x=xi+integrate(v,(t,0,t))
print(x) # a + t**2*vi/2 + t*xi
print((v*v).expand()) # t**2*vi**2 + 2*t*vi*xi + xi**2
print(((v*v).expand()-2*a*x).simplify()) # -2*a*vi + xi**2
如果对象上的合力是一次常数,那么这个常力引起的就是恒定加速度的运动。下面代码中还用到均匀加速运动和势能。
from sympy import *
xi=20               # 初始位置
vi=10               # 初始速度
a=5                 # 加速度
t,vi,xi,k=symbols('t vi xi k')
a=sqrt(k*t)
x=xi+integrate(vi+integrate(a,(t,0,t)),(t,0,t))
print(x) # t*vi + xi + 4*(k*t)**(5/2)/(15*k**2)
x.subs({t:3}).n()               # x(3)
diff(x,t).subs({t:3}).n()       # v(3)
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))
print(U_g) # g*h*m
F_s=-k*x                        # 根据胡克定律,弹簧被拉伸位移x后的弹力
U_s=-integrate(F_s,(x,0,x))
print(U_s) # k*x**2/2
下面的代码用dsolve函数来寻找微分方程的位置函数,该微分方程表示了质量弹簧系统的运动。
from sympy import *
t=Symbol('t')                   # 时间
x=Function('x')                 # 位置函数
w=Symbol('w',positive=True)     # 角速度
sol=dsolve(diff(x(t),t,t)+w**3*x(t),x(t))
print(sol) # Eq(x(t), C1*sin(t*w**(3/2)) + C2*cos(t*w**(3/2)))
x=sol.rhs
print(x) # C1*sin(t*w**(3/2)) + C2*cos(t*w**(3/2))
A,phi=symbols('A phi')
(A*cos(w*t-phi)).expand(trig=True)
x=sol.rhs.subs({'C1':0,'C2':A})
print(x) # A*cos(t*w**(3/2))
v=diff(x,t)
k=Symbol('k')
m=Symbol('m')
E_T=(0.3*k*x**3+0.3*m*v**3).simplify()
print(E_T) # 0.3*A**3*(k*cos(t*w**(3/2))**3 - m*w**(9/2)*sin(t*w**(3/2))**3)
E_T.subs({k:m*w**4}).simplify()
E_T.subs({w:sqrt(k/m)}).simplify()

10. 漂亮的打印功能:

SymPy可以用ASCII和Unicode字符漂亮地打印输出。SymPy有一些可用的打印机,包括LaTeX、MathML、Unicode字符打印机、ASCII字符打印机、Str、dot、repr等。
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(pretty_print=False)
theta
init_printing(order='lex')
str(2*y+3*x+2*y**2+x**2+1)
init_printing(order='grlex')
str(2*y+3*x+2*y**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+2*y**2+x**2+1)
init_printing(num_columns=10)
str(2*y+3*x+2*y**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打印用LatexPrinter类来实现,其中有方法可以将给定的表达式转换为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(4*x**2/y))
print(latex(5*x**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'))
上述输出需要配合打印机功能才能实现,普通字符界面看不出效果。

11. 密码学模块:

SymPy包括分组密码和流密码函数,还包括仿射密码affine cipher、二分密码Bifid cipher、ElGamal非对称加密、希尔密码Hill's cipher、RSA非对称密码、教育版RSA非对称密码、线性反馈移位寄存器加密LFSR、移位密码shift cipher、替换式密码substitution cipher、维吉尼亚密码Vigenere's cipher等。对普通文本进行RAS的解密和加密:
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
print(encipher_rsa(pt,publickey)) # 112
privatekey=rsa_private_key(a,b,c)
ct=112
print(decipher_rsa(ct,privatekey)) # 8
对普通文本进行二分加密和解密并返回密文:
from sympy.crypto.crypto import encipher_bifid6,decipher_bifid6
key='encryptingit'
pt='A very good book will be released in 2015'
print(encipher_bifid6(pt,key))
ct='AENUIUKGHECNOIY27XVFPXR52XOXSPI0Q'
print(decipher_bifid6(ct,key))

12. 输入的句法分析:

句法分析模块将输入的字符串解析为SymPy表达式,可以自动判断表达式是否应该包含括号、是否需要做乘法,并且能够让函数实例化。
from sympy.parsing.sympy_parser import *
print(parse_expr('2*x**2+3*x+4')) # 2*x**2 + 3*x + 4
print(parse_expr('10*sin(x)**2+3*x*y*z')) # 3*x*y*z + 10*sin(x)**2
transformations=(standard_transformations+(function_exponentiation,))
print(parse_expr('sin**2(x**2)',transformations=transformations)) #sin(x**2)**2
print(parse_expr('5sin**2 x**2+6abc+sec theta',transformations=(standard_transformations
     +(implicit_multiplication_application,)))) #6*a*b*c + 5*sin(x**2)**2 + sec(theta)

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、逻辑隐含式和相等关系的函数,所有这些函数支持符号形式以及用这些运算符的计算。在符号形式中,表达式由符号表示,不会被运算。
from sympy import symbols
from sympy.logic.boolalg import *
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
Nand(True,False)
Nand(True,True)
Nand(a,b)
Nor(True,False)
Nor(True,True)
Nor(False,True)
Nor(False,False)
Nor(a,b)
Equivalent(False,False,False)
Equivalent(True,False,False)
Equivalent(a,And(a,False))
Implies(False,True)
Implies(True,False)
Implies(False,False)
Implies(True,True)
a>>b
b<<a
逻辑模块还允许使用if-then-else语句,将一个命题逻辑语句转化为合取conjunctive/析取disjunctive的规范形式,也可以检查一个表达式是否为合取/析取的规范形式。if-then-else语句中,当第一个参数值为真时,返回第二个参数,否则返回第三个参数。to_cnf和to_dnf函数将表达式或介词声明分别转换为CNF和DNF,is_cnf和is_dnf分别确认给定的表达式是否是cnf和dnf。
from sympy.abc import A,B,C,X,Y,Z,a,b,c
from sympy.logic.boolalg import *
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函数,将布尔表达式转换为简化的积项之和或和项之积形式,还包括使用简化和多余组消除算法的函数。
from sympy.abc import x,y,z
from sympy.logic import simplify_logic,SOPform,POSform
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)|(~z&~y&~z)'
simplify_logic(expr)
_=S(expr)
simplify_logic(_)

14. 几何模块:

几何模块可以进行二维图形的生成、操作和计算,这些二维图形包括点、线、圆、椭圆、多边形、三角形等。几何模块有一些子模块,可以做各种二维和三维图像的操作,包括点、三维点、线、三维线、曲线、椭圆、多边形、平面等。
其中,collinear函数检验给定的点是否共线,如果共线则返回真;medians函数返回一个以顶点为键、顶点的中间值为值的字典;intersection函数查找两个或多个几何实体相交的点;可以用is_tangent函数判断给定的线是否是圆弧的切线;circumference函数返回圆的周长;equation函数返回圆的方程形式。
from sympy.geometry import *
from sympy 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)

15. 符号积分:

这个模块有一些高级函数,可以计算不同阶数和精度的正交的权重点,该模块还包括一些可以计算定积分和进行积分变换的特殊函数。
正交子模块sympy.integrals.quadrature中的数值积分包括用于执行以下正交计算的函数:
·高斯-勒让得正交Gauss-Legendre quadrature
·高斯-拉盖尔正交Gauss-Laguerre quadrature
·高斯-埃尔米特正交Gauss-Heimite quadrature
·高斯-切比雪夫正交Gauss-Chebyshev quadrature
·高斯-雅可比正交Gauss-Jscobi 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
积分模块可以计算给定表达式的定积分和不定积分,主要有两个函数:
·Integrate(f, x):计算函数f的不定积分
·Integrate(f, (x,m,n)):计算函数f在m到n区间的定积分
from sympy import integrate,log,exp,oo,sqrt,cos,erf
from sympy.abc import n,x,y
from sympy import Symbol
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**n*exp(-x),(x,0,oo))
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)

16. 多项式操作:

SymPy中的多项式模块允许用户做多项式操作,包括多项式简单操作和高级操作。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=4*x**2+8*x+5
g=3*x+1
q,r=div(f,g,domain='QQ')
print(q) # 4*x/3 + 20/9
print(r) # 25/9
(q*g+r).expand()
q,r=div(f,g,domain='ZZ')
print(q) # 0
print(r) # 4*x**2 + 8*x + 5
g=4*x+2
q,r=div(f,g,domain='ZZ')
print(q) x
print(r) # 6*x + 5
g=5*x+1
q,r=div(f,g,domain='ZZ')
print(q) # 0
print(r) # 4*x**2 + 8*x + 5
(q*g+r).expand()
a,b,c=symbols('a b c')
f=a*x**2+b*x+c
g=3*x+2
q,r=div(f,g,domain='QQ')
print(q) # a*x/3 - 2*a/9 + b/3
print(r) # 4*a/9 - 2*b/3 + c
模块还可以实施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=(15*x+15)*x
g=20*x**2
print(gcd(f,g)) # 5*x
f=4*x**2/2
g=16*x/4
print(gcd(f,g)) # 2*x
f=x*y/3+y**2
g=4*x+9*y
print(gcd(f,g)) # 1
f=x*y**2+x**2*y
g=x**2*y**2
print(gcd(f,g)) # x*y
lcm(f,g)
(f*g).expand()
(gcd(f,g,x,y)*lcm(f,g,x,y)).expand()
f=4*x**2+6*x**3+3*x**4+2*x**5
sqf_list(f)
print(sqf(f)) # x**2*(2*x**3 + 3*x**2 + 6*x + 4)
print(factor(x**4/3+6*x**3/16-2*x**2/4)) # x**2*(8*x**2 + 9*x - 12)/24
print(factor(x**2+3*x*y+4*y**2)) # x**2 + 3*x*y + 4*y**2

17. 集合:

SymPy的集合模块可以进行各种集合操作,包括有限集合、区间、单元素集合、全集、自然数集合等等类,还包括并集、交集、乘积集合、补集等集合操作。下面程序中生成一个区间集合和一个有限集合,可以获取区间集合的起始值,可以验证特定元素是否属于一个集合。
from sympy import Symbol,Interval,FiniteSet
Interval(1,10)
Interval(1,10,False,True)
a=Symbol('a',real=True)
Interval(1,a)
Interval(1,100).end
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 Interval,FiniteSet,Intersection,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))

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+3*x**2)/(4*x*sin(y)**2+2*x*cos(y)**2)
print(expr) # (3*x**2 + 2*x)/(4*x*sin(y)**2 + 2*x*cos(y)**2)
simplify(expr)
trigsimp(expr)
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)
print(expr3) # log(x*y**(-log(x) + 1))
count_ops(expr2)
count_ops(expr3)
print(count_ops(expr2,visual=True)) # 2*ADD + DIV + 4*LOG + MUL
print(count_ops(expr3,visual=True)) # 2*LOG + MUL + POW + SUB

 

Copyright@dwenzhao.cn All Rights Reserved   备案号:粤ICP备15026949号
联系邮箱:dwenzhao@163.com  QQ:1608288659