我就廢話不多說了，直接上代碼吧！

# 龍貝格法求積分import matha=0 # 積分下限b=1 # 積分上限eps=10**-5 # 精度T=[] # 復(fù)化梯形序列S=[] # Simpson序列C=[] # Cotes序列R=[] # Romberg序列def func(x): # 被積函數(shù) 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)# 計(jì)算機(jī)參考值0.6321205588print('積分結(jié)果為：{:.5f}'.format(R[-1]))

補(bǔ)充拓展：python實(shí)現(xiàn)數(shù)值分析之龍貝格求積公式

復(fù)合梯形公式的提出：

1.首先，什么是梯形公式：

梯形公式表明：f(x)在[a,b]兩點(diǎn)之間的積分（面積），近似地可以用一個(gè)梯形的面積表示。

2.顯然，這個(gè)梯形公式對于不同的f(x)而言，其代數(shù)精度不同。為了能適合更多的f(x)，我們一般使用牛頓-科特斯公式其中比較高次的公式來進(jìn)行數(shù)值求積。但高次的缺陷是當(dāng)次數(shù)大于8次，求積公式就會不穩(wěn)定。因此，我們用于數(shù)值積分的牛頓-科特斯公式通常是一次的梯形公式、二次的辛普森公式和4此的科特斯公式。

辛普森公式：

科特斯公式：

3.牛頓-科特斯公式次數(shù)高于8次不能用，但是低次公式又精度不夠。解決辦法就是使用：復(fù)合梯形求積公式。復(fù)合求積公式就是在區(qū)間[a,b]上劃分n格小區(qū)間。一個(gè)大區(qū)間[a,b]上用一次梯形公式精度不夠，那么在n個(gè)小區(qū)間都使用梯形公式，最后將小區(qū)間的和累加起來，就可以得到整個(gè)大區(qū)間[a,b]的積分近似值。

a = x0 < x1 <x2 …<xn-1 < xn =b

令Tn為將[a,b]劃分n等分的復(fù)合梯形求積公式，h =(b-a)/n為小區(qū)間的長度。h/2類似于梯形公式中的（b-a）/2

注意：這里的k+1是下標(biāo)

通過研究我們發(fā)現(xiàn)：T2n與Tn之間存在一些遞推關(guān)系。

注意：這里的k+1/2是下標(biāo)。并且其中的h/2是中的h是Tn(n等分中的h = (b-a)/n))

于是乎，我們可以一次推出T1,T2,T4,T8…T2n序列

引出這些之后，才是我們的主題：龍貝格求積公式

龍貝格求積公式的實(shí)質(zhì)是用T2n序列構(gòu)造，S2n序列，

再用S2n序列構(gòu)造C2n序列

最后用C2n序列構(gòu)造R2n序列。

編程實(shí)現(xiàn)，理解下面的幾個(gè)公式即可。

python編程代碼如下：

# coding=UTF-8# Author:winyn’’’給定一個(gè)函數(shù)，如：f(x)= x^(3/2），和積分上下限a,b，用機(jī)械求積Romberg公式求積分。’’’import numpy as npdef func(x): return x**(3/2)class Romberg: def __init__(self, integ_dowlimit, integ_uplimit): ’’’ 初始化積分上限integ_uplimit和積分下限integ_dowlimit 輸入一個(gè)函數(shù)，輸出函數(shù)在積分上下限的積分 ’’’ self.integ_uplimit = integ_uplimit self.integ_dowlimit = integ_dowlimit def calc(self): ’’’ 計(jì)算Richardson外推算法的四個(gè)序列 ’’’ t_seq1 = np.zeros(5, ’f’) s_seq2 = np.zeros(4, ’f’) c_seq3 = np.zeros(3, ’f’) r_seq4 = np.zeros(2, ’f’) # 循環(huán)生成hm間距序列 hm = [(self.integ_uplimit - self.integ_dowlimit) / (2 ** i) for i in range(0,5)] print(hm) # 循環(huán)生成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 # 多出來的點(diǎn)的累加和 for each in range(1, 2**i,2): sum =sum + hm[i]*func( self.integ_dowlimit+each * hm[i])#計(jì)算兩項(xiàng)值 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))) # 循環(huán)生成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))) # 循環(huán)生成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))) # 循環(huán)生成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龍貝格法求積分實(shí)例就是小編分享給大家的全部內(nèi)容了，希望能給大家一個(gè)參考，也希望大家多多支持好吧啦網(wǎng)。