精华内容
下载资源
问答
  • 使用光流预测的核心有两个利用多张连续回波图像(变分)或者相邻两张回波图(交叉相关,LK等)计算出光流场(速度矢量场),之后利用半拉格朗日方法进行外推。 半拉格朗日方法 三时间步半拉格朗日方法,其核心公式...

    背景介绍

    鉴于本人实在太过疲惫和劳累,一些简单的便于理解的废话就不在记录,耗费两天时间研究了一下光流法,本帖用于记录学习过程中源码的难点细节。

    使用光流预测的核心有两个利用多张连续回波图像(变分)或者相邻两张回波图(交叉相关,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,只是数值进行移动,边界进行填充:
    在这里插入图片描述
    在这里插入图片描述
    -------未完待续---------

    展开全文
  • 并以山西省霍州煤电集团公司某煤矿采区水文地质补充勘查报告矿井涌水量预测为例进行了计算,计算结果表明:非稳定流解析法的雅柯布(近似)公式对于煤层底板突水的矿井涌水量预测,较稳定流解析法和传统的类比外推法在...
  • 首先,我们将计算π和σ介子筛选质量的公式扩展到有限化学势μ的情况。 然后,我们考虑虚构-μ方法,这是一种从虚构化学势(μ=iμI)到实数(μ=μR)的外推方法。 基于纠缠Polyakov环扩展的Nambu–Jona-Lasinio...
  • 作为一种外推算法, 它在不增加计算量的前提下提高了误差的精度. 在等距基点的情况下,用计算机计算积分值通常都采用把区间逐次分半的方法进行。这样,前一次分割得到的函数值在分半以后仍可被利用,且易于编程。
  • 8.2.1 混合波场外推 8.2.2 相对误差分析 8.2.3 成像计算与并行实现 8.3 混合三维平面波合成叠前深度偏移 8.3.1 三维平面波合成与目标照明 8.3.2 因子分解波场外推 8.3.3 成像计算 8.4 共方位数据三维叠前偏移 ...
  • 作为一种外推算法,它在不增加计算量的前提下提高了误差的精度。 在等距基点的情况下,用计算机计算积分值通常都采用把区间逐次分半的方法进行。这样,前一次分割得到的函数值在分半以后仍可被利用,且易于编程。 ...

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

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

    一、变步长的梯形法

    求 \int_{a}^{b}f(x)dx

    <1> 首先在[a,b]上用梯形公式,得T_{1}

    <2>将[a,b] 二等分,步长h=\frac{b-a}{2},用复化梯形公式,得T_{2}

    <3>将[a,b]四等分,步长h=\frac{b-a}{4}  ,用复化梯形公式,得T_{4}

    \cdots \cdots

    <n>将将[a,b] 2n 等分,步长h=\frac{b-a}{2n}  ,用复化梯形公式,得T_{2n}

    问题:

    1. 每一步是否能用到上一步得计算结果?
    2. 区间划分到什么程度,才满足精度要求?

    分析:

    <1>、n等分,步长 h=\frac{b-a}{n} ,T_{n}=\frac{h}{2}[f(a)+f(b)+2\sum_{k=1}^{n-1}f(x_{k})]

    <2>、 2n 等分,步长 h^{'}=\frac{h}{2} ,记[x_{k},x_{k+1}] 的中点(图中红线)为 x_{k+\frac{1}{2}} 。

    则积分节点为:x_{0},x_{\frac{1}{2}},x_{1},x_{\frac{3}{2}},x_{2},\cdots ,x_{n-1},x_{n-\frac{1}{2}},x_{n}       (分为原节点 和 新增节点)

    T_{2n} = \frac{h^{'}}{2}[f(a)+f(b)+2\sum_{k=1}^{n-1}f(x_{k}) + 2\sum_{k=0}^{n-1}f(x_{k+\frac{1}{2}})]

    =\frac{1}{2}T_{n}+\frac{h}{2}\cdot \sum_{k=0}^{n-1}f(x_{k+\frac{1}{2}})     (注意:此h是计算T_{n}时的)

    算法总结:

    1. 计算T_{n}=\frac{h}{2}[f(a)+f(b)+2\sum_{k=1}^{n-1}f(x_{k})]
    2. 由递推公式T_{2n}=\frac{1}{2}T_{n}+\frac{h}{2}\cdot \sum_{k=0}^{n-1}f(x_{k+\frac{1}{2}})
    3. 判断  |T_{n} - T_{2n}| < \varepsilon     (\varepsilon 为误差限),成立则返回 T_{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;
    }

    二、龙贝格算法

    I - T_{2n} \approx \frac{1}{3}(T_{2n}-T_{n})       

    • I  准确值
    • T_{2n} 数值计算近似值
    • \frac{1}{3}(T_{2n}-T_{n})  为误差

    1、 I \approx T_{2n}+ \frac{1}{3}(T_{2n}-T_{n})     记为 S_{n} .

    =S_{n}            (simpson序列)         

    还存在误差

    2、simpson公式的余项,可得:

    \frac{I-S_{2n}}{I-S_{n}}\approx \frac{1}{16}\Rightarrow I-S_{2n}\approx \frac{1}{15}(S_{2n}-S_{n})

    \therefore C_{n} = S_{2n}+\frac{1}{15}(S_{2n}-S_{n})       (柯特斯序列)

    还存在误差

    3、求柯特斯公式的余项,可得:

    R_{n}=C_{2n}+\frac{1}{63}\cdot (C_{2n}-C_{n})     (龙贝格序列)

    加权平均公式:\frac{4T_{n}-T_{n}}{4-1} = S_{n}\cdots \frac{4^{2}S_{2n}-S_{n}}{4^{2}-1}=C_{n}\cdots \frac{4^{3}C_{2n}-C_{n}}{4^{3}-1}=R_{n}

    三、代码实现

    • T_{n}=\frac{h}{2}[f(a)+f(b)+2\sum_{k=1}^{n-1}f(x_{k})]
    • S_{n}=T_{2n}+ \frac{1}{3}(T_{2n}-T_{n}) = \frac{4}{3}T_{2n}-\frac{1}{3}T_{n}
    • C_{n} = S_{2n}+\frac{1}{15}(S_{2n}-S_{n}) = \frac{16}{15}S_{2n}-\frac{1}{15}S_{n}
    • R_{n}=C_{2n}+\frac{1}{63}\cdot (C_{2n}-C_{n}) = \frac{64}{63}C_{2n}-\frac{1}{63}C_{n} 

    例: f(x) = \int_{1}^{2}\frac{log(1+x^{2})}{1+x^{2}}      ,程序结果为: 

    具体代码为(没有安装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语言实现及应用

    千次阅读 2019-06-15 16:29:43
    作为一种外推算法,在不增加计算量的前提下提高了误差的精度。 龙贝格公式在变步长的求积过程中运用三个加速公式,由变步长的梯形加速,逐步得到梯形加速公式、抛物线加速公式,进而由粗糙的积分近似值迅速加工...
    • 问题背景与提出

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

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

     

    • 实验目的

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

    • 实验原理与数学模型

       

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

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

    若  变化不大时,即  。则

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

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

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

     

    • 实验内容(要点)

    1、龙贝格积分算法C语言实现

    2、龙贝格积分算法应用

    • 实验过程记录(含基本步骤、主要程序清单及异常情况记录等)
    1. 龙贝格积分算法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、主函数程序及函数子程序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语言实现及应用,对龙贝格求积公式有了一定认识,认识到其优势在于收敛速度快,积分近似值精度高。但实验中仅对该积分公式进行实现及应用,未与其他积分公式进行对比,其特性不够明显。

    参考文献

    1. 同等科等,中国石油大学出版社,计算方法.
    2. 网络资源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

     

    展开全文
  • Python 科学计算

    2018-09-20 16:59:31
    3.3.2 外推和 Spline 拟合.....................90 3.3.3 二维插值.....................................91 3.4 数值积分——integrate .................. 93 3.4.1 球的体积.......................................
  • 数值积分3.1 牛顿-柯特斯求积公式3.2 复合求积公式3.3 外推法与龙贝格公式3.4 高斯求积公式小结参考文献 1.引言 2.数值微分 3.数值积分 3.1 牛顿-柯特斯求积公式 3.2 复合求积公式 3.3 外推法与龙贝格公式 3.4 高斯...

    1. 绪论

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

    2. 数值积分

    2.1 引言

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

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

    2.2 数值积分的基本思想

    定积分的几何意义是由曲线y=f(x)y=f(x),直线y=a,y=by=a,y=bxx轴围城的曲边梯形的面积,故而不管是f(x)f(x)以何种形式给出的,只要近似计算出相应的曲边梯形的面积,就可以求出定积分的近似值,这就是数值求积的基本思想。
    矩形法:用分点a=x0<x1<<xn=ba=x_{0}<x_{1}<\dots<x_{n}=b,将区间[a,b][a,b]进行划分,记作
    Δxi=xi+1xi,f(xi)=fi\Delta{x}_{i}=x_{i+1}-x_{i},f(x_{i})=f_{i}
    若取fif_{i}[xi,xi+1][x_{i},x_{i+1}]区间内矩形的高,则得到左矩形公式
    I(f)=abf(x)dxΔx0f0+Δx1f1++Δxn1fn1+0fnI(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}}
    若取fi+1f_{i+1}[xi,xi+1][x_{i},x_{i+1}]区间内矩形的高,则得到右矩形公式
    I(f)=abf(x)dx0f0+Δx0f1++Δxn2fn1+Δxn1fnI(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)=abf(x)dxΔx02f0+Δx0+Δx12f1++Δxn2+Δxn12fn1+Δxn12fnI(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][a,b]内适当选取n+1n+1个节点xi,(i=0,1,,n)x_{i},(i=0,1,\dots,n),然后用f(xi)f(x_{i})的现行组合作为定积分的近似值,所以,数值求积公式的一般形式为
    I(f)=abf(x)dxi=0nAif(xi)=In(f)I(f)=\int_{a}^{b}f(x)dx\approx{\sum\limits_{i=0}^{n}A_{i}f(x_{i})}=I_{n}(f)
    上式中,xix_{i}称为求积节点,AiA_{i}称为求积系数,AiA_{i}的取值仅仅与节点xix_{i}的选取有关系,二并不依赖于被积函数f(x)f(x),所以上式具有通用性.

    2.3 代数精度

    代数精度(定义):若数值求积公式
    I(f)=abf(x)dxi=0nAif(xi)I(f)=\int_{a}^{b}f(x)dx\approx{\sum\limits_{i=0}^{n}A_{i}f(x_{i})}
    对被积函数f(x)=1,x,x2,,xmf(x)=1,x,x^{2},\dots,x^{m}都能精确成立,而对被积函数f(x)=xm+1f(x)=x^{m+1}不能精确成立,则称秋季公式具有mm次代数精度.
    可见,一般地,要是使得上市具有nn次代数精度,只要令它对于f(x)=1,x,x2,,xnf(x)=1,x,x^{2},\dots,x^{n}都能精确成立,也就是对给定n+1n+1各互异节点xi(i=0,1,,n)x_{i}(i=0,1,\dots,n),相应的求积系数AiA_{i}满足
    {A0+A1++An=baA0x0+A1x1++Anxn=b2a22A0x0n+A1x1n++Anxnn=bnann+1\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}
    上述线性方程的系数行列式式范德蒙行列式,当xi,(i=0,1,2,,n)x_{i},(i=0,1,2,\dots,n)互异的时候,其值式非零的,利用克莱姆法则可以唯一地求出Ai(i=0,1,2,,n)A_{i}(i=0,1,2,\dots,n),进而可以构造出求积公式.所以这样就会有这样的结论.
    定理 对任意给定n+1n+1个互异节点xi,(i=0,1,,n)x_{i},(i=0,1,\dots,n),总存在相应系数Ai,(i=0,1,,n)A_{i},(i=0,1,\dots,n),使得求积分公式具有nn次代数精度.

    2.4 插值型求积公式

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

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

    3.1 牛顿-柯特斯求积公式

    求积公式的推导
    设将积分区间[a,b][a,b]n等分,记步长h=banh=\dfrac{b-a}{n},求积节点xix_{i}为等距节点xi=a+ih(i=0,1,2,,n)x_{i}=a+ih(i=0,1,2,\dots,n),并且对应函数值f(xi)f(x_{i})是已知的,则以这些等距节点为插值节点所导出的插值型求积公式就称为是牛顿-科特斯求积公式,简称为NCN-C公式.
    在等距的情形下,我们改写求定积分公式
    I(f)=abf(x)dx(ba)abCi(n)f(xi)dxI(f)=\int_{a}^{b}f(x)dx\approx(b-a)\int_{a}^{b}C_{i}^{(n)}f(x_{i})dx
    比较差值多项式可以得到
    Ci(n)=Aiba=1baj=0,jixxixixjdxC_{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+thx=a+th,则可以得到
    dx=hdt,xxj=(tj)h,xixj=(ij)hdx=hdt,x-x_{j}=(t-j)h,x_{i}-x_{j}=(i-j)h
    所以这样就可以得到系数
    Ci(n)=1baab(xx0)(xx1)(xxi1)(xxi+1)(xxn)(xix0)(xix1)(xixi1)(xixi+1)(xixn)dxC_{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
    =hba0n(t0)(t1)(ti+1)(ti1)(tn)(i0)(i1)1(1)((ni))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
    =(1)nini!(ni)!0n(t0)(t1)(ti+1)(ti1)(tn)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
    上式系数表达的求积公式称为牛顿-科斯特公式.一般情况下采用的是系数n4n\leq{4}的牛顿-科特斯公式,在n8n\geq{8}的时候系数出现了负数,这样会造成不稳定的现象.
    n=1n=1的时候,求积公式变为梯形求积公式
    T=ba2[f(a)+f(b)]T=\dfrac{b-a}{2}[f(a)+f(b)]
    n=2n=2的时候,求积公式变为抛物线求积公式或者是辛普森求积公式
    S=ba6[f(a)+4f(a+b2)+f(b)]S=\dfrac{b-a}{6}\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right]
    n=3n=3的时候,求积公式变为3/83/8辛普森公式
    Q=ba8[f(x0)+3f(x1)+3f(x2)+f(x3)]Q=\dfrac{b-a}{8}\left[f(x_{0})+3f(x_{1})+3f(x_{2})+f(x_{3})\right]
    其中xi=a+ih,h=ba3x_{i}=a+ih,h=\dfrac{b-a}{3}
    n=4n=4的时候,求积公式变为科斯特求积公式
    C=ba90[7f(x0)+32f(x1)+12f(x2)+32f(x3)+7f(x4)]C=\dfrac{b-a}{90}\left[7f(x_{0})+32f(x_{1})+12f(x_{2})+32f(x_{3})+7f(x_{4})\right]
    其中xi=a+ih,h=ba4x_{i}=a+ih,h=\dfrac{b-a}{4}
    代数精度
    我们可以得到以下的定理:对于牛顿-科斯特求积公式,当nn为奇数的时候,它至少具有nn次代数精度;当nn为偶数的时候,它具有n+1n+1次代数精度.
    证明:
    nn为偶数的时候,设被积函数f(x)=xn+1f(x)=x^{n+1},由于f(n+1)(x)=(n+1)!f^{(n+1)}(x)=(n+1)!,所以此时的余项即为
    R[f]=abf(n+1)(x)(n+1)!j=0n(xxj)dx=abj=0n(xxj)dxR[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,xj=a+jhx=a+th,x_{j}=a+jh,所以就会得到
    R[f]=hn+20nj=0n(tj)dtR[f]=h^{n+2}\int_{0}^{n}\prod_{j=0}^{n}(t-j)dt
    n=2k,(kN+)n=2k,(k\in{N^{+}})进行变换t=s+kt=s+k,则可以得到
    R[f]=h2k+2kkj=02k(s+kj)dsR[f]=h^{2k+2}\int_{-k}^{k}\prod_{j=0}^{2k}(s+k-j)ds
    设关于ss的函数
    G(s)=j=0j=0(s+kj)=(s+k)(s+k1)s(s1)(sk+1)(sk)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)}
    =j=kk(sj)=\prod_{j=-k}^{k}(s-j)
    所以这就会得到
    G(s)=j=kk(sj)=(1)2k+1j=kk(s+j)=j=kk(sj)=G(s)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)G(s)是奇函数,所以它在对称区间[k,k][-k,k]的定积分为00,故而就会得到R[f]=0R[f]=0.即当nn为偶数的时候,对被积函数f(x)=xn+1f(x)=x^{n+1}的余项为00,故而得到对f(x)=xn+1f(x)=x^{n+1}精确成立.
    nn为奇数的时候,不妨设f(x)=xnf(x)=x^{n},很显然余项fn+1(x)=0f^{n+1}(x)=0,余项R[f]=0R[f]=0,所以当nn为奇数的时候,对被积函数xnx^{n}的余项为00,故而得到得到对f(x)=xnf(x)=x^{n}精确成立.
    误差估计
    可以得到牛顿-科特斯公式的误差为
    R[f]=abfn+1(ξ)(n+1)!j=0n(xxj)dxR[f]=\int_{a}^{b}\dfrac{f^{n+1}(\xi)}{(n+1)!}\prod_{j=0}^{n}(x-x_{j})dx
    其中ξ=a+(ba)θ(x)\xi=a+(b-a)\theta(x),由积分中值定理可以得到η(a,b)\exists{\eta\in{(a,b)}},使得
    R[f]=fn+1(η)(n+1)!abj=0n(xxj)dxR[f]=\dfrac{f^{n+1}(\eta)}{(n+1)!}\int_{a}^{b}\prod_{j=0}^{n}(x-x_{j})dx

    xj=a+jh,x=a+th,h=ban,dx=hdtx_{j}=a+jh,x=a+th,h=\frac{b-a}{n},dx=hdt
    所以就会得到
    R[f]=hn+2fn+1(η)(n+1)!0nj=0n(tj)dxR[f]=\dfrac{h^{n+2}f^{n+1}(\eta)}{(n+1)!}\int_{0}^{n}\prod_{j=0}^{n}(t-j)dx

    an=0nj=0n(tj)dxa_{n}=\int_{0}^{n}\prod_{j=0}^{n}(t-j)dx
    特别地,当n=1n=1的时候,a1=16a_{1}=-\dfrac{1}{6}就会得到梯形公式的截断误差为
    R1(f)=f(η)12(ba)3R_{1}(f)=-\dfrac{f^{''}(\eta)}{12}(b-a)^{3}
    n=2n=2的时候,a2=16024a_{2}=-\dfrac{1}{60\cdot{2^{4}}}就会得到梯形公式的截断误差为
    R2(f)=f(η)18024(ba)5R_{2}(f)=-\dfrac{f^{''}(\eta)}{180\cdot{2^{4}}}(b-a)^{5}
    收敛性和稳定性
    首先我们这样定义数值的稳定性
    定义:对于任意ϵ>0\epsilon>0,若存在δ>0\delta>0,只要δi=f(xi)fˉi,(i=0,1,2,,n)\delta_{i}=\left|f(x_{i})-\bar{f}_{i}\right|,(i=0,1,2,\dots,n),就有
    In(f)In(fˉ)ϵ\left|I_{n}(f)-I_{n}(\bar{f})\right|\leq{\epsilon}
    则称背记公式是数值稳定的.
    对于牛顿-科特斯公式来说有以下的定理:
    若牛顿-科特斯公式中的系数Ci(n)>0(i=0,1,2,,n)C_{i}^{(n)}>0(i=0,1,2,\dots,n),则求积公式是稳定的,证明如下
    由于
    Ci(n)>0,(i=0,1,2,,n),f(xi)fˉδC_{i}^{(n)}>0,(i=0,1,2,\dots,n),\left|f(x_{i})-\bar{f}\right|\leq{\delta}
    所以就会得到
    In(f)In(fˉ)=(ba)i=0n[Ci(n)(f(xi)fˉ)]δ(ba)i=0nCi(n)\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)}
    故而对于任意的ϵ>0\epsilon>0,存在δ=ϵba\delta=\dfrac{\epsilon}{b-a},只要f(xi)fˉδ\left|f(x_{i})-\bar{f}\right|\leq{\delta},就会得到
    In(f)In(fˉ)δ(ba)ϵ\left|I_{n}(f)-I_{n}(\bar{f})\right|\leq{\delta{(b-a)}}\leq{\epsilon},
    所以求积公式是数值稳定的.
    在牛顿-科特斯公式中,n8n\geq{8}的时候出现负值的系数,此时求积公式的收敛性没有保证所以在实际中很少使用到高阶的牛顿-科特斯公式,一般使用的是n=14n=1\sim{4}的牛顿-科特斯公式.

    3.2 复合求积公式

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

    3.2.1 复合梯形求积公式

    [a,b][a,b]区间nn等分,记分点为
    xi=a+ih,h=ban,(i=0,1,2,,n)x_{i}=a+ih,h=\dfrac{b-a}{n},(i=0,1,2,\dots,n)
    并在每个小区间[xi,xi+1][x_{i},x_{i+1}]内应用梯形公式,得到
    abf(x)dx=i=0n1xixi+1f(x)dxi=0n1h2[f(xi)+f(xi+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]
    =h2[f(a)+2i=1n1f(xi)+f(b)]=\dfrac{h}{2}\left[f(a)+2\sum\limits_{i=1}^{n-1}f(x_{i})+f(b)\right]

    Tn=h2[f(a)+2i=1n1f(xi)+f(b)]T_{n}=\dfrac{h}{2}\left[f(a)+2\sum\limits_{i=1}^{n-1}f(x_{i})+f(b)\right]
    上式称为复合梯形公式,下标nn表示将区间nn等分.若将区间等分为2n2n等分,在每个小区间内仍然用梯形求积公式,则可以得到T2nT_{2n}.TnT_{n}T2nT_{2n}之间的关系为
    T2n=12(Tn+Hn)T_{2n}=\dfrac{1}{2}(T_{n}+H_{n})
    等式中Hn=hi=1nf(xi+12)H_{n}=h\sum\limits_{i=1}^{n}f(x_{i+\frac{1}{2}}),xi+12x_{i+\frac{1}{2}}为区间[xi,xi+12][x_{i},x_{i+\frac{1}{2}}]的中点,即xi+12=xi+12hx_{i+\frac{1}{2}}=x_{i}+\frac{1}{2}h,根据定积分的定义可以得到
    limn,h0Tn=limh012[i=0n1f(xi)h+i=1nf(xi)h]=abf(x)dx=I(f)\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)
    可见复合型公式是收敛的,并且等式中的系数Ai>0(i=0,1,2,,n)A_{i}>0(i=0,1,2,\dots,n),故而也是稳定的.截断误差为
    Rn(f)=I(f)Tn=i=0n1[h312f(ηi)]=h212(ba)1ni=0n1f(ηi),(ηi(xi,xi+1))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}))
    由连续函数的性质以及介值定理可以得到,η(a,b)\exists{\eta\in{(a,b)}},使得
    f(η)=1ni=0n1f(ηi)f^{''}(\eta)=\dfrac{1}{n}\sum\limits_{i=0}^{n-1}f^{''}(\eta_{i})
    所以得到
    R(f)=ba12h2f(η),(η(a,b))R(f)=-\dfrac{b-a}{12}h^{2}f^{''}(\eta),(\eta\in(a,b))
    即复合型梯形公式的截断误差,Rn(f)=O(h2)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()
    

    运行的结果为(精确值为01sinxx=0.9460831\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 复合辛普森求积公式

    在每个小区间[xi,xi+1][x_{i},x_{i+1}],使用辛普森公式可以得到
    abf(x)dxi=0n1h6[f(xi)+4f(xi+12)+f(xi+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]
    =h6[f(a)+4i=0n1f(xi+12)+2i=1n1f(xi)+f(b)]=\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]
    上式中xi+12x_{i+\frac{1}{2}}为区间[xi,xi+1][x_{i},x_{i+1}]的中点,即xi+12=xi+h2x_{i+\frac{1}{2}}=x_{i}+\dfrac{h}{2}

    Sn=h6[f(a)+4i=0n1f(xi+12)+2i=1n1f(xi)+f(b)]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]
    上式称为辛普森公式,其中余项可以得到
    Rn(f)=I(f)Sn=h180(h2)4i=0n1f(4)(ηi),(ηi(xi,xi+1))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})})
    由介值定理同样可以得到η(a,b)\exists\eta\in{(a,b)},使得
    Rn(f)=ba180(h2)4f(4)(η)R_{n}(f)=-\dfrac{b-a}{180}\left(\dfrac{h}{2}\right)^{4}f^{(4)}(\eta)
    上式即为辛普森公式的截断误差,即得到Rn(f)=O(h4)R_{n}(f)=O(h^{4}).
    同样可以证明
    limnSn=abf(x)dx\lim\limits_{n\rightarrow\infty}S_{n}=\int_{a}^{b}f(x)dx
    由于各个系数Ai>0,i=0,1,2,,nA_{i}>0,i=0,1,2,\dots,n,所以这个结果是稳定的.
    为了编程方便可以将等式改写为
    Sn=h6{f(a)f(b)+i=1n[4f(xi+12)+2f(xi)]}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 复合科特斯求积公式

    同样,在每一个小区间内[xi,xi+1][x_{i},x_{i+1}],使用科斯特公式可以得到
    abf(x)dxi=0n1h90[7f(xi)+32f(xi+14)+12f(xi+12)+32f(xi+34)+7f(xi+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]}
    =h90[7f(a)+14i=1n1f(xi)+32i=0n1f(xi+14)+=\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.
    12i=0n1f(xi+12)+32i=0n1f(xi+34)+7f(b)]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]

    Cn=h90[7f(a)+14i=1n1f(xi)+32i=0n1f(xi+14)+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.
    12i=0n1f(xi+12)+32i=0n1f(xi+34)+7f(b)]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]
    显然CnC_{n}是收敛的,并且截断误差
    Rn(f)=2(ba)945(ba4)6f6(η)R_{n}(f)=-\dfrac{2(b-a)}{945}\left(\dfrac{b-a}{4}\right)^{6}f^{6}(\eta)其中η(a,b)\eta\in{(a,b)},所以Rn(f)=O(h6)R_{n}(f)=O(h^{6}),可以得到系数Ai>0A_{i}>0,所以此公式是稳定的.
    为了编程方便,我们可以将公式改写为
    Cn=h90{7f(a)7f(b)+i=1n[14f(xi)+32f(xi+14)+12f(xi+12)+32f(xi+34)]}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\}
    其中上式中xi+k4=xi+k4h,(k=0,1,2,3)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 通过精度和截断误差公式确定所分区间的最小个数

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

    3.3 开牛顿-科特斯方法

    上述的求解方法一般称作是闭牛顿-科特斯方法.它的特点是需要积分区间端点上的函数值.一些被积分函数在端点上有可以移除的奇点,使用开牛顿-科特斯方法可以更加容易处理一些,其中不使用在区间端点上的函数值.
    h=xx0h=x-x_{0},这样f(x)f(x)在区间中点ω=(x0+h2)\omega=(x_{0}+\frac{h}{2})处的泰勒展开式表示为
    f(x)=f(ω)+(xω)f(ω)+12(xω)2f(cx)f(x)=f(\omega)+(x-\omega)f^{'}(\omega)+\frac{1}{2}(x-\omega)^{2}f^{''}(c_{x})
    所以在[x0,x1][x_{0},x_{1}]之间的积分可以表示为
    x0x1f(x)dx=(x1x0)f(ω)+f(ω)x0x1(xω)dx+12x0x1f(cx)(xω)2dx\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(ω)+h324f(cx)=hf(\omega)+\frac{h^{3}}{24}f^{''}(c_{x})
    这种方法称之为中心法则,即得到积分表达式
    x0x1f(x)dx=hf(ω)+h324f(cx)\int_{x_{0}}^{x_{1}}f(x)dx=hf(\omega)+\frac{h^{3}}{24}f^{''}(c_{x})
    其中h=x1x0,ω=a+b2h=x_{1}-x_{0},\omega=\frac{a+b}{2}
    复合中心法则表达式为
    abf(x)dx=hi=1mf(ωi)+(ba)h224f(cx)\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=ban,wi=x0+h2h=\frac{b-a}{n},w_{i}=x_{0}+\frac{h}{2}

    小结

    最为重要的两个基本的法则(由牛顿-科特斯公式推导出来的结果):
    梯形法则(h=x1x0,c(x0,x1)h=x_{1}-x_{0},c\in{(x_{0},x_{1})})
    x0x1f(x)dx=h2(y0+y1)h312f(c)\int_{x_{0}}^{x_{1}}f(x)dx=\dfrac{h}{2}(y_{0}+y_{1})-\dfrac{h^{3}}{12}f^{''}(c)
    辛普森法则(h=x1x0=x2x1,c(x0,x2)h=x_{1}-x_{0}=x_{2}-x_{1},c\in{(x_{0},x_{2})})
    x0x1f(x)dx=h3(y0+4y1+y2)h590f(4)(c)\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)
    复合梯形公式
    abf(x)dx=h2(y0+ym+2k=1m1yi)(ba)h212f(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=bam,c(a,b)h=\frac{b-a}{m},c\in{(a,b)}
    复合辛普森公式
    abf(x)dx=h3(y0+ym+4i=0n1yi+12+2i=1n1yi)(ba)h4180f(4)(c)\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=bam,c(a,b)h=\frac{b-a}{m},c\in{(a,b)}

    参考文献

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

    展开全文
  • 机器学习笔记--概率与数理统计

    千次阅读 2017-09-02 16:13:13
    1 概率记作P(E),比如掷骰子,每一面的概率就是P(E) = 1/61.1 古典概率通常又叫事前概率,是指当随机事件中各种可能发生的结果及其出现的次数都可以由演绎或外推法得知,而无需经过任何统计试验即可计算各种可能发生...
  • 龙贝格积分

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

    2020-12-21 19:08:45
    作为一种外推算法,在不增加计算量的前提下提高了误差的精度。 求积步骤 算法设计 设计思想为 梯形公式经过 区间逐步分半的方法的梯形公式求积 就等于辛普森公式求积 辛普森公式经过 区间逐步分半的方法的梯形公式...
  • C语言实现龙贝格求积

    2016-01-20 11:26:31
    龙贝格求积公式也称为逐次分半加速。它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上,构造出一种加速计算积分的方法。 作为一种外推算法, 它在不增加计算量的前提下提高了误差的精度.
  • 外推法建立了变断面非满水钻井井壁竖向稳定临界深度计算公式,并分析了其影响因素;借助ABAQUS软件分析了不同阶段钻井井壁的竖向稳定性问题,计算了钻井井壁竖向稳定临界深度,分析了井壁自重、不同边界条件、井壁底部...
  • 龙贝格求积算法

    千次阅读 2018-07-15 09:23:44
     作为一种外推算法,它在不增加计算量的前提下提高了误差的精度。在等距基点的情况下,用计算机计算积分值通常都采用把区间逐次分半的方法进行。这样,前一次分割得到的函数值在分半以后仍可被利用,且易于编程。...
  • Romberg求积分算法

    千次阅读 2014-12-27 22:37:27
    // Integral-romberg方法求积分.cpp : 定义控制台应用程序的入口点。...作为一种外推算法,它在不增加计算量的前提下提高了误差的精度。 在等距基点的情况下,用计算机计算积分值通常都采用吧区间逐次分半的方
  • DEWT 用外推法求一阶常微分方程的数值解 DEWT_glg 用格拉格外推法求一阶常微分方程的数值解 第16章: 偏微分方程的数值解法 peEllip5 用五点差分格式解拉普拉斯方程 peEllip5m 用工字型差分格式解拉普拉斯方程 ...
  • MATLAB常用算法

    热门讨论 2010-04-05 10:34:28
    DEWT 用外推法求一阶常微分方程的数值解 DEWT_glg 用格拉格外推法求一阶常微分方程的数值解 第16章: 偏微分方程的数值解法 peEllip5 用五点差分格式解拉普拉斯方程 peEllip5m 用工字型差分格式解拉普拉斯方程 ...
  • 《统计预测:方法与应用》比较详尽地介绍了用于预测的定量分析... 第五章 趋势外推法  第六章 季节变动预测法  第七章 马尔可夫法  第八章 博克斯——詹金斯法  第九章 ECM模型和ARCH模型的应用  参考书目
  • 三角函数的有理逼近习题六第七章 数值积分与数值微分§7.1 插值型求积公式与代数度§7.2 复合求积公式及算法§7.3 外推原理与龙贝格算法§7.4 高斯型求积公式及其复合公式§7.5 数值微分应用:通信卫星覆盖地球...
  • +数值微分:一阶导数的数值计算及其MATLAB程序(差商求导,中心差商公式,理查森外推法求导等等,不再赘述)[ +解线性方程组的直接方法:方程组的逆矩阵解法及其MATLAB程序,三角形方程组的解法及其MATLAB程序,高斯...
  • 求根过程:是叠代过程,即由(x1,x2)→f(x1)、f(x2)、f(x中)或f(x外推)→(X1,X2),大写X1,X2就是下一轮计算的小写x1,x2,二分法、截弦法、牛顿迭代法计算公式不同,一个用中值外推,后二者用直线外推,二者用直线外推...
  • 基于人工冻结稳态温度场假设,依托Бахолдин单排管模型建立相应稳态温度场的周期形式计算模型,引入狄拉克函数,利用格林函数获得其解析解的周期表达形式,并外推获得广义单排管冻结稳态温度场解析解的矩阵...
  • 适应最近由C. Marboe和D.... 对于L = 1和自旋S = 1的最简单算子,在比较大的耦合值范围内,12回路结果的Padé外推法与可用的热力学Bethe Ansatz数据非常吻合。 附带了带有结果选择的Mathematica笔记本。
  • 统计学方法与数据分析(上下册)

    热门讨论 2013-12-29 11:32:47
    16.3外推问题 16.4多维协变量和更复杂的设计 16.5小结 补充练习 第十七章一些固定效应、随机效应和混合效应模型的方差分析 17.1引言和案例 17.2具有随机处理效应的单因子试验:随机效应模型 17.3随机效应模型...
  • 数学基础研究基本数学概念(数、几何形状、集合、函数...),及其 如何构成更复杂结构和概念的 层次结构,特别是一类 也被称为元数学概念的 基础性重要结构,用它们来形成数学语言(公式、理论、以及它们的 用来表意...

空空如也

空空如也

1 2
收藏数 33
精华内容 13
关键字:

外推法计算公式