数值积分中的辛普森方法及其误差估计

$f$在$(x_0,x_2)$上四阶可导，且在$[x_0,x_2]$上三阶导函数连续.则
\begin{equation}
\int_{x_0}^{x_2}f(x)dx=\frac{h}{3}[f(x_0)+4f(x_1)+f(x_2)]-\frac{h^5}{90}f^{(4)}(\xi)
\end{equation}其中$h=x_1-x_0=x_2-x_1$.

证明:我们仍然进行牛顿插值.设牛顿插值的插值点为
\begin{equation}
x_0,x_1,x_2
\end{equation}
则经过这三个插值点的牛顿插值多项式为
\begin{equation}
f(x_0)+(x-x_0)f[x_0,x_1]+(x-x_0)(x-x_1)f[x_0,x_1,x_2]
\end{equation}
插值余项为
\begin{equation}
(x-x_0)(x-x_1)(x-x_2)f[x,x_0,x_1,x_2]
\end{equation}
我们先来看
\begin{equation}
\int_{x_0}^{x_2}\{f(x_0)+(x-x_0)f[x_0,x_1]+(x-x_0)(x-x_1)f[x_0,x_1,x_2]\}dx
\end{equation}
我们先来算
\begin{equation}
\int_{x_0}^{x_2}f(x_0)dx=2hf(x_0)
\end{equation}
\begin{align*}
\int_{x_0}^{x_2}(x-x_0)f[x_0,x_1]dx=\frac{1}{2}(x_2-x_0)^2f[x_0,x_1]=\frac{1}{2}(2h)^2 \frac{f(x_1)-f(x_0)}{h}=2h(f(x_1)-f(x_0))
\end{align*}
\begin{align*}
\int_{x_0}^{x_2}(x-x_0)(x-x_1)f[x_0,x_1,x_2]dx&=f[x_0,x_1,x_2]\int_{x_0}^{x_2}(x-x_0)(x-x_1)dx\\&=f[x_0,x_1,x_2]\int_{x_0}^{x_2}(x^2-(x_0+x_1)x+x_0x_1)dx\\&=f[x_0,x_1,x_2][(\frac{1}{3}x_2^3-\frac{1}{2}(x_0+x_1)x_2^2+x_0x_1x_2)-(\frac{1}{3}x_0^3-\frac{1}{2}(x_0+x_1)x_0^2+x_0^2x_1)]\\&=\frac{1}{6}f[x_0,x_1,x_2](x_2-x_0)^2(2x_2+x_0-3x_1)
\end{align*}
而
\begin{align*}
f[x_0,x_1,x_2]=\frac{f(x_0)}{2h^2}+\frac{f(x_1)}{-h^2}+\frac{f(x_2)}{2h^2}
\end{align*}
可见，
\begin{align*}
\int_{x_0}^{x_2}(x-x_0)(x-x_1)f[x_0,x_1,x_2]dx=\frac{[f(x_0)-2f(x_1)+f(x_2)]h}{3}
\end{align*}
可见，
\begin{align*}
\int_{x_0}^{x_2}\{f(x_0)+(x-x_0)f[x_0,x_1]+(x-x_0)(x-x_1)f[x_0,x_1,x_2]\}dx=\frac{(f(x_0)+4f(x_1)+f(x_2))h}{3}
\end{align*}
下面我来证明
\begin{equation}\label{eq:12345}
\int_{x_0}^{x_2} (x-x_0)(x-x_1)(x-x_2)f[x,x_0,x_1,x_2]dx=\frac{-h^5}{90}f^{(4)}(\xi)
\end{equation}
令
\begin{align*}
g(t)=\int_{x_0}^t(x-x_0)(x-x_1)(x-x_2)f[x,x_0,x_1,x_2]dx
\end{align*}

则
\begin{align*}
g'(t)=(t-x_0)(t-x_1)(t-x_2)f[t,x_0,x_1,x_2]
\end{align*}
则
\begin{align*}
g'(x_0)=g'(x_1)=g'(x_2)=0,g(x_0)=0
\end{align*}
我们进行Hermite插值.由于Hermite插值是牛顿插值的极限情形，为此我们先进行牛顿插值.我们设立插值点
\begin{align*}
x_0,x_0',x_1,x_1',x_2,x_2'
\end{align*}
经过这几个插值点的牛顿插值多项式为
\begin{align*}
p(x)= g(x_0)&+(x-x_0)g[x_0,x_0']\\&+(x-x_0)(x-x_0')g[x_0,x_0',x_1]\\&+(x-x_0)(x-x_0')(x-x_1)g[x_0,x_0',x_1,x_1']\\&+(x-x_0)(x-x_0')(x-x_1)(x-x_1')g[x_0,x_0',x_1,x_1',x_2]\\&+(x-x_0)(x-x_0')(x-x_1)(x-x_1')(x-x_2)g[x_0,x_0',x_1,x_1',x_2,x_2']
\end{align*}
令$x_0'\to x_0,x_1'\to x_1,x_2'\to x_2$可得相应的Hermite插值多项式为
\begin{align*}
q(x)=g(x_0)&+(x-x_0)g[x_0,x_0]\\&+(x-x_0)^2g[x_0,x_0,x_1]\\&+(x-x_0)^2(x-x_1)g[x_0,x_0,x_1,x_1]\\&+(x-x_0)^2(x-x_1)^2g[x_0,x_0,x_1,x_1,x_2]\\&+(x-x_0)^2(x-x_1)^2(x-x_2)g[x_0,x_0,x_1,x_1,x_2,x_2]
\end{align*}
令
\begin{align*}
k(x)=g(x)-q(x)
\end{align*}
可得
\begin{align*}
k(x_0)=k(x_1)=k(x_2)=0,k'(x_0)=k'(x_1)=k'(x_2)=0
\end{align*}
于是使用好几次Rolle定理我们可得
\begin{align*}
k^{(5)}(\xi)=0
\end{align*}
即
\begin{align*}
g^{(5)}(\xi)=5!g[x_0,x_0,x_1,x_1,x_2,x_2]
\end{align*}

易得
\begin{align*}
g^{(5)}(\xi)=f^{(4)}(\xi)
\end{align*}
因此
\begin{align*}
f^{(4)}(\xi)=5!g[x_0,x_0,x_1,x_1,x_2,x_2]
\end{align*}
下面来计算
\begin{align*}
g[x_0,x_0,x_1,x_1,x_2,x_2]
\end{align*}

\begin{align*}
g[x_0,x_0,x_1,x_1,x_2,x_2]&=\frac{g[x_0,x_0,x_1,x_1,x_2]-g[x_0,x_1,x_1,x_2,x_2]}{x_0-x_2}\\
&=\frac{\frac{g[x_0,x_0,x_1,x_1]-g[x_0,x_1,x_1,x_2]}{x_0-x_2}-\frac{g[x_0,x_1,x_1,x_2]-g[x_1,x_1,x_2,x_2]}{x_0-x_2}}{x_0-x_2}\\
&=\frac{g[x_0,x_0,x_1,x_1]-2g[x_0,x_1,x_1,x_2]+g[x_1,x_1,x_2,x_2]}{(x_0-x_2)^2}\\
&=\frac{\frac{g[x_0,x_0,x_1]-g[x_0,x_1,x_1]}{x_0-x_1}-2
\frac{g[x_0,x_1,x_1]-g[x_1,x_1,x_2]}{x_0-x_2}+\frac{g[x_1,x_1,x_2]-g[x_1,x_2,x_2]}{x_1-x_2}}{(x_0-x_2)^2}\\&=\frac{g[x_0,x_0,x_1]-2g[x_0,x_1,x_1]+2g[x_1,x_1,x_2]-g[x_1,x_2,x_2]}{(x_0-x_2)^2(x_0-x_1)}\\&=\frac{3\frac{g(x_1)}{(x_0-x_1)^2}-3
\frac{g(x_1)-g(x_2)}{(x_1-x_2)^2}}{(x_0-x_2)^2(x_0-x_1)}\\&=\frac{3g(x_2)}{(x_1-x_2)^2(x_0-x_2)^2(x_0-x_1)}\\&=\frac{3g(x_2)}{-4h^5}
\end{align*}
可见，
\begin{align*}
-4h^5\frac{f^{(4)}(\xi)}{5!\times 3}=g(x_2)
\end{align*}
即
\begin{align*}
g(x_2)=\frac{-h^{5}}{90}f^{(4)}(\xi)
\end{align*}

即\ref{eq:12345}成立.

注1:在做这个问题的过程中经历了一些曲折.现在把曲折记录如下:

为了证明

\begin{equation}
\int_{x_0}^{x_2} (x-x_0)(x-x_1)(x-x_2)f[x,x_0,x_1,x_2]dx=\frac{-h^5}{90}f^{(4)}(\xi)
\end{equation}

即证明

\begin{align*}
\int_{x_0}^{x_1}(x-x_0)(x-x_1)(x-x_2)f[x,x_0,x_1,x_2]dx+\int_{x_1}^{x_2}(x-x_0)(x-x_1)(x-x_2)f[x,x_0,x_1,x_2]dx=\frac{-h^5}{90}f^{(4)}(\xi)
\end{align*}

 \begin{align*}&\int_{x_0}^{x_1}(x-x_0)(x-x_1)(x-x_2)dx=\int_{x_0}^{x_1}[x^3-(x_0+x_1+x_2)x^2+(x_0x_1+x_0x_2+x_1x_2)x-x_0x_1x_2]dx\\&=(\frac{1}{4}x_1^4-\frac{1}{3}(x_0+x_1+x_2)x_{1}^3+\frac{1}{2}(x_0x_1+x_1x_2+x_2x_0)x_{1}^2-x_0x_1x_2x_{1})\\&-(\frac{1}{4}x_0^4-\frac{1}{3}(x_0+x_1+x_2)x_{0}^3+\frac{1}{2}(x_0x_1+x_1x_2+x_2x_0)x_{0}^2-x_0x_1x_2x_{0})\\&=\frac{1}{12}[3(x_1^2+x_0^2)(x_1+x_0)(x_1-x_0)\\&+4(x_0+x_1+x_2)(x_0-x_1)(x_0^2+x_1^2+x_0x_1)+6(x_0x_1+x_1x_2+x_2x_0)(x_1-x_0)(x_1+x_0)+12x_0x_1x_2(x_0-x_1)]\\&=\frac{1}{12}[(x_1-x_0)(3x_1^3+3x_1^2x_0+3x_0^2x_1+3x_0^3-4x_0^3-4x_0x_1^2-4x_0^2x_1-4x_1x_0^2-4x_1^3\\&-4x_0x_1^2-4x_2x_0^2-4x_2x_1^2-4x_0x_1x_2+6x_0x_1^2+6x_0^2x_1+6x_1^2x_2+6x_0x_1x_2+6x_0x_1x_2+6x_0^2x_2-12x_0x_1x_2)]\\&=\frac{1}{12}(x_1-x_0)(-x_1^3-x_0^3+x_1^2x_0+x_0^2x_1+2x_0^2x_2+2x_1^2x_2-4x_0x_1x_2)\\&=\frac{-1}{12}(x_1-x_0)(x_1^3+x_0^3-x_1^2x_0-x_0^2x_1-2x_0^2x_2-2x_1^2x_2+4x_0x_1x_2)\\&=\frac{-1}{12}(x_1-x_0)[x_1^2(x_1-x_0)-x_0^2(x_1-x_0)-2x_2(x_0^2+x_1^2-2x_0x_1)]\\&=\frac{-1}{12}(x_1-x_0)^2(x_1^2-x_0^2-2x_2(x_1-x_0))=-\frac{1}{12}(x_1-x_0)^3(x_1+x_0-2x_2)\end{align*}

我们再来计算
\begin{align*}
\int_{x_1}^{x_2}(x-x_0)(x-x_1)(x-x_2)dx&=\int_{x_1}^{x_2}[x^3-(x_0+x_1+x_2)x^2+(x_0x_1+x_0x_2+x_1x_2)x-x_0x_1x_2]dx\\&=(\frac{1}{4}x_2^4-\frac{1}{3}(x_0+x_1+x_2)x_{2}^3+\frac{1}{2}(x_0x_1+x_1x_2+x_2x_0)x_{2}^2-x_0x_1x_2x_{2})\\&-(\frac{1}{4}x_1^4-\frac{1}{3}(x_0+x_1+x_2)x_{1}^3+\frac{1}{2}(x_0x_1+x_1x_2+x_2x_0)x_{1}^2-x_0x_1x_2x_{1})
\end{align*}
根据对称性，容易得到它等于
\begin{align*}
-\frac{1}{12}(x_2-x_1)^3(x_1+x_2-2x_0)
\end{align*}

于是根据加权的积分中值定理，可得

\begin{align*}
&\int_{x_0}^{x_1}(x-x_0)(x-x_1)(x-x_2)f[x,x_0,x_1,x_2]dx+\int_{x_1}^{x_2}(x-x_0)(x-x_1)(x-x_2)f[x,x_0,x_1,x_2]dx\\&=-f[c_1,x_0,x_1,x_2]\frac{1}{12}(x_2-x_1)^3(x_1+x_2-2x_0)-f[c_2,x_0,x_1,x_2]\frac{1}{12}(x_1-x_0)^3(x_0+x_1-2x_2)
\end{align*}
而
\begin{align*}
\frac{1}{12}(x_2-x_1)^3(x_1+x_2-2x_0)=-\frac{1}{12}(x_1-x_0)^3(x_0+x_1-2x_2)&=\frac{1}{4}h^4
\end{align*}

因此
\begin{align*}
\int_{x_0}^{x_2}(x-x_0)(x-x_1)(x-x_2)f[x,x_0,x_1,x_2]dx=\frac{1}{4}h^4(f[c_2,x_0,x_1,x_2]-f[c_1,x_0,x_1,x_2])
\end{align*}

即
\begin{align*}
\int_{x_0}^{x_2}(x-x_0)(x-x_1)(x-x_2)f[x,x_0,x_1,x_2]dx=\frac{1}{4}h^4(c_2-c_1)f[c_1,c_2,x_0,x_1,x_2]
\end{align*}

令$c_2-c_1=\alpha h$.下面我们来确定$\alpha$.我们设

\begin{align*}
\int_{x_0}^{x_2}(x-x_0)(x-x_1)(x-x_2)f[x,x_0,x_1,x_2]dx=\alpha h^5f''''(\xi)
\end{align*}
即
\begin{align*}
\int_{x_0}^{x_0+2h}(x-x_0)(x-x_0-h)(x-x_0-2h)f[x,x_0,x_0+h,x_0+2h]dx=\alpha h^{5}f''''(\xi)
\end{align*}

然后我就不知道怎么办了.

注2：类似地可以证明$n=4$时的Newton-Cotes公式:

$n=4$时的Newton-Cotes公式：

设$f$在$(x_0,x_4)$上6阶可导，且在$[x_0,x_4]$上5阶导函数连续，则

\begin{align*}
\int_{x_0}^{x_4}f(x)dx=\frac{2h}{45}[7f(x_0)+32f(x_1)+12f(x_2)+32f(x_3)+7f(x_4)]-\frac{8h^7}{945}f^{(6)}(\xi)
\end{align*},$\xi\in (x_0,x_4)$,$h=x_1-x_0=x_2-x_1=x_3-x_2=x_4-x_3$.

我先给出一个失败的尝试:

证明:具体的验证

$$\frac{2h}{45}[7f(x_0)+32f(x_1)+12f(x_2)+32f(x_3)+7f(x_4)]$$

这种事情我是绝对不会干的，要干也得交给机器干，否则就被算死了.关键之处都在误差的估计部分,关键部分是证明

\begin{align*}
\int_{x_0}^{x_4}(x-x_0)(x-x_1)(x-x_2)(x-x_3)(x-x_4)f[x,x_0,x_1,x_2,x_3,x_4]dx=\frac{-8h^7}{945}f^{(6)}(\xi)
\end{align*}

令
\begin{align*}
g(t)=\int_{x_0}^t(x-x_0)(x-x_1)(x-x_2)(x-x_3)(x-x_4)f[x,x_0,x_1,x_2,x_3,x_4]dx
\end{align*}
则
\begin{align*}
g'(t)=(t-x_0)(t-x_1)(t-x_2)(t-x_3)(t-x_4)f[t,x_0,x_1,x_2,x_3,x_4]
\end{align*}
则
\begin{align*}
g'(x_0)=g'(x_1)=g'(x_2)=g'(x_3)=g'(x_4)=0,g(x_0)=0
\end{align*}

由于$f$在$(x_0,x_4)$上6阶可导，在$[x_0,x_4]$上5阶导函数连续，因此我们得恰当地选取好插值点，不能全选，否则导数的阶数将过大.我们选择如下的插值点进行牛顿插值
\begin{align*}
x_0,x_0' x_1,x_1',x_3,x_3',x_4,x_4'
\end{align*}
经过这几个点的牛顿插值多项式为
\begin{align*}
g(x_0)&+(x-x_0)g[x_0,x_0']\\&+(x-x_0)(x-x_0')g[x_0,x_0',x_1]\\&+(x-x_0)(x-x_0')(x-x_1)g[x_0,x_0',x_1,x_1']\\&+(x-x_0)(x-x_0')(x-x_1)(x-x_1')g[x_0,x_0',x_1,x_1',x_{3}]\\&+(x-x_0)(x-x_0')(x-x_1)(x-x_1')(x-x_3)g[x_0,x_0',x_1,x_1',x_3,x_3']\\&+(x-x_0)(x-x_0')(x-x_1)(x-x_1')(x-x_3)(x-x_3')g[x_0,x_0',x_1,x_1',x_3,x_3',x_4]\\&+(x-x_0)(x-x_0')(x-x_1)(x-x_1')(x-x_3)(x-x_3')(x-x_4)g[x_0,x_0',x_1,x_1',x_3,x_3',x_4,x_4']
\end{align*}
令
\begin{align*}
x_0'\to x_0, x_1'\to x_1,x_3'\to x_3,x_4'\to x_4
\end{align*}
可得相应的Hermite插值多项式

\begin{align*}
q(x)= g(x_0)&+(x-x_0)g[x_0,x_0]\\&+(x-x_0)^2g[x_0,x_0,x_1]\\&+(x-x_0)^2(x-x_1)g[x_0,x_0,x_1,x_1]\\&+(x-x_0)^2(x-x_1)^2g[x_0,x_0,x_1,x_1,x_3]\\&+(x-x_0)^2(x-x_1)^2(x-x_3)g[x_0,x_0,x_1,x_1,x_3,x_3]\\&+(x-x_0)^2(x-x_1)^2(x-x_3)^2g[x_0,x_0,x_1,x_1,x_3,x_3,x_4]\\&+(x-x_0)^2(x-x_1)^2(x-x_3)^2(x-x_4)g[x_0,x_0,x_1,x_1,x_3,x_3,x_4,x_4]
\end{align*}

令
\begin{align*}
k(x)=g(x)-q(x)
\end{align*}
则
\begin{align*}
k(x_0)= k(x_1)=k(x_3)=k(x_4)=0,k'(x_0)=k'(x_1)=k'(x_3)=k'(x_4)=0
\end{align*}

于是使用几次Rolle定理，可得
\begin{align*}
k^{(7)}(\xi)=0
\end{align*}
即
\begin{align*}
g^{(7)}(\xi)=q^{(7)}(\xi)
\end{align*}
即
\begin{align*}
g^{(7)}(\xi)=7!g[x_0,x_0,x_1,x_1,x_3,x_3,x_4,x_4]
\end{align*}
由于
\begin{align*}
g^{(7)}(\xi)=f^{(6)}(\xi)
\end{align*}
因此
\begin{align*}
f^{(6)}(\xi)=7!g[x_0,x_0,x_1,x_1,x_3,x_3,x_4,x_4]
\end{align*}
下面我们来计算
\begin{align*}
g[x_0,x_0,x_1,x_1,x_3,x_3,x_4,x_4]
\end{align*}

具体的计算是十分麻烦的.
\begin{align*}
& g[x_0,x_0,x_1,x_1,x_3,x_3,x_4,x_4]=\frac{g[x_0,x_0,x_1,x_1,x_3,x_3,x_4]-g[x_0,x_1,x_1,x_3,x_3,x_4,x_4]}{x_0-x_4}\\&=\frac{\frac{g[x_0,x_0,x_1,x_1,x_3,x_3]-g[x_0,x_1,x_1,x_3,x_3,x_4]}{x_0-x_4}-\frac{g[x_0,x_1,x_1,x_3,x_3,x_4]-g[x_1,x_1,x_3,x_3,x_4,x_4]}{x_0-x_4}}{x_0-x_4}\\&=\frac{g[x_0,x_0,x_1,x_1,x_3,x_3]-2g[x_0,x_1,x_1,x_3,x_3,x_4]+g[x_1,x_1,x_3,x_3,x_4,x_4]}{(x_0-x_4)^2}\\&=\frac{\frac{g[x_0,x_0,x_1,x_1,x_3]-g[x_0,x_1,x_1,x_3,x_3]}{x_0-x_3}-2
\frac{g[x_0,x_1,x_1,x_3,x_3]-g[x_1,x_1,x_3,x_3,x_4]}{x_0-x_4}+\frac{g[x_1,x_1,x_3,x_3,x_4]-g[x_1,x_3,x_3,x_4,x_4]}{x_1-x_4}}{(x_0-x_4)^2}\\&=\frac{g[x_0,x_0,x_1,x_1,x_3]-\frac{5}{2}g[x_0,x_1,x_1,x_3,x_3]+\frac{5}{2}g[x_1,x_1,x_3,x_3,x_4]-g[x_1,x_3,x_3,x_4,x_4]}{(x_0-x_4)^2(x_0-x_3)}\\&=\frac{\frac{g[x_0,x_0,x_1,x_1]-g[x_0,x_1,x_1,x_3]}{x_0-x_3}-\frac{5}{2}\frac{g[x_0,x_1,x_1,x_3]-g[x_1,x_1,x_3,x_3]}{x_0-x_3}+\frac{5}{2}\frac{g[x_1,x_1,x_3,x_3]-g[x_1,x_3,x_3,x_4]}{x_1-x_4}-\frac{g[x_1,x_3,x_3,x_4]-g[x_3,x_3,x_4,x_4]}{x_1-x_4}}{(x_0-x_4)^2(x_0-x_3)}\\&=\frac{g[x_0,x_0,x_1,x_1]-\frac{7}{2}g[x_0,x_1,x_1,x_3]+5g[x_1,x_1,x_3,x_3]-\frac{7}{2}g[x_1,x_3,x_3,x_4]+g[x_3,x_3,x_4,x_4]}{(x_0-x_4)^2(x_0-x_3)^2}\\&=\frac{\frac{g[x_0,x_0,x_1]-g[x_0,x_1,x_1]}{x_0-x_1}-\frac{7}{2}\frac{g[x_0,x_1,x_1]-g[x_1,x_1,x_3]}{x_0-x_3}+5\frac{g[x_1,x_1,x_3]-g[x_1,x_3,x_3]}{x_1-x_3}-\frac{7}{2}\frac{g[x_1,x_3,x_3]-g[x_3,x_3,x_4]}{x_1-x_4}+\frac{g[x_3,x_3,x_4]-g[x_3,x_4,x_4]}{x_3-x_4}}{(x_0-x_4)^2(x_0-x_3)^2}\\&=\frac{g[x_0,x_0,x_1]-\frac{13}{6}g[x_0,x_1,x_1]+\frac{11}{3}g[x_1,x_1,x_3]-\frac{11}{3}g[x_1,x_3,x_3]+\frac{13}{6}g[x_3,x_3,x_4]-g[x_3,x_4,x_4]}{(x_0-x_4)^2(x_0-x_3)^2(x_0-x_1)}
\end{align*}

最后我发现$g(x_1)$无法消去，所以做到这里就失败了.正确做法请见

The remainder term of the $n+1$n+1-point closed Newton-Cotes formula