pip install numpy matplotlib scipy

import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.animation import FuncAnimationfrom matplotlib.patches import Polygonfrom scipy import integrate# ==============================================================================# 1. 定义目标函数# ==============================================================================def f(x):    """定义我们要积分的函数: f(x) = (π * cos(e^x)) / x^e"""    return np.pi * np.cos(np.exp(x)) / (x**np.e)# ==============================================================================# 2. 核心算法实现# ==============================================================================def simpson_rule(func, a, b):    """基础的辛普森法，用于计算单个区间[a, b]的积分"""    c = (a + b) / 2.0    h = (b - a) / 6.0    return h * (func(a) + 4.0 * func(c) + func(b))def adaptive_simpson_with_logging(func, a, b, tol, level, log):    """    带日志记录的自适应辛普森递归函数。
    参数:    func: 目标函数    a, b: 当前积分区间    tol: 当前区间的误差容忍度    level: 当前递归深度 (用于动画着色)    log: 一个列表，用于记录每一步的操作 (a, b, level, is_accepted)    """    c = (a + b) / 2.0    whole = simpson_rule(func, a, b)    left = simpson_rule(func, a, c)    right = simpson_rule(func, c, b)
    # 误差估计    error_estimate = abs(left + right - whole)
    # 误差阈值: 15 * tol 是一个经验系数，确保总误差在tol范围内    if error_estimate < 15.0 * tol:        # 如果误差足够小，接受该区间的积分值        log.append({'a': a, 'b': b, 'level': level, 'accepted': True})        return left + right    else:        # 否则，记录当前区间被拒绝，然后递归处理子区间        log.append({'a': a, 'b': b, 'level': level, 'accepted': False})        # 递归调用，并将误差容忍度减半        return (adaptive_simpson_with_logging(func, a, c, tol / 2, level + 1, log) +                adaptive_simpson_with_logging(func, c, b, tol / 2, level + 1, log))# ==============================================================================# 3. 动画生成逻辑# ==============================================================================def create_animation(log_data, integral_result, a, b):    """使用matplotlib创建并显示动画"""    fig, ax = plt.subplots(figsize=(12, 7))    # 配置中文字体    plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'WenQuanYi Micro Hei']    plt.rcParams['axes.unicode_minus'] = False
    # 绘制完整的函数曲线作为背景    x_full = np.linspace(a, b, 500)    y_full = f(x_full)    ax.plot(x_full, y_full, 'k-', lw=2, label=r'$f(x) = \frac{\pi \cos(e^x)}{x^e}$')
    # 设置图形标题和坐标轴    ax.set_title("自适应辛普森积分动态演示", fontsize=16)    ax.set_xlabel("x", fontsize=12)    ax.set_ylabel("f(x)", fontsize=12)    ax.set_xlim(a - 0.1, b + 0.1)    # 动态调整y轴范围以更好地显示函数    y_min, y_max = np.min(y_full), np.max(y_full)    ax.set_ylim(y_min - 0.5, y_max + 0.5)    ax.grid(True, linestyle='--', alpha=0.6)    ax.legend(loc='upper right')    # 用于存储已接受区间的多边形    accepted_patches = []
    # 用于显示当前信息的文本    info_text = ax.text(0.02, 0.95, '', transform=ax.transAxes, fontsize=11,                        verticalalignment='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))    def init():        """动画初始化函数"""        return []    def animate(frame):        """动画更新函数，每一帧对应log_data中的一个条目"""        # 清除上一帧的临时图形（不包括函数曲线和网格）        for p in accepted_patches:            p.remove()        accepted_patches.clear()        current_step_data = log_data[frame]
        # 遍历到当前帧为止的所有数据，绘制所有已接受的区间        current_integral = 0        for i in range(frame + 1):            step = log_data[i]            if step['accepted']:                verts = [                    (step['a'], 0),                     (step['a'], f(step['a'])),                     (step['b'], f(step['b'])),                     (step['b'], 0)                ]                # 用绿色半透明多边形表示已接受的区间                poly = Polygon(verts, facecolor='0.8', edgecolor='green', alpha=0.5, lw=1.5)                ax.add_patch(poly)                accepted_patches.append(poly)                # 累加已接受区间的积分值                current_integral += simpson_rule(f, step['a'], step['b'])
        # 高亮当前正在处理的区间        if not current_step_data['accepted']:            verts = [                (current_step_data['a'], 0),                 (current_step_data['a'], f(current_step_data['a'])),                 (current_step_data['b'], f(current_step_data['b'])),                 (current_step_data['b'], 0)            ]            # 用红色半透明多边形表示正在被细分的区间            poly = Polygon(verts, facecolor='red', alpha=0.4, lw=1.5)            ax.add_patch(poly)            accepted_patches.append(poly) # 临时加入，以便在下一帧清除        # 更新信息文本        status = "已接受" if current_step_data['accepted'] else "正在细分"        info_text.set_text(            f"步骤: {frame + 1} / {len(log_data)}\n"            f"当前区间: [{current_step_data['a']:.3f}, {current_step_data['b']:.3f}]\n"            f"递归深度: {current_step_data['level']}\n"            f"状态: {status}\n"            f"当前积分值: {current_integral:.10f}"        )
        return accepted_patches + [info_text]    # 创建动画    ani = FuncAnimation(fig, animate, frames=len(log_data),                        init_func=init, blit=True, interval=100, repeat=False)    return ani# ==============================================================================# 4. 主程序# ==============================================================================if __name__ == "__main__":    # 定义积分参数    INTEGRAL_LOWER_BOUND = 1    INTEGRAL_UPPER_BOUND = 4    TOLERANCE = 1e-7  # 总体误差容忍度    print(f"正在计算函数 f(x) 在区间 [{INTEGRAL_LOWER_BOUND}, {INTEGRAL_UPPER_BOUND}] 上的积分...")    print(f"设定总体误差容忍度为: {TOLERANCE}")
    # 创建一个空列表来记录算法的每一步    algorithm_log = []
    # 运行自适应辛普森积分并记录过程    integral_value = adaptive_simpson_with_logging(        f,         INTEGRAL_LOWER_BOUND,         INTEGRAL_UPPER_BOUND,         TOLERANCE,         level=0,         log=algorithm_log    )
    print("-" * 50)    print(f"自适应辛普森法计算结果: {integral_value:.10f}")    # 使用SciPy的quad函数作为精度基准    scipy_result, scipy_error = integrate.quad(f, INTEGRAL_LOWER_BOUND, INTEGRAL_UPPER_BOUND)    print(f"SciPy quad函数计算结果: {scipy_result:.10f} (估计误差: {scipy_error:.2e})")    print("-" * 50)    print("动画演示即将开始...")    # 创建并运行动画    animation = create_animation(algorithm_log, integral_value, INTEGRAL_LOWER_BOUND, INTEGRAL_UPPER_BOUND)    plt.show()