数值积分 – Notes

数值积分

牛顿-柯特斯求积公式

定义
- 已知积分区间 $[a,b]$ 和被积函数 $f(x)$ 在 $[a,b]$ 上等距的 $n+1$ 个插值节点 $x_0,x_1,\dots,x_n$ 和函数值。
- $I=\displaystyle\int_a^b f(x) \mathrm dx \approx \sum_{i=0}^n \int_a^b a_i(x)f(x_i) \mathrm dx = \sum_{i=0}^n c_i f(x_i) = (b-a) \sum_{i=0}^n C_i^{(n)} f(x_i)$ 。
- 最后 $I \approx (b-a) \displaystyle\sum_{i=0}^n C_i^{(n)} f(x_i)$ 称牛顿-柯特斯求积公式.
- 定义 $C_i^{(n)}$ 为柯特斯系数，只与插值节点数和下标有关。
- 对于任意积分，只要在区间 $[a,b]$ 选取等距的 $n+1$ 个点和函数值，分别相乘系数和相加即可得近似值。

柯特斯系数表

$n$	$C_0^{(n)}$	$C_1^{(n)}$	$C_2^{(n)}$	$C_3^{(n)}$	$C_4^{(n)}$
$1$	$1/2$	$1/2$
$2$	$1/6$	$4/6$	$1/6$
$3$	$1/8$	$3/8$	$3/8$	$1/8$
$4$	$7/90$	$16/45$	$2/15$	$16/45$	$7/90$

$C_i^{(n)}$ 有对称性， $C_i^{(n)}=C_{n-i}^{(n)}$ 。
$\displaystyle \sum_{i=0}^n C_i^{(n)}=1$ 。

误差估计
- 牛顿-柯特斯求积公式由拉格朗日公式推导而来，截断误差则也可以使用插值的误差推导。
- $n$ 次插值多项式有 $R_n(x) = \dfrac{f^{(n+1)}(\xi)}{(n+1)!}\displaystyle\prod_{i=0}^n(x-x_i)$ ，则求积公式有 $R=\displaystyle\int_a^b R_n(x)\mathrm dx$ 。
- 舍入误差有 $\epsilon = (b-a)\displaystyle\sum_{i=0}^n (|C_i^{(n)}| \Delta f| + |f| \Delta C_i^{(n)})$ 。
常用形式
- 梯形公式 $n=1$
  - $I \approx \dfrac{b-a}{2} (f(a) + f(b))$
  - $R=-\dfrac{(b-a)^3}{12}f''(\xi)$
- 辛普森公式 $n=2$
  - $I \approx \dfrac{b-a}{6}\left( f(a) + 4f\left(\dfrac{a+b}{2}\right) + f(b) \right)$
  - $R=-\dfrac{(b-a)^5}{2880}f^{(4)}(\xi)=-\dfrac{h^5}{90}f^{(4)}(\xi)$
- 柯特斯公式 $n=4$
  - $I \approx \dfrac{b-a}{90} (7f(x_0) + 32f(x_1) + 12f(x_2) + 32f(x_3) + 7f(x_4))$

复化求积公式
- 定义
  - 对于积分区间 $[a,b]$ ，等距分为 $m$ 大段，每个大段分别运用求积公式，最后累加。
  - 每一大段内分为 $n$ 个小段，即每一大段的积分值使用 $n$ 次牛顿-柯特斯求积公式求解。
  - 总共的小段数记作 $M=mn$ 。
- 常用形式
  - 复化梯形公式 $n=1$
    - $M=m$
    - $R = -\dfrac{(b-a)^3}{12m^2}f''(\xi)$
  - 复化辛普森公式 $n=2$
    - $M=2m$
    - $R = -\dfrac{(b-a)^5}{2880m^4}f^{(4)}(\xi)$
  - 通用形式
    - 若 $R$ 中包含 $(b-a)^k$ ，则复化后的 $R' = m\dfrac{R}{m^k} = \dfrac{R}{m^{k-1}}$ 。
- 递推形式
  - 复化求积公式支持增量地添加节点，而已有不需要计算原来函数值。
  - 对梯形公式，设 $T_M$ 表示有 $M=m$ 个小段的复化求积结果。
  - $T_{2^k} = \dfrac{1}{2}T_{2^{k-1}} + h\displaystyle\sum_{i=1}^{2^{k-1}} f(a+(2i-1)h)$ ，后面的 $a+(2i-1)h$ 为原来 $2^{k-1}$ 个区间的中点。
- 事后估计误差
  - 已知 $T_M,T_{2M}$ ，则 $I \approx T_{2M} + \dfrac{1}{3}(T_{2M} - T_M)$ ， $\Delta = \dfrac{1}{3}(T_{2M} - T_M)$ 为估计的 $T_{2M}$ 的误差。
  - 以上公式也可以看作加上误差来修正结果，对各常用公式都有：
    - 梯形公式： $I \approx T_{2M} + \dfrac{1}{3}(T_{2M} - T_M) = \dfrac{4}{4-1}T_{2M} - \dfrac{1}{4-1}T_M$
    - 辛普森公式： $I \approx S_{2M} + \dfrac{1}{15}(S_{2M} - S_M) = \dfrac{4^2}{4^2-1}S_{2M} - \dfrac{1}{4^2-1}S_M$
    - 柯特斯公式： $I \approx C_{2M} + \dfrac{1}{63}(C_{2M} - C_M) = \dfrac{4^3}{4^3-1}C_{2M} - \dfrac{1}{4^3-1}C_M$
龙贝格法
- 推导可得 $S_M = \dfrac{4}{4-1}T_{2M} - \dfrac{1}{4-1}T_M$ ， $C_M = \dfrac{4^2}{4^2-1}S_{2M} - \dfrac{1}{4^2-1}S_M$ 。
- 每种两个积分值可以得到一个更精确的积分值。
- 定义龙贝格序列 $R_M = \dfrac{4^3}{4^3-1}C_{2M} - \dfrac{1}{4^3-1}C_M$ 。
- 理论上可以继续基于 $R_M$ 组合，但系数的开始变化不大，对精度要求更高。
- 实际计算采用如下顺序： $T_1 \to T_2 \to S_1 \to T_4 \to S_2 \to C_1 \to T_8 \to C_4 \to S_2 \to R_1 \to \cdots$ 。
高斯求积公式
- 已知积分区间 $[a,b]$ ，被积函数 $f(x)$ 。
- $I \approx \dfrac{b-a}{2} \displaystyle\sum_{i=1}^n w_i^{(n)} f(x_i)$ ， $x_i = \dfrac{a+b}{2} + \dfrac{b-a}{2}t_i^{(n)}$ ， $t_i^{(n)}$ 为 $n$ 次勒让德多项式的 $n$ 个零点。
- $t_i^{(n)},w_i^{(n)},R_n$ 都可以查表得到。
求积公式的代数精度
- 定义
  - 求积公式 $I_n = \displaystyle\sum_{i=0}^n A_i f(x_i)$ $I_{n} = i = 0 \sum n A_{i} f (x_{i})$ 的代数精度为 $m$ $m$ ，当且仅当
    - 对次数小于等于 $m$ 的多项式的截断误差 $I-I_n=0$ ，
    - 对次数大于 $m$ 的 $I-I_n\ne 0$ 。
- 计算
  - 对于已知求积公式，只要不断判断 $1,x,x^2,x^3,\dots$ 的积分值是否精确，直到不相等。
  - 对于 $A_i$ $A_{i}$ 未知的求积公式，要求确定 $A_i$ $A_{i}$ 以得到最高代数精度
    - 先用 $1,x,x^2,x^3,\dots,x^n$ 的积分值构造方程组求解 $A_i$ ，此时至少有 $n$ 次精度。
    - 从 $x^{n+1}$ 开始判断求解的 $A_i$ 是否计算精确，若准确则代数精度可以增加。
- 性质
  - $n$ $n$ 次牛顿-柯特斯求积公式（ $n+1$ $n + 1$ 个点）的代数精度：
    - 当 $n$ 为奇数时，代数精度为 $n$ 。
    - 当 $n$ 为偶数时，代数精度为 $n+1$ 。
  - 使用 $n$ 个节点的高斯求积公式的代数精度为 $2n-1$ 。

插值