3.2 Sympy：Python中的符號數學 · SciPy Lecture Notes 中文版

# 3.2 Sympy：Python中的符號數學 > **作者** : Fabian Pedregosa **目的** * 從任意的精度評估表達式。 * 在符號表達式上進行代數運算。 * 用符號表達式進行基本的微積分任務 (極限、微分法和積分法)。 * 求解多項式和超越方程。 * 求解一些微分方程。為什么是SymPy? SymPy是符號數學的Python庫。它的目的是成為Mathematica或Maple等系統的替代品，同時讓代碼盡可能簡單并且可擴展。SymPy完全是用Python寫的，并不需要外部的庫。 Sympy文檔及庫安裝見[http://www.sympy.org/](http://www.sympy.org/) **章節內容** * SymPy第一步 * 使用SymPy作為計算器 * 練習 * 符號 * 代數運算 * 展開 * 化簡 * 微積分 * 極限 * 微分法 * 序列擴展 * 積分法 * 練習 * 方程求解 * 練習 * 線性代數 * 矩陣 * 微分方程 ## 3.2.1 SymPy第一步 ### 3.2.1.1 使用SymPy作為計算器 SymPy定義了三種數字類型：實數、有理數和整數。有理數類將有理數表征為兩個整數對: 分子和分母，因此Rational(1,2)代表1/2, Rational(5,2)代表5/2等等: In?[2]: ``` from sympy import * a = Rational(1,2) ``` In?[2]: ``` a ``` Out[2]: ``` 1/2 ``` In?[3]: ``` a*2 ``` Out[3]: ``` 1 ``` SymPy在底層使用mpmath, 這使它可以用任意精度的算術進行計算。這樣，一些特殊的常數，比如e, pi, oo (無限), 可以被作為符號處理并且可以以任意精度來評估: In?[4]: ``` pi**2 ``` Out[4]: ``` pi**2 ``` In?[5]: ``` pi.evalf() ``` Out[5]: ``` 3.14159265358979 ``` In?[6]: ``` (pi + exp(1)).evalf() ``` Out[6]: ``` 5.85987448204884 ``` 如你所見，將表達式評估為浮點數。也有一個類代表數學的無限, 稱為 oo: In?[7]: ``` oo > 99999 ``` Out[7]: ``` True ``` In?[8]: ``` oo + 1 ``` Out[8]: ``` oo ``` ### 3.2.1.2 練習 * 計算 $\sqrt{2}$ 小數點后一百位。 * 用有理數算術計算1/2 + 1/3 in rational arithmetic. ### 3.2.1.3 符號與其他計算機代數系統不同，在SymPy你需要顯性聲明符號變量: In?[4]: ``` from sympy import * x = Symbol('x') y = Symbol('y') ``` 然后你可以計算他們: In?[10]: ``` x + y + x - y ``` Out[10]: ``` 2*x ``` In?[11]: ``` (x + y)**2 ``` Out[11]: ``` (x + y)**2 ``` 符號可以使用一些Python操作符操作: +, -, *, ** (算術), &, |, ~ , >>, << (布爾邏輯). **打印** 這里我們使用下列設置打印 In?[?]: ``` sympy.init_printing(use_unicode=False, wrap_line=True) ``` ## 3.2.2 代數運算 SymPy可以進行強大的代數運算。我們將看一下最常使用的：展開和化簡。 ### 3.2.2.1 展開使用這個模塊展開代數表達式。它將試著密集的乘方和相乘: In?[13]: ``` expand((x + y)**3) ``` Out[13]: ``` x**3 + 3*x**2*y + 3*x*y**2 + y**3 ``` In?[14]: ``` 3*x*y**2 + 3*y*x**2 + x**3 + y**3 ``` Out[14]: ``` x**3 + 3*x**2*y + 3*x*y**2 + y**3 ``` 可以通過關鍵詞的形式使用更多的選項: In?[15]: ``` expand(x + y, complex=True) ``` Out[15]: ``` re(x) + re(y) + I*im(x) + I*im(y) ``` In?[16]: ``` I*im(x) + I*im(y) + re(x) + re(y) ``` Out[16]: ``` re(x) + re(y) + I*im(x) + I*im(y) ``` In?[17]: ``` expand(cos(x + y), trig=True) ``` Out[17]: ``` -sin(x)*sin(y) + cos(x)*cos(y) ``` In?[18]: ``` cos(x)*cos(y) - sin(x)*sin(y) ``` Out[18]: ``` -sin(x)*sin(y) + cos(x)*cos(y) ``` ## 3.2.2.2 化簡如果可以將表達式轉化為更簡單的形式，可以使用化簡: In?[19]: ``` simplify((x + x*y) / x) ``` Out[19]: ``` y + 1 ``` 化簡是一個模糊的術語，更準確的詞應該是：powsimp (指數化簡)、 trigsimp (三角表達式)、logcombine、radsimp一起。 **練習** * 計算$(x+y)^6$的展開。 * 化簡三角表達式$ \sin(x) / \cos(x)$ ## 3.2.3 微積分 ### 3.2.3.1 極限在SymPy中使用極限很簡單，允許語法limit(function, variable, point), 因此要計算f(x)類似$x \rightarrow 0$, 你應該使用limit(f, x, 0): In?[5]: ``` limit(sin(x)/x, x, 0) ``` Out[5]: ``` 1 ``` 你也可以計算一下在無限時候的極限: In?[6]: ``` limit(x, x, oo) ``` Out[6]: ``` oo ``` In?[7]: ``` limit(1/x, x, oo) ``` Out[7]: ``` 0 ``` In?[8]: ``` limit(x**x, x, 0) ``` Out[8]: ``` 1 ``` ### 3.2.3.2 微分法你可以使用`diff(func, var)`微分任何SymPy表達式。例如: In?[9]: ``` diff(sin(x), x) ``` Out[9]: ``` cos(x) ``` In?[10]: ``` diff(sin(2*x), x) ``` Out[10]: ``` 2*cos(2*x) ``` In?[11]: ``` diff(tan(x), x) ``` Out[11]: ``` tan(x)**2 + 1 ``` 你可以用下列方法檢查是否正確: In?[12]: ``` limit((tan(x+y) - tan(x))/y, y, 0) ``` Out[12]: ``` tan(x)**2 + 1 ``` 可以用`diff(func, var, n)`方法來計算更高的導數: In?[13]: ``` diff(sin(2*x), x, 1) ``` Out[13]: ``` 2*cos(2*x) ``` In?[14]: ``` diff(sin(2*x), x, 2) ``` Out[14]: ``` -4*sin(2*x) ``` In?[15]: ``` diff(sin(2*x), x, 3) ``` Out[15]: ``` -8*cos(2*x) ``` ### 3.2.3.3 序列展開 SymPy也知道如何計算一個表達式在一個點的Taylor序列。使用`series(expr, var)`: In?[16]: ``` series(cos(x), x) ``` Out[16]: ``` 1 - x**2/2 + x**4/24 + O(x**6) ``` In?[17]: ``` series(1/cos(x), x) ``` Out[17]: ``` 1 + x**2/2 + 5*x**4/24 + O(x**6) ``` **練習** 計算$\lim_{x\rightarrow 0} \sin(x)/x$ 計算`log(x)`對于x的導數。 ### 3.2.3.4 積分法 SymPy支持超驗基礎和特殊函數的無限和有限積分，通過`integrate()` 功能, 使用了強大的擴展的Risch-Norman算法和啟發式和模式匹配。你可以積分基本函數: In?[18]: ``` integrate(6*x**5, x) ``` Out[18]: ``` x**6 ``` In?[19]: ``` integrate(sin(x), x) ``` Out[19]: ``` -cos(x) ``` In?[20]: ``` integrate(log(x), x) ``` Out[20]: ``` x*log(x) - x ``` In?[21]: ``` integrate(2*x + sinh(x), x) ``` Out[21]: ``` x**2 + cosh(x) ``` 也可以很簡單的處理特殊函數: In?[22]: ``` integrate(exp(-x**2)*erf(x), x) ``` Out[22]: ``` sqrt(pi)*erf(x)**2/4 ``` 也可以計算一下有限積分: In?[23]: ``` integrate(x**3, (x, -1, 1)) ``` Out[23]: ``` 0 ``` In?[24]: ``` integrate(sin(x), (x, 0, pi/2)) ``` Out[24]: ``` 1 ``` In?[25]: ``` integrate(cos(x), (x, -pi/2, pi/2)) ``` Out[25]: ``` 2 ``` 不標準積分也支持: In?[26]: ``` integrate(exp(-x), (x, 0, oo)) ``` Out[26]: ``` 1 ``` In?[27]: ``` integrate(exp(-x**2), (x, -oo, oo)) ``` Out[27]: ``` sqrt(pi) ``` #### 3.2.3.5 練習 ### 3.2.4 方程求解 SymPy可以求解線性代數方程，一個或多個變量: In?[28]: ``` solve(x**4 - 1, x) ``` Out[28]: ``` [-1, 1, -I, I] ``` 如你所見，第一個參數是假設等于0的表達式。它可以解一個很大的多項式方程，也可以有能力求解多個方程，可以將各自的多個變量作為元組以第二個參數給出: In?[29]: ``` solve([x + 5*y - 2, -3*x + 6*y - 15], [x, y]) ``` Out[29]: ``` {x: -3, y: 1} ``` 也直接求解超越方程（有限的）: In?[30]: ``` solve(exp(x) + 1, x) ``` Out[30]: ``` [I*pi] ``` 多項式方程的另一個應用是`factor`。`factor`將多項式因式分解為可化簡的項，并且可以計算不同域的因式: In?[31]: ``` f = x**4 - 3*x**2 + 1 factor(f) ``` Out[31]: ``` (x**2 - x - 1)*(x**2 + x - 1) ``` In?[32]: ``` factor(f, modulus=5) ``` Out[32]: ``` (x - 2)**2*(x + 2)**2 ``` SymPy也可以解布爾方程，即，判斷一個布爾表達式是否滿足。對于這個情況，我們可以使用`satisfiable`函數: In?[33]: ``` satisfiable(x & y) ``` Out[33]: ``` {x: True, y: True} ``` 這告訴我們`(x & y)`是真，當x和y都是True的時候。如果一個表達式不是True，即它的任何參數值都無法使表達式為真，那么它將返回False: In?[34]: ``` satisfiable(x & ~x) ``` Out[34]: ``` False ``` ### 3.2.4.1 練習 * 求解系統方程$x + y = 2$, $2\cdot x + y = 0$ * 是否存在布爾值，使$(~x | y) & (~y | x)$為真? ### 3.2.5 線性代數 #### 3.2.5.1 矩陣矩陣通過Matrix類的一個實例來創建: In?[35]: ``` from sympy import Matrix Matrix([[1,0], [0,1]]) ``` Out[35]: ``` Matrix([ [1, 0], [0, 1]]) ``` 與NumPy數組不同，你也可以在里面放入符號: In?[36]: ``` x = Symbol('x') y = Symbol('y') A = Matrix([[1,x], [y,1]]) A ``` Out[36]: ``` Matrix([ [1, x], [y, 1]]) ``` In?[37]: ``` A**2 ``` Out[37]: ``` Matrix([ [x*y + 1, 2*x], [ 2*y, x*y + 1]]) ``` ### 3.2.5.2 微分方程 SymPy可以解 (一些) 常規微分。要求解一個微分方程，使用`dsolve`。首先，通過傳遞cls=Function來創建一個未定義的符號函數: In?[38]: ``` f, g = symbols('f g', cls=Function) ``` f 和 g是未定義函數。我們可以調用f(x), 并且它可以代表未知的函數: In?[39]: ``` f(x) ``` Out[39]: ``` f(x) ``` In?[40]: ``` f(x).diff(x, x) + f(x) ``` Out[40]: ``` f(x) + Derivative(f(x), x, x) ``` In?[41]: ``` dsolve(f(x).diff(x, x) + f(x), f(x)) ``` Out[41]: ``` f(x) == C1*sin(x) + C2*cos(x) ``` 關鍵詞參數可以向這個函數傳遞，以便幫助確認是否找到最適合的解決系統。例如，你知道它是獨立的方程，你可以使用關鍵詞hint=’separable’來強制`dsolve`來將它作為獨立方程來求解: In?[42]: ``` dsolve(sin(x)*cos(f(x)) + cos(x)*sin(f(x))*f(x).diff(x), f(x), hint='separable') ``` Out[42]: ``` [f(x) == -asin(sqrt(C1/cos(x)**2 + 1)) + pi, f(x) == asin(sqrt(C1/cos(x)**2 + 1)) + pi, f(x) == -asin(sqrt(C1/cos(x)**2 + 1)), f(x) == asin(sqrt(C1/cos(x)**2 + 1))] ``` **練習** * 求解Bernoulli微分方程 $x \frac{d f(x)}{x} + f(x) - f(x)^2=0$ * 使用hint=’Bernoulli’求解相同的公式。可以觀察到什么?