复化梯形求积分实例——用Python进行数值计算
作者:行行无别语只道早还乡 发布时间:2021-10-08 13:03:26
用程序来求积分的方法有很多,这篇文章主要是有关牛顿-科特斯公式。
学过插值算法的同学最容易想到的就是用插值函数代替被积分函数来求积分,但实际上在大部分场景下这是行不通的。
插值函数一般是一个不超过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个小区间。
来源:https://www.cnblogs.com/zhangte/p/6156212.html
猜你喜欢
- 记录应用程序的操作日志可以使用数据库、文本文件、XML文件等。我这里介绍的是使用 XML 文件记录操作日志。我觉得使用 XML 记录操作日志
- 问题:在安装SP4补丁的时候,老是报验证密码错误。上网查了一下资料,发现是一个小bug。按照一下操作,安装正常。SQL Server补丁安装
- 在SQL Server数据库管理中,针对分析服务Analysis Services 的性能优化必不可少,这里我们将学习到使用DMV来进行An
- 用于制作自动化微信聊天图片,通过图片生成段子视频根据一个txt文档input.txtL 一路走过来好热啊
- 创意404页面的文章我们似乎已经出过两篇了,今天hongkiat又带来了60个创意404页面.相关404页面设计文章:国外404错误页面的创
- 前言上一篇介绍了服务端流式RPC,客户端发送请求到服务器,拿到一个流去读取返回的消息序列。 客户端读取返回的流的数据。本篇将介绍客户端流式R
- 很多时候我们都需要了解下python中导入包的属性方法信息,当然dir 是最便捷的了,不过如果想知道特定的,例如以_ 开头的属性,需要写个筛
- 此文刊登在《程序员》三月期,有删改提到安全问题,首先想到应付这些问题的应该是系统管理员以及后台开发工程师们,而前端开发工程师似乎离这些问题很
- 关于JavaSctipt的兼容性,最懒的办法就是用jQuery的工具函数。尽量不要用那些什么ECMAScript之类的函数,因为很多浏览器都
- 废话不多说了,直接给大家贴js代码了,具体代码如下所示:<!DOCTYPE html><html><head&
- 本文实例讲述了Yii2框架实现登陆添加验证码功能。分享给大家供大家参考,具体如下:models中LoginForm.phppublic $v
- Python中的模块(.py文件)在创建之初会自动加载一些内建变量,__name__就是其中之一。Python模块中通常会定义很多变量和函数
- ubuntu 系统自带的 python 有多个版本,使用时难免会遇到环境变量出错,特别是当自动化运行脚本的时候。特别是近一个月来,实验室的小
- table.rows集合中是cell对象 cell.innerHTML = "<td>123</td>&q
- 1、我的源码在 /home/topsec/Documents/php-7.0.11 ,安装位置在 /usr/local/php7, php.
- 一、 文件的操作1.1创建文件格式:f = open(‘文件', ‘w')或者f = open(‘文件', ‘r
- 数据库并行访问,也就是两个或两以上用户同时访问同一数据,这也是数据库引擎如何设计和实现适度反应所面临的最大问题。设计优良、性能卓越的数据库引
- 目的对字符串的每个字符进行处理,其实每个字符(Char)就是一个长度为1的字符串。方法1.使用内建函数list()>>>
- 一、目的之前在博文SQL Server数据库最小宕机迁移方案中提到了使用了完全备份+差异备份的功能完成了数据库的转移,但是这个方法在遇到了7
- Macromedia Dreamweaver MX 2004提供了更多功能强劲的可视化设计工具、应用开