Weighted Integrals Example — Oscillatory, Algebraic, 3D & ND Integration¶
This example demonstrates weighted oscillatory quadrature (IMSL QDAWO), algebraic-weighted quadrature (IMSL QDAWS), adaptive 1D quadrature (IMSL QDAG1D), 3D rectangular quadrature with a corner singularity (IMSL QDAG3D), and N-dimensional integration (IMSL QAND):
∫₀^π x·sin(x) dx = π(oscillatory weight sin(x), ω = 1) viaintegrate_weighted_oscillatory()∫₀¹ √(x(1-x)) dx = π/8(algebraic weight, α = β = 0.5) viaintegrate_weighted_algebraic()∫₀¹ x·eˣ dx = 1(smooth 1D integrand) viaintegrate_adaptive_1d()∫₀¹∫₀¹∫₀¹ (xyz)^(-1/4) dz dy dx = (4/3)³(separable 3D singularity) viaintegrate_3d_singular()∫∫∫ e^(-x²-y²-z²) dx dy dz = π^(3/2)(3D Gaussian over [-4, 4]³) viaintegrate_nd()
Example Code¶
"""IMSL weighted and multidimensional integration examples.
Demonstrates weighted oscillatory quadrature (IMSL QDAWO), algebraic-weighted
quadrature (IMSL QDAWS), adaptive 1D quadrature (IMSL QDAG1D), 3D rectangular
quadrature (IMSL QDAG3D), and N-dimensional quadrature (IMSL QAND).
Integrals demonstrated:
- ``∫₀^π x·sin(x) dx = π`` (oscillatory weight sin(x))
via :func:`integrate_weighted_oscillatory`
- ``∫₀¹ √(x(1-x)) dx = π/8`` (algebraic weight x^0.5·(1-x)^0.5)
via :func:`integrate_weighted_algebraic`
- ``∫₀¹ x·eˣ dx = 1`` (smooth 1D integrand)
via :func:`integrate_adaptive_1d`
- ``∫₀¹∫₀¹∫₀¹ (xyz)^(-1/4) dz dy dx = (4/3)³`` (3D integrable singularity)
via :func:`integrate_3d_singular`
- ``∫∫∫ e^(-x²-y²-z²) dx dy dz = π^(3/2)`` (3D Gaussian, N-dimensional)
via :func:`integrate_nd`
Outputs:
- Table printed to stdout
- SVG plot saved to ``test_output/example_imsl_weighted_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_1d,
integrate_weighted_oscillatory,
integrate_weighted_algebraic,
integrate_3d_singular,
integrate_nd,
)
def run_demo_imsl_weighted_integrals() -> Dict[str, object]:
"""Run IMSL weighted and multidimensional integration examples.
Covers :func:`integrate_weighted_oscillatory`,
:func:`integrate_weighted_algebraic`, :func:`integrate_adaptive_1d`,
:func:`integrate_3d_singular`, and :func:`integrate_nd`.
Args:
None
Returns:
Dict[str, object]: Keys ``r1`` through ``r5`` (IntegrationResult)
and ``plot_path`` (str).
"""
# ------------------------------------------------------------------
# 1. Oscillatory weight: ∫₀^π x·sin(x) dx = π
# f(x) = x, omega = 1, weight = 'sin'
# Exact: [-x cos(x) + sin(x)]₀^π = π
# ------------------------------------------------------------------
f_osc = lambda x: x # noqa: E731
r1 = integrate_weighted_oscillatory(
f_osc, a=0.0, b=math.pi, omega=1.0, weight="sin", abs_tol=1e-10
)
exact1 = math.pi
# ------------------------------------------------------------------
# 2. Algebraic weight: ∫₀¹ √(x(1-x)) dx = π/8
# f(x) = 1, alpha = 0.5, beta = 0.5, weight_type = 1
# weight = (x-0)^0.5 · (1-x)^0.5 = √(x(1-x))
# Exact: B(3/2, 3/2) = Γ(3/2)²/Γ(3) = (√π/2)²/2 = π/8
# ------------------------------------------------------------------
f_alg = lambda x: 1.0 # noqa: E731
r2 = integrate_weighted_algebraic(
f_alg, a=0.0, b=1.0, alpha=0.5, beta=0.5, weight_type=1, abs_tol=1e-10
)
exact2 = math.pi / 8.0
# ------------------------------------------------------------------
# 3. Adaptive 1D: ∫₀¹ x·eˣ dx = 1
# Exact: [xeˣ - eˣ]₀¹ = (e - e) - (0 - 1) = 1
# ------------------------------------------------------------------
f_1d = lambda x: x * math.exp(x) # noqa: E731
r3 = integrate_adaptive_1d(f_1d, a=0.0, b=1.0, abs_tol=1e-12)
exact3 = 1.0
# ------------------------------------------------------------------
# 4. 3D rectangular with integrable singularity
# ∫₀¹∫₀¹∫₀¹ (xyz)^(-1/4) dz dy dx = (∫₀¹ t^(-1/4) dt)³ = (4/3)³
# f(z, y, x) — z is innermost argument
# ------------------------------------------------------------------
def f_3d(z: float, y: float, x: float) -> float:
"""Separable integrand with corner singularity.
Args:
z (float): Innermost integration variable.
y (float): Middle integration variable.
x (float): Outer integration variable.
Returns:
float: (x·y·z)^(-1/4) or 0 at the origin.
"""
if x <= 0.0 or y <= 0.0 or z <= 0.0:
return 0.0
return (x * y * z) ** (-0.25)
r4 = integrate_3d_singular(
f_3d,
ax=0.0, bx=1.0,
ay=0.0, by=1.0,
az=0.0, bz=1.0,
abs_tol=1e-6,
)
exact4 = (4.0 / 3.0) ** 3 # ≈ 2.3704
# ------------------------------------------------------------------
# 5. N-dimensional Gaussian: ∫∫∫ e^(-x²-y²-z²) = π^(3/2)
# via integrate_nd over [-4, 4]³ (truncation error < 1e-14)
# ------------------------------------------------------------------
def f_nd(*args: float) -> float:
"""3D Gaussian integrand.
Args:
*args (float): Coordinates (z, y, x) — innermost dimension first.
Returns:
float: e^(-(x²+y²+z²)).
"""
return math.exp(-sum(xi * xi for xi in args))
r5 = integrate_nd(
f_nd, ranges=[(-4.0, 4.0), (-4.0, 4.0), (-4.0, 4.0)], abs_tol=1e-6
)
exact5 = math.pi ** 1.5 # ≈ 5.5683
# ------------------------------------------------------------------
# Print results
# ------------------------------------------------------------------
print("\nIMSL Weighted & Multidimensional Integration Examples")
print("=" * 65)
print("\n1. Oscillatory weight — integrate_weighted_oscillatory")
print(f" Integral : ∫₀^π x·sin(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. Algebraic weight — integrate_weighted_algebraic")
print(f" Integral : ∫₀¹ √(x(1-x)) dx")
print(f" Computed : {r2.value:+.10f} Error est: {r2.error:.2e}")
print(f" Exact : {exact2:+.10f} (π/8)")
print(f" Abs diff : {abs(r2.value - exact2):.2e}")
print("\n3. Adaptive 1D — integrate_adaptive_1d")
print(f" Integral : ∫₀¹ x·eˣ dx")
print(f" Computed : {r3.value:+.10f} Error est: {r3.error:.2e}")
print(f" Exact : {exact3:+.10f} (= 1)")
print(f" Abs diff : {abs(r3.value - exact3):.2e}")
print("\n4. 3D singular rectangular — integrate_3d_singular")
print(f" Integral : ∫₀¹∫₀¹∫₀¹ (xyz)^(-1/4) dz dy dx")
print(f" Computed : {r4.value:+.10f} Error est: {r4.error:.2e}")
print(f" Exact : {exact4:+.10f} ((4/3)³)")
print(f" Abs diff : {abs(r4.value - exact4):.2e}")
print("\n5. N-dimensional Gaussian — integrate_nd (3D)")
print(f" Integral : ∫∫∫ exp(-x²-y²-z²) dx dy dz over [-4,4]³")
print(f" Computed : {r5.value:+.10f} Error est: {r5.error:.2e}")
print(f" Exact : {exact5:+.10f} (π^(3/2))")
print(f" Abs diff : {abs(r5.value - exact5):.2e}")
print("=" * 65)
# ------------------------------------------------------------------
# Plot
# ------------------------------------------------------------------
output_dir = Path("test_output")
output_dir.mkdir(parents=True, exist_ok=True)
plot_path = output_dir / "example_imsl_weighted_integrals.svg"
fig, axes = plt.subplots(2, 2, figsize=(13, 10))
# Panel 1: oscillatory integrand x·sin(x) on [0,π]
ax = axes[0, 0]
xv = np.linspace(0.0, math.pi, 500)
y_osc = xv * np.sin(xv)
ax.plot(xv, y_osc, color="#0f766e", lw=2.2, label=r"$x\,\sin(x)$")
ax.fill_between(xv, y_osc, 0, alpha=0.18, color="#0f766e")
ax.set_xlabel("x")
ax.set_ylabel(r"$x\,\sin(x)$")
ax.set_title(
r"$\int_0^\pi x\sin(x)\,dx = \pi$"
+ f"\nComputed: {r1.value:.8f}"
)
ax.axhline(0, color="gray", lw=0.7)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# Panel 2: algebraic-weighted integrand √(x(1-x)) on [0,1]
ax = axes[0, 1]
xv2 = np.linspace(0.001, 0.999, 500)
y_alg = np.sqrt(xv2 * (1.0 - xv2))
ax.plot(xv2, y_alg, color="#7c3aed", lw=2.2, label=r"$\sqrt{x(1-x)}$")
ax.fill_between(xv2, y_alg, 0, alpha=0.18, color="#7c3aed")
ax.set_xlabel("x")
ax.set_ylabel(r"$\sqrt{x(1-x)}$")
ax.set_title(
r"$\int_0^1 \sqrt{x(1-x)}\,dx = \pi/8$"
+ f"\nComputed: {r2.value:.8f}"
)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# Panel 3: adaptive 1D x·eˣ
ax = axes[1, 0]
xv3 = np.linspace(0.0, 1.0, 300)
y_1d = xv3 * np.exp(xv3)
ax.plot(xv3, y_1d, color="#b45309", lw=2.2, label=r"$x\,e^x$")
ax.fill_between(xv3, y_1d, 0, alpha=0.18, color="#b45309")
ax.set_xlabel("x")
ax.set_ylabel(r"$x\,e^x$")
ax.set_title(
r"$\int_0^1 x\,e^x\,dx = 1$"
+ f"\nComputed: {r3.value:.10f}"
)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# Panel 4: error summary bar chart
ax = axes[1, 1]
labels = [
"Oscillatory\n(1D)",
"Algebraic wt\n(1D)",
"Adaptive 1D",
"3D singular",
"ND Gaussian\n(3D)",
]
errors = [
abs(r1.value - exact1),
abs(r2.value - exact2),
abs(r3.value - exact3),
abs(r4.value - exact4),
abs(r5.value - exact5),
]
colors = ["#0f766e", "#7c3aed", "#b45309", "#dc2626", "#1d4ed8"]
ax.bar(labels, errors, color=colors, alpha=0.82, edgecolor="white", lw=0.8)
ax.set_yscale("log")
ax.set_ylabel("Absolute error |computed − exact|")
ax.set_title("Accuracy comparison across all examples")
ax.grid(True, alpha=0.3, axis="y")
for i, e in enumerate(errors):
ax.text(i, e * 2.0, f"{e:.1e}", ha="center", va="bottom", fontsize=7.5)
fig.suptitle(
"IMSL Weighted & Multidimensional Integration",
fontsize=13,
fontweight="bold",
)
fig.tight_layout()
fig.savefig(plot_path, format="svg")
plt.close(fig)
return {
"r1": r1, "r2": r2, "r3": r3, "r4": r4, "r5": r5,
"plot_path": str(plot_path),
}
if __name__ == "__main__":
run_demo_imsl_weighted_integrals()
Plot Output¶
Top-left: \(\int_0^\pi x\sin(x)\,dx = \pi\) — oscillatory-weight integrand. Top-right: \(\int_0^1 \sqrt{x(1-x)}\,dx = \pi/8\) — algebraic-weighted integrand (Beta function). Bottom-left: \(\int_0^1 x\,e^x\,dx = 1\) — smooth 1D adaptive quadrature. Bottom-right: absolute-error comparison across all five examples on a log scale.¶
Console Output¶
E:\ITDS\itds\packages\integrationdifferentiation\examples\example_imsl_weighted_integrals.py:260: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all Axes decorations.
fig.tight_layout()
IMSL Weighted & Multidimensional Integration Examples
=================================================================
1. Oscillatory weight ΓÇö integrate_weighted_oscillatory
Integral : ∫₀^π x·sin(x) dx
Computed : +3.1415926536 Error est: 3.49e-14
Exact : +3.1415926536 (π)
Abs diff : 0.00e+00
2. Algebraic weight ΓÇö integrate_weighted_algebraic
Integral : ∫₀¹ √(x(1-x)) dx
Computed : +0.3926990817 Error est: 1.44e-15
Exact : +0.3926990817 (π/8)
Abs diff : 5.55e-17
3. Adaptive 1D ΓÇö integrate_adaptive_1d
Integral : ∫₀¹ x·eˣ dx
Computed : +1.0000000000 Error est: 1.11e-14
Exact : +1.0000000000 (= 1)
Abs diff : 0.00e+00
4. 3D singular rectangular ΓÇö integrate_3d_singular
Integral : ∫₀¹∫₀¹∫₀¹ (xyz)^(-1/4) dz dy dx
Computed : +2.3703703704 Error est: 4.69e-13
Exact : +2.3703703704 ((4/3)┬│)
Abs diff : 0.00e+00
5. N-dimensional Gaussian ΓÇö integrate_nd (3D)
Integral : Γê½Γê½Γê½ exp(-x┬▓-y┬▓-z┬▓) dx dy dz over [-4,4]┬│
Computed : +5.5683277393 Error est: 9.50e-07
Exact : +5.5683279968 (π^(3/2))
Abs diff : 2.58e-07
=================================================================