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])