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

# 龍貝格法求積分
import math
a=0    # 積分下限
b=1    # 積分上限
eps=10**-5  # 精度
T=[]   # 復化梯形序列
S=[]   # Simpson序列
C=[]   # Cotes序列
R=[]   # Romberg序列
def func(x): # 被積函式
 y=math.exp(-x)
 return y
def 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,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 < 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 np


def 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),3)]
  print('C序列：' + str(list(c_seq3)))
  # 迴圈生成r_seq4
  r_seq4 = [round((4 ** 3 * c_seq3[i + 1] - c_seq3[i]) / (4 ** 3 - 1),2)]
  print('R序列：' + str(list(r_seq4)))
  return 'end'


rom = Romberg(0,1)
print(rom.calc())

以上這篇Python龍貝格法求積分例項就是小編分享給大家的全部內容了，希望能給大家一個參考，也希望大家多多支援我們。