### 龙贝格积分法的手动实现 龙贝格积分法是一种高效的数值积分方法，它通过梯形法则和理查森外推技术逐步提高积分精度。下面展示了一个完全不依赖 `numpy` 库的 Python 实现，并允许用户输入自定义被积函数。 --- ### 功能描述 1. **动态加载被积函数**：用户可通过字符串形式输入任意合法的数学表达式。 2. **核心逻辑分离**：分别实现了梯形法则、逐层细化以及最终的龙贝格表构建。 3. **无第三方库依赖**：整个实现仅使用标准库模块（如 `math`），确保轻量化。 --- ### 代码实现 以下是完整的代码实现： ```python import math def parse_function(func_str): """ 将用户输入的字符串转换为可调用的函数对象 """ def func(x): try: return eval(func_str.replace('e', 'math.e').replace('^', '**')) except Exception as e: raise ValueError(f"无法解析函数: {func_str}, 错误信息: {str(e)}") return func def trapezoidal_rule(f, a, b, n): """ 使用梯形法则计算基础积分值 """ h = (b - a) / n result = 0.5 * f(a) + 0.5 * f(b) for i in range(1, n): xi = a + i * h result += f(xi) return result * h def romberg_integration(f, a, b, tol=1e-8, max_iter=20): """ 主体龙贝格外推积分算法 """ R = [[trapezoidal_rule(f, a, b, 1)]] k = 1 while True: # 计算新的第零列项 new_R = [trapezoidal_rule(f, a, b, 2**k)] # 构建当前行剩余元素 for j in range(1, len(R)+1): term = (new_R[j-1] - R[j-1][j-1]) / ((4**j) - 1) value = new_R[j-1] + term new_R.append(value) # 判断是否满足收敛条件 if abs(new_R[-1] - R[-1][-1]) < tol or k >= max_iter: break # 更新表格数据结构 R = [row[:len(new_R)-i] for i, row in enumerate(R)] R.append(new_R) k += 1 return new_R[-1] if __name__ == "__main__": # 获取用户输入的函数表达式 func_expr = input("请输入被积函数 (例如 x^2*e^x 表示为 x**2*math.exp(x)): ") user_func = parse_function(func_expr) # 输入积分上下限 limits = list(map(float, input("请输入积分区间 [a, b], 格式为两个数空格隔开: ").split())) a, b = min(limits), max(limits) # 执行积分运算 integral_result = romberg_integration(user_func, a, b) print(f"龙贝格积分法得到的积分结果为: {integral_result:.7f}") ``` --- ### 关键功能解释 1. **函数解析器 (`parse_function`)** - 接收用户的字符串表达式并将其转化为实际可用的匿名函数[^1]。 - 支持常见数学操作符替换（如将幂次符号 `^` 替换为 Python 的 `**` 运算符）。 2. **梯形法则 (`trapezoidal_rule`)** - 提供最基础的一维积分近似方式，作为后续外推的基础[^2]。 3. **龙贝格积分主体逻辑 (`romberg_integration`)** - 循环执行梯形法则细分与理查森外推公式应用，直至达到预设误差容忍度或最大迭代次数限制[^3]。 4. **交互界面设计** - 用户能够自由指定待求解的函数及其作用域边界；程序自动适配不同类型的输入情况。 --- ### 测试案例 假设我们想对 \( y = x^2 e^x \) 在区间 [0, 1] 上进行积分，则按照提示依次填入如下内容即可获得相应答案： - 被积函数：`x**2*math.exp(x)` - 积分范围：`0 1` 预期输出应接近理论真值约等于 `1.96511`. --- ####