TWODQ Example — 2D Iterated Quadrature

This example demonstrates 2D numerical integration (IMSL TWODQ):

• Integrate f(x, y) = exp(-(x² + y²)) over [0,1]×[0,1] → ≈ 0.55777

• Integrate f(x, y) = x · y over [0,1]×[0,1] → 0.25

Example Code

"""IMSL TWODQ example: 2D iterated quadrature.

Demonstrates 2D numerical integration (IMSL TWODQ):

- Integrate ``f(x, y) = exp(-(x² + y²))`` over ``[0,1]×[0,1]``
  → ``(√π/2 · erf(1))² ≈ 0.55777``
- Integrate ``f(x, y) = x · y`` over ``[0,1]×[0,1]`` → ``0.25``

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

from __future__ import annotations

from pathlib import Path
from typing import Dict

import matplotlib.pyplot as plt
import numpy as np
from scipy.special import erf

from integrationdifferentiation import integrate_2d


def run_demo_imsl_twodq() -> Dict[str, object]:
    """Run IMSL TWODQ 2D quadrature examples.

    Args:
        None

    Returns:
        Dict[str, object]: Result dict with keys ``r1``, ``r2``
            (IntegrationResult) and ``plot_path`` (str).
    """
    # Case 1: exp(-(x² + y²)) over [0,1]×[0,1]
    def f1(y: float, x: float) -> float:
        return float(np.exp(-(x ** 2 + y ** 2)))

    r1 = integrate_2d(f1, 0.0, 1.0, lambda x: 0.0, lambda x: 1.0)
    exact1 = float((np.sqrt(np.pi) / 2.0 * erf(1.0)) ** 2)

    # Case 2: x * y over [0,1]×[0,1] → exact = 0.25
    def f2(y: float, x: float) -> float:
        return float(x * y)

    r2 = integrate_2d(f2, 0.0, 1.0, lambda x: 0.0, lambda x: 1.0)
    exact2 = 0.25

    print("\nIMSL TWODQ Example: 2D Iterated Quadrature")
    print("=" * 72)
    print(f"{'Integrand':<30} {'Result':>12} {'Error':>12} {'Exact':>12}")
    print("-" * 72)
    print(
        f"{'∬exp(-(x²+y²)) dxdy, [0,1]²':<30} {r1.value:>12.8f}"
        f" {r1.error:>12.2e} {exact1:>12.8f}"
    )
    print(
        f"{'∬xy dxdy, [0,1]²':<30} {r2.value:>12.8f}"
        f" {r2.error:>12.2e} {exact2:>12.8f}"
    )
    print("=" * 72)

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

    x_vals = np.linspace(0, 1, 60)
    y_vals = np.linspace(0, 1, 60)
    X, Y = np.meshgrid(x_vals, y_vals)

    fig, axes = plt.subplots(1, 2, figsize=(12, 5), subplot_kw={"projection": "3d"})

    # Panel 1: exp(-(x² + y²))
    Z1 = np.exp(-(X ** 2 + Y ** 2))
    ax = axes[0]
    ax.plot_surface(X, Y, Z1, cmap="Greens", alpha=0.85, edgecolor="none")
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_zlabel("f(x,y)")
    ax.set_title(
        f"$f = e^{{-(x^2+y^2)}}$\n"
        f"Result: {r1.value:.6f}  Exact: {exact1:.6f}"
    )

    # Panel 2: x * y
    Z2 = X * Y
    ax = axes[1]
    ax.plot_surface(X, Y, Z2, cmap="YlGn", alpha=0.85, edgecolor="none")
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_zlabel("f(x,y)")
    ax.set_title(
        f"$f = x \\cdot y$\n"
        f"Result: {r2.value:.6f}  Exact: {exact2:.6f}"
    )

    fig.suptitle("IMSL TWODQ: 2D Iterated Quadrature", fontsize=13, fontweight="bold")
    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)

    return {"r1": r1, "r2": r2, "plot_path": str(plot_path)}


if __name__ == "__main__":
    run_demo_imsl_twodq()

Plot Output

Left: \(f(x,y) = e^{-(x^2+y^2)}\) over \([0,1]^2\). Right: \(f(x,y) = x \cdot y\) over \([0,1]^2\). Titles show computed and exact values.

Console Output

IMSL TWODQ Example: 2D Iterated Quadrature
========================================================================
Integrand                            Result        Error        Exact
------------------------------------------------------------------------
∬exp(-(x²+y²)) dxdy, [0,1]²      0.55774629     8.29e-15   0.55774629
∬xy dxdy, [0,1]²                 0.25000000     5.54e-15   0.25000000
========================================================================