Chapt-4.数值积分
根据代数插值法，对于任意被积函数$f(x)$，都可以构造一个插值多项式$p(x)$来近似代替，两边积分$\int_a^bf(x)dx\approx\int_a^b p(x)dx$
一般地，在区间$[a,b]$上选取节点$x_i$，然后用$f(x_i)$的加权平均得到平均高度$f(\zeta)$的近似值，求积公式具有一下形式 \(\int_a^b f(x)\approx \sum\limits_{i=1}^n A_kf(x_k) \\ x_k求积节点\\ A_k求积系数，A_k的选择只和x_k有关，和f(x)无关\)
梯形/Simpson/Newton-Cotes公式
-
梯形求积公式：一次代数插值，用两点$(a,f(a)),(b,f(b))$ \(\int_a^b f(x)dx\approx \int_a^b [\frac{x-b}{a-b}f(a)+\frac{x-a}{b-a}f(b)]dx \\ = \frac{a-b}{2}[f(a)+f(b)]\)
-
Simpson求积公式：二次代数插值，用三点$(a,f(a)),(\frac{a+b}{2},f(\frac{a+b}{2})),(b,f(b))$ \(\int_a^b f(x)dx\approx \int_a^b [\frac{(x-\frac{a+b}{2})(x-b)}{(a-\frac{a+b}{2})(a-b)}f(a)\\ +\frac{(x-a)(x-b)}{(\frac{a+b}{2}-a)(\frac{a+b}{2}-b)}f(\frac{a+b}{2})\\ +\frac{(x-a)(x-\frac{a+b}{2})}{(b-a)(b-\frac{a+b}{2})}f(b)]dx \\ = \frac{b-a}{6}[f(a)+4f(\frac{a+b}{2})+f(b)]\)
-
Newton-Cotes求积公式：把积分区间$[a,b]$划分$n$等分，得到$n+1$个积分节点$x_i=a+ih,(i=0,1,\cdots,n)$，其中$h=\frac{b-a}{n}$ \(f(x)\approx p_n(x)=\sum\limits_{i=0}^{n} \frac{\omega(x)}{(x-x_i)\omega'(x_i)}f(x_i) \\ 其中\omega(x)=(x-x_0)(x-x_1)\cdots(x-x_n) \\ \omega'(x_i)=(x_i-x_0)\cdots(x_i-x_{i-1})(x_i-x_{i+1})\cdots(x_i-x_n) \\\) 求积公式表示为 \(\int_a^b f(x)dx\approx \int_a^b \sum\limits_{i=0}^{n} \frac{\omega(x)}{(x-x_i)\omega'(x_i)}f(x_i) dx \\ = \sum\limits_{i=0}^{n}[\int_a^b \frac{\omega(x)}{(x-x_i)\omega'(x_i)}f(x_i)dx] \\ = \sum\limits_{i=0}^{n}A_i f(x_i) \\ 其中 A_i=\int_a^b\frac{\omega(x)}{(x-x_i)\omega'(x_i)}dx\) 经过一堆复杂的计算，得到$A_i=(b-a)c_i^{(n)}$，其中$c_i^{(n)}$为Newton-Cotes系数，一般直接使用系数表
$n$ $c_0^{(n)}$ $c_1^{(n)}$ $c_2^{(n)}$ $c_3^{(n)}$ 1 $\frac{1}{2}$ $\frac{1}{2}$ 2 $\frac{1}{6}$ $\frac{4}{6}$ $\frac{1}{6}$ 3 $\frac{1}{8}$ $\frac{3}{8}$ $\frac{3}{8}$ $\frac{1}{8}$ 所有系数之和为1？
求积公式的代数精度
数值积分是近似方法，为了保证精度，求积公式应该对“尽可能”多的函数准确成立，代数精度就是用来描述求积公式这种“性能”的方法
-
积分近似公式具有$n$次代数精度定义：设任意一个次数不高于$n$次的代数多项式$f(x)$，对于一般求积公式，积分近似公式$\int_a^b f(x)dx\approx \sum\limits_{k=0}^m A_k f(x_k)$精确成立，而当$f(x)$是$n+1$次代数多项式时不精确成立
一般的，$n$次代数精度的求积公式$\sum\limits_{k=0}^m A_kf(x_k)$对于$f(x)=1,x,\cdots,x^n$都准确成立 \(\left\{ \begin{array}{**lr**} \sum\limits_{k=0}^m A_k x_k^0=\int_a^b xdx=b-a \\ \sum\limits_{k=0}^m A_k x_k^1=\int_a^b x^2dx=\frac{1}{2}(b^2-a^2) \\ \vdots \\ \sum\limits_{k=0}^m A_k x_k^n=\int_a^b x^ndx=\frac{1}{n+1}(b^{n+1}-a^{n+1}) \end{array} \right.\)
-
定理：梯形求积公式具有1次代数精度
定理：Simpson求积公式具有3次代数精度
定理：Newton-Cotes公式至少具有n次代数精度，当n为偶数时，代数精度至少为n+1次 \(代数插值多项式的余项公式：\\ R_n(x)=f(x)-p_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega(x)\\ 设n是偶数，f(x)=\sum\limits_{k=0}^{n+1}b_kx^k \\ \int_a^b R_n(x)dx=\int_a^b \frac{f^{(n+1)}(\xi)}{(n+1)!}\omega(x) dx \\ = \int_a^b \frac{b_{n+1}(n+1)!}{(n+1)!}\omega(x) dx = b_{n+1}\int_a^b \omega(x)dx=0 \\ \because \omega(x)=(x-x_0)\cdots(x-x_{n})关于(x_{\frac{n}{2}},0)中心对称\)
-
加权积分中值定理：如果$f(x)$在$[a,b]$上连续，$g$在$[a,b]$可积并且$\forall x\in[a,b],g(x)\ge 0$（不变号），则： \(\int_a^b f(x)g(x)dx=f(\xi)\int_a^b g(x)dx，\xi\in[a,b]\)
梯形和Simpson求积公式误差估计
-
定理：若$f(x)\in C^2[a,b]$，梯形公式误差估计为： \(\int_a^b f(x)dx-\frac{b-a}{2}[f(a)+f(b)]=-\frac{(b-a)^3}{12}f''(\eta) \\ 其中a\le \eta \le b \\ \int_a^b R(x)dx=-\int_a^b \frac{f''(\xi)}{2}(a-x)(x-b)dx \\ \xi 取值和x有关\\ 令h(x)=\frac{f''(\xi)}{2}，g(x)=(a-x)(x-b) \\ \because g(x)\ge0,x\in[a,b] \\ \therefore 使用积分中值定理 \int_a^b R(x)dx= -\frac{f''(\eta)}{2}\int_a^b (a-x)(x-b)dx\)
-
定理：若$f(x)\in C^4[a,b]$，Simpson公式误差估计 \(\int_a^b f(x)dx-\frac{b-a}{6}[f(a)+f(\frac{a+b}{2})+f(b)] =-\frac{(b-a)^5}{2880}f^{(4)}(\eta) \\ \int_a^b R(x)dx=\int_a^b \frac{f^{(4)} (\xi)}{4!}(x-x_0)(x-x_1)^2(x-x_2)dx \\ 再次使用积分中值定理\)
复合求积公式
从几何意义上观察，梯形和Simpson求积公式会产生较大的误差；对Newton-Cotes公式，n过大时稳定性得不到保证
-
复合梯形求积公式定义：把积分区间$[a,b]$划分$n$等分小区间，共$n+1$个求积节点，其中$x_k=a+kh,h=\frac{b-a}{n}$，复合梯形求积公式$T_n$表示为： \(\int_a^b f(x)dx=\int_{x_0}^{x_1}f(x)dx+\int_{x_1}^{x_2}f(x)dx +\cdots +\int_{x_{n-1}}^{x_n}f(x)dx \\ = \sum\limits_{k=0}^{n-1}\int_{x_k}^{x_{k+1}}f(x)dx \approx \sum\limits_{k=0}^{n-1}\frac{x_{k+1}-x_{k}}{2} [f(x_{k+1})-f(x_{k})] \\ =\frac{h}{2}[f(a)+f(b)+2\sum\limits_{k=1}^{n-1}f(a+kh)]\)
-
复合梯形求积公式误差
定理：若$f(x)\in C^2[a,b]$，则对复合梯形求积公式$T_n$有： \(\int_a^bf(x)dx-T_n=-\frac{b-a}{12}h^2f''(\eta)\) 应用：使用误差估计式判断$n$应该取多大才能满足需要的精度要求
定理：积分近似值$T_n,T_{2n}$有如下关系 \(T_{2n}=\frac{1}{2}[T_n+h\sum\limits_{k=1}^n f(a+(2k-1)\frac{b-a}{2n})]\)
-
复合Simpson求积公式定义：把积分区间$[a,b]$划分$n$等分，$n$为偶数，设$n=2m$，共有$2m+1$个求积节点，其中$x_k=a+kh$，$h=\frac{b-a}{n}$ \(\int_{x_{2(k-1)}}^{x_{2k}}f(x)dx=\frac{h}{3} [f(x_{(2k-2)})+f(x_{(2k-1)})+f(x_{(2k)})] \\ S_n=\int_a^b f(x)dx=\frac{h}{3}[f(a)+f(b) +4\sum\limits_{k=1}^mf(x_{2k-1})+2\sum\limits_{k=1}^mf(x_{2k})]\)
-
复合Simpson求积公式误差
定理：若$f(x)\in C^4[a,b]$，则对复合Simposon求积公式$S_n$有 \(\int_a^bf(x)dx-S_n=-\frac{b-a}{2880}(2h)^4f^{(4)}(\eta)\)
自动选取步长梯形法
复合梯形公式在求积分之前必须选定$n$，可以通过误差估计公式确定，也可以使用自动选取步长梯形法 \(\int_a^b f(x)dx-T_n=-\frac{b-a}{12}h^2f''(\eta_n) \\ \int_a^b f(x)dx-T_{2n}=-\frac{b-a}{12}(\frac{h}{2})^2f''(\eta_{2n}) \\ T_{2n}-T_n=-\frac{b-a}{12}(\frac{h}{2})^2(4f''(\eta_n)-f''(\eta_{2n})) \\ \because 设n充分大的时候，f''(\eta_n)\approx f''(\eta_{2n}) \\ \therefore T_{2n}-T_n \approx-\frac{b-a}{12}(\frac{h}{2})^2f''(\eta_{2n})\times3 =(\int_a^bf(x)dx-T_{2n})\times3 \\ \therefore 若要求\vert \int_a^b f(x)dx-T_{2n}\vert <\varepsilon \\ 只需要 \vert T_{2n}-T_{n}\vert <3\varepsilon\) 迭代计算$T_{2n}$时使用前面$T_{2n}$和$T_n$关系的定理$T_{2n}=\frac{1}{2}[T_n+h\sum\limits_{k=1}^n f(a+(2k-1)\frac{b-a}{2n})]$
Richardson外推法
梯形法的递推化：复合求积方法可以提高求积精度，实际计算时如果精度不够，通过将步长逐次分半，提高求积公式的精度，以复合梯形积分为例，导出递推公式 \(T_{2n}=\frac{1}{2}[T_n+h\sum\limits_{k=1}^n f(a+(2k-1)\frac{b-a}{2n})]\) 外推技巧：从梯形公式触发，将区间$[a,b]$逐次二分可提高求积公式精度，$n$等分时，有 \(I-T_n=-\frac{b-a}{12}h^2f''(\eta) \\ 若记T_n=T(h),则T_{2n}=T(\frac{h}{2})\) 可以证明梯形公式的余项可以泰勒展开成级数形式，得到如下定理：设$f(x)\in C^\infty[a,b]$，则有 \(T(h)=I+a_1h^2+a_2h^4+\cdots+a_lh^{2l}+\cdots \\ 系数a_l和h无关\) Richardson外推算法定义：通过提升计算$I$的近似值的误差阶，提高数值计算精度的方法；外推法是数值分析的重要技巧，只要真值和近似值的误差项可以展开成$h$的幂级数，可以使用外推法，提高精度，比如 \(T(h)=I+a_1h^2+a_2h^4+\cdots+a_lh^{2l}+\cdots \\ T(\frac{h}{2})=I+a_1\frac{h^2}{4}+a_2\frac{h^4}{16} +\cdots+a_l(\frac{h}{2})^{2l}+\cdots \\ 从T_n和T_{2n}这两个求积公式构造新的求积公式 \\ S(h)=\frac{4T(\frac{h}{2})-T(h)}{3}=I+\beta_1h^4+\beta_2h^6+\cdots \\ 系数\beta_l 和h无关\) 可以看到，$T(n)$和$T({\frac{n}{2}})$的误差阶为$O(h^2)$，$S(h)$的误差阶为$O(h^4)$
同理，进一步地用$S(n),S(\frac{n}{2})$代替$T(n),T({\frac{n}{2}})$，重复上述步骤得到新的求积公式$C(h)$，误差阶为$O(h^8)$ \(C(h)=\frac{16S(\frac{h}{2}-S(h))}{15} \\ C(h)=I+r_1h^6+r_2h^8+\cdots\)
Romberg求积法
-
Richadson外推加速方法：$T_m^{(k)}$表示序列${T_0^{(k)}}$的第$m$次加速值
-
Romberg求积公式： \(T_m^{(k)}=\frac{4^m}{4^m-1}T_{m-1}^{(k+1)} -\frac{1}{4^m-1}T_{m-1}^{(k)}\) 算法计算过程：
- 取$k=0,h=b-a$，求$T_0^{(0)}=\frac{h}{2}[f(a)+f(b)]$
- 求梯形值$T_0(\frac{b-a}{2^k})$，按递推公式计算$T_0^{(k)}$
- 求加速值，当 $\vert T_k^{(0)}-T_{k-1}^{(0)}\vert <\varepsilon$ ，终止计算
Gauss 型求积公式
对一般的机械求积公式，含有$2n+2$个待定参数，$x_k,A_k$，当$x_k$为等距节点的插值求积公式的代数精度至少为$n$次 \(\int_a^b f(x)dx\approx \sum\limits_{k=0}^n A_k f(x_k)\) 如果适当选取$x_k$，有可能使求积公式具有$2n+1$次代数精度！
设带权积分$I=\int_a^bf(x)\rho(x)dx$，其中$\rho(x)$为权函数，其求积公式为： \(\int_a^bf(x)\rho(x)dx\approx \sum\limits_{k=0}^{n}A_kf(x_k)\)
- Gauss求积公式定义：对于节点$x_k,(k=1,2,\cdots,n)$，求积公式具有$2n+1$次代数精度，称为高斯型求积公式，选取的求积节点称为高斯点
高斯点的求法
-
选定权函数$\rho(x)$
-
指定$q(x)$分别为$1,x,x^2,\cdots,x^{n-1}$，令$\int_a^b \rho(x)q(x)w_n(x)dx=0$
得到$x_i,i=1,2,\cdots,n$的方程组，求出$x_i$
-
再求出$A_k$ \(A_k=\int_a^b\frac{\rho(x)w_n(x)}{(x-x_k)w_n'(x_k)}dx\)