前言

上一篇在函数逼近中我们介绍了如何在不同的误差衡量下来逼近函数，如最佳一致逼近函数OUAP以及最小二乘展开，还介绍了利用正交多项式的零点来进行多项式插值，避免边界处的Rouge现象。这次我们来介绍数值求积的思想。

具体内容详见 《数值分析》李庆扬、王能超、易大义（华中科大出版社）第四章

敲黑板

牛顿-莱布尼兹公式的局限性

在高数中我们都学过基础的微积分知识 （微积分困难户在此 @-@），给定区间范围，给定被积函数，我们一般会利用牛顿——莱布尼兹公式，刷刷刷完成题目，可能还需要一些技巧，比如积分变量代换或者分部积分啥的，但根本的求解方法还是找到被积函数的原函数。但不可否认，高数里的定积分例子只是九牛一毛，在现实中很多函数是找不到原函数的，即不一定是初等函数，而且在定量测量实验中的函数数据也可能没有表达式，只有离散的数据点，此时也使用不了牛莱公式，GG！
回归定积分最根本的思想，无限分割，以直代曲。在每个小区间内用区间的矩形面积去代替真实的面积（以一维积分为例）。由积分中值定理，我们有∫abf(x)dx=(b−a)f(ξ)\int_a^b f(x) dx=(b-a)f(\xi)∫abf(x)dx=(b−a)f(ξ)，即区间内一定存在一点f(ξ)f(\xi)f(ξ)是区间上的平均高度。而接下来介绍的几种数值积分方法的核心思想就是尽可能准确的逼近这个平均高度。

常见的近似求法

梯形公式：以左右积分端点函数值的平均值作为平均高度：∫abf(x)dx≈(b−a)f(a)+f(b)2\int_a^b f(x) dx \approx (b-a) \frac{f(a)+f(b)}{2}∫abf(x)dx≈(b−a)2f(a)+f(b)
中距公式：以区间中间的函数值作为平均高度：∫abf(x)dx≈(b−a)f(a+b2)\int_a^b f(x) dx \approx (b-a) f(\frac{a+b}{2})∫abf(x)dx≈(b−a)f(2a+b)
对于一般情况，即只在区间上取某些节点xkx_kxk，然后使用其对应的函数值加权求和来近似得到平均高度：

∫abf(x)dx≈∑k=0nAkf(xk)(1)\int_a^b f(x) dx \approx \sum_{k=0}^n A_k f(x_k) \tag{1}∫abf(x)dx≈k=0∑nAkf(xk)(1)

可以看到式1是将积分的计算转化为特定节点上函数值的加权计算。

代数精度

衡量一个近似求积公式的精度，作出一下定义：如果一个求积公式对于次数不超过m次的多项式均能够准确成立，但对于m+1次多项式不一定准确成立，则该求积公式具有m次代数精度。

1:∑Ak=b−ax:∑Akxk=12(b2−a2)...xm:∑Akxkm=1m+1(bm+1−am+1)\begin{aligned} 1 &: \sum A_k=b-a \\ x &: \sum A_k x_k=\frac{1}{2}(b^2-a^2)\\ &... \\ x^m &: \sum A_k x_k^m=\frac{1}{m+1}(b^{m+1}-a^{m+1}) \end{aligned} 1xxm:∑Ak=b−a:∑Akxk=21(b2−a2)...:∑Akxkm=m+11(bm+1−am+1)

进一步理解，如果求积公式对于m次多项式准确成立，则必定具有m次代数精度，因为m次多项式可以转化成所有次数小于m次的多项式的线性组合。
验证梯形公式只有1次代数精度：

1:∫ab1dx=b−a=(b−a)1+12x:∫abxdx=12(b2−a2)=(b−a)a+b2x2:∫abx2dx=13(b3−a3)≠(b−a)a2+b22\begin{aligned} 1 &: \int_a^b 1 dx=b-a=(b-a) \frac{1+1}{2}\\ x &: \int_a^b x dx=\frac{1}{2}(b^2-a^2)=(b-a) \frac{a+b}{2}\\ x^2 &: \int_a^b x^2 dx=\frac{1}{3}(b^{3}-a^{3}) \neq (b-a) \frac{a^2+b^2}{2}\\ \end{aligned} 1xx2:∫ab1dx=b−a=(b−a)21+1:∫abxdx=21(b2−a2)=(b−a)2a+b:∫abx2dx=31(b3−a3)=(b−a)2a2+b2

插值型求积公式

如果给定一系列节点以及对应的函数值，求一个代数精度尽可能高的求积公式，自然而然可以想到利用多项式插值的方法，因为n次插值多项式对于次数不大于n的多项式其插值余项是0。
对于给定节点a≤x0<x1<...<xn≤ba \leq x_0 < x1 <... < x_n \leq ba≤x0<x1<...<xn≤b，作插值函数（以Lagrange为例）Ln(x)L_n(x)Ln(x)，即：

Pn(x)=∑k=0nf(xk)lk(x)P_n(x)=\sum_{k=0}^n f(x_k)l_k(x)Pn(x)=∑k=0nf(xk)lk(x)

lk(x)=∏j=0,j≠knx−xjxk−xjl_k(x)=\prod_{j=0,j \neq k}^n \frac{x-x_j}{x_k-x_j}lk(x)=∏j=0,j=knxk−xjx−xj

结合式（1）则有：

∫abf(x)dx≈∑k=0nf(xk)∫ablk(x)dx=∑k=0nAkf(xk)(2)\int_a^b f(x) dx \approx \sum_{k=0}^n f(x_k) \int_a^b l_k(x) dx=\sum_{k=0}^n A_k f(x_k) \tag{2}∫abf(x)dx≈k=0∑nf(xk)∫ablk(x)dx=k=0∑nAkf(xk)(2)

其中，Ak=∫ablk(x)dxA_k=\int_a^b l_k(x) dxAk=∫ablk(x)dx

在函数多项式插值一篇中，我们可以得到插值余项为：

R[f]=∫abfn+1(ξ)(n+1)!wn+1(x)dx(3)R[f]=\int_a^b\frac{f^{n+1}(\xi)}{(n+1)!}w_{n+1}(x)dx \tag{3}R[f]=∫ab(n+1)!fn+1(ξ)wn+1(x)dx(3)

其中，wn+1(x)=∏k=0n(x−xk)w_{n+1}(x)=\prod_{k=0}^n (x-x_k)wn+1(x)=∏k=0n(x−xk)

这里可以验证，对于n次插值型求积公式，对于次数不大于n次的多项式都准确满足，即至少具有n次插值精度。引出教材定理4.1—— 形如式1的求积公式至少有n次代数精度的充要条件为其是插值型的。

Newton-Cotes公式

将等距节点插值思想应用到多项式式插值公式中就得到了Newton-Cotes公式，设将[a,b][a,b][a,b]划分为n等分，步长为h=(b−a)/nh=(b-a)/nh=(b−a)/n，选取等距节点xk=a+kh,(k=0,1,...,n)x_k=a+kh,(k=0,1,...,n)xk=a+kh,(k=0,1,...,n)构造插值求积公式：

In=(b−a)∑k=0nCk(n)f(xk)(4)I_n=(b-a)\sum_{k=0}^n C_k^{(n)} f(x_k) \tag{4}In=(b−a)k=0∑nCk(n)f(xk)(4)

其中，Ck(n)C_k^{(n)}Ck(n)为Cotes系数，由

Ak=∫ablk(x)dx=∫x0xn∏j=0,j≠knx−xjxk−xjdx=∫0n∏j=0,j≠kn(t−j)h(k−j)hhdt=(b−a)(−1)n−knk!(n−k)!∫0n∏j=0,j≠k(t−j)dtA_k=\int_a^b l_k(x) dx=\int_{x_0}^{x_n} \prod_{j = 0,j \neq k}^n \frac{x-x_j}{x_k-x_j} dx=\int_0^n \prod_{j =0,j \neq k} ^n\frac{(t-j)h}{(k-j)h} h dt=\frac{(b-a)(-1)^{n-k}}{nk!(n-k)!}\int_0^n \prod_{j =0,j \neq k}(t-j)dtAk=∫ablk(x)dx=∫x0xnj=0,j=k∏nxk−xjx−xjdx=∫0nj=0,j=k∏n(k−j)h(t−j)hhdt=nk!(n−k)!(b−a)(−1)n−k∫0nj=0,j=k∏(t−j)dt

即Cotes系数为：Ck(n)=(−1)n−knk!(n−k)!∫0n∏j=0,j≠k(t−j)dt\large C_k^{(n)}=\frac{(-1)^{n-k}}{nk!(n-k)!}\int_0^n \prod_{j =0,j \neq k}(t-j)dtCk(n)=nk!(n−k)!(−1)n−k∫0n∏j=0,j=k(t−j)dt

可以发现，Cotes系数与函数值和区间范围均无关，只与n，k有关，可以事先求解得到系数表！

n=1时，得到梯形公式

n=2时，得到常用的Simpson公式，
S=b−a6[f(a)+4f(a+b2)+f(b)](5)S=\frac{b-a}{6}[f(a)+4f(\frac{a+b}{2})+f(b)] \tag{5}S=6b−a[f(a)+4f(2a+b)+f(b)](5)

可以验证得到偶次n阶Cotes公式具有n+1次代数精度：即Simpson公式具有3次代数精度。

当n≥8n \geq 8n≥8，Cotes系数存在正负，此时稳定性得不到保证，所以一般只适用于低阶近似，其实实际使用的很少。

几种低阶Cotes公式的余项 (h=(b−a)/nh=(b-a)/nh=(b−a)/n)：
- 一阶：R[f]=∫abf′′(ξx)2!(x−a)(x−b)dx=−112h3f′′(ξ),ξ∈[a,b]R[f]=\int_a^b \frac{f^{''}(\xi_x)}{2!}(x-a)(x-b)dx=-\frac{1}{12}h^3f^{''}(\xi), \xi \in [a,b]R[f]=∫ab2!f′′(ξx)(x−a)(x−b)dx=−121h3f′′(ξ),ξ∈[a,b]
- 二阶：R[f]=∫abf4(ξx)4!(x−a)(x−c)2(x−b)dx=−190h5f4(ξ),ξ∈[a,b]R[f]=\int_a^b \frac{f^{4}(\xi_x)}{4!}(x-a)(x-c)^2(x-b)dx=-\frac{1}{90}h^5f^{4}(\xi), \xi \in [a,b]R[f]=∫ab4!f4(ξx)(x−a)(x−c)2(x−b)dx=−901h5f4(ξ),ξ∈[a,b]

复化求积

使用Newton插值法会在边界处引起Runge现象，在阶数更高时更加明显，产生不稳定的拟合结果。 为了改善求积精度，一般会采用复化求积的方法，即在每个小区间上计算积分，再将各个局部积分加和得到总积分的近似值。如果在Newton-Cotes条件下使用，则得到分段低次合成的Newton-Cotes复化求积公式！
设将区间[a,b][a,b][a,b]分为n等份，其中间隔h=b−anh=\frac{b-a}{n}h=nb−a，xk=a+kh(k=0,1,...,n)x_k=a+kh(k=0,1,...,n)xk=a+kh(k=0,1,...,n)，在每个[xk−1,xk][x_{k-1},x_k][xk−1,xk]上使用梯形公式，得到：
∫xk−1xkf(x)dx≈xk−xk−12[f(xk−1)+f(xk)]∫abf(x)dx≈h2[f(a)+2∑k=1n−1f(xk)+f(b)],(k=1,2,...,n)\int_{x_{k-1}}^{x_k} f(x)dx \approx \frac{x_k-x_{k-1}}{2}[f(x_{k-1})+f(x_k)]\\ \int_a^b f(x)dx \approx \frac{h}{2}[f(a)+2\sum_{k=1}^{n-1} f(x_k)+f(b)],(k=1,2,...,n) ∫xk−1xkf(x)dx≈2xk−xk−1[f(xk−1)+f(xk)]∫abf(x)dx≈2h[f(a)+2k=1∑n−1f(xk)+f(b)],(k=1,2,...,n)

可以看到除了左右端点值，内部的节点值都用了两次。结合中值定理，可得到其余项形式：

R[f]=∑k=1n[−112h3f′′(ξk)]=−112h2(b−a)∑k=1nf′′(ξk)n=−112h2(b−a)f′′(ξ),ξ∈(a,b)R[f]=\sum_{k=1}^n[ -\frac{1}{12}h^3 f^{''}(\xi_k)]= -\frac{1}{12}h^2 (b-a) \frac{\sum_{k=1}^n f^{''}(\xi_k)}{n}=-\frac{1}{12}h^2 (b-a) f^{''}(\xi), \xi \in (a,b)R[f]=∑k=1n[−121h3f′′(ξk)]=−121h2(b−a)n∑k=1nf′′(ξk)=−121h2(b−a)f′′(ξ),ξ∈(a,b)
如果额外的在每个区间[xk−1,xk][x_{k-1},x_k][xk−1,xk]中取节点xk+1/2x_{k+1/2}xk+1/2，可得到复化的Simpson求积公式：

Sn=∑k=0n−1h6[f(xk)+4f(xk+1/2)+f(xk+1)]=h6[f(a)+4∑k=0n−1f(xk+1/2)+2∑k=1n−1f(xk)+f(b)]S_n=\sum_{k=0}^{n-1} \frac{h}{6} [f(x_k)+4f(x_{k+1/2})+f(x_{k+1})]= \frac{h}{6}[f(a)+4 \sum_{k=0}^{n-1}f(x_{k+1/2})+2 \sum_{k=1}^{n-1}f(x_{k})+f(b)]Sn=k=0∑n−16h[f(xk)+4f(xk+1/2)+f(xk+1)]=6h[f(a)+4k=0∑n−1f(xk+1/2)+2k=1∑n−1f(xk)+f(b)]

同理可以用相同的方法得到Simposon复化求积公式的余项：
R[S]=∑k=0n−1[−190(h2)5f(4)(ξx)]=−b−a180(h2)4∑k=1nf(4)(ξx)n=−b−a180(h2)4f(4)(ξ),ξ∈(a,b)R[S]=\sum_{k=0}^{n-1}[-\frac{1}{90}(\frac{h}{2})^5 f^{(4)}(\xi_x)]=-\frac{b-a}{180}(\frac{h}{2})^4\frac{\sum_{k=1}^n f^{(4)}(\xi_x)}{n}=-\frac{b-a}{180}(\frac{h}{2})^4 f^{(4)}(\xi), \xi \in (a,b)R[S]=k=0∑n−1[−901(2h)5f(4)(ξx)]=−180b−a(2h)4n∑k=1nf(4)(ξx)=−180b−a(2h)4f(4)(ξ),ξ∈(a,b)
定义求积公式的收敛阶数，设有复化求积公式InI_nIn，当h→0h \rightarrow 0h→0，有：I−Inhp→C(C≠0,定数)\frac{I-I_n}{h^p} \rightarrow C (C \neq 0, 定数)hpI−In→C(C=0,定数)，则称InI_nIn是p阶收敛的。由以上定义，可知复化的梯形求积公式，Simpson公式和N-C公式具有2，4，6阶收敛精度。

Gauss 求积公式

n+1个节点的插值型求积公式至少具有n阶代数精度，能不能继续提高精度的同时保证稳定性和收敛性呢？Gauss给我们带来了从整体区间上优化的方法：

对于∫baf(x)ddx≈∑k=0nAkf(xk)\int_b^a f(x)ddx \approx \sum_{k=0}^n A_k f(x_k)∫baf(x)ddx≈∑k=0nAkf(xk)，设有2n+2个未知参数Ak,xk,(k=0,1,...,n)A_k,x_k,(k=0,1,...,n)Ak,xk,(k=0,1,...,n)，选择合适的参数，则有可能使求积公式具有2n+1次代数精度。
人们将使用n+1个点且有2n+1次代数精度的求积公式称为Gauss求积公式，其中的节点称为Gauss节点。前面我们提到，n+1个节点的插值型求积公式至少具有n阶代数精度，所以Gauss求积公式也势必是插值型的，那怎么确定Gauss节点呢？下面我们给出定理：对于插值型求积公式的节点xk,(k=0,1,...,n)x_k,(k=0,1,...,n)xk,(k=0,1,...,n)是Gauss点的充要条件为，以这些节点为零点的多项式w(x)=∏k=0n(x−xk)w(x)=\prod _{k=0}^n (x-x_k)w(x)=∏k=0n(x−xk)与任意次数不超过n的多项式P(x)P(x)P(x)都正交，即∫baP(x)w(x)dx=0\int_b^a P(x)w(x) dx=0∫baP(x)w(x)dx=0。

证明如下：
- 证必要性：设P(x)P(x)P(x)是次数不超过n的多项式，则P(x)w(x)P(x)w(x)P(x)w(x)的次数不超过2n+1，则求积公式对其准确成立，有∫baP(x)w(x)dx=∑k=0nAkP(xk)w(xk)=0，其中w(xk)=0,(k=0,1,...,n)\int_b^a P(x)w(x)dx=\sum_{k=0}^n A_k P(x_k)w(x_k)=0，其中w(x_k)=0,(k=0,1,...,n)∫baP(x)w(x)dx=k=0∑nAkP(xk)w(xk)=0，其中w(xk)=0,(k=0,1,...,n)
- 证充分性：给定次数不超过2n+1的多项式f(x)f(x)f(x)，由多项式除法，f(x)f(x)f(x)可表示为：f(x)=P(x)w(x)+Q(x)f(x)=P(x)w(x)+Q(x)f(x)=P(x)w(x)+Q(x)其中，P(x),Q(x)P(x),Q(x)P(x),Q(x)都是次数不超过n的多项式
  
  由∫baP(x)w(x)dx=0\int_b^a P(x)w(x) dx=0∫baP(x)w(x)dx=0，可得
  
  ∫baf(x)dx=∫baQ(x)dx\int_b^a f(x)dx =\int_b^a Q(x) dx∫baf(x)dx=∫baQ(x)dx
  
  对于节点x0,x1,...,xkx_0,x_1,...,x_kx0,x1,...,xk,构建插值型积分公式，∫baQ(x)dx=∑k=0nAkQ(xk)\int_b^a Q(x)dx = \sum_{k=0}^n A_k Q(x_k)∫baQ(x)dx=∑k=0nAkQ(xk)对Q(x)Q(x)Q(x)准确成立。
  
  又由于w(xk)=0w(x_k)=0w(xk)=0，可得Q(xk)=f(xk)Q(x_k)=f(x_k)Q(xk)=f(xk)，从而有
  ∫baQ(x)dx=∑k=0nAkf(xk)\int_b^a Q(x)dx = \sum_{k=0}^n A_k f(x_k)∫baQ(x)dx=k=0∑nAkf(xk)，即是∫baf(x)dx=∑k=0nAkf(xk)\int_b^a f(x)dx = \sum_{k=0}^n A_k f(x_k)∫baf(x)dx=∑k=0nAkf(xk)，即以x0,x1,...,xkx_0,x_1,...,x_kx0,x1,...,xk,构建的插值型积分公式对所有次数不超过2n+1的多项式均准确成立，x0,x1,...,xkx_0,x_1,...,x_kx0,x1,...,xk是Gauss点，证毕。

基于正交多项式的Gauss求积公式

由Gauss点的性质，我们很容易联想到正交多项式，考虑区间[−1,1],ρ=1[-1,1], \rho=1[−1,1],ρ=1，基于Legendre多项式可得到∫−11f(x)dx≈∑k=0nAkf(xk)\int_{-1}^1 f(x)dx \approx \sum_{k=0}^n A_k f(x_k)∫−11f(x)dx≈k=0∑nAkf(xk)

其中，插值节点即为Legendre多项式的零点。
同理，若取区间[−1,1],ρ=11−x2\large [-1,1], \rho=\frac{1}{\sqrt{1-x^2}}[−1,1],ρ=1−x21，可得到基于Chebyshev的Gauss求积公式，特别地，其加权系数Ak=πn+1\large A_k=\frac{\pi}{n+1}Ak=n+1π，节点为cos⁡2k+12n+2π,(k=0,1,...,n)\large \cos{\frac{2k+1}{2n+2}\pi},(k=0,1,...,n)cos2n+22k+1π,(k=0,1,...,n)。

数值分析——数值积分（Newton-Cotes、复化求积、Gauss求积、正交多项式Gauss）相关推荐

自适应复化辛普森公式求积算法（C语言实现）
自适应复化辛普森公式求积算法(C语言实现) 利用复化辛普森公式求积分自适应步骤基于C语言实现的代码利用复化辛普森公式求积分自适应步骤 h为步长,a为积分下限,b为积分上限,f为积分函数,n为划分的 ...
2021-01-07 matlab数值分析数值积分与数值微分复合梯形公式复合Simpson公式
matlab数值分析数值积分与数值微分 1 复合梯形公式 function I=ftrapz(f,a,b,n) format long %显示15位双精度 h=(b-a)/n; x=linspace ...
龙贝格求积分算法例题_数值分析实习作业之龙贝格求积
北京化工大学数值分析,龙贝尔求积, 实习题五:用龙贝格求积公式法计算习题五第一题积分,是结果尽可能准确(eps = 0.000001) 主程序: %习题五:用龙贝格求积公式计算下列积分 clear c ...
基于复化梯度求积的求积步长自适应matlab实现
一.基本思想把区间逐次二分,反复利用复化求积公式进行计算,直至二分前后的两次积分近似值之差符合精度要求为止.采用误差的"事后估计"方法进行步长的自动选择. 以复化梯形公式为例,设 ...
复化科特斯公式matlab_matlab实现复化NewtonCotes公式求积分的程序应用和代码
命令窗口输入被积函数 2 1 2 0 t t e dt  . 输入应为:. Step2 :执行函数.输入形式为 mymulNewtonCotes(ft,a,b,m,n) : 其中 ft -被积函数, ...
C++实现复化辛普森公式求积分算法
1. 算法原理简介步1 将积分区间 [a,b] 分成 n 等分,分点为 xk=a+kh(k=0,1,⋯,n),其中 h=(b-a)/n. 步2 记区间 [xk,x(k+1)] 的中点为 x(k+1/ ...
变步长复化梯形法求积
#include <iostream> using namespace std; void main() { float a,b,e,result; float f(float x); f ...
数值积分的python实现——NewtonCotes、复化求积、Romberg、richardson递推
数值积分的python实现--NewtonCotes,复化求积,Romberg,richardson递推 1. 机械求积 2. Newton-Cotes公式 3. 复化求积方法 3.1 复化梯形公式 ...
角速度求积分能得到欧拉角吗_数值计算方法第六章数值积分和数值微分
写在章前:积分与微分的计算,是具有广泛应用的古典问题. 然而,在微积分教材中,只对简单的或特殊的情况,提供了函数的积分或微分的解析表达式,在理论上可以使用Newton-Leibniz公式计算.但对很多 ...

数值分析——数值积分（Newton-Cotes、复化求积、Gauss求积、正交多项式Gauss）

数值分析-数值积分

前言

敲黑板