python Romber求积
admin
2023-07-30 19:54:25
0

 

import numpy as np


# 被积函数
def fun(x):
    f = pow(x, 1.5)
    return f


# 求梯形值
def T2n(a, b, n, Tn):
    h = (b - a) / n  # 步长
    sum = 0.
    for k in range(n):
        sum += fun(a + (k + 0.5) * h)
    T2n = Tn / 2. + sum * h / 2.
    return T2n


# 求加速值
def romberg(max, a, b, epsilon):  # max为计算最大次数,a、b为积分下、上限,epsilon为限定误差
    tm = np.zeros(max, dtype=float)  # 第m行元素
    tm1 = np.zeros(max, dtype=float)  # 第m+1行元素
    tm[0] = (b - a) * (fun(a) + fun(b)) / 2.  # 初始值
    print(tm)
    k = 0
    err = 1
    while (err > epsilon and k < max - 1):  # 当误差小于预定误差,或计算次数大于最大次数时停止
        n = 2 ** k
        m = 1
        tm1[0] = T2n(a, b, n, tm[0])
        while (err > epsilon and m <= (k + 1)):  # 当误差小于预定误差,或本行全部计算完毕后停止
            tm1[m] = tm1[m - 1] + (tm1[m - 1] - tm[m - 1]) / (4. ** m - 1)
            result = tm1[m]
            err1 = abs(tm1[m] - tm[m - 1])  # 对角元素差
            err2 = abs(tm1[m] - tm1[m - 1])  # 同行前后两元素差
            err = min(err1, err2)
            m += 1
        tm = np.copy(tm1)  # 下移一行
        k += 1
        print(tm)
    return result


if __name__ == '__main__':
    f1 = romberg(6, 0, 1, 1.e-10)
    print(f1)

上一篇:贪心算法

下一篇:Selenium常用方法

相关内容

热门资讯

Mobi、epub格式电子书如... 在wps里全局设置里有一个文件关联,打开,勾选电子书文件选项就可以了。
定时清理删除C:\Progra... C:\Program Files (x86)下面很多scoped_dir开头的文件夹 写个批处理 定...
scoped_dir32_70... 一台虚拟机C盘总是莫名奇妙的空间用完,导致很多软件没法再运行。经过仔细检查发现是C:\Program...
500 行 Python 代码... 语法分析器描述了一个句子的语法结构,用来帮助其他的应用进行推理。自然语言引入了很多意外的歧义,以我们...
小程序支付时提示:appid和... [Q]小程序支付时提示:appid和mch_id不匹配 [A]小程序和微信支付没有进行关联,访问“小...
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 版本已于...
项目管理和工程管理的区别 项目管理 项目管理,顾名思义就是专注于开发和完成项目的管理,以实现目标并满足成功标准和项目要求。 工...