用程序来求积分的方法有很多，这篇文章主要是有关牛顿-科特斯公式。

学过插值算法的同学最容易想到的就是用插值函数代替被积分函数来求积分，但实际上在大部分场景下这是行不通的。

插值函数一般是一个不超过n次的多项式，如果用插值函数来求积分的话，就会引进高次多项式求积分的问题。这样会将原来的求积分问题带到另一个求积分问题：如何求n次多项式的积分，而且当次数变高时，会出现龙悲歌现象，误差反而可能会增大，并且高次的插值求积公式有可能会变得不稳定：详细原因不赘述。

牛顿-科特斯公式解决这一问题的办法是将大的插值区间分为一堆小的插值区间，使得多项式的次数不会太高。然后通过引入参数函数

将带有幂的项的取值范围固定在一个固定范围内，这样一来就将多项式带有幂的部分的求积变为一个固定的常数，只需手工算出来即可。这个常数可以直接带入多项式求积函数。

上式中x的求积分区间为[a, b]，h = (b - a)/n, 这样一来积分区间变为[0， n]，需要注意的是从这个公式可以看出一个大的区间被分为n个等长的小区间。 这一部分具体请参见任意一本有关数值计算的书！

n是一个事先确定好的值。

又因为一个大的插值区间需要被分为等长的多个小区间，并在这些小区间上分别进行插值和积分，因此此时的牛顿-科特斯公式被称为：复化牛顿-科特斯公式。

并且对于n的不同取值牛顿-科特斯有不同的名称： 当n=1时，叫做复化梯形公式，复化梯形公式也就是将每一个小区间都看为一个梯形(高为h，上底为f(t), 下底为f(t+1))。这与积分的本质：无限分隔  相同。

当n=2时，复化牛顿-科特斯公式被称为复化辛普森公式(非美国法律界著名的那个辛普森)。

我这篇文章实现的是复化梯形公式：

首先写一个函数求节点函数值求和那部分：

"""

@brief: 求和 ∑f(xk) : xk表示等距节点的第k个节点，不包括端点

xk = a + kh (k = 0, 1, 2, ...)

积分区间为[a, b]

@param: xk 积分区间的等分点x坐标集合(不包括端点)

@param: func 求积函数

@return: 返回值为集合的和

"""

def sum_fun_xk(xk, func):

return sum([func(each) for each in xk])

然后就可以写整个求积分函数了：

"""

@brief: 求func积分 :

@param: a 积分区间左端点

@param: b 积分区间右端点

@param: n 积分分为n等份(复化梯形求积分要求)

@param: func 求积函数

@return: 积分值

"""

def integral(a, b, n, func):

h = (b - a)/float(n)

xk = [a + i*h for i in range(1, n)]

return h/2 * (func(a) + 2 * sum_fun_xk(xk, func) + func(b))

相当的简单

试验：

当把大区间分为两个小区间时：

分为20个小区间时：

求的积分值就是这些彩色的梯形面积之和。

测试代码：

if __name__ == "__main__":

func = lambda x: x**2

a, b = 2, 8

n = 20

print integral(a, b, n, func)

''' 画图 '''

import matplotlib.pyplot as plt

plt.figure("play")

ax1 = plt.subplot(111)

plt.sca(ax1)

tmpx = [2 + float(8-2) /50 * each for each in range(50+1)]

plt.plot(tmpx, [func(each) for each in tmpx], linestyle = '-', color='black')

for rang in range(n):

tmpx = [a + float(8-2)/n * rang, a + float(8-2)/n * rang, a + float(8-2)/n * (rang+1), a + float(8-2)/n * (rang+1)]

tmpy = [0, func(tmpx[1]), func(tmpx[2]), 0]

c = ['r', 'y', 'b', 'g']

plt.fill(tmpx, tmpy, color=c[rang%4])

plt.grid(True)

plt.show()

注意上面代码中的n并不是上文开篇提到的公式中的n，开篇提到的n是指将每一个具体的插值区间(也就是小区间)等距插n个节点，复化梯形公式的n是固定的为1.

而代码中的n指将大区间分为n个小区间。

python与数值计算环境搭建

数值计算的编程的软件很多种,也见过一些编程绘图软件的对比. 利用Python进行数值计算,需要用到numpy(矩阵) ,scipy(公式符号), matplotlib(绘图)这些工具包. 1.Linu ...

【编程练习】收集的一些c++代码片，算法排序，读文件，写日志，快速求积分等等

写日志: class LogFile { public: static LogFile &instance(); operator FILE *() const { return m_file ...

之前写的收集的一些c++代码片，算法排序，读文件，写日志，快速求积分等等

写日志: class LogFile { public: static LogFile &instance(); operator FILE *() const { return m_file ...

牛顿插值法——用Python进行数值计算

拉格朗日插值法的最大毛病就是每次引入一个新的插值节点,基函数都要发生变化,这在一些实际生产环境中是不合适的,有时候会不断的有新的测量数据加入插值节点集, 因此,通过寻找n个插值节点构造的的插值函数与n ...

Python 3 数值计算

Python 3.4.3 (v3.4.3:9b73f1c3e601, Feb 24 2015, 22:43:06) [MSC v.1600 32 bit (Intel)] on win32Type & ...

分段二次插值——用Python进行数值计算

事实上在实际使用中,高次插值显然是很不适合的,高次插值将所有样点包涵进一个插值函数中,这是次幂高的原因.高次计算复杂,而且刚开始的一点误差会被方的很大.因此将整个区间分为若干个小区间,在每一个小区间进 ...

埃尔米特插值问题——用Python进行数值计算

当插值的要求涉及到对插值函数导数的要求时,普通插值问题就变为埃尔米特插值问题.拉格朗日插值和牛顿插值的要求较低,只需要插值函数的函数值在插值点与被插函数的值相等,以此来使得在其它非插值节点插值函数的值 ...

拉格朗日插值法——用Python进行数值计算

插值法的伟大作用我就不说了.... 那么贴代码? 首先说一下下面几点: 1. 已有的数据样本被称之为 "插值节点" 2. 对于特定插值节点,它所对应的插值函数是必定存在且唯一的(关 ...

simpson法求积分 专题练习

[xsy1775]数值积分 题意 多组询问,求\(\int_l^r\sqrt{a(1-{x^2\over b})}dx\) 分析 double f(double x) { return sqrt(a* ...

随机推荐

Reverse Core 第一部分 代码逆向技术基础

@date: 2016/10/14 笔记 记录书中较重要的知识,方便回顾 ps. 因有一些逆向基础,所以我本来就比较熟悉的知识并未详细记录 第一章 关于逆向工程 目标, ...

HDU 1052

http://acm.hdu.edu.cn/showproblem.php?pid=1052 田忌赛马本质就是一个贪心 res表示田忌的胜利场次 1.田忌最快马快于王的最快马,两个最快马比,res++ ...

SQL Sever2008 新手入门第一天安装软件

(计应154兰家才)这学期,新来了一门课程,数据库应用.刚开始什么都不懂,也不知道这东西到底是干嘛,本着路漫漫其修远兮,吾将上下而求索的精神,开始了一段求知路程.刚开始找了一个简单的绿化版sql200 ...

v2013调试无法访问此网站 localhost 拒绝了我们的连接请求

问题描述:          别人给的服务器代码,在本地部署以后调试的,localhost:8080 可以访问,localhost:2524访问不了需要改什么配置吗 解决思路:           这 ...

k8s实战之从私有仓库拉取镜像 - kubernetes

1.实战目的 从私有docker仓库拉取镜像,部署pod.上一篇中,我们搭建了私有的镜像仓库,这一篇我们将与k8s结合实战使用私有仓库. 2.登录docker 为了完成本次实战,需要登录docker, ...

HBase——HMaster启动之一(HMaster的构建)

首先,让我们来到HMaster的main方法.我们今天的流程就从这里开始. 我们需要注意,下图所示的tool的类型就是HMasterCommandLine. 接下来,让我们来到HMasterComma ...

2019 校内赛 RPG的地牢猎手(bfs+优先队列)

Problem Description Luke最近沉迷一款RPG游戏,游戏中角色可以进入地牢关卡,只要顺利走出地牢就可以获得奖励.地牢表示为n行m列的块矩阵,其中每个块只可以是障碍块.入口.出口或数 ...

zabbix3.x添加H3C网络设备详解

zabbix3.x添加H3C网络设备详解 作者:尹正杰 版权声明:原创作品,谢绝转载!否则将追究法律责任. 前言: 欢迎加入:高级运维工程师之路 598432640 相信大家在看我的文章之前,也看过其 ...

Mysql 登录及用户切换、用户权限查询

启动mysql: 方法一:net start mysql(或者是其他服务名) 方法二:在windows下启动MySQL服务 MySql安装目录:"d:\MySql\" 进入命令行输 ...

Activity设置Dialog属性点击区域外消失实现方式

通过配置:true 通过代码:setFinish ...