python 符号计算_使用Python做科学计算之符号计算

这里我们使用的是Python3 、Sympy库和Jupyter NoteBook。

预备知识：熟悉Python3语法，了解numpy。

>>> import sympy

>>> sympy.I #复数

I

>>> sympy.I ** 2

-1

>>> sympy.sqrt(-1)

I

>>> sympy.E #自然底数

E

>>> sympy.log(sympy.E)

1

>>> sympy.oo #无穷

oo

>>> 1/sympy.oo

0

>>> 1+sympy.oo

oo

>>> sympy.pi #圆周率

pi

>>> sympy.sin(sympy.pi/2)

1

>>>

二、Sympy简单计算。

Sympy有三种内置数据类型，浮点数，有理数，整数。

还有一些初等函数，示例如下：

>>> sympy.Rational(5,6)

5/6

>>> sympy.sin(12.5)

-0.0663218973512007

>>> sympy.cos(56.9)

0.938915058767187

>>> sympy.log(10,100)

log(10)/log(100)

>>> _

log(10)/log(100)

>>> a = _

>>> type(a)

>>> sympy.log(1000,10)

3

>>> sympy.log(4564,10)

log(4564)/log(10)

>>> sympy.sqrt(456)

2*sqrt(114)

>>> sympy.sqrt(4)

2

>>> sympy.factorial(9)

362880

>>>

因为在命令行下出了点问题，导致一些符号的显式并不能完美适配。

(1)函数方程与表达式求值。

>>> x = sympy.Symbol('x') #获得变量

>>> f = x**2 + x + 1 #构建表达式

>>> f

2

x + x + 1

>>> f.evalf(subs = {x:67}) #带入值，格式为 subs = {符号:数值}

4557.00000000000

>>> sympy.solve(f) #解方程

⎡ 1 √3⋅ⅈ 1 √3⋅ⅈ⎤

⎢- ─ - ────, - ─ + ────⎥

⎣ 2 2 2 2 ⎦

>>>

>>> # --------------------------------------------------

>>> df = sympy.symbols('x y') #多元表达式

>>> df

(x, y)

>>> x,y = sympy.symbols('x y')

>>> df = x * y

>>> df.evalf(subs = {x : 45 , y : 56})

2520.00000000000

>>> df.evalf(subs = {x : 1}) #带单个值不能得出值来。

x⋅y

>>> sympy.solve(df,x) #解方程，后边的参数表示对哪一个参数求解

[0]

>>> sympy.solve(df,y)

[0]

>>>

我们也可以使用Lambda函数来构建表达式。这时表达式更像一个函数。

>>> x = sympy.symbols('x') #获取变量

>>> x

x

>>> f = sympy.Lambda(x , 1 / (x**2+x+1)) #构造方程

>>> f

1

x ↦ ──────────

2

x + x + 1

>>> sympy.integrate(f(x) , x) #求积分

⎛2⋅√3⋅x √3⎞

2⋅√3⋅atan⎜────── + ──⎟

⎝ 3 3 ⎠

──────────────────────

3

>>> f(x) #这里就像数学中的函数一样使用了

1

──────────

2

x + x + 1

>>> f(2)

1/7

>>> pf = sympy.Lambda(x , sympy.integrate(f(x) , x))

>>> pf

⎛2⋅√3⋅x √3⎞

2⋅√3⋅atan⎜────── + ──⎟

⎝ 3 3 ⎠

x ↦ ──────────────────────

3

>>> sympy.diff(pf(x) , x) #求导数

4

──────────────────────

⎛ 2 ⎞

⎜⎛2⋅√3⋅x √3⎞ ⎟

3⋅⎜⎜────── + ──⎟ + 1⎟

⎝⎝ 3 3 ⎠ ⎠

>>> sympy.simplify(sympy.diff(pf(x) , x)) #表达式化简，因为符号计算没有经过化简，化简到最简形式。

1

──────────

2

x + x + 1

>>> pf(4) - pf(0)

√3⋅π 2⋅√3⋅atan(3⋅√3)

- ──── + ───────────────

9 3

>>> sympy.simplify(pf(4) - pf(0))

√3⋅(-π + 6⋅atan(3⋅√3))

──────────────────────

9

>>> sympy.simplify(pf(4) - pf(0)).evalf() #求出精确值，注意，这个精确值是相对的。因为不含未知元，故不需要传递值。

0.989661439612300

>>>

我们还可以对部分分式进行分解。

>>> x , y = sympy.symbols('x y')

>>> sympy.apart(1/((x+2)*(x+1)), x)

1 1

- ───── + ─────

x + 2 x + 1

>>> sympy.apart((x+1)/(x-1))

2

1 + ─────

x - 1

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

x + y

─────

x⋅y

>>>

(2)极限。

>>> x = sympy.Symbol('x')

>>> sympy.limit(sympy.tan(x)/x , x , 0) # x趋近于0的极限

1

>>> sympy.limit(sympy.tan(x)/x , x , sympy.oo) #x趋近于无穷的极限

⎛tan(x)⎞

lim ⎜──────⎟

x─→∞⎝ x ⎠

>>> sympy.limit(x,x,0)

0

>>> sympy.limit(x,x,sympy.oo)

∞

>>>

(3)微分。

>>> x = sympy.Symbol('x')

>>> sympy.diff(sympy.sin(x)) #单个参数时可以不指定对哪个变量求导

cos(x)

>>> sympy.diff(sympy.sin(x),x)

cos(x)

>>> sympy.diff(sympy.sin(x),x , 10) #第三个参数表示求几阶导数。

-sin(x)

>>> sympy.diff(sympy.sin(x),x , 2)

-sin(x)

>>> sympy.diff(sympy.sin(x),x , 1)

cos(x)

>>> y = sympy.Symbol('y')

>>> sympy.diff(sympy.sin(x+y) , x) #两个参数时必须指定对哪个参数求导

cos(x + y)

>>> sympy.diff(sympy.sin(x+y) , y)

cos(x + y)

>>> sympy.limit(sympy.sin(x**2+y) , x , 0)

sin(y)

>>>

(4)泰勒展开。

>>> x = sympy.Symbol('x')

>>> sympy.cos(x).series(x,0,10) #求x的函数在0处的展开，10为10阶无穷小项。

2 4 6 8

x x x x ⎛ 10⎞

1 - ── + ── - ─── + ───── + O⎝x ⎠

2 24 720 40320

>>> sympy.cos(x).series(x,1,10)

2 3 4

(x - 1) ⋅cos(1) (x - 1) ⋅sin(1) (x - 1) ⋅cos(1)

cos(1) - (x - 1)⋅sin(1) - ─────────────── + ─────────────── + ───────────────

2 6 24

5 6 7 8

(x - 1) ⋅sin(1) (x - 1) ⋅cos(1) (x - 1) ⋅sin(1) (x - 1) ⋅cos(1) (x -

- ─────────────── - ─────────────── + ─────────────── + ─────────────── - ────

120 720 5040 40320

9

1) ⋅sin(1) ⎛ 10 ⎞

─────────── + O⎝(x - 1) ; x → 1⎠

362880

>>> y = sympy.Symbol('y')

>>>

>>> sympy.cos(x+y).series(x,0,10)

2 3 4 5 6

x ⋅cos(y) x ⋅sin(y) x ⋅cos(y) x ⋅sin(y) x ⋅cos(y)

cos(y) - x⋅sin(y) - ───────── + ───────── + ───────── - ───────── - ─────────

2 6 24 120 720

7 8 9

x ⋅sin(y) x ⋅cos(y) x ⋅sin(y) ⎛ 10⎞

+ ───────── + ───────── - ───────── + O⎝x ⎠

5040 40320 362880

>>> sympy.cos(x+y).series(x,-1,4)

2 3

(x + 1) ⋅cos(y - 1) (x + 1) ⋅sin(y - 1)

cos(y - 1) - (x + 1)⋅sin(y - 1) - ─────────────────── + ─────────────────── +

2 6

⎛ 4 ⎞

O⎝(x + 1) ; x → -1⎠

>>>

(5)求和。

>>> i , m , n = sympy.symbols('i m n' , integer = True) #integer表示取值为正数

>>> sympy.summation( i , (i , 1,n) )

2

n n

── + ─

2 2

>>> sympy.summation( 1/i , (i , 1,sympy.oo) )　#参数意义为summation(符号 ， (要进行求和的符号 ， 下限，上限))

∞

>>> sympy.summation( 1/i**2 , (i , 1,sympy.oo) )

2

π

──

6

>>> sympy.summation(1/(i**2 + 1) , (i , 0,n) , (n , 0,m))

m n

____ ____

╲ ╲

╲ ╲ 1

╲ ╲ ──────

╱ ╱ 2

╱ ╱ i + 1

╱ ╱

‾‾‾‾ ‾‾‾‾

n = 0 i = 0

>>> sympy.summation(1 , (i , 0,n) , (n , 0,m))

2

m 3⋅m

── + ─── + 1

2 2

>>> sympy.summation(i**2 , (i , 0,n) , (n , 0,m))

4 3 2

m m 5⋅m m

── + ── + ──── + ─

12 3 12 6

>>>

(6)积分。

>>> x , y = sympy.symbols('x y')

>>> sympy.integrate(x**2) #普通积分

3

x

──

3

>>> sympy.integrate(x**2 , x)

3

x

──

3

>>> sympy.integrate(sympy.exp(-x**2) , (x , -1 , 1)) #定积分 ， 第二个参数为(符号，下限，上限)

√π⋅erf(1)

>>> sympy.integrate(sympy.sin(x**2) , (x , -1 , 1))

⎛√2⎞

3⋅√2⋅√π⋅fresnels⎜──⎟⋅Γ(3/4)

⎝√π⎠

───────────────────────────

4⋅Γ(7/4)

>>> sympy.integrate(sympy.sin(x) , (x , -1 , 1))

0

>>> sympy.integrate(sympy.cos(x) , (x , -1 , 1))

2⋅sin(1)

>>> sympy.integrate(sympy.cos(x) , (x , -1 , 1)).evalf()

1.68294196961579

>>> sympy.integrate(sympy.exp(-x) , (x , 0 , sympy.oo)) #求反常积分

1

>>>

(7)微分方程。

>>> x , y = sympy.symbols('x y')

>>> f = sympy.Function('f') #函数对象

>>> f #函数

f

>>> f(x).series(x,0,n=4) #关于x的函数f(x)在0处展开到4阶。

⎛ 2 ⎞│ ⎛ 3 ⎞│

2 ⎜ d ⎟│ 3 ⎜ d ⎟│

x ⋅⎜───(f(x))⎟│ x ⋅⎜───(f(x))⎟│

⎜ 2 ⎟│ ⎜ 3 ⎟│

⎛d ⎞│ ⎝dx ⎠│x=0 ⎝dx ⎠│x=0 ⎛ 4⎞

f(0) + x⋅⎜──(f(x))⎟│ + ────────────────── + ────────────────── + O⎝x ⎠

⎝dx ⎠│x=0 2 6

>>> f(x)

f(x)

>>> f(x).diff(x) # f是x的函数，并对x求一次导数

d

──(f(x))

dx

>>> f(x).diff(x,x) # f是x的函数，并对x求二次导数

2

d

───(f(x))

2

dx

>>> f(x).diff(x,x,x) # f是x的函数，并对x求三次导数

3

d

───(f(x))

3

dx

>>> f(x).diff(x,x,x,x)

4

d

───(f(x))

4

dx

>>> f(x,y).diff(x,y,x) #多元函数求导

3

∂

──────(f(x, y))

2

∂y ∂x

>>> f(x,y).diff(y,x,y)

3

∂

──────(f(x, y))

2

∂y ∂x

>>>

(8)解方程组。

>>> x , y = sympy.symbols('x y')

>>> sympy.solve([x - y - 8 , x + y] , [x,y]) #解方程，x-y-8 = 0,x+y = 0.

{x: 4, y: -4}

>>>

(9)线性代数。

1、一些复杂的符号。

>>> sympy.symbols('A:Z')

(A, B, C, D, E, F, G, H, I, J, K, L, M, N, O, P, Q, R, S, T, U, V, W, X, Y, Z)

>>> sympy.symbols('x:3')

(x₀, x₁, x₂)

>>> sympy.symbols('A:3(1:9)')

(A₀₁, A₀₂, A₀₃, A₀₄, A₀₅, A₀₆, A₀₇, A₀₈, A₁₁, A₁₂, A₁₃, A₁₄, A₁₅, A₁₆, A₁₇, A₁

₈, A₂₁, A₂₂, A₂₃, A₂₄, A₂₅, A₂₆, A₂₇, A₂₈)

>>> sympy.symbols('A1:3(1:9)')

(A₁₁, A₁₂, A₁₃, A₁₄, A₁₅, A₁₆, A₁₇, A₁₈, A₂₁, A₂₂, A₂₃, A₂₄, A₂₅, A₂₆, A₂₇, A₂

₈)

>>> sympy.symbols('A1:3(1:3)(1:3)')

(A₁₁₁, A₁₁₂, A₁₂₁, A₁₂₂, A₂₁₁, A₂₁₂, A₂₂₁, A₂₂₂)

>>> x = sympy.symbols('x:3')

>>> x

(x₀, x₁, x₂)

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

>>> [sympy.diff(f(*x) , xx) for xx in x]

⎡ ∂ ∂ ∂ ⎤

⎢───(f(x₀, x₁, x₂)), ───(f(x₀, x₁, x₂)), ───(f(x₀, x₁, x₂))⎥

⎣∂x₀ ∂x₁ ∂x₂ ⎦

>>>

2、矩阵。

>>> sympy.Matrix([[1,0] , [0,1]])

⎡1 0⎤

⎢ ⎥

⎣0 1⎦

>>> x , y = sympy.symbols('x,y')

>>> sympy.Matrix([[sympy.sin(x) , x] , [sympy.cos(y) , y]])

⎡sin(x) x⎤

⎢ ⎥

⎣cos(y) y⎦

>>> x = sympy.symbols('x:2')

>>> F = sympy.Lambda(x , sympy.Matrix([sympy.sin(x[0]) + sympy.cos(2*x[1]) , sympy.sin(x[0]*sympy.cos(x[1]))]))

>>> F

⎡sin(x₀) + cos(2⋅x₁)⎤

(x₀, x₁) ↦ ⎢ ⎥

⎣ sin(x₀⋅cos(x₁)) ⎦

>>> F(*x).jacobian(x)

⎡ cos(x₀) -2⋅sin(2⋅x₁) ⎤

⎢ ⎥

⎣cos(x₁)⋅cos(x₀⋅cos(x₁)) -x₀⋅sin(x₁)⋅cos(x₀⋅cos(x₁))⎦

>>> phi = sympy.symbols('phi') #可以使用希腊字母

>>> rot = sympy.Matrix([[sympy.cos(phi) , -sympy.sin(phi)],[sympy.sin(phi) , sympy.cos(phi)]])

>>> rot

⎡cos(φ) -sin(φ)⎤

⎢ ⎥

⎣sin(φ) cos(φ) ⎦

>>> rot.T * rot

⎡ 2 2 ⎤

⎢sin (φ) + cos (φ) 0 ⎥

⎢ ⎥

⎢ 2 2 ⎥

⎣ 0 sin (φ) + cos (φ)⎦

>>> sympy.simplify(rot.T * rot)

⎡1 0⎤

⎢ ⎥

⎣0 1⎦

>>> sympy.simplify(rot.T * rot - sympy.eye(2)) #sympy.eye生成一个单位矩阵

⎡0 0⎤

⎢ ⎥

⎣0 0⎦

>>> rot.inv() #rot.inv() 对矩阵求逆

⎡ 2 ⎤

⎢1 - sin (φ) ⎥

⎢─────────── sin(φ)⎥

⎢ cos(φ) ⎥

⎢ ⎥

⎣ -sin(φ) cos(φ)⎦

>>> rot.T

⎡cos(φ) sin(φ)⎤

⎢ ⎥

⎣-sin(φ) cos(φ)⎦

>>> sympy.simplify(rot.T - rot.inv())

⎡0 0⎤

⎢ ⎥

⎣0 0⎦

>>>

还可以以如下方式来生成矩阵。

>>> sympy.Matrix(3,4,sympy.symbols('M:3(:4)'))

⎡M₀₀ M₀₁ M₀₂ M₀₃⎤

⎢ ⎥

⎢M₁₀ M₁₁ M₁₂ M₁₃⎥

⎢ ⎥

⎣M₂₀ M₂₁ M₂₂ M₂₃⎦

>>> M = sympy.Matrix(3,4,sympy.symbols('M:3(:4)'))

>>> M[1:]

[M₀₁, M₀₂, M₀₃, M₁₀, M₁₁, M₁₂, M₁₃, M₂₀, M₂₁, M₂₂, M₂₃]

>>> M[1,:] #与numpy数组类似

[M₁₀ M₁₁ M₁₂ M₁₃]

>>> def toe(n):

... a = sympy.symbols('a:' + str(n * 2))

... f = lambda i , j : a[i - j + n - j]

... return sympy.Matrix(n,n,f)

...

>>> toe(5)

⎡a₅ a₃ a₁ a₉ a₇⎤

⎢ ⎥

⎢a₆ a₄ a₂ a₀ a₈⎥

⎢ ⎥

⎢a₇ a₅ a₃ a₁ a₉⎥

⎢ ⎥

⎢a₈ a₆ a₄ a₂ a₀⎥

⎢ ⎥

⎣a₉ a₇ a₅ a₃ a₁⎦

>>>

还可以对矩阵进行分解等操作。

>>> A = sympy.Matrix(3,3,sympy.symbols('A1:4(1:4)'))

>>> b = sympy.Matrix(3,1, sympy.symbols('b1:4'))

>>> A

⎡A₁₁ A₁₂ A₁₃⎤

⎢ ⎥

⎢A₂₁ A₂₂ A₂₃⎥

⎢ ⎥

⎣A₃₁ A₃₂ A₃₃⎦

>>> B

Traceback (most recent call last):

File "", line 1, in

NameError: name 'B' is not defined

>>> b

⎡b₁⎤

⎢ ⎥

⎢b₂⎥

⎢ ⎥

⎣b₃⎦

>>> x = A.LUsolve(b) # 分解之后字符非常多就不展示了。

>>> sympy.simplify(x) #输出还是有点问题

⎡ A₁₂⋅A₂₃⋅b₃ - A₁₂⋅A₃₃⋅b₂ - A₁₃⋅A₂₂⋅b₃ + A₁₃⋅A₃₂⋅b₂ + A₂₂⋅A₃₃⋅b₁ - A₂₃⋅A₃₂⋅b

⎢─────────────────────────────────────────────────────────────────────────────

⎢A₁₁⋅A₂₂⋅A₃₃ - A₁₁⋅A₂₃⋅A₃₂ - A₁₂⋅A₂₁⋅A₃₃ + A₁₂⋅A₂₃⋅A₃₁ + A₁₃⋅A₂₁⋅A₃₂ - A₁₃⋅A₂₂

⎢

⎢ -A₁₁⋅A₂₃⋅b₃ + A₁₁⋅A₃₃⋅b₂ + A₁₃⋅A₂₁⋅b₃ - A₁₃⋅A₃₁⋅b₂ - A₂₁⋅A₃₃⋅b₁ + A₂₃⋅A₃₁⋅

⎢─────────────────────────────────────────────────────────────────────────────

⎢A₁₁⋅A₂₂⋅A₃₃ - A₁₁⋅A₂₃⋅A₃₂ - A₁₂⋅A₂₁⋅A₃₃ + A₁₂⋅A₂₃⋅A₃₁ + A₁₃⋅A₂₁⋅A₃₂ - A₁₃⋅A₂₂

⎢

⎢ A₁₁⋅A₂₂⋅b₃ - A₁₁⋅A₃₂⋅b₂ - A₁₂⋅A₂₁⋅b₃ + A₁₂⋅A₃₁⋅b₂ + A₂₁⋅A₃₂⋅b₁ - A₂₂⋅A₃₁⋅b

⎢─────────────────────────────────────────────────────────────────────────────

⎣A₁₁⋅A₂₂⋅A₃₃ - A₁₁⋅A₂₃⋅A₃₂ - A₁₂⋅A₂₁⋅A₃₃ + A₁₂⋅A₂₃⋅A₃₁ + A₁₃⋅A₂₁⋅A₃₂ - A₁₃⋅A₂₂

₁ ⎤

────⎥

⋅A₃₁⎥

⎥

b₁ ⎥

────⎥

⋅A₃₁⎥

⎥

₁ ⎥

────⎥

⋅A₃₁⎦

>>>

(10)符号赋值。

>>> x , a = sympy.symbols('x a')

>>> f = x + a

>>> nf = f.subs(x,0) #对x赋值0

>>> nf

a

>>> f.subs(a,1)

x + 1

>>> f.subs(a , a**2) # 将a赋值为a的平方

2

a + x

>>> f.subs({x:a}) #将x赋值为a

2⋅a

>>> f.subs({x:1}) #将x赋值为1

a + 1

>>> f.subs({x:a/x , a:a**x+x})

a

─

2⋅a x

─── + a

x

>>> f.subs({x:a/x , a:a+x})

2⋅a

a + ───

x

>>> f.subs(f,f**2) #也可以对整个式子进行替换

2

(a + x)

>>> f.subs(f,f**2 + f)

2

a + x + (a + x)

>>> M = sympy.Matrix(3,3 , sympy.symbols('A:3(:3)'))

>>> M.subs(M[1,1],1)

⎡A₀₀ A₀₁ A₀₂⎤

⎢ ⎥

⎢A₁₀ 1 A₁₂⎥

⎢ ⎥

⎣A₂₀ A₂₁ A₂₂⎦

>>>

以上就是符号计算的内容，因为内容较多，又是一口气弄得，有可能有的地方有疏漏，欢迎指正。