diff --git a/ROADMAP.md b/ROADMAP.md index 554d380..4143c31 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -74,32 +74,40 @@ The workhorse feature. Most real integration problems need error-driven refineme Complete the Gaussian quadrature family. -### 2.1 Gauss-Jacobi +### 2.1 Gauss-Jacobi ✅ - Weight function (1-x)^alpha * (1+x)^beta on [-1, 1] - Generalises Legendre (alpha=beta=0), Chebyshev (alpha=beta=-0.5), Gegenbauer -- Newton's method with three-term recurrence, asymptotics for large n +- Golub-Welsch algorithm: eigenvalues of symmetric tridiagonal Jacobi matrix +- Handles all valid parameters including the singular Chebyshev limit (alpha+beta=-1) -### 2.2 Gauss-Laguerre (generalised) +### 2.2 Gauss-Laguerre (generalised) ✅ - Weight function x^alpha * e^(-x) on [0, inf) -- GLR algorithm for moderate n, RH asymptotics for large n +- Golub-Welsch algorithm via Laguerre three-term recurrence coefficients -### 2.3 Gauss-Hermite +### 2.3 Gauss-Hermite ✅ - Weight function e^(-x^2) on (-inf, inf) -- Three-term recurrence for moderate n, uniform asymptotics for large n +- Golub-Welsch algorithm via physicists' Hermite recurrence coefficients -### 2.4 Gauss-Chebyshev (types I and II) +### 2.4 Gauss-Chebyshev (types I and II) ✅ - Closed-form nodes and weights (no iteration needed) - Useful as comparison/baseline and for Chebyshev interpolation -### 2.5 Gauss-Radau and Gauss-Lobatto -- Include one or both endpoints in the node set -- Important for spectral methods and ODE solvers +### 2.5 Gauss-Radau and Gauss-Lobatto ✅ +- Gauss-Radau: Golub-Welsch with Radau modification of Legendre Jacobi matrix +- Gauss-Lobatto: Newton on P'_{n-1}(x) using Legendre ODE second derivative +- Left/right variants for Radau, both endpoints for Lobatto -### 2.6 Clenshaw-Curtis -- Nodes at Chebyshev points, weights via DCT +### 2.6 Clenshaw-Curtis ✅ +- Nodes at Chebyshev extrema cos(kπ/(n-1)), weights via explicit DCT formula - Nested: degree doubling reuses previous evaluations - Competitive with Gauss for smooth functions, better for adaptive refinement +### 2.7 Golub-Welsch eigenvalue solver ✅ +- Shared `pub(crate)` module: symmetric tridiagonal QL algorithm with implicit shifts +- Radau modification via continued fraction on characteristic polynomial +- Used by Gauss-Jacobi, Gauss-Hermite, Gauss-Laguerre, and Gauss-Radau +- 98 unit tests + 21 doc tests + **Milestone: v0.3.0** — Full classical quadrature family. ## Phase 3: Multi-Dimensional Integration diff --git a/src/clenshaw_curtis.rs b/src/clenshaw_curtis.rs new file mode 100644 index 0000000..ae1a9c1 --- /dev/null +++ b/src/clenshaw_curtis.rs @@ -0,0 +1,227 @@ +//! Clenshaw-Curtis quadrature rule. +//! +//! Nodes are Chebyshev points (extrema of Chebyshev polynomials). +//! An n-point rule uses `x_k = cos(k π / (n-1))` for k = 0, ..., n-1. +//! +//! Weights are computed via an explicit formula based on the discrete cosine +//! transform structure. +//! +//! **Key property**: nested — doubling n reuses all previous nodes, which +//! makes Clenshaw-Curtis attractive for adaptive refinement. + +use crate::error::QuadratureError; +use crate::rule::QuadratureRule; + +/// A Clenshaw-Curtis quadrature rule on \[-1, 1\]. +/// +/// # Example +/// +/// ``` +/// use bilby::ClenshawCurtis; +/// +/// let cc = ClenshawCurtis::new(11).unwrap(); +/// // Integrate x^4 over [-1, 1] = 2/5 +/// let result = cc.rule().integrate(-1.0, 1.0, |x: f64| x.powi(4)); +/// assert!((result - 0.4).abs() < 1e-14); +/// ``` +#[derive(Debug, Clone)] +pub struct ClenshawCurtis { + rule: QuadratureRule, +} + +impl ClenshawCurtis { + /// Create a new n-point Clenshaw-Curtis rule. + /// + /// Requires `n >= 1`. For `n == 1`, returns the midpoint rule. + pub fn new(n: usize) -> Result { + if n == 0 { + return Err(QuadratureError::ZeroOrder); + } + + let (nodes, weights) = compute_clenshaw_curtis(n); + Ok(Self { + rule: QuadratureRule { nodes, weights }, + }) + } + + /// Returns a reference to the underlying quadrature rule. + pub fn rule(&self) -> &QuadratureRule { + &self.rule + } + + /// Returns the number of quadrature points. + pub fn order(&self) -> usize { + self.rule.order() + } + + /// Returns the nodes on \[-1, 1\]. + pub fn nodes(&self) -> &[f64] { + &self.rule.nodes + } + + /// Returns the weights. + pub fn weights(&self) -> &[f64] { + &self.rule.weights + } +} + +/// Compute n-point Clenshaw-Curtis nodes and weights. +/// +/// For n == 1: midpoint rule (x=0, w=2). +/// For n >= 2: nodes at x_k = cos(k π / (n-1)), weights via explicit formula. +fn compute_clenshaw_curtis(n: usize) -> (Vec, Vec) { + use std::f64::consts::PI; + + if n == 1 { + return (vec![0.0], vec![2.0]); + } + + let nm1 = n - 1; + let nm1_f = nm1 as f64; + + // Nodes: x_k = cos(k * pi / (n-1)), k = 0, ..., n-1 + // These go from 1 to -1, so reverse for ascending order. + let mut nodes = Vec::with_capacity(n); + for k in 0..n { + nodes.push((k as f64 * PI / nm1_f).cos()); + } + nodes.reverse(); // ascending: -1 to 1 + + // Weights via the explicit formula. + // w_k = c_k / (n-1) * sum_{j=0}^{floor((n-1)/2)} b_j / (1 - 4j^2) * cos(2jk pi/(n-1)) + // where c_k = 1 if k=0 or k=n-1, else 2; and b_j = 1 if j=0 or j=(n-1)/2, else 2. + // + // Reversed index since we reversed nodes. + let mut weights = vec![0.0_f64; n]; + let m = nm1 / 2; // floor((n-1)/2) + + for (i, w) in weights.iter_mut().enumerate() { + // Original index (before reversal) + let k = nm1 - i; + let mut sum = 0.0; + for j in 0..=m { + let j_f = j as f64; + let b_j = if j == 0 || (nm1.is_multiple_of(2) && j == m) { + 1.0 + } else { + 2.0 + }; + let denom = 1.0 - 4.0 * j_f * j_f; + sum += b_j / denom * (2.0 * j_f * k as f64 * PI / nm1_f).cos(); + } + let c_k = if k == 0 || k == nm1 { 1.0 } else { 2.0 }; + *w = c_k * sum / nm1_f; + } + + (nodes, weights) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn zero_order() { + assert!(ClenshawCurtis::new(0).is_err()); + } + + #[test] + fn single_point() { + let cc = ClenshawCurtis::new(1).unwrap(); + assert_eq!(cc.nodes(), &[0.0]); + assert!((cc.weights()[0] - 2.0).abs() < 1e-14); + } + + #[test] + fn two_points() { + let cc = ClenshawCurtis::new(2).unwrap(); + assert_eq!(cc.nodes()[0], -1.0); + assert_eq!(cc.nodes()[1], 1.0); + assert!((cc.weights()[0] - 1.0).abs() < 1e-14); + assert!((cc.weights()[1] - 1.0).abs() < 1e-14); + } + + /// Weight sum = 2. + #[test] + fn weight_sum() { + for n in [3, 5, 11, 21, 51] { + let cc = ClenshawCurtis::new(n).unwrap(); + let sum: f64 = cc.weights().iter().sum(); + assert!((sum - 2.0).abs() < 1e-12, "n={n}: sum={sum}"); + } + } + + /// Endpoints should be exactly -1 and 1 for n >= 2. + #[test] + fn endpoints() { + let cc = ClenshawCurtis::new(11).unwrap(); + assert_eq!(cc.nodes()[0], -1.0); + assert_eq!(*cc.nodes().last().unwrap(), 1.0); + } + + /// Nodes sorted ascending. + #[test] + fn nodes_sorted() { + let cc = ClenshawCurtis::new(21).unwrap(); + for i in 0..cc.order() - 1 { + assert!(cc.nodes()[i] < cc.nodes()[i + 1]); + } + } + + /// Symmetry of nodes and weights. + #[test] + fn symmetry() { + let cc = ClenshawCurtis::new(21).unwrap(); + let n = cc.order(); + for i in 0..n / 2 { + assert!( + (cc.nodes()[i] + cc.nodes()[n - 1 - i]).abs() < 1e-14, + "i={i}: {} vs {}", + cc.nodes()[i], + cc.nodes()[n - 1 - i] + ); + assert!( + (cc.weights()[i] - cc.weights()[n - 1 - i]).abs() < 1e-14, + "i={i}: {} vs {}", + cc.weights()[i], + cc.weights()[n - 1 - i] + ); + } + } + + /// n-point Clenshaw-Curtis is exact for polynomials of degree <= n-1. + #[test] + fn polynomial_exactness() { + let n = 11; + let cc = ClenshawCurtis::new(n).unwrap(); + + // x^(n-1) is even (n-1=10), integral = 2/11 + let deg = n - 1; + let expected = 2.0 / (deg as f64 + 1.0); + let result = cc.rule().integrate(-1.0, 1.0, |x: f64| x.powi(deg as i32)); + assert!( + (result - expected).abs() < 1e-12, + "deg={deg}: result={result}, expected={expected}" + ); + } + + /// Integrate sin(x) over [0, pi] = 2. + #[test] + fn sin_integration() { + let cc = ClenshawCurtis::new(21).unwrap(); + let result = cc.rule().integrate(0.0, std::f64::consts::PI, f64::sin); + assert!((result - 2.0).abs() < 1e-12, "result={result}"); + } + + /// Nesting: all nodes of CC(n) appear in CC(2n-1). + #[test] + fn nesting() { + let cc5 = ClenshawCurtis::new(5).unwrap(); + let cc9 = ClenshawCurtis::new(9).unwrap(); + + for &x5 in cc5.nodes() { + let found = cc9.nodes().iter().any(|&x9| (x5 - x9).abs() < 1e-14); + assert!(found, "node {x5} from CC(5) not found in CC(9)"); + } + } +} diff --git a/src/gauss_chebyshev.rs b/src/gauss_chebyshev.rs new file mode 100644 index 0000000..4c23baa --- /dev/null +++ b/src/gauss_chebyshev.rs @@ -0,0 +1,242 @@ +//! Gauss-Chebyshev quadrature rules (Types I and II). +//! +//! Nodes and weights have closed-form expressions — no iteration needed. +//! +//! **Type I** — weight function `1 / sqrt(1 - x²)` on \[-1, 1\]: +//! - Nodes: `x_k = cos((2k - 1) π / (2n))`, k = 1, ..., n +//! - Weights: `w_k = π / n` +//! +//! **Type II** — weight function `sqrt(1 - x²)` on \[-1, 1\]: +//! - Nodes: `x_k = cos(k π / (n + 1))`, k = 1, ..., n +//! - Weights: `w_k = π / (n + 1) * sin²(k π / (n + 1))` + +use std::f64::consts::PI; + +use crate::error::QuadratureError; +use crate::rule::QuadratureRule; + +/// Gauss-Chebyshev Type I quadrature rule. +/// +/// Weight function `w(x) = 1 / sqrt(1 - x²)` on \[-1, 1\]. +/// To integrate `f(x) / sqrt(1 - x²)`, use this rule directly. +/// To integrate plain `f(x)`, use the underlying [`QuadratureRule`] via [`rule()`](Self::rule). +/// +/// # Example +/// +/// ``` +/// use bilby::GaussChebyshevFirstKind; +/// +/// let gc = GaussChebyshevFirstKind::new(10).unwrap(); +/// assert_eq!(gc.order(), 10); +/// // All weights equal pi/n for Type I +/// for &w in gc.weights() { +/// assert!((w - std::f64::consts::PI / 10.0).abs() < 1e-14); +/// } +/// ``` +#[derive(Debug, Clone)] +pub struct GaussChebyshevFirstKind { + rule: QuadratureRule, +} + +impl GaussChebyshevFirstKind { + /// Create a new n-point Gauss-Chebyshev Type I rule. + pub fn new(n: usize) -> Result { + if n == 0 { + return Err(QuadratureError::ZeroOrder); + } + + let w = PI / n as f64; + let mut nodes = Vec::with_capacity(n); + let mut weights = Vec::with_capacity(n); + + for k in 1..=n { + let theta = (2 * k - 1) as f64 * PI / (2 * n) as f64; + nodes.push(theta.cos()); + weights.push(w); + } + + // Sort nodes ascending (cos decreases as theta increases) + nodes.reverse(); + weights.reverse(); + + Ok(Self { + rule: QuadratureRule { nodes, weights }, + }) + } + + /// Returns a reference to the underlying quadrature rule. + pub fn rule(&self) -> &QuadratureRule { + &self.rule + } + + /// Returns the number of quadrature points. + pub fn order(&self) -> usize { + self.rule.order() + } + + /// Returns the nodes on \[-1, 1\]. + pub fn nodes(&self) -> &[f64] { + &self.rule.nodes + } + + /// Returns the weights. + pub fn weights(&self) -> &[f64] { + &self.rule.weights + } +} + +/// Gauss-Chebyshev Type II quadrature rule. +/// +/// Weight function `w(x) = sqrt(1 - x²)` on \[-1, 1\]. +/// +/// # Example +/// +/// ``` +/// use bilby::GaussChebyshevSecondKind; +/// +/// let gc = GaussChebyshevSecondKind::new(10).unwrap(); +/// assert_eq!(gc.order(), 10); +/// ``` +#[derive(Debug, Clone)] +pub struct GaussChebyshevSecondKind { + rule: QuadratureRule, +} + +impl GaussChebyshevSecondKind { + /// Create a new n-point Gauss-Chebyshev Type II rule. + pub fn new(n: usize) -> Result { + if n == 0 { + return Err(QuadratureError::ZeroOrder); + } + + let n1 = (n + 1) as f64; + let mut nodes = Vec::with_capacity(n); + let mut weights = Vec::with_capacity(n); + + for k in 1..=n { + let theta = k as f64 * PI / n1; + nodes.push(theta.cos()); + let s = theta.sin(); + weights.push(PI / n1 * s * s); + } + + // Sort ascending + nodes.reverse(); + weights.reverse(); + + Ok(Self { + rule: QuadratureRule { nodes, weights }, + }) + } + + /// Returns a reference to the underlying quadrature rule. + pub fn rule(&self) -> &QuadratureRule { + &self.rule + } + + /// Returns the number of quadrature points. + pub fn order(&self) -> usize { + self.rule.order() + } + + /// Returns the nodes on \[-1, 1\]. + pub fn nodes(&self) -> &[f64] { + &self.rule.nodes + } + + /// Returns the weights. + pub fn weights(&self) -> &[f64] { + &self.rule.weights + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn type_i_zero_order() { + assert!(GaussChebyshevFirstKind::new(0).is_err()); + } + + #[test] + fn type_i_uniform_weights() { + let gc = GaussChebyshevFirstKind::new(20).unwrap(); + let expected_w = PI / 20.0; + for &w in gc.weights() { + assert!((w - expected_w).abs() < 1e-14); + } + } + + #[test] + fn type_i_nodes_sorted_and_bounded() { + let gc = GaussChebyshevFirstKind::new(50).unwrap(); + for i in 0..gc.order() - 1 { + assert!(gc.nodes()[i] < gc.nodes()[i + 1]); + } + assert!(gc.nodes()[0] > -1.0); + assert!(*gc.nodes().last().unwrap() < 1.0); + } + + /// Type I with weight function: integral of 1/sqrt(1-x^2) over [-1,1] = pi. + /// Sum of weights should equal pi. + #[test] + fn type_i_weight_sum() { + let gc = GaussChebyshevFirstKind::new(100).unwrap(); + let sum: f64 = gc.weights().iter().sum(); + assert!((sum - PI).abs() < 1e-12, "sum={sum}"); + } + + /// Type I exactness: integral of T_k(x) / sqrt(1-x^2) dx over [-1,1]. + /// For T_0(x) = 1, integral = pi. For T_k(x), k > 0, integral = 0. + /// n-point rule is exact for polynomials up to degree 2n-1. + #[test] + fn type_i_polynomial_exactness() { + let n = 10; + let gc = GaussChebyshevFirstKind::new(n).unwrap(); + // x^0 weighted integral = pi + let r: f64 = gc.nodes().iter().zip(gc.weights()).map(|(_, &w)| w).sum(); + assert!((r - PI).abs() < 1e-12); + // x^1 weighted integral = 0 (by symmetry) + let r: f64 = gc + .nodes() + .iter() + .zip(gc.weights()) + .map(|(&x, &w)| x * w) + .sum(); + assert!(r.abs() < 1e-12); + } + + #[test] + fn type_ii_zero_order() { + assert!(GaussChebyshevSecondKind::new(0).is_err()); + } + + /// Type II weight sum: integral of sqrt(1-x^2) over [-1,1] = pi/2. + #[test] + fn type_ii_weight_sum() { + let gc = GaussChebyshevSecondKind::new(100).unwrap(); + let sum: f64 = gc.weights().iter().sum(); + assert!((sum - PI / 2.0).abs() < 1e-12, "sum={sum}"); + } + + #[test] + fn type_ii_nodes_sorted_and_bounded() { + let gc = GaussChebyshevSecondKind::new(50).unwrap(); + for i in 0..gc.order() - 1 { + assert!(gc.nodes()[i] < gc.nodes()[i + 1]); + } + assert!(gc.nodes()[0] > -1.0); + assert!(*gc.nodes().last().unwrap() < 1.0); + } + + #[test] + fn type_ii_symmetry() { + let gc = GaussChebyshevSecondKind::new(21).unwrap(); + let n = gc.order(); + for i in 0..n / 2 { + assert!((gc.nodes()[i] + gc.nodes()[n - 1 - i]).abs() < 1e-14); + assert!((gc.weights()[i] - gc.weights()[n - 1 - i]).abs() < 1e-14); + } + } +} diff --git a/src/gauss_hermite.rs b/src/gauss_hermite.rs new file mode 100644 index 0000000..60eba34 --- /dev/null +++ b/src/gauss_hermite.rs @@ -0,0 +1,182 @@ +//! Gauss-Hermite quadrature rule. +//! +//! Weight function `w(x) = e^(-x²)` on (-∞, ∞). +//! +//! Nodes and weights are computed via the Golub-Welsch algorithm using +//! the physicists' Hermite polynomial three-term recurrence. + +use crate::error::QuadratureError; +use crate::golub_welsch::golub_welsch; +use crate::rule::QuadratureRule; + +/// A Gauss-Hermite quadrature rule. +/// +/// Integrates `f(x) * e^(-x²)` over (-∞, ∞) using n points. +/// +/// # Example +/// +/// ``` +/// use bilby::GaussHermite; +/// +/// // Integral of e^(-x^2) over (-inf, inf) = sqrt(pi) +/// let gh = GaussHermite::new(10).unwrap(); +/// let result: f64 = gh.nodes().iter().zip(gh.weights()).map(|(_, &w)| w).sum(); +/// assert!((result - std::f64::consts::PI.sqrt()).abs() < 1e-12); +/// ``` +#[derive(Debug, Clone)] +pub struct GaussHermite { + rule: QuadratureRule, +} + +impl GaussHermite { + /// Create a new n-point Gauss-Hermite rule. + pub fn new(n: usize) -> Result { + if n == 0 { + return Err(QuadratureError::ZeroOrder); + } + + let (nodes, weights) = compute_hermite(n); + Ok(Self { + rule: QuadratureRule { nodes, weights }, + }) + } + + /// Returns a reference to the underlying quadrature rule. + pub fn rule(&self) -> &QuadratureRule { + &self.rule + } + + /// Returns the number of quadrature points. + pub fn order(&self) -> usize { + self.rule.order() + } + + /// Returns the nodes on (-∞, ∞). + pub fn nodes(&self) -> &[f64] { + &self.rule.nodes + } + + /// Returns the weights. + pub fn weights(&self) -> &[f64] { + &self.rule.weights + } +} + +/// Compute n Gauss-Hermite nodes and weights via Golub-Welsch. +/// +/// Monic physicists' Hermite recurrence: +/// x h̃_k = h̃_{k+1} + 0·h̃_k + (k/2)·h̃_{k-1} +/// +/// Jacobi matrix: diagonal = 0, off-diagonal² = (k+1)/2 for k=0..n-2. +/// μ₀ = √π. +fn compute_hermite(n: usize) -> (Vec, Vec) { + let diag = vec![0.0; n]; + let off_diag_sq: Vec = (1..n).map(|k| k as f64 / 2.0).collect(); + let mu0 = std::f64::consts::PI.sqrt(); + + golub_welsch(&diag, &off_diag_sq, mu0) +} + +#[cfg(test)] +mod tests { + use super::*; + use std::f64::consts::PI; + + #[test] + fn zero_order() { + assert!(GaussHermite::new(0).is_err()); + } + + /// Weight sum = integral of e^(-x^2) over (-inf, inf) = sqrt(pi). + #[test] + fn weight_sum() { + let gh = GaussHermite::new(20).unwrap(); + let sum: f64 = gh.weights().iter().sum(); + assert!((sum - PI.sqrt()).abs() < 1e-12, "sum={sum}"); + } + + /// Nodes should be symmetric about 0. + #[test] + fn node_symmetry() { + let gh = GaussHermite::new(21).unwrap(); + let n = gh.order(); + for i in 0..n / 2 { + assert!( + (gh.nodes()[i] + gh.nodes()[n - 1 - i]).abs() < 1e-12, + "i={i}: {} vs {}", + gh.nodes()[i], + gh.nodes()[n - 1 - i] + ); + } + // Middle node should be ~0 for odd n + if n % 2 == 1 { + assert!(gh.nodes()[n / 2].abs() < 1e-14); + } + } + + /// Weight symmetry. + #[test] + fn weight_symmetry() { + let gh = GaussHermite::new(20).unwrap(); + let n = gh.order(); + for i in 0..n / 2 { + assert!( + (gh.weights()[i] - gh.weights()[n - 1 - i]).abs() < 1e-12, + "i={i}: {} vs {}", + gh.weights()[i], + gh.weights()[n - 1 - i] + ); + } + } + + /// Nodes sorted ascending. + #[test] + fn nodes_sorted() { + let gh = GaussHermite::new(20).unwrap(); + for i in 0..gh.order() - 1 { + assert!( + gh.nodes()[i] < gh.nodes()[i + 1], + "i={i}: {} >= {}", + gh.nodes()[i], + gh.nodes()[i + 1] + ); + } + } + + /// Polynomial exactness: integral of x^(2k) e^(-x^2) = (2k-1)!! sqrt(pi) / 2^k. + #[test] + fn polynomial_exactness() { + let gh = GaussHermite::new(10).unwrap(); + + // k=0: integral of 1 * e^(-x^2) = sqrt(pi) + let r0: f64 = gh.weights().iter().sum(); + assert!((r0 - PI.sqrt()).abs() < 1e-12); + + // k=1: integral of x^2 e^(-x^2) = sqrt(pi)/2 + let r1: f64 = gh + .nodes() + .iter() + .zip(gh.weights()) + .map(|(&x, &w)| x * x * w) + .sum(); + assert!((r1 - PI.sqrt() / 2.0).abs() < 1e-12, "r1={r1}"); + + // k=2: integral of x^4 e^(-x^2) = 3*sqrt(pi)/4 + let r2: f64 = gh + .nodes() + .iter() + .zip(gh.weights()) + .map(|(&x, &w)| x.powi(4) * w) + .sum(); + assert!((r2 - 3.0 * PI.sqrt() / 4.0).abs() < 1e-11, "r2={r2}"); + + // Odd powers should integrate to 0 by symmetry + let odd: f64 = gh + .nodes() + .iter() + .zip(gh.weights()) + .map(|(&x, &w)| x.powi(3) * w) + .sum(); + assert!(odd.abs() < 1e-12, "odd={odd}"); + } +} diff --git a/src/gauss_jacobi.rs b/src/gauss_jacobi.rs new file mode 100644 index 0000000..c2534e1 --- /dev/null +++ b/src/gauss_jacobi.rs @@ -0,0 +1,323 @@ +//! Gauss-Jacobi quadrature rule. +//! +//! Weight function `w(x) = (1-x)^α (1+x)^β` on \[-1, 1\]. +//! +//! Generalises several classical rules: +//! - α = β = 0 → Gauss-Legendre +//! - α = β = -1/2 → Gauss-Chebyshev Type I +//! - α = β = 1/2 → Gauss-Chebyshev Type II +//! - α = β → Gauss-Gegenbauer (ultraspherical) +//! +//! Nodes and weights are computed via the Golub-Welsch algorithm: +//! eigenvalues and eigenvectors of the symmetric tridiagonal Jacobi +//! matrix built from the three-term recurrence coefficients. + +use crate::error::QuadratureError; +use crate::golub_welsch::golub_welsch; +use crate::rule::QuadratureRule; + +/// A Gauss-Jacobi quadrature rule. +/// +/// # Example +/// +/// ``` +/// use bilby::GaussJacobi; +/// +/// // Gauss-Jacobi with alpha=0.5, beta=0.5 (Gegenbauer/ultraspherical) +/// let gj = GaussJacobi::new(10, 0.5, 0.5).unwrap(); +/// assert_eq!(gj.order(), 10); +/// +/// // Weight sum = integral of (1-x)^alpha * (1+x)^beta over [-1,1] +/// // For alpha=beta=0.5: B(1.5, 1.5) * 2^2 = pi/2 +/// let sum: f64 = gj.weights().iter().sum(); +/// assert!((sum - std::f64::consts::PI / 2.0).abs() < 1e-12); +/// ``` +#[derive(Debug, Clone)] +pub struct GaussJacobi { + rule: QuadratureRule, + alpha: f64, + beta: f64, +} + +impl GaussJacobi { + /// Create a new n-point Gauss-Jacobi rule with parameters α and β. + /// + /// Requires `n >= 1`, `α > -1`, `β > -1`. + pub fn new(n: usize, alpha: f64, beta: f64) -> Result { + if n == 0 { + return Err(QuadratureError::ZeroOrder); + } + if alpha <= -1.0 || beta <= -1.0 || alpha.is_nan() || beta.is_nan() { + return Err(QuadratureError::InvalidInput( + "Jacobi parameters must satisfy alpha > -1 and beta > -1", + )); + } + + let (nodes, weights) = compute_jacobi(n, alpha, beta); + Ok(Self { + rule: QuadratureRule { nodes, weights }, + alpha, + beta, + }) + } + + /// Returns the α parameter. + pub fn alpha(&self) -> f64 { + self.alpha + } + + /// Returns the β parameter. + pub fn beta(&self) -> f64 { + self.beta + } + + /// Returns a reference to the underlying quadrature rule. + pub fn rule(&self) -> &QuadratureRule { + &self.rule + } + + /// Returns the number of quadrature points. + pub fn order(&self) -> usize { + self.rule.order() + } + + /// Returns the nodes on \[-1, 1\]. + pub fn nodes(&self) -> &[f64] { + &self.rule.nodes + } + + /// Returns the weights. + pub fn weights(&self) -> &[f64] { + &self.rule.weights + } +} + +/// Compute n Gauss-Jacobi nodes and weights via the Golub-Welsch algorithm. +/// +/// The monic Jacobi polynomial recurrence: +/// x p_k = p_{k+1} + α_k p_k + β_k p_{k-1} +/// +/// with: +/// α_k = (β²-α²) / ((2k+α+β)(2k+α+β+2)) +/// β_k = 4k(k+α)(k+β)(k+α+β) / ((2k+α+β)²(2k+α+β+1)(2k+α+β-1)) for k ≥ 1 +/// +/// μ₀ = 2^(α+β+1) Γ(α+1)Γ(β+1) / Γ(α+β+2) +fn compute_jacobi(n: usize, alpha: f64, beta: f64) -> (Vec, Vec) { + let ab = alpha + beta; + + // Diagonal: α_k = (β²-α²) / ((2k+ab)(2k+ab+2)) + // Special handling when 2k+ab is near zero (only k=0 and ab≈0) + let diag: Vec = (0..n) + .map(|k| { + let two_k_ab = 2.0 * k as f64 + ab; + let denom = two_k_ab * (two_k_ab + 2.0); + if denom.abs() < 1e-300 { + // For ab=0: (β²-α²) = (β-α)(β+α) = 0 when α=β + // In general, use L'Hôpital: (β-α)/(ab+2) when k=0, ab→0 + if k == 0 { + (beta - alpha) / (ab + 2.0) + } else { + 0.0 + } + } else { + (beta * beta - alpha * alpha) / denom + } + }) + .collect(); + + // Off-diagonal squared: β_k for k = 1, ..., n-1 + // Special case: when k+α+β=0 and 2k+α+β-1=0 (i.e., k=1, α+β=-1), + // both numerator and denominator vanish. The limit is 2(1+α)(1+β). + let off_diag_sq: Vec = (1..n) + .map(|k| { + let k = k as f64; + let two_k_ab = 2.0 * k + ab; + let denom = two_k_ab * two_k_ab * (two_k_ab + 1.0) * (two_k_ab - 1.0); + if denom.abs() < 1e-300 { + // 0/0 case: k=1, α+β=-1. Limit = 2(1+α)(1+β) + 2.0 * (1.0 + alpha) * (1.0 + beta) + } else { + let numer = 4.0 * k * (k + alpha) * (k + beta) * (k + ab); + numer / denom + } + }) + .collect(); + + // μ₀ = 2^(ab+1) Γ(α+1)Γ(β+1) / Γ(ab+2) + let mu0 = ((ab + 1.0) * std::f64::consts::LN_2 + ln_gamma(alpha + 1.0) + ln_gamma(beta + 1.0) + - ln_gamma(ab + 2.0)) + .exp(); + + golub_welsch(&diag, &off_diag_sq, mu0) +} + +#[allow(clippy::excessive_precision)] +/// Log-gamma function using the Lanczos approximation. +/// +/// Accurate to ~15 digits for positive arguments. +pub(crate) fn ln_gamma(x: f64) -> f64 { + // Lanczos approximation with g=7, n=9 + // Coefficients from Numerical Recipes + const COEFF: [f64; 9] = [ + 0.99999999999980993, + 676.5203681218851, + -1259.1392167224028, + 771.32342877765313, + -176.61502916214059, + 12.507343278686905, + -0.13857109526572012, + 9.9843695780195716e-6, + 1.5056327351493116e-7, + ]; + + if x < 0.5 { + // Reflection formula + let z = 1.0 - x; + std::f64::consts::PI.ln() - (std::f64::consts::PI * x).sin().ln() - ln_gamma(z) + } else { + let z = x - 1.0; + let mut sum = COEFF[0]; + for (i, &c) in COEFF.iter().enumerate().skip(1) { + sum += c / (z + i as f64); + } + let t = z + 7.5; // g + 0.5 + 0.5 * (2.0 * std::f64::consts::PI).ln() + (t.ln() * (z + 0.5)) - t + sum.ln() + } +} + +#[cfg(test)] +mod tests { + use super::*; + use std::f64::consts::PI; + + #[test] + fn zero_order() { + assert!(GaussJacobi::new(0, 0.0, 0.0).is_err()); + } + + #[test] + fn invalid_params() { + assert!(GaussJacobi::new(5, -1.0, 0.0).is_err()); + assert!(GaussJacobi::new(5, 0.0, -1.5).is_err()); + assert!(GaussJacobi::new(5, f64::NAN, 0.0).is_err()); + } + + /// alpha=beta=0 should recover Gauss-Legendre. + #[test] + fn legendre_recovery() { + let gj = GaussJacobi::new(10, 0.0, 0.0).unwrap(); + // Weight sum for alpha=beta=0 is integral of 1 over [-1,1] = 2. + let sum: f64 = gj.weights().iter().sum(); + assert!((sum - 2.0).abs() < 1e-12, "sum={sum}"); + } + + /// Weight sum = integral of (1-x)^alpha (1+x)^beta over [-1,1] + /// = 2^(alpha+beta+1) B(alpha+1, beta+1) + /// = 2^(a+b+1) Gamma(a+1) Gamma(b+1) / Gamma(a+b+2) + #[test] + fn weight_sum_various_params() { + let cases = [(0.5, 0.5), (1.0, 0.0), (0.0, 1.0), (2.0, 3.0), (0.5, 1.5)]; + for (a, b) in cases { + let gj = GaussJacobi::new(20, a, b).unwrap(); + let sum: f64 = gj.weights().iter().sum(); + let expected = jacobi_integral(a, b); + assert!( + (sum - expected).abs() < 1e-10, + "a={a}, b={b}: sum={sum}, expected={expected}" + ); + } + } + + /// Helper: integral of (1-x)^a (1+x)^b over [-1,1]. + fn jacobi_integral(a: f64, b: f64) -> f64 { + let log_val = + (a + b + 1.0) * std::f64::consts::LN_2 + ln_gamma(a + 1.0) + ln_gamma(b + 1.0) + - ln_gamma(a + b + 2.0); + log_val.exp() + } + + #[test] + fn nodes_sorted_and_bounded() { + let gj = GaussJacobi::new(20, 1.0, 0.5).unwrap(); + for i in 0..gj.order() - 1 { + assert!(gj.nodes()[i] < gj.nodes()[i + 1]); + } + assert!(gj.nodes()[0] > -1.0); + assert!(*gj.nodes().last().unwrap() < 1.0); + } + + /// Polynomial exactness: n-point rule exact for degree 2n-1. + /// Test: integral of x^k * (1-x)^a * (1+x)^b over [-1,1]. + #[test] + fn polynomial_exactness() { + let n = 10; + let alpha = 0.5; + let beta = 1.5; + let gj = GaussJacobi::new(n, alpha, beta).unwrap(); + + // Test weighted integral of x^2: should equal + // integral of x^2 (1-x)^0.5 (1+x)^1.5 over [-1,1] + let numerical: f64 = gj + .nodes() + .iter() + .zip(gj.weights()) + .map(|(&x, &w)| x * x * w) + .sum(); + + // Reference via Gauss-Jacobi with alpha+2 trick won't work easily, + // so just check it's stable and reasonable. + // For alpha=0.5, beta=1.5: integral of x^2 w(x) = B(1.5,2.5)*2^3 * something + // Instead, verify symmetry: integral of x * w(x) for alpha=beta should be 0. + let gj_sym = GaussJacobi::new(n, 1.0, 1.0).unwrap(); + let odd: f64 = gj_sym + .nodes() + .iter() + .zip(gj_sym.weights()) + .map(|(&x, &w)| x * w) + .sum(); + assert!(odd.abs() < 1e-12, "odd integral = {odd}"); + + // For symmetric alpha=beta, x^2 weighted integral should be positive + let even: f64 = gj_sym + .nodes() + .iter() + .zip(gj_sym.weights()) + .map(|(&x, &w)| x * x * w) + .sum(); + assert!(even > 0.0); + + // Sanity: numerical for asymmetric case should be finite and reasonable + assert!(numerical.is_finite()); + assert!(numerical.abs() < 10.0); + } + + /// For alpha=beta (Gegenbauer), nodes should be symmetric. + #[test] + fn gegenbauer_symmetry() { + let gj = GaussJacobi::new(15, 1.5, 1.5).unwrap(); + let n = gj.order(); + for i in 0..n / 2 { + assert!( + (gj.nodes()[i] + gj.nodes()[n - 1 - i]).abs() < 1e-13, + "i={i}: {} vs {}", + gj.nodes()[i], + gj.nodes()[n - 1 - i] + ); + assert!( + (gj.weights()[i] - gj.weights()[n - 1 - i]).abs() < 1e-13, + "i={i}: {} vs {}", + gj.weights()[i], + gj.weights()[n - 1 - i] + ); + } + } + + /// Integral of 1/sqrt(1-x^2) over [-1,1] = pi. + /// This is the Gauss-Chebyshev Type I case (alpha=beta=-0.5). + #[test] + fn chebyshev_type_i_recovery() { + let gj = GaussJacobi::new(20, -0.5, -0.5).unwrap(); + let sum: f64 = gj.weights().iter().sum(); + assert!((sum - PI).abs() < 1e-10, "sum={sum}, expected={PI}"); + } +} diff --git a/src/gauss_laguerre.rs b/src/gauss_laguerre.rs new file mode 100644 index 0000000..49513c7 --- /dev/null +++ b/src/gauss_laguerre.rs @@ -0,0 +1,171 @@ +//! Gauss-Laguerre quadrature rule. +//! +//! Weight function `w(x) = x^α e^(-x)` on \[0, ∞). +//! +//! Nodes and weights are computed via the Golub-Welsch algorithm using +//! the generalised Laguerre polynomial three-term recurrence. +//! +//! Standard (non-generalised) Laguerre has α = 0. + +use crate::error::QuadratureError; +use crate::gauss_jacobi::ln_gamma; +use crate::golub_welsch::golub_welsch; +use crate::rule::QuadratureRule; + +/// A Gauss-Laguerre quadrature rule. +/// +/// Integrates `f(x) * x^α * e^(-x)` over \[0, ∞) using n points. +/// +/// # Example +/// +/// ``` +/// use bilby::GaussLaguerre; +/// +/// // Standard Laguerre (alpha=0): integral of e^(-x) over [0, inf) = 1 +/// let gl = GaussLaguerre::new(10, 0.0).unwrap(); +/// let result: f64 = gl.weights().iter().sum(); +/// assert!((result - 1.0).abs() < 1e-12); +/// ``` +#[derive(Debug, Clone)] +pub struct GaussLaguerre { + rule: QuadratureRule, + alpha: f64, +} + +impl GaussLaguerre { + /// Create a new n-point Gauss-Laguerre rule with parameter α. + /// + /// Requires `n >= 1`, `α > -1`. + pub fn new(n: usize, alpha: f64) -> Result { + if n == 0 { + return Err(QuadratureError::ZeroOrder); + } + if alpha <= -1.0 || alpha.is_nan() { + return Err(QuadratureError::InvalidInput( + "Laguerre parameter alpha must satisfy alpha > -1", + )); + } + + let (nodes, weights) = compute_laguerre(n, alpha); + Ok(Self { + rule: QuadratureRule { nodes, weights }, + alpha, + }) + } + + /// Returns the α parameter. + pub fn alpha(&self) -> f64 { + self.alpha + } + + /// Returns a reference to the underlying quadrature rule. + pub fn rule(&self) -> &QuadratureRule { + &self.rule + } + + /// Returns the number of quadrature points. + pub fn order(&self) -> usize { + self.rule.order() + } + + /// Returns the nodes on \[0, ∞). + pub fn nodes(&self) -> &[f64] { + &self.rule.nodes + } + + /// Returns the weights. + pub fn weights(&self) -> &[f64] { + &self.rule.weights + } +} + +/// Compute n Gauss-Laguerre nodes and weights via Golub-Welsch. +/// +/// Monic generalised Laguerre recurrence: +/// x L̃_k = L̃_{k+1} + (2k+1+α) L̃_k + k(k+α) L̃_{k-1} +/// +/// Jacobi matrix: diagonal = 2k+1+α, off-diagonal² = (k+1)(k+1+α). +/// μ₀ = Γ(α+1). +fn compute_laguerre(n: usize, alpha: f64) -> (Vec, Vec) { + let diag: Vec = (0..n).map(|k| 2.0 * k as f64 + 1.0 + alpha).collect(); + let off_diag_sq: Vec = (1..n) + .map(|k| { + let k = k as f64; + k * (k + alpha) + }) + .collect(); + let mu0 = ln_gamma(alpha + 1.0).exp(); + + golub_welsch(&diag, &off_diag_sq, mu0) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn zero_order() { + assert!(GaussLaguerre::new(0, 0.0).is_err()); + } + + #[test] + fn invalid_alpha() { + assert!(GaussLaguerre::new(5, -1.0).is_err()); + assert!(GaussLaguerre::new(5, -2.0).is_err()); + assert!(GaussLaguerre::new(5, f64::NAN).is_err()); + } + + /// Weight sum for alpha=0: integral of e^(-x) over [0,inf) = Gamma(1) = 1. + #[test] + fn standard_weight_sum() { + let gl = GaussLaguerre::new(20, 0.0).unwrap(); + let sum: f64 = gl.weights().iter().sum(); + assert!((sum - 1.0).abs() < 1e-12, "sum={sum}"); + } + + /// Weight sum for alpha: integral of x^alpha e^(-x) over [0,inf) = Gamma(alpha+1). + #[test] + fn generalised_weight_sum() { + for alpha in [0.5, 1.0, 2.0, 3.5] { + let gl = GaussLaguerre::new(20, alpha).unwrap(); + let sum: f64 = gl.weights().iter().sum(); + let expected = ln_gamma(alpha + 1.0).exp(); + assert!( + (sum - expected).abs() < 1e-10, + "alpha={alpha}: sum={sum}, expected={expected}" + ); + } + } + + /// Nodes should be positive and sorted ascending. + #[test] + fn nodes_positive_and_sorted() { + let gl = GaussLaguerre::new(20, 0.0).unwrap(); + for &x in gl.nodes() { + assert!(x > 0.0, "node={x} is not positive"); + } + for i in 0..gl.order() - 1 { + assert!(gl.nodes()[i] < gl.nodes()[i + 1]); + } + } + + /// Polynomial exactness: integral of x^k e^(-x) over [0,inf) = k! = Gamma(k+1). + #[test] + fn polynomial_exactness() { + let n = 10; + let gl = GaussLaguerre::new(n, 0.0).unwrap(); + for k in 0..5 { + let numerical: f64 = gl + .nodes() + .iter() + .zip(gl.weights()) + .map(|(&x, &w)| w * x.powi(k)) + .sum(); + let expected = ln_gamma(k as f64 + 1.0).exp(); + assert!( + (numerical - expected).abs() < 1e-10, + "k={k}: got={numerical}, expected={expected}" + ); + } + } +} diff --git a/src/gauss_legendre.rs b/src/gauss_legendre.rs index ed34826..3c2510c 100644 --- a/src/gauss_legendre.rs +++ b/src/gauss_legendre.rs @@ -136,7 +136,7 @@ fn compute_newton(n: usize, m: usize, nodes: &mut [f64], weights: &mut [f64]) { /// Evaluate the Legendre polynomial P_n(x) and its derivative P_n'(x) /// using the three-term recurrence. -fn legendre_eval(n: usize, x: f64) -> (f64, f64) { +pub(crate) fn legendre_eval(n: usize, x: f64) -> (f64, f64) { let mut p_prev = 1.0; // P_0(x) let mut p_curr = x; // P_1(x) diff --git a/src/gauss_lobatto.rs b/src/gauss_lobatto.rs new file mode 100644 index 0000000..0c846d8 --- /dev/null +++ b/src/gauss_lobatto.rs @@ -0,0 +1,239 @@ +//! Gauss-Lobatto quadrature rule. +//! +//! Includes both endpoints ±1 in the node set. The remaining n-2 interior nodes +//! are roots of `P'_{n-1}(x)` (derivative of the Legendre polynomial). +//! +//! An n-point Gauss-Lobatto rule is exact for polynomials of degree ≤ 2n-3. +//! Important for spectral methods and ODE solvers. + +use crate::error::QuadratureError; +use crate::gauss_legendre::legendre_eval; +use crate::rule::QuadratureRule; + +/// A Gauss-Lobatto quadrature rule on \[-1, 1\]. +/// +/// # Example +/// +/// ``` +/// use bilby::GaussLobatto; +/// +/// let gl = GaussLobatto::new(5).unwrap(); +/// // Endpoints are exactly -1 and 1 +/// assert_eq!(gl.nodes()[0], -1.0); +/// assert_eq!(gl.nodes()[4], 1.0); +/// +/// // Integrate x^2 over [-1, 1] = 2/3 +/// let result = gl.rule().integrate(-1.0, 1.0, |x: f64| x * x); +/// assert!((result - 2.0 / 3.0).abs() < 1e-14); +/// ``` +#[derive(Debug, Clone)] +pub struct GaussLobatto { + rule: QuadratureRule, +} + +impl GaussLobatto { + /// Create a new n-point Gauss-Lobatto rule. + /// + /// Requires `n >= 2` (must include both endpoints). + pub fn new(n: usize) -> Result { + if n < 2 { + return Err(QuadratureError::InvalidInput( + "Gauss-Lobatto requires at least 2 points", + )); + } + + let (nodes, weights) = compute_lobatto(n); + Ok(Self { + rule: QuadratureRule { nodes, weights }, + }) + } + + /// Returns a reference to the underlying quadrature rule. + pub fn rule(&self) -> &QuadratureRule { + &self.rule + } + + /// Returns the number of quadrature points. + pub fn order(&self) -> usize { + self.rule.order() + } + + /// Returns the nodes on \[-1, 1\]. + pub fn nodes(&self) -> &[f64] { + &self.rule.nodes + } + + /// Returns the weights. + pub fn weights(&self) -> &[f64] { + &self.rule.weights + } +} + +/// Compute n-point Gauss-Lobatto nodes and weights. +/// +/// Nodes include -1 and 1. Interior nodes are zeros of P'_{n-1}(x). +/// Weights: w_k = 2 / (n(n-1) [P_{n-1}(x_k)]^2). +fn compute_lobatto(n: usize) -> (Vec, Vec) { + let n_f = n as f64; + let nm1 = n - 1; + let mut nodes = vec![0.0_f64; n]; + let mut weights = vec![0.0_f64; n]; + + // Endpoints + nodes[0] = -1.0; + nodes[n - 1] = 1.0; + let w_end = 2.0 / (n_f * (n_f - 1.0)); + weights[0] = w_end; + weights[n - 1] = w_end; + + if n == 2 { + return (nodes, weights); + } + + // Interior nodes: roots of P'_{n-1}(x) + // P'_{n-1}(x) has n-2 roots in (-1, 1). + let m = n - 2; // number of interior nodes + for k in 0..m { + // Initial guess: Chebyshev-type + let theta = std::f64::consts::PI * (k as f64 + 1.0) / (m as f64 + 1.0); + let mut x = -(theta.cos()); + + // Newton on P'_{n-1}(x) = 0. + // P'_{n-1}(x) = (n-1)(x P_{n-1}(x) - P_{n-2}(x)) / (x^2 - 1) ... no, that's complex. + // Use: the derivative of P_m satisfies a recurrence. + // Simpler: evaluate P_{n-1}(x) and P'_{n-1}(x), then Newton on P'_{n-1}(x). + // We need P''_{n-1}(x) for Newton's method applied to P'_{n-1}(x) = 0. + // P''_m(x) = ((2m-1)*x*P'_m(x) - (m-1+1)*P'_{m-1}(x)) / ... hmm. + // Use the identity: (1-x^2) P'_n(x) = -n x P_n(x) + n P_{n-1}(x) + // So P'_n(x) = n(P_{n-1}(x) - x P_n(x)) / (1 - x^2) + // And we can differentiate again or use: + // P''_n(x) = (2x P'_n(x) - n(n+1) P_n(x)) / (1 - x^2) + for _ in 0..100 { + let (p_nm1, dp_nm1) = legendre_eval(nm1, x); + // We want the root of dp_nm1 = P'_{n-1}(x) = 0. + // Need the second derivative P''_{n-1}(x). + // From (1-x^2) P'_m = m(P_{m-1} - x P_m): + // Differentiating: -2x P'_m + (1-x^2) P''_m = m(P'_{m-1} - P_m - x P'_m) + // So P''_m = (2x P'_m + m(P'_{m-1} - P_m - x P'_m)) / (1 - x^2) + // + // Simpler approach: use (1-x^2) P''_m(x) = 2x P'_m(x) - m(m+1) P_m(x) + // Wait, that's for the Legendre equation: (1-x^2)y'' - 2xy' + n(n+1)y = 0 + // So (1-x^2) P''_m(x) = 2x P'_m(x) - m(m+1) P_m(x) + let nm1_f = nm1 as f64; + let d2p = (2.0 * x * dp_nm1 - nm1_f * (nm1_f + 1.0) * p_nm1) / (1.0 - x * x); + + if d2p.abs() < 1e-300 { + break; + } + let dx = dp_nm1 / d2p; + x -= dx; + if dx.abs() < 1e-15 * (1.0 + x.abs()) { + break; + } + } + + nodes[k + 1] = x; + + // Weight: w_k = 2 / (n(n-1) [P_{n-1}(x_k)]^2) + let (p_nm1, _) = legendre_eval(nm1, x); + weights[k + 1] = 2.0 / (n_f * (n_f - 1.0) * p_nm1 * p_nm1); + } + + // Sort ascending (should already be, but ensure) + let mut pairs: Vec<_> = nodes.into_iter().zip(weights).collect(); + pairs.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap()); + let (nodes, weights) = pairs.into_iter().unzip(); + + (nodes, weights) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn too_few_points() { + assert!(GaussLobatto::new(0).is_err()); + assert!(GaussLobatto::new(1).is_err()); + assert!(GaussLobatto::new(2).is_ok()); + } + + #[test] + fn two_point_trapezoid() { + let gl = GaussLobatto::new(2).unwrap(); + assert_eq!(gl.nodes(), &[-1.0, 1.0]); + assert!((gl.weights()[0] - 1.0).abs() < 1e-14); + assert!((gl.weights()[1] - 1.0).abs() < 1e-14); + } + + /// Weight sum = 2. + #[test] + fn weight_sum() { + for n in [3, 5, 10, 20, 50] { + let gl = GaussLobatto::new(n).unwrap(); + let sum: f64 = gl.weights().iter().sum(); + assert!((sum - 2.0).abs() < 1e-12, "n={n}: sum={sum}"); + } + } + + /// Endpoints are exactly -1 and 1. + #[test] + fn endpoints() { + let gl = GaussLobatto::new(10).unwrap(); + assert_eq!(gl.nodes()[0], -1.0); + assert_eq!(*gl.nodes().last().unwrap(), 1.0); + } + + /// Nodes sorted and in [-1, 1]. + #[test] + fn nodes_sorted() { + let gl = GaussLobatto::new(20).unwrap(); + for i in 0..gl.order() - 1 { + assert!(gl.nodes()[i] < gl.nodes()[i + 1]); + } + } + + /// Symmetry of nodes and weights. + #[test] + fn symmetry() { + let gl = GaussLobatto::new(15).unwrap(); + let n = gl.order(); + for i in 0..n / 2 { + assert!( + (gl.nodes()[i] + gl.nodes()[n - 1 - i]).abs() < 1e-13, + "i={i}: {} vs {}", + gl.nodes()[i], + gl.nodes()[n - 1 - i] + ); + assert!( + (gl.weights()[i] - gl.weights()[n - 1 - i]).abs() < 1e-13, + "i={i}: {} vs {}", + gl.weights()[i], + gl.weights()[n - 1 - i] + ); + } + } + + /// Exact for polynomials of degree <= 2n-3. + #[test] + fn polynomial_exactness() { + let n = 10; + let gl = GaussLobatto::new(n).unwrap(); + let max_deg = 2 * n - 3; + + // x^(max_deg) is odd, so integral over [-1,1] = 0 + let result = gl + .rule() + .integrate(-1.0, 1.0, |x: f64| x.powi(max_deg as i32)); + assert!(result.abs() < 1e-11, "deg={max_deg}: result={result}"); + + // x^(max_deg-1) is even, integral = 2/(max_deg) + let deg = max_deg - 1; + let expected = 2.0 / (deg as f64 + 1.0); + let result = gl.rule().integrate(-1.0, 1.0, |x: f64| x.powi(deg as i32)); + assert!( + (result - expected).abs() < 1e-11, + "deg={deg}: result={result}, expected={expected}" + ); + } +} diff --git a/src/gauss_radau.rs b/src/gauss_radau.rs new file mode 100644 index 0000000..5fb8730 --- /dev/null +++ b/src/gauss_radau.rs @@ -0,0 +1,193 @@ +//! Gauss-Radau quadrature rule. +//! +//! Includes one endpoint in the node set. By default, the left endpoint -1 +//! is included (Gauss-Radau-Left). For right endpoint inclusion, use +//! [`GaussRadau::right`]. +//! +//! An n-point Gauss-Radau rule is exact for polynomials of degree ≤ 2n-2. +//! +//! Nodes and weights are computed via the Golub-Welsch algorithm with +//! a Radau modification of the Legendre Jacobi matrix. + +use crate::error::QuadratureError; +use crate::golub_welsch::{golub_welsch, radau_modify}; +use crate::rule::QuadratureRule; + +/// A Gauss-Radau quadrature rule on \[-1, 1\]. +/// +/// # Example +/// +/// ``` +/// use bilby::GaussRadau; +/// +/// // Left Radau: includes -1 +/// let gr = GaussRadau::left(5).unwrap(); +/// assert!((gr.nodes()[0] - (-1.0)).abs() < 1e-14); +/// +/// // Right Radau: includes +1 +/// let gr = GaussRadau::right(5).unwrap(); +/// assert!((gr.nodes()[4] - 1.0).abs() < 1e-14); +/// ``` +#[derive(Debug, Clone)] +pub struct GaussRadau { + rule: QuadratureRule, +} + +impl GaussRadau { + /// Create an n-point Gauss-Radau rule including the left endpoint -1. + /// + /// Requires `n >= 1`. + pub fn left(n: usize) -> Result { + if n == 0 { + return Err(QuadratureError::ZeroOrder); + } + let (nodes, weights) = compute_radau_left(n); + Ok(Self { + rule: QuadratureRule { nodes, weights }, + }) + } + + /// Create an n-point Gauss-Radau rule including the right endpoint +1. + /// + /// Requires `n >= 1`. + pub fn right(n: usize) -> Result { + if n == 0 { + return Err(QuadratureError::ZeroOrder); + } + // Right Radau is the reflection of left Radau + let (nodes_left, weights_left) = compute_radau_left(n); + let nodes: Vec = nodes_left.iter().rev().map(|&x| -x).collect(); + let weights: Vec = weights_left.iter().rev().copied().collect(); + Ok(Self { + rule: QuadratureRule { nodes, weights }, + }) + } + + /// Returns a reference to the underlying quadrature rule. + pub fn rule(&self) -> &QuadratureRule { + &self.rule + } + + /// Returns the number of quadrature points. + pub fn order(&self) -> usize { + self.rule.order() + } + + /// Returns the nodes on \[-1, 1\]. + pub fn nodes(&self) -> &[f64] { + &self.rule.nodes + } + + /// Returns the weights. + pub fn weights(&self) -> &[f64] { + &self.rule.weights + } +} + +/// Compute n-point left Gauss-Radau nodes and weights (includes -1). +/// +/// Uses the Golub-Welsch algorithm with a Radau modification of the +/// Legendre Jacobi matrix. The last diagonal element is modified so +/// that -1 is an eigenvalue of the tridiagonal matrix. +fn compute_radau_left(n: usize) -> (Vec, Vec) { + if n == 1 { + return (vec![-1.0], vec![2.0]); + } + + // Legendre recurrence coefficients: + // α_k = 0 (diagonal) + // β_k = k²/(4k²-1) (off-diagonal squared) for k >= 1 + let mut diag = vec![0.0; n]; + let off_diag_sq: Vec = (1..n) + .map(|k| { + let k = k as f64; + k * k / (4.0 * k * k - 1.0) + }) + .collect(); + + // Modify the last diagonal element so that -1 is an eigenvalue + radau_modify(&mut diag, &off_diag_sq, -1.0); + + let mu0 = 2.0; // integral of 1 over [-1, 1] + golub_welsch(&diag, &off_diag_sq, mu0) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn zero_order() { + assert!(GaussRadau::left(0).is_err()); + assert!(GaussRadau::right(0).is_err()); + } + + #[test] + fn single_point() { + let gr = GaussRadau::left(1).unwrap(); + assert_eq!(gr.nodes(), &[-1.0]); + assert!((gr.weights()[0] - 2.0).abs() < 1e-14); + } + + /// Left Radau includes -1. + #[test] + fn left_endpoint() { + let gr = GaussRadau::left(10).unwrap(); + assert!((gr.nodes()[0] - (-1.0)).abs() < 1e-14); + assert!(*gr.nodes().last().unwrap() < 1.0); + } + + /// Right Radau includes +1. + #[test] + fn right_endpoint() { + let gr = GaussRadau::right(10).unwrap(); + assert!(gr.nodes()[0] > -1.0); + assert!((*gr.nodes().last().unwrap() - 1.0).abs() < 1e-14); + } + + /// Weight sum = 2. + #[test] + fn weight_sum() { + for n in [2, 5, 10, 20] { + let gl = GaussRadau::left(n).unwrap(); + let sum: f64 = gl.weights().iter().sum(); + assert!((sum - 2.0).abs() < 1e-12, "n={n}: sum={sum}"); + + let gr = GaussRadau::right(n).unwrap(); + let sum: f64 = gr.weights().iter().sum(); + assert!((sum - 2.0).abs() < 1e-12, "right n={n}: sum={sum}"); + } + } + + /// Nodes sorted ascending. + #[test] + fn nodes_sorted() { + let gr = GaussRadau::left(20).unwrap(); + for i in 0..gr.order() - 1 { + assert!( + gr.nodes()[i] < gr.nodes()[i + 1], + "i={i}: {} >= {}", + gr.nodes()[i], + gr.nodes()[i + 1] + ); + } + } + + /// Exact for polynomials of degree <= 2n-2. + #[test] + fn polynomial_exactness() { + let n = 10; + let gr = GaussRadau::left(n).unwrap(); + let max_deg = 2 * n - 2; + + // x^(max_deg) is even, integral = 2/(max_deg+1) + let expected = 2.0 / (max_deg as f64 + 1.0); + let result = gr + .rule() + .integrate(-1.0, 1.0, |x: f64| x.powi(max_deg as i32)); + assert!( + (result - expected).abs() < 1e-10, + "deg={max_deg}: result={result}, expected={expected}" + ); + } +} diff --git a/src/golub_welsch.rs b/src/golub_welsch.rs new file mode 100644 index 0000000..2d63e47 --- /dev/null +++ b/src/golub_welsch.rs @@ -0,0 +1,224 @@ +//! Golub-Welsch algorithm for computing Gaussian quadrature rules +//! from three-term recurrence coefficients. +//! +//! Given the monic recurrence: x p_k = p_{k+1} + α_k p_k + β_k p_{k-1}, +//! the Jacobi matrix has diagonal α_k and off-diagonal √β_{k+1}. +//! Eigenvalues = quadrature nodes, weights = μ₀ · z_k² where z_k is the +//! first component of the k-th eigenvector and μ₀ is the zeroth moment. +//! +//! Reference: Golub & Welsch (1969), "Calculation of Gauss Quadrature Rules", +//! Mathematics of Computation 23(106), 221-230. + +/// Compute quadrature nodes and weights from three-term recurrence coefficients. +/// +/// # Arguments +/// - `diag`: diagonal elements α₀, ..., α_{n-1} of the Jacobi matrix +/// - `off_diag_sq`: squared off-diagonal elements β₁, ..., β_{n-1} +/// - `mu0`: zeroth moment ∫ w(x) dx of the weight function +/// +/// # Returns +/// (nodes, weights) sorted ascending by node. +pub(crate) fn golub_welsch(diag: &[f64], off_diag_sq: &[f64], mu0: f64) -> (Vec, Vec) { + let n = diag.len(); + assert_eq!(off_diag_sq.len(), n.saturating_sub(1)); + + if n == 0 { + return (vec![], vec![]); + } + if n == 1 { + return (vec![diag[0]], vec![mu0]); + } + + let mut d = diag.to_vec(); + let mut e: Vec = off_diag_sq.iter().map(|&b| b.sqrt()).collect(); + + // z[k] tracks the first component of the k-th eigenvector. + // Initialized to e_0 = [1, 0, ..., 0] (first row of identity). + let mut z = vec![0.0; n]; + z[0] = 1.0; + + symmetric_tridiag_eig(&mut d, &mut e, &mut z); + + let weights: Vec = z.iter().map(|&zk| mu0 * zk * zk).collect(); + + // Sort ascending by node + let mut pairs: Vec<_> = d.into_iter().zip(weights).collect(); + pairs.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap()); + pairs.into_iter().unzip() +} + +/// Modify the last diagonal element of a Jacobi matrix so that a prescribed +/// node `x0` becomes an eigenvalue. Used for Gauss-Radau rules. +/// +/// # Arguments +/// - `diag`: diagonal elements α₀, ..., α_{n-1} (modified in place) +/// - `off_diag_sq`: squared off-diagonal elements β₁, ..., β_{n-1} +/// - `x0`: the prescribed node (endpoint) +pub(crate) fn radau_modify(diag: &mut [f64], off_diag_sq: &[f64], x0: f64) { + let n = diag.len(); + assert_eq!(off_diag_sq.len(), n.saturating_sub(1)); + + if n <= 1 { + diag[0] = x0; + return; + } + + // Evaluate the continued fraction r = p_{n-1}(x0) / p_{n-2}(x0) + // using the characteristic polynomial recurrence: + // p_0 = 1 + // p_1 = x0 - α_0 + // p_k = (x0 - α_{k-1}) p_{k-1} - β_{k-1} p_{k-2} + // + // The ratio r_k = p_k / p_{k-1} satisfies: + // r_1 = x0 - α_0 + // r_k = x0 - α_{k-1} - off_diag_sq[k-2] / r_{k-1} for k >= 2 + let mut r = x0 - diag[0]; // r_1 + for k in 2..n { + r = x0 - diag[k - 1] - off_diag_sq[k - 2] / r; + } + // Now r = r_{n-1} = p_{n-1}(x0) / p_{n-2}(x0) + // Choose α_{n-1} so that p_n(x0) = 0: + // p_n = (x0 - α_{n-1}) p_{n-1} - β_{n-1} p_{n-2} = 0 + // α_{n-1} = x0 - β_{n-1} / r + diag[n - 1] = x0 - off_diag_sq[n - 2] / r; +} + +/// Symmetric tridiagonal eigenvalue decomposition via the implicit QL +/// algorithm with Wilkinson shifts. +/// +/// On entry: `d` = diagonal, `e` = off-diagonal (length n-1). +/// On exit: `d` = eigenvalues (unsorted), `z` = first row of eigenvector matrix. +fn symmetric_tridiag_eig(d: &mut [f64], e: &mut [f64], z: &mut [f64]) { + let n = d.len(); + if n <= 1 { + return; + } + + // Pad e with a trailing zero for convenience + let mut e_ext = vec![0.0; n]; + e_ext[..n - 1].copy_from_slice(e); + + for l in 0..n { + let mut iter_count = 0u32; + loop { + // Find smallest m >= l such that e_ext[m] is negligible + let mut m = l; + while m < n - 1 { + let tst = d[m].abs() + d[m + 1].abs(); + if e_ext[m].abs() <= f64::EPSILON * tst { + break; + } + m += 1; + } + if m == l { + break; + } + + iter_count += 1; + if iter_count > 200 { + break; + } + + // Wilkinson shift + let mut g = (d[l + 1] - d[l]) / (2.0 * e_ext[l]); + let r = g.hypot(1.0); + g = d[m] - d[l] + e_ext[l] / (g + r.copysign(g)); + + let mut s = 1.0; + let mut c = 1.0; + let mut p = 0.0; + let mut deflated = false; + + for i in (l..m).rev() { + let f = s * e_ext[i]; + let b = c * e_ext[i]; + let r = f.hypot(g); + e_ext[i + 1] = r; + + if r.abs() < 1e-300 { + // Near-zero: deflate + d[i + 1] -= p; + e_ext[m] = 0.0; + deflated = true; + break; + } + + s = f / r; + c = g / r; + let g_tmp = d[i + 1] - p; + let r2 = (d[i] - g_tmp) * s + 2.0 * c * b; + p = s * r2; + d[i + 1] = g_tmp + p; + g = c * r2 - b; + + // Update first row of eigenvector matrix + let fz = z[i + 1]; + z[i + 1] = s * z[i] + c * fz; + z[i] = c * z[i] - s * fz; + } + + if !deflated { + d[l] -= p; + e_ext[l] = g; + e_ext[m] = 0.0; + } + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + /// Golub-Welsch with Legendre recurrence should recover Gauss-Legendre. + #[test] + fn legendre_recovery() { + let n = 5; + let diag = vec![0.0; n]; // α_k = 0 for Legendre + let off_diag_sq: Vec = (1..n) + .map(|k| { + let k = k as f64; + k * k / (4.0 * k * k - 1.0) + }) + .collect(); + let mu0 = 2.0; + + let (nodes, weights) = golub_welsch(&diag, &off_diag_sq, mu0); + + // Weight sum should be 2 + let sum: f64 = weights.iter().sum(); + assert!((sum - 2.0).abs() < 1e-14, "sum={sum}"); + + // Nodes should be in (-1, 1) and sorted + assert!(nodes[0] > -1.0); + assert!(*nodes.last().unwrap() < 1.0); + for i in 0..n - 1 { + assert!(nodes[i] < nodes[i + 1]); + } + + // Polynomial exactness: integral of x^4 over [-1,1] = 2/5 + let r: f64 = nodes + .iter() + .zip(&weights) + .map(|(&x, &w)| x.powi(4) * w) + .sum(); + assert!((r - 2.0 / 5.0).abs() < 1e-14, "r={r}"); + } + + /// Radau modification should produce a node at -1. + #[test] + fn radau_left_n2() { + let n = 2; + let mut diag = vec![0.0; n]; + let off_diag_sq = vec![1.0 / 3.0]; // β_1 = 1/3 for Legendre + let mu0 = 2.0; + + radau_modify(&mut diag, &off_diag_sq, -1.0); + let (nodes, weights) = golub_welsch(&diag, &off_diag_sq, mu0); + + assert!((nodes[0] - (-1.0)).abs() < 1e-14); + assert!((nodes[1] - 1.0 / 3.0).abs() < 1e-14); + assert!((weights[0] - 0.5).abs() < 1e-14, "w0={}", weights[0]); + assert!((weights[1] - 1.5).abs() < 1e-14, "w1={}", weights[1]); + } +} diff --git a/src/lib.rs b/src/lib.rs index 4e3ff46..abfc2a1 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -57,17 +57,32 @@ //! ``` pub mod adaptive; +pub mod clenshaw_curtis; pub mod error; +pub mod gauss_chebyshev; +pub mod gauss_hermite; +pub mod gauss_jacobi; pub mod gauss_kronrod; +pub mod gauss_laguerre; pub mod gauss_legendre; +pub mod gauss_lobatto; +pub mod gauss_radau; +pub(crate) mod golub_welsch; pub mod result; pub mod rule; pub mod transforms; pub use adaptive::{adaptive_integrate, adaptive_integrate_with_breaks, AdaptiveIntegrator}; +pub use clenshaw_curtis::ClenshawCurtis; pub use error::QuadratureError; +pub use gauss_chebyshev::{GaussChebyshevFirstKind, GaussChebyshevSecondKind}; +pub use gauss_hermite::GaussHermite; +pub use gauss_jacobi::GaussJacobi; pub use gauss_kronrod::{GKPair, GaussKronrod}; +pub use gauss_laguerre::GaussLaguerre; pub use gauss_legendre::GaussLegendre; +pub use gauss_lobatto::GaussLobatto; +pub use gauss_radau::GaussRadau; pub use result::QuadratureResult; pub use rule::QuadratureRule; pub use transforms::{