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

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

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

插值函数一般是一个不超过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进行数值计算就是小编分享给大家的全部内容了,希望能给大家一个参考,也希望大家多多支持我们。

(0)

相关推荐

  • python用quad、dblquad实现一维二维积分的实例详解

    背景: python函数库scipy的quad.dblquad实现一维二维积分的范例.需要注意dblquad的积分顺序问题. 代码: import numpy as np from scipy import integrate def half_circle(x): """ 原心:(1,0),半径为1 半圆函数:(x-1)^2+y^2 = 1 """ return (1-(x-1)**2)**0.5 """ 梯形法求

  • Python 做曲线拟合和求积分的方法

    这是一个由加油站油罐传感器测量的油罐高度数据和出油体积,根据体积和高度的倒数,用截面积来描述油罐形状,求出拟合曲线,再用标准数据,求积分来验证拟合曲线效果和误差的一个小项目. 主要的就是首先要安装Anaconda  python库,然后来运用这些数学工具. ###最小二乘法试验### import numpy as np import pymysql from scipy.optimize import leastsq from scipy import integrate ###绘图,看拟合效

  • 利用python求积分的实例

    python的numpy库集成了很多的函数.利用其中的函数可以很方便的解决一些数学问题.本篇介绍如何使用python的numpy来求解积分. 代码如下: # -*- coding: utf-8 -*- import numpy as np from scipy.integrate import quad,dblquad,nquad def main(): print quad(lambda x:np.exp(-x),0,np.inf) '''求积分,np.inf代表正无穷. 结果第一个数值代表运

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

    用程序来求积分的方法有很多,这篇文章主要是有关牛顿-科特斯公式. 学过插值算法的同学最容易想到的就是用插值函数代替被积分函数来求积分,但实际上在大部分场景下这是行不通的. 插值函数一般是一个不超过n次的多项式,如果用插值函数来求积分的话,就会引进高次多项式求积分的问题.这样会将原来的求积分问题带到另一个求积分问题:如何求n次多项式的积分,而且当次数变高时,会出现龙悲歌现象,误差反而可能会增大,并且高次的插值求积公式有可能会变得不稳定:详细原因不赘述. 牛顿-科特斯公式解决这一问题的办法是将大的插

  • Python龙贝格法求积分实例

    我就废话不多说了,直接上代码吧! # 龙贝格法求积分 import math a=0 # 积分下限 b=1 # 积分上限 eps=10**-5 # 精度 T=[] # 复化梯形序列 S=[] # Simpson序列 C=[] # Cotes序列 R=[] # Romberg序列 def func(x): # 被积函数 y=math.exp(-x) return y def Romberg(a,b,eps,func): h = b - a T.append(h * (func(a) + func(

  • 使用Python中的reduce()函数求积的实例

    编写一个prod()函数,可以接受一个list并利用reduce()求积. from functools import reduce def prod(x,y): return x * y L = reduce(prod,[3,5,7,9]) print(L) 打印结果如下: 以上这篇使用Python中的reduce()函数求积的实例就是小编分享给大家的全部内容了,希望能给大家一个参考,也希望大家多多支持我们.

  • 通过实例了解python property属性

    这篇文章主要介绍了通过实例了解python property属性,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友可以参考下 1. 什么是property属性 一种用起来像是使用的实例属性一样的特殊属性,可以对应于某个方法 # ############### 定义 ############### class Foo: def func(self): pass # 定义property属性 @property def prop(self): pass # ##

  • 通过实例解析Python调用json模块

    这篇文章主要介绍了通过实例解析Python调用json模块,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友可以参考下 介绍 今天介绍一种数据格式,json.Json是JavaScript Object Notation的缩写,区别于txt.csv,json编码格式更加灵活,在工作也会经常遇到.在Python中要读写json是十分方便的,只需要调用json模块. 使用 直接导入模块 import json 两个读写数据的函数: json.dumps() 和

  • 关于初始种子自动选取的区域生长实例(python+opencv)

    算法中,初始种子可自动选择(通过不同的划分可以得到不同的种子,可按照自己需要改进算法),图分别为原图(自己画了两笔为了分割成不同区域).灰度图直方图.初始种子图.区域生长结果图. 另外,不管时初始种子选择还是区域生长,阈值选择很重要. import cv2 import numpy as np import matplotlib.pyplot as plt #初始种子选择 def originalSeed(gray, th): ret, thresh = cv2.cv2.threshold(gr

  • 通过实例了解Python str()和repr()的区别

    这篇文章主要介绍了通过实例了解Python str()和repr()的区别,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友可以参考下 区别 其实用处就是最大的区别了:str()主要用来为终端用户输出一些信息,而repr()主要用来调试:同时后者的目标是为了消除一些歧义(例如浮点数的精度问题),前者主要为了可读. 使用 In [12]: s = 'abc' In [13]: print(str(s)) abc In [14]: print(2.0/11) 0

  • 通过实例解析python描述符原理作用

    这篇文章主要介绍了通过实例解析python描述符原理作用,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友可以参考下 本质上看,描述符是一个类,只不过它定义了另一个类中属性的访问方式.换句话说,一个类可以将属性管理全权委托给描述符类. 描述符类基于以下三种特殊方法,换句话说,这三种方法组成了描述符协议: __set__(self, obj, type = None): 在设置属性时,将调用这一方法. __get__(self, obj, value): 在读

随机推荐