數值積分

誤差分析

 

分析數值積分的工具很自然地還是用泰勒展開式,我們令
(式一) : \(g(x) = \int_{x_k}^x f(t)dt \) ,
則我們有 \(g'(x) = f(x)\), \(g^{(2)}(x) = f'(x)\), \(g^{(3)}(x) = f^{(2)}(x)\)。把 \(g(x)\) 對 \(x_k\) 做展開

 

\[g(x) = g(x_k) + g'(x_k)(x-x_k) +\frac{1}{2}g^{(2)}(x_k)(x-x_k)^2 + \frac{1}{3!}g^{(3)}(x_k)(x-x_k)^3+\dots\]

將(式一)帶入,並且用 \(x_{k+1}\) 替換 \(x\) , \(h\) 替換 \(x_{k+1}-x_k\),我們得到

\[\int_{x_k}^{x_{k+1}} f(x) dx = hf(x_k) + \frac{h^2}{2}f'(x_k)+\frac{h^3}{3!}f^{(2)}(x_k)+ \dots\]

我們用同樣的方法,把 \(x\) 用 \(x_{k-1} \) 替換,我們得到類似的式子
\[ \int_{x_k}^{x_{k-1}} f(x) dx = -hf(x_k) + \frac{h^2}{2}f'(x_k)-\frac{h^3}{3!}f^{(2)}(x_k)+\dots\]

整合上述兩式我們得到

\[ \int_{x_{k-1}}^{x_{k+1}} f(x) dx = 2hf(x_k) + \frac{2h^3}{3!}f^{(3)}(x_k)+\frac{2h^5}{5!}f^{(4)}(x_k)+\dots\]

若我們改變間距為原來的兩倍,令 \( x_{mk} = \frac{x_k+x_{k+1}}{2} \),得到
\[\int_{x_k}^{x_{k+1}} f(x) dx - hf(x_{mk}) = \frac{h^3}{24}f^{(2)}(x_{mk}) + h^{(4)}f(x_{mk}) + \dots = O(h^3) \]

因此,當我們計算 \(\int_a^b f(x)dx\) 時,得到
\[| \int_a^b f(x) dx-\sum hf(x_{mk}) | = n*O(h^3) = \frac{b-a}{h}O(h^3) = O(h^2) \]

當我們用類似的方法將 \(f_{k-1}\), \(f_{k+1}\) 對 \(x_k\) 展開,並且代入Simpson's rule,我們得到

\[\int_{x_{k-1}}^{x_{k+1}} f(x) dx-\frac{h}{3}(f_{k-1}+f_{k}+f_{k+1}) = -\frac{h^5}{90}f^{(4)}(x_k)+ \dots = O(h^5)\]

因此,用 Simpon's rule 計算數值積分時的 truncation error 是 \(O(h^4)\)。

問題:

試用類似的手法計算梯形法的truncation error,證明梯形法和中點法有相同的誤差階數。 

 

這裡有兩件事情可能已經引起讀者的興趣了。 

- 數值積分得到的 truncation error 的次數,比數值微分高。
- 用不同階數的函數去逼近 \(f(x)\) 竟然得到相同階數的 truncation error。


對第一點的觀察是,我們用常數函數去取代 \(f(x)\) 計算出來的誤差階數是 2,而我們用2 次函數去取代 \(f(x)\) 所得的誤差階數是4,這種誤差階數高於逼近函數的階數的現象是一個通則。也因此,我們在數值解的解題習慣上會喜歡積分方程甚於微分方程,雖然積分方程所對應的解往往是一種弱解的形式,但是在計算上會比較準確。 

對第二點的觀察是,在 midpoint rule 和 trapezoidal rule 中,雖然代表的逼近函數階數不同,但他們的誤差階數卻是相同的。這裡宣示了一件特別的現象,就是當我們做積分時,sampling 的點若選取的好,我們會得到比預期更好的 truncation error order。在 Gauss Quadrature 這一節中我們會詳細介紹這個特性。

當我們知道truncation error 的階次之後,由於我們並不會認真去算(或是實務上無法計算)\(f(x)\) 的高階導數,因此 truncation error 的領導係數是多少,是我們需要去估計的。估計的方式,和數值微分一樣。我們透過更改 step size 的大小,觀察兩者計算結果的差別來估計 truncation error 的領導係數。以 Simpson's rule 為例,我們令數值積分為 \(I_S(a,b,h) \simeq \int_a^b f(x) dx\),則

\[ \int_a^b f(x) dx - I_S(a,b,h) = O(h^4) = C h^4 \]
同理

\[ \int_a^b f(x) dx - I_S(a,b,\frac{h}{2}) = O(\frac{h}{2}^4) = C \frac{h^4}{16} \]

因此
\[| I_S(a,b,\frac{h}{2})-I_S(a,b,h) | = \frac{15}{16}h^4|C| \]

我們得到
\[|C| \simeq \frac{16}{15}h^{-4} | I_S(a,b,\frac{h}{2})-I_S(a,b,h) |\]。

這一個式子右端是我們透過數值積分計算出來的結果,有了這個領導係數的估計值,我們便可以知道
(EA:式一):\(| \int_a^b f(x) dx - I_S(a,b,h) | \simeq \frac{1}{15} | I_S(a,b,\frac{h}{2})-I_S(a,b,h) | \)。

這是一個很棒的結果,雖然我們不知道真解的積分值,但透過這個公式,我們觀察 \(I_S(a,b,h)\) 和 \(I_S(a,b,\frac{h}{2}) 的誤差,就可以控制數值解和真解的距離。若我們使用 trapezoidal rule 來計算積分,我們使用下列誤差估計公式來控制我們的計算結果。

 

(EA:式二):\(| \int_a^b f(x) dx - I_T(a,b,h) | \simeq \frac{1}{3} | I_T(a,b,\frac{h}{2})-I_T(a,b,h) | \)。

問題:

請計算為什麼梯形法估計誤差時的對應係數是 \(\frac{1}{3}\)?