Numerical Integration Libraries: Domain Explainer#

Purpose: Help technical decision makers understand Python numerical integration libraries and choose the right tool for computing definite integrals in scientific, engineering, and financial applications.

Audience: Software engineers, data scientists, researchers without deep numerical analysis expertise.

What This Solves#

The Problem#

Imagine you’re building a Monte Carlo option pricing model. You need to compute:

∫₀^∞ e^(-x²/2) * payoff(x) dx

Manual approach: Implement trapezoid rule, tune step size, check convergence. Time: Hours to days debugging precision/performance issues.

Library approach: scipy.integrate.quad(lambda x: exp(-x**2/2) * payoff(x), 0, inf) Time: 3 minutes (1 line of code + validation).

This is the numerical integration problem: Computing definite integrals when analytic solutions don’t exist or are impractical.

Who Encounters This#

• Scientific Computing: Physics simulations, computational chemistry, climate models
• Financial Engineering: Option pricing, risk models, VaR calculations
• Data Science: Statistical inference, Bayesian posteriors, expectation computation
• Engineering: Structural analysis, control systems, signal processing

Why It Matters#

Business Impact:

• Speed: Seconds vs hours for integral computation (10,000x speedup typical)
• Accuracy: Adaptive error control (1e-6 to 1e-12 typical, arbitrary precision available)
• Reliability: Battle-tested algorithms (QUADPACK from 1970s-1980s still standard)

Technical Impact:

• High-dimensional integration (6-50+ dimensions via Monte Carlo)
• Infinite intervals, singularities, oscillatory integrands handled automatically
• Multi-platform (NumPy/SciPy runs everywhere Python runs)

Accessible Analogies#

What Is Numerical Integration?#

Analogy: Measuring the Area Under a Curve

Imagine tracing a mountain profile on graph paper:

• Symbolic integration: Find the exact formula for the area (calculus class)
• Numerical integration: Count grid squares under the curve (approximation)

Problem: Most real curves don’t have simple formulas. You MUST approximate.

Numerical integration is like having a smart assistant that:

1. Samples the curve at strategic points (not uniform grid)
2. Fits smooth shapes between points (not just rectangles)
3. Estimates error and refines automatically (adaptive)
4. Handles edge cases (infinite curves, sharp corners)

Why Libraries Matter#

Without libraries: Implement adaptive quadrature yourself

• Research algorithms (Gauss-Kronrod, Clenshaw-Curtis, etc.)
• Handle numerical stability (rounding errors, overflow)
• Debug convergence failures
• Optimize performance

Time: Weeks to months (even for experts)

With libraries: quad(f, a, b) in SciPy

• Pre-built algorithms (QUADPACK Fortran, decades of optimization)
• Automatic error estimation
• Production-tested edge cases
• 1 line of code

Time: Minutes

When You Need This#

Clear Decision Criteria#

You NEED numerical integration if: ✅ Computing expectations, probabilities, statistical posteriors ✅ Solving differential equations (via integral equations) ✅ Volume/area calculations in computational geometry ✅ Transform computations (Fourier, Laplace, etc.) ✅ Physics simulations (partition functions, path integrals)

You DON’T need this if: ❌ Symbolic solution exists and is tractable (use SymPy) ❌ Data is already discrete (use NumPy.sum) ❌ Simple geometric shapes (use formulas: πr², 4/3πr³, etc.)

Concrete Use Case Examples#

Monte Carlo Option Pricing:

• Problem: Compute expected payoff under Black-Scholes SDE
• Solution: scipy.integrate.quad for 1D integrals, vegas for multi-asset (high-D)
• ROI: Minutes to implement vs hours debugging trapezoid rule

Bayesian Posterior Normalization:

• Problem: Compute normalizing constant ∫ p(x|data) dx
• Solution: scipy.integrate.quad (1-3D) or vegas (>3D parameters)
• ROI: Accurate posteriors for inference

Computational Geometry Volume:

• Problem: Volume of intersection between 3D meshes
• Solution: quadpy with tetrahedron schemes (FEM-style)
• ROI: Robust volume calculation for CAD/graphics

Trade-offs#

What You’re Choosing Between#

1. Adaptive vs Fixed-Order Quadrature#

Adaptive (scipy.quad):

• Pros: Automatic error control, handles singularities, reliable
• Cons: Variable runtime, harder to parallelize

When: Default choice (95% of use cases)

Fixed-Order (quadpy Gauss schemes):

• Pros: Predictable cost, easy to parallelize, ultra-high-order available
• Cons: Manual error estimation, may over/under-sample

When: FEM/FVM (fixed integration on standard elements), GPU acceleration

2. Deterministic vs Monte Carlo#

Deterministic (scipy.quad, quadpy):

• Pros: Reproducible, converges faster for low dimensions (1-6D)
• Cons: Exponential cost scaling with dimension (curse of dimensionality)

When: 1-6 dimensions, need reproducibility

Monte Carlo (vegas):

• Pros: Dimension-independent scaling, handles 20-50+ dimensions
• Cons: Stochastic (different results each run), slower convergence (O(1/√N))

When: >6 dimensions (physics, finance with many assets)

3. Standard vs Arbitrary Precision#

Standard Precision (scipy, quadpy with float64):

• Pros: Fast (hardware floating point), sufficient for engineering (1e-12 error)
• Cons: Limited to ~15 decimal digits

When: 99% of applications (engineering, science, finance)

Arbitrary Precision (mpmath.quad):

• Pros: 1-10,000 digits, verification of symbolic results
• Cons: 10-100x slower (pure Python), complex numbers less efficient

When: Theory/research, numerical verification, catastrophic cancellation cases

Cost Considerations#

Why Cost Matters Here#

Numerical integration libraries are free (BSD/MIT licenses). Cost is developer time and compute time.

Pricing Models#

Libraries: $0 (open-source: scipy, quadpy, mpmath, vegas)

Developer Time:

• Setup: 1-4 hours (install scipy, learn API)
• Integration: 1-8 hours (adapt to specific use case)
• Maintenance: 2-10 hours/year (update dependencies)

Compute Time:

• 1D integrals: Microseconds to milliseconds (scipy.quad)
• Multi-D integrals: Seconds to minutes (adaptive subdivision)
• High-D Monte Carlo: Minutes to hours (vegas, depends on samples)

ROI Analysis#

Example: Monte Carlo Pricer (10,000 options/day):

Manual trapezoid rule:

• Development: 2 weeks × $5K/week = $10K
• Runtime: 10 ms/option × 10K = 100 seconds/day
• Debugging: 20 hours × $100/hour = $2K (first year)

scipy.integrate.quad:

• Development: 4 hours × $100/hour = $400
• Runtime: 0.1 ms/option × 10K = 1 second/day (100x faster via adaptive)
• Debugging: Minimal (library stability)

Savings: $11.6K first year, 99% runtime reduction

Implementation Reality#

Realistic Timeline Expectations#

Prototype (1-2 days):

• Install scipy/quadpy
• Test on sample integrals
• Validate accuracy vs analytic solutions
• Team: 1 developer/scientist

Production MVP (1-2 weeks):

• Integrate with existing codebase
• Error handling (singularities, convergence failures)
• Performance tuning (tolerance settings)
• Testing and validation
• Team: 1 developer + 1 domain expert

Optimized Production (1-3 months):

• Parallel execution (multi-core, GPU if applicable)
• Batch processing pipelines
• Monitoring and logging
• Team: 1-2 developers

Team Skill Requirements#

Minimum (Using scipy.quad):

• Python: Basic (functions, NumPy arrays)
• Numerical Methods: None (library abstracts algorithms)
• Domain Expertise: Medium (validate results)
• Training Time: 1-2 days

Typical (Using quadpy or vegas):

• Python: Moderate (classes, generators)
• Numerical Methods: Basic (understand convergence, error estimation)
• Domain Expertise: High
• Training Time: 1-2 weeks

Common Pitfalls#

Pitfall 1: “Adaptive means automatic accuracy”

• Reality: Must set tolerances (epsabs, epsrel) for scipy.quad
• Fix: Start with 1e-6 relative error, tighten if needed

Pitfall 2: “Monte Carlo is always better for high-D”

• Reality: Smooth 6D integrals may be faster with nested scipy.nquad
• Fix: Benchmark both (vegas vs nquad) for your problem

Pitfall 3: “One library solves everything”

• Reality: scipy (1-6D), quadpy (FEM domains), mpmath (precision), vegas (>6D)
• Fix: Choose based on dimensionality and precision needs

Pitfall 4: “Singularities will cause failures”

• Reality: scipy.quad handles most singularities via weight functions
• Fix: Use quad(..., weight='alg-loga', wvar=(alpha, beta)) for algebraic singularities

Key Takeaways for Decision Makers#

Top 3 Decisions to Make#

Decision 1: Dimensionality

• Rule: 1-6D → scipy.integrate (default), >6D → vegas (Monte Carlo)

Decision 2: Precision

• Rule: Standard (float64, ~1e-12) → scipy/quadpy, Arbitrary (50+ digits) → mpmath

Decision 3: Integration Domain

• Rule: Rectangles/standard intervals → scipy, Exotic (triangles, spheres) → quadpy

Budget Guidance#

Setup (One-Time):

• Developer time: 1-4 weeks × $5K/week = $5K-20K
• Infrastructure: $0 (runs on standard hardware, no GPU required)
• Total: $5K-20K

Ongoing (Per Year):

• Maintenance: 10-20 hours × $100/hour = $1K-2K
• Total: $1K-2K/year

ROI:

• Manual implementation: $10K-50K + weeks of debugging
• Library-based: $5K-20K + minimal debugging
• Payback: First project (weeks to months)

Questions to Ask Vendors/Consultants#

Technical Questions:

1. “Which library do you recommend: scipy, quadpy, or vegas? Why?” (Tests dimensionality understanding)
2. “How do you handle singularities?” (Tests scipy.quad weight functions knowledge)
3. “What’s the expected error for our integrals?” (Should quote 1e-6 to 1e-12 for scipy, O(1/√N) for vegas)

Business Questions:

1. “What’s the runtime for our problem size?” (Should give ballpark: µs for 1D, seconds for 6D, minutes for high-D)
2. “How do we validate accuracy?” (Should mention analytic test cases, convergence plots)
3. “What’s the maintenance burden?” (scipy: minimal updates, quadpy/vegas: more frequent)

Red Flags:

• ❌ Claims “any dimensionality, any precision” with one library (no such thing)
• ❌ Recommends custom implementation over scipy (reinventing the wheel)
• ❌ Doesn’t mention curse of dimensionality for >6D

Green Flags:

• ✅ Recommends scipy.quad first, vegas only if >6D
• ✅ Discusses tolerance settings (epsabs, epsrel) for error control
• ✅ Mentions QUADPACK heritage (scipy.quad = proven Fortran code)

Glossary#

Numerical Integration: Computing definite integrals via approximation (when symbolic integration is intractable)

Adaptive Quadrature: Automatic refinement of sampling points to meet error tolerance (scipy.quad’s approach)

Gauss-Kronrod: Nested quadrature rule for error estimation (7-15 point scheme in scipy.quad)

Monte Carlo Integration: Stochastic sampling for high-dimensional integrals (vegas library)

Curse of Dimensionality: Exponential cost growth for deterministic methods in high dimensions (why vegas needed for >6D)

QUADPACK: Fortran library (1970s-1980s) underlying scipy.integrate.quad (battle-tested, reliable)

Arbitrary Precision: Computation with 50-10,000 decimal digits (mpmath.quad)

scipy.integrate: Python’s standard numerical integration module (part of SciPy scientific stack)

quadpy: Modern high-order quadrature library (150+ schemes, exotic domains)

vegas: Adaptive Monte Carlo integration (high-dimensional integrals)

Singularity: Point where integrand is infinite or undefined (e.g., 1/x at x=0)

Tolerance: Acceptable error for numerical result (epsabs=1e-8 means ±1e-8 absolute error)

Further Reading#

Non-Technical:

• “Numerical Integration” (Wikipedia) - Accessible overview
• “The QUADPACK Story” (Piessens et al.) - History of scipy’s core algorithm

Technical:

• scipy.integrate documentation: https://docs.scipy.org/doc/scipy/reference/integrate.html
• quadpy documentation: https://github.com/sigma-py/quadpy
• vegas documentation: https://vegas.readthedocs.io/

Books:

• “Numerical Recipes” (Press et al.) - Classic numerical methods reference
• “Numerical Integration” (Davis & Rabinowitz) - Comprehensive theory

Community:

• SciPy mailing list: Scientific Python discussions
• Computational Science Stack Exchange: Numerical integration Q&A

S1 Approach: Rapid Discovery#

What S1 Discovers#

WHAT numerical integration libraries exist in Python for computing definite integrals?

S1 is an ecosystem scan: library positioning, algorithm coverage, performance characteristics, and specialized capabilities.

S1 Content Format#

For each library, document:

• Maturity: PyPI downloads, GitHub activity, production usage
• Algorithm Coverage: Which methods implemented (adaptive quadrature, Gauss-Kronrod, Romberg, Monte Carlo, etc.)
• Precision: Single/double/arbitrary precision support
• Dimensionality: 1D, 2D, 3D, or N-dimensional integrals
• Special Features: Singularity handling, infinite intervals, oscillatory integrands
• Performance: Speed benchmarks for typical integrals
• Best for: Quick positioning statement

What S1 Excludes#

❌ Installation instructions ❌ Code examples ❌ Configuration guides ❌ API documentation ❌ Usage tutorials

S1 helps you choose, not use

Reading Time#

5-10 minutes for complete ecosystem scan

mpmath.quad#

Positioning: Arbitrary-precision numerical integration in pure Python

Maturity#

• Stars: ~1.3K (mpmath/mpmath)
• Downloads: ~20M/month (PyPI)
• Maintainer: Fredrik Johansson, mpmath team
• Production: Symbolic computing, high-precision scientific calculations
• First Release: 2007

Algorithm Coverage#

1D Integration:

• Tanh-sinh quadrature (double exponential method)
• Gauss-Legendre quadrature
• Adaptive subdivision

Multi-dimensional:

• Nested 1D integration (via quad in loops)

Precision#

• Arbitrary precision: 1 to 10,000+ decimal digits
• Error estimates: Yes (adaptive tolerance control)
• Working precision: User-configurable via mp.dps

Dimensionality#

• 1D: Primary focus
• 2D+: Manual nesting (no built-in)

Special Features#

✅ Arbitrary precision: Main selling point (mpmath backend) ✅ Infinite intervals: (-∞, ∞), [a, ∞), (-∞, b] ✅ Oscillatory integrands: Tanh-sinh handles many cases well ✅ Singularities: Adaptive subdivision + endpoint handling ✅ Symbolic integration: Interop with SymPy

Performance#

1D integral (100 digits precision):

• ~10-100 ms (orders of magnitude slower than scipy)

Double precision equivalent:

• ~1-10 ms (still 10-100x slower than scipy.quad)

Trade-off: Precision over speed

API Complexity#

Simple: Similar to scipy.quad

from mpmath import mp, quad
mp.dps = 50  # 50 decimal places
result = quad(lambda x: mp.exp(-x**2), [0, mp.inf])

Best For#

• ✅ High-precision calculations (50-10,000 digits)
• ✅ Numerical verification (check symbolic results)
• ✅ Avoiding rounding errors (critical precision requirements)
• ✅ Research/theory (exact results, not engineering)

Not Ideal For#

• ❌ Speed-critical - 10-100x slower than scipy
• ❌ Multi-dimensional - no built-in support
• ❌ Large-scale - pure Python overhead
• ❌ Standard precision - scipy faster for float64

Dependencies#

Required: None (pure Python with gmpy2 acceleration optional) Optional: gmpy2 (faster arbitrary precision arithmetic)

License#

BSD-3-Clause (permissive, commercial-friendly)

quadpy#

Positioning: Modern high-order quadrature library with extensive algorithm collection

Maturity#

• Stars: ~800 (sigma-py/quadpy)
• Downloads: ~50K/month (PyPI)
• Maintainer: Nico Schlömer (active)
• Production: Research, engineering simulations, FEM/FVM codes
• First Release: 2017

Algorithm Coverage#

1D Integration (70+ schemes):

• Gauss-Legendre (arbitrary order)
• Clenshaw-Curtis
• Gauss-Kronrod (nested)
• Tanh-sinh (double exponential)
• Newton-Cotes (trapezoid, Simpson, Boole, etc.)

2D/3D Integration (150+ schemes):

• Triangle: Dunavant, Lyness-Jespersen, Strang-Fix
• Quadrilateral: Gauss product, Stroud
• Tetrahedron: Hammer-Stroud, Keast
• Hexahedron: Gauss product, Stroud
• Sphere: Lebedev, spherical designs
• Ball: Stroud, Hammer-Stroud

Multi-dimensional:

• Product rules (tensor products of 1D)
• Sparse grids (Smolyak construction)

Precision#

• Double precision: Standard
• Arbitrary precision: Via mpmath (all schemes)
• Error estimates: Some schemes (nested rules)

Dimensionality#

• 1D: Extensive (70+ schemes)
• 2D: Triangles, quads, disks, circles (150+ schemes)
• 3D: Tetrahedra, hexahedra, spheres, balls (100+ schemes)
• N-D: Product rules, sparse grids

Special Features#

✅ High-order schemes: Up to 50th order for many domains ✅ Exotic domains: Spheres, balls, pyramids, wedges, simplices ✅ Nested rules: Gauss-Kronrod for error estimation ✅ Arbitrary precision: Full mpmath support ✅ Tabulated weights: Pre-computed optimal nodes/weights

Performance#

1D Gauss-Legendre (order 10):

• ~5-10 µs (pre-computed weights)

2D Triangle (Dunavant order 10):

• ~20-50 µs (dozens of evaluation points)

High-order overhead: Weight table lookup + evaluation

API Complexity#

Moderate: Requires choosing scheme explicitly

val = quadpy.c1.gauss_legendre(n=10).integrate(f, [0.0, 1.0])

Benefit: Fine control over scheme selection

Best For#

• ✅ High-order accuracy (10th-50th order Gauss rules)
• ✅ Exotic domains (spheres, simplices, non-rectangular)
• ✅ FEM/FVM integration (standard element shapes)
• ✅ Arbitrary precision (mpmath backend)
• ✅ Algorithm research (extensive scheme catalog)

Not Ideal For#

• ❌ Adaptive integration - mostly fixed-order schemes
• ❌ High-dimensional (>6D) - no specialized N-D methods
• ❌ General-purpose - requires domain knowledge
• ❌ Simple use cases - scipy.integrate easier

Dependencies#

Required: numpy, sympy, ndim Optional: mpmath (arbitrary precision)

License#

MIT (permissive, commercial-friendly)

S1 Recommendation: Quick Choice Guide#

Decision Tree#

1-3 Dimensions, Standard Precision#

→ scipy.integrate.quad (default choice)

• Adaptive, accurate, fast
• Handles singularities and infinite intervals
• Battle-tested, stable API
• dblquad/tblquad for 2D/3D

High-Order Accuracy Needed (FEM/FVM, specialized domains)#

→ quadpy

• 10th-50th order Gauss schemes
• Triangles, tetrahedra, spheres
• Arbitrary precision via mpmath
• Requires domain knowledge

Arbitrary Precision (50+ digits)#

→ mpmath.quad

• 1-10,000 decimal places
• Pure Python (slower but precise)
• Verification and theory work

High-Dimensional (>6D)#

→ vegas

• Adaptive Monte Carlo
• Scales to 20-50+ dimensions
• Physics simulations (Feynman integrals, QCD)
• Statistical uncertainty acceptable

Quick Comparison#

Typical Workflow#

1. Try scipy.integrate.quad first (covers 80% of use cases)
2. Need >6D? → vegas
3. Need >50 digits? → mpmath.quad
4. Need exotic domain or ultra-high-order? → quadpy

Not Covered in S1#

• Symbolic integration: SymPy (not numerical)
• GPU acceleration: cupy, JAX (separate ecosystem)
• Specialized oscillatory: scipy.integrate doesn’t specialize, quadpy has some schemes
• Parallel frameworks: Dask integration (advanced)

Next Steps#

• S2: Deep dive into algorithm details, benchmarks, API patterns
• S3: Need-driven selection (your specific use case)
• S4: Strategic recommendation with code examples

scipy.integrate#

Positioning: Industry standard numerical integration in SciPy scientific computing stack

Maturity#

• Stars: ~13K (scipy/scipy repo)
• Downloads: 50M+/month (PyPI - scipy package)
• Maintainer: SciPy development community (NumFOCUS sponsored)
• Production: Used in scientific computing, engineering, finance, research worldwide
• First Release: 2001 (original), 2006 (modern scipy.integrate)

Algorithm Coverage#

1D Integration:

• quad: Adaptive Gauss-Kronrod quadrature (QUADPACK Fortran)
• romberg: Romberg integration (Richardson extrapolation)
• fixed_quad: Fixed-order Gaussian quadrature
• simpson: Simpson’s rule
• trapezoid: Trapezoidal rule

Multi-dimensional:

• dblquad: Double integrals (nested 1D)
• tblquad: Triple integrals
• nquad: N-dimensional integrals (arbitrary nesting)

Special:

• quad_vec: Vectorized quad (multiple integrands simultaneously)
• quadrature: Adaptive Gaussian quadrature (pure Python)

Precision#

• Double precision: Standard (float64)
• Arbitrary precision: Limited (via quad with mpmath)
• Error estimates: Yes (absolute/relative error)

Dimensionality#

• 1D: quad, romberg, simpson, trapezoid, fixed_quad
• 2D: dblquad
• 3D: tblquad
• N-D: nquad (up to ~6D practical)

Special Features#

✅ Singularities: Weight functions for algebraic/logarithmic singularities ✅ Infinite intervals: (-∞, ∞), [a, ∞), (-∞, b] ✅ Oscillatory integrands: Not specialized (use general adaptive) ✅ Vectorized: quad_vec for multiple integrands ✅ Symbolic integration: No (use SymPy)

Performance#

Typical 1D integral (smooth function, 0-1):

• quad: ~10-50 µs (adaptive, high accuracy)
• fixed_quad: ~5-10 µs (fixed order)
• simpson: ~20-40 µs (fixed samples)

Higher dimensions: Exponential scaling (~100x per dimension)

API Complexity#

Simple: One-liner for basic integrals

result, error = quad(lambda x: x**2, 0, 1)

Advanced: Weight functions, limits, tolerances

Best For#

• ✅ General-purpose 1D-3D integration (default choice)
• ✅ High accuracy (adaptive error control)
• ✅ Standard scientific Python stack (ships with SciPy)
• ✅ Production code (battle-tested, stable API)

Not Ideal For#

• ❌ High-dimensional (>6D) - exponential scaling (use Monte Carlo)
• ❌ Arbitrary precision - limited mpmath support (use mpmath.quad)
• ❌ Highly oscillatory - no specialized methods (use specialized libs)
• ❌ GPU acceleration - CPU only

Dependencies#

Required: numpy, scipy Optional: mpmath (arbitrary precision)

License#

BSD-3-Clause (permissive, commercial-friendly)

vegas#

Positioning: Adaptive Monte Carlo integration for high-dimensional integrals

Maturity#

• Stars: ~50 (gplepage/vegas)
• Downloads: ~20K/month (PyPI)
• Maintainer: G. Peter Lepage (Cornell)
• Production: High-energy physics, statistical mechanics, QCD calculations
• First Release: 2013

Algorithm Coverage#

Monte Carlo Methods:

• VEGAS algorithm (adaptive importance sampling)
• Adaptive stratified sampling
• Variance reduction techniques

Precision#

• Double precision: Standard
• Statistical uncertainty: Yes (Monte Carlo error estimates)
• Convergence: O(1/√N) where N = sample count

Dimensionality#

• 1D-20D+: Scales well to high dimensions
• Sweet spot: 5-20 dimensions
• Practical limit: ~50-100 dimensions (depends on integrand smoothness)

Special Features#

✅ High-dimensional: Main strength (>6D) ✅ Adaptive sampling: Learns important regions automatically ✅ Variance reduction: Importance sampling, stratification ✅ Parallel: Multi-core support (via nitn iterations) ✅ Warm-up: Training phase to optimize sampling grid

Performance#

6D integral (10^6 samples):

• ~1-10 seconds (depends on integrand complexity)

Scaling: O(N) with sample count (dimension-independent)

Accuracy: Improves as 1/√N (need 100x more samples for 10x accuracy)

API Complexity#

Moderate: Requires tuning iterations and sample counts

import vegas
integ = vegas.Integrator([[0, 1], [0, 1], [0, 1]])
result = integ(f, nitn=10, neval=1e4)

Best For#

• ✅ High-dimensional (>6D) - only practical method
• ✅ Smooth integrands - adaptive sampling effective
• ✅ Physics simulations - QCD, Feynman integrals
• ✅ Statistical uncertainty acceptable - Monte Carlo variance

Not Ideal For#

• ❌ Low-dimensional (<6D) - scipy.quad faster and more accurate
• ❌ High accuracy - need 10^9+ samples for 0.01% error
• ❌ Discontinuous integrands - slow convergence
• ❌ Real-time - stochastic, variable runtime

Dependencies#

Required: numpy, cython Optional: None

License#

GPLv3+ (copyleft - consider for commercial use)

S2 Approach: Comprehensive Analysis#

What S2 Discovers#

HOW do numerical integration algorithms work, and how do they compare?

S2 is a deep dive: algorithm internals, performance benchmarks, API patterns, edge case handling.

S2 Content Format#

For each library/algorithm:

• Algorithm Details: Mathematical foundation, convergence properties
• Performance Benchmarks: Timing across dimensions, accuracy vs cost trade-offs
• API Patterns: Code examples, parameter tuning
• Error Handling: Singularities, infinite intervals, oscillatory integrands
• Comparison: Head-to-head on common test integrals

What S2 Excludes#

❌ Production deployment (S4) ❌ Use case specific recommendations (S3) ❌ Superficial surveys (S1 already covered that)

S2 helps you understand how libraries work internally

Reading Time#

30-60 minutes for complete algorithmic understanding

S2 Recommendation: Algorithm Selection#

Key Insights from Deep Dive#

1. scipy.integrate.quad (QUADPACK) Strengths#

✅ Adaptive error control: Automatically refines where needed ✅ Battle-tested: 40+ years of production use ✅ Singularity handling: Weight functions for algebraic/log singularities ✅ Infinite intervals: Transforms to finite [0, 1]

When: Default choice for 1-3D integrals

2. quadpy High-Order Advantage#

✅ Fixed-cost predictability: Gauss-50 always uses 50 points (no subdivision unpredictability) ✅ FEM/FVM integration: Standard elements (triangles, tetrahedra) have optimal schemes ✅ Parallel-friendly: Fixed evaluation points → easy GPU/multi-core

When: High-order accuracy needed, FEM applications, GPU acceleration

3. mpmath.quad Precision Trade-off#

✅ Arbitrary precision: 50-10,000 digits ✅ Catastrophic cancellation: Avoids subtraction errors in float64

⚠️ Cost: 10-100x slower than scipy (pure Python vs Fortran)

When: Verification, theory, extreme precision requirements

4. vegas Monte Carlo Sweet Spot#

✅ Dimension-independent: O(N) cost regardless of D (5D = 20D = 50D) ✅ Adaptive sampling: Learns important regions

⚠️ Convergence: O(1/√N) slower than deterministic O(1/N^k)

When: >6 dimensions (only practical option)

Algorithm Decision Matrix#

Production Recommendations#

Scientific Python stack users: scipy.integrate (already installed) FEM/FVM developers: quadpy (optimal element integration) High-precision researchers: mpmath.quad (verification) Physics simulations: vegas (Feynman integrals, partition functions)

Next: S3 Need-Driven Discovery#

Match these algorithms to specific use cases (option pricing, Bayesian inference, etc.)

scipy.integrate: Algorithm Deep Dive#

Core Algorithm: QUADPACK#

scipy.integrate.quad wraps QUADPACK Fortran library (1970s-1980s Piessens, de Doncker, et al.)

Gauss-Kronrod Quadrature#

Idea: Nested quadrature rules for error estimation without wasting function evaluations

7-15 Point Rule (scipy default):

• 7-point Gauss-Legendre quadrature (initial estimate)
• 15-point Kronrod extension (refined estimate)
• Error = |15-point - 7-point|

Benefit: Reuses all 7 Gauss points in 15-point rule (only 8 new evaluations for error)

Adaptive Subdivision#

1. Compute G-K on [a, b]
2. If error < tolerance → done
3. Else: split [a, b] → [a, (a+b)/2], [(a+b)/2, b]
4. Recurse on each subinterval
5. Stop when global error < tolerance or max subdivisions reached

Key insight: Focus evaluations where integrand is complex (non-smooth regions)

Performance Characteristics#

Smooth integrands:

• 20-50 function evaluations typical
• Exponential convergence (error ∝ exp(-n))

Oscillatory/non-smooth:

• 100-1000+ evaluations
• Subdivision concentrates on difficult regions

Infinite intervals:

• Transform: ∫₀^∞ f(x) dx → ∫₀¹ f(t/(1-t)) / (1-t)² dt
• Maps infinite to finite [0, 1]

API Patterns#

from scipy.integrate import quad

# Basic: integral of x^2 from 0 to 1
result, error = quad(lambda x: x**2, 0, 1)
# result = 0.333..., error ~ 1e-14

# Infinite interval
result, error = quad(lambda x: np.exp(-x**2), 0, np.inf)

# Singularity: ∫₀¹ x^(-0.5) dx (diverges at x=0)
result, error = quad(lambda x: x**(-0.5), 0, 1,
                      weight='alg', wvar=(-0.5, 0))

Convergence Analysis#

Target accuracy: 1e-6 to 1e-12 typical (engineering/science)

Achieving 1e-8:

• Smooth: ~30 evaluations
• Oscillatory: ~200 evaluations
• Singularity: ~50-100 evaluations (with weight functions)

When QUADPACK Fails#

High-dimensional (>3D): Exponential cost Highly oscillatory: No specialized methods (generic adaptive) Discontinuous: Subdivision helps but slow convergence Arbitrary precision: Limited mpmath support

→ Use alternative libraries for these cases

S3 Approach: Need-Driven Discovery#

What S3 Discovers#

WHY would I use numerical integration for MY specific use case?

S3 is problem-first: concrete use cases, requirements, library matching.

S3 Content Format#

For each use case:

• Use Case Description: Real-world problem
• Technical Requirements: Dimensions, precision, performance
• Library Recommendation: Which library + why
• Alternative Approaches: When NOT to use numerical integration
• Example Code: Minimal working example

What S3 Excludes#

❌ General algorithm theory (S2 covered that) ❌ Production patterns (S4 strategic) ❌ Library surveys (S1 rapid)

S3 helps you answer: “Is numerical integration right for MY problem?”

Reading Time#

10-20 minutes to find and evaluate your use case

Use Case: Monte Carlo Option Pricing#

Problem Description#

Financial engineering: Price European/exotic options when closed-form solutions don’t exist

Example: Barrier option with path-dependent payoff

• Must integrate over stochastic paths (high-dimensional)
• No Black-Scholes closed form

Technical Requirements#

Dimensionality: 10-100 time steps × multiple assets = 10-1000 dimensions Precision: 0.1-1% accuracy (financial tolerance) Performance: Seconds to minutes for portfolio (1000s of options) Constraints: Reproducible pricing (regulatory), confidence intervals

Library Recommendation#

Primary: vegas (adaptive Monte Carlo)

• Handles 10-100D easily
• Statistical error estimates (confidence intervals)
• Adaptive importance sampling learns option payoff regions

Alternative: scipy.integrate.nquad (nested 1D)

• Only for 2-4 assets (6-12D max)
• Deterministic (reproducible)
• Slower convergence than vegas for >6D

Why Not Symbolic/Analytic?#

❌ SymPy integration: Path-dependent payoffs have no closed form ❌ Closed-form: Black-Scholes only for vanilla European (not exotic)

Example Code#

import vegas
import numpy as np

def barrier_payoff(x):
    """Asian barrier call: max(0, avg(S_t) - K) if never hit barrier"""
    # x = [S_1, S_2, ..., S_T] normalized paths
    S = 100 * np.exp(np.cumsum(x))  # geometric Brownian motion
    if np.any(S < 80):  # barrier hit
        return 0
    return max(0, np.mean(S) - 100)  # Asian call payoff

# 20-dimensional integral (20 time steps)
integ = vegas.Integrator([[- 0.5, 0.5]] * 20)
result = integ(barrier_payoff, nitn=10, neval=1e5)
print(f"Option value: ${result.mean:.2f} ± ${result.sdev:.2f}")

Performance Expectations#

10D integral (10 time steps):

• 10⁵ samples: ~1 second, ±1% error
• 10⁶ samples: ~10 seconds, ±0.3% error

100D integral (100 time steps):

• 10⁶ samples: ~10 seconds, ±1% error (dimension-independent!)
• 10⁷ samples: ~100 seconds, ±0.3% error

Trade-offs#

Monte Carlo (vegas): ✅ Scales to 100D+ ✅ Statistical error estimates ⚠️ Stochastic (different results each run)

Deterministic (scipy.nquad): ✅ Reproducible ❌ Exponential cost (unusable >6D)

Verdict: Monte Carlo is only viable approach for real options (typically 20-100D)

S3 Recommendation: Use Case Matching#

Quick Use Case Lookup#

Financial Engineering#

• Option pricing (10-100D): vegas
• Risk integrals (VaR, 1-3D): scipy.integrate.quad
• Expectation computation (1-5D): scipy or vegas (depends on D)

Bayesian Inference#

• Posterior normalization (1-3D): scipy.integrate.quad
• High-D posteriors (>5D): vegas or MCMC (sampling, not integration)
• Evidence computation: scipy (low-D) or vegas (high-D)

Physics Simulations#

• Feynman path integrals (20-100D): vegas (standard in HEP)
• Partition functions (10-50D): vegas
• 1D integrals (quantum mechanics): scipy.integrate.quad

Computational Geometry#

• Volume of meshes (FEM): quadpy (tetrahedron integration)
• Area under parametric curves: scipy.integrate.quad
• 3D CAD intersections: quadpy (triangle/tet schemes)

Engineering#

• Structural analysis (FEM): quadpy (element integration)
• Control theory (Lyapunov functions, 1-3D): scipy.integrate.quad
• Signal processing (Fourier integrals, 1D): scipy.integrate.quad

Decision Heuristic#

if dimensions <= 3:
    use scipy.integrate.quad/dblquad/tblquad
elif dimensions <= 6 and integrand smooth:
    try scipy.integrate.nquad first
    if too slow, use vegas
elif dimensions > 6:
    use vegas (Monte Carlo only practical option)

if arbitrary_precision_needed:
    use mpmath.quad (1D only, slow)

if FEM_element_integration:
    use quadpy (standard element schemes)

When NOT to Use Numerical Integration#

❌ Symbolic solution exists: Use SymPy (exact answer) ❌ Data already discrete: Use NumPy.sum (no integral needed) ❌ High-D sampling problem: Use MCMC (Markov Chain Monte Carlo) instead of integration ❌ Optimization problem: Use scipy.optimize (not integration)

Next: S4 Strategic Selection#

Production patterns, dependency management, testing strategies

S4 Approach: Strategic Selection#

What S4 Delivers#

HOW do I deploy numerical integration in production?

S4 is production-ready: dependency management, error handling, testing, monitoring.

S4 Content Format#

• Final Recommendation: Specific library for specific context
• Dependency Management: pip/conda, version pinning
• Error Handling: Convergence failures, singularities, edge cases
• Testing Strategy: Validation against analytic solutions
• Code Examples: Production-ready patterns
• Monitoring: Performance tracking, accuracy verification

What S4 Excludes#

❌ Algorithm theory (S2) ❌ Use case exploration (S3) ❌ Library survey (S1)

S4 provides: Copy-pasteable production patterns

Reading Time#

20-40 minutes for complete production deployment guide

Production Deployment: scipy.integrate#

Final Recommendation#

Use scipy.integrate.quad for 1-3D integrals in production

✅ Mature (40+ years QUADPACK heritage) ✅ Stable API (scipy 0.x → 1.x compatibility) ✅ Minimal dependencies (numpy, scipy already in scientific stack) ✅ Excellent documentation and community support

Dependency Management#

# pip
pip install scipy>=1.7.0

# conda (recommended for scientific stack)
conda install scipy>=1.7.0

# Poetry/requirements.txt
scipy>=1.7.0,<2.0.0  # Pin major version

Version strategy:

• scipy 1.7+ (Python 3.7+)
• Pin major version (avoid breaking changes)
• Update minor versions regularly (bug fixes, performance)

Error Handling Patterns#

from scipy.integrate import quad, IntegrationWarning
import warnings

def safe_integrate(f, a, b, **kwargs):
    """Production-safe integration with error handling"""
    try:
        with warnings.catch_warnings():
            warnings.simplefilter("error", IntegrationWarning)
            result, error = quad(f, a, b, **kwargs)

        # Verify error is acceptable
        if error > kwargs.get('epsrel', 1e-8) * abs(result):
            warnings.warn(f"Integration error {error:.2e} exceeds tolerance")

        return result, error

    except IntegrationWarning as e:
        # Convergence failure - subdivide or increase limit
        return quad(f, a, b, limit=1000, **kwargs)

    except Exception as e:
        # Log and re-raise with context
        logging.error(f"Integration failed for f={f.__name__}, bounds=[{a}, {b}]: {e}")
        raise

Testing Strategy#

Unit tests: Validate against analytic solutions

import numpy as np
from scipy.integrate import quad

def test_polynomial_integral():
    """∫₀¹ x² dx = 1/3"""
    result, error = quad(lambda x: x**2, 0, 1)
    assert np.isclose(result, 1/3, rtol=1e-8)
    assert error < 1e-10

def test_gaussian_integral():
    """∫₀^∞ e^(-x²) dx = √π/2"""
    result, error = quad(lambda x: np.exp(-x**2), 0, np.inf)
    assert np.isclose(result, np.sqrt(np.pi)/2, rtol=1e-8)

Integration tests: Verify convergence on production integrals

def test_convergence_tolerance():
    """Ensure error estimates are accurate"""
    f = lambda x: np.sin(1/x) if x != 0 else 0  # tricky integrand

    # Different tolerances should converge
    r1, e1 = quad(f, 0.01, 1, epsrel=1e-6)
    r2, e2 = quad(f, 0.01, 1, epsrel=1e-10)

    assert abs(r1 - r2) < e1  # Coarse result within error of fine result

Performance Monitoring#

import time
import logging

def monitored_integrate(f, a, b, **kwargs):
    """Integration with performance logging"""
    start = time.time()
    result, error = quad(f, a, b, **kwargs)
    elapsed = time.time() - start

    # Log metrics
    logging.info(f"Integration: {elapsed*1000:.2f}ms, error={error:.2e}")

    # Alert on slow integrals
    if elapsed > 1.0:  # > 1 second
        logging.warning(f"Slow integration: {elapsed:.2f}s for {f.__name__}")

    return result, error

Production Patterns#

Batch Processing#

def integrate_batch(funcs, a, b):
    """Integrate multiple functions in parallel"""
    from multiprocessing import Pool

    with Pool() as pool:
        args = [(f, a, b) for f in funcs]
        results = pool.starmap(quad, args)

    return results

Caching for Repeated Integrals#

from functools import lru_cache

@lru_cache(maxsize=128)
def cached_integrate(f_key, a, b, epsrel=1e-8):
    """Cache integration results for repeated queries"""
    # Note: f must be hashable (use string key for function ID)
    f = FUNCTION_REGISTRY[f_key]
    result, error = quad(f, a, b, epsrel=epsrel)
    return result

Common Pitfalls#

Pitfall 1: Ignoring Error Estimates#

❌ Bad: result, _ = quad(f, 0, 1) # Ignore error ✅ Good: Check error < tolerance

Pitfall 2: Default Tolerance Too Loose#

❌ Bad: Use default (1.49e-8) for financial calculations ✅ Good: Specify epsrel=1e-10 or tighter

Pitfall 3: Infinite Interval Without Care#

❌ Bad: quad(f, 0, np.inf) for slow-decaying f ✅ Good: Test finite truncation first, then infinite

Deployment Checklist#

• scipy >= 1.7.0 installed
• Error handling wrapper implemented
• Unit tests for analytic test cases
• Integration tests for production integrals
• Performance monitoring/logging
• Tolerance settings documented
• Convergence failure handling tested

Next Steps#

For multi-dimensional (2D-3D), use dblquad/tblquad with same patterns. For high-D (>6D), switch to vegas (separate deployment guide).

S4 Strategic Recommendation#

Production-Ready Summary#

Primary Recommendation: scipy.integrate#

For 95% of use cases (1-3D integrals, standard precision):

from scipy.integrate import quad

# Production pattern
result, error = quad(f, a, b, epsrel=1e-10, limit=100)
assert error < 1e-9, f"Integration failed: error={error}"

Why scipy.integrate:

• ✅ Installed everywhere (part of scientific Python stack)
• ✅ 40+ years battle-tested (QUADPACK Fortran)
• ✅ Stable API (scipy 0.x → 1.x compatibility)
• ✅ Excellent docs, community support

Dependencies: numpy, scipy (likely already installed)

Alternative: quadpy (Specialized)#

For FEM/FVM, high-order schemes (triangles, tetrahedra, exotic domains):

import quadpy

# Triangle integration (FEM)
scheme = quadpy.t2.get_good_scheme(10)  # 10th order
val = scheme.integrate(f, triangle_vertices)

When to use:

• FEM/FVM element integration
• Need 20th-50th order accuracy
• Exotic domains (spheres, pyramids)

Dependencies: numpy, sympy, quadpy

Alternative: mpmath.quad (High-Precision)#

For arbitrary precision (50-10,000 digits):

from mpmath import mp, quad

mp.dps = 100  # 100 decimal places
result = quad(f, [0, 1])

When to use:

• Numerical verification of symbolic results
• Catastrophic cancellation in float64
• Theory/research

Trade-off: 10-100x slower than scipy

Alternative: vegas (High-Dimensional)#

For >6D integrals (only practical option):

import vegas

integ = vegas.Integrator([[0, 1]] * 20)  # 20D
result = integ(f, nitn=10, neval=1e6)

When to use:

• Monte Carlo option pricing (10-100D)
• Physics simulations (Feynman integrals)
• Bayesian evidence (high-D posteriors)

Trade-off: O(1/√N) convergence (need 100x more samples for 10x accuracy)

Deployment Decision Matrix#

Production Checklist#

scipy.integrate (recommended default):

• scipy >= 1.7.0 installed
• Error handling (IntegrationWarning)
• Unit tests vs analytic solutions
• Tolerance settings documented (epsrel, epsabs)
• Performance monitoring

quadpy:

• Scheme selection strategy (order vs cost)
• Element geometry validation
• Arbitrary precision config (if using mpmath)

mpmath:

• Precision requirements documented (mp.dps)
• Performance testing (pure Python overhead)
• Validate gmpy2 acceleration installed

vegas:

• Sample count tuning (nitn, neval)
• Warm-up strategy (adaption iterations)
• Statistical error validation
• Seed management (reproducibility)

Final Advice#

Start simple: scipy.integrate.quad covers 95% of needs

Scale up: Only adopt specialized libraries when scipy.integrate limits reached

• High-D (>6D) → vegas
• High-precision (>15 digits) → mpmath
• FEM/exotic domains → quadpy

Production patterns:

• Error handling mandatory (convergence failures happen)
• Unit tests against analytic solutions
• Monitor performance (slow integrals indicate problems)
• Document tolerance settings (epsrel, epsabs)

Avoid premature optimization: scipy.quad is fast enough for most applications (µs-ms per integral)