基于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)

x 2 \displaystyle \sqrt{x^{2}} x2

 sympy.sqrt(y ** 2)

y \displaystyle y 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 ⁡ ( 3.14159265358979 n ) \displaystyle \cos{\left(3.14159265358979 n \right)} cos(3.14159265358979n)

sympy.cos(n2 * np.pi)

cos ⁡ ( 3.14159265358979 n ) \displaystyle \cos{\left(3.14159265358979 n \right)} cos(3.14159265358979n)

sympy.cos(n3 *np.pi)

cos ⁡ ( 3.14159265358979 n ) \displaystyle \cos{\left(3.14159265358979 n \right)} cos(3.14159265358979n)

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

sympy.cos(n1 * pi)

cos ⁡ ( π n ) \displaystyle \cos{\left(\pi n \right)} cos(πn)

sympy.cos(n2 * pi)

( − 1 ) n \displaystyle \left(-1\right)^{n} (−1)n

sympy.cos(n3 *pi)

− 1 \displaystyle -1 −1

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

oo

∞ \displaystyle \infty ∞

pi

π \displaystyle \pi π

sympy.E

e \displaystyle e e

sympy.EulerGamma

γ \displaystyle \gamma γ

sympy.I

i \displaystyle i i

sympy.oo

∞ \displaystyle \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 \displaystyle 0.2999999999999999888977698 0.2999999999999999888977698

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

0.3 \displaystyle 0.3 0.3

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

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

浮点错误。

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

r1

2 3 \displaystyle \frac{2}{3} 32

5/6

0.8333333333333334 \displaystyle 0.8333333333333334 0.8333333333333334

x/y

x y \displaystyle \frac{x}{y} yx

定义符号函数

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

f = sympy.Function("f")

type(f)

sympy.core.function.UndefinedFunction

f(x)

f ( x ) \displaystyle f{\left(x \right)} f(x)

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

g ( x , y , z ) \displaystyle g{\left(x,y,z \right)} g(x,y,z)

g.free_symbols

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

g ( x , y , z ) \displaystyle g{\left(x,y,z \right)} g(x,y,z)

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

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

( x ↦ x 2 + 2 x ) \displaystyle \left( x \mapsto x^{2} + 2 x \right) (x↦x2+2x)

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

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

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

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

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

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

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

f(x)

x 2 + 2 x \displaystyle x^{2} + 2 x x2+2x

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

expre

5 x 3 + x 2 + 1 \displaystyle 5 x^{3} + x^{2} + 1 5x3+x2+1

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

expre.args

( 1 , x 2 , 5 x 3 ) \displaystyle \left( 1, \ x^{2}, \ 5 x^{3}\right) (1, x2, 5x3)

expre.args[2]

5 x 3 \displaystyle 5 x^{3} 5x3

expre.args[2].args[1]

x 3 \displaystyle x^{3} x3

expre.args[2].args[0]

5 \displaystyle 5 5

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

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

expr

2 x 2 − x ( x + 1 ) − 2 x \displaystyle 2 x^{2} - x \left(x + 1\right) - 2 x 2x2−x(x+1)−2x

sympy.simplify(expr)

x ( x − 3 ) \displaystyle x \left(x - 3\right) x(x−3)

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

expr2

2 sin ⁡ ( x ) cos ⁡ ( x ) \displaystyle 2 \sin{\left(x \right)} \cos{\left(x \right)} 2sin(x)cos(x)

sympy.simplify(expr2)

sin ⁡ ( 2 x ) \displaystyle \sin{\left(2 x \right)} sin(2x)

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

化简函数

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

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

sympy.powsimp 用幂指数化简

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

sympy.ratsimp 提取公分母

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

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

sympy.expand(expr)

x 2 + 9 x + 20 \displaystyle x^{2} + 9 x + 20 x2+9x+20

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

sin ⁡ ( x + y ) \displaystyle \sin{\left(x + y \right)} sin(x+y)

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

sin ⁡ ( x ) cos ⁡ ( y ) + sin ⁡ ( y ) cos ⁡ ( x ) \displaystyle \sin{\left(x \right)} \cos{\left(y \right)} + \sin{\left(y \right)} \cos{\left(x \right)} sin(x)cos(y)+sin(y)cos(x)

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

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

log ⁡ ( a ) + log ⁡ ( b ) \displaystyle \log{\left(a \right)} + \log{\left(b \right)} log(a)+log(b)

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

i e b sin ⁡ ( a ) + e b cos ⁡ ( a ) \displaystyle i e^{b} \sin{\left(a \right)} + e^{b} \cos{\left(a \right)} iebsin(a)+ebcos(a)

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

a x b x \displaystyle a^{x} b^{x} axbx

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

e a x e − b x \displaystyle e^{a x} e^{- b x} eaxe−bx

因子、汇集和组合

sympy.factor(x**2 - 1)

( x − 1 ) ( x + 1 ) \displaystyle \left(x - 1\right) \left(x + 1\right) (x−1)(x+1)

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

x ( sin ⁡ ( z ) + cos ⁡ ( y ) ) \displaystyle x \left(\sin{\left(z \right)} + \cos{\left(y \right)}\right) x(sin(z)+cos(y))

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

log ⁡ ( a b ) \displaystyle \log{\left(\frac{a}{b} \right)} log(ba)

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

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

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

x 2 y z + x + y \displaystyle x^{2} y z + x + y x2yz+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))

( sin ⁡ ( x ) + cos ⁡ ( x ) ) ( − sin ⁡ ( y ) + cos ⁡ ( y ) ) \displaystyle \left(\sin{\left(x \right)} + \cos{\left(x \right)}\right) \left(- \sin{\left(y \right)} + \cos{\left(y \right)}\right) (sin(x)+cos(x))(−sin(y)+cos(y))

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

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

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

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

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

− 1 x + 2 + 1 x + 1 \displaystyle - \frac{1}{x + 2} + \frac{1}{x + 1} −x+21+x+11

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

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

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

1 x + 1 \displaystyle \frac{1}{x + 1} x+11

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

y x y + 4 \displaystyle \frac{y}{x y + 4} xy+4y

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

1 x + 4 y \displaystyle \frac{1}{x + 4 y} x+4y1

代换 Substitutions

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

2 y \displaystyle 2 y 2y

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

sin ⁡ ( y e y ) \displaystyle \sin{\left(y e^{y} \right)} sin(yey)

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

cos ⁡ ( y e y ) \displaystyle \cos{\left(y e^{y} \right)} cos(yey)

当然也可以带入值

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

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

expr.subs(values)

13.3 \displaystyle 13.3 13.3

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

sympy.N(pi, 50)

3.1415926535897932384626433832795028841971693993751 \displaystyle 3.1415926535897932384626433832795028841971693993751 3.1415926535897932384626433832795028841971693993751

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

x + 0.3183098862 \displaystyle x + 0.3183098862 x+0.3183098862

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

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

d d x f ( x ) \displaystyle \frac{d}{d x} f{\left(x \right)} dxdf(x)

sympy.diff(f, x, x)

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

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

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

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

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

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

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

∂ 5 ∂ y 2 ∂ x 3 g ( x , y ) \displaystyle \frac{\partial^{5}}{\partial y^{2}\partial x^{3}} g{\left(x,y \right)} ∂y2∂x3∂5g(x,y)

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

expr.diff(x)

4 x 3 + 3 x 2 + 2 x + 1 \displaystyle 4 x^{3} + 3 x^{2} + 2 x + 1 4x3+3x2+2x+1

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

expr.diff(x)

y cos ⁡ ( x 2 ) cos ⁡ ( x y ) − sin ⁡ ( x 2 ) sin ⁡ ( x y ) 2 \displaystyle y \cos{\left(\frac{x}{2} \right)} \cos{\left(x y \right)} - \frac{\sin{\left(\frac{x}{2} \right)} \sin{\left(x y \right)}}{2} ycos(2x)cos(xy)−2sin(2x)sin(xy)

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

expr.diff(x).doit()

2 x π polygamma ⁡ ( 0 , 1 2 − x 2 ) 2 Γ ( 1 2 − x 2 ) + 2 x π log ⁡ ( 2 ) Γ ( 1 2 − x 2 ) \displaystyle \frac{2^{x} \sqrt{\pi} \operatorname{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)} 2Γ(21−2x)2xπ polygamma(0,21−2x)+Γ(21−2x)2xπ log(2)

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

d d x e cos ⁡ ( x ) \displaystyle \frac{d}{d x} e^{\cos{\left(x \right)}} dxdecos(x)

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

− e cos ⁡ ( x ) sin ⁡ ( x ) \displaystyle - e^{\cos{\left(x \right)}} \sin{\left(x \right)} −ecos(x)sin(x)

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

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

 sympy.integrate(f)

∫ f ( x ) d x \displaystyle \int f{\left(x \right)}\, dx ∫f(x)dx

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

∫ a b f ( x ) d x \displaystyle \int\limits_{a}^{b} f{\left(x \right)}\, dx a∫bf(x)dx

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

π 2 \displaystyle \frac{\sqrt{\pi}}{2} 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 − x 2 2 + x 4 24 + O ( x 6 ) \displaystyle 1 - \frac{x^{2}}{2} + \frac{x^{4}}{24} + O\left(x^{6}\right) 1−2x2+24x4+O(x6)

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

1 − x + x 2 − x 3 + x 4 − x 5 + O ( x 6 ) \displaystyle 1 - x + x^{2} - x^{3} + x^{4} - x^{5} + O\left(x^{6}\right) 1−x+x2−x3+x4−x5+O(x6)

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

 expr.series(x, n=4)

1 − x y + x 2 ( y 2 − 1 2 ) + x 3 ( − 5 y 3 6 + y 2 ) + O ( x 4 ) \displaystyle 1 - x y + x^{2} \left(y^{2} - \frac{1}{2}\right) + x^{3} \left(- \frac{5 y^{3}}{6} + \frac{y}{2}\right) + O\left(x^{4}\right) 1−xy+x2(y2−21)+x3(−65y3+2y)+O(x4)

expr.series(y, n=4)

cos ⁡ ( x ) − x y cos ⁡ ( x ) + x 2 y 2 cos ⁡ ( x ) − 5 x 3 y 3 cos ⁡ ( x ) 6 + O ( y 4 ) \displaystyle \cos{\left(x \right)} - x y \cos{\left(x \right)} + x^{2} y^{2} \cos{\left(x \right)} - \frac{5 x^{3} y^{3} \cos{\left(x \right)}}{6} + O\left(y^{4}\right) cos(x)−xycos(x)+x2y2cos(x)−65x3y3cos(x)+O(y4)

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

1 \displaystyle 1 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 ⁡ ( x ) \displaystyle - \sin{\left(x \right)} −sin(x)

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

cos ⁡ ( x ) \displaystyle \cos{\left(x \right)} cos(x)

以上代码用于说名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

( 1 2 , − 1 ) \displaystyle \left( \frac{1}{2}, \ -1\right) (21, −1)

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

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

n \displaystyle n n

∑ n = 1 ∞ 1 n 2 \displaystyle \sum_{n=1}^{\infty} \frac{1}{n^{2}} n=1∑∞n21

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

π 2 6 \displaystyle \frac{\pi^{2}}{6} 6π2

x = sympy.Symbol("x")

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

e x − 1 \displaystyle e^{x} - 1 ex−1

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

wz.doit()

x ( e x x − 1 x ) \displaystyle x \left(\frac{e^{x}}{x} - \frac{1}{x}\right) x(xex−x1)

wz.doit().simplify()

e x − 1 \displaystyle e^{x} - 1 ex−1

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

方程

 x = sympy.Symbol("x")

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

[ − 3 , 1 ] \displaystyle \left[ -3, \ 1\right] [−3, 1]

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

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

[ − b + − 4 a c + b 2 2 a , − b + − 4 a c + b 2 2 a ] \displaystyle \left[ \frac{- b + \sqrt{- 4 a c + b^{2}}}{2 a}, \ - \frac{b + \sqrt{- 4 a c + b^{2}}}{2 a}\right] [2a−b+−4ac+b2 , −2ab+−4ac+b2 ]

eq1 = x + 2 * y - 1

eq2 = x - y + 1

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

[ { x : − 1 3 , y : 2 3 } ] \displaystyle \left[ \left\{ x : - \frac{1}{3}, \ y : \frac{2}{3}\right\}\right] [{x:−31, y:32}]

#也可以获取结果
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])

[ 1 2 ] \displaystyle \left[\begin{matrix}1\\2\end{matrix}\right] [12]

np.array([1,2])

array([1, 2])

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

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

[ a b c d ] \displaystyle \left[\begin{matrix}a & b\\c & d\end{matrix}\right] [acbd]

 M * M

[ a 2 + b c a b + b d a c + c d b c + d 2 ] \displaystyle \left[\begin{matrix}a^{2} + b c & a b + b d\\a c + c d & b c + d^{2}\end{matrix}\right] [a2+bcac+cdab+bdbc+d2]

M.norm()

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

M.rank()

2 \displaystyle 2 2

M.trace()

a + d \displaystyle a + d a+d

M.transpose()

[ a c b d ] \displaystyle \left[\begin{matrix}a & c\\b & d\end{matrix}\right] [abcd]

M^T

---------------------------------------------------------------------------NameError                                 Traceback (most recent call last)<ipython-input-193-0be5c4114203> in <module>
----> 1 M^TNameError: name 'T' is not defined

M.T

[ a c b d ] \displaystyle \left[\begin{matrix}a & c\\b & d\end{matrix}\right] [abcd]

M.H

[ a ‾ c ‾ b ‾ d ‾ ] \displaystyle \left[\begin{matrix}\overline{a} & \overline{c}\\\overline{b} & \overline{d}\end{matrix}\right] [abcd]

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)

[ b 1 − p ( − b 1 q + b 2 ) − p q + 1 − b 1 q + b 2 − p q + 1 ] \displaystyle \left[\begin{matrix}b_{1} - \frac{p \left(- b_{1} q + b_{2}\right)}{- p q + 1}\\\frac{- b_{1} q + b_{2}}{- p q + 1}\end{matrix}\right] [b1−−pq+1p(−b1q+b2)−pq+1−b1q+b2]

x1 = M.inv() * b

x1

[ b 1 − p q + 1 − b 2 p − p q + 1 − b 1 q − p q + 1 + b 2 − p q + 1 ] \displaystyle \left[\begin{matrix}\frac{b_{1}}{- p q + 1} - \frac{b_{2} p}{- p q + 1}\\- \frac{b_{1} q}{- p q + 1} + \frac{b_{2}}{- p q + 1}\end{matrix}\right] [−pq+1b1−−pq+1b2p−−pq+1b1q+−pq+1b2]

x ==x1

False

rt = x.T*x1
rt

[ ( − b 1 + b 2 p ) ( b 1 − p ( − b 1 q + b 2 ) − p q + 1 ) p q − 1 + ( − b 1 q + b 2 ) ( b 1 q − b 2 ) ( − p q + 1 ) ( p q − 1 ) ] \displaystyle \left[\begin{matrix}\frac{\left(- b_{1} + b_{2} p\right) \left(b_{1} - \frac{p \left(- b_{1} q + b_{2}\right)}{- p q + 1}\right)}{p q - 1} + \frac{\left(- b_{1} q + b_{2}\right) \left(b_{1} q - b_{2}\right)}{\left(- p q + 1\right) \left(p q - 1\right)}\end{matrix}\right] [pq−1(−b1+b2p)(b1−−pq+1p(−b1q+b2))+(−pq+1)(pq−1)(−b1q+b2)(b1q−b2)]

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)

[ 25 76 29 76 ] \displaystyle \left[\begin{matrix}\frac{25}{76}\\\frac{29}{76}\end{matrix}\right] [76257629]

[ b 1 b 2 ] \displaystyle \left[\begin{matrix}b_{1}\\b_{2}\end{matrix}\right] [b1b2]

b[0]

b 1 \displaystyle b_{1} b1

b[1]

b 2 \displaystyle b_{2} b2

x1.subs(ics)

[ 25 76 29 76 ] \displaystyle \left[\begin{matrix}\frac{25}{76}\\\frac{29}{76}\end{matrix}\right] [76257629]

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

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

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