Python数值积分实战：复合梯形vs复合辛普森公式性能对比（附可视化代码）

Python数值积分实战：复合梯形vs复合辛普森公式性能对比（附可视化代码）

在工程计算、物理模拟乃至金融建模中，我们常常会遇到一个看似简单却至关重要的任务：计算一个函数曲线下的面积，也就是积分。很多时候，我们面对的函数没有简单的解析原函数，或者其表达式复杂到难以手工处理。这时，数值积分就成了我们手中不可或缺的“瑞士军刀”。对于刚接触这个领域的朋友，可能会被各种公式弄得眼花缭乱——梯形法、辛普森法、龙贝格法，到底该选哪个？而对于有经验的老手，在面对一个具体的新问题时，可能也需要快速评估：在有限的算力和精度要求下，哪种方法能以最小的代价给出最可靠的结果？

今天，我们就聚焦于数值积分中最经典、应用最广泛的两种方法：复合梯形公式和复合辛普森公式。我们不打算停留在教科书式的理论推导，而是直接切入实战，用Python作为我们的实验平台。我们将通过一系列精心设计的对比实验，直观地展示这两种方法在处理不同特性函数（如平滑的多项式、振荡的三角函数）时的精度差异、收敛速度以及计算效率。更重要的是，我会提供完整的、可直接复用的代码，并教你如何生成可视化图表，让抽象的误差和性能对比变得一目了然。无论你是正在完成数值分析作业的学生，还是需要在项目中快速集成可靠积分模块的工程师，这篇文章都将为你提供清晰的路径和实用的工具。

1. 核心概念与算法原理速览

在深入代码之前，我们有必要快速回顾一下这两种方法的核心理念。理解其背后的思想，能帮助我们在实际应用中做出更明智的选择，而不仅仅是机械地调用函数。

复合梯形公式的思想最为直观：它将整个积分区间 [a, b] 分割成 n 个等宽的小区间。在每个小区间上，我们用连接区间两端点的直线（即梯形）来近似代替原函数曲线，然后计算所有这些梯形面积之和作为积分的近似值。它的公式简洁明了：

T_n = h/2 * [f(a) + 2 * Σ_{i=1}^{n-1} f(x_i) + f(b)]

其中，h = (b - a)/n 是步长，x_i = a + i*h。你可以把它想象成用一系列“台阶”去拟合曲线，当台阶足够细密时，拟合效果就会变好。

复合辛普森公式则更进了一步。它要求将区间分割成偶数 n 个子区间。在每两个相邻的子区间（即一个宽度为 2h 的大区间）上，它不是用直线，而是用一个二次抛物线来拟合函数曲线。这个抛物线会穿过区间两个端点和中点的函数值，从而能更好地捕捉函数的弯曲特性。其公式为：

S_n = h/3 * [f(a) + 4 * Σ_{i=1,3,5...}^{n-1} f(x_i) + 2 * Σ_{j=2,4,6...}^{n-2} f(x_j) + f(b)]

从理论上讲，如果被积函数足够光滑（即高阶导数连续且有限），复合辛普森公式的误差阶为 O(h^4)，而复合梯形公式的误差阶仅为 O(h^2)。这意味着，当我们将步长 h 减半时，梯形法的误差大约减少到原来的1/4，而辛普森法的误差能减少到原来的约1/16！这是一个巨大的优势。

为了更清晰地对比两者的计算模式，我们可以看下面的表格：

提示：误差阶 O(h^k) 是一个渐进概念，它描述了当步长 h 趋近于0时，误差减少的速度。k 越大，方法对于光滑函数的收敛速度越快。但这并不意味着在任何情况下辛普森法都绝对更好，对于不光滑或奇异的函数，高阶方法可能反而不稳定。

2. 从零搭建Python积分工具箱

理论说得再多，不如动手写一行代码。让我们抛开复杂的库，从最基础的原理开始，构建我们自己的数值积分函数。这不仅能加深理解，也让你拥有完全可控、可定制的工具。

首先，我们导入必要的库。NumPy 负责高效的数组运算，Matplotlib 则是我们进行可视化的得力助手。

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid') # 使用更美观的绘图样式

2.1 实现复合梯形公式

复合梯形公式的实现非常直接，几乎就是公式的直译。我们只需要注意数组切片的高效性。

def composite_trapezoidal(f, a, b, n):
    """
    使用复合梯形公式计算定积分近似值。

    参数
    ----------
    f : function
        被积函数，应能接受数组输入。
    a, b : float
        积分下限和上限。
    n : int
        子区间数量。

    返回
    -------
    float
        积分的近似值。
    """
    if n <= 0:
        raise ValueError("区间数 n 必须为正整数。")
    h = (b - a) / n
    x = np.linspace(a, b, n+1)  # 生成 n+1 个节点
    y = f(x)
    # 公式核心：首尾项权重为1，中间项权重为2
    T = h * (0.5*y[0] + np.sum(y[1:-1]) + 0.5*y[-1])
    return T

2.2 实现复合辛普森公式

复合辛普森公式的实现需要一点技巧，主要是对奇数下标和偶数下标的函数值进行不同的加权求和。利用 NumPy 的数组切片，我们可以写得既简洁又高效。

def composite_simpson(f, a, b, n):
    """
    使用复合辛普森公式计算定积分近似值。

    参数
    ----------
    f : function
        被积函数，应能接受数组输入。
    a, b : float
        积分下限和上限。
    n : int
        子区间数量，必须为偶数。

    返回
    -------
    float
        积分的近似值。
    """
    if n <= 0 or n % 2 != 0:
        raise ValueError("区间数 n 必须为正偶数。")
    h = (b - a) / n
    x = np.linspace(a, b, n+1)
    y = f(x)
    # 核心加权求和：奇数下标乘4，偶数下标乘2（不包括首尾）
    S = (h / 3) * (y[0] + y[-1] + 4 * np.sum(y[1:-1:2]) + 2 * np.sum(y[2:-2:2]))
    return S

注意：composite_simpson 函数中，y[1:-1:2] 选取了所有奇数索引的节点（从第1个到倒数第2个，步长为2），y[2:-2:2] 选取了所有偶数索引的内部节点（从第2个到倒数第3个，步长为2）。这种写法比用循环快得多。

有了这两个核心函数，我们的工具箱就准备好了。接下来，让我们用它们去解决一些实际问题，并看看它们的表现究竟如何。

3. 实战对比：当方法遇到不同性格的函数

数值方法的性能高度依赖于被积函数的“性格”。一个在多项式上表现完美的方法，面对剧烈振荡或存在奇点的函数时可能会惨败。我们设计三个具有代表性的测试案例，来全面考察我们的两种方法。

3.1 案例一：平滑多项式 - 基本功测试

我们首先选择一个非常平滑的函数：f(x) = x^4 + 2*x^2 + 1，在区间 [0, 2] 上积分。这个函数的高阶导数仍然连续且有限，是数值方法的“理想客户”。我们可以轻松求出其精确积分值为 (1/5)*2^5 + (2/3)*2^3 + 2 ≈ 14.9333...。

让我们编写一个对比函数，它能够计算不同区间数 n 下的近似值、绝对误差，并绘制误差随 n 变化的收敛图。

def compare_methods_on_polynomial():
    """在平滑多项式函数上对比两种方法"""
    f = lambda x: x**4 + 2*x**2 + 1
    a, b = 0, 2
    exact = (1/5)*b**5 + (2/3)*b**3 + b - ((1/5)*a**5 + (2/3)*a**3 + a)

    # 测试不同的区间数
    n_values = np.array([4, 8, 16, 32, 64, 128, 256])
    errors_trap = []
    errors_simp = []

    for n in n_values:
        I_trap = composite_trapezoidal(f, a, b, int(n))
        I_simp = composite_simpson(f, a, b, int(n))
        errors_trap.append(abs(I_trap - exact))
        errors_simp.append(abs(I_simp - exact))

    # 输出结果表格
    print(f"精确积分值: {exact:.10f}")
    print("\n不同区间数下的积分结果与误差对比：")
    print("-" * 70)
    print(f"{'n':>6} | {'梯形法近似值':>18} | {'梯形法误差':>12} | {'辛普森法近似值':>18} | {'辛普森法误差':>12}")
    print("-" * 70)
    for i, n in enumerate(n_values):
        I_trap = composite_trapezoidal(f, a, b, int(n))
        I_simp = composite_simpson(f, a, b, int(n))
        print(f"{n:6d} | {I_trap:18.10f} | {abs(I_trap-exact):12.4e} | {I_simp:18.10f} | {abs(I_simp-exact):12.4e}")

    # 绘制误差收敛图
    plt.figure(figsize=(10, 6))
    plt.loglog(n_values, errors_trap, 'o-', label='复合梯形公式误差', linewidth=2, markersize=8)
    plt.loglog(n_values, errors_simp, 's-', label='复合辛普森公式误差', linewidth=2, markersize=8)

    # 绘制参考斜率线，展示理论误差阶
    ref_h = 1.0 / n_values
    plt.loglog(n_values, 10 * ref_h**2, 'k--', label='斜率 O(1/n²)', alpha=0.7)
    plt.loglog(n_values, 0.1 * ref_h**4, 'k:', label='斜率 O(1/n⁴)', alpha=0.7)

    plt.xlabel('区间数量 n (对数坐标)', fontsize=12)
    plt.ylabel('绝对误差 (对数坐标)', fontsize=12)
    plt.title('平滑多项式函数上的误差收敛对比', fontsize=14)
    plt.legend(fontsize=11)
    plt.grid(True, which="both", ls="--", alpha=0.5)
    plt.tight_layout()
    plt.show()

# 运行测试
compare_methods_on_polynomial()

运行这段代码，你会立刻从输出的表格和图表中获得深刻印象。对于这个平滑函数，即使只用很少的区间（如n=8），辛普森法给出的结果误差已经非常小（可能达到1e-5量级），而梯形法需要更多的区间（如n=64）才能达到相近的精度。在双对数坐标的误差图上，两条线的斜率清晰地印证了理论：梯形法误差线的斜率接近-2（对应 O(1/n²)），而辛普森法误差线的斜率接近-4（对应 O(1/n⁴)）。这意味着为了将误差降低一个数量级，梯形法需要将区间数增加约3倍，而辛普森法只需要增加不到2倍。

3.2 案例二：振荡三角函数 - 稳定性挑战

现在我们来点更有挑战性的：一个高频振荡的函数，例如 f(x) = sin(10*x) + cos(5*x)，在区间 [0, π] 上积分。这类函数在工程中非常常见，比如信号处理。其精确积分值可以通过解析计算得到：∫[0,π] sin(10x)+cos(5x) dx = [-cos(10x)/10 + sin(5x)/5]_0^π = (1/5)*(sin(5π)-sin(0)) - (1/10)*(cos(10π)-cos(0)) = 0 - (1/10)*(1 - 1) = 0。

积分值为0，这对数值方法是一个有趣的考验，因为正负误差容易抵消，也可能因为微小的偏差而产生显著相对误差。

def compare_methods_on_oscillatory():
    """在高频振荡函数上对比两种方法"""
    f = lambda x: np.sin(10*x) + np.cos(5*x)
    a, b = 0, np.pi
    exact = 0.0  # 解析解

    n_values = np.array([10, 20, 40, 80, 160, 320, 640])
    errors_trap = []
    errors_simp = []

    for n in n_values:
        # 确保辛普森法的n为偶数
        n_even = int(n) if int(n) % 2 == 0 else int(n) + 1
        I_trap = composite_trapezoidal(f, a, b, int(n))
        I_simp = composite_simpson(f, a, b, n_even)
        errors_trap.append(abs(I_trap - exact))
        errors_simp.append(abs(I_simp - exact))

    # 可视化：绘制函数图像及两种方法的近似积分区域（以n=20为例）
    plt.figure(figsize=(14, 5))

    # 子图1：函数图像与近似积分区域
    plt.subplot(1, 2, 1)
    x_fine = np.linspace(a, b, 1000)
    y_fine = f(x_fine)
    plt.plot(x_fine, y_fine, 'b-', label='f(x) = sin(10x)+cos(5x)', linewidth=1.5)

    n_vis = 20
    x_trap = np.linspace(a, b, n_vis + 1)
    y_trap = f(x_trap)
    plt.fill_between(x_trap, 0, y_trap, alpha=0.3, color='orange', label='梯形法近似面积 (n=20)')
    plt.scatter(x_trap, y_trap, color='red', s=20, zorder=5)

    x_simp = np.linspace(a, b, n_vis + 1) # 节点相同
    # 为了展示辛普森的抛物线拟合，我们绘制更密的点
    for i in range(0, n_vis, 2):
        xs = np.linspace(x_simp[i], x_simp[i+2], 50)
        # 在每两个区间上用二次插值（这里简化用三次样条示意辛普森的思想）
        from scipy.interpolate import CubicSpline
        cs = CubicSpline(x_simp[i:i+3], y_trap[i:i+3])
        plt.plot(xs, cs(xs), 'g--', linewidth=1, alpha=0.7)
    plt.title('振荡函数与积分近似 (n=20)', fontsize=12)
    plt.xlabel('x')
    plt.ylabel('f(x)')
    plt.legend(loc='upper right')
    plt.grid(True, alpha=0.3)

    # 子图2：误差收敛对比
    plt.subplot(1, 2, 2)
    plt.loglog(n_values, errors_trap, 'o-', label='复合梯形公式误差', linewidth=2)
    plt.loglog(n_values, errors_simp, 's-', label='复合辛普森公式误差', linewidth=2)
    plt.xlabel('区间数量 n (对数坐标)')
    plt.ylabel('绝对误差 (对数坐标)')
    plt.title('振荡函数上的误差收敛对比')
    plt.legend()
    plt.grid(True, which="both", ls="--", alpha=0.5)

    plt.tight_layout()
    plt.show()

    # 打印部分结果
    print("振荡函数积分（精确值=0）误差对比：")
    print("-" * 60)
    print(f"{'n':>6} | {'梯形法误差':>15} | {'辛普森法误差':>15}")
    print("-" * 60)
    for i, n in enumerate(n_values[:4]): # 只看前几个
        n_even = int(n) if int(n) % 2 == 0 else int(n) + 1
        I_trap = composite_trapezoidal(f, a, b, int(n))
        I_simp = composite_simpson(f, a, b, n_even)
        print(f"{n:6d} | {abs(I_trap):15.6e} | {abs(I_simp):15.6e}")

compare_methods_on_oscillatory()

运行这个测试，你会发现一些有趣的现象。对于振荡函数，两种方法的收敛速度可能不如在平滑多项式上那么理想，误差曲线可能不会严格遵循理论斜率。辛普森法的优势依然存在，但可能没那么悬殊。左图直观展示了原因：梯形法用直线去拟合快速波动的曲线，在每个波峰波谷都会丢失大量面积信息。而辛普森法用抛物线拟合，能更好地跟随曲线的弯曲趋势，因此在相同区间数下，其面积近似（绿色虚线围成的区域）更接近真实曲线下的净面积（正负相抵后接近0）。

3.3 案例三：存在快速变化区域的函数 - 局部适应性考验

最后，我们看一个在局部区域变化剧烈的函数，比如 f(x) = exp(-50*(x-0.5)**2) + 0.1，这是一个中心在0.5处非常尖锐的高斯峰加上一个基底。在 [0, 1] 区间上积分。这个函数在大部分区域几乎为常数0.1，但在x=0.5附近有一个很窄的尖峰。精确积分值没有简单表达式，我们可以用 SciPy 的高精度积分函数 quad 作为参考基准。

from scipy import integrate

def compare_methods_on_peak():
    """在存在尖锐峰值的函数上对比两种方法"""
    f = lambda x: np.exp(-50 * (x - 0.5)**2) + 0.1
    a, b = 0, 1

    # 使用高精度方法计算近似精确值
    exact, _ = integrate.quad(f, a, b, epsabs=1e-14, epsrel=1e-14)
    print(f"参考精确值 (SciPy quad): {exact:.12f}")

    # 测试从稀疏到密集的区间划分
    n_values = np.array([5, 10, 20, 40, 80, 160, 320, 640])
    errors_trap = []
    errors_simp = []

    for n in n_values:
        n_even = int(n) if int(n) % 2 == 0 else int(n) + 1
        I_trap = composite_trapezoidal(f, a, b, int(n))
        I_simp = composite_simpson(f, a, b, n_even)
        errors_trap.append(abs(I_trap - exact))
        errors_simp.append(abs(I_simp - exact))

    # 可视化
    plt.figure(figsize=(12, 10))

    # 子图1：函数形态及不同n下的近似
    plt.subplot(2, 2, 1)
    x_fine = np.linspace(a, b, 1000)
    y_fine = f(x_fine)
    plt.plot(x_fine, y_fine, 'k-', label='f(x): 尖锐高斯峰', linewidth=2)
    for n in [5, 10, 20]:
        x_nodes = np.linspace(a, b, n+1)
        y_nodes = f(x_nodes)
        plt.plot(x_nodes, y_nodes, 'o--', label=f'n={n} 采样点', markersize=4, linewidth=1, alpha=0.7)
    plt.xlabel('x')
    plt.ylabel('f(x)')
    plt.title('函数形态与不同密度采样')
    plt.legend()
    plt.grid(True, alpha=0.3)

    # 子图2：误差收敛对比
    plt.subplot(2, 2, 2)
    plt.loglog(n_values, errors_trap, 'o-', label='梯形法误差', linewidth=2)
    plt.loglog(n_values, errors_simp, 's-', label='辛普森法误差', linewidth=2)
    plt.xlabel('区间数量 n')
    plt.ylabel('绝对误差')
    plt.title('误差收敛对比 (对数坐标)')
    plt.legend()
    plt.grid(True, which="both", ls="--", alpha=0.5)

    # 子图3：绝对误差随n变化（线性坐标）
    plt.subplot(2, 2, 3)
    plt.plot(n_values, errors_trap, 'o-', label='梯形法误差')
    plt.plot(n_values, errors_simp, 's-', label='辛普森法误差')
    plt.xlabel('区间数量 n')
    plt.ylabel('绝对误差')
    plt.title('误差下降趋势 (线性坐标)')
    plt.legend()
    plt.grid(True, alpha=0.3)

    # 子图4：达到指定精度所需的区间数对比
    plt.subplot(2, 2, 4)
    target_errors = [1e-2, 1e-3, 1e-4, 1e-5]
    n_required_trap = []
    n_required_simp = []
    for target in target_errors:
        # 简单地从已有数据中查找（实际应用中可用更智能的搜索）
        for i, err in enumerate(errors_trap):
            if err < target:
                n_required_trap.append(n_values[i])
                break
        for i, err in enumerate(errors_simp):
            if err < target:
                n_required_simp.append(n_values[i])
                break
    width = 0.35
    x = np.arange(len(target_errors))
    plt.bar(x - width/2, n_required_trap, width, label='梯形法所需 n', color='orange')
    plt.bar(x + width/2, n_required_simp, width, label='辛普森法所需 n', color='green')
    plt.xlabel('目标误差')
    plt.ylabel('所需最小区间数 n')
    plt.title('达到指定精度所需计算量对比')
    plt.xticks(x, [f'1e-{int(-np.log10(te))}' for te in target_errors])
    plt.legend()
    plt.grid(True, axis='y', alpha=0.3)

    plt.tight_layout()
    plt.show()

compare_methods_on_peak()

这个案例的图表非常具有启发性。左上图显示，当采样点很稀疏（n=5或10）时，采样点很可能完全错过那个尖锐的峰，导致两种方法都严重低估了积分值。只有当区间划分足够细密，能够“捕捉”到峰值时，精度才会快速提升。误差收敛图（右上和左下）显示，在n较小时，误差下降缓慢甚至出现平台；当n大到足以解析峰值后，误差才开始以理论速度下降。最说明问题的是右下角的柱状图：为了达到 1e-4 的精度，梯形法可能需要超过600个区间，而辛普森法可能只需要不到100个。这生动地展示了高阶方法在计算效率上的巨大优势，尤其是在处理复杂函数时。

4. 进阶话题：误差分析与自适应策略

经过前面的实战，我们已经对两种方法的性能有了直观感受。但作为一名严谨的实践者，我们还需要更深入地理解误差的来源，并探索如何让算法更智能。

4.1 误差来源与理论估计

数值积分的误差主要来自两个方面：截断误差和舍入误差。

• 截断误差：源于我们用简单几何图形（梯形或抛物线）代替复杂曲线所带来的近似误差。这正是我们之前讨论的 O(h²) 和 O(h⁴) 的理论误差阶。
• 舍入误差：由于计算机浮点数精度有限，在大量加减乘除运算中累积的微小误差。通常，当步长 h 非常小时，计算节点数量剧增，舍入误差可能会开始显现甚至主导。

对于复合梯形公式，其截断误差有一个实用的估计式： E_T ≈ - (b-a) * h² * f''(ξ) / 12 其中 ξ 是 [a, b] 区间内的某一点。这意味着误差与步长的平方成正比，也与函数二阶导数的平均值有关。如果函数曲率很大（二阶导绝对值大），误差就大。

对于复合辛普森公式，误差估计为： E_S ≈ - (b-a) * h⁴ * f⁽⁴⁾(ξ) / 180 误差与步长的四次方成正比，与函数的四阶导数有关。

注意：这些误差估计式中的 ξ 是未知的，因此它们主要用于理论分析和指导步长选择，不能用于给出精确的误差值。在实际中，我们常采用理查德森外推或运行两种不同步长的计算并比较结果来估计误差。

4.2 实现一个简单的自适应积分器

固定步长的方法有个明显缺点：它在函数平缓区域浪费计算资源，在变化剧烈区域又可能精度不足。自适应积分的思想是“把钱花在刀刃上”，在需要的地方自动加密网格。

这里我们实现一个基于复合辛普森公式的自适应算法（递归版本）的核心思路：

1. 用辛普森公式计算整个区间 [a, b] 的积分近似值 S。
2. 将区间对半分为两个子区间 [a, c] 和 [c, b]，分别计算积分 S_left 和 S_right。
3. 如果 |S - (S_left + S_right)| < tolerance（容差），我们认为当前精度足够，返回 S_left + S_right 作为结果。
4. 否则，说明精度不够，我们递归地对左右两个子区间分别应用同样的过程，并将容差减半（或按区间长度比例分配）。

def adaptive_simpson(f, a, b, tol=1e-6, max_depth=20, depth=0):
    """
    递归自适应辛普森积分法。

    参数
    ----------
    f : function
        被积函数。
    a, b : float
        积分区间。
    tol : float
        目标误差容限。
    max_depth : int
        最大递归深度，防止栈溢出。
    depth : int
        当前递归深度（内部使用）。

    返回
    -------
    float
        积分的近似值。
    """
    if depth > max_depth:
        print(f"警告：达到最大递归深度 {max_depth}。在区间 [{a:.3f}, {b:.3f}] 可能未收敛。")
        return composite_simpson(f, a, b, 2)  # 退回简单辛普森

    c = (a + b) / 2
    S_ab = composite_simpson(f, a, b, 2)          # 整个区间用1个抛物线（2个子区间）
    S_ac = composite_simpson(f, a, c, 2)          # 左半区间
    S_cb = composite_simpson(f, c, b, 2)          # 右半区间
    S_total = S_ac + S_cb

    # 误差估计：|S_ab - S_total| / 15 是一个常用的误差估计量
    error_est = abs(S_ab - S_total) / 15

    if error_est < tol * (b - a) / (1.0 if depth == 0 else 1.0): # 根节点用绝对容差，子节点用相对
        return S_total
    else:
        # 递归细化
        left_result = adaptive_simpson(f, a, c, tol/2, max_depth, depth+1)
        right_result = adaptive_simpson(f, c, b, tol/2, max_depth, depth+1)
        return left_result + right_result

# 测试自适应积分器
def test_adaptive():
    """测试自适应积分器在尖锐峰值函数上的表现"""
    f = lambda x: np.exp(-50 * (x - 0.5)**2) + 0.1
    a, b = 0, 1
    exact, _ = integrate.quad(f, a, b)

    tol = 1e-5
    result_adaptive = adaptive_simpson(f, a, b, tol=tol)
    result_fixed_n100 = composite_simpson(f, a, b, 100)
    result_fixed_n500 = composite_simpson(f, a, b, 500)

    print("自适应辛普森积分测试 (函数: exp(-50*(x-0.5)^2)+0.1)")
    print("=" * 60)
    print(f"参考精确值 (SciPy): {exact:.12f}")
    print(f"自适应结果 (tol={tol}): {result_adaptive:.12f}, 误差: {abs(result_adaptive-exact):.2e}")
    print(f"固定 n=100 结果: {result_fixed_n100:.12f}, 误差: {abs(result_fixed_n100-exact):.2e}")
    print(f"固定 n=500 结果: {result_fixed_n500:.12f}, 误差: {abs(result_fixed_n500-exact):.2e}")

    # 可视化自适应过程（记录递归调用的区间）
    intervals_visited = []
    def record_adaptive(f, a, b, tol=1e-6, max_depth=20, depth=0):
        intervals_visited.append((a, b, depth))
        if depth > max_depth:
            return composite_simpson(f, a, b, 2)
        c = (a + b) / 2
        S_ab = composite_simpson(f, a, b, 2)
        S_ac = composite_simpson(f, a, c, 2)
        S_cb = composite_simpson(f, c, b, 2)
        S_total = S_ac + S_cb
        error_est = abs(S_ab - S_total) / 15
        if error_est < tol * (b - a):
            return S_total
        else:
            left = record_adaptive(f, a, c, tol/2, max_depth, depth+1)
            right = record_adaptive(f, c, b, tol/2, max_depth, depth+1)
            return left + right

    _ = record_adaptive(f, a, b, tol=1e-5)
    # 绘制函数和自适应区间
    plt.figure(figsize=(10, 6))
    x_fine = np.linspace(a, b, 1000)
    y_fine = f(x_fine)
    plt.plot(x_fine, y_fine, 'k-', linewidth=2, label='f(x)')
    colors = plt.cm.viridis(np.linspace(0, 1, max([iv[2] for iv in intervals_visited])+1))
    for (a_i, b_i, d) in intervals_visited:
        plt.axvspan(a_i, b_i, alpha=0.2, color=colors[d], edgecolor=None)
    plt.xlabel('x')
    plt.ylabel('f(x)')
    plt.title(f'自适应辛普森积分区间细化 (容差={tol})')
    # 创建深度图例
    from matplotlib.patches import Patch
    legend_elements = [Patch(facecolor=colors[i], alpha=0.5, label=f'递归深度 {i}') for i in range(0, 5, 2)]
    plt.legend(handles=legend_elements + [plt.Line2D([0], [0], color='k', linewidth=2, label='f(x)')])
    plt.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.show()

test_adaptive()

运行这个测试，你会看到自适应积分器如何智能地将计算资源集中在函数变化剧烈的峰值区域（x=0.5附近），而在平坦区域使用较粗的网格。它用比固定步长 n=500 少得多的函数求值次数，达到了相近甚至更高的精度。这在实际计算中，尤其是对复杂函数进行高精度积分时，能显著节省计算时间。

5. 性能基准测试与选择指南

在文章的最后，我们从一个更工程化的角度出发，系统地对比一下两种方法的计算性能，并给出一些实用的选择建议。

我们设计一个简单的基准测试，统计计算不同函数、不同精度要求下，两种方法所需的函数求值次数和计算时间。函数求值次数是衡量数值积分算法成本的一个关键指标，因为在实际应用中，函数 f(x) 本身可能非常复杂和耗时。

import time

def benchmark_integration(f, a, b, exact, method, n):
    """基准测试单个方法和区间数"""
    start_time = time.perf_counter()
    if method == 'trap':
        result = composite_trapezoidal(f, a, b, n)
        evals = n + 1
    elif method == 'simp':
        n = n if n % 2 == 0 else n + 1
        result = composite_simpson(f, a, b, n)
        evals = n + 1
    else:
        raise ValueError("方法必须是 'trap' 或 'simp'")
    end_time = time.perf_counter()
    error = abs(result - exact)
    elapsed = end_time - start_time
    return error, evals, elapsed

def run_comprehensive_benchmark():
    """运行综合性能基准测试"""
    # 定义三个测试函数及其精确积分值（在区间[0,1]上）
    test_cases = [
        ("平滑多项式", lambda x: x**4 + 2*x**2 + 1, 1.0 + 2.0/3.0 + 1.0/5.0),
        ("振荡函数", lambda x: np.sin(10*x) + np.cos(5*x), (1/5)*np.sin(5) - (1/10)*(np.cos(10)-1)),
        ("尖锐峰值", lambda x: np.exp(-50*(x-0.5)**2) + 0.1, integrate.quad(lambda x: np.exp(-50*(x-0.5)**2)+0.1, 0, 1)[0])
    ]

    target_errors = [1e-4, 1e-6, 1e-8]
    methods = ['trap', 'simp']
    results = {case[0]: {method: {str(target): {'n': None, 'evals': None, 'time': None, 'error': None} for target in target_errors} for method in methods} for case in test_cases}

    # 为每个测试用例和每个目标误差，寻找满足条件的最小n
    for case_name, func, exact in test_cases:
        print(f"\n{'='*60}")
        print(f"测试函数: {case_name}")
        print(f"{'='*60}")
        for target in target_errors:
            for method in methods:
                n = 4 if method == 'simp' else 2  # 起始n
                found = False
                while n <= 100000:  # 设置上限防止无限循环
                    error, evals, elapsed = benchmark_integration(func, 0, 1, exact, method, n)
                    if error < target:
                        results[case_name][method][str(target)]['n'] = n
                        results[case_name][method][str(target)]['evals'] = evals
                        results[case_name][method][str(target)]['time'] = elapsed
                        results[case_name][method][str(target)]['error'] = error
                        found = True
                        break
                    n *= 2  # 指数增加n以快速搜索
                if not found:
                    results[case_name][method][str(target)]['n'] = '>100k'

    # 打印结果摘要表格
    for case_name in results:
        print(f"\n**{case_name} 性能对比摘要**")
        print("-" * 80)
        header = f"{'目标误差':<12} | {'方法':<10} | {'所需 n':<8} | {'函数求值次数':<12} | {'计算时间(秒)':<15} | {'实际误差':<12}"
        print(header)
        print("-" * 80)
        for target in target_errors:
            for method in methods:
                res = results[case_name][method][str(target)]
                method_name = '梯形法' if method == 'trap' else '辛普森法'
                n_str = str(res['n']) if res['n'] is not None else 'N/A'
                evals_str = str(res['evals']) if res['evals'] is not None else 'N/A'
                time_str = f"{res['time']:.6f}" if res['time'] is not None else 'N/A'
                error_str = f"{res['error']:.2e}" if res['error'] is not None else 'N/A'
                print(f"{target:<12.0e} | {method_name:<10} | {n_str:<8} | {evals_str:<12} | {time_str:<15} | {error_str:<12}")

    # 可视化：为第一个测试用例绘制性能对比图
    case_name = "平滑多项式"
    func, exact = test_cases[0][1], test_cases[0][2]
    n_range = np.logspace(1, 3, 20, dtype=int)  # 从10到1000
    times_trap, errors_trap = [], []
    times_simp, errors_simp = [], []

    for n in n_range:
        # 梯形法
        err, evals, t = benchmark_integration(func, 0, 1, exact, 'trap', n)
        errors_trap.append(err)
        times_trap.append(t)
        # 辛普森法 (确保n为偶数)
        n_simp = n if n % 2 == 0 else n + 1
        err, evals, t = benchmark_integration(func, 0, 1, exact, 'simp', n_simp)
        errors_simp.append(err)
        times_simp.append(t)

    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
    # 左图：误差 vs 计算时间
    ax = axes[0]
    ax.loglog(times_trap, errors_trap, 'o-', label='复合梯形公式', linewidth=2)
    ax.loglog(times_simp, errors_simp, 's-', label='复合辛普森公式', linewidth=2)
    ax.set_xlabel('计算时间 (秒, 对数坐标)')
    ax.set_ylabel('绝对误差 (对数坐标)')
    ax.set_title('计算时间 vs 精度 (平滑多项式)')
    ax.legend()
    ax.grid(True, which="both", ls="--", alpha=0.5)

    # 右图：误差 vs 函数求值次数
    ax = axes[1]
    # 注意：这里为了简化，用n近似代表求值次数（对于梯形法是n+1，辛普森法也是n+1）
    ax.loglog(n_range, errors_trap, 'o-', label='复合梯形公式', linewidth=2)
    ax.loglog([n if n%2==0 else n+1 for n in n_range], errors_simp, 's-', label='复合辛普森公式', linewidth=2)
    ax.set_xlabel('区间数量 n (近似求值次数, 对数坐标)')
    ax.set_ylabel('绝对误差 (对数坐标)')
    ax.set_title('计算成本 vs 精度 (平滑多项式)')
    ax.legend()
    ax.grid(True, which="both", ls="--", alpha=0.5)

    plt.tight_layout()
    plt.show()

run_comprehensive_benchmark()

基准测试的结果表格和图表会清晰地告诉你：在大多数情况下，尤其是对光滑函数和高精度要求，复合辛普森公式在计算效率上具有压倒性优势。为了达到相同的精度（例如 1e-8），辛普森法所需的函数求值次数和计算时间通常比梯形法少一个数量级甚至更多。

然而，这并不意味着辛普森法是万能的。根据我的经验，在以下场景中，你可能需要重新考虑选择：

• 函数不光滑或存在奇点：如果函数有间断点、尖点或导数不存在，高阶方法（如辛普森法）的理论误差阶可能不成立，其表现可能反而不如更稳健的梯形法，甚至可能发散。这时，可能需要先将积分区间在奇点处拆分，或者考虑使用专门处理奇异积分的方法。
• 超低精度快速估算：如果你只需要一个数量级正确的粗略估计，梯形法实现简单，几行代码就能搞定，心理负担小。
• 函数求值极其昂贵：虽然辛普森法通常需要更少的求值次数，但它的逻辑稍复杂。如果每次函数求值都需要调用一个耗时数秒的仿真程序，那么两种方法的总时间差异可能被掩盖，此时代码的简洁性和可维护性可能成为更重要的考量。
• 自适应积分：如我们前面实现的，自适应策略可以同时发挥两种方法的优点。通常，自适应算法底层会使用一个像辛普森法这样的高阶方法进行局部积分和误差估计。

所以，我的建议是：对于一般的平滑函数，优先使用复合辛普森公式。在实现时，记得检查区间数 n 是否为偶数。如果对函数的性态不确定，或者需要处理可能存在问题的积分，可以先尝试梯形法作为基准，再逐步切换到更高效的方法。对于复杂的实际问题，利用 SciPy 库中成熟的积分函数（如 scipy.integrate.quad）通常是更稳妥和高效的选择，它们内部已经集成了自适应、误差控制等多种高级策略。我们自己动手实现这些经典方法，最大的价值在于理解其原理，从而能在关键时刻做出正确的判断和选择。