Singular Integrals Example — Adaptive Singular-Point & Cauchy PV Integration¶
This example demonstrates adaptive integration with interior singularity points (IMSL QDAGP), Cauchy principal-value integration (IMSL QDAWC), and 2D quadrature over a rectangular domain with integrable edge singularities (IMSL QDAG2D):
∫₋₁¹ ln|x| dx = -2(log singularity at interior point x = 0) viaintegrate_adaptive_singular_points()PV ∫₋₁¹ 1/x dx = 0(Cauchy principal value, f = 1, c = 0) viaintegrate_cauchy_principal_value()∫₀¹ ∫₀¹ 1/√(x·y) dy dx = 4(separable, integrable edge singularities) viaintegrate_2d_singular()
Example Code¶
"""IMSL singular and principal-value integration examples.
Demonstrates adaptive integration with interior singularity points (IMSL QDAGP),
Cauchy principal-value integration (IMSL QDAWC), and 2D quadrature over a
rectangular domain with integrable endpoint singularities (IMSL QDAG2D).
Integrals demonstrated:
- ``∫₋₁¹ ln|x| dx = -2`` (log singularity at interior point x = 0)
via :func:`integrate_adaptive_singular_points`
- ``PV ∫₋₁¹ 1/x dx = 0`` (Cauchy principal value)
via :func:`integrate_cauchy_principal_value`
- ``∫₀¹ ∫₀¹ 1/√(x·y) dy dx = 4`` (integrable edge singularities)
via :func:`integrate_2d_singular`
Outputs:
- Table printed to stdout
- SVG plot saved to ``test_output/example_imsl_singular_integrals.svg``
"""
from __future__ import annotations
import math
from pathlib import Path
from typing import Dict
import matplotlib.pyplot as plt
import numpy as np
from integrationdifferentiation import (
integrate_adaptive_singular_points,
integrate_cauchy_principal_value,
integrate_2d_singular,
)
def run_demo_imsl_singular_integrals() -> Dict[str, object]:
"""Run IMSL singular and principal-value integration examples.
Covers :func:`integrate_adaptive_singular_points`,
:func:`integrate_cauchy_principal_value`, and
:func:`integrate_2d_singular`.
Args:
None
Returns:
Dict[str, object]: Keys ``r1`` through ``r3`` (IntegrationResult)
and ``plot_path`` (str).
"""
# ------------------------------------------------------------------
# 1. Adaptive integration with interior singularity
# ∫₋₁¹ ln|x| dx — log singularity at interior point x = 0
# Exact: 2 * ∫₀¹ ln(x) dx = 2 * (-1) = -2
# ------------------------------------------------------------------
f_log = lambda x: math.log(abs(x)) if x != 0.0 else -1e15 # noqa: E731
r1 = integrate_adaptive_singular_points(
f_log, a=-1.0, b=1.0, points=[0.0], abs_tol=1e-9
)
exact1 = -2.0
# ------------------------------------------------------------------
# 2. Cauchy principal value
# PV ∫₋₁¹ 1/x dx = 0 (f(x) = 1, singularity at c = 0)
# The function computes ∫ f(x)/(x-c) dx → PV ∫ 1/x dx
# ------------------------------------------------------------------
f_one = lambda x: 1.0 # noqa: E731
r2 = integrate_cauchy_principal_value(
f_one, a=-1.0, b=1.0, c=0.0, abs_tol=1e-9
)
exact2 = 0.0
# ------------------------------------------------------------------
# 3. 2D quadrature over [0,1]² with edge singularities
# ∫₀¹ ∫₀¹ 1/√(x·y) dy dx = (∫₀¹ x⁻¹/² dx)² = 2² = 4
# f(y, x) — inner variable y is the first argument
# ------------------------------------------------------------------
f_2d = lambda y, x: 1.0 / math.sqrt(x * y) if (x > 0 and y > 0) else 0.0 # noqa: E731
r3 = integrate_2d_singular(
f_2d, ax=0.0, bx=1.0, ay=0.0, by=1.0, abs_tol=1e-6
)
exact3 = 4.0
# ------------------------------------------------------------------
# Print results
# ------------------------------------------------------------------
print("\nIMSL Singular & Principal-Value Integration Examples")
print("=" * 65)
print("\n1. Interior log singularity — integrate_adaptive_singular_points")
print(f" Integral : ∫₋₁¹ ln|x| dx")
print(f" Computed : {r1.value:+.10f} Error est: {r1.error:.2e}")
print(f" Exact : {exact1:+.10f}")
print(f" Abs diff : {abs(r1.value - exact1):.2e}")
print("\n2. Cauchy principal value — integrate_cauchy_principal_value")
print(f" Integral : PV ∫₋₁¹ 1/x dx")
print(f" Computed : {r2.value:+.10f} Error est: {r2.error:.2e}")
print(f" Exact : {exact2:+.10f}")
print(f" Abs diff : {abs(r2.value - exact2):.2e}")
print("\n3. 2D with edge singularities — integrate_2d_singular")
print(f" Integral : ∫₀¹∫₀¹ 1/√(x·y) dy dx")
print(f" Computed : {r3.value:+.10f} Error est: {r3.error:.2e}")
print(f" Exact : {exact3:+.10f}")
print(f" Abs diff : {abs(r3.value - exact3):.2e}")
print("=" * 65)
# ------------------------------------------------------------------
# Plot
# ------------------------------------------------------------------
output_dir = Path("test_output")
output_dir.mkdir(parents=True, exist_ok=True)
plot_path = output_dir / "example_imsl_singular_integrals.svg"
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
# Panel 1: ln|x| on [-1, 1] with singularity marked
ax = axes[0]
x_neg = np.linspace(-1.0, -0.02, 300)
x_pos = np.linspace(0.02, 1.0, 300)
ax.plot(x_neg, np.log(np.abs(x_neg)), color="#1d4ed8", lw=2.0, label=r"$\ln|x|$")
ax.plot(x_pos, np.log(np.abs(x_pos)), color="#1d4ed8", lw=2.0)
ax.fill_between(x_neg, np.log(np.abs(x_neg)), 0, alpha=0.15, color="#1d4ed8")
ax.fill_between(x_pos, np.log(np.abs(x_pos)), 0, alpha=0.15, color="#1d4ed8")
ax.axvline(0.0, color="#ef4444", lw=1.5, ls="--", label="singularity x=0")
ax.set_xlabel("x")
ax.set_ylabel(r"$\ln|x|$")
ax.set_title(
r"$\int_{-1}^{1} \ln|x|\,dx = -2$" + f"\nComputed: {r1.value:.6f}"
)
ax.set_ylim(-5, 0.5)
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
# Panel 2: 1/x with PV region
ax = axes[1]
ax.plot(x_neg, 1.0 / x_neg, color="#7c3aed", lw=2.0, label=r"$f(x) = 1$")
ax.plot(x_pos, 1.0 / x_pos, color="#7c3aed", lw=2.0)
ax.fill_between(x_neg, 1.0 / x_neg, 0, alpha=0.15, color="#7c3aed")
ax.fill_between(x_pos, 1.0 / x_pos, 0, alpha=0.15, color="#db2777",
label="areas cancel (PV = 0)")
ax.axvline(0.0, color="#ef4444", lw=1.5, ls="--", label="Cauchy pt c=0")
ax.set_xlabel("x")
ax.set_ylabel(r"$1/(x - c)$ with $c = 0$")
ax.set_title(
r"$\mathrm{PV}\int_{-1}^{1}\frac{1}{x}\,dx = 0$"
+ f"\nComputed: {r2.value:.2e}"
)
ax.set_ylim(-20, 20)
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
# Panel 3: 2D integrand colourmap
ax = axes[2]
xv = np.linspace(0.02, 1.0, 200)
yv = np.linspace(0.02, 1.0, 200)
X, Y = np.meshgrid(xv, yv)
Z = 1.0 / np.sqrt(X * Y)
c = ax.pcolormesh(X, Y, Z, cmap="magma_r", shading="auto", vmax=20)
fig.colorbar(c, ax=ax, fraction=0.046, pad=0.04)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title(
r"$\int_0^1\!\int_0^1 \frac{1}{\sqrt{xy}}\,dy\,dx = 4$"
+ f"\nComputed: {r3.value:.6f}"
)
ax.set_aspect("equal")
fig.suptitle(
"IMSL Singular & Principal-Value Integration",
fontsize=13,
fontweight="bold",
)
fig.tight_layout()
fig.savefig(plot_path, format="svg")
plt.close(fig)
return {"r1": r1, "r2": r2, "r3": r3, "plot_path": str(plot_path)}
if __name__ == "__main__":
run_demo_imsl_singular_integrals()
Plot Output¶
Left: \(\int_{-1}^{1} \ln|x|\,dx = -2\) with log singularity marked at \(x = 0\). Centre: Cauchy principal-value integral \(\mathrm{PV}\int_{-1}^{1}\frac{1}{x}\,dx = 0\) showing cancellation of positive and negative areas. Right: colourmap of \(1/\sqrt{x y}\) over \([0,1]^2\) with the computed result annotated.¶
Console Output¶
IMSL Singular & Principal-Value Integration Examples
=================================================================
1. Interior log singularity ΓÇö integrate_adaptive_singular_points
Integral : ∫₋₁¹ ln|x| dx
Computed : -2.0000000000 Error est: 3.33e-15
Exact : -2.0000000000
Abs diff : 8.88e-16
2. Cauchy principal value ΓÇö integrate_cauchy_principal_value
Integral : PV ∫₋₁¹ 1/x dx
Computed : +0.0000000000 Error est: 7.27e-13
Exact : +0.0000000000
Abs diff : 2.22e-16
3. 2D with edge singularities ΓÇö integrate_2d_singular
Integral : ∫₀¹∫₀¹ 1/√(x·y) dy dx
Computed : +4.0000000000 Error est: 1.92e-12
Exact : +4.0000000000
Abs diff : 1.78e-15
=================================================================