时间:2021-05-22
我就废话不多说了,直接上代码吧!
# 龙贝格法求积分import matha=0 # 积分下限b=1 # 积分上限eps=10**-5 # 精度T=[] # 复化梯形序列S=[] # Simpson序列C=[] # Cotes序列R=[] # Romberg序列def func(x): # 被积函数 y=math.exp(-x) return ydef Romberg(a,b,eps,func): h = b - a T.append(h * (func(a) + func(b)) / 2) ep=eps+1 m=0 while(ep>=eps): m=m+1 t=0 for i in range(2**(m-1)-1): t=t+func(a+(2*(i+1)-1)*h/2**m)*h/2**m t=t+T[-1]/2 T.append(t) if m>=1: S.append((4**m*T[-1]-T[-2])/(4**m-1)) if m>=2: C.append((4**m*S[-1]-S[-2])/(4**m-1)) if m>=3: R.append((4**m*C[-1]-C[-2])/(4**m-1)) if m>4: ep=abs(10*(R[-1]-R[-2]))Romberg(a,b,eps,func)# print(T)# print(S)# print(C)# print(R)# 计算机参考值0.6321205588print("积分结果为:{:.5f}".format(R[-1]))补充拓展:python实现数值分析之龙贝格求积公式
复合梯形公式的提出:
1.首先,什么是梯形公式:
梯形公式表明:f(x)在[a,b]两点之间的积分(面积),近似地可以用一个梯形的面积表示。
2.显然,这个梯形公式对于不同的f(x)而言,其代数精度不同。为了能适合更多的f(x),我们一般使用牛顿-科特斯公式其中比较高次的公式来进行数值求积。但高次的缺陷是当次数大于8次,求积公式就会不稳定。因此,我们用于数值积分的牛顿-科特斯公式通常是一次的梯形公式、二次的辛普森公式和4此的科特斯公式。
辛普森公式:
科特斯公式:
3.牛顿-科特斯公式次数高于8次不能用,但是低次公式又精度不够。解决办法就是使用:复合梯形求积公式。复合求积公式就是在区间[a,b]上划分n格小区间。一个大区间[a,b]上用一次梯形公式精度不够,那么在n个小区间都使用梯形公式,最后将小区间的和累加起来,就可以得到整个大区间[a,b]的积分近似值。
a = x0 < x1 <x2 …<xn-1 < xn =b
令Tn为将[a,b]划分n等分的复合梯形求积公式,h =(b-a)/n为小区间的长度。h/2类似于梯形公式中的(b-a)/2
注意:这里的k+1是下标
通过研究我们发现:T2n与Tn之间存在一些递推关系。
注意:这里的k+1/2是下标。并且其中的h/2是中的h是Tn(n等分中的h = (b-a)/n))
于是乎,我们可以一次推出T1,T2,T4,T8…T2n序列
引出这些之后,才是我们的主题:龙贝格求积公式
龙贝格求积公式的实质是用T2n序列构造,S2n序列,
再用S2n序列构造C2n序列
最后用C2n序列构造R2n序列。
编程实现,理解下面的几个公式即可。
python编程代码如下:
# coding=UTF-8# Author:winyn'''给定一个函数,如:f(x)= x^(3/2),和积分上下限a,b,用机械求积Romberg公式求积分。'''import numpy as npdef func(x): return x**(3/2)class Romberg: def __init__(self, integ_dowlimit, integ_uplimit): ''' 初始化积分上限integ_uplimit和积分下限integ_dowlimit 输入一个函数,输出函数在积分上下限的积分 ''' self.integ_uplimit = integ_uplimit self.integ_dowlimit = integ_dowlimit def calc(self): ''' 计算Richardson外推算法的四个序列 ''' t_seq1 = np.zeros(5, 'f') s_seq2 = np.zeros(4, 'f') c_seq3 = np.zeros(3, 'f') r_seq4 = np.zeros(2, 'f') # 循环生成hm间距序列 hm = [(self.integ_uplimit - self.integ_dowlimit) / (2 ** i) for i in range(0,5)] print(hm) # 循环生成t_seq1 fa = func(self.integ_dowlimit) fb = func(self.integ_uplimit) t0 = (1 / 2) * (self.integ_uplimit - self.integ_dowlimit) * (fa+fb) t_seq1[0] = t0 for i in range(1, 5): sum = 0 # 多出来的点的累加和 for each in range(1, 2**i,2): sum =sum + hm[i]*func( self.integ_dowlimit+each * hm[i])#计算两项值 temp1 = 1 / 2 * t_seq1[i - 1] temp2 =sum temp = temp1 + temp2 # 求t_seql的1-4位 t_seq1[i] = temp print('T序列:'+ str(list(t_seq1))) # 循环生成s_seq2 s_seq2 = [round((4 * t_seq1[i + 1] - t_seq1[i]) / 3,6) for i in range(0, 4)] print('S序列:' + str(list(s_seq2))) # 循环生成c_seq3 c_seq3 = [round((4 ** 2 * s_seq2[i + 1] - s_seq2[i]) / (4 ** 2 - 1),6) for i in range(0, 3)] print('C序列:' + str(list(c_seq3))) # 循环生成r_seq4 r_seq4 = [round((4 ** 3 * c_seq3[i + 1] - c_seq3[i]) / (4 ** 3 - 1),6) for i in range(0, 2)] print('R序列:' + str(list(r_seq4))) return 'end'rom = Romberg(0, 1)print(rom.calc())以上这篇Python龙贝格法求积分实例就是小编分享给大家的全部内容了,希望能给大家一个参考,也希望大家多多支持。
声明:本页内容来源网络,仅供用户参考;我单位不保证亦不表示资料全面及准确无误,也不保证亦不表示这些资料为最新信息,如因任何原因,本网内容或者用户因倚赖本网内容造成任何损失或损害,我单位将不会负任何法律责任。如涉及版权问题,请提交至online#300.cn邮箱联系删除。
本文实例讲述了python实现数值积分的Simpson方法。分享给大家供大家参考。具体如下:#coding=utf-8#simpson法计算积分,数值积分,效果
本文实例为大家分享了C语言求解定积分的具体方法,供大家参考,具体内容如下题目要求:求下面函数的定积分:思路:求一个函数的定积分,其实就是求它的面积,如对函数求积
首先是对一元函数求积分,使用Scipy下的integrate函数:fromscipyimportintegratedefg(x):return(1-x**2)*
12月3日消息《》获悉,西贝创始人贾国龙在首届中国餐饮品牌节上表示,西贝已经决定上市,目前正在寻找合适的时间,选择资本投资西贝。据悉,此前,贾国龙曾对外界宣称:
本文实例讲述了Python实现的北京积分落户数据分析。分享给大家供大家参考,具体如下:北京积分落户状况获取数据(爬虫/文件下载)—>分析(维度—指标)从公司维度