<dd id="sga4y"><optgroup id="sga4y"></optgroup></dd>
<nav id="sga4y"><code id="sga4y"></code></nav>
  • <optgroup id="sga4y"></optgroup>

    SymPy-符號運算好幫手?

    SymPy是Python的數學符號計算庫,用它可以進行數學公式的符號推導。為了調用方便,下面所有的實例程序都假設事先從sympy庫導入了所有內容:

    >>> from sympy import *
    

    封面上的經典公式?

    本書的封面上的公式:

    e^{\mathrm{i} \pi} + 1 = 0

    叫做歐拉恒等式,其中e是自然指數的底,i是虛數單位, \pi 是圓周率。此公式被譽為數學最奇妙的公式,它將5個基本數學常數用加法、乘法和冪運算聯系起來。下面用SymPy驗證一下這個公式。

    載入的符號中,E表示自然指數的底,I表示虛數單位,pi表示圓周率,因此上述的公式可以直接如下計算:

    >>> E**(I*pi)+1
    0
    

    歐拉恒等式可以下面的公式進行計算,

    e^{{\mathrm{i}}x}=\cos x+ {\mathrm{i}} \sin x

    為了用SymPy求證上面的公式,我們需要引入變量x。在SymPy中,數學符號是Symbol類的對象,因此必須先創建之后才能使用:

    >>> x = Symbol('x')
    

    expand函數可以將公式展開,我們用它來展開E**(I*pi)試試看:

    >>> expand( E**(I*x) )
    exp(I*x)
    

    沒有成功,只是換了一種寫法而已。這里的exp不是math.exp或者numpy.exp,而是sympy.exp,它是一個類,用來表述自然指數函數。

    expand函數有關鍵字參數complex,當它為True時,expand將把公式分為實數和虛數兩個部分:

    >>> expand(exp(I*x), complex=True)
    I*exp(-im(x))*sin(re(x)) + cos(re(x))*exp(-im(x))
    

    這次得到的結果相當復雜,其中sin, cos, re, im都是sympy定義的類,re表示取實數部分,im表示取虛數部分。顯然這里的運算將符號x當作復數了。為了指定符號x必須是實數,我們需要如下重新定義符號x:

    >>> x = Symbol("x", real=True)
    >>> expand(exp(I*x), complex=True)
    I*sin(x) + cos(x)
    

    終于得到了我們需要的公式。那么如何證明它呢。我們可以用泰勒多項式展開:

    >>> tmp = series(exp(I*x), x, 0, 10)
    >>> pprint(tmp)
               2      3    4      5     6      7      8        9
              x    I*x    x    I*x     x    I*x      x      I*x
    1 + I*x - -- - ---- + -- + ---- - --- - ---- + ----- + ------ + O(x**10)
              2     6     24   120    720   5040   40320   362880
    

    series是泰勒展開函數,pprint將公式用更好看的格式打印出來。下面分別獲得tmp的實部和虛部,分別和cos(x)和sin(x)的展開公式進行比較:

    >>> pprint(re(tmp))
                        2    4     6      8
                       x    x     x      x
    1 + re(O(x**10)) - -- + -- - --- + -----
                       2    24   720   40320
    
    >>> pprint( series( cos(x), x, 0, 10) )
         2    4     6      8
        x    x     x      x
    1 - -- + -- - --- + ----- + O(x**10)
        2    24   720   40320
    
    >>> pprint(im(tmp))
                        3     5     7       9
                       x     x     x       x
    x + im(O(x**10)) - -- + --- - ---- + ------
                       6    120   5040   362880
    
    >>> pprint(series(sin(x), x, 0, 10))
         3     5     7       9
        x     x     x       x
    x - -- + --- - ---- + ------ + O(x**10)
        6    120   5040   362880
    

    球體體積?

    用SciPy數值積分一節我們介紹了如何使用數值定積分計算球體的體積,而SymPy的符號積分函數integrate則可以幫助我們進行符號積分。integrate可以進行不定積分:

    >>> integrate(x*sin(x), x)
    -x*cos(x) + sin(x)
    

    如果指定x的取值范圍的話,integrate則進行定積分運算:

    >>> integrate(x*sin(x), (x, 0, 2*pi))
    -2*pi
    

    為了計算球體體積,首先讓我們來看看如何計算圓形面積,假設圓形的半徑為r,則圓上任意一點的Y坐標函數為:

    y(x) = \sqrt{r^2 - x^2}

    因此我們可以直接對上述函數在-r到r區間上進行積分得到半圓面積,注意這里我們使用symbols函數一次創建多個符號:

    >>> x, y, r = symbols('x,y,r')
    >>> 2 * integrate(sqrt(r*r-x**2), (x, -r, r))
    2*Integral((r**2 - x**2)**(1/2), (x, -r, r))
    

    很遺憾,integrate函數沒有計算出結果,而是直接返回了我們輸入的算式。這是因為SymPy不知道r是大于0的,如下重新定義r,就可以得到正確答案了:

    >>> r = symbols('r', positive=True)
    >>> circle_area = 2 * integrate(sqrt(r**2-x**2), (x, -r, r))
    >>> circle_area
    pi*r**2
    

    接下來對此面積公式進行定積分,就可以得到球體的體積,但是隨著X軸坐標的變化,對應的切面的的半徑會發生變化,現在假設X軸的坐標為x,球體的半徑為r,則x處的切面的半徑為可以使用前面的公式y(x)計算出。

    _images/mayavi2_sphere.png

    球體體積的雙重定積分示意圖

    因此我們需要對circle_area中的變量r進行替代:

    >>> circle_area = circle_area.subs(r, sqrt(r**2-x**2))
    >>> circle_area
    pi*(r**2 - x**2)
    

    用subs進行算式替換

    subs函數可以將算式中的符號進行替換,它有3種調用方式:

    • expression.subs(x, y) : 將算式中的x替換成y
    • expression.subs({x:y,u:v}) : 使用字典進行多次替換
    • expression.subs([(x,y),(u,v)]) : 使用列表進行多次替換

    請注意多次替換是順序執行的,因此:

    expression.sub([(x,y),(y,x)])
    

    并不能對兩個符號x,y進行交換。

    然后對circle_area中的變量x在區間-r到r上進行定積分,得到球體的體積公式:

    >>> integrate(circle_area, (x, -r, r))
    4*pi*r**3/3
    

    首頁目錄

    上一篇文章

    SciPy-數值計算庫

    下一篇文章

    matplotlib-繪制精美的圖表

    97av