數值積分

基礎數值積分

對於一個從 \(R^1 \to R^1\) 的函數,在 \([a,b]\) 區間上的定積分,就是把 \(a,b\) 線段上切割成許多細小的格點 \(x_k\),我們很容易想像將每一個小區間 \([x_k,x_{k+1}]\) 上的積分值加總起來就是 \(\int_a^b f(x) dx\)。當函數 \(f (x)\) 在 \([x_k,x_{k+1}]\) 上夠平滑時,我們可以用 \(f(x_{mk})(x_{k+1}-x_k)\) 來代替這個小區間的積分,其中 \( x_{mk} = \frac{x_k+x_{k+1}}{2}\)這個方法稱為中點法。我們也可以用這個區段頭尾相接的梯形面積來代表,稱為平行四邊形法。大一微積分裡我們通常會教到下列方法:

midpoint rule

\(\int_{x_k}^{x_{k+1}} f(x) dx \simeq h f_{mk}\),其中 \(h = x_{k+1}-x_k\), \(f_{mk} = f(x_{mk}) \)

trapezoidal rule

\(\int_{x_k}^{x_{k+1}} f(x) dx \simeq \frac{h}{2}(f_k+f_{k+1})\),其中 \(f_k = f(x_k) \)

Simpson's rule

\(\int_{x_{k-1}}^{x_{k+1}} f(x) dx \simeq \frac{h}{3}(f_{k-1}+4f_k+f_{k+1})\),其中 \(h = \frac{x_{k+1}-x_{k-1}}{2}\)


要理解 midpoint rule 和 trapezoidal rule 很容易從圖形中想像並理解如何得到上述的公式,但是 Simpson's rule 的係數,就沒有那麼直覺。我們可以這樣看,在 midpoint rule 中我們用一個常數函數來取代區間 \([x_k,x_{k+1}]\) 裡的 \(f(x)\),在 trapezoidal rule 中,我們用一個線型函數來取代區間 \([x_k,x_{k+1}]\) 裡的 \(f(x)\),而 Simpson's rule 呢?直覺上猜想是用二次函數來取代,意思是我們要找一個二次多項式 \(p_{(2)}(x)\) 通過 \((x_{k-1},f_{k-1})\), \((x_k,f_k)\) 和 \((x_{k+1},f_{k+1})\) 三個點。然後用這個多項式在\([x_{k-1},x_{k+1}]\) 上的積分來逼近 \(f (x)\) 的積分。經過這樣的計算我們便能得到 Simpson's rule 的對應係數。

若我們能逼近每一個 \([x_k,x_{k+1}]\) 區間上的積分值,將這些區間上的值加總起來就是在 \([a,b]\) 區間上函數 \(f(x)\) 積分的逼近值。

下列兩個 Python 程式分別對應上述的 trapezoidal rule 數值積分和 Simpson's rule 數值積分。

def trapz(f,a,b,N):

    """
    Approximate the definite integral of f(x) from a to b by the 
    composite trapezoidal rule, using N subintervals

    For example:
    trapz(lambda x: x**2+2*x+1,0.0,3.0,100)
    """
    return (b-a)*(float(f(a)+f(b))/2+sum([float(f(a+(b-a)*k/N)) for k in xrange(1,N)]))/N

 

因為頭尾相接的緣故,我們只有最前項和最後項除上2,中間項直接相加。程式中特別加註float() 函式的部份,是避免程式使用到整數乘法,造成計算上的錯誤。

def simpson(f, a, b, N):
    """
    Approximate the definite integral of f(x) from a to b by the
    Simpson's rule, using N subintervals, N must be even integer.

    For example:
    simpson(lambda x: x**2+2*x+1,0.0,3.0,100)
    """
    if N%2==0:
        n = N
    else:
        n = N+1
        
    h = float(b - a) / n;
    S = float(f(a))
 
    for i in range(1, n, 2):
        x = a + h * i
        S += 4 * f(x)
 
    for i in range(2, n-1, 2):
        x = a + h * i
        S += 2 * f(x)
 
    S += f(b)
    S = h*S/3
 
    return S

 

上列是Simpson 法的Python 程式,使用的方法很簡單,如果你的函數不複雜,可以直接用說明範例的方式來計算積分值。若函數較為複雜,我們用以下範例來執行。

>>> def myfun(x):
.... if x < 0:
.... return 1-x
.... else:
.... return x**2
>>> simpson(lambda x: myfun(x),-1,2,1000)
4.1676666679999999