Python 符號運算套件：SymPy

作者：王一哲 日期：2018/11/29

前言

SymPy 是 Python 的符號運算套件，可以幫助我們計算函數的微分、積分，而且套件裡已經內建了常用的常數和函數，例如圓周率

$\pi =3.14159265358979$、尤拉常數

$E=2.71828182845905$、正弦函數

$\sin$。如果已經在電腦上安裝好 Python，接下來只要在文字介面中執行下的指令即可安裝 SymPy

pip install sympy

如果是 Windows 的用戶，可以使用 Windows PowerShell(系統管理員) 替這部電腦所有的使用者安裝，如果是 Linux 的用戶，在執行指令時可能需要在最前面加上 sudo 暫時取得管理者權限。如果只想要試用看看，可以到 SymPy 的官方網站按右上角的 Online Shell 使用線上版的 SymPy Live。

SymPy 官方網站

SymPy Live

基本運算

如果想要引入 SymPy 函式庫，需要在程式碼的最上方加上

from __future__ import division
from sympy import *

其中第一行是為了引入 SymPy 的除法型式，第二行才是引入 SymPy，兩行的順序不能顛倒，以下的例子會省略這兩行。

SymPy 的數字格式有 3 種：整數(Integer)、浮點數(Float)、有理數(Rational)，其中有理數是以分數的型式顯示，例如

print("Rational(2, 3) =", Rational(2, 3))
print("Rational(2) / 3) =", Rational(2) / 3)
print("2 / 3 =", 2 / 3)

輸出的結果為

Rational(2, 3) = 2/3
Rational(2) / 3) = 2/3
2 / 3 = 0.6666666666666666

我們也可以把數字代入內建的函數中，輸出時會保留數學上的表達方式，如果要把數值計算出來，則要加上 evalf()，例如

$$\sin \left(\frac{\pi }{3}\right)=\frac{\sqrt{3}}{2}\approx 0.866025403784439$$

$$e^{-\mathrm{\infty }}=\exp [-\mathrm{\infty }]=0$$

程式碼為

print("sin(pi / 3) =", sin(pi / 3))
print("sin(pi / 3) =", sin(pi / 3).evalf())
print("exp(-oo) =", exp(-oo))

輸出的結果為

sin(pi / 3) = sqrt(3)/2
sin(pi / 3) = 0.866025403784439
exp(-oo) = 0

我們也可以自訂函數，例如

$$f(x)=x^2-2x-8$$

$$f(x)=(x-4)(x+2)=0\Rightarrow x=4\vee -2$$

$$f(2)=2^2-2\times 2-8=-8$$

程式碼為

x = symbols('x', communtative = True)
f = symbols('f', cls = Function)
f = x**2 - 2*x - 8
print("f(x) =", f)
print("roots of f(x):", roots(f))
print("roots of f(2):", f.evalf(subs = {x:2}))

我們用到以下的函式

1. symbols: 定義所用的符號，communtative = True 將符號設定為可交換，cls = Function 是將這些符號設定為函數的代號。
2. roots: 求多項式函數的根。
3. subs = {x:2}: 利用字典格式將函數中的 x 用 2 代入。

輸出的結果為

f(x) = x**2 - 2*x - 8
roots of f(x): {4: 1, -2: 1}
roots of f(2): -8.00000000000000

我們也可以用 SymPy 解聯立方程式，使用的函式為 solve，語法為

solve((方程式1, 方程式2), (變數1, 變數2))

以下列的二元一次聯立方程式為例

$$\begin{array}{c}x+2y=8\\ 2x-y=6\end{array}\Rightarrow x=4,y=2$$

程式碼為

x, y = symbols('x y', communtative = True)
print(solve((x + 2*y - 8, 2*x - y - 6), (x, y)))

輸出為

{x: 4, y: 2}

微分

如果要用 SymPy 計算微分，所用的函式為 diff，語法為

函數名稱.diff()
diff(函數名稱, 變數)

兩種寫法的效果相同，但我比較傾向使用第二種寫法，能夠一看就知道是對誰微分。以多項式函數的微分為例：

x = symbols('x', communtative = True)
f = symbols('f', cls = Function)
f = x**2 - 2*x - 8
print("f'(x) =", f.diff())
print("f'(x) =", diff(f, x))

輸出為

f'(x) = 2*x - 2

如果以簡諧運動的位置

$x(t)$、速度

$v(t)$、加速度

$a(t)$ 為例：

$$x(t)=R\cos (\omega t+\varphi )$$

$$v(t)=-R\omega \sin (\omega t+\varphi )$$

$$x(t)=-R{\omega }^2\cos (\omega t+\varphi )$$
程式碼的寫法如下

R, w, t, phi = symbols('R w t phi', commutative = True)
x, v, a = symbols('x v a', cls = Function)
x = R * cos(w * t + phi)
print("x(t) =", x)
v = diff(x, t)
a = diff(v, t)
print("v(t) =", v)
print("a(t) =", a)

輸出為

x(t) = R*cos(phi + t*w)
v(t) = -R*w*sin(phi + t*w)
a(t) = -R*w**2*cos(phi + t*w)

積分

如果要用 SymPy 計算積分，所用的函式為 integrate，語法為

integrate(函數名稱, 變數)
integrate(函數名稱, (變數, 積分下限, 積分上限))

以下以一些常見的函數為例：

$$\int \sin (x)dx=-\cos (x)+c$$

$${\int }_0^{\pi }\sin (x)dx={[-\cos (x)]}_0^{\pi }=-(-1)-(-1)=2$$

$$\int e^{x}dx=e^x+c$$

$${\int }_0^1{e}^{x}dx=e^1-e^0=e-1\approx 1.71828182845905$$

$$\int \frac{dx}{x}=\ln x+c$$

$${\int }_2^3\frac{dx}{x}=\ln 3-\ln 2=\ln \left(\frac{3}{2}\right)\approx 0.405465108108164$$

程式碼為

x = symbols('x', communtative = True)
print(integrate(sin(x), x))
print(integrate(sin(x), (x, 0, pi)))
print(integrate(exp(x), x))
print(integrate(exp(x), (x, 0, 1)))
print(integrate(exp(x), (x, 0, 1)).evalf())
print(integrate(1 / x, x))
print(integrate(1 / x, (x, 2, 3)))
print(integrate(1 / x, (x, 2, 3)).simplify())
print(integrate(1 / x, (x, 2, 3)).evalf())

輸出為

-cos(x)
2
exp(x)
-1 + E
1.71828182845905
log(x)
-log(2) + log(3)
-log(2) + log(3)
0.405465108108164

在倒數第2行的程式碼中，我們使用了 simplify()，用來化簡數學式子，但是在 Python Shell 中看不出效果，如果在 SymPy Live 中會將 -log(2) + log(3) 合併為 log(3 / 2)。還有一點，在SymPy 當中

$\log ={\log }_e=\ln$。

接下來我們試著處理更複雜的函數，例如前一篇文章〈利用殼層定理推導電量均勻分布的球殼產生的電場〉當中的積分式

$$E_{out}=\int d{E}_r=\frac{kQ}{4r^{2}R}{\int }_{r-R}^{r+R}\frac{r^2+s^2-R^2}{s^2}ds=\frac{kQ}{r^2}$$

$$E_{in}=\int d{E}_r=\frac{kQ}{4r^{2}R}{\int }_{R-r}^{R+r}\frac{r^2+s^2-R^2}{s^2}ds=0$$

使用的程式碼如下

k, Q, r, R, s = symbols('k Q r R s', commutative = True)
dE = k * Q / (4 * r**2 * R) * ((r**2 + s**2 - R**2) / s**2)
E_out = integrate(dE, (s, r - R, r + R), manual = True)
E_in = integrate(dE, (s, R - r, R + r), manual = True).simplify()
print("E_out =", E_out)
print("E_in =", E_in)

我們在積分指令最後加上 manual = True，可以將計算結果盡量顯示為人類手算積分的表達方式。另外在 E_in 這行最後加上 simplify() 將計算結果化簡，否則答案不會是0，會變成

E_in = Q*k/(2*r**2) - Q*k*(R - r - (-R + r)*(R + r)/(R - r))/(4*R*r**2)

結語

這是最近幾天為了做符號運算找到的使用方法，只運用到 SymPy 一點點的功能而已，有興趣的同學可以再上網搜尋更多的教學及範例。

參考資料

1. SymPy 官方網站: https://www.sympy.org/en/index.html
2. SymPy Live: https://live.sympy.org/
3. SymPy Core: https://docs.sympy.org/latest/modules/core.html
4. Calculus: https://docs.sympy.org/latest/tutorial/calculus.html
5. Symbolic Integrals: https://docs.sympy.org/latest/modules/integrals/integrals.html

---

tags:Python