复化梯形求积分——用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个小区间。

时间: 2024-10-06 15:50:25

复化梯形求积分——用Python进行数值计算的相关文章

python与数值计算环境搭建

数值计算的编程的软件很多种,也见过一些编程绘图软件的对比. 利用Python进行数值计算,需要用到numpy(矩阵) ,scipy(公式符号), matplotlib(绘图)这些工具包. 1.Linux系统中一般会带有Python.可以用命令查看是否安装Python $ python Python 2.7.5 (default, Feb 11 2014, 07:46:25) [GCC 4.8.2 20140120 (Red Hat 4.8.2-13)] on linux2 Type "help&

[算法][求积分][复合辛普森公式]

1 //这里f()为被积函数,输入a,b为积分上下限, 2 //eps为计算精度[这里要注意假收敛,一般设小一点好] 3 #include <iostream> 4 #include <cmath> 5 #define eps 1e-6 6 using namespace std; 7 double f(double x){ 8 return sqrt(1+cos(x)*cos(x)); 9 }//被积函数 10 double Sn(double a,double b,double

Python 3 数值计算

Python 3.4.3 (v3.4.3:9b73f1c3e601, Feb 24 2015, 22:43:06) [MSC v.1600 32 bit (Intel)] on win32Type "copyright", "credits" or "license()" for more information.>>> 17 /3 #典型的除法返回一个浮点数5.666666666666667>>> 17 //

Python实现数值计算----牛顿插值法

拉格朗日插值法的最大毛病就是每次引入一个新的插值节点,基函数都要发生变化,这在一些实际生产环境中是不合适的,有时候会不断的有新的测量数据加入插值节点集, 因此,通过寻找n个插值节点构造的的插值函数与n+1个插值节点构造的插值函数之间的关系,形成了牛顿插值法.推演牛顿插值法的方式是归纳法,也就是计算Ln(x)- Ln+1(x),并且从n=1开始不断的迭代来计算n+1时的插值函数. 牛顿插值法的公式是: 注意:在程序中我用W 代替  计算牛顿插值函数关键是要计算差商,n阶差商的表示方式如下:   关

用python做数值计算

http://sebug.net/paper/books/scipydoc/scipy_intro.html http://www.cnblogs.com/weilq/p/3432817.html https://eclipse.org/downloads/packages/eclipse-classic-422/junosr2 https://www.ics.uci.edu/~pattis/common/handouts/introtopythonineclipse/

Python实现数值计算----分段二次插值

事实上在实际使用中,高次插值显然是很不适合的,高次插值将所有样点包涵进一个插值函数中,这是次幂高的原因.高次计算复杂,而且刚开始的一点误差会被方的很大.因此将整个区间分为若干个小区间,在每一个小区间进行插值这样更好,实现容易,也方便在一些嵌入式设备上使用.有不少需要插值方法的场景是在嵌入式的应用中. 我以等距节点的二次插值为例,以每三个节点为一个子区间. 等距节点二次插值很好写,由于每个区间只有三个插值节点,计算差商也不必使用拉格朗日插值中使用的递归,直接列表达式也很简单(实际上等距节点二次插值

复化梯形求积法

#include<stdio.h>#include<string.h>#include<math.h>#include<conio.h>#include<stdlib.h>#define epsilon 0.00001float f(float x){ return (1.0/(1.0+x*x));}float computeT(float a,float b){ float T=0,h=(b-a)/2; T=h*(f(a)+2*T+f(b))/

猴子分桃—Python

def f(): for i in range(3120,4000): flag = 1 k=i for j in range(5): if i%5==1: i=(i//5)*4 else: flag=-1 break if flag==1: print(k) f() 原文地址:https://www.cnblogs.com/Python-T/p/9426974.html

第六章实验报告(2)

C程序设计实验报告 实验项目: 6.4.2.2.利用复化梯形公式计算定积分 6.4.2.3.编计算Ackerman函数 6.4.3.1.编写计算x的y次幂的递归函数getpower(int x,int y).并在主程序中实现输入输出. 6.4.3.2.编写计算学生年龄的递归函数 6.4.3.3.编写递归函数下hi先Ackerma函数 姓名:钟俊敏    实验地点:教学楼514教室     实验时间:5月16日 6.4.2.2.利用复化梯形公式计算定积分 ● 掌握c语言中定义函数的方法 ●掌握通过