Improper Integrals Example — Infinite-Domain Quadrature¶
This example demonstrates integrationdifferentiation.integrate_adaptive_infinite(),
integrationdifferentiation.integrate_weighted_fourier(), and
integrationdifferentiation.integrate_adaptive_endpoint_singular()
(IMSL QDAGI / QDAWF / QDAGS):
\(\int_0^\infty e^{-x^2}\,dx = \sqrt{\pi}/2 \approx 0.88623\)
\(\int_0^\infty \frac{\sin x}{x}\,dx = \pi/2 \approx 1.57080\) (Dirichlet integral via weighted Fourier)
\(\int_0^1 \frac{\ln x}{\sqrt{x}}\,dx = -4\) (endpoint singularity at \(x = 0\))
Example Code¶
"""IMSL QDAGI example: improper integrals over infinite domains.
Demonstrates :func:`integrationdifferentiation.integrate_adaptive_infinite`
and :func:`integrationdifferentiation.integrate_adaptive_endpoint_singular`
(IMSL QDAGI / QDAGS):
- ``∫₀^∞ exp(-x²) dx = √π/2 ≈ 0.88623``
- ``∫₀^∞ sin(x)/x dx = π/2 ≈ 1.57080`` (Dirichlet integral)
- ``∫₀^1 ln(x)/√x dx = −4`` (endpoint singularity at x=0)
Outputs:
- Table printed to stdout
- SVG plot saved to ``test_output/example_imsl_improper_integrals.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,
integrate_adaptive_endpoint_singular,
integrate_weighted_fourier,
)
def run_demo_imsl_improper_integrals() -> Dict[str, object]:
"""Run three improper integral examples and visualise convergence.
Args:
None
Returns:
Dict[str, object]: Keys ``results`` (list of dicts with value/exact/error),
``plot_path`` (str).
"""
# ---------------------------------------------------------------- case 1
# ∫₀^∞ exp(-x²) dx = √π / 2
def f1(x: float) -> float:
return float(np.exp(-float(x) ** 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 dx = π/2 — use weighted Fourier quadrature for precision
# Write as ∫_{ε}^∞ (1/x) * sin(1*x) dx ≈ π/2 (missing strip [0,ε] ≈ ε)
def one_over_x(x: float) -> float:
x = float(x)
return 1.0 / x if x > 1e-12 else 1.0
r2 = integrate_weighted_fourier(one_over_x, a=1e-10, omega=1.0, weight="sin")
exact2 = float(np.pi / 2.0)
# ---------------------------------------------------------------- case 3
# ∫₀^1 ln(x)/√x dx = -4 (endpoint singularity: ln(x)→-∞ as x→0+)
def f3(x: float) -> float:
x = float(x)
if x <= 0.0:
return 0.0
return float(np.log(x) / np.sqrt(x))
r3 = integrate_adaptive_endpoint_singular(f3, a=0.0, b=1.0)
exact3 = -4.0
cases = [
{"label": "∫₀^∞ exp(-x²) dx", "result": r1, "exact": exact1},
{"label": "∫₀^∞ sin(x)/x dx", "result": r2, "exact": exact2},
{"label": "∫₀^1 ln(x)/√x dx", "result": r3, "exact": exact3},
]
print("\nIMSL QDAGI Example: Improper Integrals")
print("=" * 78)
print(f"{'Integral':<26} {'Result':>12} {'Error':>12} {'Exact':>12} {'|Δ|':>12}")
print("-" * 78)
for c in cases:
r = c["result"]
delta = abs(r.value - c["exact"])
print(
f"{c['label']:<26} {r.value:>12.8f} {r.error:>12.2e}"
f" {c['exact']:>12.8f} {delta:>12.2e}"
)
print("=" * 78)
# -------------------------------------------------------- convergence data
# Show how partial integrals ∫₀^T converge as T grows for case 1 and 2
T_vals = np.linspace(0.5, 15.0, 60)
partial1 = []
for T in T_vals:
res = integrate_adaptive_infinite(f1, a=0.0, b=float(T))
partial1.append(res.value)
partial2 = []
def sinc_safe(x: float) -> float:
x = float(x)
return 1.0 if x < 1e-12 else float(np.sin(x) / x)
for T in T_vals:
res = integrate_adaptive(sinc_safe, 1e-8, float(T))
partial2.append(res.value)
# ----------------------------------------------------------------- plots
output_dir = Path("test_output")
output_dir.mkdir(parents=True, exist_ok=True)
plot_path = output_dir / "example_imsl_improper_integrals.svg"
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
# --- panel 1: exp(-x²) integrand
ax = axes[0]
x1 = np.linspace(0, 4, 400)
ax.plot(x1, np.exp(-x1 ** 2), color="#059669", lw=2.0, label=r"$e^{-x^2}$")
ax.fill_between(x1, np.exp(-x1 ** 2), alpha=0.2, color="#059669")
ax.axhline(0, color="black", lw=0.6)
ax.set_xlabel("x")
ax.set_title(r"$\int_0^\infty e^{-x^2}\,dx = \sqrt{\pi}/2$")
ax.legend()
ax.grid(True, alpha=0.3)
ax.annotate(
f"Numerical: {r1.value:.6f}\nExact: {exact1:.6f}",
xy=(0.97, 0.95), xycoords="axes fraction", ha="right", va="top",
fontsize=8, bbox=dict(boxstyle="round,pad=0.3", fc="#ecfdf5", ec="#a7f3d0"),
)
# --- panel 2: sin(x)/x integrand
ax = axes[1]
x2 = np.linspace(1e-4, 20, 800)
ax.plot(x2, np.sin(x2) / x2, color="#2563eb", lw=2.0, label=r"$\sin(x)/x$")
ax.axhline(0, color="black", lw=0.6)
ax.set_xlabel("x")
ax.set_title(r"$\int_0^\infty \frac{\sin x}{x}\,dx = \pi/2$")
ax.legend()
ax.grid(True, alpha=0.3)
ax.annotate(
f"Numerical: {r2.value:.6f}\nExact: {exact2:.6f}",
xy=(0.97, 0.95), xycoords="axes fraction", ha="right", va="top",
fontsize=8, bbox=dict(boxstyle="round,pad=0.3", fc="#eff6ff", ec="#bfdbfe"),
)
# --- panel 3: partial integral convergence
ax = axes[2]
ax.plot(T_vals, partial1, color="#059669", lw=2.0, label=r"$\int_0^T e^{-x^2}$")
ax.axhline(exact1, color="#059669", lw=1.0, ls="--", alpha=0.6, label=f"exact={exact1:.5f}")
ax.plot(T_vals, partial2, color="#2563eb", lw=2.0, label=r"$\int_{0^+}^T \sin(x)/x$")
ax.axhline(exact2, color="#2563eb", lw=1.0, ls="--", alpha=0.6, label=f"exact={exact2:.5f}")
ax.set_xlabel("Upper limit T")
ax.set_ylabel("Partial integral")
ax.set_title("Convergence of partial integrals")
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
fig.suptitle("IMSL QDAGI: Improper Integrals", fontsize=13, fontweight="bold")
fig.tight_layout()
fig.savefig(plot_path, format="svg")
plt.close(fig)
return {"results": cases, "plot_path": str(plot_path)}
if __name__ == "__main__":
run_demo_imsl_improper_integrals()
Plot Output¶
Left: \(e^{-x^2}\) integrand with shaded area. Centre: \(\sin(x)/x\) integrand. Right: convergence of partial integrals \(\int_0^T\) as the upper limit \(T\) increases.¶
Console Output¶
IMSL QDAGI Example: Improper Integrals
==============================================================================
Integral Result Error Exact |Δ|
------------------------------------------------------------------------------
∫₀^∞ exp(-x²) dx 0.88622693 7.10e-09 0.88622693 0.00e+00
∫₀^∞ sin(x)/x dx 1.57079633 3.83e-09 1.57079633 1.18e-10
∫₀^1 ln(x)/√x dx -4.00000000 1.35e-13 -4.00000000 8.53e-14
==============================================================================