计算方法-数值积分
问题定义
求 \(I(f)=\int_a^bf(x)\mathrm{d}x\),但是 \(f(x)\) 的原函数很难求、\(f(x)\) 未明确给出(可能只是一张离散的数表),同时对计算精度没有过高的要求。
机械求积公式
尝试用给定的一系列节点的函数值 \(f(x_0), f(x_1), \cdots\) 来线性组合出积分的近似值,即
\[ Q(f)=\sum_{j=0}^nf(x_j)H_j \]
其中,\(x_i(i=0, 1, \cdots, n)\) 是给出函数值的 \((n+1)\) 个节点,\(H_j\) 是与前面那 \((n+1)\) 个节点有关,而与 \(f(x)\) 本身无关的一系列系数。显然,机械求积公式的关键是从给出的 \(x_i\) 确定 \(H_j\)。
代数精度
如果一个机械求积公式,它用在所有不超过 \(r\) 次的多项式(即 \(f(x)\) 是不超过 \(r\) 次的多项式)上时能精确成立,而对于 \((r+1)\) 次多项式至少有一个不能精确成立,称这个求积公式有 \(r\) 次代数精度。
例如,梯形公式 \(I_1(f)=\frac{b-a}{2}(f(a)+f(b))\) 有 1 次代数精度:
- \(f(x)=1\) 时,精确成立。
- \(f(x)=x\) 时,\(\int_a^bx\mathrm{d}x=\frac{b^2-a^2}{2}=\frac{b-a}{2}(b+a)\),精确成立。
- \(f(x)=x^2\) 时,不精确成立。
插值型求积公式
将 \(f(x)\) 在 \(x_0, x_1, \cdots, x_n\) 处 Lagrange 插值,得到 \(L_n(x)=\sum_{j=0}^nl_j(x)f(x_j)\),取
\[ H_j=\int\limits_a^bl_j(x)\mathrm{d}x \]
来作为机械求积公式的系数,可以证明这个公式至少有 \(n\) 次代数精度。事实上,如果一个机械求积公式有 \(n\) 次及以上的代数精度,它必然是插值型的。
等距节点的 Newton-Cotes 公式
当求积节点等距时,将区间 \([a, b]\) 分成 \(n\) 等分,记 \(x_0=a\),\(x_i=a+ih\),\(i=0, 1, \cdots, n\)。
记 \(x=a+th\),则
\[ \begin{aligned} H_j &= \int\limits_a^bl_j(x)\mathrm{d}x \\ &= \int\limits_a^b\prod_{i=0,i\neq j}^n\frac{(x-x_i)}{(x_j-x_i)}\mathrm{d}x \\ &= \int\limits_a^b\prod_{i=0, i\neq j}^n\frac{t-i}{j-i}\mathrm{d}(a+th) \\ &= \frac{(-1)^{n-j}h}{j!(n-j)!}\int\limits_a^b\prod_{i=0, i\neq j}^n(t-i)\mathrm{d}t \end{aligned} \]
因为 \(b-a=nh\),将上式除以 \((b-a)\),得到
\[ C_j=\frac{H_j}{b-a}=\frac{(-1)^{n-j}}{nj!(n-j)!}\int\limits_a^b\prod_{i=0, i\neq j}^n(t-i)\mathrm{d}t \]
称之为「Cotes 系数」。由此可以得到等距节点的 Newton-Cotes 公式:
\[ Q_n(f)=(b-a)\sum_{j=0}^nC_jf(x_j) \]
特殊情况
梯形公式:\(n=1, x_0=a, x_1=b\),
\[ Q_1(f)=\frac{b-a}{2}(f(a)+f(b)) \]
Simpson 公式(抛物线公式):\(n=2, x_0=a, x_2=b, x_1=\frac{a+b}{2}\),
\[ Q_2(f)=\frac{b-a}{6}(f(a)+4f(\frac{a+b}{2})+f(b)) \]
Cotes 公式:\(n=4\),
\[ Q_4(f)=\frac{2h}{45}(7f(a)+32f(a+h)+12f(a+2h)+32f(a+3h)+7f(b)) \]
收敛性
Newton-Cotes 公式并不总是收敛于积分的真值。
数值稳定性
设 \(f(x_i)\) 的计算值为 \(\tilde{f}(x_i)\),且 \(|f(x_i)-\tilde{f}(x_i)|\leqslant\varepsilon\),则
\[ |\sum_{j=0}^nH_jf(x_j)-\sum_{j=0}^nH_j\tilde{f}(x_j)|=|\sum_{j=0}^nH_j(f(x_j)-\tilde{f}(x_j))|\leqslant\varepsilon\sum_{j=0}^n|H_j| \]
如果 \(H_j\) 全是正数,有 \(|H_j|=H_j\),那么
\[ \varepsilon\sum_{j=0}^n|H_j|=\varepsilon\sum_{j=0}^nH_j=(b-a)\varepsilon \]
这个误差不会因为 \(n\) 的增大而增大,数值稳定。
如果 \(H_j\) 有正有负,则 \(\varepsilon\sum_{j=0}^n|H_j|\) 是可能随 \(n\) 变大而无限增长的,数值不稳定。
只有 \(n\leqslant7\) 和 \(n=9\) 的 Newton-Cotes 公式 \(H_j\) 是全正的。
代数精度与余项
当 \(n\) 为奇数时,NC 公式有 \(n\) 次代数精度。设 \(f\in C^{n+1}[a, b]\),则总是 \(\exists\xi\in(a, b)\),使得
\[ I(f)-Q(f)=E(f)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\int\limits_a^bp_{n+1}(x)\mathrm{d}x \]
当 \(n\) 为偶数时,NC 公式有 \((n+1)\) 次代数精度。设 \(f\in C^{n+2}[a, b]\),则总是 \(\exists\xi\in(a, b)\),使得
\[ I(f)-Q(f)=E(f)=\frac{f^{(n+2)}(\xi)}{(n+2)!}\int\limits_a^bxp_{n+1}(x)\mathrm{d}x \]
常用 NC 公式的余项(也就是误差):
- 梯形公式(\(n=1,h=b-a\)):\(-\frac{h^3}{12}f''(\eta)\);
- Simpson 公式(\(n=2,h=\frac{b-a}{2}\)):\(-\frac{h^5}{90}f^{(4)}(\eta)\);
- Cotes 公式(\(n=4,h=\frac{b-a}{4}\)):\(-\frac{8h^7}{945}f^{(6)}(\eta)\)。
复化的 Newton-Cotes 公式
所谓「复化」,是指将区间 \([a, b]\) 进行 \(n\) 等分后,在每个小区间上用一些简单的求积公式(例如梯形公式或 Simpson 公式),然后进行求和。
复化的梯形公式
记 \(h=\frac{b-a}{n}\),\(x_i=x_0+ih\)。
\[ \begin{aligned} T_n&=\sum_{i=0}^{n-1}\frac{h}{2}(f(x_i)+f(x_{i+1})) \\ &= \frac{h}{2}(f(a)+2\sum_{i=1}^{n-1}f(x_i)+f(b)) \end{aligned} \]
如果在 \(T_n\) 的基础上,将每个区间平分,即 \([x_i, x_{i+1}]\to[x_i, x_{i+\frac{1}{2}}], [x_{i+\frac{1}{2}}, x_{i+1}]\),得到 \(2n\) 等分的复化梯形公式 \(T_{2n}=\frac{h}{4}\sum_{i=0}^{n-1}(f(x_i)+2f(x_{i+\frac{1}{2}})+f(x_{i+1}))\)。
记 \(U_n=h\sum_{i=0}^{n-1}f(x_{i+\frac{1}{2}})\),有
\[ T_{2n}=\frac{1}{2}(T_n+U_n) \]
这提供了一种计算高阶复化梯形公式的方法。例如,要计算 \(T_8\),只要按 \(T_1\to T_2\to T_4\to T_8\) 的顺序算就可以了。
复化的 Simpson 公式
记 \(h=\frac{b-a}n\),\(x_i=x_0+ih\)。
\[ \begin{aligned} S_n &=\frac{h}6\sum_{i=0}^{n-1}(f(x_i)+4f(x_{i+\frac{1}{2}})+f(x_{i+1})) \\ &= \frac{h}6(f(a)+4\sum_{i=0}^{n-1}f(x_{i+\frac{1}{2}})+2\sum_{i=1}^{n-1}f(x_i)+f(b)) \\ &=\frac{1}{3}T_n+\frac{2}{3}U_n=\frac{4T_{2n}-T_n}{3} \end{aligned} \]
即,用高一阶的复化梯形公式就可以计算出复化 Simpson 公式。这种计算方式比较方便。
复化的 Cotes 公式
\[ C_n=\frac{h}{90}(7f(a)+32\sum_{i=0}^{n-1}f(x_{i+\frac{1}{4}})+12\sum_{i=0}^{n-1}f(x_{i+\frac{1}{2}})+32\sum_{i=0}^{n-1}f(x_{i+\frac{3}{4}})+14\sum_{i=1}^{n-1}f(x_i)+7f(b)) \]
使用 \(S_{2n}\) 和 \(S_n\) 表示:
\[ C_n=\frac{4^2S_{2n}-S_{n}}{4^2-1} \]