基于Python3的科学运算与常用算法-第4 符号运算

本系列是使用Ubuntu 64位操作系统上，anaconda环境下，采用jupyter notebook 实现python3科学计算的文章，应用于青岛科技大学信息学院科学计算与数据挖掘等多个课程中python教学。需要安装的模块，numpy ，matplotlib，scipy，sympy，pandas等科学计算常用模块。由于匆忙成稿，个人能力所限，显得粗糙，并且有几个地方是错误的。主要参考了Numerical Python 2nd Edition

科学计算包括数值计算和符号计算两种计算 。计算机能够对数值进行一系列运算是人所共知的事，但计算机也能够对含未知量的式子直接进行推导、演算则并不是人人皆知。数值计算和符号计算本来应该是并存的两种计算，是计算的平行的两个部分，决不能厚此薄彼，因此这两种计算都是一样重要的。利用计算机对一个函数进行求导、积分，这早己成为事实。

import sympy

sympy.init_printing()#调用及时显示系统。但对于jupyter似乎不需要。

 x = sympy.Symbol("x")#创建符号

符号名是一个字符串，可以选择包含类似Latex标记

使符号名在ipython的显示系统中显示良好。

符号对象的名称在创建时设置。符号可以采用不同方式创建，

例如，使用sympy.symbol、sympy.symbols和sympy。

通常需要将sympy符号与python变量关联起来

同一个名称或与符号名称紧密对应的名称。例如,

要创建一个名为x的符号，并将其绑定到同名的python变量，

我们可以使用symbol类的构造函数并传递包含该符号的字符串

名称作为第一个参数

 y = sympy.Symbol("y", real=True)

y.is_real

True

 x.is_real is None

True

另一方面，如果我们要使用is real来查询先前定义的符号

x，它没有显式地指定为real，因此可以同时表示real和

非实变量，结果是没有.但如果被显式定义，则为false

y.is_real is None

False

sympy.Symbol("z", imaginary=True).is_real

False

 x = sympy.Symbol("x")

y = sympy.Symbol("y", positive=True)

 sympy.sqrt(x ** 2)

$\sqrt{x^2}$

 sympy.sqrt(y ** 2)

$y$

#以上运算可以清晰的发现定义的清晰度，可以影响运算结果。

n1 = sympy.Symbol("n")

n2 = sympy.Symbol("n", integer=True)

n3 = sympy.Symbol("n", odd=True)

import numpy as np

sympy.cos(n1 * np.pi)

$\cos \left(3.14159265358979n\right)$

sympy.cos(n2 * np.pi)

$\cos \left(3.14159265358979n\right)$

sympy.cos(n3 *np.pi)

$\cos \left(3.14159265358979n\right)$

 from sympy import I, pi, oo #然而在sympy中不建议采用sympy之外的常数

sympy.cos(n1 * pi)

$\cos \left(\pi n\right)$

sympy.cos(n2 * pi)

${\left(-1\right)}^n$

sympy.cos(n3 *pi)

$-1$

可以看到采用sympy内部的常数，虚数，π，无穷大，更利于公式简化推导。

oo

$\infty$

pi

$\pi$

sympy.E

$e$

sympy.EulerGamma

$\gamma$

sympy.I

$i$

sympy.oo

$\infty$

创建符号时，符号的属性可以用下列属性标定：

real ,imaginary 属性判定 is_real, is_imaginary

positive, negative 属性判定 is_positive, is_negative

integer 属性判定 is_integer

odd, even is_odd, 属性判定 is_even

prime 属性判定 is_prime

finite, infinite 属性判定 is_finite, is_infinite

a, b, c = sympy.symbols("a, b, c", negative=True)

d, e, f = sympy.symbols("d, e, f", positive=True)

i = sympy.Integer(19)

type(i)

sympy.core.numbers.Integer

i.is_Integer, i.is_real, i.is_odd

(True, True, True)

 "%.25f" % 0.3  # create a string representation with 25 decimals，把0.3表示成有25个小数的浮点数。

'0.2999999999999999888977698'

sympy.Float(0.3, 25)

$0.2999999999999999888977698$

sympy.Float('0.3', 25)#请同学们注意的是，这个才是精确的表示浮点数0.3

$0.3$

但是，请注意，要正确地将0.3表示为浮点对象，必须

从字符串“0.3”而不是python float 0.3初始化它，后者已经包含

浮点错误。

r1 = sympy.Rational(2, 3)#分式的表示方式

r1

$\frac{2}{3}$

5/6

$0.8333333333333334$

x/y

$\frac{x}{y}$

定义符号函数

x,y,z = sympy.symbols("x,y,z")

f = sympy.Function("f")

type(f)

sympy.core.function.UndefinedFunction

f(x)

$f\left(x\right)$

g = sympy.Function("g")(x,y,z)

g

$g\left(x,y,z\right)$

g.free_symbols

{ x , y , z } \displaystyle \left\{x, y, z\right\} {x,y,z}

g

$g\left(x,y,z\right)$

定义函数符号，下面可以定义函数的表达式

h = sympy.Lambda(x,x**2+2*x)

h

$\left(x\mapsto x^2+2x\right)$

#lambda函数定义一类函数只有当你的情况完全满足这四个标准时, 我才会说你可以使用 lambda 表达式:
所谓的匿名函数

你所要做的操作是不重要的:函数不值得一个名称

使用 lambda 表达式比你所能想到的函数名称让代码更容易理解

你很确定还没有一个函数能满足你的需求

你团队的每个人都了解 lambda 表达式, 并且都同意使用它们

如果上述四条中的任何一条都不符合你的情况, 我建议用 def 来写一个新的函数, (如果可能)接受一个在 Python 中已经存在且能满足你需求的函数.

def f(x):
    return x**2+2*x

f(x)

$x^2+2x$

#很多时候，我们并不需要函数表示，而是需要如下的表达式的(柄,感觉叫句柄有点意思)
expre = 1+x**2+5*x**3

expre

$5x^3+x^2+1$

#expre,有自己的完善的参数体系，这个体系呈现树形结构。

expre.args

$\left(1,x^2,5x^3\right)$

expre.args[2]

$5x^3$

expre.args[2].args[1]

$x^3$

expre.args[2].args[0]

$5$

#化简操作，是所有表达式运算过程中最常用的函数

expr = 2*(x**2-x) - x*(x+1)

expr

$2x^2-x\left(x+1\right)-2x$

sympy.simplify(expr)

$x\left(x-3\right)$

expr2 =2*sympy.cos(x)*sympy.sin(x)

expr2

$2\sin \left(x\right)\cos \left(x\right)$

sympy.simplify(expr2)

$\sin \left(2x\right)$

化简成功的关键是对所有的函数和常数采用sympy定义的函数或常数。

化简函数

sympy.simplify 调用不同的方法得到化简结果

sympy.trigsimp 尝试三角恒等式化简结果

sympy.powsimp 用幂指数化简

sympy.compsimp 通过组合化简多项式

sympy.ratsimp 提取公分母

下面是展开函数的使用说明：

expr=(x+4)*(x+5)

sympy.expand(expr)

$x^2+9x+20$

sympy.sin(x+y).expand()

$\sin \left(x+y\right)$

sympy.sin(x+y).expand(trig= True)

$\sin \left(x\right)\cos \left(y\right)+\sin \left(y\right)\cos \left(x\right)$

a, b = sympy.symbols("a, b", positive=True)# 使用化简必须保证，所有的数学符号，函数都是sympy产生的

sympy.log(a * b).expand(log=True)

$\log \left(a\right)+\log \left(b\right)$

sympy.exp(I*a + b).expand(complex=True)

$i{e}^b\sin \left(a\right)+e^b\cos \left(a\right)$

sympy.expand((a * b)**x, power_base=True)

$a^x{b}^x$

sympy.exp((a-b)*x).expand(power_exp=True)

$e^{ax}{e}^{-bx}$

因子、汇集和组合

sympy.factor(x**2 - 1)

$\left(x-1\right)\left(x+1\right)$

 sympy.factor(x * sympy.cos(y) + sympy.sin(z) * x)

$x\left(\sin \left(z\right)+\cos \left(y\right)\right)$

sympy.logcombine(sympy.log(a) - sympy.log(b))

$\log \left(\frac{a}{b}\right)$

sympy.trigsimp, sympy.powsimp, and sympy.logcombine,都可以产生简单的化简

expr = x + y + x**2 * y * z

expr.collect(x) #按照x的指数大小排列

$x^{2}yz+x+y$

expr = sympy.cos(x + y) + sympy.sin(x - y)

 expr.expand(trig=True).collect([sympy.cos(x), 
                                 sympy.sin(x)]).collect(sympy.cos(y) - sympy.sin(y))

$\left(\sin \left(x\right)+\cos \left(x\right)\right)\left(-\sin \left(y\right)+\cos \left(y\right)\right)$

expr.expand(trig=True).
collect([sympy.cos(x), sympy.sin(x)]).collect(sympy.cos(y) - sympy.sin(y))

  File "<ipython-input-96-fefb96c05190>", line 1
    expr.expand(trig=True).
                           ^
SyntaxError: invalid syntax

expr.expand(trig=True).collect([
    sympy.cos(x), sympy.sin(x)]).collect(sympy.cos(y) - sympy.sin(y))

$\left(\sin \left(x\right)+\cos \left(x\right)\right)\left(-\sin \left(y\right)+\cos \left(y\right)\right)$

事实上，在那个地方回车是有讲究的，如果愿意全部试试。

 sympy.apart(1/(x**2 + 3*x + 2), x)#非常有用的分解

$-\frac{1}{x+2}+\frac{1}{x+1}$

 sympy.together(1 / (y * x + y) + 1 / (1+x))

$\frac{y+1}{y\left(x+1\right)}$

sympy.cancel(y / (y * x + y)) #约分化简

$\frac{1}{x+1}$

sympy.cancel(y / (y * x + 4)) 

$\frac{y}{xy+4}$

sympy.cancel(y / (y * x + 4*y**2)) 

$\frac{1}{x+4y}$

代换 Substitutions

(x + y).subs(x, y)

$2y$

sympy.sin(x * sympy.exp(x)).subs(x, y)

$\sin \left(y{e}^y\right)$

 sympy.sin(x * z).subs({z: sympy.exp(y), x: y, sympy.sin: sympy.cos})

$\cos \left(y{e}^y\right)$

当然也可以带入值

expr = x * y + z**2 *x

values = {x: 1.25, y: 0.4, z: 3.2}

expr.subs(values)

$13.3$

赋值过程中，追求精度，必须要采用，sympy.N()函数，或者用evalf()函数、

sympy.N(pi, 50)

$3.1415926535897932384626433832795028841971693993751$

 (x + 1/pi).evalf(10)

$x+0.3183098862$

f = sympy.Function('f')(x)

sympy.diff(f, x)               # equivalent to f.diff(x)

$\frac{d}{dx}f\left(x\right)$

sympy.diff(f, x, x)

$\frac{d^2}{d{x}^2}f\left(x\right)$

sympy.diff(f, x, 3)   # equivalent to sympy.diff(f, x, x, x)

$\frac{d^3}{d{x}^3}f\left(x\right)$

g = sympy.Function('g')(x, y)

 g.diff(x, y)          # equivalent to sympy.diff(g, x, y)

$\frac{{\partial }^2}{\partial y\partial x}g\left(x,y\right)$

 g.diff(x, 3, y, 2)    # equivalent to sympy.diff(g, x, x, x, y, y)

$\frac{{\partial }^5}{\partial y^2\partial x^3}g\left(x,y\right)$

expr = x**4 + x**3 + x**2 + x + 1

expr.diff(x)

$4x^3+3x^2+2x+1$

expr = sympy.sin(x * y) * sympy.cos(x / 2)

expr.diff(x)

$y\cos \left(\frac{x}{2}\right)\cos \left(xy\right)-\frac{\sin \left(\frac{x}{2}\right)\sin \left(xy\right)}{2}$

 expr2 = sympy.special.polynomials.hermite(x, 0)

expr.diff(x).doit()

$\frac{2^x\sqrt{\pi }\mathrm{polygamma}\left(0,\frac{1}{2}-\frac{x}{2}\right)}{2\Gamma \left(\frac{1}{2}-\frac{x}{2}\right)}+\frac{2^x\sqrt{\pi }\log \left(2\right)}{\Gamma \left(\frac{1}{2}-\frac{x}{2}\right)}$

 d = sympy.Derivative(sympy.exp(sympy.cos(x)), x)

d

$\frac{d}{dx}{e}^{\cos \left(x\right)}$

d.doit()#执行上一步运算

$-e^{\cos \left(x\right)}\sin \left(x\right)$

 a, b, x, y = sympy.symbols("a, b, x, y")

 f = sympy.Function("f")(x)

 sympy.integrate(f)

$\int f\left(x\right)dx$

sympy.integrate(f, (x, a, b))

$\underset{a}{\overset{b}{\int }}f\left(x\right)dx$

sympy.integrate(sympy.exp(-x**2), (x, 0, oo))

$\frac{\sqrt{\pi }}{2}$

级数

 x, y = sympy.symbols("x, y")

 f = sympy.Function("f")(x)

sympy.series(f, x)

KaTeX parse error: Undefined control sequence: \substack at position 87: …ght)} \right|_{\̲s̲u̲b̲s̲t̲a̲c̲k̲{ x=0 }} + \fra…

 x0 = sympy.Symbol("{x_0}")

f.series(x, x0, n=3)

KaTeX parse error: Undefined control sequence: \substack at position 124: …ght)} \right|_{\̲s̲u̲b̲s̲t̲a̲c̲k̲{ \xi_{1}={x_0}…

f.series(x, x0, n=2).removeO() #移除掉截断项，参与后续运算。

KaTeX parse error: Undefined control sequence: \substack at position 99: …ght)} \right|_{\̲s̲u̲b̲s̲t̲a̲c̲k̲{ \xi_{1}={x_0}…

f.series(x, x0, n=3).removeO()

KaTeX parse error: Undefined control sequence: \substack at position 117: …ght)} \right|_{\̲s̲u̲b̲s̲t̲a̲c̲k̲{ \xi_{1}={x_0}…

sympy.cos(x).series()

$1-\frac{x^2}{2}+\frac{x^4}{24}+O\left(x^6\right)$

(1/(1+x)).series()

$1-x+x^2-x^3+x^4-x^5+O\left(x^6\right)$

 expr = sympy.cos(x) / (1 + sympy.sin(x * y))

 expr.series(x, n=4)

$1-xy+x^2\left(y^2-\frac{1}{2}\right)+x^3\left(-\frac{5y^3}{6}+\frac{y}{2}\right)+O\left(x^4\right)$

expr.series(y, n=4)

$\cos \left(x\right)-xy\cos \left(x\right)+x^2{y}^2\cos \left(x\right)-\frac{5x^3{y}^3\cos \left(x\right)}{6}+O\left(y^4\right)$

sympy.limit(sympy.sin(x) / x, x, 0) #极限

$1$

f = sympy.Function('f')

x, h = sympy.symbols("x, h")

diff_limit = (f(x + h) - f(x))/h

sympy.limit(diff_limit.subs(f, sympy.cos), h, 0)

$-\sin \left(x\right)$

sympy.limit(diff_limit.subs(f, sympy.sin), h, 0)

$\cos \left(x\right)$

以上代码用于说名sin函数的极限由来，没有太强的意义

 expr = (x**2 - 3*x) / (2*x - 2)

p = sympy.limit(expr/x, x, sympy.oo)

q = sympy.limit(expr - p*x, x, sympy.oo)

p,q

$\left(\frac{1}{2},-1\right)$

 n = sympy.symbols("n", integer=True)

 x = sympy.Sum(1/(n**2), (n, 1, oo))

n

$n$

x

$\sum _{n=1}^{\infty }\frac{1}{n^2}$

x.doit() #有的函数必须用doit（）函数来完成，最终的运算

$\frac{{\pi }^2}{6}$

x = sympy.Symbol("x")

sympy.Sum((x)**n/(sympy.factorial(n)), (n, 1, oo)).doit().simplify()

$e^x-1$

wz = sympy.Sum((x)**n/(sympy.factorial(n)), (n, 1, oo))

wz.doit()

$x\left(\frac{e^x}{x}-\frac{1}{x}\right)$

wz.doit().simplify()

$e^x-1$

涉及到各类化简，再加上doit（）命令，在函数化简的路上，要多尝试！

方程

 x = sympy.Symbol("x")

 sympy.solve(x**2 + 2*x - 3)

$\left[-3,1\right]$

 a, b, c = sympy.symbols("a, b, c")

sympy.solve(a * x**2 + b * x + c, x)

$\left[\frac{-b+\sqrt{-4ac+b^2}}{2a},-\frac{b+\sqrt{-4ac+b^2}}{2a}\right]$

eq1 = x + 2 * y - 1

eq2 = x - y + 1

 sympy.solve([eq1, eq2], [x, y], dict=True)

$\left[\left\{x:-\frac{1}{3},y:\frac{2}{3}\right\}\right]$

#也可以获取结果
sols = sympy.solve([eq1, eq2], [x, y], dict=True)

[eq1.subs(sol).simplify() == 0 and eq2.subs(sol).simplify() == 0 for sol in sols]

[True]

线性代数

sympy.Matrix([1,2])

$\left[\begin{array}{c}1\\ 2\end{array}\right]$

np.array([1,2])

array([1, 2])

 a, b, c, d = sympy.symbols("a, b, c, d")

 M = sympy.Matrix([[a, b], [c, d]])

M

$\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]$

 M * M

$\left[\begin{array}{cc}{a}^2+bc& ab+bd\\ ac+cd& bc+d^2\end{array}\right]$

M.norm()

$\sqrt{{\left|a\right|}^2+{\left|b\right|}^2+{\left|c\right|}^2+{\left|d\right|}^2}$

M.rank()

$2$

M.trace()

$a+d$

M.transpose()

$\left[\begin{array}{cc}a& c\\ b& d\end{array}\right]$

M^T

---------------------------------------------------------------------------

NameError                                 Traceback (most recent call last)

<ipython-input-193-0be5c4114203> in <module>
----> 1 M^T


NameError: name 'T' is not defined

M.T

$\left[\begin{array}{cc}a& c\\ b& d\end{array}\right]$

M.H

$\left[\begin{array}{cc}\overline{a}& \overline{c}\\ \overline{b}& \overline{d}\end{array}\right]$

transpose/T 转置
adjoint/H 转置加共轭
trace 求迹
det 行列式
inv 矩阵的逆
LUdecomposition LU分解
LUsolve LU求解方程
QRdecomposition 分解
diagonalize 对角化 D = P^(-1) MP
norm 范数
nullspace Compute a set of vectors that span the null space of a Matrix.
rank 计算矩阵的秩
singular_values 奇异值
solve Solve a linear system of equations in the form Mx = b.

example： $ x+py = b_{1} $,$ qx+y = b_{2} $
我们可以把上面的方程用矩阵表示为
$ MX =b $

p, q = sympy.symbols("p, q")

M = sympy.Matrix([[1, p], [q, 1]])

b = sympy.Matrix(sympy.symbols("b_1, b_2"))

x = M.LUsolve(b)

x

$\left[\begin{array}{c}{b}_1-\frac{p\left(-b_{1}q+b_2\right)}{-pq+1}\\ \frac{-b_{1}q+b_2}{-pq+1}\end{array}\right]$

x1 = M.inv() * b

x1

$\left[\begin{array}{c}\frac{b_1}{-pq+1}-\frac{b_{2}p}{-pq+1}\\ -\frac{b_{1}q}{-pq+1}+\frac{b_2}{-pq+1}\end{array}\right]$

x ==x1

False

rt = x.T*x1
rt

$\left[\begin{array}{c}\frac{\left(-b_1+b_{2}p\right)\left(b_1-\frac{p\left(-b_{1}q+b_2\right)}{-pq+1}\right)}{pq-1}+\frac{\left(-b_{1}q+b_2\right)\left(b_{1}q-b_2\right)}{\left(-pq+1\right)\left(pq-1\right)}\end{array}\right]$

rt.collect(p)#调用化简函数都会报错，原因在于多变量矩阵没有该类化简的属性

---------------------------------------------------------------------------

AttributeError                            Traceback (most recent call last)

<ipython-input-232-b3d07a6ec5ae> in <module>
----> 1 rt.collect(p)


~/.conda/envs/qutip-env/lib/python3.7/site-packages/sympy/matrices/matrices.py in __getattr__(self, attr)
   2139         else:
   2140             raise AttributeError(
-> 2141                 "%s has no attribute %s." % (self.__class__.__name__, attr))
   2142 
   2143     def __len__(self):


AttributeError: MutableDenseMatrix has no attribute collect.

ics={b[0]:3,b[1]:4,p:7,q:11}
x.subs(ics)

$\left[\begin{array}{c}\frac{25}{76}\\ \frac{29}{76}\end{array}\right]$

b

$\left[\begin{array}{c}{b}_1\\ b_2\end{array}\right]$

b[0]

$b_1$

b[1]

$b_2$

x1.subs(ics)

$\left[\begin{array}{c}\frac{25}{76}\\ \frac{29}{76}\end{array}\right]$

由此可以看到，直接比较是会得到荒谬的答案的，然而带入数据后，发现其实是相等的。

然而，如果解方程mx=b是目标，计算矩阵的逆比执行LU分解的难度更大，

那么使用lu分解更有效。这对于较大的方程系统，尤其重要