Romberg数值积分代码

import numpy as np


def f(x):
    """定义被积函数 f(x) = ln(1 + x) / (1 + x^2)"""
    return np.log(1 + x) / (1 + x ** 2)


def trapezoidal_rule(a, b, n, func):
    """梯形法则实现"""
    h = (b - a) / n
    x = np.linspace(a, b, n + 1)
    y = func(x)
    return (h / 2) * (y[0] + 2 * np.sum(y[1:n]) + y[n])


def romberg_integration(a, b, func, max_steps=4):
    """
    龙贝格积分实现
    参数：
    a, b: 积分上下限
    func: 被积函数
    max_steps: 最大迭代次数 (对应 n = 2^k)
    返回：龙贝格表格和最终近似值
    """
    # 初始化龙贝格表格
    R = np.zeros((max_steps, max_steps))

    # 计算 T_n (梯形法则)
    for k in range(max_steps):
        n = 2 ** k
        R[k, 0] = trapezoidal_rule(a, b, n, func)

    # 理查森外推
    for m in range(1, max_steps):
        for k in range(m, max_steps):
            R[k, m] = (4 ** m * R[k, m - 1] - R[k - 1, m - 1]) / (4 ** m - 1)

    # 打印表格
    print("龙贝格积分表格：")
    for i in range(max_steps):
        row = [f"{R[i, j]:.6f}" if j <= i else " " for j in range(max_steps)]
        print(f"k={i}: {row}")

    return R[max_steps - 1, max_steps - 1]


# 测试代码
if __name__ == "__main__":
    a, b = 0, 1  # 积分区间 [0, 1]
    max_steps = 4  # 迭代到 k=3 (n=1, 2, 4, 8)

    result = romberg_integration(a, b, f, max_steps)
    exact = (np.pi / 8) * np.log(2)  # 精确值
    print(f"\n龙贝格积分近似值: {result:.6f}")
    print(f"精确值: {exact:.6f}")
    print(f"误差: {abs(result - exact):.6f}")