Python龙贝格法求积分实例
我就废话不多说了,直接上代码吧!
#龙贝格法求积分 importmath a=0#积分下限 b=1#积分上限 eps=10**-5#精度 T=[]#复化梯形序列 S=[]#Simpson序列 C=[]#Cotes序列 R=[]#Romberg序列 deffunc(x):#被积函数 y=math.exp(-x) returny defRomberg(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 foriinrange(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) ifm>=1: S.append((4**m*T[-1]-T[-2])/(4**m-1)) ifm>=2: C.append((4**m*S[-1]-S[-2])/(4**m-1)) ifm>=3: R.append((4**m*C[-1]-C[-2])/(4**m-1)) ifm>4: ep=abs(10*(R[-1]-R[-2])) Romberg(a,b,eps,func) #print(T) #print(S) #print(C) #print(R) #计算机参考值0.6321205588 print("积分结果为:{:.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 令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编程代码如下: 以上这篇Python龙贝格法求积分实例就是小编分享给大家的全部内容了,希望能给大家一个参考,也希望大家多多支持毛票票。 声明:本文内容来源于网络,版权归原作者所有,内容由互联网用户自发贡献自行上传,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任。如果您发现有涉嫌版权的内容,欢迎发送邮件至:czq8825#qq.com(发邮件时,请将#更换为@)进行举报,并提供相关证据,一经查实,本站将立刻删除涉嫌侵权内容。
#coding=UTF-8
#Author:winyn
'''
给定一个函数,如:f(x)=x^(3/2),和积分上下限a,b,用机械求积Romberg公式求积分。
'''
importnumpyasnp
deffunc(x):
returnx**(3/2)
classRomberg:
def__init__(self,integ_dowlimit,integ_uplimit):
'''
初始化积分上限integ_uplimit和积分下限integ_dowlimit
输入一个函数,输出函数在积分上下限的积分
'''
self.integ_uplimit=integ_uplimit
self.integ_dowlimit=integ_dowlimit
defcalc(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)foriinrange(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
foriinrange(1,5):
sum=0
#多出来的点的累加和
foreachinrange(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)foriinrange(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)foriinrange(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)foriinrange(0,2)]
print('R序列:'+str(list(r_seq4)))
return'end'
rom=Romberg(0,1)
print(rom.calc())