复化梯形求积分实例——用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


猜你喜欢
- 自带验证器1、UniqueValidator用于验证(唯一)unique=True的字段,常用参数:queryset: required,用
- strSourceFile = Server.MapPath(dataXML&"/Advertisement/"
- 在开始课程之前,我要求学生们填写一份调查表,这个调查表反映了它们对Python中一些概念的理解情况。一些话题("if/
- 本文实例讲述了Python使用win32 COM实现Excel的写入与保存功能。分享给大家供大家参考,具体如下:很久之前通过东拼西凑实现过使
- 1.环境准备1.MySQL 安装路径:/usr/local2.CentOS 6.2 服务器(64 bit)3.MySQL 5.7.28 下载
- 某大师曾说过,像了解自己的老婆 一样了解自己管理的数据库,个人认为包含了两个方面的了解:1,在稳定性层面来说,更多的是关注高可用、读写分离、
- int(10)int(20)分别代表什么意思储备知识在设计数据库表的时候,经常需要设计一个id字段,它的类型一般都是整型int,经常会遇到i
- 一个简单的验证码爬取程序本文介绍了在Python2.7环境下爬取网站验证码:思路就是获取验证码对应的url,然后发起requst请求,读取该
- 有一个表tb_3a_huandan_detail,每天有300W左右的数据。查询太慢了,网上了解了一下,可以做表分区。由于数据较大,所以决定
- 音乐流媒体服务的兴起使得音乐无处不在。我们在上下班的时候听音乐,锻炼身体,工作或者只是放松一下。这些服务的一个关键特性是播放列表,通常按流派
- 本文实例讲述了Python wxPython库消息对话框MessageDialog用法。分享给大家供大家参考,具体如下:消息对话框即我们平时
- 内容有些啰嗦,内容记载了当时遇到了bug以及解决问题的思路。业务场景简述:前端做配置化组件,通过url内的唯一标识,请求后端sql 哪取页面
- creatdoc.asp<!DOCTYPE HTML PUBLIC "-//W3C/DTD&n
- 目录前言Tips - django版本区别路由匹配无名分组&有名分组无名分组有名分组小提示反向解析路由不涉及分组的反向解析有名分组&
- 准备工作首先,我们需要确保以下几点:你已经安装了MySQL数据库,并且可以正常连接。你已经配置好了MyBatis的环境,并且可以成功执行单条
- 分组取TOP数据是T-SQL中的常用查询, 如学生信息管理系统中取出每个学科前3名的学生。这种查询在SQL Server 2005之前,写起
- 1. zip() 函数的介绍1.1 功能zip() 函数用于将可迭代的对象作为参数,将对象中对应的元素打包成一个个元组,然后返回由这些元组组
- Jones向量假设光波沿z轴传播,那么其三个方向的电场分量可以表示为Jones矩阵能够保证二维列向量形状不变的运算有无穷多种,但最符合我们直
- 在开发过程中,我们经常会将日期时间的毫秒数存放到数据库,但是它对应的时间看起来就十分不方便,我们可以使用一些函数将毫秒转换成date格式。
- 何为质数: 只能被1 和 自身 整除的数;方法: 利用js中求模, 看是否有余数. ---> 3%2 = 1; 5%2 = 3....