Python 龙贝格/Romberg算法
admin
2023-07-30 19:52:53
0

龙贝格求积公式也称为逐次分半加速法。是数值计算方法之一,用以求解数值积分。是在梯形公式、辛普森公式和柯特斯公式之间关系的基础上,构造出一种加速计算积分的方法。 作为一种外推算法,在不增加计算量的前提下提高了误差的精度。

在等距基点的情况下,用计算机计算积分值通常都采用把区间逐次分半的方法进行。这样,前一次分割得到的函数值在分半以后仍可被利用,且易于编程。

# 用于储存函数,返回函数值
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])

相关内容

热门资讯

Mobi、epub格式电子书如... 在wps里全局设置里有一个文件关联,打开,勾选电子书文件选项就可以了。
定时清理删除C:\Progra... C:\Program Files (x86)下面很多scoped_dir开头的文件夹 写个批处理 定...
scoped_dir32_70... 一台虚拟机C盘总是莫名奇妙的空间用完,导致很多软件没法再运行。经过仔细检查发现是C:\Program...
小程序支付时提示:appid和... [Q]小程序支付时提示:appid和mch_id不匹配 [A]小程序和微信支付没有进行关联,访问“小...
500 行 Python 代码... 语法分析器描述了一个句子的语法结构,用来帮助其他的应用进行推理。自然语言引入了很多意外的歧义,以我们...
pycparser 是一个用... `pycparser` 是一个用 Python 编写的 C 语言解析器。它可以用来解析 C 代码并构...
微信小程序使用slider实现... 众所周知哈,微信小程序里面的音频播放是没有进度条的,但最近有个项目呢,客户要求音频要有进度条控制,所...
65536是2的几次方 计算2... 65536是2的16次方:65536=2⁶ 65536是256的2次方:65536=256 6553...
Apache Doris 2.... 亲爱的社区小伙伴们,我们很高兴地向大家宣布,Apache Doris 2.0.0 版本已于...
项目管理和工程管理的区别 项目管理 项目管理,顾名思义就是专注于开发和完成项目的管理,以实现目标并满足成功标准和项目要求。 工...