Symbolic Mathematics in Python with SymPy

• 前言
• Using SymPy as a calculator
  • 1. SymPy 的基本運算
    • 除法
    • Rational class
    • Float class
    • SymPy 中的常數
    • 表達特定位數的浮點數
    • SymPy 中的無限大
• Defining symbols and symbolic operations
  • 定義符號
  • 定義多個符號
  • 展開式
  • 因式分解
  • x次方由小到大排列
  • 帶入值
  • 簡化
• Solving equations symbolically
  • 解方程式
    • 字典形式
  • 解方程式組
    • 解的形式
• Plotting using SymPy
• Calculus with SymPy
  • 極限
    • 極限的實值
    • 負極限
  • 微分
    • 微分的結果
    • 帶入
  • Finding the integrals of functions
    • 積分
    • 定積分
  • 序列展開
  • lambafy()
• Linear algebra with SymPy
  • 對角
  • 矩陣運算
    • 相加
    • 相乘
  • 反矩陣
  • 找特徵值
  • 找特徵向量

前言

• SymPy 是一款將編碼轉為可執行數學演算的 Python 套件。以其易執行的可做為的替代品，例如 Mathematica, Maple, Sage 等等。
• SymPy 基本上全部建置在 Python 裡並不需要載額外封包。

Using SymPy as a calculator

import sympy as sp

1. SymPy 的基本運算

除法

u = sp.Integer(24)
u/8, u/7
## Output: (3, 24/7)

Rational class

用分數的形式呈現

a = sp.Rational(1, 2)
a, a*2
## Output: (1/2, 1)

Float class

用浮點數的形式呈現

sp.Float(3.53) + sp.Float(3)
## Output: 6.53

在一般的 Python 中，浮點數的運算會有誤差，但在 SymPy 中，浮點數的運算會是精確的。

3.53 + 3
## Output: 6.529999999999999

SymPy 中的常數

SymPy 中的常數是以 sp.常數名稱 的形式呈現，例如 sp.pi、sp.E、sp.I 等等。
如果要將常數轉為浮點數，可以使用 evalf() 函數。

sp.pi**2, sp.pi.evalf()
## Output: (pi**2, 3.14159265358979)
(sp.pi + sp.exp(1)).evalf()
## Output: 5.85987448204884

表達特定位數的浮點數

sqrt2 = sp.sqrt(2).evalf(10)
sqrt2
## Output: 1.414213562

SymPy 中的無限大

sp.oo

Defining symbols and symbolic operations

定義符號

可以使用 Symbol() 函數來定義符號

x = sp.Symbol('x')

定義多個符號

定義多個符號時，可以使用 symbols() 函數

x,y = sp.symbols('x,y')

展開式

p = (x + 2)*(x + 3)
p2 = sp.expand(p)
## Output: x**2 + 5*x + 6

因式分解

sp.factor(p2)
## Output: (x + 2)*(x + 3)

x次方由小到大排列

sp.init_printing(order='rev-lex')
p2
## Output: 6 + 5*x + x**2

sp.init_printing(order='lex')
x = sp.Symbol('x')
y = sp.Symbol('y')
p3 = x*x + x*y + x*y + y*y
p3
## Output: x**2 + 2*x*y + y**2

帶入值

p3.subs({x:1, y:2}) # We use dictionary to specify the values of x and y
## Output: 9

簡化

用simplify()函數來簡化

sp.simplify(p3.subs({x:1-y})) 
# The result turns out to be 1 because the other terms of the expression cancel each other.
## Output: 1

Solving equations symbolically

解方程式

solve()函數可以解方程式

x = sp.Symbol('x')
expr = x - 5 - 7
sp.solve(expr)
## Output: [12]

字典形式

用dict=True得到字典形式

x = sp.Symbol('x')
expr = x**2 + 5*x + 4
sp.solve(expr), sp.solve(expr, dict=True)
## Output: ([-4, -1], [{x: -4}, {x: -1}])

解方程式組

x = sp.Symbol('x')
y = sp.Symbol('y')
expr1 = 2*x + 3*y - 6
expr2 = 3*x + 2*y - 12

sp.solve((expr1, expr2), dict=True)
## Output: [{x: 24/5, y: -6/5}]

解的形式

nsolve() -> 帶著符號的解
nroots() -> 數值解

expr1 = x**3 + 2*x**2 -1
sp.solve(expr1), sp.nroots(expr1, n=3) # n is the number of digits to calculate

Plotting using SymPy

可以用plot()函數來畫圖
加上設定title、xlabel、ylabel等

x=sp.Symbol('x')
sp.plot(2*x + 3, (x, -5, 5), title='A Line', xlabel='x', ylabel='2x+3');

Calculus with SymPy

極限

x = sp.Symbol('x')
l = sp.Limit(1/x, x, sp.oo) # (function, variable, limit)

極限的實值

l.doit()
## Output: 0

負極限

sp.Limit(1/x, x, 0, dir='-').doit()
## Output: -oo

微分

t = sp.Symbol('t')
St = 5*t**2 + 2*t + 8
sp.Derivative(St, t)
## Output: (d/dt)(5*t**2 + 2*t + 8, t)

微分的結果

sp.Derivative(St, t).doit()
## Output: 10*t + 2

我們也可以直接使用 diff() 函數來進行微分

sp.diff(St, t)
## Output: 10*t + 2

帶入

t1 = sp.Symbol('t1')
d.doit().subs({t:t1}), d.doit().subs({t:1})
## Output: (10*t1 + 2, 12)

Finding the integrals of functions

積分

x = sp.Symbol('x')
sp.Integral(sp.sin(x), x)
## Output: Integral(sin(x), x)dx

積分的結果

sp.Integral(sp.sin(x), x).doit()
## Output: -cos(x)

同樣的，我們也可以直接使用 integrate() 函數來進行積分

sp.integrate(sp.sin(x), x), sp.integrate(sp.log(x), x)
## Output: (-cos(x), x*log(x) - x)

定積分

要計算定積分，可以使用 Integral() 函數，並且指定積分的範圍

# sp.Integral(sp.sin(x), (x, 0, sp.pi / 2)).doit()
sp.integrate(sp.sin(x), (x, 0, sp.pi / 2))
## Output: 1

反常積分也可以使用

sp.integrate(sp.exp(-x), (x, 0, sp.oo))
## Output: 1

序列展開

我們可以使用 series() 函數來展開序列

sp.series(sp.cos(x), x, 0, 10)
## Output: 1 - x**2/2 + x**4/24 - x**6/720 + x**8/40320 + O(x**10)

也可以移除 O ($x^{10}$) 的部分

sp.series(sp.cos(x), x, 0, 10).removeO()
## Output: 1 - x**2/2 + x**4/24 - x**6/720 + x**8/40320

lambafy()

lambafy() 函數可以將符號函數轉換成數值函數

expr = sp.series(sp.cos(x), x, 0, 10).removeO()
f = sp.utilities.lambdify(x, expr, 'numpy')
f(0)
## Output: 1.0

Linear algebra with SymPy

建立矩陣

m1 = sp.Matrix([[a,b,c], [x, y, z], [a, b, c]])
m1
## Output: Matrix([
## [a, b, c],
## [x, y, z],
## [a, b, c]])

對角

用元素建立對角矩陣

m2 = sp.diag(a,b,c)
m2
## Output: Matrix([
## [a, 0, 0],
## [0, b, 0],
## [0, 0, c]])

矩陣運算

相加

3*m1 + m2

相乘

m1*m2.T

反矩陣

m2.inv()

找特徵值

m2.eigenvals()
## Output: {a: 1, b: 1, c: 1}

找特徵向量

m2.eigenvects()