數值積分

Adaptive Quadrature


從Romberg Integration 的程式中我們可以看到,給定 \([a,b]\) 和 \(f(x)\),我們在不知道真解的情況下可以透過兩個不同的step size 的計算結果,來推估誤差。從泰勒展開式的分析中,我們知道函數 \(f(x)\) 的平滑度是影響truncation error 的主因。因此,不同平滑度的函數,要達到同樣的精準度需要的step size 也不同。若一個函數在不同區段上的平滑度不同,我們可以透過在平滑的區域使用較大的 \(h\) ,在變化率較大的區域使用較小的 \(h\) 達到使用最少的計算量卻得到好的精準度的計算結果,這樣的數值積分方法稱為Adaptive Quadrature。要做到這件事其實不難,從Romberg Integration 的公式(RI:式3)可以看到,程式有能力自行估計在該區間上的誤差。因此,當誤差達到預期,就停止計算,若誤差過大,就把 \(h\) 變小,進行更細緻的計算。

所以,Adaptive Quadrature 的概念如下,假設我們使用Simpson's rule 來計算數值積分,給定一個區間 \([a,b]\),我們會得到一個只用三個函數值所計算出來的積分值,\(h = b-1\)。然後我們取\(a,b \) 的中點 \(c\),把區間切割為 \([a,c]\) 和 \([c,b]\)。並且分別計算這兩個積分所加總的積分值與先前 \([a,b]\) 上的積分值的差距。透過(EA:式一)的公式我們知道數值解和真解的差距為這兩次計算的差距除上15。因此當這個量小於我們可以忍受的誤差時,我們停止計算程序;若不是我們就繼續在每一個小區段上進行重複的切割,再計算數值積分,直到每一個區段上計算出來的誤差都小於我們可以忍受的誤差。以Simpson's rule 計算的Adaptive Quadrature 程式如下,若要換成別的計算方式,對應的誤差估計式也需要同步更換,例如,以trapezoidal rule 計算時我們用(EA:式二)。

def adapt_simpson(f,a,b,tol=10**(-9)):
    """
    The adaptive quadrature of f(x) on [a,b] interval with the tolerent 
    tol by Simpson's rule

    f : function of f(x)
    [a,b] : the interval
    tol : tolerance

    For example:
    S, err, nodes = adapt_simpson(lambda x: 3*x**2+2*x+1,0,2,10**(-12))

    S : the numerical integration result
    nodes : the nodes for computation
    err : the estimated error
    """
    S = simpson(f,a,b,1)
    S, err, nodes = recursive_simpson(f,a,b,S,tol)

    return S, err, nodes

def recursive_simpson(f,a,b,S,tol):
    """
    Recursive implementation of adaptive Simpson's rule

    f: function of f(x)
    [a,b] : the interval of integration
    S : the previous integration result
    tol : the tolerance
    
    This is a subfunction of adapt_simpson.
    """
    c = float(a+b)/2
    SL = simpson(f,a,c,1)
    SR = simpson(f,c,b,1)
    Sn = SL + SR
    err = abs(Sn-S)/15.0
    if err <= tol:
        S = SL+SR +err
        nodes = [a,c,b]
        return S, err, nodes
    else:
        SL, err1, nodes1 = recursive_simpson(f,a,c,SL,tol/2.0)
        SR, err2, nodes2 = recursive_simpson(f,c,b,SR,tol/2.0)
        err = err1 + err2
        nodes = nodes1[0:-1]
        nodes.extend(nodes2)
    
    return SL+SR, err, nodes

 

adapt_simpson 是執行的主程式,recursive_simpson 是遞迴的計算式,使用方法如下: 

>>> S, err, nodes = adapt_simpson(lambda x: np.sin(x),0,2,10**(-12))
>>> S,err
(1.4161468365475729, 2.1530912903272581e-13)
>>> len(nodes)
461
>>> from pylab import *
>>> x = np.linspace(0,2,100)
>>> y = np.sin(x)
>>> x1 = np.array(nodes)
>>> y1 = np.sin(x1)
>>> plot(x,y,'r')
[<matplotlib.lines.Line2D object at 0x1d341d0>]
>>> plot(x1,y1,'b*')
[<matplotlib.lines.Line2D object at 0x60bd290>]
>>> fig1 = figure(1)
>>> fig1.show()
>>> fig1.savefig('adap_quadrature_sin.png')

這裡我們積\([0,2]\) 區間的 \(sin()\) 函數,顯示 S, err 得到積分值為1.4161468365475729 並且估計誤差小於 \(10^{-12}\)。nodes 是 \([0,2]\) 區間中彈性調整使用到的取值點,共有461 個。接下來的指令是把函數以及取點的位置畫在同一張圖上,方便讀者觀察取點位置與函數平滑度的關係。