计算方法-数值积分

问题定义

\(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} \]


计算方法-数值积分
https://eupho.me/a558a3e4.html
作者
Lambert Swizzer
发布于
2023年5月12日
许可协议