七、数值微分与数值积分

1 数值微分

1.1 差商型求导公式

待补充

两点公式（n=1）

过节点 x0,x1=x0+hx_0,x_1=x_0+hx0,x1=x0+h 的Lagrange插值多项式为：

L1(x)=x−x1−hf(x0)+x−x0hf(x1)L_1(x)=\frac{x-x_1}{-h}f(x_0)+\frac{x-x_0}{h}f(x_1)L1(x)=−hx−x1f(x0)+hx−x0f(x1)

故有：

{f′(x0)≈L1(x0)=f(x1)−f(x0)hf′(x1)≈L1′(x1)=f(x1)−f(x0)h(7−9)\begin{cases}f'(x_0)\approx L_1(x_0)=\frac{f(x_1)-f(x_0)}{h}\\f'(x_1)\approx L_1'(x_1)=\frac{f(x_1)-f(x_0)}{h}\end{cases}\quad\quad\quad\quad\quad\quad\quad(7-9){f′(x0)≈L1(x0)=hf(x1)−f(x0)f′(x1)≈L1′(x1)=hf(x1)−f(x0)(7−9)

截断误差为：

待补充

三点公式

过节点 xi=x0+ih(i=0,1,2)x_i=x_0+ih\ (i=0,1,2)xi=x0+ih (i=0,1,2) , f(x)f(x)f(x) 的Lagrange插值多项式为：

L2(x)=(x−x1)(x−x2)2h2f(x0)+(x−x0)(x−x2)−h2f(x1)f(x2)+(x−x0)(x−x1)2h2f(x2)L_2(x)=\frac{(x-x_1)(x-x_2)}{2h^2}f(x_0)+\frac{(x-x_0)(x-x_2)}{-h^2}f(x_1)f(x_2)+\frac{(x-x_0)(x-x_1)}{2h^2}f(x_2)L2(x)=2h2(x−x1)(x−x2)f(x0)+−h2(x−x0)(x−x2)f(x1)f(x2)+2h2(x−x0)(x−x1)f(x2)

求导得

$$$$

1.2 插值型求导公式

待补充

下面给出等距节点情况下，求一阶导数的几个常用公式。

两点公式（n=1）

过节点 x0,x1=x0+hx_0,x_1=x_0+hx0,x1=x0+h 的Lagrange插值多项式为：

L1(x)=x−x1−hf(x0)+x−x0hf(x1)L_1(x)=\frac{x-x_1}{-h}f(x_0)+\frac{x-x_0}{h}f(x_1)L1(x)=−hx−x1f(x0)+hx−x0f(x1)

故有：

{f′(x0)≈L1′(x0)=f(x1)−f(x0)hf′(x1)≈L1′(x1)=f(x1)−f(x0)h(7−9)\begin{cases}f'(x_0)\approx L_1'(x_0)=\frac{f(x_1)-f(x_0)}{h}\\f'(x_1)\approx L_1'(x_1)=\frac{f(x_1)-f(x_0)}{h}\end{cases}\quad\quad\quad\quad\quad(7-9){f′(x0)≈L1′(x0)=hf(x1)−f(x0)f′(x1)≈L1′(x1)=hf(x1)−f(x0)(7−9)

截断误差（？）为：

{R1′(x0)=−h2f′′(ξ0)R1′(x1)=h2f′′(ξ1)(ξ0,ξ1∈(x0,x1)(7−10)\begin{cases}R_1'(x_0)=-\frac{h}{2}f''(\xi_0)\\R'_1(x_1)=\frac{h}{2}f''(\xi_1)\end{cases}\quad(\xi_0,\xi_1\in(x_0,x_1)\quad\quad\quad\quad\quad\quad(7-10){R1′(x0)=−2hf′′(ξ0)R1′(x1)=2hf′′(ξ1)(ξ0,ξ1∈(x0,x1)(7−10)

三点公式

过节点xi=x0+ih(i=0,1,2),f(x)x_i=x_0+ih\ (i=0,1,2), f(x)xi=x0+ih (i=0,1,2),f(x)的Lagrange插值多项式（？）为：

待补充

分别代入x0,x1,x2x_0,x_1,x_2x0,x1,x2，得三点公式：

{f′(x0)≈12h[−3f(x0)+4f(x1)−f(x2)]f′(x1)≈12h[−f(x0)+f(x2)]f′(x2)≈12h[f(x0)−4f(x1)+3f(x2)](7−11)\begin{cases}f'(x_0)\approx\frac{1}{2h}[-3f(x_0)+4f(x_1)-f(x_2)]\\f'(x_1)\approx\frac{1}{2h}[-f(x_0)+f(x_2)]\\f'(x_2)\approx\frac{1}{2h}[f(x_0)-4f(x_1)+3f(x_2)]\end{cases}\quad\quad\quad\quad\quad\quad(7-11)⎩⎪⎨⎪⎧f′(x0)≈2h1[−3f(x0)+4f(x1)−f(x2)]f′(x1)≈2h1[−f(x0)+f(x2)]f′(x2)≈2h1[f(x0)−4f(x1)+3f(x2)](7−11)

待补充

2 牛顿-科茨（Newton-Cotes）求积公式

2.1 数值积分的基本思想

定积分的定义表明，函数f(x)f(x)f(x)在[a,b][a,b][a,b]上的积分值是和式的极限：

I(f)=∫abf(x)dx=lim⁡λ→0∑i=1nf(ξi)ΔxiI(f)=\int_a^bf(x)\text{d}x=\lim\limits_{\lambda\to0}\sum\limits_{i=1}^nf(\xi_i)\Delta x_iI(f)=∫abf(x)dx=λ→0limi=1∑nf(ξi)Δxi

其中λ=max⁡1≤i≤nΔxi\lambda=\max\limits_{1\le i\le n}\Delta x_iλ=1≤i≤nmaxΔxi，积分和式∑i=1nf(ξi)Δxi\sum\limits_{i=1}^nf(\xi_i)\Delta x_ii=1∑nf(ξi)Δxi是一些点上函数值的线性组合。构造数值积分方法的基本思想就是用被积函数在积分区间[a,b][a,b][a,b]上的某些节点xk(k=0,1,⋯,n)x_k\ (k=0,1,\cdots,n)xk (k=0,1,⋯,n)处的函数值的线性组合作为定积分的近似值，即：

I(f)=∫abf(x)dx≈∑k=0nAkf(xk)(7−16)I(f)=\int_a^bf(x)\text{d}x\approx\sum\limits_{k=0}^nA_kf(x_k)\quad\quad\quad\quad\quad(7-16)I(f)=∫abf(x)dx≈k=0∑nAkf(xk)(7−16)

其中{Ak}\{A_k\}{Ak}称为求积系数，它只与节点{xk}\{x_k\}{xk}有关，与被积函数f(x)f(x)f(x)无关。

由节点{xk}\{x_k\}{xk}和求积系数{Ak}\{A_k\}{Ak}的不同取法，可以得到不同的求积公式。最直接自然的一种想法是用f(x)f(x)f(x)在[a,b][a,b][a,b]上的插值多项式φn(x)\varphi_n(x)φn(x)代替f(x)f(x)f(x)，以φn(x)\varphi_n(x)φn(x)在区间[a,b][a,b][a,b]上的积分值作为所求积分I(f)I(f)I(f)的近似值，即：

I(f)≈∫abφn(x)dxI(f)\approx\int_a^b\varphi_n(x)\text{d}xI(f)≈∫abφn(x)dx

这样得到的求积公式称为插值型求积公式，通常采用Lagrange插值。

设 [a,b][a,b][a,b] 上有 n+1n+1n+1 个互异节点 x0,x1,…,xnx_0,x_1,\dots,x_nx0,x1,…,xn，f(x)f(x)f(x) 的 nnn 次Lagrange插值多项式为：

LN(x)=∑k=0nlk(x)f(xk)L_N(x)=\sum\limits_{k=0}^nl_k(x)f(x_k)LN(x)=k=0∑nlk(x)f(xk)

其中lk(x)=∏j=0,j≠knx−xixk−xjl_k(x)=\prod\limits_{j=0,j\ne k}^n\frac{x-x_i}{x_k-x_j}lk(x)=j=0,j=k∏nxk−xjx−xi，插值型求积公式为：

I(f)≈∫abLn(x)dx=∑k=0n(∫ablk(x)dx)f(xk)(7−17)I(f)\approx\int_a^bL_n(x)\text{d}x=\sum\limits_{k=0}^n(\int_a^bl_k(x)\text{d}x)f(x_k)\quad\quad\quad\quad(7-17)I(f)≈∫abLn(x)dx=k=0∑n(∫ablk(x)dx)f(xk)(7−17)

显然，求积系数Ak=∫ablk(x)dx(k=0,1,⋯,n)A_k=\int_a^bl_k(x)\text{d}x\ (k=0,1,\cdots,n)Ak=∫ablk(x)dx (k=0,1,⋯,n)仅与节点{xk}\{x_k\}{xk}有关，与被积函数f(x)f(x)f(x)的具体形式无关。求积公式（7-17）

待补充

2.2 Newton-Cotes公式

如果取等距节点{xk}\{x_k\}{xk}，则插值型求积公式更简便易算：

设xk=a+kh(k=0,1,⋯,n),h=b−anx_k=a+kh\ (k=0,1,\cdots,n),\ h=\frac{b-a}{n}xk=a+kh (k=0,1,⋯,n), h=nb−a，由式（7-17）：

Ak=∫ablk(x)dx=∫ab(∏j=0,j≠knx−xjxk−xj)dxA_k=\int_a^bl_k(x)\text{d}x=\int_a^b(\prod\limits_{j=0,j\ne k}^n\frac{x-x_j}{x_k-x_j})\text{d}xAk=∫ablk(x)dx=∫ab(j=0,j=k∏nxk−xjx−xj)dx

待补充

令

CK(n)=(−1)n−kk!(n−k)!n∫0n∏j=0,j≠kn(t−j)dt(7−19)C_K^{(n)}=\frac{(-1)^{n-k}}{k!(n-k)!n}\int_0^n\prod\limits_{j=0,j\ne k}^n(t-j)dt\quad\quad\quad\quad\quad(7-19)CK(n)=k!(n−k)!n(−1)n−k∫0nj=0,j=k∏n(t−j)dt(7−19)

则Ak=(b−a)Ck(n)A_k=(b-a)C_k^{(n)}Ak=(b−a)Ck(n)，求积公式（7-17）可简化为：

I(f)=∫abf(x)dx≈(b−a)∑k=0nCK(n)f(xk)(7−20)I(f)=\int_a^bf(x)d_x\approx(b-a)\sum\limits_{k=0}^nC_K^{(n)}f(x_k)\quad\quad\quad\quad\quad(7-20)I(f)=∫abf(x)dx≈(b−a)k=0∑nCK(n)f(xk)(7−20)

式（7-20）称为n阶Newton-Cotes公式，简记为N-C公式，{Ck(n)}\{C_k^{(n)}\}{Ck(n)}称为Cotes系数。

N-C公式的截断误差为：

Rn(f)=∫abf(n+1)(ξ)(n+1)!∏j=0n(x−xj)dx=hn+2(n+1)!∫0nf(n+1)(ξ)∏j=0n(t−j)dt(7−21)R_n(f)=\int_a^b\frac{f^{(n+1)}(\xi)}{(n+1)!}\prod\limits_{j=0}^n(x-x_j)dx\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\\=\frac{h^{n+2}}{(n+1)!}\int_0^nf^{(n+1)}(\xi)\prod\limits_{j=0}^n(t-j)dt\quad\quad\quad\quad\quad(7-21)Rn(f)=∫ab(n+1)!f(n+1)(ξ)j=0∏n(x−xj)dx=(n+1)!hn+2∫0nf(n+1)(ξ)j=0∏n(t−j)dt(7−21)

因为CknC_k^{n}Ckn仅与等分区间数n有关，与积分空间及被积函数均无关，故可事先造好Cotes系数表（？表7-2），这样使用N-C公式计算定积分近似值更加方便。例如，n=4时，N-C公式为：

待补充

式（7-22）是常用的求积公式，称为Cotes公式。

特别地，当n=1时，C0(1)=C1(1)=12C_0^{(1)}=C_1^{(1)}=\frac{1}{2}C0(1)=C1(1)=21，因此有：

I(f)=∫abf(x)dx≈b−a2[f(a)+f(b)](7−23)I(f)=\int_a^bf(x)\text{d}x\approx\frac{b-a}{2}[f(a)+f(b)]\quad\quad\quad\quad\quad(7-23)I(f)=∫abf(x)dx≈2b−a[f(a)+f(b)](7−23)

这是以过点(a,f(a)),(b,f(b))(a,f(a)),(b,f(b))(a,f(a)),(b,f(b))的直线代替曲线y=f(x)y=f(x)y=f(x)，以梯形面积近似曲边梯形面积。故式（7-23）又称为梯形公式。

当n=2时，N-C公式为：

I(f)=∫abf(x)dx=b−a6[f(a)+4f(a+b2)+f(b)](7−24)I(f)=\int_a^bf(x)\text{d}x=\frac{b-a}{6}[f(a)+4f(\frac{a+b}{2})+f(b)]\quad\quad\quad\quad(7-24)I(f)=∫abf(x)dx=6b−a[f(a)+4f(2a+b)+f(b)](7−24)

这是以过曲线上三点(a,f(a)),(a+b2,f(a+b2)),(b,f(b))(a,f(a)),(\frac{a+b}{2},f(\frac{a+b}{2})),(b,f(b))(a,f(a)),(2a+b,f(2a+b)),(b,f(b))的抛物线代替曲线y=f(x)y=f(x)y=f(x)求曲边梯形面积的近似值，因而此公式称为抛物线求积公式，又称辛普森（Simpson）公式。

2.3 误差估计

首先引入一种衡量数值积分公式近似程度的概念。

求积公式的代数精确度

定义7.1

若当f(x)f(x)f(x)为任意次数不高于m的多项式时，求积公式

∫abf(x)dx≈∑k=0nAkf(xk)\int_a^bf(x)dx\approx\sum\limits_{k=0}^nA_kf(x_k)∫abf(x)dx≈k=0∑nAkf(xk)

均精确成立，而对某个m+1次多项式，公式不精确成立，则称该求积公式具有m次代数精确度。

可以验证，梯形公式具有一次代数精确度。因为若f(x)f(x)f(x)是一次多项式，则：

R1(x)=f(x)−L1(x)=f′′(ξ)2(x−a)(x−b)=0R_1(x)=f(x)-L_1(x)=\frac{f''(\xi)}{2}(x-a)(x-b)=0R1(x)=f(x)−L1(x)=2f′′(ξ)(x−a)(x−b)=0

从而有：

∫abf(x)dx=∫abL1(x)dx=b−a2[f(a)+f(b)]\int_a^bf(x)\text{d}x=\int_a^bL_1(x)dx=\frac{b-a}{2}[f(a)+f(b)]∫abf(x)dx=∫abL1(x)dx=2b−a[f(a)+f(b)]

若取f(x)=x2f(x)=x^2f(x)=x2，则：

∫abf(x)dx=∫abx2dx=13(b3−a3)≠b−a2(a2+b2)\int_a^bf(x)dx=\int_a^bx^2dx=\frac{1}{3}(b^3-a^3)\ne\frac{b-a}{2}(a^2+b^2)∫abf(x)dx=∫abx2dx=31(b3−a3)=2b−a(a2+b2)

一般地，由n次插值多项式导出的求积公式至少有n次代数精确度。对N-C公式，更有以下结论：

定理7.1

2n阶Newton-Cotes公式至少具有2n+1次代数精确度。

容易证明，求积公式（7-16）具有m次代数精确度的充要条件是它对f(x)=1,x,⋯,xmf(x)=1,x,\cdots,x^mf(x)=1,x,⋯,xm都能精确成立，但对f(x)=xm+1f(x)=x^{m+1}f(x)=xm+1，式（7-16）不准确成立。这一结论给出了判别一个求积公式的代数精确度的方法。

待补充 215

利用余项公式及积分中值定理，可以导出低阶N-C公式的误差估计。

梯形公式和辛普森公式的误差估计

由式（7-21），梯形公式的截断误差为：

R1(f)=∫abf′′(ξx)2(x−a)(x−b)dx,ξx∈(a,b)R_1(f)=\int_a^b\frac{f''(\xi_x)}{2}(x-a)(x-b)\text{d}x,\quad\xi_x\in(a,b)R1(f)=∫ab2f′′(ξx)(x−a)(x−b)dx,ξx∈(a,b)

因为(x-a)(x-b)在区间[a,b][a,b][a,b]内不变号，如果f(x)f(x)f(x)二阶连续可导，则由反常积分中值定理（？），存在ξ∈(a,b)\xi\in(a,b)ξ∈(a,b)，使得：

R1(f)=f′′(ξ)2∫ab(x−a)(x−b)dx=f′′(ξ)2[x33−(a+b)x22+abx]ab=−f′′(ξ)12(b−a)3R_1(f)=\frac{f''(\xi)}{2}\int_a^b(x-a)(x-b)\text{d}x\\=\frac{f''(\xi)}{2}[\frac{x^3}{3}-(a+b)\frac{x^2}{2}+abx]_a^b\\=-\frac{f''(\xi)}{12}(b-a)^3R1(f)=2f′′(ξ)∫ab(x−a)(x−b)dx=2f′′(ξ)[3x3−(a+b)2x2+abx]ab=−12f′′(ξ)(b−a)3

定理7.2

（？）这个C是啥意思

若f(x)∈C(2)[a,b]f(x)\in C^{(2)}[a,b]f(x)∈C(2)[a,b]，则梯形公式的误差为：

R1(f)=∫abf(x)dx−b−a2[f(a)+f(b)]=−(b−a)312f′′(ξ),ξ∈(a,b)(7−25)R_1(f)=\int_a^bf(x)\text{d}x-\frac{b-a}{2}[f(a)+f(b)]=-\frac{(b-a)^3}{12}f''(\xi),\quad\xi\in(a,b)\quad\quad\quad\quad\quad(7-25)R1(f)=∫abf(x)dx−2b−a[f(a)+f(b)]=−12(b−a)3f′′(ξ),ξ∈(a,b)(7−25)

借助于f(x)f(x)f(x)的Hermite插值与定理7.1（？），可导出Simpson公式的误差估计：

定理7.3

若f(x)∈C(4)[a,b]f(x)\in C^{(4)}[a,b]f(x)∈C(4)[a,b]，则Simpson公式的误差为：

R2(f)=∫abf(x)dx−b−a6[f(a)+4f(a+b2)+f(b)]=−(b−a)52880f(4)(ξ)=−h590f(4)(ξ),ξ∈(a,b)(7−26)R_2(f)=\int_a^bf(x)\text{d}x-\frac{b-a}{6}[f(a)+4f(\frac{a+b}{2})+f(b)]=-\frac{(b-a)^5}{2880}f^{(4)}(\xi)=-\frac{h^5}{90}f^{(4)}(\xi),\quad\xi\in(a,b)\quad\quad\quad(7-26)R2(f)=∫abf(x)dx−6b−a[f(a)+4f(2a+b)+f(b)]=−2880(b−a)5f(4)(ξ)=−90h5f(4)(ξ),ξ∈(a,b)(7−26)

其中h=(b−a)2h=\frac{(b-a)}{2}h=2(b−a)为步长。

C(2)[a,b]C^{(2)}[a,b]C(2)[a,b]应该是对[a,b][a,b][a,b]二等分的含义。

3 复化求积公式

前面导出的误差估计式表明，用N-C公式计算积分近似值时，步长越小，截断误差越小（？）

待补充

3.1 复化梯形公式

用分段线性插值函数近似被积函数，等于把积分区间分成若干小区间，在每个小区间上以梯形面积近似曲边梯形面积，即用梯形公式求小区间上积分的近似值。这样求得的近似值显然比用梯形公式计算精度高。定积分存在定理表明，只要被积函数连续，当小区间长度趋于零时，小梯形面积之和趋于曲边梯形面积的准确值，即定积分的准确值。

将积分区间 [a,b][a,b][a,b] 分成 nnn 等份，记h=b−an,xk=a+kh(k=0,1,⋯,n)h=\frac{b-a}{n},\ x_k=a+kh\ (k=0,1,\cdots,n)h=nb−a, xk=a+kh (k=0,1,⋯,n)，在每个小区间[xk,xk+1](k=0,1,⋯,n−1)[x_k,x_{k+1}]\ (k=0,1,\cdots,n-1)[xk,xk+1] (k=0,1,⋯,n−1)上用梯形公式并求和，得：

∫abf(x)dx=∑k=0n−1∫kk+1f(x)dx≈∑k=0n−1h2[f(xk)+f(xk+1)]\int_a^bf(x)\text{d}x=\sum\limits_{k=0}^{n-1}\int_k^{k+1}f(x)\text{d}x\approx\sum\limits_{k=0}^{n-1}\frac{h}{2}[f(x_k)+f(x_{k+1})]∫abf(x)dx=k=0∑n−1∫kk+1f(x)dx≈k=0∑n−12h[f(xk)+f(xk+1)]

整理得：

∫abf(x)dx≈h2[f(a)+f(b)+2∑k=1n−1f(xk)]=Tn(7−28)\int_a^bf(x)\text{d}x\approx\frac{h}{2}[f(a)+f(b)+2\sum\limits_{k=1}^{n-1}f(x_k)]=T_n\quad\quad\quad\quad\quad(7-28)∫abf(x)dx≈2h[f(a)+f(b)+2k=1∑n−1f(xk)]=Tn(7−28)

式（7-28）称为复化梯形公式。

如果f(x)∈C(2)[a,b]f(x)\in C^{(2)}[a,b]f(x)∈C(2)[a,b]（？），在小区间[xk,xk+1][x_k,x_{k+1}][xk,xk+1]上，梯形公式的截断误差为：

∫xkxk+1dx−h2[f(xk)+f(xk+1)]=−h312f′′(ξk),ξk∈(xk,xk+1)\int_{x_k}^{x_{k+1}}dx-\frac{h}{2}[f(x_k)+f(x_{k+1})]=-\frac{h^3}{12}f''(\xi_k),\quad\xi_k\in(x_k,x_{k+1})∫xkxk+1dx−2h[f(xk)+f(xk+1)]=−12h3f′′(ξk),ξk∈(xk,xk+1)
？？？？ 12 是哪来的

因此

RT(f)=∫abf(x)dx−Tn=−h312∑k=0n−1f′′(ξk)R_T(f)=\int_a^bf(x)dx-T_n=-\frac{h^3}{12}\sum\limits_{k=0}^{n-1}f''(\xi_k)RT(f)=∫abf(x)dx−Tn=−12h3k=0∑n−1f′′(ξk)

因为f′′(x)f''(x)f′′(x)在[a,b]连续，由介值定理，存在ξ∈(a,b)\xi\in(a,b)ξ∈(a,b)，使得：

f′′(ξ)=1n∑k=0n−1f′′(ξk)f''(\xi)=\frac{1}{n}\sum\limits_{k=0}^{n-1}f''(\xi_k)f′′(ξ)=n1k=0∑n−1f′′(ξk)

从而有：

RT(f)=∫abf(x)dx−Tn=−h312nf′′(ξ)=−b−a12h2f′′(ξ),ξ∈(a,b)(7−29)R_T(f)=\int_a^bf(x)dx-T_n=-\frac{h^3}{12}nf''(\xi)=-\frac{b-a}{12}h^2f''(\xi),\quad\xi\in(a,b)\quad\quad\quad\quad\quad\quad\quad(7-29)RT(f)=∫abf(x)dx−Tn=−12h3nf′′(ξ)=−12b−ah2f′′(ξ),ξ∈(a,b)(7−29)

这就是复化梯形公式的截断误差。

3.2 复化Simpson公式

如果用分段二次插值函数近似被积函数，即在小区间上用Simpson公式计算积分近似值，就导出复化Simpson公式。

将区间[a,b]分成n=2m等份，分点为xk=a+kh(k=0,1,⋯,n),h=b−anx_k=a+kh\ (k=0,1,\cdots,n),h=\frac{b-a}{n}xk=a+kh (k=0,1,⋯,n),h=nb−a，在每个小区间[x2k−2,x2k][x_{2k-2},x_{2k}][x2k−2,x2k]上，用Simpson公式求积分，则有：

待补充

整理后得到：

∫abf(x)dx≈h3[f(a)+f(b)+2∑k=1m−1f(x2k)+4∑k=1mf(x2k−1)](7−30)\int_a^bf(x)\text{d}x\approx\frac{h}{3}[f(a)+f(b)+2\sum\limits_{k=1}^{m-1}f(x_{2k})+4\sum\limits_{k=1}^mf(x_{2k-1})]\quad\quad\quad\quad\quad(7-30)∫abf(x)dx≈3h[f(a)+f(b)+2k=1∑m−1f(x2k)+4k=1∑mf(x2k−1)](7−30)

式（7-30）称为复化Simpson公式。

4 龙贝格（Romberg）求积公式

4.1 理查逊（Richardson）外推法

设用某种数值方法求量III的近似值，一般地，近似值是步长hhh的函数，记为I1(h)I_1(h)I1(h)，相应的误差为：

I−I1(h)=a1hp1+a2hp2+⋯+akhpk+⋯(7−36)I-I_1(h)=a_1h^{p_1}+a_2h^{p_2}+\cdots+a_kh^{p_k}+\cdots\quad\quad\quad\quad\quad\quad(7-36)I−I1(h)=a1hp1+a2hp2+⋯+akhpk+⋯(7−36)

其中0<p1<p2<⋯<pk<⋯,ai(i=1,2,⋯)0<p_1<p_2<\cdots<p_k<\cdots,\ a_i(i=1,2,\cdots)0<p1<p2<⋯<pk<⋯, ai(i=1,2,⋯)与h无关。若用αh\alpha hαh代替式（7-36）中的h，则得：

I−I1(αh)=a1αp1hp1+a2αp2hp2+⋯+akαpkhpk+⋯(7−37)I-I_1(\alpha h)=a_1\alpha^{p_1}h^{p_1}+a_2\alpha^{p_2}h^{p_2}+\cdots+a_k\alpha^{p_k}h^{p_k}+\cdots\quad\quad\quad\quad\quad(7-37)I−I1(αh)=a1αp1hp1+a2αp2hp2+⋯+akαpkhpk+⋯(7−37)

式（7-37）减去式（7-36）乘以αp1\alpha^{p_1}αp1，得：

I−I1(αh)−αp1[I−I1(h)]=a2(αp2−αp1)hp2+a3(αp3−αp1)hp3+⋯+ak(αpk−αp1)hpk+⋯I-I_1(\alpha h)-\alpha^{p_1}[I-I_1(h)]\\=a_2(\alpha^{p_2}-\alpha^{p_1})h^{p_2}+a_3(\alpha^{p_3}-\alpha^{p_1})h^{p_3}+\cdots+a_k(\alpha^{p_k}-\alpha^{p_1})h^{p_k}+\cdotsI−I1(αh)−αp1[I−I1(h)]=a2(αp2−αp1)hp2+a3(αp3−αp1)hp3+⋯+ak(αpk−αp1)hpk+⋯

取α\alphaα满足∣α∣≠1|\alpha|\ne1∣α∣=1，以1−αp11-\alpha^{p_1}1−αp1除上式两端，得：

I−I1(αh)−αp1I1(h)1−αp1=b2hp2+b3hp3+⋯+bkhp+k+⋯(7−38)I-\frac{I_1(\alpha h)-\alpha^{p_1}I_1(h)}{1-\alpha^{p_1}}=b_2h^{p_2}+b_3h^{p_3}+\cdots+b_kh^{p+k}+\cdots\quad\quad\quad\quad\quad(7-38)I−1−αp1I1(αh)−αp1I1(h)=b2hp2+b3hp3+⋯+bkhp+k+⋯(7−38)

其中bi=ai(αpi−αp1)/(1−αp1)(i=2,3,…)b_i=a_i(\alpha^{p_i}-\alpha^{p_1})/(1-\alpha^{p_1})\ (i=2,3,\dots)bi=ai(αpi−αp1)/(1−αp1) (i=2,3,…)仍与h无关。令：

I2(h)=I1(αh)−αp1I1(h)1−αp1I_2(h)=\frac{I_1(\alpha h)-\alpha^{p_1}I_1(h)}{1-\alpha^{p_1}}I2(h)=1−αp1I1(αh)−αp1I1(h)

由式（7-38），以I2(h)I_2(h)I2(h)作为I的近似值，其误差至少为O(hp2)O(h^{p_2})O(hp2)（？），因此I2(h)I_2(h)I2(h)收敛于III的速度比I1(h)I_1(h)I1(h)快。不断重复以上做法，可以得到一个函数序列：

Im(h)=Im−1(αh)=αPm−1Im−1(h)1−αPm−1(m=2,3,⋯)I_m(h)=\frac{I_{m-1}(\alpha h)=\alpha^{P_{m-1}}I_{m-1}(h)}{1-\alpha^{P_{}m-1}}\quad(m=2,3,\cdots)Im(h)=1−αPm−1Im−1(αh)=αPm−1Im−1(h)(m=2,3,⋯)

以Im(h)I_m(h)Im(h)近似I，误差为I−Im(h)=O(hPm)I-I_m(h)=O(h^{P_m})I−Im(h)=O(hPm)。随着m的增大，收敛速度越来越快，这就是Richardson外推法。将此方法用于梯形值序列{T2k}\{T_{2^k}\}{T2k}，可以构造出一种收敛很快，计算简便的数值积分方法。

4.2 Romberg求积公式

设 I=∫abf(x)dxI=\int_a^bf(x)\text{d}xI=∫abf(x)dx ，记 2k2^k2k 等分区间 [a,b][a,b][a,b] 用复化梯形公式所求得的近似值为 T2k=T0(k)=I1(b−a2k)(k=0,1,2,…)T_{2^k}=T_0^{(k)}=I_1(\frac{b-a}{2^k})\ (k=0,1,2,\dots)T2k=T0(k)=I1(2kb−a) (k=0,1,2,…) 。由Euler-Maclaurin公式

R(f,T0(k))=a1(b−a2k)2+a2(b−a2k)4+⋯+ai(b−a2k)2i+⋯R(f,T_0^{(k)})=a_1(\frac{b-a}{2^{k}})^2+a_2(\frac{b-a}{2^{k}})^4+\cdots+a_i(\frac{b-a}{2^k})^{2i}+\cdotsR(f,T0(k))=a1(2kb−a)2+a2(2kb−a)4+⋯+ai(2kb−a)2i+⋯

取α=12\alpha=\frac{1}{2}α=21，由Richardson外推公式：

I2(b−a2k)=I1(b−a2k+1)−(12)2I1(b−a2k)1−(12)2=4T0(k+1)−T0(k)3I_2(\frac{b-a}{2^k})=\frac{I_1(\frac{b-a}{2^{k+1}})-(\frac{1}{2})^2I_1(\frac{b-a}{2^k})}{1-(\frac{1}{2})^2}=\frac{4T_0^{(k+1)}-T_0^{(k)}}{3}I2(2kb−a)=1−(21)2I1(2k+1b−a)−(21)2I1(2kb−a)=34T0(k+1)−T0(k)

记I2=(b−a2k)=T1(k)I_2=(\frac{b-a}{2^k})=T_1^{(k)}I2=(2kb−a)=T1(k)，则：

T1(k)=4T0(k+1)−T0(k)3(k=0,1,2⋯)(7−39)T_1^{(k)}=\frac{4T_0^{(k+1)}-T_0^{(k)}}{3}\quad(k=0,1,2\cdots)\quad\quad\quad\quad(7-39)T1(k)=34T0(k+1)−T0(k)(k=0,1,2⋯)(7−39)

且有：

R(f,T1(k))=O((b−a2k)4)R(f,T_1^{(k)})=O((\frac{b-a}{2^k})^4)R(f,T1(k))=O((2kb−a)4)

事实上，容易验证，{T1(k)}\{T_1^{(k)}\}{T1(k)}就是Simpson值序列{S2k}\{S_{2^k}\}{S2k}

待补充

记Tm(k)=Im+1(b−a2k)T_m^{(k)}=I_{m+1}(\frac{b-a}{2^k})Tm(k)=Im+1(2kb−a)，则上式可记为：

Tm(k)=4mTm−1(k+1)−Tm−1(k)4m−1(m=1,2,⋯;k=0,1,2,⋯)(7−40)T_m^{(k)}=\frac{4^mT_{m-1}^{(k+1)}-T_{m-1}^{(k)}}{4^m-1}\quad(m=1,2,\cdots;k=0,1,2,\cdots)\quad\quad\quad\quad\quad(7-40)Tm(k)=4m−14mTm−1(k+1)−Tm−1(k)(m=1,2,⋯;k=0,1,2,⋯)(7−40)

5 高斯（Gauss）型求积公式

5.1 一般理论

Gauss型求积公式

Gauss型求积公式的基本思想是：在节点数n固定时，适当地选取节点{xk}\{x_k\}{xk}与求积系数{Ak}\{A_k\}{Ak}，使求积公式具有最高的代数精确度。

由定义（7.1），求积公式

∫abρ(x)f(x)dx≈∑k=1nAkf(xk)(7−41)\int_a^b\rho(x)f(x)\text{d}x\approx\sum\limits_{k=1}^nA_kf(x_k)\quad\quad\quad\quad\quad(7-41)∫abρ(x)f(x)dx≈k=1∑nAkf(xk)(7−41)

其中ρ(x)≥0\rho(x)\ge0ρ(x)≥0为权函数，具有m次代数精确度，充要条件是对f(x)=xl(l=0,1,⋯,m)f(x)=x^l\ (l=0,1,\cdots,m)f(x)=xl (l=0,1,⋯,m)式（7-41）精确成立，即：

∫abρ(x)xldx≈∑k=1nAkxkl(7−42)\int_a^b\rho(x)x^l\text{d}x\approx\sum\limits_{k=1}^nA_kx_k^l\quad\quad\quad\quad\quad(7-42)∫abρ(x)xldx≈k=1∑nAkxkl(7−42)

令：

μl=∫abρ(x)xldx(l=0,1,⋯,m)(7−43)\mu_l=\int_a^b\rho(x)x^l\text{d}x\quad(l=0,1,\cdots,m)\quad\quad\quad\quad\quad(7-43)μl=∫abρ(x)xldx(l=0,1,⋯,m)(7−43)

将式（7-43）代入式（7-42），得节点{xk}\{x_k\}{xk}与求积系数{Ak}\{A_k\}{Ak}应满足方程组：

{A1+A2+⋯+An=μ0A1xa+A2x2+⋯+Anxn=μ1A1x1m+A2x2m+⋯+Anxnm=μm(7−44)\begin{cases}A_1+A_2+\cdots+A_n=\mu_0\\A_1x_a+A_2x_2+\cdots+A_nx_n=\mu_1\\A_1x_1^m+A_2x_2^m+\cdots+A_nx_n^m=\mu_m\end{cases}\quad\quad\quad\quad\quad(7-44)⎩⎪⎨⎪⎧A1+A2+⋯+An=μ0A1xa+A2x2+⋯+Anxn=μ1A1x1m+A2x2m+⋯+Anxnm=μm(7−44)

式（7-44）含有2n个未知数，m+1个方程，可以证明，当m+1≤2nm+1\le 2nm+1≤2n时，方程组（7-44）有解。这表明n个节点的求积公式的代数精确度可达到2n-1。但若取2n次多项式：

P2n(x)=(x−x1)2(x−x2)2⋯(x−xn)2P_{2n}(x)=(x-x_1)^2(x-x_2)^2\cdots(x-x_n)^2P2n(x)=(x−x1)2(x−x2)2⋯(x−xn)2

显然，对任意n个互异节点x1,x2,⋯,xnx_1,x_2,\cdots,x_nx1,x2,⋯,xn，都有P2n(xk)=0(k=1,2,⋯,n)P_{2n}(x_k)=0\ (k=1,2,\cdots,n)P2n(xk)=0 (k=1,2,⋯,n)，从而有∑k=1nP2n(xk)=0\sum\limits_{k=1}^nP_{2n}(x_k)=0k=1∑nP2n(xk)=0。然而：

∫abρ(x)P2n(x)dx>0\int_a^b\rho(x)P_{2n}(x)\text{d}x>0∫abρ(x)P2n(x)dx>0

因此，n个节点的求积公式的最高代数精确度为2n-1（？）

定义 7.2

如果一组节点x1,x2,⋯,xn∈[a,b]x_1,x_2,\cdots,x_n\in[a,b]x1,x2,⋯,xn∈[a,b]能使求积公式：

∫abρ(x)f(x)dx≈∑k=1nAkf(xk)(7−41)\int_a^b\rho(x)f(x)\text{d}x\approx\sum\limits_{k=1}^nA_kf(x_k)\quad\quad\quad\quad\quad(7-41)∫abρ(x)f(x)dx≈k=1∑nAkf(xk)(7−41)

具有2n-1次代数精确度，则称这组节点{xk}\{x_k\}{xk}为Gauss点，公式（7-41）称为带权函数ρ(x)\rho(x)ρ(x)的Gauss型求积公式。

Gauss点与正交多项式的关系