adaptiveSimpson function

double adaptiveSimpson(double f(double), double a, double b, { double tol = 1e-12, int maxRecursion = 20, })

Adaptive Simpson integrator.

High-accuracy adaptive Simpson that recursively subdivides intervals until the requested tolerance is reached. Suitable for smooth integrands and robust enough for production workloads.

Implementation

double adaptiveSimpson(
  double Function(double) f,
  double a,
  double b, {
  double tol = 1e-12,
  int maxRecursion = 20,
}) {
  double simpson(double fa, double fm, double fb, double h) =>
      (fa + 4 * fm + fb) * h / 6.0;

  double recurse(
    double a,
    double b,
    double fa,
    double fb,
    double fm,
    double whole,
    int depth,
  ) {
    final m = 0.5 * (a + b);
    final lm = 0.5 * (a + m);
    final rm = 0.5 * (m + b);
    final flm = f(lm);
    final frm = f(rm);
    final left = simpson(fa, flm, fm, (m - a));
    final right = simpson(fm, frm, fb, (b - m));
    final delta = left + right - whole;
    if (delta.abs() <= 15 * tol || depth >= maxRecursion) {
      return left + right + delta / 15.0;
    }
    return recurse(a, m, fa, fm, flm, left, depth + 1) +
        recurse(m, b, fm, fb, frm, right, depth + 1);
  }

  final fa = f(a);
  final fb = f(b);
  final m = 0.5 * (a + b);
  final fm = f(m);
  final whole = (fa + 4 * fm + fb) * (b - a) / 6.0;
  return recurse(a, b, fa, fb, fm, whole, 0);
}