SymPy 是一個可以進行數學代數運算的 Python 模組,他的發展目標是一個完整的電腦代數系統,而此模組完全都是使用 Python 寫成的,所以不需要依賴其他的函式庫。

準備 SymPy 環境

在 Ubuntu 或 Debian Linux 中若要安裝 SymPy,可以直接使用 apt 安裝:

sudo apt-get install python-sympy

裝完之後,接著就可以執行 Python 了:

python

進入 Python 的互動式模式(interactive mode)之後,就會看到三個大於符號(>>>),這就是 Python 的主要命令提示字元(primary prompt)。

sympy1

如果是只要使用 SymPy 的功能,可以使用 isympy 這個指令,這是專門為了 SymPy 而設計的互動式介面,其實跟 Python 的互動式環境一樣,但是 isympy 會在一開始設定一些基本的變數與環境,讓使用起來更方便。

isympy

藍色的 In 就是輸入指令的地方,輸出則是紅色的 Out,後面中括弧中的數字是自動編上去的行號。

如果不想那麼麻煩安裝 SymPy,也可以直接使用 SymPy Live,直接在網頁上運行 SymPy,其實更方便,開啓網頁後,左下方就可以直接輸入 Python 指令了。

sympy_live

而 SymPy Live 的環境就跟上面的 isympy 差不多,一開始就會把一些基本的變數與環境設定好。這種線上版本的排版比一般終端機還要漂亮,我們在這裡執行相同的指令,與上面的 isympy 比較一下。

isympy2

因為 SymPy Live 是用 MathJax 來顯示數學式子的,所以看起來很漂亮,幾乎跟 LaTeX 的排版差不多。若是沒什麼特別需求,用這個就很棒了,省去安裝的麻煩。

開始使用 SymPy

接著就可以輸入 Python 指令了,在使用之前要先載入 sympy 這個 Python 模組,輸入:

>>> from sympy import *

輸入完這行指令之後,記得按下 Enter 鍵,正常的話應該沒有什麼訊息出現。這個指令是將 sympy 中所有的東西都載入進來,如果一開始執行了這行指令,之後所有的 import 步驟都可以省略,而如果不想一次載入所有的東西也可以,那就按照下面的 import 作法,需要使用時才載入。

接著就可以來使用 SymPy 了,首先來介紹一些基本的使用方式。在 SymPy 中的數值分為三種類型,分別為浮點數(Float)、有理數(Rational)與整數(Integer)。有理數是由兩個整數相除所得到的數值,例如 Rational(1, 2) 則代表 1/2:

>>> from sympy import Rational
>>> a = Rational(1, 2)

>>> a
1/2

>>> a*2
1

>>> Rational(2)**50/Rational(10)**50
1/88817841970012523233890533447265625

在使用浮點數與整數運算時(尤其是除法)要注意所使用的數值是 SymPy 數值或是原生的 Python 數值,若是使用一般原生的 Python 數值,相除的結果會是浮點數,若在 isympy 中則會使用 __future__ 中的 division 進行除法運算。

>>> 1/2
0.5

若是在舊版的 Python 2 中,division 沒有被 import,則結果會被 truncated:

>>> 1/2
0

這兩種狀況都是使用 Python 原生的數值,不是 SymPy 數值,在大部分的狀況下,在使用 SymPy 時都會使用 SymPy 數值,所以其實可以使用這樣的方式簡化程式碼:

>>> R = Rational
>>> R(1, 2)
1/2

>>> R(1)/2
1/2

最後一個指令中 R(1) 會得到一個 Rational 數值,再除以原生的 Python 數值 2,得到的結果也是一個 Rational 數值。

在數學上有一些常數,例如圓周率(pi)與尤拉常數(e),在 SymPy 中也有支援,在運算時他會以符號的方式進行運算,也就是說 1 + pi 不會計算成浮點數 4.14,而會保持 1 + pi 這個形式,若是要算出實際的浮點數,可以使用 evalf() 這個函數。

>>> from sympy import pi, E
>>> pi**2
pi**2

>>> pi.evalf()
3.14159265358979

>>> (pi + E).evalf()
5.85987448204884

數學上的無限大(infinity)在 SymPy 中也有,表示方式為兩個 o(看起來就像是無限大的寫法)。

>>> from sympy import oo
>>> oo
oo

>>> oo + 1
oo

符號(Symbol)與代數運算

SymPy 與其他的 CAS(Computer Algebra System)不同,在 SymPy 中所有的變數都要經過宣告才能使用:

>>> from sympy import Symbol
>>> x = Symbol('x')
>>> y = Symbol('y')

第二行是將 Python 變數 x 指定為 SymPy Symbol。

在 sympy.abc 中有預先定義一些 Symbol(包涵一些以希臘字母命名的 Symbol)。

>>> from sympy.abc import x, theta

若要建立多個 Symbol 可以使用 symbols 與 var 兩個函數,兩個函數作用相同,都可以接受 range notation,但 var 會自動將建立的 Symbol 加入命名空間(namespace)。

>>> from sympy import symbols, var
>>> a, b, c = symbols('a,b,c')
>>> d, e, f = symbols('d:f')
>>> var('g:h')
(g, h)
>>> var('g:2')
(g0, g1)

建立了 Symbol 之後,就可以使用它們進行代數運算了。

>>> x + y + x - y
2*x

>>> (x + y)**2
(x + y)**2

>>> ((x + y)**2).expand()
x**2 + 2*x*y + y**2

Symbol 亦可被替換為其他的 Symbol 或是數值。

>>> ((x + y)**2).subs(x, 1)
(y + 1)**2

>>> ((x + y)**2).subs(x, y)
4*y**2

>>> ((x + y)**2).subs(x, 1 - y)
1

apart(expr, x) 函數可以處理 partial fraction decomposition。

>>> from sympy import apart
>>> from sympy.abc import x, y, z
>>> 1/( (x + 2)*(x + 1) )
       1
---------------
(x + 1)*(x + 2)

>>> apart(1/( (x + 2)*(x + 1) ), x)
    1       1
- ----- + -----
  x + 2   x + 1

>>> (x + 1)/(x - 1)
x + 1
-----
x - 1

>>> apart((x + 1)/(x - 1), x)
      2
1 + -----
    x - 1

together(expr, x) 則與 apart 相反,將所有的項結合在一起。

>>> from sympy import together
>>> together(1/x + 1/y + 1/z)
x*y + x*z + y*z
---------------
     x*y*z

>>> together(apart((x + 1)/(x - 1), x), x)
x + 1
-----
x - 1

>>> together(apart(1/( (x + 2)*(x + 1) ), x), x)
       1
---------------
(x + 1)*(x + 2)

微積分(Calculus)

SymPy 也可以用於微積分的計算,這裡介紹在微積分中常用的運算。

SymPy 提供了ㄧ個 limit(function, variable, point) 函數可以計算函數的極限值。下面的這個指令可以計算 sin(x) 函數在 x 趨近於零的極限值。

>>> from sympy import limit, Symbol, sin, oo
>>> x = Symbol('x')
>>> limit(sin(x)/x, x, 0)
1

亦可計算無限大的極限值:

>>> limit(x, x, oo)
oo

>>> limit(1/x, x, oo)
0

>>> limit(x**x, x, 0)
1

再舉一個更複雜的例子:

>>> limit((1+1/x)**x, x, oo)
E

>>> limit((x**3-1)/(x**2-1), x, 1)
3/2

>>> limit(E**x/(1-x**3), x, 0)
1

這裡的 E 就是數學上的尤拉常數。

diff(func, var) 這個函數可以計算微分(differentiation):

>>> from sympy import diff, Symbol, sin, tan
>>> x = Symbol('x')
>>> diff(sin(x), x)
cos(x)

>>> diff(sin(2*x), x)
2*cos(2*x)

>>> diff(tan(x), x)
   2
tan (x) + 1

我們可以用下面的指令確認這個答案是正確的:

>>> from sympy import limit
>>> from sympy.abc import delta
>>> limit((tan(x + delta) - tan(x))/delta, delta, 0)
   2
tan (x) + 1

更高階的微分可以使用 diff(func, var, n)

>>> diff(sin(2*x), x, 1)
2*cos(2*x)

>>> diff(sin(2*x), x, 2)
-4*sin(2*x)

>>> diff(sin(2*x), x, 3)
-8*cos(2*x)

SymPy 的 Symbol 可以用 .series(var, point, order) 轉成泰勒展開式(Taylor series):

>>> from sympy import Symbol, cos
>>> x = Symbol('x')
>>> cos(x).series(x, 0, 10)
     2    4     6      8
    x    x     x      x      / 10
1 - -- + -- - --- + ----- + Ox  /
    2    24   720   40320
>>> (1/cos(x)).series(x, 0, 10)
     2      4       6        8
    x    5*x    61*x    277*x     / 10
1 + -- + ---- + ----- + ------ + Ox  /
    2     24     720     8064

另外一個例子:

>>> from sympy import Integral, pprint
>>> y = Symbol("y")
>>> e = 1/(x + y)
>>> s = e.series(x, 0, 5)
>>> print(s)
1/y - x/y**2 + x**2/y**3 - x**3/y**4 + x**4/y**5 + O(x**5)
>>> pprint(s)
          2    3    4
1   x    x    x    x     / 5
- - -- + -- - -- + -- + Ox /
y    2    3    4    5
    y    y    y    y

summation(f, (i, a, b)) 這個函數可以計算 f 函數從 ab 的加總,也就是:

\[ \sum_{i=a}^b f(i) \]

summation 函數如果沒辦法計算出總和,就會顯示原本的加總公式。下面是一些範例:

>>> from sympy import summation, oo, symbols, log
>>> i, n, m = symbols('i n m', integer=True)
>>> summation(2*i - 1, (i, 1, n))
 2
n
>>> summation(1/2**i, (i, 0, oo))
2
>>> summation(1/log(n)**n, (n, 2, oo))
  oo
 ___
   `
       -n
  /   log (n)
 /__,
n = 2
>>> summation(i, (i, 0, n), (n, 0, m))      3    2
m    m    m
-- + -- + -
6    2    3
>>>  from sympy.abc import x
>>>  from sympy import factorial
>>>  summation(x**n/factorial(n), (n, 0, oo))
 x
e

在積分方面,SymPy 支援定積分(definite integration)與不定積分(indefinite integration),而積分的計算是使用 integrate() 函數。

>>> from sympy import integrate, erf, exp, sin, log, oo, pi, sinh, symbols
>>> x, y = symbols('x,y')
>>> integrate(6*x**5, x)
 6
x
>>> integrate(sin(x), x)
-cos(x)
>>> integrate(log(x), x)
x*log(x) - x
>>> integrate(2*x + sinh(x), x)
 2
x  + cosh(x)

更複雜的例子:

>>> integrate(exp(-x**2)*erf(x), x)
  ____    2
/ pi *erf (x)
--------------
      4

計算定積分:

>>> integrate(x**3, (x, -1, 1))
0
>>> integrate(sin(x), (x, 0, pi/2))
1
>>> integrate(cos(x), (x, -pi/2, pi/2))
2

瑕積分(improper integral)也同樣支援:

>>> integrate(exp(-x), (x, 0, oo))
1
>>> integrate(log(x), (x, 0, 1))
-1

複數(Complex numbers)

SymPy 也可以處理複數,Symbol 亦可以指定屬性(attribute),例如:real、positive 或 complex 等,指定屬性會影響其運算的結果。

>>> from sympy import Symbol, exp, I
>>> x = Symbol("x") # a plain x with no attributes
>>> exp(I*x).expand()
 I*x
e
>>> exp(I*x).expand(complex=True)
   -im(x)               -im(x)
I*e      *sin(re(x)) + e      *cos(re(x))
>>> x = Symbol("x", real=True)
>>> exp(I*x).expand(complex=True)
I*sin(x) + cos(x)

三角函數(Trigonometric)

>>> from sympy import asin, asinh, cos, sin, sinh, symbols, I
>>> x, y = symbols('x,y')
>>> sin(x + y).expand(trig=True)
sin(x)*cos(y) + sin(y)*cos(x)
>>> cos(x + y).expand(trig=True)
-sin(x)*sin(y) + cos(x)*cos(y)
>>> sin(I*x)
I*sinh(x)
>>> sinh(I*x)
I*sin(x)
>>> asinh(I)
I*pi
----
 2
>>> asinh(I*x)
I*asin(x)
>>> sin(x).series(x, 0, 10)
     3     5     7       9
    x     x     x       x       / 10
x - -- + --- - ---- + ------ + Ox  /
    6    120   5040   362880

>>> sinh(x).series(x, 0, 10)
     3     5     7       9
    x     x     x       x       / 10
x + -- + --- + ---- + ------ + Ox  /
    6    120   5040   362880

>>> asin(x).series(x, 0, 10)
     3      5      7       9
    x    3*x    5*x    35*x     / 10
x + -- + ---- + ---- + ----- + Ox  /
    6     40    112     1152

>>> asinh(x).series(x, 0, 10)
     3      5      7       9
    x    3*x    5*x    35*x     / 10
x - -- + ---- - ---- + ----- + Ox  /
    6     40    112     1152

Spherical Harmonics

>>> from sympy import Ylm
>>> from sympy.abc import theta, phi
>>> Ylm(1, 0, theta, phi)
  ___
/ 3 *cos(theta)
----------------
        ____
    2*/ pi

>>> Ylm(1, 1, theta, phi)
   ___  I*phi
-/ 6 *e     *sin(theta)
------------------------
            ____
        4*/ pi

>>> Ylm(2, 1, theta, phi)
   ____  I*phi
-/ 30 *e     *sin(theta)*cos(theta)
------------------------------------
                  ____
              4*/ pi

階乘(factorials)與 Gamma 函數

>>> from sympy import factorial, gamma, Symbol
>>> x = Symbol("x")
>>> n = Symbol("n", integer=True)

>>> factorial(x)
x!

>>> factorial(n)
n!

>>> gamma(x + 1).series(x, 0, 3) # i.e. factorial(x)
                      /          2     2
                    2 |EulerGamma    pi |    / 3
1 - EulerGamma*x + x *|----------- + ---| + Ox /
                           2         12/

Zeta 函數

>>> from sympy import zeta
>>> zeta(4, x)
zeta(4, x)

>>> zeta(4, 1)
  4
pi
---
 90

>>> zeta(4, 2)
       4
     pi
-1 + ---
      90

>>> zeta(4, 3)
         4
  17   pi
- -- + ---
  16    90

多項式(polynomials)

>>> from sympy import assoc_legendre, chebyshevt, legendre, hermite
>>> chebyshevt(2, x)
   2
2*x  - 1

>>> chebyshevt(4, x)
   4      2
8*x  - 8*x  + 1

>>> legendre(2, x)
   2
3*x    1
---- - -
 2     2

>>> legendre(8, x)
      8         6         4        2
6435*x    3003*x    3465*x    315*x     35
------- - ------- + ------- - ------ + ---
  128        32        64       32     128

>>> assoc_legendre(2, 1, x)
        __________
       /    2
-3*x*/  - x  + 1

>>> assoc_legendre(2, 2, x)
     2
- 3*x  + 3

>>> hermite(3, x)
   3
8*x  - 12*x

微分方程(Differential Equations)

>>> from sympy import Function, Symbol, dsolve
>>> f = Function('f')
>>> x = Symbol('x')
>>> f(x).diff(x, x) + f(x)
        2
       d
f(x) + ---(f(x))
         2
       dx

>>> dsolve(f(x).diff(x, x) + f(x), f(x))
f(x) = C1*sin(x) + C2*cos(x)

代數方程式(Algebraic Equations)

>>> from sympy import solve, symbols
>>> x, y = symbols('x,y')
>>> solve(x**4 - 1, x)
[-1, 1, -I, I]

>>> solve([x + 5*y - 2, -3*x + 6*y - 15], [x, y])
{x: -3, y: 1}

線性代數(Linear Algebra)

建立矩陣可以使用 Matrix 類別:

>>> from sympy import Matrix, Symbol
>>> Matrix([[1, 0], [0, 1]])
[1  0]
[    ]
[0  1]

矩陣中亦可包含 Symbol:

>>> x = Symbol('x')
>>> y = Symbol('y')
>>> A = Matrix([[1, x], [y, 1]])
>>> A
[1  x]
[    ]
[y  1]

>>> A**2
[x*y + 1    2*x  ]
[                ]
[  2*y    x*y + 1]

Pattern Matching

使用 .match() 方法並配合 Wild 類別可以對 Symbol 進行比對,例如抓取某一項的係數等,這個方法會傳回一個 dictionary,例如:

>>> from sympy import Symbol, Wild
>>> x = Symbol('x')
>>> p = Wild('p')
>>> (5*x**2).match(p*x**2)
{p: 5}

>>> q = Wild('q')
>>> (x**2).match(p*x**q)
{p: 1, q: 2}

如果沒有找到指定的 pattern,則會傳回 None:

>>> print (x + 1).match(p**x)
None

亦可使用 Wild 的 exclude 參數排除指定的比對結果:

>>> p = Wild('p', exclude=[1, x])
>>> print (x + 1).match(x + p) # 1 is excluded
None
>>> print (x + 1).match(p + 1) # x is excluded
None
>>> print (x + 1).match(x + 2 + p) # -1 is not excluded
{p_: -1}

輸出顯示(Printing)

SymPy 的運算結果可以用好幾種不同的方式輸出,一般預設的方式是使用 str(expression) 的方式,例如:

>>> from sympy import Integral
>>> from sympy.abc import x
>>> print x**2
x**2
>>> print 1/x
1/x
>>> print Integral(x**2, x)
Integral(x**2, x)

pprint() 函數(pretty printing)可以將結果以 ascii-art 的方式「畫」出來,例如:

>>> from sympy import Integral, pprint
>>> from sympy.abc import x
>>> pprint(x**2)
 2
x
>>> pprint(1/x)
1
-
x
>>> pprint(Integral(x**2, x))
  /
 |
 |  2
 | x  dx
 |
/

若您的系統有安裝 unicode 字型,則 pprint() 函數預設就會使用它,亦可使用 use_unicode 參數自行設定。

>>> pprint(Integral(x**2, x), use_unicode=True)

  2
 x  dx

下面示範將 pprint() 設定為預設的輸出函數:

$ python
Python 2.5.2 (r252:60911, Jun 25 2008, 17:58:32)
[GCC 4.3.1] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> from sympy import init_printing, var, Integral
>>> init_printing(use_unicode=False, wrap_line=False, no_global=True)
>>> var("x")
x
>>> x**3/3
 3
x
--
3
>>> Integral(x**2, x) #doctest: +NORMALIZE_WHITESPACE
  /
 |
 |  2
 | x  dx
 |
/

使用 Python 輸出:

>>> from sympy.printing.python import python
>>> from sympy import Integral
>>> from sympy.abc import x
>>> print python(x**2)
x = Symbol('x')
e = x**2
>>> print python(1/x)
x = Symbol('x')
e = 1/x
>>> print python(Integral(x**2, x))
x = Symbol('x')
e = Integral(x**2, x)

使用 LaTeX 輸出:

>>> from sympy import Integral, latex
>>> from sympy.abc import x
>>> latex(x**2)
x^{2}
>>> latex(x**2, mode='inline')
$x^{2}$
>>> latex(x**2, mode='equation')
begin{equation}x^{2}end{equation}
>>> latex(x**2, mode='equation*')
begin{equation*}x^{2}end{equation*}
>>> latex(1/x)
frac{1}{x}
>>> latex(Integral(x**2, x))
int x^{2}, dx

使用 MathML 輸出:

>>> from sympy.printing.mathml import mathml
>>> from sympy import Integral, latex
>>> from sympy.abc import x
>>> print mathml(x**2)
<apply><power><ci>x</ci><cn>2</cn></power></apply>
>>> print mathml(1/x)
<apply><power><ci>x</ci><cn>-1</cn></power></apply>

使用 Pyglet:

>>> from sympy import Integral, preview
>>> from sympy.abc import x
>>> preview(Integral(x**2, x))

如果系統上有安裝 pyglet 的話,就會跳出這個視窗:

pngview1

參考資料