数值分析：数值积分

• 数值分析：数值积分
  • 数值积分定义
    • 近似定积分
    • 数值求积公式
      • 例1
      • 例2
  • 插值型求积公式
  • 牛顿科特斯（Newton-Cotes）求积公式
    • 梯形公式
    • 辛普森（Simpson）公式
    • 科特斯（Cotes）公式
      • 例
  • 复化求积公式
    • 复化梯形公式
    • 复化辛普森公式
      • 例
    • 龙贝格（Romberg）求积公式
    • 复化Cotes公式
  • 高斯型求积公式
  • 差商型求积公式
    • 一阶向前差商
    • 一阶向后差商
    • 一阶中心差商
    • 二阶中心差商

数值积分定义

大多数定积分问题都不能通过解析法，即计算牛顿一莱布尼茨公式$\int_{a}^{b}{f(x)}\,{\rm d}x = F(b) - F(a)$来求解，数值积分方法就是将$f(x)$用简单函数近似替代（比如泰勒展开），利用离散数据求函数的数值积分，从而避免了寻求原函数的困难，并为计算机求积分提供了可行性。

近似定积分

左矩形

$\int_{a}^{b}{f(x)}\,{\rm d}x = f(a)(b-a) + \dfrac{(b-a)^2}{2} f'(η) { \ \ \ } η∈(a,b)$

右矩形

$\int_{a}^{b}{f(x)}\,{\rm d}x = f(b)(b-a) - \dfrac{(b-a)^2}{2} f'(η) { \ \ \ } η∈(a,b)$

中矩形

$\int_{a}^{b}{f(x)}\,{\rm d}x = f(\dfrac{a+b}{2})(b-a) -\dfrac{(b-a)^3}{24} f''(η), { \ \ \ } η∈(a,b)$

数值求积公式

例1

确定下式的代数精度

• $\int_{-1}^{1}f(x){\rm d}x ≈ \dfrac{f(-1) + 2f(0) + f(1)}{2}$
• $\int_{-1}^{1}f(x){\rm d}x ≈ f(-\dfrac{1}{\sqrt{3}}) + f(\dfrac{1}{\sqrt{3}})$

解

记

• $I(f) = \int_{-1}^{1}f(x){\rm d}x$
• $I_1(f) = \dfrac{f(-1) + 2f(0) + f(1)}{2}$
• $I_2(f) = f(-\dfrac{1}{\sqrt{3}}) + f(\dfrac{1}{\sqrt{3}})$

$f$	$I(f)$	$I_1(f)$	$I_2(f)$
$1$	$\int^{1}_{-1}1{\rm d}x=2$	$\frac{1+2+1}{2}=2$	$1+1=2$
$x$	$\int^{1}_{-1}x{\rm d}x=0$	$\frac{-1+0+1}{2}=0$	$\frac{-1+1}{\sqrt{3}}=0$
$x^2$	$\int^{1}_{-1}x^2{\rm d}x=\frac{2}{3}$	$\frac{1+0+1}{2}=1$	$\frac{1+1}{3}=\frac{2}{3}$
$x^3$	$\int^{1}_{-1}x^3{\rm d}x=0$	INVALID	$\frac{1+1}{3\sqrt{3}}=0$
$x^4$	$\int^{1}_{-1}x^4{\rm d}x=\frac{2}{5}$	INVALID	$\frac{1+1}{9}=\frac{2}{9}$

综上可知，$I_1(f)$为1次代数精度，$I_2(f)$为3次代数精度。

例2

试确定$\int_{0}^{h} f(x){\rm d}x ≈ A_0f(0) + A_1f(h) + A_2f'(0)$的系数$A_0,A_1,A_2$，使该多项式具有尽可能高的代数精度。

解

$f$	$\int_{0}^{h} f(x){\rm d}x$	$A_0f(0) + A_1f(h) + A_2f'(0)$
$1$	$h$	$A_0+A_1$
$x$	$\frac{h^2}{2}$	$0+A_1h+A_2$
$x^2$	$\frac{h^3}{3}$	$0+A_1h^2+0$
$x^3$	$\frac{h^4}{4}$	$0+A_1h^3+0$

从而得$\begin{cases} h = A_0+A_1 \\ \frac{h^2}{2} = A_1h+A_2 \\ \frac{h^3}{3} = A_1h^2 \end{cases}$，解得$\begin{cases} A_0 = \frac{2}{3}h \\ A_1 = \frac{1}{3}h \\ A_2 = \frac{1}{6}h^2 \end{cases}$

$f$	$\int_{0}^{h} f(x){\rm d}x$	$A_0f(0) + A_1f(h) + A_2f'(0)$
$x^3$	$\frac{h^4}{4}$	$A_1h^3=\frac{1}{3}h^4$

故多项式最高可达2次代数精度，此时系数为$\begin{cases} A_0 = \frac{2}{3}h \\ A_1 = \frac{1}{3}h \\ A_2 = \frac{1}{6}h^2 \end{cases}$

插值型求积公式

拉格朗日（Lagrange）插值多项式：
$L_n(x) = \sum\limits_{i=0}^{n} l_i(x)f(x_i)$

• $l_i(x) = \prod\limits_{^{j=0}_{j≠i}}^{n} \dfrac{x-x_j}{x_i-x_j}$
• $R_n(x) = f(x) - L_n(x) = \dfrac{f^{(n+1)}(ξ)}{(n+1)!} ω_{n+1}(x), { \ \ \ } ξ∈(a,b)$

牛顿科特斯（Newton-Cotes）求积公式

当插值节点等距时，得到得插值型求积公式即为牛顿科特斯（Newton-Cotes）求积公式。

$\int_a^b f(x) \,{\rm d}x ≈ \sum\limits_{i=0}^{n} A_i f(x_i) = (b-a)\sum\limits_{i=0}^{n} C_i^{(n)}f(x_i)$

• $A_i = \int_{a}^{b}{ l_i(x) }\,{\rm d}x \xrightarrow[(h=\frac{b-a}{n},\,i=0,1,\cdots,n)]{x_i=a+ih} \, \xrightarrow[]{x=a+th} \dfrac{(-1)^{n-i}h}{i!(n-i)!} \int_0^n{ \prod\limits_{^{j=0}_{j≠i}}^n (t-j) }\,{\rm d}t$
  • $\sum_{i=0}^{n} A_i = b-a$
• $C_i^{(n)}=\dfrac{1}{b-a}A_i$，称为Cotes系数
  • $\sum_{i=0}^{n} C_i^{(n)} = 1$

梯形公式

当牛顿科特斯（Newton-Cotes）求积公式的$n=1$，即2个节点时

辛普森（Simpson）公式

当牛顿科特斯（Newton-Cotes）求积公式的$n=2$，即3个节点时

科特斯（Cotes）公式

当牛顿科特斯（Newton-Cotes）求积公式的$n=4$，即5个节点时

$\int_a^b f(x) {\rm d}x ≈ \dfrac{b-a}{90} \left[ 7f(x_0) + 32f(x_1) + 12f(x_2) + 32f(x_3) + 7f(x_4) \right] = C$

• $C_{0}^{(4)} = \frac{7}{90}$；$C_{1}^{(4)} = \frac{32}{90}$；$C_{2}^{(4)} = \frac{12}{90}$；$C_{3}^{(4)} = \frac{32}{90}$；$C_{4}^{(4)} = \frac{7}{90}$；

例

用Simpson公式求$y(x) = \int_0^1 (1-3xt)y(t) {\,{\rm d}t} +(1-3x)$

解

记$f(t) = (1-3xt)y(t)$，由Simpson公式得

$\int_0^1 (1-3xt)y(t) {\,{\rm d}t} ≈ \dfrac{1}{6}\left[ f(0)+4f(\dfrac{1}{2})+f(1) \right] = \dfrac{y(0) + 4(1-\dfrac{3}{2}x)y(\frac{1}{2}) + (1-3x)y(1)}{6}$

因此近似解

$\tilde{y}(x) = \dfrac{y(0)+4(1-\dfrac{3}{2}x)y(\frac{1}{2})+(1-3x)y(1)}{6} + (1-3x)$

取$x=0,\frac{1}{2},1$得

$\begin{cases} \tilde{y}(0) = \frac{1}{6} \left[ \tilde{y}(0) + 4\tilde{y}(\frac{1}{2}) + \tilde{y}(1) \right] + 1 \\\\ \tilde{y}(\frac{1}{2}) = \frac{1}{6} \left[ \tilde{y}(0) + \tilde{y}(\frac{1}{2}) - \frac{1}{2}\tilde{y}(1) \right] - \frac{1}{2} \\\\ \tilde{y}(1) = \frac{1}{6} \left[ \tilde{y}(0) - 2\tilde{y}(\frac{1}{2}) - 2\tilde{y}(1) \right] - 2 \end{cases}$

解得 $\begin{cases} \tilde{y}(0) = \frac{2}{3} \\ \tilde{y}(\frac{1}{2}) = -\frac{1}{3} \\ \tilde{y}(1) = -\frac{4}{3} \end{cases}$，带入$\tilde{y}(x)$得

$\tilde{y}(x) = \dfrac{2}{3} - 2x$

复化求积公式

通过增加求积节点的个数来提高计算精度是不可行的（高阶Newton-Cotes公式稳定性得不到保证）。

将区间分成若干个小区间，每个区间的结果加起来得到整个区间的求积公式。

复化梯形公式

复化辛普森公式

例

计算积分$I = \int_{0}^{\frac{\pi}{2}} \sin x \,{\rm d}x$，当使用梯形公式时，应将区间$[0, \frac{\pi}{2}]$分成多少等份，才能使误差不超过$\frac{1}{2}×10^{-3}$，若取同样的求积节点，使用复化辛普森公式时截断误差是多少。

解

• $f(x) = \sin x$
• $f'(x) = \cos x$
• $f''(x) = -\sin x$
• $f^{(3)}(x) = -\cos x$
• $f^{(4)}(x) = \sin x$

$\begin{array}{l} |R(T_n)| \\\\ \\\\ \\\\ \end{array} \begin{array}{l} = |-\dfrac{b-a}{12}h^2 f''(η)| \\\\ ≤ \dfrac{\pi}{24} \cdot (\dfrac{\pi}{2n})^2 \cdot 1 \\\\ ≤ \dfrac{1}{2}×10^{-3} \end{array}$

即$n^2 ≥ \dfrac{\pi^3}{48} × 10^3$，取$n=26$

当$n=26$时，使用复化辛普森公式有

$\begin{array}{l} |R(S_n)| \\\\ \\\\ \\\\ \end{array} \begin{array}{l} = |-\dfrac{b-a}{2880}h^4 f^{(4)}(η)| \\\\ ≤ -\dfrac{\pi}{5760} \cdot (\dfrac{\pi}{26})^4 \cdot 1 \\\\ ≤ 0.7266303 × 10^{−9} \end{array}$

龙贝格（Romberg）求积公式

略

复化Cotes公式

略

高斯型求积公式

正交多项式

使求积公式具有$2n+1$次代数精度的节点$x_0,x_1,x_2,\dots,x_n$（$n+1$个节点）称为Gauss点，此时的插值型求积公式称为Gauss型求积公式。

当且仅当$x_0,x_1,x_2,\dots,x_n$是$[a,b]$上以$ρ(x)$为权函数的$n+1$次正交多项式$P_{n+1}(x)$的零点，$x_0,x_1,x_2,\dots,x_n$为Gauss点。

• Gauss-Legendre求积公式
• Gauss-Laguerre求积公式
• Gauss-Hermite求积公式

差商型求积公式

一阶向前差商

• $O(h) = -\dfrac{h}{2}f''(x_0+θh), { \ \ \ } 0≤θ≤1$

一阶向后差商

• $O(h) = \dfrac{h}{2}f''(x_0-θh), { \ \ \ } 0≤θ≤1$

一阶中心差商

• $O(h^2) = -\dfrac{h^2}{6}f'''(x_0+θh), { \ \ \ } -1≤θ≤1$

二阶中心差商

• $O(h^2) = -\dfrac{h^2}{12}f^{(4)}(x_0+θh), { \ \ \ } -1<θ<1$