Recurrence Rules Example — Quadrature Rule Construction¶
This example demonstrates the full workflow for constructing Gauss quadrature rules from three-term recurrence coefficients and the Fejér Type-1 rule, including a round-trip validation:
Obtain Legendre recurrence coefficients via
recurrence_coefficients_classical()(IMSL RECCF).Build Gauss-Legendre rules (n = 8 and n = 16) via
gauss_quadrature_from_recurrence()(IMSL GQRCF) using the Golub–Welsch algorithm.Build a Fejér Type-1 rule (n = 8) via
fejer_quadrature_rule()(IMSL FQRUL).Round-trip: recover recurrence coefficients from the GL-8 rule via
recurrence_from_quadrature()(IMSL RECQR) and verify near-zero round-trip error.Compare accuracy of GL-8, Fejér-8, and GL-16 on
∫₋₁¹ eˣ dx = e − 1/e ≈ 2.35040.
Example Code¶
"""IMSL quadrature rule construction from recurrence coefficients.
Demonstrates classical recurrence coefficients (IMSL RECCF), Gauss quadrature
rule construction from recurrence (IMSL GQRCF), Fejér Type-1 rule (IMSL FQRUL),
and reverse round-trip from quadrature rule to recurrence (IMSL RECQR).
Workflow:
1. Obtain Legendre recurrence coefficients via
:func:`recurrence_coefficients_classical`.
2. Build Gauss-Legendre rules via :func:`gauss_quadrature_from_recurrence`.
3. Build Fejér Type-1 rules via :func:`fejer_quadrature_rule`.
4. Round-trip: recover recurrence from the GL rule via
:func:`recurrence_from_quadrature`.
5. Compare accuracy of GL(8), Fejér(8), and GL(16) on
``∫₋₁¹ eˣ dx = e − 1/e ≈ 2.35040``.
Outputs:
- Node/weight tables printed to stdout
- SVG plot saved to ``test_output/example_imsl_recurrence_rules.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 (
recurrence_coefficients_classical,
gauss_quadrature_from_recurrence,
recurrence_from_quadrature,
fejer_quadrature_rule,
)
def _apply_rule(rule, f) -> float: # type: ignore[no-untyped-def]
"""Apply a GaussQuadratureRule to integrate f on [-1, 1].
Args:
rule (GaussQuadratureRule): Quadrature rule with ``points`` and
``weights`` attributes.
f (Callable): Integrand function.
Returns:
float: Weighted sum approximating the integral.
"""
return float(np.dot(rule.weights, np.array([f(xi) for xi in rule.points])))
def run_demo_imsl_recurrence_rules() -> Dict[str, object]:
"""Run IMSL recurrence coefficient and quadrature rule construction examples.
Covers :func:`recurrence_coefficients_classical`,
:func:`gauss_quadrature_from_recurrence`, :func:`fejer_quadrature_rule`,
and :func:`recurrence_from_quadrature`.
Args:
None
Returns:
Dict[str, object]: Keys ``gl8``, ``fej8``, ``gl16`` (GaussQuadratureRule),
``alpha_rt``, ``beta_rt`` (round-trip recurrence), and
``plot_path`` (str).
"""
n_small = 8
n_large = 16
# ------------------------------------------------------------------
# 1. Classical Legendre recurrence coefficients
# ------------------------------------------------------------------
alpha8, beta8 = recurrence_coefficients_classical(n_small, weight_type="legendre")
# ------------------------------------------------------------------
# 2. Build GL-8 via Golub-Welsch
# ------------------------------------------------------------------
gl8 = gauss_quadrature_from_recurrence(alpha8, beta8)
# ------------------------------------------------------------------
# 3. GL-16 from recurrence
# ------------------------------------------------------------------
alpha16, beta16 = recurrence_coefficients_classical(n_large, weight_type="legendre")
gl16 = gauss_quadrature_from_recurrence(alpha16, beta16)
# ------------------------------------------------------------------
# 4. Fejér Type-1 rule (n=8)
# ------------------------------------------------------------------
fej8 = fejer_quadrature_rule(n_small)
# ------------------------------------------------------------------
# 5. Round-trip: rule → recurrence
# ------------------------------------------------------------------
alpha_rt, beta_rt = recurrence_from_quadrature(gl8.points, gl8.weights)
# ------------------------------------------------------------------
# 6. Compare accuracy on ∫₋₁¹ eˣ dx = e − 1/e
# ------------------------------------------------------------------
f_exp = math.exp
exact = math.exp(1.0) - math.exp(-1.0) # ≈ 2.35040...
val_gl8 = _apply_rule(gl8, f_exp)
val_fej8 = _apply_rule(fej8, f_exp)
val_gl16 = _apply_rule(gl16, f_exp)
err_gl8 = abs(val_gl8 - exact)
err_fej8 = abs(val_fej8 - exact)
err_gl16 = abs(val_gl16 - exact)
# ------------------------------------------------------------------
# Print results
# ------------------------------------------------------------------
print("\nIMSL Recurrence Coefficients & Quadrature Rule Construction")
print("=" * 70)
print(f"\nLegendre recurrence coefficients (n={n_small}):")
print(f" {'k':>3} {'alpha_k':>14} {'beta_k':>14}")
print(f" {'-'*3} {'-'*14} {'-'*14}")
for k in range(n_small):
print(f" {k:>3} {alpha8[k]:>14.8f} {beta8[k]:>14.8f}")
print(f"\n{n_small}-point Gauss-Legendre nodes and weights (from recurrence):")
print(f" {'Node':>16} {'Weight':>16}")
print(f" {'-'*16} {'-'*16}")
for xi, wi in zip(gl8.points, gl8.weights):
print(f" {xi:>16.10f} {wi:>16.10f}")
print(f"\n{n_small}-point Fejér Type-1 nodes and weights:")
print(f" {'Node':>16} {'Weight':>16}")
print(f" {'-'*16} {'-'*16}")
for xi, wi in zip(fej8.points, fej8.weights):
print(f" {xi:>16.10f} {wi:>16.10f}")
print(f"\nRound-trip recurrence from GL-{n_small} rule:")
print(f" Max |α_orig − α_rt| = {np.max(np.abs(alpha8 - alpha_rt)):.2e}")
print(f" Max |β_orig − β_rt| = {np.max(np.abs(beta8 - beta_rt)):.2e}")
print(f"\nAccuracy on ∫₋₁¹ eˣ dx = e − 1/e ≈ {exact:.10f}:")
print(f" {'Rule':<14} {'Computed':>14} {'|Error|':>12}")
print(f" {'-'*14} {'-'*14} {'-'*12}")
print(f" {'GL-8':<14} {val_gl8:>14.10f} {err_gl8:>12.3e}")
print(f" {'Fejer-8':<14} {val_fej8:>14.10f} {err_fej8:>12.3e}")
print(f" {'GL-16':<14} {val_gl16:>14.10f} {err_gl16:>12.3e}")
print("=" * 70)
# ------------------------------------------------------------------
# Plot
# ------------------------------------------------------------------
output_dir = Path("test_output")
output_dir.mkdir(parents=True, exist_ok=True)
plot_path = output_dir / "example_imsl_recurrence_rules.svg"
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# Panel 1: GL-8 nodes and Fejér-8 nodes on [-1, 1]
ax = axes[0]
y_gl = np.zeros(n_small) + 0.15
y_fej = np.zeros(n_small) - 0.15
ax.scatter(gl8.points, y_gl, s=60, color="#1d4ed8", zorder=5, label=f"GL-{n_small} nodes")
ax.scatter(fej8.points, y_fej, s=60, color="#dc2626", marker="^", zorder=5,
label=f"Fejér-{n_small} nodes")
ax.axhline(0.15, color="#1d4ed8", lw=0.5, alpha=0.4)
ax.axhline(-0.15, color="#dc2626", lw=0.5, alpha=0.4)
for xi in gl8.points:
ax.axvline(xi, color="#1d4ed8", lw=0.6, alpha=0.25)
for xi in fej8.points:
ax.axvline(xi, color="#dc2626", lw=0.6, alpha=0.25)
ax.set_xlim(-1.15, 1.15)
ax.set_ylim(-0.5, 0.5)
ax.set_yticks([])
ax.set_xlabel("Node position on [-1, 1]")
ax.set_title(f"Quadrature node locations\n(n = {n_small})")
ax.legend(fontsize=9, loc="upper left")
ax.grid(True, alpha=0.3, axis="x")
# Panel 2: weight comparison bar chart
ax = axes[1]
x_pos = np.arange(n_small)
width = 0.38
ax.bar(x_pos - width / 2, gl8.weights, width=width, color="#1d4ed8",
alpha=0.8, label=f"GL-{n_small} weights")
ax.bar(x_pos + width / 2, fej8.weights, width=width, color="#dc2626",
alpha=0.8, label=f"Fejér-{n_small} weights")
ax.set_xticks(x_pos)
ax.set_xticklabels([f"k={i+1}" for i in range(n_small)], fontsize=8)
ax.set_xlabel("Quadrature point index")
ax.set_ylabel("Weight")
ax.set_title(f"Quadrature weights comparison\n(n = {n_small})")
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, axis="y")
# Panel 3: error comparison bar chart (log scale)
ax = axes[2]
rule_labels = [f"GL-{n_small}", f"Fejér-{n_small}", f"GL-{n_large}"]
rule_errors = [err_gl8, err_fej8, err_gl16]
rule_colors = ["#1d4ed8", "#dc2626", "#059669"]
bars = ax.bar(rule_labels, rule_errors, color=rule_colors, alpha=0.85,
edgecolor="white", lw=0.8)
ax.set_yscale("log")
ax.set_ylabel(r"$|$computed $-$ exact$|$ (log scale)")
ax.set_title(
r"Accuracy on $\int_{-1}^{1} e^x\,dx = e - 1/e$"
+ f"\nExact ≈ {exact:.6f}"
)
for bar, err in zip(bars, rule_errors):
ax.text(
bar.get_x() + bar.get_width() / 2.0,
err * 3.0,
f"{err:.2e}",
ha="center", va="bottom", fontsize=9,
)
ax.grid(True, alpha=0.3, axis="y")
fig.suptitle(
"IMSL Recurrence Coefficients & Quadrature Rule Construction",
fontsize=13,
fontweight="bold",
)
fig.tight_layout()
fig.savefig(plot_path, format="svg")
plt.close(fig)
return {
"gl8": gl8,
"fej8": fej8,
"gl16": gl16,
"alpha_rt": alpha_rt,
"beta_rt": beta_rt,
"plot_path": str(plot_path),
}
if __name__ == "__main__":
run_demo_imsl_recurrence_rules()
Plot Output¶
Left: positions of GL-8 and Fejér-8 nodes on \([-1, 1]\). Centre: bar-chart comparison of weights for each quadrature point. Right: absolute-error bar chart (log scale) for GL-8, Fejér-8, and GL-16 applied to \(\int_{-1}^{1} e^x\,dx = e - 1/e\).¶
Console Output¶
E:\ITDS\itds\packages\integrationdifferentiation\examples\example_imsl_recurrence_rules.py:221: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all Axes decorations.
fig.tight_layout()
IMSL Recurrence Coefficients & Quadrature Rule Construction
======================================================================
Legendre recurrence coefficients (n=8):
k alpha_k beta_k
--- -------------- --------------
0 0.00000000 2.00000000
1 0.00000000 0.33333333
2 0.00000000 0.26666667
3 0.00000000 0.25714286
4 0.00000000 0.25396825
5 0.00000000 0.25252525
6 0.00000000 0.25174825
7 0.00000000 0.25128205
8-point Gauss-Legendre nodes and weights (from recurrence):
Node Weight
---------------- ----------------
-0.9602898565 0.1012285363
-0.7966664774 0.2223810345
-0.5255324099 0.3137066459
-0.1834346425 0.3626837834
0.1834346425 0.3626837834
0.5255324099 0.3137066459
0.7966664774 0.2223810345
0.9602898565 0.1012285363
8-point Fejér Type-1 nodes and weights:
Node Weight
---------------- ----------------
-0.9807852804 0.0669829457
-0.8314696123 0.2229879330
-0.5555702330 0.3241525191
-0.1950903220 0.3858766022
0.1950903220 0.3858766022
0.5555702330 0.3241525191
0.8314696123 0.2229879330
0.9807852804 0.0669829457
Round-trip recurrence from GL-8 rule:
Max |α_orig − α_rt| = 6.83e-16
Max |β_orig − β_rt| = 4.44e-16
Accuracy on ∫₋₁¹ eˣ dx = e − 1/e ≈ 2.3504023873:
Rule Computed |Error|
-------------- -------------- ------------
GL-8 2.3504023873 8.882e-16
Fejer-8 2.3504023937 6.367e-09
GL-16 2.3504023873 0.000e+00
======================================================================