數值積分

Romberg Integration

這一節的重點與 Richardson's extrapolation 類似,就是我們可以透過已知的數值積分 \(I(a,b,h)\) 和 \(I(a,b,\frac{h}{2})\) 進一步提昇 truncation error 的階數。首先我們以 trapezoidal rule 為例,我們以 \(I_{T,2}(a,b,\frac{h}{2})\) 符號代表2 階的 trapezoidal rule 計算的積分結果,下標的 T 代表 trapezoidal rule,2 代表 2 階。

(RI:式1):

\[I_{T,2}(a,b,\frac{h}{2}) = \frac{h}{2} \left[ \frac{f(a)+f(b)}{2} + \sum_{k=1}^{2N-1}f(x_{k/2}) \right] \\= \frac{h}{2} \left[ \frac{f(a)+f(b)}{2} + \sum_{k=1}^{N-1} f(x_{2k/2}) + \sum_{k=0}^{N-1}f(x_{(2k+1)/2}) \right] \\= \frac{1}{2} \left[ I_{T,2}(a,b,h)  + h\sum_{k=0}^{N-1}f(x_{(2k+1)/2}) \right]\]

透過這個算法,我們從已計算的 \(I_{T,2}(a,b,h)\) 結果利用上式的最後一項去計算 \(I_{T,2}(a,b,\frac{h}{2}) \) 的結果,避免一些計算上的浪費。當我們有了 \(I_{T,2}(a,b,h)\) 和 \(I_{T,2}(a,b,\frac{h}{2}) \) 的值,我們知道他們都是 \(O(h^2)\) 的誤差,配合類似Richarson's extrapolation 的技巧,我們可以把 truncation error 的階數提昇。所以更高階的數值積分式為
\[ \frac{2^{2}I_{T,2}(a,b,h)-I_{T,2}(a,b,2h)}{2^2-1}\\
=\frac{1}{3}\left[2h\left(f(a)+f(b)+2\sum_{k=1}^{N-1} f(x_k))-h(f(a)+f(b)+2\sum_{m=1}^{N/2-1} f(x_{2m})\right)\right] \\

=\frac{h}{3}\left(f(a)+f(b)+4\sum_{m=1}^{N/2}f(x_{2m-1})+2\sum_{m=1}^{N/2-1}f(x_{2m})\right) \]

這個式子剛好是Simpson's rule,因此提昇後數值積分不是3 階,而是4 階。由於Simpson's rule 的高階誤差表示剛好是4, 6, 8 這種偶數的間距,所以我們從2 階的trapezoidal rule 出發,得到下列的提升公式

(RI:式2):

\[I_{T,2(k+1)}(a,b,h) = \frac{2^{2k}I_{T,2k}(a,b,h)-I_{T,2k}(a,b,2h)}{2^{2k}-1}\]

和先前的數值微分一樣,知道truncation error 的階數,只是告訴我們當 \(h\) 變化時可以預期它往真值靠近的收斂速率,要知道計算當時的真正誤差,我們需要知道 truncation error term 的領導係數。要估計這個領導係數並且估計與真值的誤差,我們也是透過比較兩個不同 step size 的 \(h\),因此我們得到下列的誤差估計式

(RI:式3):
\[|E_{T,2(k+1)} (2^{-k}h)| \cong \frac{|I_{T,2k}(2^{-k}h)-I_{T,2k}(2^{-(k-1)}h)|}{2^{2k}-1}\]

這個提昇的過程稱為Romberg integration,我們可以從下列這個圖表看到,往下是 \(h\) 的step size 縮小,我們用(RI:式1)來有效率計算更小step size 的值,往右是把truncation error 的階數增加,我們用(RI:式2)來計算。

Table of Romberg's Integration

以下是Python 的Romberg integration table 程式:

import numpy as np
def romberg(f,a,b,tol = 10**(-10),K=10):
    """
    Construct Romberg table to find the definite integral of f over [a,b]
    
    f : the input function
    [a,b] : the interval
    tol : the tolerance(default is 10**(-10))
    K : the maximal step of iteration

    For example:
    (S,R,err) = romberg(lambda x: 2*x**3+3*x+2,0,1,10**(-8),15)
    
    S : the numerical integration of f on [a,b]
    R : the list of Romberg table
    err : the estimatation error
    
    If the rusult reach the tolerance, then we will break the loop 
    before K iteration.
    """
    h = float(b-a)
    N = 1
    R = [[float(f(a)+f(b))*h/2]]
    for k in xrange(2,K+1):
        h /= 2
        N *= 2
        x = np.array([a+i*h for i in xrange(1,N,2)])
        Rt = [R[-1][-1]/2 + h*sum(f(x))]
        tmp = 1.0
        for n in xrange(k-1):
              tmp *= 4
              Rt.append((tmp*Rt[n]-R[-1][n])/(tmp-1))

        R.append(Rt)
        err = abs(R[-1][-2]-R[-2][-1])/(tmp-1)
        if err < tol:
            break

    S = R[-1][-1]
    return S,R,err

 

上列的程式會計算 \(f(x)\) 在 \([a,b]\) 區間的積分,tol 是設定我們期待的計算誤差,若不給值預設為10 的-10次方,\(K\) 設定為我們希望計算多少次提昇的步驟,預設為10次。若提早達到計算誤差,該程式會提前結束。