Numerical Integration

import numpy as np
from typing import Callable, Tuple


class NumericalIntegrator:
    """Comprehensive numerical integration toolkit."""

    @staticmethod
    def trapezoidal(f: Callable, a: float, b: float,
                    n: int = 1000) -> Tuple[float, float]:
        """Composite trapezoidal rule. Returns (result, error_estimate)."""
        x = np.linspace(a, b, n + 1)
        y = f(x)
        h = (b - a) / n
        result = h * (y[0] / 2 + np.sum(y[1:-1]) + y[-1] / 2)

        # Error estimate via Richardson extrapolation
        x2 = np.linspace(a, b, 2 * n + 1)
        y2 = f(x2)
        h2 = h / 2
        result2 = h2 * (y2[0] / 2 + np.sum(y2[1:-1]) + y2[-1] / 2)

        error = abs(result2 - result) / 3  # Richardson estimate
        return result, error

    @staticmethod
    def simpson(f: Callable, a: float, b: float,
                n: int = 1000) -> Tuple[float, float]:
        """Composite Simpson's rule (n forced even). Returns (result, error)."""
        if n % 2 != 0:
            n += 1
        x = np.linspace(a, b, n + 1)
        y = f(x)
        h = (b - a) / n
        result = h / 3 * (y[0] + 4 * np.sum(y[1::2]) +
                          2 * np.sum(y[2:-1:2]) + y[-1])

        # Error estimate
        n2 = 2 * n
        x2 = np.linspace(a, b, n2 + 1)
        y2 = f(x2)
        h2 = (b - a) / n2
        result2 = h2 / 3 * (y2[0] + 4 * np.sum(y2[1::2]) +
                            2 * np.sum(y2[2:-1:2]) + y2[-1])

        error = abs(result2 - result) / 15  # Richardson estimate
        return result, error

    @staticmethod
    def gauss_legendre(f: Callable, a: float, b: float,
                       n: int = 5) -> float:
        """n-point Gauss-Legendre quadrature."""
        nodes, weights = np.polynomial.legendre.leggauss(n)
        t = 0.5 * (b - a) * nodes + 0.5 * (b + a)
        w = 0.5 * (b - a) * weights
        return np.sum(w * f(t))

    @staticmethod
    def romberg(f: Callable, a: float, b: float,
                max_level: int = 10) -> np.ndarray:
        """Romberg integration using Richardson extrapolation."""
        R = np.zeros((max_level + 1, max_level + 1))

        # Trapezoidal estimates at various levels
        n = 1
        h = b - a
        R[0, 0] = h / 2 * (f(a) + f(b))

        for i in range(1, max_level + 1):
            n = 2 ** i
            h = (b - a) / n
            # Add new midpoints
            x_new = np.arange(1, n, 2) * h + a
            R[i, 0] = 0.5 * R[i-1, 0] + h * np.sum(f(x_new))

            # Richardson extrapolation
            for j in range(1, i + 1):
                R[i, j] = R[i, j-1] + (R[i, j-1] - R[i-1, j-1]) / (4**j - 1)

        return R

    @staticmethod
    def monte_carlo(f: Callable, domain: list,
                    n_samples: int = 100000) -> Tuple[float, float]:
        """Monte Carlo integration over a rectangular domain."""
        a, b = domain[0], domain[1]
        volume = b - a

        x = np.random.uniform(a, b, n_samples)
        values = f(x)
        result = volume * np.mean(values)
        std_err = volume * np.std(values) / np.sqrt(n_samples)

        return result, std_err


def demonstrate_integration():
    """Full demonstration of all integration methods."""
    print("=" * 60)
    print("NUMERICAL INTEGRATION METHOD COMPARISON")
    print("=" * 60)

    # Test 1: Smooth function
    f1 = lambda x: np.exp(-x**2)
    exact1 = np.sqrt(np.pi)  # ∫₋∞^∞ e^(-x²) dx = √π
    # But we integrate from -5 to 5 (close enough)

    from scipy.integrate import quad
    exact1_quad, _ = quad(f1, -5, 5)

    print(f"\n∫₋₅⁵ e^(-x²) dx")
    print(f"scipy quad: {exact1_quad:.12f}")

    integrator = NumericalIntegrator()

    for n in [10, 50, 200]:
        trap, err = integrator.trapezoidal(f1, -5, 5, n)
        simp, err = integrator.simpson(f1, -5, 5, n)
        print(f"  n={n:>3}: Trapezoidal {trap:.10f} (err {err:.2e}) | "
              f"Simpson {simp:.10f} (err {err:.2e})")

    # Gauss-Legendre
    for n in [3, 5, 10]:
        gl = integrator.gauss_legendre(f1, -5, 5, n)
        print(f"  Gauss-Legendre (n={n}): {gl:.10f}")

    # Romberg
    R = integrator.romberg(f1, -5, 5, max_level=6)
    print(f"\nRomberg integration:")
    for i in range(min(7, R.shape[0])):
        print(f"  Level {i}: {R[i, i]:.12f}")

    # Test 2: Oscillatory function
    print(f"\n{'=' * 60}")
    f2 = lambda x: np.sin(100 * x) / (1 + x**2)
    mc_result, mc_err = integrator.monte_carlo(f2, [-10, 10], 1_000_000)
    quad_result, _ = quad(f2, -10, 10)
    print(f"∫₋₁₀¹⁰ sin(100x)/(1+x²) dx")
    print(f"scipy quad:  {quad_result:.10f}")
    print(f"Monte Carlo: {mc_result:.10f} (±{mc_err:.2e})")


demonstrate_integration()