数值微分与数值积分(一)

最新推荐文章于 2025-04-28 15:27:03 发布

Mobtgzhang

最新推荐文章于 2025-04-28 15:27:03 发布

阅读量2.8k

点赞数 10

分类专栏：数值分析文章标签：数值积分算法

本文链接：https://blog.csdn.net/Zhang_Pro/article/details/105849743

版权

数值分析专栏收录该内容

5 篇文章

订阅专栏

数值积分与数值微分

1. 绪论
2. 数值积分
3. 常见的几种数值积分求值公式
小结
参考文献

1. 绪论

数值积分，用于求定积分的近似值。在数值分析中，数值积分是计算定积分数值的方法和理论。在数学分析中，给定函数的定积分的计算不总是可行的。许多定积分不能用已知的积分公式得到精确值。数值积分是利用黎曼积分等数学定义，用数值逼近的方法近似计算给定的定积分值。这就是计算机求解数值积分的方法。
数值微分，即根据函数在一些离散点的函数值，推算它在某点的导数或某高阶导数的近似值。通常用差商代替微商，或用一能近似代替该函数的较简单的函数（如多项式、样条函数）的相应导数作为所求导数的近似值。这就是计算机求解数值微分的方法。

2. 数值积分

2.1 引言

在工程技术和科学研究中，经常会计算定积分。很多数学问题中，例如求微分方程和积分方程的求解等也需要计算定积分。所以在区间 $[a, b]$ 内求定积分
$I(f)=\int_{a}^{b}f(x)dx$
是一个比较传统但是应用广泛的问题。从理论上讲，要计算连续函数 $f (x)$ 在 $[a, b]$ 内的定积分，只要找到 $f (x)$ 的原函数 $F (x)$ ,则由牛顿-莱布尼茨公式
$\int_{a}^{b}f(x)dx=F(b)-F(a)$
很容易求解出对应的积分函数的值，但是在实际过程中被积函数的形式比较复杂，例如下面的几种情形：

被积函数 $f (x)$ 的原函数不能使用初等函数的有限形式表示,例如 $\int_{0}^{1}\frac{\sin{x}}{x}dx$ , $\int_{0}^{1}e^{-x^{2}}dx$ 等等；
被积函数 $f (x)$ 是用函数表的形式或者图形的形式给出的，没有解析表达式；
虽然被积函数的原函数 $f (x)$ 能够使用初等函数表示，但是由于表达式过于复杂，计算定积分的值会非常困难。

2.2 数值积分的基本思想

定积分的几何意义是由曲线 $y = f (x)$ ,直线 $y = a, y = b$ 与 $x$ 轴围城的曲边梯形的面积,故而不管是 $f (x)$ 以何种形式给出的，只要近似计算出相应的曲边梯形的面积，就可以求出定积分的近似值，这就是数值求积的基本思想。
矩形法:用分点 $a=x_{0}<x_{1}<\dots<x_{n}=b$ ,将区间 $[a, b]$ 进行划分,记作
$\Delta{x}_{i}=x_{i+1}-x_{i},f(x_{i})=f_{i}$
若取 $f_{i}$ 为 $x_{i},x_{i+1}]$ 区间内矩形的高,则得到左矩形公式
$I(f)=\int_{a}^{b}f(x)dx{\approx}\Delta{x}_{0}{\cdot}f_{0}+\Delta{x}_{1}{\cdot}f_{1}+\dots+\Delta{x}_{n-1}{\cdot}f_{n-1}+0\cdot{f_{n}}$
若取 $f_{i+1}$ 为 $x_{i},x_{i+1}]$ 区间内矩形的高,则得到右矩形公式
$I(f)=\int_{a}^{b}f(x)dx{\approx}0{\cdot}f_{0}+\Delta{x}_{0}{\cdot}f_{1}+\dots+\Delta{x}_{n-2}{\cdot}f_{n-1}+\Delta{x}_{n-1}\cdot{f_{n}}$
取上述两式的平均值,得到梯形公式
$I(f)=\int_{a}^{b}f(x)dx{\approx}\dfrac{\Delta{x_{0}}}{2}f_{0}+\dfrac{\Delta{x_{0}}+\Delta{x_{1}}}{2}f_{1}+\dots+\dfrac{\Delta{x_{n-2}}+\Delta{x_{n-1}}}{2}f_{n-1}+\frac{\Delta{x_{n-1}}}{2}f_{n}$
上述三个公式的共同特点式将定积分的计算转换为计算各点函数值的线性组合,只是各个函数值的系数不同而已.一般地,在 $[a, b]$ 内适当选取 $n + 1$ 个节点 $x_{i},(i=0,1,\dots,n)$ ,然后用 $f(x_{i})$ 的现行组合作为定积分的近似值,所以,数值求积公式的一般形式为
$I(f)=\int_{a}^{b}f(x)dx\approx{\sum\limits_{i=0}^{n}A_{i}f(x_{i})}=I_{n}(f)$
上式中, $x_{i}$ 称为求积节点, $A_{i}$ 称为求积系数, $A_{i}$ 的取值仅仅与节点 $x_{i}$ 的选取有关系,二并不依赖于被积函数 $f (x)$ ,所以上式具有通用性.

2.3 代数精度

代数精度(定义)：若数值求积公式
$I(f)=\int_{a}^{b}f(x)dx\approx{\sum\limits_{i=0}^{n}A_{i}f(x_{i})}$
对被积函数 $f(x)=1,x,x^{2},\dots,x^{m}$ 都能精确成立,而对被积函数 $f(x)=x^{m+1}$ 不能精确成立,则称秋季公式具有 $m$ 次代数精度.
可见，一般地，要是使得上市具有 $n$ 次代数精度，只要令它对于 $f(x)=1,x,x^{2},\dots,x^{n}$ 都能精确成立,也就是对给定 $n + 1$ 各互异节点 $x_{i}(i=0,1,\dots,n)$ ,相应的求积系数 $A_{i}$ 满足
$\begin{cases} A_{0}+A_{1}+\dots+A_{n}=b-a\\ A_{0}x_{0}+A_{1}x_{1}+\dots+A_{n}x_{n}=\dfrac{b^{2}-a^{2}}{2}\\ \dots\\ A_{0}x_{0}^{n}+A_{1}x_{1}^{n}+\dots+A_{n}x_{n}^{n}=\dfrac{b^{n}-a^{n}}{n+1}\\ \end{cases}$
上述线性方程的系数行列式式范德蒙行列式，当 $x_{i},(i=0,1,2,\dots,n)$ 互异的时候,其值式非零的,利用克莱姆法则可以唯一地求出 $A_{i}(i=0,1,2,\dots,n)$ ,进而可以构造出求积公式.所以这样就会有这样的结论.
定理对任意给定 $n + 1$ 个互异节点 $x_{i},(i=0,1,\dots,n)$ ,总存在相应系数 $A_{i},(i=0,1,\dots,n)$ ,使得求积分公式具有 $n$ 次代数精度.

2.4 插值型求积公式

计算一种求定积分的方法可以使用插值多项式来构造数值求积公式.
设 $f (x)$ 在互异的结点 $x_{i},(i=0,1,2,\dots,n)$ 处的函数值为 $f(x_{i}),(i=0,1,2,\dots,n)$ ,则可以构造拉格朗日插值多项式
$L_{n}(x)=\sum\limits_{i=0}^{n}l_{i}(x)f(x_{i})=\sum\limits_{i=0}^{n}\left(\prod_{j=0,j\neq{i}}^{n}\dfrac{x-x_{i}}{x_{i}-x_{j}}\right)f(x_{i})$
由于代数差值多项式 $L_{n}(x)$ 的原函数是很容易求出来的,所以可以取
$\int_{a}^{b}f(x)dx\approx{\int_{a}^{b}L_{n}(x)dx}=\int_{a}^{b}\sum\limits_{i=0}^{n}l_{i}(x)f(x_{i})dx=\sum\limits_{i=0}^{n}\left(\int_{a}^{b}l_{i}(x)dx\right)f(x_{i})$
所以这样就会得到求积公式
$\int_{a}^{b}f(x)dx=\sum\limits_{i=0}^{n}A_{i}f(x_{i})$
等式中
$A_{i}=\int_{a}^{b}l_{i}(x)dx=\int_{a}^{b}\prod_{j=0,j\neq{i}}^{n}\dfrac{x-x_{i}}{x_{i}-x_{j}}dx$
所以由上式确定的求积公式称为插值型求积公式,其中余项为
$R[f]=\int_{a}^{b}\dfrac{f^{(n+1)}(\xi)}{(n+1)!}\prod_{j=0}^{n}(x-x_{j})dx$
当被积函数取的次数不超过 $n$ 次多项式时,因为 $f^{n+1}(x)=0$ ,所以余项 $R [f] = 0$ ,这说明求积公式对一切次数不超过 $n$ 次的多项式都精确成立,所以含有 $n + 1$ 个互异节点 $x_{i}(i=0,1,2,\dots,n)$ 的插值型求积公式至少具有 $n$ 次代数精度.
特别地,若求积公式至少有 $n$ 次代数精度,则它对于 $n$ 次插值基函数 $l_{i}(x)$ 是准确成立的,即有
$\int_{a}^{b}l_{i}(x)dx=\sum\limits_{j=0}^{n}A_{j}l_{j}=A_{i}$
所以就会有以下的一个基本的结论:求积公式为插值型求积公式的充分必要条件是它至少具有 $n$ 次代数精度.

3. 常见的几种数值积分求值公式

3.1 牛顿-柯特斯求积公式

求积公式的推导
设将积分区间 $[a, b]$ n等分,记步长 $h=\dfrac{b-a}{n}$ ,求积节点 $x_{i}$ 为等距节点 $x_{i}=a+ih(i=0,1,2,\dots,n)$ ,并且对应函数值 $f(x_{i})$ 是已知的,则以这些等距节点为插值节点所导出的插值型求积公式就称为是牛顿-科特斯求积公式,简称为 $N - C$ 公式.
在等距的情形下,我们改写求定积分公式
$I(f)=\int_{a}^{b}f(x)dx\approx(b-a)\int_{a}^{b}C_{i}^{(n)}f(x_{i})dx$
比较差值多项式可以得到
$C_{i}^{(n)}=\dfrac{A_{i}}{b-a}=\dfrac{1}{b-a}\prod_{j=0,j\neq{i}}\dfrac{x-x_{i}}{x_{i}-x_{j}}dx$
为了进一步简化计算,做变换 $x = a + t h$ ,则可以得到
$dx=hdt,x-x_{j}=(t-j)h,x_{i}-x_{j}=(i-j)h$
所以这样就可以得到系数
$C_{i}^{(n)}=\frac{1}{b-a}\int_{a}^{b}\dfrac{(x-x_{0})(x-x_{1})\dots(x-x_{i-1})(x-x_{i+1})\dots(x-x_{n})}{(x_{i}-x_{0})(x_{i}-x_{1})\dots(x_{i}-x_{i-1})(x_{i}-x_{i+1})\dots(x_{i}-x_{n})}dx$
$=\frac{h}{b-a}\int_{0}^{n}\dfrac{(t-0)(t-1)\dots(t-i+1)(t-i-1)\dots(t-n)}{(i-0)(i-1)\dots\cdot{1}\cdot(-1)\dots(-(n-i))}dx$
$=\frac{(-1)^{n-i}}{ni!(n-i)!}\int_{0}^{n}(t-0)(t-1)\dots(t-i+1)(t-i-1)\dots(t-n)dx$
上式系数表达的求积公式称为牛顿-科斯特公式.一般情况下采用的是系数 $n\leq{4}$ 的牛顿-科特斯公式,在 $n\geq{8}$ 的时候系数出现了负数,这样会造成不稳定的现象.
当 $n = 1$ 的时候,求积公式变为梯形求积公式
$T=\dfrac{b-a}{2}[f(a)+f(b)]$
当 $n = 2$ 的时候,求积公式变为抛物线求积公式或者是辛普森求积公式
$S=\dfrac{b-a}{6}\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right]$
当 $n = 3$ 的时候,求积公式变为 $3 / 8$ 辛普森公式
$Q=\dfrac{b-a}{8}\left[f(x_{0})+3f(x_{1})+3f(x_{2})+f(x_{3})\right]$
其中 $x_{i}=a+ih,h=\dfrac{b-a}{3}$
当 $n = 4$ 的时候,求积公式变为科斯特求积公式
$C=\dfrac{b-a}{90}\left[7f(x_{0})+32f(x_{1})+12f(x_{2})+32f(x_{3})+7f(x_{4})\right]$
其中 $x_{i}=a+ih,h=\dfrac{b-a}{4}$
代数精度
我们可以得到以下的定理:对于牛顿-科斯特求积公式,当 $n$ 为奇数的时候,它至少具有 $n$ 次代数精度;当 $n$ 为偶数的时候,它具有 $n + 1$ 次代数精度.
证明:
当 $n$ 为偶数的时候,设被积函数 $f(x)=x^{n+1}$ ,由于 $f^{(n+1)}(x)=(n+1)!$ ,所以此时的余项即为
$R[f]=\int_{a}^{b}\dfrac{f^{(n+1)}(x)}{(n+1)!}\prod_{j=0}^{n}(x-x_{j})dx=\int_{a}^{b}\prod_{j=0}^{n}(x-x_{j})dx$
进行变换 $x=a+th,x_{j}=a+jh$ ,所以就会得到
$R[f]=h^{n+2}\int_{0}^{n}\prod_{j=0}^{n}(t-j)dt$
设 $n=2k,(k\in{N^{+}})$ 进行变换 $t = s + k$ ,则可以得到
$R[f]=h^{2k+2}\int_{-k}^{k}\prod_{j=0}^{2k}(s+k-j)ds$
设关于 $s$ 的函数
$G(s)=\prod_{j=0}^{j=0}(s+k-j)=(s+k)(s+k-1)\dots{s(s-1)}\dots{(s-k+1)(s-k)}$
$=\prod_{j=-k}^{k}(s-j)$
所以这就会得到
$G(-s)=\prod_{j=-k}^{k}(-s-j)=(-1)^{2k+1}\prod_{j=-k}^{k}(s+j)=-\prod_{j=-k}^{k}(s-j)=-G(s)$
所以 $G (s)$ 是奇函数,所以它在对称区间 $[- k, k]$ 的定积分为 $0$ ,故而就会得到 $R [f] = 0$ .即当 $n$ 为偶数的时候,对被积函数 $f(x)=x^{n+1}$ 的余项为 $0$ ,故而得到对 $f(x)=x^{n+1}$ 精确成立.
当 $n$ 为奇数的时候,不妨设 $f(x)=x^{n}$ ,很显然余项 $f^{n+1}(x)=0$ ,余项 $R [f] = 0$ ,所以当 $n$ 为奇数的时候,对被积函数 $x^{n}$ 的余项为 $0$ ,故而得到得到对 $f(x)=x^{n}$ 精确成立.
误差估计
可以得到牛顿-科特斯公式的误差为
$R[f]=\int_{a}^{b}\dfrac{f^{n+1}(\xi)}{(n+1)!}\prod_{j=0}^{n}(x-x_{j})dx$
其中 $\xi=a+(b-a)\theta(x)$ ,由积分中值定理可以得到 $\exists{\eta\in{(a,b)}}$ ,使得
$R[f]=\dfrac{f^{n+1}(\eta)}{(n+1)!}\int_{a}^{b}\prod_{j=0}^{n}(x-x_{j})dx$
而
$x_{j}=a+jh,x=a+th,h=\frac{b-a}{n},dx=hdt$
所以就会得到
$R[f]=\dfrac{h^{n+2}f^{n+1}(\eta)}{(n+1)!}\int_{0}^{n}\prod_{j=0}^{n}(t-j)dx$
令
$a_{n}=\int_{0}^{n}\prod_{j=0}^{n}(t-j)dx$
特别地,当 $n = 1$ 的时候, $a_{1}=-\dfrac{1}{6}$ 就会得到梯形公式的截断误差为
$R_{1}(f)=-\dfrac{f^{''}(\eta)}{12}(b-a)^{3}$
当 $n = 2$ 的时候, $a_{2}=-\dfrac{1}{60\cdot{2^{4}}}$ 就会得到梯形公式的截断误差为
$R_{2}(f)=-\dfrac{f^{''}(\eta)}{180\cdot{2^{4}}}(b-a)^{5}$
收敛性和稳定性
首先我们这样定义数值的稳定性
定义:对于任意 $\epsilon>0$ ,若存在 $\delta>0$ ,只要 $\delta_{i}=\left|f(x_{i})-\bar{f}_{i}\right|,(i=0,1,2,\dots,n)$ ,就有
$\left|I_{n}(f)-I_{n}(\bar{f})\right|\leq{\epsilon}$
则称背记公式是数值稳定的.
对于牛顿-科特斯公式来说有以下的定理:
若牛顿-科特斯公式中的系数 $C_{i}^{(n)}>0(i=0,1,2,\dots,n)$ ,则求积公式是稳定的,证明如下
由于
$C_{i}^{(n)}>0,(i=0,1,2,\dots,n),\left|f(x_{i})-\bar{f}\right|\leq{\delta}$
所以就会得到
$\left|I_{n}(f)-I_{n}(\bar{f})\right|=\left|(b-a)\sum\limits_{i=0}^{n}\left[C_{i}^{(n)}\left(f(x_{i})-\bar{f}\right)\right]\right|\leq{\delta{(b-a)}}\sum\limits_{i=0}^{n}C_{i}^{(n)}$
故而对于任意的 $\epsilon>0$ ,存在 $\delta=\dfrac{\epsilon}{b-a}$ ,只要 $\left|f(x_{i})-\bar{f}\right|\leq{\delta}$ ,就会得到
$\left|I_{n}(f)-I_{n}(\bar{f})\right|\leq{\delta{(b-a)}}\leq{\epsilon}$ ,
所以求积公式是数值稳定的.
在牛顿-科特斯公式中, $n\geq{8}$ 的时候出现负值的系数,此时求积公式的收敛性没有保证所以在实际中很少使用到高阶的牛顿-科特斯公式,一般使用的是 $n=1\sim{4}$ 的牛顿-科特斯公式.

3.2 复合求积公式

从梯形公式、辛普森公式以及科特斯公式求积公式的余项式可以看出,数值求积的误差除了与被积函数有关系,还与积分区间的长度 $(b - a)$ 有关系,积分区间越小,求积公式的截断误差就会越小,故而在求定积分的时候,通常把积分区间等分成若干小区间,在每个小区间内采用次数不高的求积公式,然后再将其加起来得到整个区间内的求积公式,这就是符合求积公式的基本思想.

3.2.1 复合梯形求积公式

将 $[a, b]$ 区间 $n$ 等分,记分点为
$x_{i}=a+ih,h=\dfrac{b-a}{n},(i=0,1,2,\dots,n)$
并在每个小区间 $x_{i},x_{i+1}]$ 内应用梯形公式,得到
$\int_{a}^{b}f(x)dx=\sum\limits_{i=0}^{n-1}\int_{x_{i}}^{x_{i+1}}f(x)dx\approx\sum\limits_{i=0}^{n-1}\dfrac{h}{2}\left[f(x_{i})+f(x_{i+1})\right]$
$=\dfrac{h}{2}\left[f(a)+2\sum\limits_{i=1}^{n-1}f(x_{i})+f(b)\right]$
令
$T_{n}=\dfrac{h}{2}\left[f(a)+2\sum\limits_{i=1}^{n-1}f(x_{i})+f(b)\right]$
上式称为复合梯形公式,下标 $n$ 表示将区间 $n$ 等分.若将区间等分为 $2 n$ 等分,在每个小区间内仍然用梯形求积公式,则可以得到 $T_{2n}$ . $T_{n}$ 与 $T_{2n}$ 之间的关系为
$T_{2n}=\dfrac{1}{2}(T_{n}+H_{n})$
等式中 $H_{n}=h\sum\limits_{i=1}^{n}f(x_{i+\frac{1}{2}})$ , $x_{i+\frac{1}{2}}$ 为区间 $[x_{i},x_{i+\frac{1}{2}}]$ 的中点,即 $x_{i+\frac{1}{2}}=x_{i}+\frac{1}{2}h$ ,根据定积分的定义可以得到
$\lim\limits_{n\rightarrow\infty,h\rightarrow{0}}T_{n}=\lim\limits_{h\rightarrow{0}}\dfrac{1}{2}\left[\sum\limits_{i=0}^{n-1}f(x_{i})h+\sum\limits_{i=1}^{n}f(x_{i})h\right]=\int_{a}^{b}f(x)dx=I(f)$
可见复合型公式是收敛的,并且等式中的系数 $A_{i}>0(i=0,1,2,\dots,n)$ ,故而也是稳定的.截断误差为
$R_{n}(f)=I(f)-T_{n}=\sum\limits_{i=0}^{n-1}\left[-\dfrac{h^{3}}{12}f^{''}(\eta_{i})\right]=-\dfrac{h^{2}}{12}(b-a)\dfrac{1}{n}\sum\limits_{i=0}^{n-1}f^{''}(\eta_{i}),(\eta_{i}\in(x_{i},x_{i+1}))$
由连续函数的性质以及介值定理可以得到, $\exists{\eta\in{(a,b)}}$ ,使得
$f^{''}(\eta)=\dfrac{1}{n}\sum\limits_{i=0}^{n-1}f^{''}(\eta_{i})$
所以得到
$R(f)=-\dfrac{b-a}{12}h^{2}f^{''}(\eta),(\eta\in(a,b))$
即复合型梯形公式的截断误差, $R_{n}(f)=O(h^{2})$ .
用程序写出来可以得到

import math
def integral(func,a,b,delta = 0.1):
    n = int((b-a)/delta)
    h = (b-a)/n
    v_sum = h*(func(a)+func(b))/2
    for k in range(1,n-1):
        xi = a + k*delta
        v_sum = v_sum + h*func(xi)
    return v_sum
def main():
    def sinx(x):
        if x == 0:
            return 1
        return math.sin(x)/x
    print(integral(sinx,0,1))
if __name__ == "__main__":
    main()

运行的结果为(精确值为 $\int_{0}^{1}\frac{\sin{x}}{x}=0.9460831\dots$ )

0.9458320718669053

第二种通过递推公式的方法可以得到

import math
def integral(func,a,b,delta = 0.1):
    n = int((b-a)/delta)
    h = (b-a)/n
    vT = (func(a)+func(b))*h/2
    vH = 0
    for k in range(1,n):
        xi = a + k*h
        xip = xi + h/2
        vH = vH + h*func(xip)
        vT = vT + h*func(xi)
    xip = xi + h/2
    vH = vH + h*func(xip)
    return (vH+vT)/2
def main():
    def sinx(x):
        if x == 0:
            return 1
        return math.sin(x)/x
    print(integral(sinx,0,1))
if __name__ == "__main__":
    main()

得到的结果为

0.9388524984416825

3.2.2 复合辛普森求积公式

在每个小区间 $x_{i},x_{i+1}]$ ,使用辛普森公式可以得到
$\int_{a}^{b}f(x)dx\approx\sum\limits_{i=0}^{n-1}\dfrac{h}{6}\left[f(x_{i})+4f(x_{i+\frac{1}{2}})+f(x_{i+1})\right]$
$=\dfrac{h}{6}\left[f(a)+4\sum\limits_{i=0}^{n-1}f(x_{i+\frac{1}{2}})+2\sum\limits_{i=1}^{n-1}f(x_{i})+f(b)\right]$
上式中 $x_{i+\frac{1}{2}}$ 为区间 $x_{i},x_{i+1}]$ 的中点,即 $x_{i+\frac{1}{2}}=x_{i}+\dfrac{h}{2}$
记
$S_{n}=\dfrac{h}{6}\left[f(a)+4\sum\limits_{i=0}^{n-1}f(x_{i+\frac{1}{2}})+2\sum\limits_{i=1}^{n-1}f(x_{i})+f(b)\right]$
上式称为辛普森公式,其中余项可以得到
$R_{n}(f)=I(f)-S_{n}=-\dfrac{h}{180}\left(\dfrac{h}{2}\right)^{4}\sum\limits_{i=0}^{n-1}f^{(4)}(\eta_{i}),(\eta_{i}\in{(x_{i},x_{i+1})})$
由介值定理同样可以得到 $\exists\eta\in{(a,b)}$ ,使得
$R_{n}(f)=-\dfrac{b-a}{180}\left(\dfrac{h}{2}\right)^{4}f^{(4)}(\eta)$
上式即为辛普森公式的截断误差,即得到 $R_{n}(f)=O(h^{4})$ .
同样可以证明
$\lim\limits_{n\rightarrow\infty}S_{n}=\int_{a}^{b}f(x)dx$
由于各个系数 $A_{i}>0,i=0,1,2,\dots,n$ ,所以这个结果是稳定的.
为了编程方便可以将等式改写为
$S_{n}=\dfrac{h}{6}\left\{f(a)-f(b)+\sum\limits_{i=1}^{n}\left[4f(x_{i+\frac{1}{2}})+2f(x_{i})\right]\right\}$
程序求解上述表达式可以得到

import math
def integral(func,a,b,delta = 0.1):
    n = int((b-a)/delta)
    h = (b-a)/n
    v_sum = h*(func(a)-func(b))/6
    for k in range(1,n+1):
        xi = a + k*h
        xip = xi+h/2
        v_sum = v_sum + h*(4*func(xip)+2*func(xi))/6
    return v_sum
def main():
    def sinx(x):
        if x == 0:
            return 1
        return math.sin(x)/x
    print(integral(sinx,0,1))
if __name__ == "__main__":
    main()

通过程序求解可以得到

0.9345186746707335

3.2.3 复合科特斯求积公式

同样,在每一个小区间内 $x_{i},x_{i+1}]$ ,使用科斯特公式可以得到
$\int_{a}^{b}f(x)dx\approx{\sum\limits_{i=0}^{n-1}\dfrac{h}{90}\left[7f(x_{i})+32f(x_{i+\frac{1}{4}})+12f(x_{i+\frac{1}{2}})+32f(x_{i+\frac{3}{4}})+7f(x_{i+1})\right]}$
$=\dfrac{h}{90}\left[7f(a)+14\sum\limits_{i=1}^{n-1}f(x_{i})+32\sum\limits_{i=0}^{n-1}f(x_{i+\frac{1}{4}})+\right.$
$12\sum\limits_{i=0}^{n-1}f(x_{i+\frac{1}{2}})+32\sum\limits_{i=0}^{n-1}f(x_{i+\frac{3}{4}})+7f(b)\left.\right]$
令
$C_{n}=\dfrac{h}{90}\left[7f(a)+14\sum\limits_{i=1}^{n-1}f(x_{i})+32\sum\limits_{i=0}^{n-1}f(x_{i+\frac{1}{4}})+\right.$
$12\sum\limits_{i=0}^{n-1}f(x_{i+\frac{1}{2}})+32\sum\limits_{i=0}^{n-1}f(x_{i+\frac{3}{4}})+7f(b)\left.\right]$
显然 $C_{n}$ 是收敛的,并且截断误差
$R_{n}(f)=-\dfrac{2(b-a)}{945}\left(\dfrac{b-a}{4}\right)^{6}f^{6}(\eta)$ 其中 $\eta\in{(a,b)}$ ,所以 $R_{n}(f)=O(h^{6})$ ,可以得到系数 $A_{i}>0$ ,所以此公式是稳定的.
为了编程方便,我们可以将公式改写为
$C_{n}=\dfrac{h}{90}\left\{7f(a)-7f(b)+\sum\limits_{i=1}^{n}\left[14f(x_{i})+32f(x_{i+\frac{1}{4}})+12f(x_{i+\frac{1}{2}})+32f(x_{i+\frac{3}{4}})\right]\right\}$
其中上式中 $x_{i+\frac{k}{4}}=x_{i}+\frac{k}{4}h,(k=0,1,2,3)$
用程序求解上述复合积分求解公式可以得到

import math
def integral(func,a,b,delta = 0.1):
    n = int((b-a)/delta)
    h = (b-a)/n
    v_sum = 7*h*(func(a)-func(b))/90
    for k in range(1,n+1):
        xi = a + k*h
        xia = xi+h/4
        xib = xi+2*h/4
        xic = xi+3*h/4
        v_sum = v_sum + h*(14*func(xi)+32*func(xia)+12*func(xib)+32*func(xic))/90
    return v_sum
def main():
    def sinx(x):
        if x == 0:
            return 1
        return math.sin(x)/x
    print(integral(sinx,0,1))
if __name__ == "__main__":
    main()

得到的结果为

0.9314371161293896

3.2.4 通过精度和截断误差公式确定所分区间的最小个数

根据截断误差公式和精度可以列出不等式
$\left|R_{n}(f)\right|=\left|(b-a)\cdot{V}\cdot{h^{m}}\right|\leq{\epsilon}$
其中 $V$ 是对应的参数, $m$ 是关于 $n$ 的正整数,求解出 $n$ 这样就可以得到所分区间的最小 $n$ 的值.

3.3 开牛顿-科特斯方法

上述的求解方法一般称作是闭牛顿-科特斯方法.它的特点是需要积分区间端点上的函数值.一些被积分函数在端点上有可以移除的奇点,使用开牛顿-科特斯方法可以更加容易处理一些，其中不使用在区间端点上的函数值.
令 $h=x-x_{0}$ ,这样 $f (x)$ 在区间中点 $\omega=(x_{0}+\frac{h}{2})$ 处的泰勒展开式表示为
$f(x)=f(\omega)+(x-\omega)f^{'}(\omega)+\frac{1}{2}(x-\omega)^{2}f^{''}(c_{x})$
所以在 $x_{0},x_{1}]$ 之间的积分可以表示为
$\int_{x_{0}}^{x_{1}}f(x)dx=(x_{1}-x_{0})f(\omega)+f^{'}(\omega)\int_{x_{0}}^{x_{1}}(x-\omega)dx+\frac{1}{2}\int_{x_{0}}^{x_{1}}f^{''}(c_{x})(x-\omega)^{2}dx$
$=hf(\omega)+\frac{h^{3}}{24}f^{''}(c_{x})$
这种方法称之为中心法则,即得到积分表达式
$\int_{x_{0}}^{x_{1}}f(x)dx=hf(\omega)+\frac{h^{3}}{24}f^{''}(c_{x})$
其中 $h=x_{1}-x_{0},\omega=\frac{a+b}{2}$
复合中心法则表达式为
$\int_{a}^{b}f(x)dx=h\sum\limits_{i=1}^{m}f(\omega_{i})+\frac{(b-a)h^{2}}{24}f^{''}(c_{x})$
其中 $h=\frac{b-a}{n},w_{i}=x_{0}+\frac{h}{2}$

小结

最为重要的两个基本的法则(由牛顿-科特斯公式推导出来的结果):
梯形法则( $h=x_{1}-x_{0},c\in{(x_{0},x_{1})}$ )
$\int_{x_{0}}^{x_{1}}f(x)dx=\dfrac{h}{2}(y_{0}+y_{1})-\dfrac{h^{3}}{12}f^{''}(c)$
辛普森法则( $h=x_{1}-x_{0}=x_{2}-x_{1},c\in{(x_{0},x_{2})}$ )
$\int_{x_{0}}^{x_{1}}f(x)dx=\dfrac{h}{3}\left(y_{0}+4y_{1}+y_{2}\right)-\frac{h^{5}}{90}f^{(4)}(c)$
复合梯形公式
$\int_{a}^{b}f(x)dx=\dfrac{h}{2}(y_{0}+y_{m}+2\sum\limits_{k=1}^{m-1}y_{i})-\dfrac{(b-a)h^{2}}{12}f^{''}(c)$
其中 $h=\frac{b-a}{m},c\in{(a,b)}$
复合辛普森公式
$\int_{a}^{b}f(x)dx=\dfrac{h}{3}\left(y_{0}+y_{m}+4\sum\limits_{i=0}^{n-1}y_{i+\frac{1}{2}}+2\sum\limits_{i=1}^{n-1}y_{i}\right)-\frac{(b-a)h^{4}}{180}f^{(4)}(c)$
其中 $h=\frac{b-a}{m},c\in{(a,b)}$