\qquad 使用牛顿-莱布尼兹公式直接去求一个函数 f ( x ) f(x) f(x) 在区间 [ a , b ] [a,b] [a,b] 上的定积分值 I I I，在数值方法上比较困难，因为需要先求出被积函数的原函数 F ( x ) F(x) F(x)。

I = ∫ a b f ( x ) d x = F ( b ) − F ( a ) \qquad\qquad I=\displaystyle\int_a^bf(x)dx=F(b)-F(a) I=∫ab​f(x)dx=F(b)−F(a)

\qquad 通过积分中值定理，在积分区间 [ a , b ] [a,b] [a,b] 上存在某个点 ξ \xi ξ，使得

I = ∫ a b f ( x ) d x = ( b − a ) f ( ξ ) , ξ ∈ ( a , b ) \qquad\qquad I=\displaystyle\int_a^bf(x)dx=(b-a)f(\xi),\quad\xi\in(a,b) I=∫ab​f(x)dx=(b−a)f(ξ),ξ∈(a,b)
\qquad

1. 机械求积

\qquad 根据积分中值定理，只需要合理选择 ξ ∈ ( a , b ) \xi\in(a,b) ξ∈(a,b) 的值，就可以求出定积分，例如最基本的梯形公式和中矩形公式。
\qquad

用红色虚线矩形的面积，代替 f ( x ) f(x) f(x) 与 x x x 轴所围的面积。

• 梯形公式：　　取 f ( ξ ) ≈ f ( a ) + f ( b ) 2 f(\xi)\approx\dfrac{f(a)+f(b)}{2} f(ξ)≈2f(a)+f(b)​

∫ a b f ( x ) d x = b − a 2 [ f ( a ) + f ( b ) ] \qquad\qquad\displaystyle\int_a^bf(x)dx=\dfrac{b-a}{2}[f(a)+f(b)] ∫ab​f(x)dx=2b−a​[f(a)+f(b)]

• 中矩形公式：　取 f ( ξ ) ≈ f ( a + b 2 ) f(\xi)\approx f\left(\dfrac{a+b}{2}\right) f(ξ)≈f(2a+b​)

∫ a b f ( x ) d x = ( b − a ) f ( a + b 2 ) \qquad\qquad\displaystyle\int_a^bf(x)dx=(b-a)f\left(\dfrac{a+b}{2}\right) ∫ab​f(x)dx=(b−a)f(2a+b​)

\qquad 更一般地，可以在积分区间 [ a , b ] [a,b] [a,b] 上选取一些求积节点 x k x_k xk​，用 f ( x k ) f(x_k) f(xk​) 的加权平均（权值为 A k A_k Ak​）来近似 f ( ξ ) f(\xi) f(ξ) 的值，也就是采用“机械求积”方法：

f ( ξ ) ≈ ∑ k = 0 n A k f ( x k ) \qquad\qquad f(\xi)\approx\displaystyle\sum_{k=0}^nA_kf(x_k) f(ξ)≈k=0∑n​Ak​f(xk​) ⟹ ∫ a b f ( x ) d x = ( b − a ) ∑ k = 0 n A k f ( x k ) \ \ \ \Longrightarrow\displaystyle\int_a^bf(x)dx=(b-a)\displaystyle\sum_{k=0}^nA_kf(x_k)    ⟹∫ab​f(x)dx=(b−a)k=0∑n​Ak​f(xk​)

梯形公式：　 n = 2 , A 0 = 1 , x 0 = a , A 1 = 1 , x 1 = b n=2,A_0=1,x_0=a,A_1=1,x_1=b n=2,A0​=1,x0​=a,A1​=1,x1​=b
中矩形公式： n = 1 , A 0 = 1 , x 0 = a + b 2 n=1,A_0=1,x_0=\dfrac{a+b}{2} n=1,A0​=1,x0​=2a+b​

2. Newton-Cotes公式

\qquad 将积分区间 [ a , b ] [a,b] [a,b] 划分为 n n n 等分，步长为 h = b − a n h=\dfrac{b-a}{n} h=nb−a​，选取等距求积节点 x k = a + k h x_k=a+kh xk​=a+kh 构造出牛顿-科特斯 (Newton-Cotes) \text{(Newton-Cotes)} (Newton-Cotes) 公式：

I n = ( b − a ) ∑ k = 0 n C k ( n ) f ( x k ) \qquad\qquad I_n=(b-a)\displaystyle\sum_{k=0}^nC_k^{(n)}f(x_k) In​=(b−a)k=0∑n​Ck(n)​f(xk​)

\qquad 其中， C k ( n ) C_k^{(n)} Ck(n)​ 为科特斯系数，若 x = a + t h x=a+th x=a+th，则

C k ( n ) = h b − a ∫ 0 n ∏ j = 0 j ≠ k n t − j k − j d t = ( − 1 ) n − k n k ! ( n − k ) ! ∫ 0 n ∏ j = 0 j ≠ k n ( t − j ) d t \qquad\qquad C_k^{(n)}=\dfrac{h}{b-a}\displaystyle\int_0^n\displaystyle\prod^n_{j=0\atop j\ne k}\dfrac{t-j}{k-j}dt=\dfrac{(-1)^{n-k}}{nk!(n-k)!}\displaystyle\int_0^n\displaystyle\prod_{j=0\atop j\ne k}^n(t-j)dt Ck(n)​=b−ah​∫0n​j=kj=0​∏n​k−jt−j​dt=nk!(n−k)!(−1)n−k​∫0n​j=kj=0​∏n​(t−j)dt

\qquad

• 当 n = 1 n=1 n=1 时， C 0 ( 1 ) = C 1 ( 1 ) = 1 2 C_0^{(1)}=C_1^{(1)}=\dfrac{1}{2} C0(1)​=C1(1)​=21​，也就是梯形公式（单区间， x 0 = a , x 1 = b x_0=a,x_1=b x0​=a,x1​=b）

  T = I 1 = b − a 2 [ f ( a ) + f ( b ) ] \qquad\qquad T=I_1=\dfrac{b-a}{2}[f(a)+f(b)] T=I1​=2b−a​[f(a)+f(b)]

• 当 n = 2 n=2 n=2 时， C 0 ( 2 ) = 1 6 , C 1 ( 2 ) = 4 6 , C 2 ( 2 ) = 1 6 C_0^{(2)}=\dfrac{1}{6},C_1^{(2)}=\dfrac{4}{6},C_2^{(2)}=\dfrac{1}{6} C0(2)​=61​,C1(2)​=64​,C2(2)​=61​，也就是辛普森 (Simpson) \text{(Simpson)} (Simpson)公式（ 2 2 2 等分区间）

  S = I 2 = b − a 6 [ f ( a ) + 4 f ( a + b 2 ) + f ( b ) ] \qquad\qquad S=I_2=\dfrac{b-a}{6}\left[f(a)+4f\left(\dfrac{a+b}{2}\right)+f(b)\right] S=I2​=6b−a​[f(a)+4f(2a+b​)+f(b)]

• 当 n = 4 n=4 n=4 时，也就是科特斯公式（ 4 4 4 等分区间）

  C = I 4 = b − a 90 [ 7 f ( a ) + 32 f ( x 1 ) + 12 f ( x 2 ) + 32 f ( x 3 ) + 7 f ( b ) ] \qquad\qquad C=I_4=\dfrac{b-a}{90}\left[7f(a)+32f(x_1)+12f(x_2)+32f(x_3)+7f(b)\right] C=I4​=90b−a​[7f(a)+32f(x1​)+12f(x2​)+32f(x3​)+7f(b)]
  \qquad 其中， x 0 = a , x 4 = b , x k = a + k h , h = b − a 4 x_0=a,x_4=b,x_k=a+kh,h=\dfrac{b-a}{4} x0​=a,x4​=b,xk​=a+kh,h=4b−a​
  　
  【注意】当 n ≥ 8 n\ge8 n≥8 时，科特斯系数 C k ( n ) C_k^{(n)} Ck(n)​ 出现负值，求积计算变得不稳定。

\qquad

3. 复化求积方法

\qquad 由于高阶的牛顿-科特斯公式不稳定，无法通过提高阶数 n n n 来提高求积精度。

\qquad 为了提高精度，可以采用复化求积方法：将积分区间 [ a , b ] [a,b] [a,b] 分成若干个子区间 [ x k , x k + 1 ] [x_k,x_{k+1}] [xk​,xk+1​]，在每个子区间上采用低阶 ( n < 8 ) (n<8) (n<8) 的牛顿-科特斯公式。
\qquad

3.1 复化梯形公式

\qquad 将积分区间 [ a , b ] [a,b] [a,b] 进行 n n n 等分，步长为 h = b − a n h=\dfrac{b-a}{n} h=nb−a​，等距求积节点 x k = a + k h , k = 0 , 1 , ⋯ , n x_k=a+kh,k=0,1,\cdots,n xk​=a+kh,k=0,1,⋯,n，在每个子区间 [ x k , x k + 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 上采用梯形公式：

T ( k ) = x k + 1 − x k 2 [ f ( x k ) + f ( x k + 1 ) ] = h 2 [ f ( x k ) + f ( x k + 1 ) ] \qquad\qquad T^{(k)}=\dfrac{x_{k+1}-x_k}{2}[f(x_k)+f(x_{k+1})]=\dfrac{h}{2}[f(x_k)+f(x_{k+1})] T(k)=2xk+1​−xk​​[f(xk​)+f(xk+1​)]=2h​[f(xk​)+f(xk+1​)]

\qquad 记 n n n 等分时求得积分的结果为 T n T_n Tn​，那么

∫ a b f ( x ) d x ≈ T n = ∑ k = 0 n − 1 T ( k ) = h 2 ∑ k = 0 n − 1 [ f ( x k ) + f ( x k + 1 ) ] \qquad\qquad\displaystyle\int_a^bf(x)dx\approx T_n=\displaystyle\sum_{k=0}^{n-1}T^{(k)}=\dfrac{h}{2}\displaystyle\sum_{k=0}^{n-1}[f(x_k)+f(x_{k+1})] ∫ab​f(x)dx≈Tn​=k=0∑n−1​T(k)=2h​k=0∑n−1​[f(xk​)+f(xk+1​)]

误差为 R n ( f ) = − b − a 12 h 2 f ′ ′ ( η ) , η ∈ ( a , b ) R_n(f)=-\dfrac{b-a}{12}h^2f^{''}(\eta),\quad\eta\in(a,b) Rn​(f)=−12b−a​h2f′′(η),η∈(a,b) 为 h 2 h^2 h2 阶

\qquad

3.2 复化辛普森公式

\qquad 将积分区间 [ a , b ] [a,b] [a,b] 进行 n n n 等分，步长为 h = b − a n h=\dfrac{b-a}{n} h=nb−a​，等距求积节点 x k = a + k h , k = 0 , 1 , ⋯ , n x_k=a+kh,k=0,1,\cdots,n xk​=a+kh,k=0,1,⋯,n，在每个子区间 [ x k , x k + 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 上采用辛普森公式：

S ( k ) = x k + 1 − x k 6 [ f ( x k ) + 4 f ( x k + 1 2 ) + f ( x k + 1 ) ] = h 6 [ f ( x k ) + 4 f ( x k + 1 2 ) + f ( x k + 1 ) ] \qquad\qquad\begin{aligned}S^{(k)}&=\dfrac{x_{k+1}-x_k}{6}\left[f(x_k)+4f\left(x_{k+\frac{1}{2}}\right)+f(x_{k+1})\right]\\ &=\dfrac{h}{6}\left[f(x_k)+4f\left(x_{k+\frac{1}{2}}\right)+f(x_{k+1})\right]\end{aligned} S(k)​=6xk+1​−xk​​[f(xk​)+4f(xk+21​​)+f(xk+1​)]=6h​[f(xk​)+4f(xk+21​​)+f(xk+1​)]​

\qquad 记 n n n 等分时求得积分的结果为 S n S_n Sn​，那么

∫ a b f ( x ) d x ≈ S n = ∑ k = 0 n − 1 S ( k ) = h 6 ∑ k = 0 n − 1 [ f ( x k ) + 4 f ( x k + 1 2 ) + f ( x k + 1 ) ] \qquad\qquad\displaystyle\int_a^bf(x)dx\approx S_n=\displaystyle\sum_{k=0}^{n-1}S^{(k)}=\dfrac{h}{6}\displaystyle\sum_{k=0}^{n-1}\left[f(x_k)+4f\left(x_{k+\frac{1}{2}}\right)+f(x_{k+1})\right] ∫ab​f(x)dx≈Sn​=k=0∑n−1​S(k)=6h​k=0∑n−1​[f(xk​)+4f(xk+21​​)+f(xk+1​)]

误差 R n ( f ) = − b − a 180 ( h 2 ) 4 f ( 4 ) ( η ) , η ∈ ( a , b ) R_n(f)=-\dfrac{b-a}{180}\left(\dfrac{h}{2}\right)^4f^{(4)}(\eta),\quad\eta\in(a,b) Rn​(f)=−180b−a​(2h​)4f(4)(η),η∈(a,b) 为 h 4 h^4 h4 阶

实现代码

例：求积分 I = ∫ 0 1 sin ⁡ ( x ) x d x I=\displaystyle\int_0^1\dfrac{\sin(x)}{x}dx I=∫01​xsin(x)​dx， 7 7 7 位有效数字准确值为 0.9460831 0.9460831 0.9460831

import numpy as npdef trapezoid(f,a,b):    # trapezoid(n=1)return (b-a)*(f(b)+f(a))/2def simpson(f,a,b):    # simpson(n=2)t = (a+b)/2val = (b-a)*(f(a)+4*f(t)+f(b))/6return valdef compositeInt(f,a,b,n,mode='trapezoid'):h = (b-a)/nval = 0eps = 1e-10if mode=='trapezoid':for k in range(n):xk = a+k*h+epsxk1 = xk+hval += trapezoid(f,xk,xk1)        if mode=='simpson':for k in range(n):xk = a+k*h+epsxk1 = xk+hval += simpson(f,xk,xk1)            return valif __name__ == '__main__':a = 0b = 1f = lambda x: np.sin(x)/x    print('====== trapezoid: newton-cotes(1st-order) ======')val = compositeInt(f, a, b, 8)print('trapezoid: %.7f' % (val))     print('====== simpson: newton-cotes(2nd-order) ======')val = compositeInt(f, a, b, 4,'simpson')print('simpson: %.7f' % (val))

运行结果：
====== trapezoid: newton-cotes(1st-order) ======
　　　trapezoid: 0.9456909　　（2位有效数字）
====== simpson: newton-cotes(2nd-order) ======
　　　simpson: 0.9460833　　 （6位有效数字）
　　　
( 1 ) \qquad(1) (1) 由于 4 4 4 等分的复化辛普森方法在每个子区间上还需要再计算一个二分点，因此区间的等分数实际上是 8 8 8 等分。
( 2 ) \qquad(2) (2) 然而， 4 4 4 等分的复化辛普森方法的精度高于 8 8 8 等分的复化梯形方法。由误差公式也可计算出 ∣ R 8 ( f ) ∣ ≤ 0.000434 \vert R_8(f)\vert\le0.000434 ∣R8​(f)∣≤0.000434 以及 ∣ R 4 ( f ) ∣ ≤ 0.271 × 1 0 − 6 \vert R_4(f)\vert\le0.271\times10^{-6} ∣R4​(f)∣≤0.271×10−6。
\qquad

4. 求积公式的递推化

\qquad 复化求积方法如果要提高精度，就必然增加求积子区间的数量。例如，将 [ a , b ] [a,b] [a,b] 进行 n n n 等分时，求积节点的数量为 n + 1 n+1 n+1；若将每个求积区间 [ x k , x k + 1 ] [x_k,x_{k+1}] [xk​,xk+1​] 再一分为二，求积节点的数量变为 2 n + 1 2n+1 2n+1，每次等分之后，求积区间的数量就要翻倍。如果仍然在每个子区间上使用复化公式，无疑会极大地增加计算量。

\qquad 为了提高计算效率，可以将等分前后的积分值与进行比较，从而形成“递推关系”。
\qquad

4.1 复化梯形公式的递推化

\qquad 为了提高复化求积公式的精度，可以选择将每个求积子区间 [ x k , x k + 1 ] [x_k,x_{k+1}] [xk​,xk+1​] 一分为二的方式，也就在区间 [ x k , x k + 1 ] [x_k,x_{k+1}] [xk​,xk+1​] 增加一个等距的中间点 x k + 1 2 = 1 2 ( x k + x k + 1 ) x_{k+\frac{1}{2}}=\dfrac{1}{2}(x_k+x_{k+1}) xk+21​​=21​(xk​+xk+1​)。

\qquad

• 假设将 [ a , b ] [a,b] [a,b] 等分 n n n 个点 a , x 1 , ⋯ , x n − 2 , b a,x_1,\cdots,x_{n-2},b a,x1​,⋯,xn−2​,b，步长 h = b − a n h=\dfrac{b-a}{n} h=nb−a​，考虑复化梯形公式：

1 ) \qquad1) 1) 以求积子区间 [ x k , x k + 1 ] [x_k,x_{k+1}] [xk​,xk+1​] 为例， n n n 等分时的梯形公式为

T 1 = h 2 [ f ( x k ) + f ( x k + 1 ) ] \qquad\qquad T_1=\dfrac{h}{2}[f(x_k)+f(x_{k+1})] T1​=2h​[f(xk​)+f(xk+1​)]

2 ) \qquad2) 2) 进行 2 n 2n 2n 等分时，需要将每个求积子区间 [ x k , x k + 1 ] [x_k,x_{k+1}] [xk​,xk+1​] 一分为二，变成了两个等距的子区间 [ x k , x k + 1 2 ] [x_k,x_{k+\frac{1}{2}}] [xk​,xk+21​​] 和 [ x k + 1 2 , x k + 1 ] [x_{k+\frac{1}{2}},x_{k+1}] [xk+21​​,xk+1​]，在两个子区间上分别采用梯形公式。此时， 2 n 2n 2n 等分时的复化梯形公式为是：

T 2 = h 4 [ f ( x k ) + f ( x k + 1 2 ) ] + h 4 [ f ( x k + 1 2 ) + f ( x k + 1 ) ] = h 4 [ f ( x k ) + 2 f ( x k + 1 2 ) + f ( x k + 1 ) ] = h 4 [ f ( x k ) + f ( x k + 1 ) ] + h 2 f ( x k + 1 2 ) = 1 2 T 2 + h 2 f ( x k + 1 2 ) \qquad\qquad\begin{aligned}T_2&=\dfrac{h}{4}[f(x_k)+f(x_{k+\frac{1}{2}})]+\dfrac{h}{4}[f(x_{k+\frac{1}{2}})+f(x_{k+1})]\\ &=\dfrac{h}{4}[f(x_k)+2f(x_{k+\frac{1}{2}})+f(x_{k+1})]\\ &=\dfrac{h}{4}[f(x_k)+f(x_{k+1})]+\dfrac{h}{2}f(x_{k+\frac{1}{2}})\\ &=\dfrac{1}{2}T_2+\dfrac{h}{2}f(x_{k+\frac{1}{2}})\end{aligned} T2​​=4h​[f(xk​)+f(xk+21​​)]+4h​[f(xk+21​​)+f(xk+1​)]=4h​[f(xk​)+2f(xk+21​​)+f(xk+1​)]=4h​[f(xk​)+f(xk+1​)]+2h​f(xk+21​​)=21​T2​+2h​f(xk+21​​)​

• 考虑整个求积区间 [ a , b ] [a,b] [a,b]， n n n 等分时的复化梯形求积为：

T n = h 2 ∑ k = 0 n − 1 [ f ( x k ) + f ( x k + 1 ) ] \qquad\qquad T_n=\dfrac{h}{2}\displaystyle\sum_{k=0}^{n-1}[f(x_k)+f(x_{k+1})] Tn​=2h​k=0∑n−1​[f(xk​)+f(xk+1​)]

\qquad 按照上述规则， 2 n 2n 2n 等分时的复化梯形求积为：

T 2 n = h 4 ∑ k = 0 n − 1 [ f ( x k ) + f ( x k + 1 ) ] + h 2 ∑ k = 0 n − 1 f ( x k + 1 2 ) \qquad\qquad T_{2n}=\dfrac{h}{4}\displaystyle\sum_{k=0}^{n-1}[f(x_k)+f(x_{k+1})]+\dfrac{h}{2}\displaystyle\sum_{k=0}^{n-1}f(x_{k+\frac{1}{2}}) T2n​=4h​k=0∑n−1​[f(xk​)+f(xk+1​)]+2h​k=0∑n−1​f(xk+21​​)

\qquad
\qquad
\qquad 对比这两个式子，可得到 n n n 等分与 2 n 2n 2n 等分时采用复化梯形公式求积结果的关系：

T 2 n = 1 2 T n + h 2 ∑ k = 0 n − 1 f ( x k + 1 2 ) \qquad\qquad T_{2n}=\dfrac{1}{2}T_n+\dfrac{h}{2}\displaystyle\sum_{k=0}^{n-1}f(x_{k+\frac{1}{2}}) T2n​=21​Tn​+2h​k=0∑n−1​f(xk+21​​)

\qquad 也就是说，从 n n n 等分到 2 n 2n 2n 等分，不需要再在 2 n 2n 2n 个新的求积区间上使用复化梯形公式重新计算，只需要在 n n n 等分求积结果 T n T_n Tn​ 的基础上，加上 n n n 个新的求积节点的值 f ( x k + 1 2 ) f(x_{k+\frac{1}{2}}) f(xk+21​​) 即可。

实现代码

仍以求积分 I = ∫ 0 1 sin ⁡ ( x ) x d x I=\displaystyle\int_0^1\dfrac{\sin(x)}{x}dx I=∫01​xsin(x)​dx为例：

import numpy as np
def trapezoid(f,a,b):               # trapezoid(n=1)return (b-a)*(f(b)+f(a))/2def compositeInt_trape(f,a,b,n):    # recursive composite integrationeps = 1e-10val = trapezoid(f,a+eps,b)print('T1: %.7f' % (val))i = 1while i<=np.log2(n):h = (b-a)/(2**i)t = 0for k in range(1,2**(i-1)+1):xkc = a + (2*k-1)*ht += f(xkc)*h    val = val/2 + tprint('T%d: %.7f' % (2**i,val))i += 1return val    if __name__ == '__main__':a = 0b = 1f = lambda x: np.sin(x)/x    n = 3           # 区间等分次数（3次等分后分为8个等距区间）val = compositeInt_trape(f, a, b, 2**n)print('recursive trapezoid: %.7f' % (val))

运行结果：
　　　T1: 0.9207355
　　　T2: 0.9397933
　　　T4: 0.9445135
　　　T8: 0.9456909
　　recursive trapezoid: 0.9456909
　　　
\qquad 本例中等分了 3 3 3 次之后的子区间数为 2 3 = 8 2^3=8 23=8，和第 3 3 3 节例子中复化梯形公式的结果相同，但是计算量减少了。　　　
\qquad

4.2 龙贝格(Romberg)算法

Romberg \qquad\text{Romberg} Romberg算法，是另一种减少计算量的算法，也是基于递推的方式。

• 分析梯形公式的误差
  \qquad
  R n ( f ) = − b − a 12 h 2 f ′ ′ ( η ) , η ∈ ( a , b ) \qquad R_n(f)=-\dfrac{b-a}{12}h^2f^{''}(\eta),\quad\eta\in(a,b) Rn​(f)=−12b−a​h2f′′(η),η∈(a,b)
  \qquad
  \qquad 可以得到： { I − T n = − b − a 12 h 2 f ′ ′ ( η ) , η ∈ ( a , b ) I − T 2 n = − b − a 12 ( h 2 ) 2 f ′ ′ ( η ˉ ) , η ˉ ∈ ( a , b ) \left\{\begin{aligned} I-T_n&=-\dfrac{b-a}{12}h^2f^{''}(\eta),\qquad\eta\in(a,b) \\ \\ I-T_{2n}&=-\dfrac{b-a}{12}\left(\dfrac{h}{2}\right)^2f^{''}(\bar\eta),\qquad\bar\eta\in(a,b) \end{aligned} \right. ⎩ ⎨ ⎧​I−Tn​I−T2n​​=−12b−a​h2f′′(η),η∈(a,b)=−12b−a​(2h​)2f′′(ηˉ​),ηˉ​∈(a,b)​

\qquad
\qquad 若假设 f ′ ′ ( η ) ≈ f ′ ′ ( η ˉ ) f^{''}(\eta)\approx f^{''}(\bar\eta) f′′(η)≈f′′(ηˉ​)，那么
\qquad
I − T 2 n I − T n ≈ 1 4 ⟹ I − T 2 n ≈ 1 3 ( T 2 n − T n ) \qquad\qquad\dfrac{I-T_{2n}}{I-T_n}\approx\dfrac{1}{4}\Longrightarrow I-T_{2n}\approx\dfrac{1}{3}(T_{2n}-T_n) I−Tn​I−T2n​​≈41​⟹I−T2n​≈31​(T2n​−Tn​)
\qquad
\qquad 从而可以通过递推关系得到积分值估计： T ‾ = T 2 n + 1 3 ( T 2 n − T n ) = 4 3 T 2 n − 1 3 T n \overline T=T_{2n}+\dfrac{1}{3}(T_{2n}-T_n)=\dfrac{4}{3}T_{2n}-\dfrac{1}{3}T_n T=T2n​+31​(T2n​−Tn​)=34​T2n​−31​Tn​

\qquad 又由　 T n = h 2 ∑ k = 0 n − 1 [ f ( x k ) + f ( x k + 1 ) ] T_n=\dfrac{h}{2}\displaystyle\sum_{k=0}^{n-1}[f(x_k)+f(x_{k+1})] Tn​=2h​k=0∑n−1​[f(xk​)+f(xk+1​)]
\qquad 　和　 T 2 n = h 4 ∑ k = 0 n − 1 [ f ( x k ) + f ( x k + 1 ) ] + h 2 ∑ k = 0 n − 1 f ( x k + 1 2 ) T_{2n}=\dfrac{h}{4}\displaystyle\sum_{k=0}^{n-1}[f(x_k)+f(x_{k+1})]+\dfrac{h}{2}\displaystyle\sum_{k=0}^{n-1}f(x_{k+\frac{1}{2}}) T2n​=4h​k=0∑n−1​[f(xk​)+f(xk+1​)]+2h​k=0∑n−1​f(xk+21​​)

\qquad 可得到积分值估计值 T ‾ \overline T T 实际上就是辛普森公式的积分值 S n S_n Sn​：

T ‾ = 4 3 T 2 n − 1 3 T n = h 6 ∑ k = 0 n − 1 [ f ( x k + f ( x k + 1 ) ] + 2 h 3 ∑ k = 0 n − 1 f ( x k + 1 2 ) = h 6 ∑ k = 0 n − 1 [ f ( x k ) + 4 f ( x k + 1 2 ) + f ( x k + 1 ) ] = S n \qquad\qquad\begin{aligned}\overline T&=\dfrac{4}{3}T_{2n}-\dfrac{1}{3}T_n\\ &=\dfrac{h}{6}\displaystyle\sum_{k=0}^{n-1}[f(x_k+f(x_{k+1})]+\dfrac{2h}{3}\displaystyle\sum_{k=0}^{n-1}f(x_{k+\frac{1}{2}})\\ &=\dfrac{h}{6}\displaystyle\sum_{k=0}^{n-1}[f(x_k)+4f(x_{k+\frac{1}{2}})+f(x_{k+1})]=S_n\end{aligned} T​=34​T2n​−31​Tn​=6h​k=0∑n−1​[f(xk​+f(xk+1​)]+32h​k=0∑n−1​f(xk+21​​)=6h​k=0∑n−1​[f(xk​)+4f(xk+21​​)+f(xk+1​)]=Sn​​
\qquad

• 依次类推，考虑（2等分）辛普森公式的误差

I − S 2 n I − S n ≈ 1 16 ⟹ I − S 2 n ≈ 1 16 ( S 2 n − S n ) \qquad\qquad\dfrac{I-S_{2n}}{I-S_n}\approx\dfrac{1}{16}\Longrightarrow I-S_{2n}\approx\dfrac{1}{16}(S_{2n}-S_n) I−Sn​I−S2n​​≈161​⟹I−S2n​≈161​(S2n​−Sn​)
\qquad
\qquad 从而可以通过递推关系得到积分值估计： S ‾ = 16 15 S 2 n − 1 15 S n \overline S= \dfrac{16}{15}S_{2n}-\dfrac{1}{15}S_n S=1516​S2n​−151​Sn​

\qquad 验证可知积分值估计值 S ‾ \overline S S 就是科特斯公式的积分值： C n = 16 15 S 2 n − 1 15 S n C_n= \dfrac{16}{15}S_{2n}-\dfrac{1}{15}S_n Cn​=1516​S2n​−151​Sn​
\qquad

• 依次类推，考虑（4等分）科特斯公式的误差

I − C 2 n I − C n ≈ 1 64 ⟹ I − C 2 n ≈ 1 64 ( C 2 n − C n ) \qquad\qquad\dfrac{I-C_{2n}}{I-C_n}\approx\dfrac{1}{64}\Longrightarrow I-C_{2n}\approx\dfrac{1}{64}(C_{2n}-C_n) I−Cn​I−C2n​​≈641​⟹I−C2n​≈641​(C2n​−Cn​)
\qquad
\qquad 从而可以通过递推关系得到积分值估计： C ‾ = 64 63 C 2 n − 1 63 C n \overline C=\dfrac{64}{63}C_{2n}-\dfrac{1}{63}C_n C=6364​C2n​−631​Cn​

\qquad 也就得到了龙贝格 ( Romberg ) (\text{Romberg}) (Romberg)公式：

R n = 64 63 C 2 n − 1 63 C n \qquad\qquad R_n=\dfrac{64}{63}C_{2n}-\dfrac{1}{63}C_n Rn​=6364​C2n​−631​Cn​

\qquad 使用龙贝格公式时，需要首先将求积区间 [ a , b ] [a,b] [a,b] 等分 2 3 = 8 2^3=8 23=8 次。在这些子区间上，采用复化梯形公式进行求积运算。其具体过程可以表示成：

k k k	h h h	T 2 k T_{2^k} T2k​	S 2 k − 1 S_{2^{k-1}} S2k−1​	C 2 k − 2 C_{2^{k-2}} C2k−2​	R 2 k − 3 R_{2^{k-3}} R2k−3​
0 0 0	b − a b-a b−a	T 1 = 0.92073549 T_1=0.92073549 T1​=0.92073549
1 1 1	b − a 2 \frac{b-a}{2} 2b−a​	T 2 = 0.93979328 T_2=0.93979328 T2​=0.93979328	S 1 S_1 S1​
2 2 2	b − a 4 \frac{b-a}{4} 4b−a​	T 4 = 0.94451352 T_4=0.94451352 T4​=0.94451352	S 2 S_2 S2​	C 1 C_1 C1​
3 3 3	b − a 8 \frac{b-a}{8} 8b−a​	T 8 = 0.94569086 T_8=0.94569086 T8​=0.94569086	S 4 S_4 S4​	C 2 C_2 C2​	R 1 R_1 R1​

表中， k k k 表示等分次数， h h h 表示子区间长度，等分后的子区间数为 2 k 2^k 2k
　　　 T 2 k T_{2^k} T2k​ 表示 k k k 次等分时，采用复化梯形公式求积分的值（详见4.1节代码运行结果）
从第 4 4 4 列开始，就是采用龙贝格递推算法的加速过程，不需要再进行复化求积。

4.3 理查森(Richardson)外推算法

\qquad 龙贝格算法实际上是理查森 (Richardson) \text{(Richardson)} (Richardson)外推算法的一种特殊情况。

梯形公式的余项，可以展开成级数形式( α l \alpha_l αl​ 与 h h h 无关)：
T ( h ) = I + α 1 h 2 + α 2 h 4 + ⋯ + α l h 2 l + ⋯ T(h)=I+\alpha_1h^2+\alpha_2h^4+\cdots+\alpha_lh^{2l}+\cdots T(h)=I+α1​h2+α2​h4+⋯+αl​h2l+⋯
T ( h 2 ) = I + α 1 h 2 4 + α 2 h 4 16 + ⋯ + α l ( h 2 ) 2 l + ⋯ T(\frac{h}{2})=I+\alpha_1\frac{h^2}{4}+\alpha_2\frac{h^4}{16}+\cdots+\alpha_l(\frac{h}{2})^{2l}+\cdots T(2h​)=I+α1​4h2​+α2​16h4​+⋯+αl​(2h​)2l+⋯
⟹ T 1 ( h ) = 4 T ( h 2 ) − T ( h ) 3 = I + β 1 h 4 + β 2 h 6 + ⋯ \Longrightarrow T_1(h)=\dfrac{4T(\frac{h}{2})-T(h)}{3}=I+\beta_1h^4+\beta_2h^6+\cdots ⟹T1​(h)=34T(2h​)−T(h)​=I+β1​h4+β2​h6+⋯ （辛普森公式序列）
T 1 ( h 2 ) = I + β 1 h 4 16 + β 2 h 6 64 + ⋯ \qquad T_1(\frac{h}{2})=I+\beta_1\frac{h^4}{16}+\beta_2\frac{h^6}{64}+\cdots T1​(2h​)=I+β1​16h4​+β2​64h6​+⋯
⟹ T 2 ( h ) = 16 T 1 ( h 2 ) − T 1 ( h ) 15 = I + γ 1 h 6 + γ 2 h 8 + ⋯ \qquad\Longrightarrow T_2(h)=\dfrac{16T_1(\frac{h}{2})-T_1(h)}{15}=I+\gamma_1h^6+\gamma_2h^8+\cdots ⟹T2​(h)=1516T1​(2h​)−T1​(h)​=I+γ1​h6+γ2​h8+⋯ （科特斯公式序列）
⋯ ⋯ \qquad\qquad\cdots\cdots ⋯⋯

\qquad 将龙贝格算法进行扩展，并进行公式化，可得到理查森 (Richardson) \text{(Richardson)} (Richardson)递推关系：

T m ( k ) = 4 m 4 m − 1 T m − 1 ( k + 1 ) − 1 4 m − 1 T m − 1 ( k ) , k = 1 , 2 , ⋯ \qquad\qquad T_m^{(k)}=\dfrac{4^m}{4^m-1}T_{m-1}^{(k+1)}-\dfrac{1}{4^m-1}T_{m-1}^{(k)}\ ,\quad k=1,2,\cdots Tm(k)​=4m−14m​Tm−1(k+1)​−4m−11​Tm−1(k)​ ,k=1,2,⋯

\qquad 特别地，当选择 m = 3 m=3 m=3 时，即为龙贝格算法：

( 1 ) \qquad(1) (1)　 m = 1 m=1 m=1 时， T 1 ( k ) = 4 3 T 0 ( k + 1 ) − 1 3 T 0 ( k ) T_1^{(k)}=\dfrac{4}{3}T_0^{(k+1)}-\dfrac{1}{3}T_0^{(k)} T1(k)​=34​T0(k+1)​−31​T0(k)​

( 2 ) \qquad(2) (2)　 m = 2 m=2 m=2 时， T 2 ( k ) = 16 15 T 1 ( k + 1 ) − 1 15 T 1 ( k ) T_2^{(k)}=\dfrac{16}{15}T_1^{(k+1)}-\dfrac{1}{15}T_1^{(k)} T2(k)​=1516​T1(k+1)​−151​T1(k)​

( 3 ) \qquad(3) (3)　 m = 3 m=3 m=3 时， T 3 ( k ) = 64 63 T 2 ( k + 1 ) − 1 63 T 2 ( k ) T_3^{(k)}=\dfrac{64}{63}T_2^{(k+1)}-\dfrac{1}{63}T_2^{(k)} T3(k)​=6364​T2(k+1)​−631​T2(k)​

\qquad 当 m > 3 m>3 m>3 时，即为理查森 (Richardson) \text{(Richardson)} (Richardson)外推算法，如下表所示：

k k k	h h h	T 0 ( k ) ( T 2 k ) T_0^{(k)}\ (T_{2^k}) T0(k)​ (T2k​)	T 1 ( k ) ( S 2 k − 1 ) T_1^{(k)}\ (S_{2^{k-1}}) T1(k)​ (S2k−1​)	T 2 ( k ) ( C 2 k − 2 ) T_2^{(k)}\ (C_{2^{k-2}}) T2(k)​ (C2k−2​)	T 3 ( k ) ( R 2 k − 3 ) T_3^{(k)}\ (R_{2^{k-3}}) T3(k)​ (R2k−3​)	T 4 ( k ) T_4^{(k)} T4(k)​
0 0 0	b − a b-a b−a	T 0 ( 1 ) : T 1 = 0.92073549 T_0^{(1)}:T_1=0.92073549 T0(1)​:T1​=0.92073549
1 1 1	b − a 2 \frac{b-a}{2} 2b−a​	T 0 ( 2 ) : T 2 = 0.93979328 T_0^{(2)}:T_2=0.93979328 T0(2)​:T2​=0.93979328	T 1 ( 1 ) : S 1 T_1^{(1)}:S_1 T1(1)​:S1​
2 2 2	b − a 4 \frac{b-a}{4} 4b−a​	T 0 ( 3 ) : T 4 = 0.94451352 T_0^{(3)}:T_4=0.94451352 T0(3)​:T4​=0.94451352	T 1 ( 2 ) : S 2 T_1^{(2)}:S_2 T1(2)​:S2​	T 2 ( 1 ) : C 1 T_2^{(1)}:C_1 T2(1)​:C1​
3 3 3	b − a 8 \frac{b-a}{8} 8b−a​	T 0 ( 4 ) : T 8 = 0.94569086 T_0^{(4)}:T_8=0.94569086 T0(4)​:T8​=0.94569086	T 1 ( 3 ) : S 4 T_1^{(3)}:S_4 T1(3)​:S4​	T 2 ( 2 ) : C 2 T_2^{(2)}:C_2 T2(2)​:C2​	T 3 ( 1 ) : R 1 T_3^{(1)}:R_1 T3(1)​:R1​
4 4 4	b − a 16 \frac{b-a}{16} 16b−a​	T 0 ( 5 ) : T 16 = 0.94598503 T_0^{(5)}:T_{16}=0.94598503 T0(5)​:T16​=0.94598503	T 1 ( 4 ) : S 8 T_1^{(4)}:S_8 T1(4)​:S8​	T 2 ( 3 ) : C 4 T_2^{(3)}:C_4 T2(3)​:C4​	T 3 ( 2 ) : R 2 T_3^{(2)}:R_2 T3(2)​:R2​	T 4 ( 4 ) T_4^{(4)} T4(4)​
⋮ \vdots ⋮	⋮ \vdots ⋮	⋮ \vdots ⋮	⋮ \vdots ⋮	⋮ \vdots ⋮	⋮ \vdots ⋮	⋮ \vdots ⋮

从第 4 4 4 列开始，每个表项的值满足其左侧、左上侧的递推关系。

实现代码

仍以求积分 I = ∫ 0 1 sin ⁡ ( x ) x d x I=\displaystyle\int_0^1\dfrac{\sin(x)}{x}dx I=∫01​xsin(x)​dx为例：

import numpy as npdef trapezoid(f,a,b):    # trapezoid(n=1)return (b-a)*(f(b)+f(a))/2def compositeInt_trape(f,a,b,n):    #recursiveeps = 1e-10num = np.log2(n).astype(int)val = np.zeros(num+1)val[0] = trapezoid(f,a+eps,b)i = 1while i<=np.log2(n):h = (b-a)/(2**i)t = 0for k in range(1,2**(i-1)+1):xkc = a + (2*k-1)*ht += f(xkc)*h    val[i] = val[i-1]/2 + ti += 1return valdef romberg(f,a,b):# 只需要等分3次val_t = compositeInt_trape(f, a, b, 2**3)    t1 = np.concatenate((np.ones(1),val_t[0:-1]))t2 = (4*val_t - t1)/3val_s = t2[1:]    t1 = np.concatenate((np.ones(1),val_s[0:-1]))t2 = (16*val_s - t1)/15val_c = t2[1:]    t1 = np.concatenate((np.ones(1),val_c[0:-1]))t2 = (64*val_c - t1)/63val_r = t2[1]    print(val_t[-1::-1].flatten())print(val_s[-1::-1].flatten())print(val_c[-1::-1].flatten())print('[%.8f]' % (val_r))    return val_rdef richardson(f,a,b,k):n = k+2val_t = compositeInt_trape(f, a, b, 2**n)print(val_t[-1::-1].flatten())    for i in range(1,n):t0 = 4**it1 = np.concatenate((np.ones(1),val_t[0:-1]))t2 = (t0*val_t - t1)/(t0-1)val_t = t2[1:]print(val_t[-1::-1].flatten())print('[%.8f]' % (val_t[1]))    return val_t[1]if __name__ == '__main__':a = 0b = 1f = lambda x: np.sin(x)/x    val = romberg(f, a, b)print('romberg: %.7f' % (val))    val = richardson(f, a, b, 2) # 比Romberg多等分一次print('richardson: %.7f' % (val))

运行结果：
[0.94569086 0.94451352 0.93979328 0.92073549]
[0.94608331 0.94608693 0.94614588]
[0.94608307 0.946083 ]
[0.94608307]
romberg: 0.9460831

只需要进行 2 3 = 8 2^3=8 23=8 等分，并分别计算复合梯形公式的结果 T 1 , T 2 , T 4 , T 8 T_1,T_2,T_4,T_8 T1​,T2​,T4​,T8​ ( T 8 T_8 T8​ 只具有2位有效数字)，就可以递推出具有更高精度的结果 T 3 ( 1 ) T_3^{(1)} T3(1)​ (具有7位有效数字)。

[0.94598503 0.94569086 0.94451352 0.93979328 0.92073549]
[0.94608309 0.94608331 0.94608693 0.94614588]
[0.94608307 0.94608307 0.946083 ]
[0.94608307 0.94608307]
[0.94608307]
richardson: 0.9460831

参考：
数值微分的python实现

数值积分的python实现——NewtonCotes、复化求积、Romberg、richardson递推相关推荐

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

   数值分析-数值积分 前言 敲黑板 牛顿-莱布尼兹公式的局限性 常见的近似求法 代数精度 插值型求积公式 Newton-Cotes公式 复化求积 Gauss 求积公式 基于正交多项式的Gauss求积公式 ...

2. 数值分析复化求积matlab,MATLAB数值分析实验二(复合梯形、辛普森和龙贝格求积,以及二重积分计算等)...

   1.理解如何在计算机上使用数值方法计算定积近似值; 2.学会复合梯形.复合Simpson和龙贝格求积分公式的编程与应用. 3.探索二重积分在矩形区域的数值积分方法. 佛山科学技术学院 实 验 报 告 ...

3. 变步长梯形公式数值积分的Python程序

   基本代码: from sympy import *def f(t):f = 2000*log(140000/(140000-2100*t))-9.8*treturn fx = symbols('x') ...

4. 梯形公式的数值积分的Python程序

   代码 例子 下面的代码的意图是改变积分区域的梯形块数,得到不同精度的积分值. 以对函数f(t) = 2000×ln(140000/(140000-2100×t))-9.8×t 进行积分为例,积分区域为 ...

5. python ftp文件夹文件递归上传推送

   python ftp文件夹文件递归上传推送 posted on 2018-10-16 17:05 秦瑞It行程实录 阅读(...) 评论(...) 编辑 收藏 转载于:https://www.cnbl ...

6. 2.3 基本算法之递归变递推 1188 菲波那契数列(2) python

   http://noi.openjudge.cn/ch0203/1760/ """2.3 基本算法之递归变递推 1188 菲波那契数列(2)--3分 http://ybt. ...

7. 2.3 基本算法之递归变递推 放苹果 python

   http://noi.openjudge.cn/ch0203/666/ """ 2.3 基本算法之递归变递推 666 放苹果 http://noi.openjudge.c ...

8. 钉钉一行代码_利用Python快速搭建钉钉和邮件数据推送系统

   前面的文章我们写到了利用Python实现钉钉和邮件的数据推送,在数据处理这一块实现了对mysql和odps的数据获取和处理,可以满足常规业务大部分数据场景需求,在一家初创公司数据基础建设还不完善的时候 ...

9. 蓝桥杯python省赛冲刺篇2——常用算法的详细解析及对应蓝桥杯真题：打表模拟法、递推递归法、枚举法、贪心算法、差分与前缀和

   注意:加了题目链接 目录 注意:加了题目链接 一.打表模拟法 介绍 1. 算式问题 题目描述 解析与代码演示 2. 求值 题目描述 解析与代码演示 3. 既约分数 题目描述 解析与代码演示 4. 天干 ...

最新文章

1. Intellij IDEA使用教程（超详细）
2. js中深拷贝和浅拷贝问题
3. Java 异常处理的误区和经验总结--转载
4. 转]Window, Linux动态链接库的分析对比
5. docker+jenkins+maven+svn
6. Markdown编译器插入公式的数学符号及字体颜色、背景
7. (转) Spring 3 报org.aopalliance.intercept.MethodInterceptor问题解决方法
8. java处理请求的流程_Java Spring mvc请求处理流程详解
9. Navicat打开保存的查询
10. 是什么意思网络语_网络语“随薪锁欲”是什么意思？
11. 无线路由器挖洞方法大比拼：白盒 or 黑盒？
12. 自我总结（五）---(学习j2ee)
13. Docker下安装Anaconda
14. WCF入门（七）——异常处理1
15. layui动态设置checkbox选中状态
16. 开源的shell工具finalShell
17. itunes安装后不能用，双击后等很长时间，提示：ITUNES 驱动程序缺少用于导入和刻录的CD与DVD注册的设置...
18. 六、量子纠错码的构成法
19. python高级函数_python高级之函数
20. kali之beef的使用

热门文章

1. 〖Python自动化办公篇⑳〗 - python实现邮件自动化 - 发送html邮件和带附件的邮件
2. 看伟大的领袖如何激励行动有感
3. 天气预报接口，在自己的网页上加入天气预报吧
4. “你家娃为什么这么爱看书？”只用3招，孩子秒变小书迷
5. 2020起重机械指挥考试题库及起重机械指挥考试试题
6. 使用springboot+netty处理tcp/ip服务端编程
7. html 人物介绍 轮播,jQuery卡通人物介绍卡牌轮播切换代码
8. 【MMU篇】初见MMU和TLB
9. Mac 之 Command Line Tools
10. 经纬恒润受邀出席首届中国汽车芯片应用创新拉力赛新闻发布会