複合梯形公式的提出:
1.首先,什麼是梯形公式:
![](https://img.laitimes.com/img/9ZDMuAjOiMmIsIjOiQnIsICM38CXlZHbvN3cpR2Lc1TPB10QGtWUCpEMJ9CXsxWam9CXwADNvwVZ6l2c052bm9CXUJDT1wkNhVzLcRnbvZ2Lc1TPB9EMZRUT5tmaNBDOsJGcohVYsR2MMBjVtJWd0ckW65UbM5WOHJWa5kHT20ESjBjUIF2LcRHelR3LcJzLctmch1mclRXY39zN2gzN1YzM1EDNyATM4EDMy8CX0Vmbu4GZzNmLn9Gbi1yZtl2Lc9CX6MHc0RHaiojIsJye.jpg)
梯形公式表明: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),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())