龙贝格求积算法python实现

import numpy as np


def trapezoid(a, b, n, func):
    """
    复化梯形公式求函数func在区间[a,b]上的积分值
    n是等分的区间数目
    """
    x = np.linspace(a, b, num=n + 1)
    y = func(x)
    h = (b - a) / (2 * n)
    return h * (y[0] + 2 * np.sum(y[1:-1]) + y[-1])


def romberg(a, b, func, eps, max_iter=100):
    """
    龙贝格求积算法
    求函数func在区间[a,b]上的积分值
    eps:误差
    max_iter:最大迭代上限，应该大于等于2
    """
    previous = [trapezoid(a, b, 1, func)]
    for i in range(1, max_iter):
        current = [trapezoid(a, b, 2 ** i, func)]
        for k in range(1, i + 1):
            tmp = ((4 ** k) * current[k - 1] - previous[k - 1]) / (4 ** k - 1)
            current.append(tmp)
        if np.abs(current[-1] - previous[-1]) < eps:
            return current[-1]
        previous = current
    return previous[-1]


if __name__ == '__main__':
    # x^2+x^3+exp(x)在[0,2]上的积分值：13.055722765711728
    print(
        romberg(0, 2, lambda x: x ** 2 + x ** 3 + np.exp(x), 1e-5)
    )