# 机械求积法

梯形公式：

$$\int_a^b f(x)dx\approx \frac12(b-a)[f(a)+f(b)]$$

中矩形公式：

$$\int_a^b f(x)dx\approx (b-a)f(\frac{a+b}{2})$$

辛普森公式：

$$\int_a^bf(x)dx\approx\frac{1}{6}(b-a)[f(a)+4f(\frac{a+b}{2})+f(b)]$$

# 代数精度

使求积公式准确成立的多项式的次数，可以作为一种衡量求积公式精确度程度的标准

若 $\int_a^b f(X)dx\approx \sum_{i=0}^nA_if(x_i)$ 对于一切次数小于等于 m 的多项式是准确的，而对于次数为 $m+1$ 的某一多项式并不准确，则称该求积公式具有 $m$ 次代数精度

推论：若 $\int_a^b f(X)dx\approx \sum_{i=0}^nA_if(x_i)$ 具有 $m$ 次代数精度，则这个公式对于一切次数小于等于 $m$ 的多项式是准确的

证明

若 $\int_a^b f(X)dx\approx \sum_{i=0}^nA_if(x_i)$ 对函数 $f(x)$ 和 $g(x)$ 准确成立，则它对于函数 $\alpha f(x)+\beta g(x)$ 准确成立

$$\begin{aligned} \int_a^b[\alpha f(x)+\beta g(x)]dx&=\int_a^b\alpha f(x)dx+\int_a^b \beta g(x)dx\\ &=\alpha \sum_{i=0}^nA_if(x_i)+\beta \sum_{i=0}^nA_ig(x_i)=\sum_{i=0}^n A_i[\alpha f(x_i)+\beta g(x_i)] \end{aligned}$$

• 梯形公式具有一次代数精度
• 辛普森具有三次代数精度

对于任意给定的 $n+1$ 个互异的节点 $a=x_0<x_1<...<x_n=b$ 总存在系数 $A_0,A_1,...,A_n$，使得 $\int_a^b f(X)dx\approx \sum_{i=0}^nA_if(x_i)$ 的代数精度最低为 $n$

证明

设以 $x_0,x_1,...,x_n$ 为插值节点的 $f(x)$ 的 $Lagrange$ 插值多项式为

$$p_n(x)=\sum_{i=0}^nf(x_i)l_i(x)$$

则

$$f(x)=p_n(x)+\frac{f^{(n+1)}(\xi)}{(n+1)!}w(x)$$

所以

$$\int_a^bf(x)dx=\sum_{i=0}^nf(x_i)\int_a^bl_i(x)dx+\frac{1}{(n+1)!}\int_a^bf^{(n+1)}(\xi)w(x)dx$$

记

$$R(f)=\frac{1}{(n+1)!}\int_a^bf^{(n+1)}(\xi)w(x)dx$$

对于任何一个次数不高于 $n$ 次的多项式 $f(x)$，$R(f)=0$

# 求积公式的构造方法

# 利用代数精度作为标准构造求积公式

因为如果代数精度为 $m$，那么对于一切次数小于等于 $m$ 的多项式是准确的，所以可以带入 $1,x,x^2,...$ 进行解方程求解

# 插值法构造求积公式

用 $p_n(x)$ 代替 $f(x)$ 从而的获得求积公式

$$\int_a^bf(x)dx\approx \int_a^bp_n(x)dx$$

$$\int_a^bp_n(x)dx=\int_a^b [\sum_{i=0}^nf(x_i)l_i(x)]dx=\sum_{i=0}^nf(x_i)\int_a^bl_i(x)dx$$

显然，求积系数为

$$A_i=\int_a^bl_i(x)dx$$

定义：求积公式的求积系数为 $A_i=\int_a^bl_i(x)dx$ ，则此求积公式是插值型的

定理：求积公式 $\int_a^bf(x)dx\approx \sum_{i=0}^n A_if(x_i)$ 至少具有 $n$ 次代数精度的充要条件是它是插值型的

证明

$\Large\Rightarrow$

因为 $\int_a^bf(x)dx\approx \sum\limits_{i=0}^n A_if(x_i)$ 至少具有 $n$ 次精度，所以它对所有次数小于等于 $n$ 次的多项式准确成立

则以 $x_0,x_1,..,x_n$ 为插值节点的 $Lagrange$ 插值基函数 $l_i(x)$ 是个 $n$ 次多项式，准确成立

$$\int_a^bf(x)dx= \sum\limits_{i=0}^n A_il_i(x)$$

因为

$$l_i(x)=\left\{\begin{matrix} 1 & i=j\\ 0 & i\neq j \end{matrix}\right.$$

则

$$\int_a^bf(x)dx= A_i$$

$\Large\Leftarrow$

因为对于任意次数小于等于 $n$ 的多项式 $f(x)$ ，其 $n$ 次插值多项式 $p(x)$ 就是其自身（因为 $f^{(n+1)}(x)=0$）

所以插值型的求积公式对于一切次数小于等于 $n$ 的多项式准确成立

# Newton-Cotes 求积公式

如果将积分区间 $[a,b]$ 分成 $n$ 等份，其求积节点 $x_i$ 为 $x_i=a+ih,i=0,1,2,...,n$，$h=\frac{b-a}{n}$ 表示各等分小区间的宽度

$$\begin{aligned} \int_a^b f(x)\,dx &= \int_a^b p_n(x)\,dx = \sum_{i=0}^n \left[f(x_i) \int_a^b l_i(x)\,dx\right] \\ &= (b-a)\sum_{i=0}^n \left[\frac{\int_a^b l_i(x)\,dx}{b-a} \cdot f(x_i)\right] \end{aligned}$$

科特斯系数：

$$C_i=\frac{\int_a^bl_i(x)dx}{b-a}$$

因为 $x=a+th$

$$x=a\rightarrow t=0\\ x=b\rightarrow t=n\\ dx=hdt\\ b-a=nh$$

所以 $C_i$ 可以化简为

$$\begin{aligned} C_i&=\frac{1}{nh}\int_0^n [\prod _{j=0,j\neq i}^n\frac{x-x_j}{x_i-x_j}]dx\\ &=\frac{1}{nh}\int_0^n [\prod _{j=0,j\neq i}^n\frac{(a+th)-(a+jh)}{(a+ih)-(a+jh)}]dx\\ &=\frac{1}{n} \int_0^n \prod _{j=0,j\neq i}^n\frac{t-j}{i-j}dt\\ &=\frac{(-1)^{n-i}}{n\times i!\times(n-i)!}\int_0^n[\prod _{j=0,j\neq i}^n(t-h)]dt \end{aligned}$$

# 牛顿 - 科特斯公式的稳定性

考虑计算 $f(x_i)$ 时产生的舍入误差 $\varepsilon_i$

$$\int_a^b f(x)dx\approx (b-a)\sum_{i=0}^n C_i[f(x_i)+\varepsilon_i]$$

则由舍入误差引起的积分误差为：

$$\begin{aligned} \varepsilon &= (b-a)\sum_{i=0}^n C_i[f(x_i)+\varepsilon_i] - (b-a)\sum_{i=0}^n C_if(x_i) \\ &= (b-a)\sum_{i=0}^n C_i\varepsilon_i \end{aligned}$$

记 $\varepsilon_{max}=\max|\varepsilon_i|$

$$\begin{aligned} |\varepsilon|=&|(b-a)\sum_{i=0}^nC_i\varepsilon_i|\leq (b-a)\sum_{i=0}^n|C_i\varepsilon_i|=(b-a)\sum_{i=0}^n|C_i||\varepsilon_i|\\ &\leq (b-a)\varepsilon_{max}\sum_{i=0}^n|C_i| \end{aligned}$$

当 $n\leq 7$ 时， $C_i$ 皆为正，故：$\sum_{i=0}^n|C_i|=\sum_{i=0}^n C_i=0$

$$|\varepsilon|\leq (b-a)\varepsilon_{max}$$

$\rightarrow$ 由舍入误差引起的积分误差有上界

当 $n\geq 8$ 时， $C_i$ 出现负数， 故：$\sum_{i=0}^n|C_i|$ 随着 $n$ 的增大而不断增大

$\rightarrow$ 由舍入误差引起的积分误差无上界

** 总结：** 高阶的牛顿 - 科特斯求积公式不是数值稳定的，所以不宜采用

# 公式的分类

梯形公式

几何意义：用四边形梯形的面积近似替代 $y=f(x)$ 围成的曲边梯形的面积

辛普森公式

几何意义：用抛物线 $y=p_2(x)$ 围成的曲边梯形面积近似替代 $y=f(x)$ 围成的曲边梯形的面积

科特斯公式

由四次 Lagrange 插值式推导而得得求积公式

$$\int_a^b f(x)dx\approx (b-a)\times\frac{7f(x_0)+32f(x_1)+12f(x_2)+32f(x_3)+7f(x_4)}{90}$$

# 公式的余项

# 牛顿 - 科特斯公式的余项

$$R=\int_a^b\frac{f^{(n+1)}(\xi)}{(n+1)!}w(x)dx$$

令 $x=a+th$

$$\begin{aligned} R&=\int_0^n\frac{f^{(n+1)}(\xi)}{(n+1)!}(a+th-a-0h)(a+th-a-1h)...(a+th-a-nh)hdt\\ &=\frac{h^{n+2}}{(n+1)!}\int_0^nf^{(n+1)}(\xi)\prod_{j=0}^n(t-j)dt \end{aligned}$$

定理：当 $n$ 为偶数时，$n$ 阶 Newton-Cotes 求积公式得代数精度至少是 $n+1$

证明

记 $n+1$ 次多项式为

$$q_{n+1}(x)=\sum_{k=0}^{n+1}b_kx^k$$

其 $n+1$ 阶导数为

$$q_{n+1}(x)^{(n+1)}=b_{n+1}(n+1)!$$

则

$$\begin{aligned} R&=\frac{h^{n+2}}{(n+1)!}\int_0^nb_{n+1}(n+1)!\prod_{j=0}^n(t-j)dt\\ &=b_{n+1}h^{n+2}\int_0^n(t-0)(t-1)....(t-n)dt \end{aligned}$$

令 $n = 2k,u=t-k$

$$R=b_{n+1}h^{n+2}\int_{-k}^k (u+k)(u+k-1)...u(u-k+1)(u-k)du$$

令

$$\begin{aligned} H(u)&=(u+k)(u+k-1)...(u-k+1)(u-k)\\ H(-u)&=(-u+k)(-u+k-1)...(-u-k+1)(-u-k)\\ &=(-1)^{2k+1}H(u)=-H(u) \end{aligned}$$

所以 $H(u)$ 是奇函数 $\int_{-k}^kH(u)=0$

$$R=0$$

# 梯形公式的余项

设 $f(x)$ 在 $[a,b]$ 上由二阶连续导数

$$\begin{aligned} R_T&=\frac{h^{1+2}}{(1+1)!}\int_0^1f^{(1+1)}(\xi)(t-0)(t-1)dt \\&=\frac{h^3}{2}\int_0^1f''(\xi)t(t-1)dt\\ &=-\frac{h^3}{12}f''(\eta) \end{aligned}$$

其中 $h=b-a$

# 辛普森公式的余项

设 $f(x)$ 在 $[a,b]$ 上有四阶连续导数

因为辛普森公式具有三次代数精度，所以只需将 $f(x)$ 表示成某个三次插值多项式与其余项的和，然后再求积分即可知道辛普森公式的余项

设 $p_3(x)$ 是三次多项式且满足

$$p_3(a)=f(a),p_3(\frac{a+b}2)=f(\frac{a+b}{2}),p_3(b)=f(b)\\ p_3'(\frac{a+b}{2})=f'(\frac{a+b}{2})$$

三次 Hermite 插值多项式的截断误差为：

$$R_3(x) = \frac{f^{(4)}(\xi)}{4!}(x-a)\left(x-\frac{a+b}{2}\right)^2(x-b)$$

所以

$$\begin{aligned} R&=\int_a^bf(x)dx-\int_a^bp_3(x)dx\\ &=\int_a^bf(x)dx-\frac{1}{6}(b-a)[f(a)+4f(\frac{a+b}{2})+f(b)]\\ &=\int_a^b\frac{f^{(4)}(\xi)}{4!}(x-a)\left(x-\frac{a+b}{2}\right)^2(x-b)dx \end{aligned}$$

因为 $x\in[a,b]$，所以 $(x-a)(x-b)\leq 0$，所以 $(x-a)\left(x-\frac{a+b}{2}\right)^2(x-b)$ 不变号

由积分中值定理得

$$\begin{aligned} R &=\frac{f^{(4)}(\eta)}{4!}\int_a^b(x-a)\left(x-\frac{a+b}{2}\right)^2(x-b)dx \\&=-\frac{1}{90}(\frac{b-a}{2})^5f^{(4)}(\eta) \end{aligned}$$

# 复化求积

# 复化梯形

$$\int_a^bf(x)dx=\frac{h}{2}[f(a)+2\sum_{i=0}^{n-1}f(x_k)+f(b)]$$

# 复化辛普森公式

$$\int_a^bf(x)dx=\frac{h}{3}[f(a)+4\sum_{k=1}^mf(x_{2k-1})+2\sum_{k=1}^{m-1}f(x_{2k})+f(b)]$$

# 复化 Cotes 公式

# 复化求积的截断误差

# 复化梯形的截断误差

$$R_{T}=-\frac{1}{12}nh^3f''(\eta),\eta\in[a,b]$$

# 复化辛普森公式的截断误差

$$R_S=-\frac{(n/2)h^5}{90}f^{(4)}(\eta),\eta\in[a,b]$$

# 步长的自动选取

小区间 $[x_i,x_{i+1}]$，记小区间中点： x_{i+\frac{1}{2}}=\frac{x_i+x_{i+1}}

区间二分前、后的积分值分别记为 T_{i1},T_

$$\begin{aligned} T_{i1}&=\frac{h}{2}[f(x_i)+f(x_{i+1})]\\ T_{i2}&=\frac{h/2}{2}[f(x_i)+f(x_{i+\frac{1}{2}})]+\frac{h/2}{2}[f(x_{i+\frac{1}{2}})+f(x_{i+1})]\\ &=\frac{1}{2}T_{i1}+\frac{h}{2}f(x_{i+\frac{1}{2}}) \\T_{2n}&=\sum_{i=0}^{n-1}T_{i2}\\ &=\frac{1}{2}\sum_{i=0}^{n-1}T_{i1}+\frac{h}{2}\sum_{i=0}^{n-1}f(x_{i+\frac{1}{2}})\\ &=\frac{1}{2}T_n+\frac{h}{2}\sum_{i=0}^{n-1}f(x_{i+\frac{1}{2}}) \end{aligned}$$

# 龙贝格（Romberg）求积公式

复化梯形求积，积分区间分为 $n$ 等分时的截断误差为：

$$I-T_n=-\frac{h^2}{12}(b-a)f''(\eta)=-\frac{h^2}{12}\int_a^bf''(x)dx=-\frac{h^2}{12}[f'(b)-f'(a)]$$

其中 h=\frac{b-a}

显然

$$I-T_n\propto h^2$$

所以

$$I-T_{2n}=\frac{1}{4}(I-T_n)$$

化简可得

$$I-T_{2n}\approx\frac{1}{3}(T_{2n}-T_n)$$

叠加上误差，可得到比 $T_{2n}$ 更精确的积分值

$$\bar T=T_{2n}+\frac{1}{3}(T_{2n}-T_n)=\frac{4}{3}T_{2n}-\frac{1}{3}T_n$$

经化简得

$$\frac{4}{3} T_{2n} - \frac{1}{3} T_n = S_{2n}$$

化简过程

$$\begin{aligned} \frac{4}{3} T_{2n} &= \frac{4}{3} \cdot \frac{h/2}{2} \left[ f(a) + 2 \sum_{i=0}^{n-1} f(x_{i+\frac{1}{2}}) + 2 \sum_{i=1}^{n-1} f(x_i) + f(b) \right] \\&= \frac{h}{6} \left[ 2f(a) + 4 \sum_{i=0}^{n-1} f(x_{i+\frac{1}{2}}) + 4 \sum_{i=1}^{n-1} f(x_i) + 2f(b) \right] \\\frac{1}{3} T_n &= \frac{1}{3} \cdot \frac{h}{2} \left[ f(a) + 2 \sum_{i=1}^{n-1} f(x_i) + f(b) \right]\\ &= \frac{h}{6} \left[ f(a) + 2 \sum_{i=1}^{n-1} f(x_i) + f(b) \right]\\ \bar{T} &= \frac{4}{3} T_{2n} - \frac{1}{3} T_n\\ &= \frac{h}{6} \left[ 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) \right]\\ &= S_{2n} \end{aligned}$$

根据上述公式，我们可以总结出

$$T_{2n}=\frac{T_n}{2}+\frac{h}{2}\sum_{i=1}^{n-1}f(x_{i+\frac{1}{2}})$$

其中 $h$ 表示二分前的步长，$f(x_{i+\frac{1}{2}})$ 表示二分节点的值