【科學計算】數值積分

數值積分法

在計算數學中，不可避免地需要對函數在給定區間上求定積分。即計算

\[\int_{a}^{b}f(x)\mathrm{d}x. \]

根據定積分的定義，定積分是積分上和與積分下和共同極限，因此只要函數可積，且區間被拆分得足夠細小，那么積分就可以近似表示為

\[\int_{a}^{b} f(x)\mathrm{d}x=\lim_{n\to \infty}\sum_{i=1}^{n}\frac{b-a}{n}\cdot f(a+\frac{b-a}{n}i). \]

因此，如果我們知道函數在足夠多的\(n\)個點處的取值，就可以用積分的定義近似計算其數值積分。不過這樣做的缺點是運算量大，且誤差難以估計。本文要介紹的Newton-Cotes積分和Gauss積分法，是兩種常用的數值積分方法。

目錄

• 數值積分法
  • 一、Newton-Cotes公式積分
    • 簡單情形：梯形積分公式與拋物線積分公式
    • Newton-Cotes公式
    • 誤差分析
    • 復化求積與逐次分半
  • 二、代數精度
    • 代數精度的定義
    • Newton-Cotes公式的代數精度
    • 證明拋物線公式的誤差估計
  • 三、Gauss求積公式
    • Gauss求積的目標
    • Gauss-Legendre求積公式

一、Newton-Cotes公式積分

利用Newton-Cotes公式對函數進行積分，是基於之前我們說過的插值理論的。在插值理論中我們說過，如果已知函數在不同的\(n+1\)個點處的取值，就可以構造一個\(n\)次多項式作為原函數的近似。而多項式的積分是很容易計算的，這樣就簡化了定積分的計算。

為方便接下來的討論，先回顧插值理論中的一些結果，設\(\omega(x)=\prod_{i=0}^{n}(x-x_i)\)。

首先，過\(n+1\)個不同的節點給出的插值多項式是唯一的，利用Lagrange插值基函數可以給出過\(x_0,x_1,\cdots,x_n\)點的插值函數為

\[P_n(x)=\sum_{i=0}^{n}\left(\prod_{j\ne i}\frac{x-x_i}{x_j-x_i} \right)f(x_i)=\sum_{i=0}^{\infty}\frac{\omega(x)}{(x-x_i)\omega'(x_i)}f(x_i) \]

對任何\(x\in [a,b]\)，用\(P_n(x)\)作為\(f(x)\)的近似的誤差估計是

\[R(x)=f(x)-p_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega(x). \]

特別當\(n=1\)時和\(n=2\)時，有

\[R_1(x)=\frac{f^{(2)}(\xi)}{2}(x-x_0)(x-x_1),\\ R_2(x)=\frac{f^{(3)}(\xi)}{6}(x-x_0)(x-x_1)(x-x_2). \]

簡單情形：梯形積分公式與拋物線積分公式

梯形積分公式基於兩個點：\((a,f(a))\)和\((b,f(b))\)，利用這兩個點構造出的插值多項式最簡單，是一條直線，可以很容易地寫出其表達式：

\[P_1(x)=f(a)+\frac{f(b)-f(a)}{b-a}(x-a), \]

對其進行積分，就得到

\[\int_{a}^{b}P_1(x)\mathrm{d}x=(b-a)f(a)+\frac{f(b)-f(a)}{b-a}\frac{(b-a)^2}{2}=(b-a)\left(\frac{1}{2}f(a)+\frac{1}{2}f(b) \right). \]

直觀上來看，這就是以\((a,f(a))\)和\((b,f(b))\)到\(x\)軸的線段為平行邊圍成的梯形的面積，因而將這個求積公式叫做梯形求積公式，它是插值多項式積分的最簡單情形。

如果我們需要稍復雜一些，增加一個點，最好的選擇是增加\(a,b\)的中點。此時，插值多項式是一個二次多項式，即一條拋物線，可以寫出其表達式為：

\[\begin{aligned} &\quad P_2(x) \\ &=\frac{(x-a)(x-\frac{a+b}{2})}{(b-a)(b-\frac{a+b}{2})}f(b)+\frac{(x-a)(x-b)}{(\frac{a+b}{2}-a)(\frac{a+b}{2}-b)}f\left(\frac{a+b}{2} \right)+\frac{(x-\frac{a+b}{2})(x-b)}{(a-\frac{a+b}{2})(a-b)}f(a)\\ &=\frac{2}{(b-a)^2}\bigg[(x-a)\left(x-\frac{a+b}{2} \right)f(b)+\\ & \qquad 2(x-a)(x-b)f\left(\frac{a+b}{2} \right)+\left(x-\frac{a+b}{2}(x-b)\right)f(a) \bigg] \end{aligned} \]

對其進行積分，得到的結果是

\[\int_{a}^{b}P_2(x)\mathrm{d}x=(b-a)\left[\frac{1}{6}f(a)+\frac{4}{6}f\left(\frac{a+b}{2} \right)+\frac{1}{6}f(b) \right]. \]

如果對比梯形求積公式，會發現，它們共同的特點是：都包含\((b-a)\)這一項，同時另一個因子都是函數值的加權平均。如果將插值多項式推廣到更高的階數，是否仍有這樣的結論？下一Part的結論將給出肯定的答復。

Newton-Cotes公式

這一部分的內容，就是對梯形求積公式和拋物線求積公式的一個推廣，它的要求是：在\([a,b]\)上選擇的分點是等間距的，即

\[x_i=a+ih,\quad i=0,1,2,\cdots,n,h=\frac{b-a}{n}. \]

（下面的推導過程可以略看，與計算過程無關）根據前言部分的結論，我們有

\[P_n(x)=\sum_{i=0}^{n}\frac{\omega(x)}{(x-x_i)\omega'(x_i)}f(x_i), \]

於是

\[\int_{a}^{b}P_n(x)\mathrm{d}x=\int_{a}^{b}\left(\sum_{i=0}^{n}\frac{\omega(x)}{(x-x_i)\omega'(x_i)}f(x_i) \right)\mathrm{d}x=\sum_{i=0}^{n}\left(\int_{a}^{b}\frac{\omega(x)}{(x-x_i)\omega'(x_i)} \right)f(x_i). \]

令

\[A_i=\int_{a}^{b}\frac{\omega(x)}{(x-x_i)\omega'(x_i)}, \]

則有

\[\int_{a}^{b}P_n(x)\mathrm{d}x=\sum_{i=0}^{n}A_if(x_i). \]

下面要計算\(A_i\)，就要用到節點等間距的要求，用\(a+ih\)替代\(x_i\)，並令\(x=a+th\)，這樣\(t\)就在\([0,n]\)中連續取值，並且

\[\omega(x)=\prod_{i=0}^{n}(x-x_i)=\prod_{i=0}^{n}[(a+ih)-(a+th)]=h^{n+1}\cdot t(t-1)\cdots(t-n),\\ \omega'(x_i)=\prod_{j\ne i}^{n}(x_i-x_j)=h^n (-1)^{n-i}[i!(n-1)!]. \]

代入就得到

\[\begin{aligned} A_i&=\int_{a}^{b}\frac{\omega(x)}{(x-x_i)\omega'(x_i)}\mathrm{d}x\\ &\xlongequal{x=a+th}\int_{0}^{n}\frac{h^{n+1}t(t-1)\cdots(t-n)}{(-1)^{n-i}h^{n}[i!(n-1)!]h(t-i)}h\mathrm{d}t \\ &=\frac{(-1)^{n-i}h}{i!(n-1)!}\int_{0}^{n}\frac{t(t-1)\cdots(t-n)}{t-i}\mathrm{d}t\\ &=(nh)\frac{(-1)^{n-i}}{n\cdot i!(n-i)!}\int_{0}^{n}\frac{t(t-1)\cdots(t-n)}{t-i}\mathrm{d}t. \end{aligned} \]

（略看部分結束）觀察\(A_i\)式，會發現\((nh)\)就是區間長度\(b-a\)，再令剩余部分為\(c_{i}^{(n)}\)，就得到積分式為：

\[\int_{a}^{b}P_n(x)\mathrm{d}x=(b-a)\sum_{i=0}^{n}c_i^{(n)}f(x_i). \]

當然，\(c_i^{(n)}\)的計算有一些繁瑣，但是注意到\(c_i^{(n)}\)的表達式中只有\(n,i\)而沒有\(h\)，這就說明\(c_i^{(n)}\)的計算不依賴於區間，只依賴於節點個數與節點的排位，因此可以事先計算出來以供查閱。這樣\(c_i^{(n)}\)被稱為Newton-Cotes系數，有表如下：

\[\begin{array}{c|cc} \hline & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 \\ \hline 1 & \frac{1}{2} & \frac{1}{2} \\ \hline 2 & \frac{1}{6} & \frac{4}{6} & \frac{1}{6} \\ \hline 3 & \frac{1}{8} & \frac{3}{8} & \frac{3}{8} & \frac{1}{8} \\ \hline 4 & \frac{7}{90} & \frac{32}{90} & \frac{12}{90} & \frac{32}{90} & \frac{7}{90} \\ \hline 5 & \frac{19}{288} & \frac{75}{288} & \frac{50}{288} & \frac{50}{288} & \frac{75}{288} & \frac{19}{288} \\ \hline 6 & \frac{41}{840} & \frac{216}{840} & \frac{27}{840} & \frac{572}{840} & \frac{27}{840} & \frac{216}{840} & \frac{41}{840} \\ \hline 7 & \frac{751}{17280} & \frac{3577}{17280} & \frac{1323}{17280} & \frac{2989}{17280} & \frac{2989}{17280} & \frac{1323}{17280} & \frac{3577}{17280} & \frac{751}{17280} \\ \hline 8 & \frac{989}{28350} & \frac{5888}{28350} & \frac{-928}{28350} & \frac{10496}{28350} & \frac{-4540}{28350} & \frac{10496}{28350} & \frac{-925}{28350} & \frac{5888}{28350} & \frac{989}{28350}\\ \hline \end{array} \]

可以發現，當\(n\le 7\)時所有的\(c_{i}^{(n)}\)都是正數，而當\(n\ge 8\)時Cotes系數出現了負數，這是我們不使用過多結點進行Newton-Cotes公式計算數值積分的原因，因為出現了負數會導致誤差估計難以確定，這是我們下一Part將討論的內容。

現在我們用Cotes公式計算\(\sin x^2\)在\([-1,1]\)上的積分，不妨取\(n=2\)和\(n=3\)，插值多項式的圖示如下（黑色為原函數，紅色為插值多項式）：

當\(n=2\)時，計算得

\[\int_{-1}^{1}\sin(x^2)\mathrm{d}x\approx 2\left[\frac{1}{6}\sin (-1)^2+\frac{4}{6}\sin(0)^2+\frac{1}{6}\sin(1^2) \right]=0.5609807, \]

當\(n=3\)時，計算得

\[\int_{-1}^{1}\sin(x^2)\mathrm{d}x\approx 2\left[\frac{1}{8}\sin(-1)^2+\frac{3}{8}\sin(-\frac{1}{3})^2+\frac{3}{8}\sin(\frac{1}{3})^2+\frac{1}{8}\sin 1^2 \right]=0.5870594. \]

而積分的更精確值是\(0.6205366\)，從圖上也可以看出紅色線要比黑色線普遍低一些，這代表對於此函數，Newton-Cotes公式積分的結果會比精確值低。此時，如果選擇\(n=1\)，結果會怎樣？

誤差分析

首先，在使用Newton-Cotes公式計算積分時，常常不會選擇太多的節點，因為\(n\ge 8\)時，Cotes系數可能會出現負數。為什么對於負數Cotes系數時一般不選擇，這是考慮到計算\(f(x_i)\)時可能出現的偏誤。一般說來，即便\(f\)已知，計算\(f(x_i)\)也會有一定的舍入誤差、截斷誤差等，這一般會有一定的絕對誤差限\(\varepsilon\)，不妨假設計算出來的有誤差的值為\(\tilde f(x_i)\)。如果Cotes系數全是正數，那么結合\(\sum_i c_i^{(n)}=1\)，有

\[\begin{aligned} &\quad \left|(b-a)\sum_{i=1}^{n}c_i^{(n)}f(x_i)-(b-a)\sum_{i=1}^{n}c_i^{(n)}\tilde f(x_i) \right| \\ &=(b-a)\left| \sum_{i=1}^{n}c_i^{(n)}[f(x_i)-\tilde f(x_i)] \right|\\ &\le (b-a)\sum_{i=1}^{n}|c_i^{(n)}||f(x_i)-\tilde f(x_i)|\\ &\le (b-a)\varepsilon \sum_{i=1}^{n}|c_i^{(n)}|\\ &=(b-a)\varepsilon. \end{aligned} \]

這使得誤差處於可以控制的范圍內，這個性質稱為數值穩定性。如果Cotes系數存在負數，這個純由計算誤差導致的誤差就可能變得難以控制。

在討論誤差估計之前，先回顧一下積分中值定理。

積分中值定理：若函數\(f(x)\)、\(g(x)\)在閉區間\([a,b]\)上連續，且\(g(x)\)在\([a,b]\)上不變號，則存在一點\(\varepsilon\in[a,b]\)，使得

\[\int_{a}^{b}f(x)g(x)\mathrm{d}x=f(\varepsilon)\int_{a}^{b}g(x)\mathrm{d}x. \]

現在，我們討論Newton-Cotes公式的誤差估計，方便起見，我們只討論比較簡單的梯形公式與拋物線公式。首先對於梯形求積公式，有如下的定理：若\(f(x)\in C^{2}[a,b]\)，則梯形求積公式有誤差估計

\[R(f)=\int_{a}^{b}f(x)\mathrm{d}x-\frac{b-a}{2}(f(a)+f(b))=-\frac{(b-a)^3}{12}f^{(2)}(\eta),\quad \eta\in(a,b). \]

顯然

\[\begin{aligned} R(f)&=\int_{a}^{b}[f(x)-P_1(x)]\mathrm{d}x\\ &=\int_{a}^{b} R_1(x)\mathrm{d}x\\ &=\int_{a}^{b}\frac{f^{(2)}(\xi)}{2}(x-a)(x-b)\mathrm{d}x. \end{aligned} \]

要注意，這里\(\xi\)是\(x\)的函數\(\xi(x)\)，所以不能貿然將其提到積分限以外，而應運用積分中值定理，因\((x-a)(x-b)\)在\([a,b]\)上非正，有

\[\int_{a}^{b}\frac{f^{(2)}(\xi(x))}{2}(x-a)(x-b)\mathrm{d}x=\frac{f^{(2)}(\eta)}{2}\int_{a}^{b}(x-a)(x-b)\mathrm{d}x. \]

對於下面的積分，我們常用以下變換：\(x=a+th\)，\(h=\frac{b-a}{n}\)。這里\(n=1\)，所以

\[\int_{a}^{b}(x-a)(x-b)\mathrm{d}x=h^3\int_{0}^{1}t(t-1)\mathrm{d}t=-\frac{(b-a)^3}{6}, \]

所以

\[R(f)=-\frac{(b-a)^3f^{(2)}(\eta)}{12}. \]

對於拋物線求積公式，有與梯形求積公式類似的、但又有所不同的誤差估計：若\(f(x)\in C^4[a,b]\)，則

\[R(f)=-\frac{(b-a)^5}{2880}f^{(4)}(\eta). \]

這是因為此時\(n\)為偶數。這里的證明，需要在下一部分給出代數精度這一概念時才能給出。

復化求積與逐次分半

在插值理論中，我們知道存在某些函數，即使節點無限加密，也難以對原始函數很好的擬合（尤其是負次數函數），即Runge現象。為了克服Runge現象，往往使用分段插值以縮小誤差。同理，如果擬合效果差，用多項式代替求積分的效果也自然不好，我們可能也會對分段插值得到的分段函數進行分段求積，這便是復化求積公式。

這里不對復化梯形公式與復化拋物線求積公式過多介紹，直接給出公式。對於梯形求積公式，把\([a,b]\)區間\(n\)等分后得到的節點\(x_k=a+kh\)，有

\[\begin{aligned} I&\approx \sum_{k=0}^{n-1}\frac{x_{k+1}-x_k}{2}(f(x_{k+1})+f(x_{k}))\\ & =\frac{h}{2}\left[f(a)+f(b)+2\sum_{k=1}^{n-1}f(a+kh)\right]\\ &\xlongequal{def}T_n. \end{aligned} \]

對於拋物線求積公式，把\([a,b]\)區間\(2n\)等分后得到的節點為\(x_k=a+kh\)，這里\(h=(b-a)/2n\)，於是

\[\begin{aligned} I&\approx \sum_{k=1}^{n}\frac{h}{3}[f(x_{2k-2})+4f(x_{2k-1})+f(x_{2k})] \\ &=\frac{h}{3}\left[ f(a)+f(b)+4\sum_{k=1}^{n}f(x_{2k-1})+2\sum_{k=1}^{n-1}f(x_{2k}) \right]\\ &\xlongequal{def} S_n. \end{aligned} \]

其誤差估計為：

\[R(f,T_n)=-\frac{b-a}{12}h^2f^{(2)}(\eta),\quad \eta\in(a,b),\\ R(f,S_n)=-\frac{b-a}{180}h^4f^{(4)}(\eta),\quad \eta\in(a,b). \]

對於梯形求積公式，有

\[\begin{aligned} R(f,T_n)&=\sum_{k=0}^{n-1}\left(-\frac{h^3}{12}f^{(2)}(\eta_k) \right),\quad \eta_k\in(x_{k},x_{k+1})\\ &=-\frac{h^3}{12}\sum_{k=0}^{n-1}f^{(2)}(\eta_k)\\ &=-\frac{h^3}{12}nf^{(2)}(\eta),\quad \eta\in(a,b)\\ &=-\frac{b-a}{12}h^2f^{(2)}(\eta). \end{aligned} \]

這里第三個等號是因為\(f^{(2)}(x)\)作為連續函數的性質。

對於拋物線求積公式，有

\[\begin{aligned} R(f,S_n)&=\sum_{k=0}^{n-1}\left(-\frac{(2h)^5}{2880}f^{(4)}(\eta_k) \right),\quad \eta_k\in(x_k,x_{k+1})\\ &=-\frac{(2h)^5}{2880}\sum_{k=0}^{n-1}f^{(4)}(\eta_k)\\ &=-\frac{2h^5}{180}nf^{(4)}(\eta),\quad \eta\in(a,b)\\ &=-\frac{b-a}{180}h^{4}f^{(4)}(\eta). \end{aligned} \]

接下來介紹逐次分半法，它是一種計算復化求積的方式。雖然有了梯形求積公式的誤差估計，但是這是一種事前估計，即為了達到我們預期的誤差限，往往會選擇超過實際所需非常多的節點，這是一種計算量上的浪費。而逐次分半法，每次對節點加密一倍，從而獲得更精確的積分結果，這樣將得到一個積分數列\(T_1,T_2,T_4,\cdots\)，當相鄰兩個積分結果小於一定值時，就不繼續計算，將最終結果作為返回值。

逐次分半法還有一個顯著優點是，重復的節點不需要重復計算。不妨設將\([a,b]\)四等分，節點是\(a,x_1,x_2,x_3,b\)，看下面的計算過程：

\[T_1=(b-a)\left(\frac{f(a)}{2}+\frac{f(b)}{2} \right),\\ T_2=\frac{(b-a)}{2}\left(\frac{f(a)}{2}+\frac{f(b)}{2}+f(x_2) \right)=\frac{T_1}{2}+\frac{b-a}{2}f(x_2),\\ T_4=\frac{b-a}{4}\left(\frac{f(a)}{2}+\frac{f(b)}{2}+f(x_1)+f(x_2)+f(x_3) \right)=\frac{T_2}{2}+\frac{b-a}{4}(f(x_1)+f(x_3)),\\ \cdots \]

在這樣的計算過程下，重復的結點沒有加入計算，所以每次只需要計算新加入的節點即可。

一般將停步的定值（即最終相鄰兩個積分值的差）設定為\(3\varepsilon\)，\(\varepsilon\)是誤差限。

二、代數精度

在Newton-Cotes公式中，最終將求積公式的形式表現為

\[\int_{a}^{b}f(x)\mathrm{d}x\approx \sum_{i=0}^{n}A_if(x_i).\int_{a}^{b}f(x)\mathrm{d}x\approx \sum_{i=0}^{n}A_if(x_i). \]

這也是我們要介紹的求積公式的一般形式，注意各個\(A_i\)是不依賴於被積函數\(f\)的。這里介紹代數精度概念，一方面是作為評判的標准，另一方面是為下面Gauss求積公式做准備。

代數精度的定義

前面所敘述的求積，都是以某種方式將被積函數\(f(x)\)用\(n\)次多項式\(P_n(x)\)代替，對多項式求積。特別在Newton-Cotes公式求積分時，使用的是插值多項式，而我們知道如果\(f(x)\)是不高於\(n\)次的多項式，則\(P_n(x)\equiv f(x)\)，這時候積分的結果是精確的；但如果\(f(x)\)是高於\(n\)次的多項式，那么\(f^{(n+1)}(x)\)不會恆等於\(0\)，用\(P_n(x)\)替代\(f(x)\)求積分就會存在誤差。代數精度這個概念，就是用來表征某種積分法是相當於用幾次多項式來替代求積的。

現給出代數精度的定義。對一個一般求積公式

\[\int_{a}^{b}f(x)\mathrm{d}x\approx \sum_{k=0}^{n}A_kf(x_k), \]

其中\(A_k\)是不依賴於函數\(f(x)\)的常數。如果求積公式中的\(f(x)\)是任意一個次數不高於\(m\)次多項式時，求積公式的約等號可改為精確成立的等號，但存在\(m+1\)次多項式\(f(x)\)使求積公式不能精確成立，就稱求積公式具有\(m\)次代數精度。

Newton-Cotes公式的代數精度

在引入Gauss求積公式之前，先對Newton-Cotes公式的代數精度進行討論。

由於Newton-Cotes公式使用\(n\)次插值多項式進行替代求值，所以對於有\(n+1\)個節點\(x_0,x_1,\cdots,x_n\)的Cotes公式，Newton-Cotes求積至少應當具有\(n\)次代數精度，因為當\(f(x)\in\mathcal P_n(\mathbb{R})\)時\(f(x)\equiv P_n(x)\)。比如，梯形求積公式對於一次函數的積分是精確的，而拋物線求積公式對於二次函數的積分也是精確的。

但值得注意的是，Newton-Cotes公式對於\(n\)為偶數的情形，可以具有\(n+1\)次代數精度。

記任意\(n+1\)次多項式為

\[q_{n+1}(x)=\sum_{k=0}^{n+1}b_{k}x^{k},\quad q_{n+1}^{(n+1)}(x)=(n+1)!b_{n+1}, \]

用\(P_n(x)\)代表其\(n\)次插值多項式，於是其誤差為

\[R_n(x)=\frac{q_{n+1}^{(n+1)}(\xi)}{(n+1)!}\omega_n(x)=b_{n+1}\omega_n(x), \]

在前面，我們已經得到了這樣的結果：

\[\int_{a}^{b}\omega_n(x)\mathrm{d}x=h^{n+2}\int_{0}^{n}t(t-1)\cdots(t-n)\mathrm{d}t, \]

接下來的目標是證明上面的積分為\(0\)，這樣就代表用\(P_n(x)\)替代\(q_{n+1}(x)\)的積分是精確的。令\(n=2k\)並引進變換\(u=t-k\)，有

\[\begin{aligned} &\quad \int_{0}^{n}t(t-1)\cdots(t-n)\mathrm{d}t\\ &=\int_{0}^{2k}t(t-1)\cdots(t-k)(t-k-1)\cdots(t-2k+1)(t-2k)\mathrm{d}t\\ &=\int_{k}^{k}(u+k)(u+k-1)\cdots u(u-1)\cdots(u-k-1)(u-k)\mathrm{d}u, \end{aligned} \]

令積分號里面的部分為\(H(u)\)，即\(H(u)=(u+k)\cdots u(u-1)\cdots(u-k)\)，有

\[\begin{aligned} H(-u)&=(-u+k)(-u+k+1)\cdots(-u)(-u-1)\cdots(-u-k) \\ &=(-1)^{2k+1}H(u)\\ &=-H(u), \end{aligned} \]

從而\(H(u)\)是一個奇函數，這說明積分結果是\(0\)，所以

\[\int_{a}^{b}\omega_n(x)\mathrm{d}x=0, \]

即\(P_n(x)\)代替\(q_{n+1}(x)\)使用Newton-Cotes求積公式是精確的。

證明拋物線公式的誤差估計

之前，我們給出了梯形求積公式與拋物線求積公式的誤差估計，並證明了梯形求積公式的誤差估計，但沒有給出拋物線公式的誤差估計。這是因為拋物線求積公式具有三次代數精度，這使得它可以用更高次的插值多項式逼近。

再回顧拋物線公式的誤差估計：

\[R(f)=-\frac{(b-a)^5}{2880}f^{(4)}(\eta),\quad \eta\in(a,b). \]

由於可以構造三次多項式，故除了\(f(a),f(\frac{a+b}{2}),f(b)\)以外，還可以增加一個條件，這里選擇的是\(f'(\frac{a+b}{2})\)。利用Newton插值法，可知\(f'(\frac{a+b}{2})\)實際上是在\((\frac{a+b}{2})\)處選擇了重節點。

構造一個三次插值多項式\(P_3(x)\)滿足

\[P_3(a)=f(a),\quad P_3(b)=f(b),\\ P_3\left(\frac{a+b}{2} \right)=f\left(\frac{a+b}{2} \right),\\ P_3'\left(\frac{a+b}{2} \right)=f'\left(\frac{a+b}{2} \right). \]

從而插值函數的誤差是

\[f(x)-P_3(x)=\frac{f^{(4)}(\xi(x))}{4!}(x-a)\left(x-\frac{a+b}{2} \right)^2(x-b), \]

易知后面的部分在\([a,b]\)上不變號。積分，由積分中值定理，並引入變換\(x=a+t(b-a)\)，就有

\[\begin{aligned} &\quad \int_{a}^{b}\frac{f^{(4)}(\xi(x))}{4!}(x-a)\left(x-\frac{a+b}{2} \right)^2(x-b)\mathrm{d}x\\ &=\frac{f^{(4)}(\eta)}{4!}\int_{a}^{b}(x-a)\left(x-\frac{a+b}{2} \right)^2(x-b)\mathrm{d}x\\ &=\frac{(b-a)^5f^{(4)}(\eta)}{4!}\int_0^1 t\left(t-\frac{1}{2}\right)^2 (t-1)\mathrm{d}t\\ &=-\frac{(b-a)^5f^{(4)}(\eta)}{4!\cdot 120}. \end{aligned} \]

容易計算最后一個積分值為\(-1/120\)，且\(4!\cdot 120=2880\)，故結論得證。

三、Gauss求積公式

Gauss型求積公式是對Newton-Cotes求積公式的改進。Newton-Cotes求積公式基於插值多項式推出，選擇了最方便的等距節點，並且根據\(A_i\)中與區間無關的部分編制了Cotes系數表，從而使計算變得方便。

Gauss型求積公式則是為了更精確的積分結果，選擇更復雜的節點與對應的系數\(A_i\)。如何獲得更精確的積分結果，就是我們接下來需要考慮的問題。

Gauss求積的目標

先明確Gauss型求積公式的目標：在節點數\(n\)給定時，最大化代數精度。要最大化代數精度，首先就要知道代數精度最高能達到多少，這就成為一個規划問題：

\[\max m,\\ \mathrm{s.t.}\forall q\in\mathcal P_{m}(\mathbb{R}),\quad \sum_{k=1}^{n}A_kq(x_k)=\int_{a}^{b}q(x)\mathrm{d}x. \]

寫出\(m\)次多項式的一般形式：

\[q_m(x)=\sum_{i=0}^{m}a_ix^{i}, \]

於是

\[\int_{a}^{b}q_m(x)\mathrm{d}x=\sum_{i=0}^{m}\int_{a}^{b}a_ix^{i}\mathrm{d}x=\sum_{i=0}^{m}a_i\int_{a}^{b}x^{i}\mathrm{d}x\xlongequal{def}\sum_{i=0}^{m}a_i\mu_i, \]

這里\(\int_{a}^{b}x_i\mathrm{d}x\)是定值，故使用\(\mu_i\)替代，從而我們得到一個等式：

\[\sum_{i=0}^{m}a_i\mu_i=\sum_{k=1}^{n}A_kq(x_k)=\sum_{k=1}^{n}\sum_{i=0}^{m}A_ka_ix_k^{i}=\sum_{i=0}^{m}a_i\sum_{k=1}^{n}A_kx_k^{i}. \]

注意這個等式中的幾個量：

• \(a_i\)是任意實數。因為一組\(a_i\)唯一決定一個\(q_m(x)\)，所以\(a_i\)任意取值等式都應該成立。
• \(\mu_i\)是定值。因為\(\mu_i=\int_{a}^{b}x^{i}\mathrm{d}x\)。
• \(A_k\)和\(x_k\)是待取的定值，即取一組\(A_k,x_k\)使得對任何\(a_i\)都有等式成立。

由於\(a_i\)是任意實數，所以可以將等式展開成方程組：\(\mu_i=\sum_{k=1}^{n}A_kx_{k}^{i},i=0,1,\cdots,m\)，即

\[\begin{array}{rl} A_1+A_2+\cdots+A_n&=\mu_0,\\ A_1x_1+A_2x_2+\cdots+A_nx_n&= \mu_1,\\ A_1x_1^2+A_2x_2^2+\cdots+A_nx_n^2&=\mu_2,\\ \cdots & \cdots \\ A_1x_1^{m}+A_1x_2^{m}+\cdots+A_nx_n^{m}&=\mu_m. \end{array} \]

因未知數個數為\(2n\)，為使這個方程組成立，必有\(m\le 2n-1\)，且如果存在一組\(A_k,x_k\)使得\(m=2n-1\)，則\(2n-1\)就是\(n\)個互異節點能達到的最大代數精度。事實上，這樣的\(A_k,x_k\)必定是可以取到的。

接下來討論如何選擇節點，即解上面的非線性方程組。為方便起見，我們將定積分的上下限控制為\([-1,1]\)，這樣有

\[\int_{a}^{b}q(x)\mathrm{d}x\xlongequal{x=\frac{b+a}{2}+\frac{b-a}{2}t}\frac{b-a}{2}\int_{-1}^{1}q\left(\frac{b+a}{2}+\frac{b-a}{2}t \right)\mathrm{d}t. \]

因為我們假定\(q(x)\)是\(2n-1\)次多項式，所以用一個含\(x_1,\cdots,x_n\)的\(n\)次多項式，這里就取為\(\omega(x)=\prod_{i=1}^{n}(x-x_i)\)去除\(q(x)\)，就有

\[q(x)=p(x)\omega(x)+r(x), \]

這里多項式除法的存在唯一性保證了\(p(x),r(x)\)都是唯一的\(n-1\)次多項式，且

\[\int_{-1}^{1}q(x)\mathrm{d}x=\int_{-1}^{1}p(x)\omega(x)\mathrm{d}x+\int_{-1}^{1}r(x)\mathrm{d}x. \]

現在，最關鍵的一步在於，我們希望對任何不超過\(n-1\)次的多項式\(p(x)\)，都有

\[\int_{-1}^{1}p(x)\omega(x)\mathrm{d}x=0,\quad (1), \]

這等價於以下\(n\)個算式：

\[\int_{-1}^{1}(x-x_1)\cdots(x-x_n)=0,\\ \int_{-1}^{1}x(x-x_1)\cdots(x-x_n)=0,\\ \cdots \\ \int_{-1}^{1}x^{n-1}(x-x_1)\cdots(x-x_n)=0. \]

一般從中解出\(x_1,x_2,\cdots,x_n\)的值，可以達到降維的效果，求解難度也大幅下降；並且，從這里可以看出，\(x_1,\cdots,x_n\)的值是不依賴於函數\(f\)的。這時候

\[\int_{-1}^{1}p(x)\mathrm{d}x=\int_{-1}^{1}r(x)\mathrm{d}x, \]

再利用求積公式對\(f(x)=x^{i}\)，\(i=0,1,\cdots,n-1\)精確成立，代入已經求得的\(x_1,\cdots,x_n\)的值，就可以解出\(A_1,\cdots,A_n\)。事實上，可以解得

\[A_k=\int_{0}^{1}\frac{\omega(x)}{(x-x_k)\omega'(x_k)}\mathrm{d}x. \]

這相當於把一個\(2n\)個未知數的非線性方程，轉化成一個含\(n\)個未知數的非線性方程和一個\(n\)個未知數的線性方程組求解。

Gauss-Legendre求積公式

在之前的討論中，我們使用\(\omega(x)\)，利用\(x^{i}\omega(x)\)的積分在\([-1,1]\)上為\(0\)來求解\(x_1,\cdots,x_n\)。而Gauss-Legendre求積公式，則將\(\omega(x)\)與Legendre多項式建立聯系來作為除函數，這樣做的好處是，對任何\(i=0,1,\cdots,n-1\)，因Legendre多項式的性質都有

\[\int_{-1}^{1}x^{i}\frac{\mathrm{d}^n}{\mathrm{d}x^n}(x^2-1)^{n}\mathrm{d}x\xlongequal{def}\int_{-1}^{1}x^{i}P(x)\mathrm{d}x=0, \]

由Legendre多項式的性質，如果取它的\(n\)個零點作為\(x_1,\cdots,x_n\)，則

\[\frac{\mathrm{d}^n}{\mathrm{d}x^{n}}(x^2-1)^{n}=(x-x_1)(x-x_2)\cdots(x-x_n). \]

並且

\[A_k=\frac{2(1-x_{k}^2)}{[nP_{n-1}(x_k)]^2},k=1,2,\cdots,n. \]

這些系數被制成表，需要使用時可以直接查閱。

\[\begin{array}{|c|c|c:c|c|c|} \hline n & x_k & A_k & n & x_k & A_k \\ \hline 2 & \pm 0.5773503 & 1.0000000 & 8 & \pm 0.9602899 & 0.1012285 \\ 4 & \pm 0.8611363 & 0.3478548 & & \pm 0.7966665 & 0.2223810 \\ & \pm 0.3399810 & 0.6521452 & & \pm 0.5255324 & 0.3137066 \\ 6 & \pm 0.9324695 & 0.1713245 & & \pm 0.1834346 & 0.3626838 \\ & \pm 0.6612094 & 0.3607616 \\ & \pm 0.2386192 & 0.4679139 \\ \hline \end{array} \]

在實際應用Gauss-Legendre公式求積時，一般步驟是先將其作積分限變換轉化為\([-1,1]\)上的積分，再根據公式代入求解。

現使用Gauss-Legendre公式計算\(\sin(x^2)\)在\([-1,1]\)上的積分，之前已經計算過其積分的准確值大約是\(0.6205366\)。當\(n=4\)時，有

\[\int_{-1}^{1}\sin (x^2)\mathrm{d}x\approx \sum_{k=1}^{4}A_k\sin(x_k^2)=0.6203309. \]

相比Newton-Cotes公式得到的結果，Gauss-Legendre積分顯然更為精確。

盡管Gauss-Legendre求積有着更高的精度，但是隨着\(n\)的增大，節點\(x_k\)的變化並沒有什么規律，所以\(f(x_k)\)的值不能重復利用。