Gauss Rules Example — Gauss Quadrature Rules

This example demonstrates Gauss quadrature rule generation (IMSL GQRUL):

• 5-point Gauss-Legendre rule: nodes and weights on [-1, 1]

• 5-point Gauss-Hermite rule: nodes and weights on (-∞, ∞)

• Use the 5-point GL rule to integrate exp(-x) on [-1, 1]

Example Code

"""IMSL GQRUL example: Gauss quadrature rules.

Demonstrates Gauss quadrature rule generation (IMSL GQRUL):

- 5-point Gauss-Legendre rule: nodes and weights on ``[-1, 1]``
- 5-point Gauss-Hermite rule: nodes and weights on ``(-∞, ∞)``
- Use the 5-point GL rule to integrate ``exp(-x)`` on ``[-1, 1]``

Outputs:
- Table printed to stdout
- SVG plot saved to ``test_output/demo_imsl_gauss_rules.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 gauss_quadrature_rule


def run_demo_imsl_gauss_rules() -> Dict[str, object]:
    """Run IMSL GQRUL Gauss quadrature rule examples.

    Args:
        None

    Returns:
        Dict[str, object]: Result dict with keys ``gl_rule``, ``gh_rule``
            (GaussQuadratureRule) and ``plot_path`` (str).
    """
    n = 5

    gl = gauss_quadrature_rule(n, weight_type="legendre")
    gh = gauss_quadrature_rule(n, weight_type="hermite")

    # Use GL rule to integrate exp(-x) on [-1, 1]
    # ∫_{-1}^{1} exp(-x) dx = e - 1/e ≈ 2.35040
    f = lambda x: np.exp(-x)  # noqa: E731
    gl_estimate = float(np.dot(gl.weights, np.array([f(xi) for xi in gl.points])))
    exact = float(np.exp(1.0) - np.exp(-1.0))

    print("\nIMSL GQRUL Example: Gauss Quadrature Rules")
    print("=" * 60)
    print(f"\n5-point Gauss-Legendre rule on [-1, 1]:")
    print(f"  {'Node':>14}  {'Weight':>14}")
    print(f"  {'-'*14}  {'-'*14}")
    for xi, wi in zip(gl.points, gl.weights):
        print(f"  {xi:>14.8f}  {wi:>14.8f}")
    print(f"\n  ∫_{{-1}}^{{1}} exp(-x) dx via GL rule : {gl_estimate:.8f}")
    print(f"  Exact (e - 1/e)                  : {exact:.8f}")
    print(f"  Error                            : {abs(gl_estimate - exact):.2e}")

    print(f"\n5-point Gauss-Hermite rule on (-∞, ∞):")
    print(f"  {'Node':>14}  {'Weight':>14}")
    print(f"  {'-'*14}  {'-'*14}")
    for xi, wi in zip(gh.points, gh.weights):
        print(f"  {xi:>14.8f}  {wi:>14.8f}")
    print(f"\n  (Weights include exp(-x²) factor; sum of weights ≈ √π)")
    print(f"  Sum of GH weights: {gh.weights.sum():.8f}  √π ≈ {np.sqrt(np.pi):.8f}")
    print("=" * 60)

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

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

    # Panel 1: Gauss-Legendre
    ax = axes[0]
    ax.bar(gl.points, gl.weights, width=0.04, color="#059669", alpha=0.85,
           label="GL weights", zorder=3)
    ax.scatter(gl.points, np.zeros(n), color="#047857", s=60, zorder=5,
               label="GL nodes", marker="|", linewidths=2)
    for xi, wi in zip(gl.points, gl.weights):
        ax.text(xi, wi + 0.005, f"{xi:.3f}", ha="center", va="bottom", fontsize=7.5)
    ax.set_xlabel("Node position")
    ax.set_ylabel("Weight")
    ax.set_title(f"{n}-point Gauss-Legendre Rule\non [-1, 1]")
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)
    ax.set_xlim(-1.3, 1.3)

    # Panel 2: Gauss-Hermite
    ax = axes[1]
    ax.bar(gh.points, gh.weights, width=0.15, color="#10b981", alpha=0.85,
           label="GH weights", zorder=3)
    ax.scatter(gh.points, np.zeros(n), color="#065f46", s=60, zorder=5,
               label="GH nodes", marker="|", linewidths=2)
    for xi, wi in zip(gh.points, gh.weights):
        ax.text(xi, wi + 0.005, f"{xi:.3f}", ha="center", va="bottom", fontsize=7.5)
    ax.set_xlabel("Node position")
    ax.set_ylabel("Weight")
    ax.set_title(f"{n}-point Gauss-Hermite Rule\non (-∞, ∞)")
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)

    fig.suptitle("IMSL GQRUL: Gauss Quadrature Rules", fontsize=13, fontweight="bold")
    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)

    return {"gl_rule": gl, "gh_rule": gh, "plot_path": str(plot_path)}


if __name__ == "__main__":
    run_demo_imsl_gauss_rules()

Plot Output

Left: 5-point Gauss-Legendre nodes and weights on \([-1, 1]\). Right: 5-point Gauss-Hermite nodes and weights on \((-\infty, \infty)\). Node positions are labelled above each bar.

Console Output

IMSL GQRUL Example: Gauss Quadrature Rules
============================================================

5-point Gauss-Legendre rule on [-1, 1]:
            Node          Weight
  --------------  --------------
     -0.90617985      0.23692689
     -0.53846931      0.47862867
      0.00000000      0.56888889
      0.53846931      0.47862867
      0.90617985      0.23692689

  ∫_{-1}^{1} exp(-x) dx via GL rule : 2.35040239
  Exact (e - 1/e)                  : 2.35040239
  Error                            : 8.25e-10

5-point Gauss-Hermite rule on (-∞, ∞):
            Node          Weight
  --------------  --------------
     -2.02018287      0.01995324
     -0.95857246      0.39361932
      0.00000000      0.94530872
      0.95857246      0.39361932
      2.02018287      0.01995324

  (Weights include exp(-x²) factor; sum of weights ≈ √π)
  Sum of GH weights: 1.77245385  √π ≈ 1.77245385
============================================================