■ 前言

在 Composite Numerical Integration 中给出了三种复合数值积分方法它们分别是：

Newton-Cotes Formulas

01三种方法的公式及其Python程序

1.复合梯形方法

假设 $f \in C^2 \left[ {a,b} \right]$， $h = \left( {b - a} \right)/n$， $x_j = a + j \cdot h$， $j = 0,1, \cdots ,n$：

其中存在 $\mu \in \left( {a,b} \right)$，使得上述积分的误差等于：
$E = {{b - 2} \over {12}}h^2 f''\left( \mu \right)$

在 Laplace数值逆运算的讨论 中给出了复合梯形方法的Python程序实现：

def trapz(f, a, b, N=50):
    x = linspace(a, b, N+1)
    y = f(x)
    y_right = y[1:]
    y_left = y[:-1]
    dx = (b-a) / N
    T = dx/2 * sum(y_right + y_left)
    return T

2.Simpson1/3法则

对于积分函数的要求以及对区间的划分与前面复合梯形方法是一样的，只是对于区间划分的数量 $n$要求必须是偶数。

积分误差与区间长度 $h$之间是五次方的关系。

def simpson(f, a, b, N=50):
    x = linspace(a, b, N+1, endpoint=True)
    y = f(x)
    y4 = y[1:-1:2]
    y2 = y[2:-2:2]

    dx = (b-a) / N
    T = dx/3*(sum(y2)*2 + sum(y4)*4 + y[0] + y[-1])
    return T

3.Simpson3/8法则

对于区间的划分与前面相同，只是对区间划分的数量 $n$为3的整数倍数。那么：

def simpson3_8(f, a, b, N=50):
    x = linspace(a, b, N+1, endpoint=True)
    y = f(x)
    y0 = y[0:-1:3]
    y3 = y[3::3]
    y1 = y[1:-1:3]
    y2 = y[2:-1:3]

    dx = (b-a) / N
    T = dx*3/8*(sum(y0) + sum(y1)*3 + sum(y2)*3 + sum(y3))
    return T

4.对比实验

N=6时，所得到的结果如下：

printf(simpson1_3(sin, 0, pi, 6))
printf(simpson3_8(sin, 0, pi, 6))
printf(trapz(sin, 0, pi/2, 6)