python Romber求积

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)