# Python 符號運算套件：SymPy
> 作者：王一哲 日期：2018/11/29
<br></br>
## 前言
SymPy 是 Python 的符號運算套件，可以幫助我們計算函數的微分、積分，而且套件裡已經內建了常用的常數和函數，例如圓周率$\pi = 3.14159265358979$、尤拉常數$E = 2.71828182845905$、正弦函數$\sin$。如果已經在電腦上安裝好 Python，接下來只要在文字介面中執行下的指令即可安裝 SymPy
```shell
pip install sympy
```
如果是 Windows 的用戶，可以使用 **Windows PowerShell(系統管理員)** 替這部電腦所有的使用者安裝，如果是 Linux 的用戶，在執行指令時可能需要在最前面加上 **sudo** 暫時取得管理者權限。如果只想要試用看看，可以到 SymPy 的官方網站按右上角的 **Online Shell** 使用線上版的 SymPy Live。
<br></br>
<img height="100%" width="100%" src="https://i.imgur.com/VWQZt0T.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">SymPy 官方網站</div>
<br></br>
<img height="100%" width="100%" src="https://i.imgur.com/JOF93JE.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">SymPy Live</div>
<br></br>
## 基本運算
如果想要引入 SymPy 函式庫，需要在程式碼的最上方加上
```python
from __future__ import division
from sympy import *
```
其中第一行是為了引入 SymPy 的除法型式，第二行才是引入 SymPy，兩行的順序不能顛倒，以下的例子會省略這兩行。
SymPy 的數字格式有 3 種：整數(Integer)、浮點數(Float)、有理數(Rational)，其中有理數是以分數的型式顯示，例如
```python
print("Rational(2, 3) =", Rational(2, 3))
print("Rational(2) / 3) =", Rational(2) / 3)
print("2 / 3 =", 2 / 3)
```
輸出的結果為
```python
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$$
程式碼為
```python
print("sin(pi / 3) =", sin(pi / 3))
print("sin(pi / 3) =", sin(pi / 3).evalf())
print("exp(-oo) =", exp(-oo))
```
輸出的結果為
```python
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$$
程式碼為
```python
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 代入。
輸出的結果為
```python
f(x) = x**2 - 2*x - 8
roots of f(x): {4: 1, -2: 1}
roots of f(2): -8.00000000000000
```
我們也可以用 SymPy 解聯立方程式，使用的函式為 **solve**，語法為
```python
solve((方程式1, 方程式2), (變數1, 變數2))
```
以下列的二元一次聯立方程式為例
$$\begin{matrix}x + 2y = 8 \\ 2x - y = 6 \end{matrix} \Rightarrow x = 4, y = 2$$
程式碼為
```python
x, y = symbols('x y', communtative = True)
print(solve((x + 2*y - 8, 2*x - y - 6), (x, y)))
```
輸出為
```python
{x: 4, y: 2}
```
<br></br>
## 微分
如果要用 SymPy 計算微分，所用的函式為 **diff**，語法為
```python
函數名稱.diff()
diff(函數名稱, 變數)
```
兩種寫法的效果相同，但我比較傾向使用第二種寫法，能夠一看就知道是對誰微分。以多項式函數的微分為例：
```python
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))
```
輸出為
```python
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)$$
程式碼的寫法如下
```python
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)
```
輸出為
```python
x(t) = R*cos(phi + t*w)
v(t) = -R*w*sin(phi + t*w)
a(t) = -R*w**2*cos(phi + t*w)
```
<br></br>
## 積分
如果要用 SymPy 計算積分，所用的函式為 **integrate**，語法為
```python
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$$
程式碼為
```python
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())
```
輸出為
```python
-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$。
接下來我們試著處理更複雜的函數，例如前一篇文章〈[利用殼層定理推導電量均勻分布的球殼產生的電場](https://hackmd.io/s/SJHCp1OR7)〉當中的積分式
$$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$$
使用的程式碼如下
```python
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，會變成
```python
E_in = Q*k/(2*r**2) - Q*k*(R - r - (-R + r)*(R + r)/(R - r))/(4*R*r**2)
```
<br></br>
## 結語
這是最近幾天為了做符號運算找到的使用方法，只運用到 SymPy 一點點的功能而已，有興趣的同學可以再上網搜尋更多的教學及範例。
<br></br>
## 參考資料
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`
Sign in with Wallet
Connect another wallet