Quadrature Comparison Example — Accuracy vs. Evaluations¶
This example compares quadrature strategies from
integrationdifferentiation for integrating
\(\int_{-3}^{3} e^{-x^2}\,dx \approx \sqrt{\pi} \approx 1.77245\):
Gauss-Legendre rules with 4, 8, 16, and 32 nodes — IMSL GQRUL (
integrationdifferentiation.gauss_quadrature_rule())Adaptive quadrature at tolerances \(10^{-4}\), \(10^{-8}\), \(10^{-12}\) — IMSL QDAG (
integrationdifferentiation.integrate_adaptive())Non-adaptive 10-point Gauss-Kronrod — IMSL QDNG (
integrationdifferentiation.integrate_nonadaptive())
Example Code¶
"""IMSL quadrature comparison: accuracy vs. function evaluations.
Demonstrates several quadrature strategies from
:mod:`integrationdifferentiation` for integrating ``exp(-x²)`` over
``[-3, 3]``:
- Gauss-Legendre rules with 4, 8, 16, and 32 nodes (IMSL GQRUL)
- Adaptive quadrature (IMSL QDAG)
- Non-adaptive Gauss-Kronrod rule (IMSL QDAG nonadaptive)
Exact answer: ``∫_{-3}^{3} exp(-x²) dx = erf(3) * √π ≈ 1.77245``
Outputs:
- Accuracy table printed to stdout
- SVG plot saved to ``test_output/example_imsl_quadrature_comparison.svg``
"""
from __future__ import annotations
from pathlib import Path
from typing import Dict, List
import matplotlib.pyplot as plt
import numpy as np
from integrationdifferentiation import (
gauss_quadrature_rule,
integrate_adaptive,
integrate_nonadaptive,
)
def _apply_gl_rule(n: int, a: float, b: float) -> tuple[float, int]:
"""Apply an n-point Gauss-Legendre rule to integrate exp(-x²) on [a,b].
Args:
n (int): Number of Gauss-Legendre nodes.
a (float): Lower integration limit in Pa units — here dimensionless x.
b (float): Upper integration limit.
Returns:
tuple[float, int]: (integral estimate, number of function evaluations).
"""
rule = gauss_quadrature_rule(n, weight_type="legendre")
# Affine transform: x = 0.5*(b-a)*t + 0.5*(b+a), t in [-1,1]
mid = 0.5 * (b + a)
half = 0.5 * (b - a)
val = half * float(
np.sum(rule.weights * np.exp(-(mid + half * rule.points) ** 2))
)
return val, n
def run_demo_imsl_quadrature_comparison() -> Dict[str, object]:
"""Compare Gauss-Legendre, adaptive, and non-adaptive quadrature accuracy.
Integrates ``exp(-x²)`` on ``[-3, 3]`` using multiple methods and
reports absolute error vs. number of function evaluations.
Args:
None
Returns:
Dict[str, object]: Keys ``rows`` (list of result dicts),
``exact`` (float), ``plot_path`` (str).
"""
from scipy.special import erf
a, b = -3.0, 3.0
exact = float(erf(3.0) * np.sqrt(np.pi))
def f(x):
return np.exp(-np.asarray(x, dtype=float) ** 2)
rows: List[dict] = []
# Gauss-Legendre rules: 4, 8, 16, 32 nodes
for n in (4, 8, 16, 32):
val, nfev = _apply_gl_rule(n, a, b)
rows.append({"method": f"Gauss-Legendre n={n}", "value": val, "n_fev": nfev, "error": abs(val - exact)})
# Adaptive quadrature (QDAG)
for tol in (1e-4, 1e-8, 1e-12):
r = integrate_adaptive(f, a, b, abs_tol=tol, rel_tol=tol)
rows.append({"method": f"Adaptive (tol={tol:.0e})", "value": r.value, "n_fev": r.n_fev, "error": abs(r.value - exact)})
# Non-adaptive (fixed Gauss-Kronrod)
r_na = integrate_nonadaptive(f, a, b)
rows.append({"method": "Non-adaptive (G-K)", "value": r_na.value, "n_fev": r_na.n_fev, "error": abs(r_na.value - exact)})
print("\nIMSL Quadrature Comparison: ∫_{-3}^{3} exp(-x²) dx")
print(f"Exact value: {exact:.10f}")
print("=" * 72)
print(f"{'Method':<28} {'Result':>14} {'|Error|':>14} {'n_fev':>8}")
print("-" * 72)
for row in rows:
print(f"{row['method']:<28} {row['value']:>14.10f} {row['error']:>14.3e} {row['n_fev']:>8}")
print("=" * 72)
# --------------------------------------------------------------- plot
output_dir = Path("test_output")
output_dir.mkdir(parents=True, exist_ok=True)
plot_path = output_dir / "example_imsl_quadrature_comparison.svg"
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))
# Panel 1: integrand and GL nodes for n=8
x_plot = np.linspace(a, b, 400)
rule8 = gauss_quadrature_rule(8, weight_type="legendre")
mid, half = 0.5 * (b + a), 0.5 * (b - a)
nodes8 = mid + half * rule8.points
ax1.plot(x_plot, np.exp(-x_plot ** 2), color="#1d4ed8", lw=2.0, label=r"$e^{-x^2}$")
ax1.fill_between(x_plot, np.exp(-x_plot ** 2), alpha=0.15, color="#1d4ed8")
ax1.scatter(nodes8, np.exp(-nodes8 ** 2), color="#dc2626", zorder=5, s=60, label="GL-8 nodes")
ax1.axhline(0, color="black", lw=0.6)
ax1.set_xlabel("x")
ax1.set_ylabel(r"$e^{-x^2}$")
ax1.set_title(r"Integrand $e^{-x^2}$ on $[-3, 3]$ with 8-pt GL nodes")
ax1.legend()
ax1.grid(True, alpha=0.3)
# Panel 2: error vs n_fev (convergence plot)
gl_rows = [r for r in rows if r["method"].startswith("Gauss-Legendre")]
ad_rows = [r for r in rows if r["method"].startswith("Adaptive")]
na_rows = [r for r in rows if r["method"].startswith("Non-adaptive")]
ax2.loglog(
[r["n_fev"] for r in gl_rows],
[max(r["error"], 1e-16) for r in gl_rows],
"o-", color="#059669", lw=2.0, ms=7, label="Gauss-Legendre",
)
ax2.loglog(
[r["n_fev"] for r in ad_rows],
[max(r["error"], 1e-16) for r in ad_rows],
"s--", color="#d97706", lw=2.0, ms=7, label="Adaptive (QDAG)",
)
ax2.loglog(
[r["n_fev"] for r in na_rows],
[max(r["error"], 1e-16) for r in na_rows],
"^", color="#7c3aed", ms=10, label="Non-adaptive (G-K)",
)
ax2.set_xlabel("Number of function evaluations")
ax2.set_ylabel("Absolute error")
ax2.set_title("Quadrature accuracy vs. evaluations")
ax2.legend()
ax2.grid(True, alpha=0.3, which="both")
fig.suptitle("IMSL Quadrature Method Comparison", fontsize=13, fontweight="bold")
fig.tight_layout()
fig.savefig(plot_path, format="svg")
plt.close(fig)
return {"rows": rows, "exact": exact, "plot_path": str(plot_path)}
if __name__ == "__main__":
run_demo_imsl_quadrature_comparison()
Plot Output¶
Left: integrand \(e^{-x^2}\) on \([-3,3]\) with the 8-point Gauss-Legendre sample nodes marked. Right: log–log convergence plot — absolute error vs. number of function evaluations for each method.¶
Console Output¶
IMSL Quadrature Comparison: Γê½_{-3}^{3} exp(-x┬▓) dx
Exact value: 1.7724146965
========================================================================
Method Result |Error| n_fev
------------------------------------------------------------------------
Gauss-Legendre n=4 1.3852665803 3.871e-01 4
Gauss-Legendre n=8 1.7688308634 3.584e-03 8
Gauss-Legendre n=16 1.7724146933 3.241e-09 16
Gauss-Legendre n=32 1.7724146965 3.331e-15 32
Adaptive (tol=1e-04) 1.7724146965 6.661e-16 63
Adaptive (tol=1e-08) 1.7724146965 6.661e-16 63
Adaptive (tol=1e-12) 1.7724146965 2.220e-16 147
Non-adaptive (G-K) 1.7722370822 1.776e-04 10
========================================================================