Python 龙贝格/Romberg算法
龙贝格求积公式也称为逐次分半加速法。是数值计算方法之一,用以求解数值积分。是在梯形公式、辛普森公式和柯特斯公式之间关系的基础上,构造出一种加速计算积分的方法。 作为一种外推算法,在不增加计算量的前提下提高了误差的精度。
在等距基点的情况下,用计算机计算积分值通常都采用把区间逐次分半的方法进行。这样,前一次分割得到的函数值在分半以后仍可被利用,且易于编程。
# 用于储存函数,返回函数值 def fun(x): return 1 / (x + 1) # 存放求积分范围 Min = float(input('请输入求积分下限:')) Max = float(input('请输入求积分上限:')) Mid = Max - Min # 存放积分结果的精度要求 precision = 10 ** int('-'+input('请输入求积精度10^-(输入值):')) # 用于存放 T S C R 的计算结果 T = [[j for j in [0, 0, 0, 0]] for i in range(4)] # 用于计行数 i = 1 T_1 = Mid * (fun(Max) + fun(Min)) / 2 T[0][0] = T_1 T[1][0] = T_1 / 2 + fun((Max - Min) / 2) * Mid / 2 T[1][1] = (T[1][0] * 4 ** i - T[0][0]) / (4 ** i - 1) # 精度达不到设定值时继续执行 # 当行数小于四时,根据行数选择计算 T S C R 中的哪几个 while abs(T[i][i]-T[i][i-1]) > precision: i += 1 Sum = 0 if i == 4 : break for j in range(2 ** (i - 1)): Sum += fun(Min + (2 * j + 1) * Mid / 2 ** i) T[i][0] = T[i - 1][0] / 2 + Sum * Mid / 2 ** i for j in range(1, i + 1): T[i][j] = (T[i][j - 1] * 4 ** i - T[i - 1][j - 1]) / (4 ** i - 1) # 当行数大于四时,计算全部的 T S C R while abs(T[i-1][-1] - T[i-1][-2]) > precision: Sum = 0 T.append([]) for j in range(2 ** (i - 1)): Sum += fun(Min + (2 * j + 1) * Mid / 2 ** i) T[i].append(T[i - 1][0] / 2 + Sum * Mid / 2 ** i) for j in range(1, 4): T[i].append((T[i][j - 1] * 4 ** i - T[i - 1][j - 1]) / (4 ** i - 1)) i += 1 print('函数在积分区间[{0},{1}]上,误差小于10^{2}的积分为:'.format(Min,Max,precision)) print(T[-1][-1])