高斯積分

前情提要: 數值積分初探

前情提要: 以內插多項式來做數值積分

我們想要找到一些插值點使得以內插多項式來做數值積分會最準. 這就是高斯積分.

Goal: 任意給定一可積分函數 $f(x)$, $x\in[-1, 1]$, 我們想要算 $\int^1_{-1} f(x) dx$.

複習一下內插多項式的積分

給定任意 $n+1$ 個不重複的點在 $[-1, 1]$ 區間, 記為 $x_0, x_1, \cdots, x_n.$ 我們可以逼近函數 $f(x)$ 的積分 $$ \int^1_{-1} f(x)dx \approx \sum^n_{i=0} w_i f_i, $$ 其中 $f_i = f(x_i)$, 並且 $$ w_i = \int^1_{-1} \ell_i(x)dx = \int^1_{-1} \left(\prod^n_{k=0, k\ne i} \frac{x-x_k}{x_i-x_k}\right) dx. $$

一個重要的 remark: 如果 $p(x)$ 是個 $n$ 次或以下的多項式, 那這個積分是等號!

也就是說 $$ \int^1_{-1} p(x)dx = \sum^n_{i=0} w_i f_i, $$ where $p(x)$ is a polynomial with degree less than or equals to $n$.

Idea:

高斯積分的想法就是, 我想要找到特別的點 $x_i$ 以及特別的權重 $w_i$ 使得我能算準的多項式次數最高.

如果我們回想一下內插多項式, 給定點然後求 $n+1$ 個權重, 也就是有 $n+1$ 個未知數, 那顯然它所能滿足的方程式就是 $n+1$ 個. 如果要把多項式算準的話從常數($0$次)到 $n$次共可以列出剛好 $n+1$個方程式. 所以最多能把 $n$ 次多項式算準.

現在假設點($x_i$)跟權重($w_i$)都是待決定的, 那我就有 $2(n+1)$ 個未知數, 能夠把多項式從常數($0$次)到 $2n+1$次都算準. 所以理論上, $n+1$ 個點的高斯積分可以將 $2n+1$次的多項式算準, 完全沒有誤差.

那如果這這樣的想法直接列式的話就會得到以下 $2n+1$ 個方程式 $$ \int^1_{-1} x^m dx=\frac{1 - (-1)^{m+1}}{m+1} = \sum^n_{i=0} w_i x^m_i, \quad 0\le m\le 2n+1. $$ 然後解以上的非線性聯立方程式就可以得到所有的點 $x_i$ 及權重 $w_i$.

看起來有點困難, 不過好消息是其實不難解.

我們需要引進正交多項式 (orthogonal polynomials)

在線性代數裡我們知道兩個向量如果正交則其內積等於零. 函數空間裡我們也可以做一樣的定義, 我們定義兩個函數的內積 $$ (f, g) = \int^1_{-1} f(x)g(x) dx. $$

Definition

$f(x)$ and $g(x)$ are orthogonal if $(f,g)=0$.

我們知道 $n$ 次多項式有一組最基本的基底 ${1, x, x^2, \cdots, x^n}$, 將這組基底以 Gram-Schimidt 做正交化就可以得到一組正交基底. 這組正交基底就是所謂的 Legendre polynomials.

我們將 Legendre polynomials 這組正交基底記為 $P_i(x)$, 其前三個如下: $$ P_0(x) = 1, \quad P_1(x) = x, \quad P_2(x)=\frac{1}{2}(3x^2-1). $$

Legendre polynomials 有以下兩個性質:

1. $P_i(x)$ 是一個 $i$次多項式
2. $(P_i(x), P_j(x)) = \delta_{ij} = \begin{cases} 0, \quad i\ne j, \\
   1, \quad i=j. \end{cases}$
   • 對任何多項式 $Q(x)$ 如果其次數小於 $i$, 則有 $(P_i(x), Q(x))=0$

高斯積分的點與權重:

任意給定一個 $2n+1$ 次多項式 $p(x)$, 我們可以透過多項式的除法把它分解成以下這種樣子 $$ p(x) = P_{n+1}(x)Q(x) + R(x), $$ 其中 $Q(x)$ 以及 $R(x)$ 都是至多 $n$ 次的多項式. 透過這樣的分解我們發現 $$ \int^1_{-1} p(x)dx = \int^1_{-1} P_{n+1}(x)Q(x) + R(x)dx = \int^1_{-1} R(x)dx, $$ 其中 $P_{n+1}$ 與 $Q$ 相乘後的積分等於零是用到上面正交多項式的性質.

以上這積分如果以數值積分來表示的話就是 $$ \int^1_{-1} p(x)dx \approx \sum^n_{i=0} w_i p(x_i) = \sum^n_{i=0} w_i \left(P_{n+1}(x_i)Q(x_i) + R(x_i)\right). $$ 我們希望數值積分也會有 “$P_{n+1}$ 與 $Q$ 相乘的積分等於零” 這件事. 那一個最簡單的取法就是我們取 $x_i$ 是 $P_{n+1}(x)$ 這個 $n+1$ 次多項式所有的根. 由於 $P_{n+1}$ 剛好有 $n+1$ 個不重複的根, 這樣我們就有 $x_i$ 這些點了, 而且 $P_{n+1}(x_i)=0$ 表示 $$ \sum^n_{i=0} w_i P_{n+1}(x_i)Q(x_i)=0, $$ 所以 $$ \int^1_{-1} p(x)dx \approx \sum^n_{i=0} w_i R(x_i), $$ 也就是說我們只需要把 $R(x)$ 這個 $n$ 次多項式算準就好. 那我們就可以用內插多項式算權重的方法把權重都寫下來 $$ w_i = \int^1_{-1} \left(\prod^n_{k=0, k\ne i} \frac{x-x_k}{x_i-x_k}\right) dx. $$

如此, 我們就把點跟權重都求出來了.

Summary:

最後總結一下:

• $n+1$ 個點的高斯積分可以將 $2n+1$ 次多項式算準無誤差

• 高斯點是正交多項式 $P_{n+1}(x)$ 的 $n+1$ 個根

• 求出高斯點後我們以內插多項式的做法來求出權重

舉個簡單的例子來試試

假設我們只想要用兩個點來估算積分, 那 Legendre polynomial 的 $P_2(x)=\frac{1}{2}(3x^2-1)$.

由此可算出兩個高斯點, 也就是 $P_2$ 的兩個根是 $$ x_0 = -\frac{1}{\sqrt{3}}, \quad x_1 = \frac{1}{\sqrt{3}}. $$ 接著權重也可以簡單的算出來, 我們得到. $$ w_0 = w_1 = 1. $$ 因此, 兩個點的高斯積分公式是 $$ \int^1_{-1} f(x)dx = f\left(-\frac{1}{\sqrt{3}}\right) + f\left(\frac{1}{\sqrt{3}}\right), $$ 可以將 3次(含)以下的多項式算準無誤差.

再延伸一下

對於更廣義有權重的積分 $$ \int^1_{-1} w(x)f(x)dx, $$ 我們只要找出相對應的正交多項式, 我們一樣可以求出其高斯積分點以及權重.

舉例來說, $w(x)=\frac{1}{1+x^2}$

其相對應的內積以及積分為 $$ (f,g) = \int^1_{-1} \frac{1}{1+x^2}f(x)g(x)dx, \quad \int^1_{-1} \frac{1}{1+x^2}f(x)dx. $$ 其正交多項式為 Chebyshev polynomials, $T_i(x)$, 前三項為: $$ T_0(x) = 1, \quad T_1(x) =x, \quad T_2(x) = 2x^2-1. $$ 比較特別的是其高斯點, 稱為 Gauss-Chebyshev points, 可以直接求出來: $$ x_i = \cos\left(\frac{(i+\frac{1}{2})\pi}{n+1}\right), \quad 0\le i\le n. $$ 而且權重也可以直接求出來: $$ w_i = \frac{\pi}{n+1}. $$