數值積分

Gauss Quadrature

我們在介紹基礎midpoint rule 和trapezoidal rule 時我們發現前者是用常數函數來逼近真實函數,後者是用一次函數來逼近真實函數,話說使用越高階的函數來逼近真實函數會得到較好的計算結果,然而這兩種算法所對應的數值公式卻是完全相同的(當然,對應的數值解也相同)。因此,我們有一個結果就是甚選計算用的點,可以增加計算的truncation error次數。

Gauss Quadrature 有下列幾種主要形式:

  1. \(\int_{a}^{b} f(x)dx\)
  2. \(\int_{-\infty}^{\infty} f(x) dx\)
  3. \(\int_{0}^{\infty} e^{-x} f(x) dx\)
  4. \(\int_{-1}^1 \frac{1}{ \sqrt{1-x^2} } f(x) dx \)
  5. \(\int_{-1}^{1} \sqrt{1-x^{2}} f(x)dx\)

無論是上列何種形式,我們的目標是選取最適合的\(w_i\) 和最適合的 \(x_i\) 使得 \(\sum w_{i} f(x_{i})\) 最靠近我們的積分值。

Gauss-Legendre Integration

假設 \(f(x)\) 是一個 3 次以下的多項式,並且我們試圖用兩個點來計算 \(f(x)\) 在 \([-1,1]\) 上的積分。

  1. 若 \(f(x) = 1\) : \(w_1 f(x_1) + w_2 f(x_2) =\int_{-1}^1 f dx = 2\)
  2. 若 \(f(x) = x\) : \(w_1 f(x_1) + w_2 f(x_2) = w_1 x_1 + w_2 x_2 = \int_{-1}^1 x dx = 0\)
  3. 若 \(f(x) = x^2\) : \(w_1 f(x_1) + w_2 f(x_2) = w_1 x_1^2 + w_2 x_2^2 = \int_{-1}^1 x^2 dx = \frac{2}{3} \)
  4. 若\(f(x) = x^3\) : \(w_1 f(x_1) + w_2 f(x_2) = w_1 x_1^3 + w_2 x_2^3 = \int_{-1}^1 x^3 dx = 0 \)

上列我們有四個未知數,四個方程式,經過求解後得到 \(w_1 = w_2 = 1\), \(x_1 = -\frac{1}{\sqrt{3}}\) 以及 \(x_2 = \frac{1}{\sqrt{3}}\)。當我們使用 \(I[-1,1] = f(-\frac{1}{\sqrt{3}}) + f(\frac{1}{\sqrt{3}})\) 時,對應的函數 \(f(x)\) 在3 階以下,我們都可以沒有誤差的計算出來。

我們可以仿造上面的方式,增加取值的點數 \(N\),因為每一個取值的點會對應一個權重,所以我們會有 \(2N\) 個變數,當我們使用 \(2N\) 個方程式時,對應的多項式階數是 \(2N-1\)。寫成矩陣表示 \(Aw=b\) 是一個由 \(x_{j}^{i}\) 構成的 \(2N\) -by- \(N\) 的矩陣\(A\),其中 \(A[i,j] = x_{j}^{i}\),\(w[j] = w_{j}\) 以及 \(b[j]=\frac{1-(-1)^{j}}{j}\),其中 \(i = 0,\dots,2N-1\),\(j = 0,\dots,N-1\)。\(A\) 矩陣的特性是只要 \(x_{j}\) 點沒有重複,\(rank(A) = N\)。意思是 \(A^{T}A\) 會是一個可逆矩陣,當我們知道每一個 \(x_{j}\) 的值時,我們可以求得 \(w\)。問題是,如何讓 \(b\) 落在 \(A\) 的 column space ,避免我們算出的是最小平方法的解呢?答案是當 \(x_{j}\) 是 Legendre polynomial 的根時,\(b\) 會落在 \(A\) 的column space 裡。有興趣深入的同學,請參考Math Word 關於Legendre-Gauss Quadrature 的介紹。Legendre polynomial 的公式為
\[L_{N}(x) = \sum_{i=0}^{floor(N/2)} (-1)^{i}\frac{(2N-2i)!x^{N-2i}}{2^{N}i!(N-i)!(N-2i)!}\]
其中floor() 函數是往下找最靠近的整數。對 Legendre polynomial 的認識是,他是一系列在 \([-1,1]\) 區間上的多項式正交基底。我們可以用
\[L_{N}(x) = \frac{1}{N}((2N-1)xL_{N-1}(x) - (N-1)L_{N-2}(x))\]
疊代公式來計算Legendre polynomial。

def Legendre_poly(N):
    """
    The generator of the coefficient of Legendre polynomial of degree N

    For example:
    p = Legendre_poly(3)
    """
    if N <=0:
        c = np.array([1])
    elif N==1:
        c = np.array([1,0])
    else:
        a = Legendre_poly(N-1)
        a = list(a)
        a.append(0)
        a = np.array(a)

        c = list(Legendre_poly(N-2))
        b = [0,0]
        b.extend(c)
        b = np.array(b)

        c = ((2*N-1)*a+(N-1)*b)/float(N)

    return c

def Gauss_Legendre(N):
    """
    The generator of the points and the weightings of Gauss-Lagendre Integration of order N

    For exammple:
    x,w = Gauss_Legendre(3)
    """
    if N<0:
        print "Gauss-Legendre polynomial of negative order??"
        return
    
    c = Legendre_poly(N)
    p = np.poly1d(c)
    x = p.r
    
    A = np.matrix(np.zeros([N,N]),dtype='complex')
    b = np.matrix(np.zeros([N,1]),dtype='complex')
    A[0,:] = 1
    b[0,0] = 2
    for i in xrange(1,N):
        for j in xrange(N):
            A[i,j] = x[j]**i
        if i%2 <>0:
            b[i,0] = 0
        else:
            b[i,0] = float(2)/(i+1)

    w = la.solve(A,b)
                  
    return x,w

def Gauss_Legendre_quadrature(f,a,b,N):
    """
    The Gauss-Legendre Integration of f(x) over [a,b] with N grid points.

    For example:
    S = Gauss_Legendre_quadrature(lambda x: 3*x**5+4*x**3-2*x+1,-2,3,10)
    """

    x,w = Gauss_Legendre(N)
    x = ((b-a)*x+a+b)/2
    y = f(x)

    S = np.dot(y,w)*(b-a)/2.0
    return float(S.real)

 

上述的程式分三部分,第一部分是Legendre_poly(),這個函式會計算出對應N 階Legendre polynomial 的係數,第二部分是Gauss_Legendre(),會計算出對應N 階Legendre polynomial 的根以及權重,第三部分是Gauss_Legendre_quadrature() 會算出Gauss-Legendre Integration數值積分。其中要特別注意的是第二段的程式中,矩陣A 與b 我們特別設定dtype 為commplex。如果少了這段,Python 在餵值到矩陣 \(A\) 的時候會自動把虛數的部份省略,這樣會造成計算結果的錯誤。由於Gauss-Legendre Integration 的預設區間為 \([-1,1]\) 對於定義在一般區間 \([a,b]\) 的函數,我們要使用變數變換的方式來計算其積分。

問題1:

為什麼做Gauss Legengre quadrature 時我們要找的是Legendre polynomial 的根?

問題2:

Legendre polynomial 彼此正交的特性為什麼可以保證數值積分的結果更靠近真實解?請用高頻與低頻的概念解釋。

Gauss-Hermite Integration

Gauss-Hermite Integration 是計算 \(\int_{-\infty}^{\infty} e^{-x^2}f(x) dx\) 這種形式的積分,可以想像當 \(f(x)\) 的定義域是正負 \(\infty\) 時,我們一定會懷疑,要怎麼選取有限的點來計算寬度是無窮大的值域的積分呢?透過 \(e^{-x^2}\) 在兩端點收斂到 \(0\) 的特性,我們可以計算這類的數值積分。Hermite polynomial 的特性是

\( \int_{-\infty}^{\infty} e^{-x^2}H_m(x) H_n(x)dx = \delta_{m,n} \),

這正交的特性可以讓我們對無窮support 的函數 \(f(x)\) 做正交基底的展開。

我們選取的取值點是Hermite polynomial 的根,Hermite polynomial 的公式是

\(H_{N}(x) = \sum_{i=0}^{floor(N/2)}\frac{(-1)^{i}}{i!} N(N-1)...(N-2i+1)(2x)^{N-2i}\)


我們用以下的遞迴公式來計算一般的Hermite polynomial

\(H_{N}(x) = 2xH_{N-1}(x)-H'_{N-1}(x)\)


若我們要計算 \(\int_{-\infty}^{\infty} f(x) dx\),我們先假設 \(g(x) = e^{x^2}f(x)\) 然後計算 \(\int_{-\infty}^{\infty} e^{-x^2}g(x) dx\) 即可。

def Hermite_poly(N):
    """
    The generator of the coefficient of Hermite polynomial of degree N

    For example:
    c = Hermite_poly(3)
    """
    if N <=0:
        c = np.array([1])
    else:
        c = np.array([2,0])
        for i in xrange(N-1):
            a = list(c)
            a.append(0)
            a = np.array(a)

            p = np.poly1d(c)
            b = [0,0]
            b.extend(list(p.deriv(1)))
            b = np.array(b)
            c = 2*a-b
    return c

def Gauss_Hermite(N):
    """
    The generator of the points and the weightings of Gauss-Hermite Integration of order N

    For exammple:
    x,w = Gauss_Hermite(3)
    """
    if N<0:
        print "Gauss-Hermite polynomial of negative order??"
        return
    
    c = Hermite_poly(N)
    p = np.poly1d(c)
    x = p.r
    
    A = np.matrix(np.zeros([N,N]),dtype='complex')
    b = np.matrix(np.zeros([N,1]),dtype='complex')
    A[0,:] = 1
    b[0,0] = np.sqrt(np.pi)
    for i in xrange(1,N):
        for j in xrange(N):
            A[i,j] = x[j]**i
        if i%2 <>0:
            b[i,0] = 0
        else:
            b[i,0] = (float(i-1)/2.0)*b[i-2,0]

    w = la.solve(A,b)
                  
    return x,w.real

 

上列程式有兩部分,第一部分Hermite_poly() 會造出N 階Hermite polynomial 的係數,第二部分Gauss_Hermite() 會產生出Gauss-Hermite Integration 的取值點和權重。使用範例如下: 

>>> x,w = Gauss_Hermite(6)
>>> def myfun(x):
return np.sin(2*x)+x**2-2

>>> y = myfun(x)
>>> S = np.dot(y,w)
>>> float(S.real)
-2.6586807763582732
>>> S
matrix([[-2.65868078]])

這裡我們令 \(f(x) = sin(2x)+x^{2}-2\),我們要計算的是 \(\int_{-\infty}^{\infty} e^{-x^2} f(x) dx\),計算結果是矩陣形式,我們用float(S.real) 取得答案。

Gauss-Laguerre Integration

Gauss-Hermite Integration 是計算 \(\int_{0}^{\infty} e^{-x} f(x) dx\) 形式的積分,我們的取值點是Laguerre polynomial 的根,Laguerre 的公式是


\(L_{n}(x) = \sum_{k=0}^{n} {n \choose k} \frac{(-1)^{k}x^{k}}{k!}\)

 

Laguerre polynomial 的特性是

\( \int_{-\infty}^{\infty} e^{-x}L_m(x) L_n(x) dx = \delta_m,n \)

有了這樣的正交基底,我們仿造上述Hermite polynomial 或是 Legendre polynomial 的討論,便可以處理 \(\int_{0}^{\infty} f(x) dx\) 的積分問題。若我們要計算 \(\int_{0}^{\infty} f(x) dx\) 的積分,我們令 \(g(x) = e^{x}f(x)\) 再使用此方法即可。

Gauss-Chebyshev Integration

Gauss-Chebyshev Integration 是處理 \(\int_{-1}^{1} \frac{1}{\sqrt{1-x^2}} f(x) dx\)  形式的積分,可以想見當要積分的函數在端點會爆開時,使用Gauss-Chebyshev 形式的積分會得到較好的結果。對應取值的點是Chebyshev polynomial 的根,公式如下

\(x_i = cos(((2i+1)\pi)/2N)\), for \(i = 0 , \dots, N-1\) 

對應的 \(w_i\) 公式如下

\(w_i = \pi/N\), for \(i = 0, \dots , N-1\)
另外一個形式的Gauss-Chebyshev Integration 是 \(\int_{-1}^{1} \sqrt{1-x^2} f(x) dx\),對應的根為

\(x_i = cos(i\frac{\pi}{N+1})\), \(w_{N,i} = \frac{\pi}{N+1}sin^2(i\frac{\pi}{N+1})\), for \(i = 1,\dots ,N\)