复化梯形求积分实例——用Python进行数值计算

所属分类: 脚本专栏 / python 阅读数: 1024
收藏 0 赞 0 分享

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

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

插值函数一般是一个不超过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实现按学生年龄排序的实际问题详解

这篇文章主要给大家介绍了关于Python实现按学生年龄排序实际问题的相关资料,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友们下面跟着小编来一起学习学习吧。
收藏 0 赞 0 分享

Python开发的HTTP库requests详解

Requests是用Python语言编写,基于urllib,采用Apache2 Licensed开源协议的HTTP库。它比urllib更加方便,可以节约我们大量的工作,完全满足HTTP测试需求。Requests的哲学是以PEP 20 的习语为中心开发的,所以它比urllib更加P
收藏 0 赞 0 分享

Python网络爬虫与信息提取(实例讲解)

下面小编就为大家带来一篇Python网络爬虫与信息提取(实例讲解)。小编觉得挺不错的,现在就分享给大家,也给大家做个参考。一起跟随小编过来看看吧
收藏 0 赞 0 分享

在python3环境下的Django中使用MySQL数据库的实例

下面小编就为大家带来一篇在python3环境下的Django中使用MySQL数据库的实例。小编觉得挺不错的,现在就分享给大家,也给大家做个参考。一起跟随小编过来看看吧
收藏 0 赞 0 分享

Python 3.x读写csv文件中数字的方法示例

在我们日常开发中经常需要对csv文件进行读写,下面这篇文章主要给大家介绍了关于Python 3.x读写csv文件中数字的相关资料,文中通过示例代码介绍的非常详细,对大家具有一定的参考学习价值,需要的朋友们下面跟着小编来一起学习学习吧。
收藏 0 赞 0 分享

Python实现解析Bit Torrent种子文件内容的方法

这篇文章主要介绍了Python实现解析Bit Torrent种子文件内容的方法,结合实例形式分析了Python针对Torrent文件的读取与解析相关操作技巧与注意事项,需要的朋友可以参考下
收藏 0 赞 0 分享

Python实现文件内容批量追加的方法示例

这篇文章主要介绍了Python实现文件内容批量追加的方法,结合实例形式分析了Python文件的读写相关操作技巧,需要的朋友可以参考下
收藏 0 赞 0 分享

Python简单实现自动删除目录下空文件夹的方法

这篇文章主要介绍了Python简单实现自动删除目录下空文件夹的方法,涉及Python针对文件与目录的读取、判断、删除等相关操作技巧,需要的朋友可以参考下
收藏 0 赞 0 分享

简单学习Python多进程Multiprocessing

这篇文章主要和大家一起简单的学习Python多进程Multiprocessing ,具有一定的参考价值,感兴趣的小伙伴们可以参考一下
收藏 0 赞 0 分享

Python导入模块时遇到的错误分析

这篇文章主要给大家详细解释了在Python处理导入模块的时候出现错误以及具体的情况分析,非常的详尽,有需要的小伙伴可以参考下
收藏 0 赞 0 分享
查看更多