Newton-Cotes求积公式（一） | 一般形式 +公式稳定性 +截断误差与代数精度

1. Newton-Cotes公式的一般形式

对于I=∫abf(x)dxI=\int_a^bf(x)dxI=∫abf(x)dx，将区间[a,b][a,b][a,b]分成n等份，步长为：
h=(b−a)nh=\frac{(b-a)}{n} h=n(b−a)
故求积节点xk=a+kh(k=0,1,2,⋯,n)x_k=a+kh(k=0,1,2,\cdots,n)xk=a+kh(k=0,1,2,⋯,n)，通过这n+1n+1n+1个节点构造一个n次代数多项式pn(x)≅f(x)p_n(x) \cong f(x)pn(x)≅f(x)。令
x=a+th(t=0,1,2,⋯,k,⋯,n)x=a+th \quad (t=0,1,2,\cdots,k,\cdots,n) x=a+th(t=0,1,2,⋯,k,⋯,n)
作积分变换
dx=hdtdx=hdt dx=hdt

x=a,t=0;x=b,t=nx=a,t=0; \quad x=b,t=n x=a,t=0;x=b,t=n

由Ak=∫ablk(x)dxA_k=\int_a^bl_k(x)dxAk=∫ablk(x)dx计算求积系数Ak(k=0,1,2,⋯,n)A_k(k=0,1,2,\cdots,n)Ak(k=0,1,2,⋯,n)，即：
Ak=∫ablk(x)dx=∫ab∏j=0,j≠knx−xjxk−xjdx=∫0n∏j=0,j≠kn(a+th)−(a+jh)(a+kh)−(a+jh)hdt=∫0n∏j=0,j≠knt−jk−jhdt=h∫0n∏j=0,j≠knt−jk−jdtA_k=\int_a^bl_k(x)dx=\int_a^b\prod^n_{j=0,j\neq k}\frac{x-x_j}{x_k-x_j}dx=\int_0^n\prod^n_{j=0,j\neq k}\frac{(a+th)-(a+jh)}{(a+kh)-(a+jh)}hdt \\ =\int_0^n \prod^n_{j=0,j\neq k}\frac{t-j}{k-j}hdt=h\int_0^n\prod_{j=0,j\neq k}^n\frac{t-j}{k-j}dt Ak=∫ablk(x)dx=∫abj=0,j=k∏nxk−xjx−xjdx=∫0nj=0,j=k∏n(a+kh)−(a+jh)(a+th)−(a+jh)hdt=∫0nj=0,j=k∏nk−jt−jhdt=h∫0nj=0,j=k∏nk−jt−jdt
因为
∏j=0,j≠kn(k−j)=k(k−1)⋯(k−(k−1))⋅(k−(k+1))⋯(k−n)=k(k−1)⋯1⋅(−1)(−2)⋯(−(n−k))=k!(−1)n−k(n−k)!\prod^n_{j=0, j\neq k}(k-j)=k(k-1)\cdots(k-(k-1))·(k-(k+1))\cdots (k-n)\\ =k(k-1)\cdots 1·(-1)(-2)\cdots (-(n-k))=k!(-1)^{n-k}(n-k)! j=0,j=k∏n(k−j)=k(k−1)⋯(k−(k−1))⋅(k−(k+1))⋯(k−n)=k(k−1)⋯1⋅(−1)(−2)⋯(−(n−k))=k!(−1)n−k(n−k)!
所以
Ak=h(−1)n−kk!(n−k)!∫0n∏j=0,j≠kn(t−j)dtA_k=\frac{h(-1)^{n-k}}{k!(n-k)!}\int_0^n\prod_{j=0,j\neq k}^n(t-j)dt Ak=k!(n−k)!h(−1)n−k∫0nj=0,j=k∏n(t−j)dt
令
Ck=Akb−a=Aknh=(−1)n−knk!(n−k)!∫0n∏j=0,j≠kn(t−j)dt(1)C_k=\frac{A_k}{b-a}=\frac{A_k}{nh}=\frac{(-1)^{n-k}}{nk!(n-k)!}\int_0^n\prod_{j=0,j\neq k}^n(t-j)dt \tag{1} Ck=b−aAk=nhAk=nk!(n−k)!(−1)n−k∫0nj=0,j=k∏n(t−j)dt(1)
称CkC_kCk为Cotes（柯特斯）系数。故
∫abf(x)dx≅(b−a)∑k=0nCkf(xk)(2)\int_a^bf(x)dx\cong (b-a)\sum_{k=0}^nC_kf(x_k) \tag{2} ∫abf(x)dx≅(b−a)k=0∑nCkf(xk)(2)
（1）式即为n阶的Newton-Cotes公式的一般形式。

如果给定n，则CkC_kCk可由（1）式计算出来。下标给出了n=1∼10n=1\sim 10n=1∼10的Cotes系数。应该注意

1）表中的系数ckc_kck要乘倍率K，才是公式（1）和（2）中的CkC_kCk。例如，当n=3时，C0=18,C1=38,C2=38,C3=18C_0=\frac{1}{8},C_1=\frac{3}{8},C_2=\frac{3}{8},C_3=\frac{1}{8}C0=81,C1=83,C2=83,C3=81。

2）当n大于等于6时，表中未列出的系数可根据对称性得到。例如，当n=6时，C6=41840C_6=\frac{41}{840}C6=84041。

根据上述规则，可以写出n=4n=4n=4的Newton-Cotes公式如下：
∫abf(x)dx≅190(b−a)[7f(x0)+32f(x1)+12f(x2)+32f(x3)+7f(x4)]\int_a^bf(x)dx \cong \frac{1}{90}(b-a)[7f(x_0)+32f(x_1)+12f(x_2)+32f(x_3)+7f(x_4)] ∫abf(x)dx≅901(b−a)[7f(x0)+32f(x1)+12f(x2)+32f(x3)+7f(x4)]
其中，x0=a,x4=bx_0=a,x_4=bx0=a,x4=b.

2. Newton-Cotes公式的稳定性

在Newton-Cotes公式（2）中，设计算函数值f(xk)(k=1,2,⋯,n)f(x_k)(k=1,2,\cdots,n)f(xk)(k=1,2,⋯,n)所产生的舍入误差为ϵk(k=1,2,⋯,n)\epsilon_k(k=1,2,\cdots,n)ϵk(k=1,2,⋯,n)，则用Newton-Cotes公式计算积分的误差为：
η=(b−a)∑k=0nCkf(xk)−(b−a)∑k=0nCk[f(xk)+ϵk]=(b−a)∑j=0nCkϵk\eta=(b-a)\sum_{k=0}^nC_kf(x_k)-(b-a)\sum_{k=0}^nC_k[f(x_k)+\epsilon_k]=(b-a)\sum_{j=0}^nC_k\epsilon_k η=(b−a)k=0∑nCkf(xk)−(b−a)k=0∑nCk[f(xk)+ϵk]=(b−a)j=0∑nCkϵk
令ϵ=max0≤k≤n∣ϵk∣\epsilon=max_{0\leq k\leq n}|\epsilon_k|ϵ=max0≤k≤n∣ϵk∣，当Cotes系数全为正时(n≤7,n=9)(n\leq 7,n=9)(n≤7,n=9)，则
∣η∣≤(b−a)∑k=0n∣Ckϵk∣≤(b−a)ϵ∑k=0n∣Ck∣=(b−a)ϵ|\eta|\leq (b-a)\sum_{k=0}^n|C_k\epsilon_k|\leq (b-a)\epsilon \sum_{k=0}^n|C_k|=(b-a)\epsilon ∣η∣≤(b−a)k=0∑n∣Ckϵk∣≤(b−a)ϵk=0∑n∣Ck∣=(b−a)ϵ
所以，Newton-Cotes公式是数值稳定的。如果Cotes系数有正有负，则∑k=0n∣Ck∣>1\sum_{k=0}^n|C_k|>1∑k=0n∣Ck∣>1

随着n的增大，该值变得越来越大，并且对∣η∣|\eta|∣η∣的影响越来越大，容易出现数值不稳定现象。因此，高阶(n≥8)(n\geq 8)(n≥8)的Newton-Cotes公式在实际中很少应用。

3. 截断误差与代数精度

截断误差

Newton-Cotes公式是插值型求积公式，故其余项R是Lagrange插值多项式余项的积分：
R=∫abf(x)dx−∑i=0nAif(xi)=∫abf(n+1)(ξ)(n+1)!∏i=0n(x−xi)dx,ξ∈[a,b]R=\int_a^bf(x)dx-\sum_{i=0}^nA_if(x_i)=\int_a^b\frac{f^{(n+1)}(\xi)}{(n+1)!}\prod_{i=0}^n(x-x_i)dx, \quad \xi\in[a,b] R=∫abf(x)dx−i=0∑nAif(xi)=∫ab(n+1)!f(n+1)(ξ)i=0∏n(x−xi)dx,ξ∈[a,b]
即
R=∫abf(n+1)(ξ)(n+1)!∏i=0n(x−xi)dx,ξ∈[a,b]R=\int_a^b\frac{f^{(n+1)}(\xi)}{(n+1)!}\prod_{i=0}^n(x-x_i)dx, \quad \xi \in [a,b] R=∫ab(n+1)!f(n+1)(ξ)i=0∏n(x−xi)dx,ξ∈[a,b]
做积分变换x=a+th,dx=hdt,x=a,t=0;x=b,t=nx=a+th,dx=hdt,x=a,t=0;x=b,t=nx=a+th,dx=hdt,x=a,t=0;x=b,t=n，有：
R=∫0af(n+1)(ξ)(n+1)!(a+th−a−0h)(a+th−a−1h)⋯(a+th−a−nh)hdtR=\int_0^a\frac{f^{(n+1)}(\xi)}{(n+1)!}(a+th-a-0h)(a+th-a-1h)\cdots(a+th-a-nh)hdt R=∫0a(n+1)!f(n+1)(ξ)(a+th−a−0h)(a+th−a−1h)⋯(a+th−a−nh)hdt
故
R=hn+2(n+1)!∫0nf(n+1)(ξ)∏k=0n(t−k)dt,ξ∈[a,b]R=\frac{h^{n+2}}{(n+1)!}\int_0^nf^{(n+1)}(\xi)\prod_{k=0}^n(t-k)dt, \quad \xi\in [a,b] R=(n+1)!hn+2∫0nf(n+1)(ξ)k=0∏n(t−k)dt,ξ∈[a,b]
在推导几个低阶Newton-Cotes公式的截断误差时，会用到广义积分中值定理。即：

定理1：设函数g(x)g(x)g(x)即f(x)g(x)f(x)g(x)f(x)g(x)在[a,b]上可积，函数f(x)f(x)f(x)在[a,b][a,b][a,b]上连续，且m≤f(x)≤Mm\leq f(x)\leq Mm≤f(x)≤M（m，M为常数），函数g(x)g(x)g(x)在[a,b][a,b][a,b]上不变号，即始终有g(x)≤0g(x)\leq 0g(x)≤0或g(x)≥0g(x) \geq 0g(x)≥0，则
∫abf(x)g(x)dx=f(ξ)∫abg(x)dx,ξ∈[a,b]\int_a^bf(x)g(x)dx=f(\xi)\int_a^bg(x)dx, \quad \xi \in[a,b] ∫abf(x)g(x)dx=f(ξ)∫abg(x)dx,ξ∈[a,b]

代数精度

定理2：具有n+1个求积节点的Newton-Cotes公式在n为偶数时，其代数精度为n+1次。