您好,欢迎来到保捱科技网。
搜索
您的当前位置:首页复化梯形求积分实例——用Python进行数值计算

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

来源:保捱科技网
复化梯形求积分实例——⽤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进⾏数值计算就是⼩编分享给⼤家的全部内容了,希望能给⼤家⼀个参考,也希望⼤家多多⽀持。

因篇幅问题不能全部显示,请点此查看更多更全内容

Copyright © 2019- baoaiwan.cn 版权所有 赣ICP备2024042794号-3

违法及侵权请联系:TEL:199 18 7713 E-MAIL:2724546146@qq.com

本站由北京市万商天勤律师事务所王兴未律师提供法律服务