每日算法 — 使用java实现辛普森法则：自适应数值积分与误差控制

数值积分是科学计算的核心工具，当函数没有解析原函数或仅有点采样数据时，数值方法成为求定积分的唯一途径。辛普森法则（Simpson’s Rule）通过用抛物线拟合局部区间，将梯形法则的误差从 $O(h^2)$ 降至 $O(h^4)$，在工程与物理仿真中广泛应用。

算法原理

基本辛普森公式

在区间 $[a, b]$ 上，用经过三点 $(a, f(a))$、$((a+b)/2, f((a+b)/2))$、$(b, f(b))$ 的抛物线近似函数曲线：

$$
\int_a^b f(x)\,dx \approx \frac{b-a}{6}\left[f(a) + 4f\left(\frac{a+b}{2}\right) + f(b)\right]
$$

复合辛普森法则

将 $[a, b]$ 等分为 $n$（偶数）个子区间，步长 $h=(b-a)/n$，则：

$$
S_n = \frac{h}{3}\left[f(x_0) + 4\sum_{i=1,3,5…}^{n-1}f(x_i) + 2\sum_{i=2,4,6…}^{n-2}f(x_i) + f(x_n)\right]
$$

自适应辛普森

通过递归比较 $S(a,b)$ 与 $S(a,m)+S(m,b)$ 的差值，仅在误差较大区域细分，兼顾精度与效率。

Java实现

import java.util.function.DoubleFunction;

/**
 * 辛普森数值积分器
 * 支持基本辛普森、复合辛普森与自适应辛普森三种模式
 */
public class SimpsonIntegrator {

    /**
     * 基本辛普森公式：单区间三点抛物线拟合
     * @param f 被积函数
     * @param a 左端点
     * @param b 右端点
     * @return 积分近似值
     */
    public static double basicSimpson(DoubleFunction<Double> f, double a, double b) {
        double mid = (a + b) / 2.0;
        return (b - a) / 6.0 * (f.apply(a) + 4 * f.apply(mid) + f.apply(b));
    }

    /**
     * 复合辛普森法则：等距细分求和
     * @param f 被积函数
     * @param a 左端点
     * @param b 右端点
     * @param n 子区间数（必须为偶数）
     * @return 积分近似值
     */
    public static double compositeSimpson(DoubleFunction<Double> f, double a, double b, int n) {
        if (n % 2 != 0) {
            throw new IllegalArgumentException("n必须为偶数");
        }
        double h = (b - a) / n;
        double sum = f.apply(a) + f.apply(b);
        // 奇数点系数为4
        for (int i = 1; i < n; i += 2) {
            sum += 4 * f.apply(a + i * h);
        }
        // 偶数点系数为2
        for (int i = 2; i < n; i += 2) {
            sum += 2 * f.apply(a + i * h);
        }
        return sum * h / 3.0;
    }

    /**
     * 自适应辛普森递归实现
     * @param f 被积函数
     * @param a 左端点
     * @param b 右端点
     * @param eps 误差容限
     * @param maxDepth 最大递归深度
     * @return 满足精度要求的积分值
     */
    public static double adaptiveSimpson(
            DoubleFunction<Double> f, double a, double b, double eps, int maxDepth) {
        return adaptiveSimpsonRecursive(f, a, b, eps, maxDepth,
                f.apply(a), f.apply((a + b) / 2.0), f.apply(b), basicSimpson(f, a, b));
    }

    private static double adaptiveSimpsonRecursive(
            DoubleFunction<Double> f, double a, double b, double eps, int depth,
            double fa, double fm, double fb, double whole) {
        double m = (a + b) / 2.0;
        double lm = (a + m) / 2.0;
        double rm = (m + b) / 2.0;
        double flm = f.apply(lm);
        double frm = f.apply(rm);
        double left = (m - a) / 6.0 * (fa + 4 * flm + fm);
        double right = (b - m) / 6.0 * (fm + 4 * frm + fb);
        double delta = left + right - whole;

        if (depth <= 0 || Math.abs(delta) <= 15 * eps) {
            return left + right + delta / 15.0; // Richardson外推修正
        }
        return adaptiveSimpsonRecursive(f, a, m, eps / 2, depth - 1, fa, flm, fm, left)
             + adaptiveSimpsonRecursive(f, m, b, eps / 2, depth - 1, fm, frm, fb, right);
    }

    public static void main(String[] args) {
        // 测试1: ∫₀^π sin(x) dx = 2
        DoubleFunction<Double> sin = Math::sin;
        System.out.println("=== 测试 ∫₀^π sin(x) dx = 2 ===");
        System.out.printf("基本辛普森: %.10f%n", basicSimpson(sin, 0, Math.PI));
        System.out.printf("复合辛普森(n=100): %.10f%n", compositeSimpson(sin, 0, Math.PI, 100));
        System.out.printf("自适应辛普森: %.10f%n", adaptiveSimpson(sin, 0, Math.PI, 1e-10, 20));

        // 测试2: ∫₀¹ e^x dx = e - 1 ≈ 1.718281828...
        DoubleFunction<Double> exp = Math::exp;
        System.out.println("\n=== 测试 ∫₀¹ e^x dx = e-1 ===");
        System.out.printf("基本辛普森: %.10f%n", basicSimpson(exp, 0, 1));
        System.out.printf("复合辛普森(n=100): %.10f%n", compositeSimpson(exp, 0, 1, 100));
        System.out.printf("自适应辛普森: %.10f%n", adaptiveSimpson(exp, 0, 1, 1e-12, 20));

        // 测试3: ∫₀¹ 4/(1+x²) dx = π (蒙特卡洛经典对比)
        DoubleFunction<Double> piFunc = x -> 4.0 / (1.0 + x * x);
        System.out.println("\n=== 测试 ∫₀¹ 4/(1+x²) dx = π ===");
        System.out.printf("基本辛普森: %.10f%n", basicSimpson(piFunc, 0, 1));
        System.out.printf("复合辛普森(n=1000): %.10f%n", compositeSimpson(piFunc, 0, 1, 1000));
        System.out.printf("自适应辛普森: %.10f%n", adaptiveSimpson(piFunc, 0, 1, 1e-12, 20));
    }
}

误差分析与复杂度

自适应辛普森通过递归细分将计算资源集中在函数曲率大的区域，对于振荡函数（如 $sin(1/x)$ 在 $0$ 附近）可节省数十倍的函数求值次数。

扩展思考

1. 龙贝格积分：在复合辛普森基础上进行Richardson外推，进一步将误差降至 $O(h^{2k})$。
2. 高斯-勒让德积分：通过选取最优节点（Legendre多项式零点），$n$ 点公式可达到 $2n-1$ 次代数精度。
3. 反常积分处理：对于无穷区间或瑕积分，可通过变量替换（如 $x=t/(1-t)$）将有界区间映射到 $[0,1]$。

辛普森法则以其简洁性与高效性，至今仍是工程计算和物理仿真的首选积分工具之一。