背景介绍
鉴于本人实在太过疲惫和劳累，一些简单的便于理解的废话就不在记录，耗费两天时间研究了一下光流法，本帖用于记录学习过程中源码的难点细节。
使用光流预测的核心有两个利用多张连续回波图像（变分）或者相邻两张回波图（交叉相关，LK等）计算出光流场（速度矢量场），之后利用半拉格朗日方法进行外推。
半拉格朗日方法
三时间步半拉格朗日方法，其核心公式如下



但这个有点难理解，但看一维的图就很好理解



它的核心就是迭代求Am（位移矢量），如果理解困难的话，其实质上就是通过未来迭代（猜）出当前时刻的Am找到过去的假设标量值（回波值）不变的点，延伸到二维就是两个方向迭代，ok实在是so easy~
源码中lk直接调用cv的lk，没什么好说的，需要注意的是img要转成灰度图，主要是交叉相关:
    m, n = R.shape[1:]
    X, Y = np.meshgrid(np.arange(n), np.arange(m))
    print(X)
    def f(v):
        XYW = [Y + v[1], X + v[0]]
        R_w = ip.map_coordinates(
            R[-2, :, :], XYW, mode="constant", cval=np.nan, order=0, prefilter=False
        )

        mask = np.logical_and(np.isfinite(R[-1, :, :]), np.isfinite(R_w))

        return -np.corrcoef(R[-1, :, :][mask], R_w[mask])[0, 1]

    options = {"initial_simplex": (np.array([(0, 1), (1, 0), (1, 1)]))}
    result = op.minimize(f, (1, 1), method="Nelder-Mead", options=options)

    return np.stack([-result.x[0] * np.ones((m, n)), -result.x[1] * np.ones((m, n))])

首先是np.meshgrid插出坐标，其次核心公式是ip.map_coordinates，这个很难理解，查阅scipy官方文档：scipy.ndimage.map_coordinates这个函数看看,output的shape是根据坐标插值的，如果坐标不变那么output还是原来的shape，只是数值进行移动，边界进行填充：





-------未完待续---------

        龙贝格求积公式也称为逐次分半加速法。它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。 作为一种外推算法，它在不增加计算量的前提下提高了误差的精度。

        在等距基点的情况下，用计算机计算积分值通常都采用把区间逐次分半的方法进行。这样，前一次分割得到的函数值在分半以后仍可被利用,且易于编程。

一、变步长的梯形法

求 

<1> 首先在上用梯形公式，得

<2>将 二等分，步长，用复化梯形公式，得

<3>将四等分，步长  ，用复化梯形公式，得



<n>将将 2n 等分，步长  ，用复化梯形公式，得

问题：

每一步是否能用到上一步得计算结果？
	
区间划分到什么程度，才满足精度要求？
分析：

<1>、n等分，步长  ，



<2>、 2n 等分，步长  ，记 的中点(图中红线)为  。

则积分节点为：       （分为原节点 和 新增节点）



     (注意：此是计算时的)

算法总结：

计算
	
由递推公式
	
判断       （ 为误差限），成立则返回 并输出。
变步长的梯形法代码（C++)：


//变步长梯形积分
double T2n(double(*f)(double), double a, double b, double error) {
	int n = 1;   
	double h = (b - a) / n;
	double Tn = (f(a) + f(b)) * h / 2.0;  //计算 n=1 时的积分，粗略值
	double T2n;     //步长变为一半后的积分值
	bool key = false;
	do {    //至少循环一次
		double Hn = 0.0;
		for (int k = 0; k < n; k++) {
			double x = a + (k + 0.5) * h;
			Hn += h * f(x);
		}
		T2n = (Tn + Hn) / 2.0;
		//判断误差
		if (fabs(Tn - T2n) < error) {
			key = true;
		}
		else {
			Tn = T2n;
			n *= 2;
			h /= 2;
		}
 
	} while (!key);
	return T2n;
}

二、龙贝格算法

       

  准确值
	
 数值计算近似值
	
  为误差
1、      记为  .

            (序列)         

还存在误差

2、求公式的余项，可得：



       （柯特斯序列）

还存在误差

3、求柯特斯公式的余项，可得：

     （龙贝格序列）



加权平均公式：

三、代码实现


	

	

	
 
例：       ,程序结果为： 



具体代码为(没有安装armadillo矩阵计算库的点这里)：

//#include<iostream>
//#include<armadillo>
//#include<iomanip>

#include"romberg.h"

using namespace std;
using namespace arma;

int main()
{
	const double eps = 1e-6;
	int size = 50;
	Romberg.zeros(size, 4);
	Coefficient << 4 / 3.0 << 1 / 3.0 << endr
		<< 16 / 15.0 << 1 / 15.0 << endr
		<< 64 / 63.0 << 1 / 63.0 << endr;

	double a =1, b=2;
	cout << "请输入积分下限, a = ";
	cin >> a;
	cout  << "请输入积分上限, b = ";
	cin >> b;
	cout << endl;
	double result = getSn(a, b, eps);
	cout << "积分结果为：" << setprecision(8) << result << endl;
}

头文件为：

#pragma once
#include<iostream>
#include<math.h>
#include<armadillo>
#include<iomanip>
#include<vector>

using namespace std;
using namespace arma;

mat Romberg;                            //用于存放T、S、 C、R
mat Coefficient;                        //存放 计算系数
// 记录坐标值
vector<double> x_vec;
vector<double> y_vec;

// 积分函数
double fun(double x) {
	double y = log(1 + pow(x, 2)) / (1 + pow(x, 2));
	return y;
}

//存储对应坐标
void GreatXY(double a, double b, int n) {
	x_vec.clear();
	y_vec.clear();
	double h = (b - a) / n;
	for (int k = 0; k <= n; k++) {
		double x = a + h * k;
		double y = fun(x);
		x_vec.push_back(x);
		y_vec.push_back(y);
	}
}

/*变步梯形积分
   a:积分下限
   b:积分上限
   n:节点数
   */
double getTn(double a, double b, int n)
{
	GreatXY(a, b, n);
	double h = (b - a) / n;
	double Tn = h / 2.0 * (fun(a) + fun(b));
	for (int k = 1; k < n; k++) {
		double x = a + k * h;
		Tn += h * fun(x);
	}
	return Tn;
}

void FillMatrix(double &a, double &b) {
	//填充Tn
	for (int k = 0; k <= 4; k++) { //先计算k =0,1,2,3,4的情况，即将S1、S2计算出来，在|S1-S2|不符合误差情况时在更新矩阵。
		Romberg(k, 0) = getTn(a, b, pow(2, k));
	}

	//填充Sn(i=1)、Cn(i=2)、Rn(i=3)       //编程思想的行列，从零开始
	for (int i = 1; i <= 3; i++) {       //列
		for (int j = i; j <= 4; j++) {    //行
			Romberg(j, i) = Romberg(j, i - 1)*Coefficient(i - 1, 0) - Romberg(j - 1, i - 1)*Coefficient(i - 1, 1);
		}
	}
}

double getSn(double a, double b,double eps){
	FillMatrix(a, b);
	//判断误差
	if (fabs(Romberg(4, 3) - Romberg(3, 3)) < eps) {
		return Romberg(4, 3);
	}
	else 
	{
		for (int k = 5;; k++) {
			Romberg(k, 0) = getTn(a, b, pow(2, k));     //T
			Romberg(k, 1) = 4 / 3.0 *Romberg(k, 0) - 1 / 3.0*Romberg(k - 1, 0);           //S
			Romberg(k, 2) = 16 / 15.0*Romberg(k, 1) - 1 / 15.0*Romberg(k - 1, 1);      //C
			Romberg(k, 3) = 64 / 63.0 *Romberg(k, 2) - 1 / 63.0*Romberg(k - 1, 2);     //R
			//再次比较
			if (fabs(Romberg(k, 3) - Romberg(k-1, 3)) < eps) {
				return Romberg(k, 3);
				break;
			}

			if (k >= 50) {
				cout << "循环次数太多，请更换程序算法";
				return -1;
				break;
			}
		}

	}
}

 

问题背景与提出
龙贝格求积公式也称为逐次分半加速法。是数值计算方法之一，用以求解数值积分。是在梯形公式、辛普森公式和柯特斯公式之间关系的基础上，构造出一种加速计算积分的方法。 作为一种外推算法，在不增加计算量的前提下提高了误差的精度。

龙贝格公式在变步长的求积过程中运用三个加速公式，由变步长的梯形加速法，逐步得到梯形加速公式、抛物线加速公式，进而由粗糙的积分近似值迅速加工成精度较高的积分近似值。本实验通过龙贝格公式的C语言实现，两个简单示例，认识龙贝格公式的高精度积分近似。

 

实验目的
    通过对龙贝格求积公式的学习，在梯形公式、辛普森公式、柯特斯公式的基础上，实现龙贝格积分公式的C语言实现，理解龙贝格积分的的积分近似值的加工精度高的效果，进一步提高对龙贝格积分的认识。

实验原理与数学模型
    由

  可以看到，用  和  的线性组合可以得到一个比 和  都好的近似公式，因此用  作为  的近似计算方式，通过验算，可以得到即复化Simpson公式。复化Simpson公式要优于复化梯形公式。

再由复化Simpson公式的截断误差公式





若  变化不大时，即  。则

所以得到上式表明，用  和  的线性组合可以得到一个比 和  都好的近似公式。通过验算，上式右端项是  。即

当然复化Cotes公式要优于复化Simpson公式。

类似于前面的推导，可以得到关于  的线性组合公式令上式的右端项为  。即称  为龙贝格（Romberg）值。可以猜想，龙贝格（Romberg）值  要优于Cotes值  。因此用龙贝格（Romberg）值计算积分，将有较快的收敛速度。 [1] 

 

实验内容（要点）
1、龙贝格积分算法C语言实现

2、龙贝格积分算法应用

实验过程记录（含基本步骤、主要程序清单及异常情况记录等）
龙贝格积分算法C语言实现
#include <stdio.h>

#include <stdlib.h>

#include <math.h>

 

double f(double x);

double rom(double a,double b,double ep)

 

{

    double h,s1,c1,r1,s2,c2,r2,t1,t2,s;

    int i,k,n;

    h=(b-a)/2.0;

    t2=(f(a)+f(b))*h;

    s2=0;

    c2=0;

    r2=0;

    n=1;

    k=0;

    do{

        t1=t2;

        s1=s2;

        c1=c2;

        r1=r2;

        s=0.0;

 

        for(i=1;i<=n;i++)

        {

            s=s+f(a+(2*i-1)*h);

        }

        t2=t1/2+s*h;

        s2=(4*t2-t1)/3;

        c2=(16*s2-s1)/15;

        r2=(64*c2-c1)/63;

        n=n*2;

        k=k+1;

        h=h/2;

        printf("R(%d)=%f\n",k,r2);

    }

    while(fabs(r2-r1)>=ep);

 

    return 0;

 

}

龙贝格积分算法应用
1、主函数程序及函数子程序1

double f(double x)

{

    double f1;

    if(x!=0) f1=sin(x)/x;

    return f1;

}

int main()

{

    double a,b,ep;

    a=0;

    b=1;

    ep=0.0000005;

    rom(a,b,ep);

}

运行结果如下：



     2、主函数程序及函数子程序2

     double f(double y)

{

    return 1/y;

}

int main()

{

    double a,b,ep;

    a=1;

    b=3;

    ep=0.0000005;

    rom(a,b,ep);

}

运行结果如下：



实验问题及解决方法：

参考相关龙贝格实现方法，主要是一些语法错误。

实验总结
通过对龙贝格积分算法的学习，龙贝格积分算法的C语言实现及应用，对龙贝格求积公式有了一定认识，认识到其优势在于收敛速度快，积分近似值精度高。但实验中仅对该积分公式进行实现及应用，未与其他积分公式进行对比，其特性不够明显。

参考文献

同等科等，中国石油大学出版社，计算方法.
	
网络资源https://baike.baidu.com/item/%E9%BE%99%E8%B4%9D%E6%A0%BC%E6%B1%82%E7%A7%AF%E5%85%AC%E5%BC%8F/3882832?fr=aladdin
 

数值积分与数值微分
1. 绪论
2. 数值积分
2.1 引言
2.2 数值积分的基本思想
2.3 代数精度
2.4 插值型求积公式
3. 常见的几种数值积分求值公式
3.1 牛顿-柯特斯求积公式
3.2 复合求积公式
3.2.1 复合梯形求积公式
3.2.2 复合辛普森求积公式
3.2.3 复合科特斯求积公式
3.2.4 通过精度和截断误差公式确定所分区间的最小个数
3.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+th$,则可以得到

$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$等分.若将区间等分为$2n$等分,在每个小区间内仍然用梯形求积公式,则可以得到$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)}$
参考文献
[1] 数值计算

[2] 华章数学译丛24 数值计算