1. 原理与公式
2. Python代码实现
# coding=utf-8import numpy as npimport math# 二分法梯形公式# func:需要积分的函数# x_min: 积分下限# x_max: 积分上限# epoch: 二分次数def compute_Tn(func, x_min=0, x_max=1, epoch=10):Tn_list = []Tn = 0h0 = x_max - x_min # 积分区间的长度,即初始步长h = h0 # 每次迭代计算更新h步长x_half_list = np.array([0]) # 二分点列表for k in range(epoch + 1):if k == 0: # 初始步长h=h0,取两个端点Tn = 0.5 * h * (func(x_min) + func(x_max)) # T1Tn_list.append({"T_%d" % 2 ** k: Tn.item(), "k": k})x_half_list = np.array([(x_min + x_max) / 2]) # 计算二分点else:Tn = 0.5 * Tn + 0.5 * h * np.sum(func(x_half_list)) # 上一轮的T2n = 0.5*Tn + 0.5*h*二分点处的函数值之和Tn_list.append({"T_%d" % 2 ** k: Tn.item(), "k": k})h = 0.5 * h # 更新步长h为原来的一半x_half_list = np.linspace(0, 2 ** k, 2 ** k, endpoint=False) # 计算下一轮所需的二分点，一共有2^k个点,0,1,2,...2^k-1x_half_list = h0 * (2 * x_half_list + 1) / (2 ** (k + 1)) + x_min # X_(k+1/2)=a + (b-a)*(2n+1)/2^(k+1)return Tn_list# 由梯形公式计算辛普森公式def compute_Sn(Tn_list):Sn_list = []for i in range(len(Tn_list) - 1):Sn = list(Tn_list[i + 1].values())[0] * 4 / 3 - list(Tn_list[i].values())[0] / 3k = list(Tn_list[i + 1].values())[1]Sn_list.append({"S_%d" % 2 ** (k - 1): Sn, "k": k})return Sn_list# 由辛普森公式计算柯特斯公式def compute_Cn(Tn_list=None, Sn_list=None):if Sn_list is None:Sn_list = compute_Sn(Tn_list)Cn_list = []for i in range(len(Sn_list) - 1):Cn = list(Sn_list[i + 1].values())[0] * 16 / 15 - list(Sn_list[i].values())[0] / 15k = list(Sn_list[i + 1].values())[1]Cn_list.append({"C_%d" % 2 ** (k - 2): Cn, "k": k})return Cn_list# 由柯特斯公式计算龙贝格公式def compute_Rn(Tn_list=None, Cn_list=None):if Cn_list is None:Cn_list = compute_Cn(Tn_list)Rn_list = []for i in range(len(Cn_list) - 1):Rn = list(Cn_list[i + 1].values())[0] * 64 / 63 - list(Cn_list[i].values())[0] / 63k = list(Cn_list[i + 1].values())[1]Rn_list.append({"R_%d" % 2 ** (k - 3): Rn, "k": k})return Rn_list# 例题测试用函数def sinx_div_x(x):y = np.sin(x) / xif not isinstance(y, np.ndarray):y = np.array([y])y[np.where(np.isnan(y))] = 1return yif __name__ == '__main__':Tn_list = compute_Tn(sinx_div_x, x_min=0, x_max=1, epoch=24)print(Tn_list)Sn_list = compute_Sn(Tn_list)print(Sn_list)Cn_list = compute_Cn(Sn_list=Sn_list)print(Cn_list)Rn_list = compute_Rn(Cn_list=Cn_list)print(Rn_list)
{'T_1': 0.9207354924039483, 'k': 0}{'T_2': 0.9397932848061772, 'k': 1}{'T_4': 0.9445135216653896, 'k': 2}{'T_8': 0.9456908635827013, 'k': 3}{'T_16': 0.945985029934386, 'k': 4}{'T_32': 0.9460585609627681, 'k': 5}{'T_64': 0.9460769430600631, 'k': 6}{'T_128': 0.946081538543152, 'k': 7}{'T_256': 0.946082687411347, 'k': 8}{'T_512': 0.9460829746282348, 'k': 9}{'T_1024': 0.9460830464324467, 'k': 10}{'T_2048': 0.946083064383499, 'k': 11}{'T_4096': 0.946083068871262, 'k': 12}{'T_8192': 0.9460830699932028, 'k': 13}{'T_16384': 0.946083070273688, 'k': 14}{'T_32768': 0.9460830703438092, 'k': 15}{'T_65536': 0.9460830703613395, 'k': 16}{'T_131072': 0.9460830703657221, 'k': 17}{'T_262144': 0.9460830703668177, 'k': 18}{'T_524288': 0.9460830703670917, 'k': 19}{'T_1048576': 0.9460830703671603, 'k': 20}{'T_2097152': 0.9460830703671774, 'k': 21}{'T_4194304': 0.9460830703671814, 'k': 22}{'T_8388608': 0.9460830703671824, 'k': 23}{'T_16777216': 0.946083070367183, 'k': 24}{'S_1': 0.9461458822735869, 'k': 1}{'S_2': 0.9460869339517937, 'k': 2}{'S_4': 0.9460833108884718, 'k': 3}{'S_8': 0.9460830853849478, 'k': 4}{'S_16': 0.9460830713055621, 'k': 5}{'S_32': 0.9460830704258281, 'k': 6}{'S_64': 0.9460830703708483, 'k': 7}{'S_128': 0.9460830703674121, 'k': 8}{'S_256': 0.9460830703671974, 'k': 9}{'S_512': 0.9460830703671841, 'k': 10}{'S_1024': 0.9460830703671832, 'k': 11}{'S_2048': 0.946083070367183, 'k': 12}{'S_4096': 0.946083070367183, 'k': 13}{'S_8192': 0.9460830703671832, 'k': 14}{'S_16384': 0.946083070367183, 'k': 15}{'S_32768': 0.946083070367183, 'k': 16}{'S_65536': 0.946083070367183, 'k': 17}{'S_131072': 0.946083070367183, 'k': 18}{'S_262144': 0.946083070367183, 'k': 19}{'S_524288': 0.9460830703671832, 'k': 20}{'S_1048576': 0.9460830703671832, 'k': 21}{'S_2097152': 0.9460830703671828, 'k': 22}{'S_4194304': 0.9460830703671828, 'k': 23}{'S_8388608': 0.9460830703671831, 'k': 24}{'C_1': 0.9460830040636741, 'k': 2}{'C_2': 0.946083069350917, 'k': 3}{'C_4': 0.9460830703513795, 'k': 4}{'C_8': 0.9460830703669365, 'k': 5}{'C_16': 0.9460830703671791, 'k': 6}{'C_32': 0.946083070367183, 'k': 7}{'C_64': 0.946083070367183, 'k': 8}{'C_128': 0.9460830703671831, 'k': 9}{'C_256': 0.9460830703671832, 'k': 10}{'C_512': 0.9460830703671832, 'k': 11}{'C_1024': 0.946083070367183, 'k': 12}{'C_2048': 0.9460830703671831, 'k': 13}{'C_4096': 0.9460830703671833, 'k': 14}{'C_8192': 0.946083070367183, 'k': 15}{'C_16384': 0.9460830703671831, 'k': 16}{'C_32768': 0.9460830703671831, 'k': 17}{'C_65536': 0.9460830703671831, 'k': 18}{'C_131072': 0.9460830703671831, 'k': 19}{'C_262144': 0.9460830703671833, 'k': 20}{'C_524288': 0.9460830703671832, 'k': 21}{'C_1048576': 0.9460830703671828, 'k': 22}{'C_2097152': 0.9460830703671829, 'k': 23}{'C_4194304': 0.9460830703671831, 'k': 24}{'R_1': 0.9460830703872224, 'k': 3}{'R_2': 0.9460830703672599, 'k': 4}{'R_4': 0.9460830703671834, 'k': 5}{'R_8': 0.946083070367183, 'k': 6}{'R_16': 0.946083070367183, 'k': 7}{'R_32': 0.946083070367183, 'k': 8}{'R_64': 0.9460830703671831, 'k': 9}{'R_128': 0.9460830703671832, 'k': 10}{'R_256': 0.9460830703671832, 'k': 11}{'R_512': 0.946083070367183, 'k': 12}{'R_1024': 0.9460830703671831, 'k': 13}{'R_2048': 0.9460830703671833, 'k': 14}{'R_4096': 0.946083070367183, 'k': 15}{'R_8192': 0.9460830703671831, 'k': 16}{'R_16384': 0.9460830703671831, 'k': 17}{'R_32768': 0.9460830703671831, 'k': 18}{'R_65536': 0.9460830703671831, 'k': 19}{'R_131072': 0.9460830703671833, 'k': 20}{'R_262144': 0.9460830703671832, 'k': 21}{'R_524288': 0.9460830703671828, 'k': 22}{'R_1048576': 0.9460830703671829, 'k': 23}{'R_2097152': 0.9460830703671831, 'k': 24}Process finished with exit code 0
