Python 符號運算套件：SymPy

前言

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^{-\infty} = \exp[-\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{matrix}x + 2y = 8 \\ 2x - y = 6 \end{matrix} \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 + \phi)$$
$$v(t) = -R \omega \sin(\omega t + \phi)$$
$$x(t) = -R \omega^2 \cos(\omega t + \phi)$$

程式碼的寫法如下

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 = \left [ -\cos(x) \right ]_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 dE_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 dE_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

---

HackMD 版本連結：https://hackmd.io/s/Syj4ncoCm