QDAG Example — Adaptive Quadrature

This example demonstrates general-purpose adaptive integration (IMSL QDAG / QDAGI):

• Integrate exp(-x²) from 0 to ∞ → √π / 2 ≈ 0.88623

• Integrate sin(x)/x from 0 to π → Si(π) ≈ 1.85194

Example Code

"""IMSL QDAG example: univariate adaptive quadrature.

Demonstrates general-purpose adaptive integration (IMSL QDAG / QDAGI):

- Integrate ``exp(-x²)`` from ``0`` to ``∞`` → ``√π / 2 ≈ 0.88623``
- Integrate ``sin(x)/x`` from ``0`` to ``π`` → ``Si(π) ≈ 1.85194``

Outputs:
- Table printed to stdout
- SVG plot saved to ``test_output/demo_imsl_qdag.svg``
"""

from __future__ import annotations

from pathlib import Path
from typing import Dict

import matplotlib.pyplot as plt
import numpy as np

from integrationdifferentiation import integrate_adaptive, integrate_adaptive_infinite


def run_demo_imsl_qdag() -> Dict[str, object]:
    """Run IMSL QDAG adaptive quadrature examples.

    Args:
        None

    Returns:
        Dict[str, object]: Result dict with keys ``results`` (list of result
            objects) and ``plot_path`` (str).
    """
    def f1(x: float) -> float:
        return float(np.exp(-x ** 2))

    def f2(x: float) -> float:
        return float(np.sinc(x / np.pi))  # sinc(x/π) = sin(x)/x

    # Case 1: exp(-x²) on [0, ∞) — exact = sqrt(π)/2
    r1 = integrate_adaptive_infinite(f1, a=0.0, b=np.inf)
    exact1 = float(np.sqrt(np.pi) / 2.0)

    # Case 2: sin(x)/x on [0, π] — avoid division by zero at x=0
    def f2_safe(x: float) -> float:
        return 1.0 if x == 0.0 else float(np.sin(x) / x)

    r2 = integrate_adaptive(f2_safe, 0.0, np.pi)
    # Si(π) computed via scipy
    from scipy.special import sici
    exact2 = float(sici(np.pi)[0])

    results = [r1, r2]

    print("\nIMSL QDAG Example: Adaptive Quadrature")
    print("=" * 72)
    print(f"{'Function':<28} {'Result':>12} {'Error':>12} {'Exact':>12} {'n_fev':>6}")
    print("-" * 72)
    print(
        f"{'∫exp(-x²) dx, 0→∞':<28} {r1.value:>12.8f} {r1.error:>12.2e}"
        f" {exact1:>12.8f} {r1.n_fev:>6}"
    )
    print(
        f"{'∫sin(x)/x dx, 0→π':<28} {r2.value:>12.8f} {r2.error:>12.2e}"
        f" {exact2:>12.8f} {r2.n_fev:>6}"
    )
    print("=" * 72)

    output_dir = Path("test_output")
    output_dir.mkdir(parents=True, exist_ok=True)
    plot_path = output_dir / "demo_imsl_qdag.svg"

    fig, axes = plt.subplots(1, 2, figsize=(12, 5))

    # Panel 1: exp(-x²)
    ax = axes[0]
    x1 = np.linspace(0, 3.5, 400)
    y1 = np.exp(-x1 ** 2)
    ax.plot(x1, y1, color="#059669", linewidth=2.0, label=r"$f(x) = e^{-x^2}$")
    ax.fill_between(x1, y1, alpha=0.18, color="#059669")
    ax.axhline(0, color="black", linewidth=0.8)
    ax.set_xlabel("x")
    ax.set_ylabel("f(x)")
    ax.set_title(r"$\int_0^\infty e^{-x^2}\,dx = \sqrt{\pi}/2 \approx 0.88623$")
    ax.legend()
    ax.grid(True, alpha=0.3)
    ax.annotate(
        f"Result: {r1.value:.6f}\nError:  {r1.error:.2e}",
        xy=(0.98, 0.95),
        xycoords="axes fraction",
        ha="right",
        va="top",
        fontsize=9,
        bbox=dict(boxstyle="round,pad=0.3", fc="#ecfdf5", ec="#a7f3d0"),
    )

    # Panel 2: sin(x)/x
    ax = axes[1]
    x2 = np.linspace(1e-6, np.pi, 400)
    y2 = np.sin(x2) / x2
    ax.plot(x2, y2, color="#10b981", linewidth=2.0, label=r"$f(x) = \sin(x)/x$")
    ax.fill_between(x2, y2, alpha=0.18, color="#10b981")
    ax.axhline(0, color="black", linewidth=0.8)
    ax.set_xlabel("x")
    ax.set_ylabel("f(x)")
    ax.set_title(r"$\int_0^\pi \frac{\sin x}{x}\,dx = \mathrm{Si}(\pi) \approx 1.85194$")
    ax.legend()
    ax.grid(True, alpha=0.3)
    ax.annotate(
        f"Result: {r2.value:.6f}\nError:  {r2.error:.2e}",
        xy=(0.98, 0.95),
        xycoords="axes fraction",
        ha="right",
        va="top",
        fontsize=9,
        bbox=dict(boxstyle="round,pad=0.3", fc="#ecfdf5", ec="#a7f3d0"),
    )

    fig.suptitle("IMSL QDAG: Adaptive Quadrature Examples", fontsize=13, fontweight="bold")
    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)

    return {"results": results, "plot_path": str(plot_path)}


if __name__ == "__main__":
    run_demo_imsl_qdag()

Plot Output

Left: \(\int_0^\infty e^{-x^2}\,dx = \sqrt{\pi}/2\). Right: \(\int_0^\pi \frac{\sin x}{x}\,dx = \mathrm{Si}(\pi)\). Both panels annotate the computed result and error estimate.

Console Output

IMSL QDAG Example: Adaptive Quadrature
========================================================================
Function                           Result        Error        Exact  n_fev
------------------------------------------------------------------------
∫exp(-x²) dx, 0→∞              0.88622693     7.10e-09   0.88622693    135
∫sin(x)/x dx, 0→π              1.85193705     2.06e-14   1.85193705     21
========================================================================