axiolid_numeric/
quad.rs

1//! Numerical integration of a scalar function over a finite interval.
2//!
3//! - [`integrate`]: globally adaptive Gauss-Kronrod (7-point Gauss inside a
4//!   15-point Kronrod rule), as QUADPACK `qag` (Piessens et al., 1983) and
5//!   OCCT `math_KronrodSingleIntegration`. The panel with the largest error
6//!   estimate is bisected until the total estimate meets the tolerance. No
7//!   node lies on a panel end, so integrable endpoint singularities
8//!   (`x^-1/2`, `ln x`) are handled by refinement toward them.
9//! - [`GaussLegendre`]: a fixed `n`-point rule, exact for polynomials of
10//!   degree `2n - 1`, for callers that assemble their own integrals (energy
11//!   terms over spline spans) and know the integrand's degree.
12//!
13//! The error estimate is QUADPACK's: per panel `|K15 - G7|`, rescaled by the
14//! integrand's variation and floored at the rounding level. It is an
15//! estimate, not a proof, and is conservative on smooth integrands; the
16//! tests check it against the true error on closed forms, including
17//! endpoint singularities.
18
19use std::cmp::Ordering;
20use std::collections::BinaryHeap;
21
22use crate::error::{finite, tolerance, NumericError, NumericResult, Status};
23
24/// Kronrod nodes on `[0, 1]`, descending; odd indices are the Gauss nodes.
25const XGK: [f64; 8] = [
26    0.991_455_371_120_812_6,
27    0.949_107_912_342_758_5,
28    0.864_864_423_359_769_1,
29    0.741_531_185_599_394_4,
30    0.586_087_235_467_691_1,
31    0.405_845_151_377_397_2,
32    0.207_784_955_007_898_5,
33    0.0,
34];
35/// Kronrod weights matching [`XGK`].
36const WGK: [f64; 8] = [
37    0.022_935_322_010_529_22,
38    0.063_092_092_629_978_55,
39    0.104_790_010_322_250_2,
40    0.140_653_259_715_525_9,
41    0.169_004_726_639_267_9,
42    0.190_350_578_064_785_4,
43    0.204_432_940_075_298_9,
44    0.209_482_141_084_727_8,
45];
46/// 7-point Gauss weights for `XGK[1], XGK[3], XGK[5], XGK[7]`.
47const WG: [f64; 4] = [
48    0.129_484_966_168_869_7,
49    0.279_705_391_489_276_7,
50    0.381_830_050_505_118_9,
51    0.417_959_183_673_469_4,
52];
53
54/// Accuracy request and budget for [`integrate`].
55#[derive(Debug, Clone, Copy, PartialEq)]
56pub struct IntegrationOptions {
57    /// Stop when the error estimate is at most this.
58    pub absolute_tolerance: f64,
59    /// Or at most this times `|integral|`.
60    pub relative_tolerance: f64,
61    /// Most panels the interval may be split into.
62    pub max_panels: usize,
63}
64
65impl IntegrationOptions {
66    /// Tolerances with a budget of 2000 panels (30 000 evaluations).
67    pub fn new(absolute_tolerance: f64, relative_tolerance: f64) -> Self {
68        Self {
69            absolute_tolerance,
70            relative_tolerance,
71            max_panels: 2000,
72        }
73    }
74}
75
76/// The value of an integral and how far it can be trusted.
77#[derive(Debug, Clone, Copy, PartialEq)]
78#[must_use]
79pub struct Integral {
80    /// The integral estimate.
81    pub value: f64,
82    /// Estimate of `|value - exact|`.
83    pub error_estimate: f64,
84    /// Integrand evaluations performed.
85    pub evaluations: usize,
86    /// Panels in the final partition.
87    pub panels: usize,
88    /// Whether the tolerance was met; otherwise the estimate says how
89    /// close it came.
90    pub status: Status,
91}
92
93impl Integral {
94    /// The integral if the tolerance was met.
95    ///
96    /// # Errors
97    ///
98    /// [`NumericError::NotConverged`] with the error estimate reached.
99    pub fn converged(self) -> NumericResult<Self> {
100        if self.status == Status::Converged {
101            Ok(self)
102        } else {
103            Err(NumericError::NotConverged {
104                name: "integrate",
105                error_estimate: self.error_estimate,
106            })
107        }
108    }
109}
110
111#[derive(Debug, Clone, Copy)]
112struct Panel {
113    a: f64,
114    b: f64,
115    value: f64,
116    error: f64,
117}
118
119impl PartialEq for Panel {
120    fn eq(&self, other: &Self) -> bool {
121        self.error.total_cmp(&other.error) == Ordering::Equal
122    }
123}
124impl Eq for Panel {}
125impl PartialOrd for Panel {
126    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
127        Some(self.cmp(other))
128    }
129}
130impl Ord for Panel {
131    fn cmp(&self, other: &Self) -> Ordering {
132        self.error.total_cmp(&other.error)
133    }
134}
135
136/// Integrate `f` over `[lower, upper]` by adaptive Gauss-Kronrod (G7/K15).
137///
138/// `lower > upper` integrates in reverse (negated value). The integrand is
139/// evaluated only strictly inside the interval.
140///
141/// # Errors
142///
143/// Refuses non-finite bounds or tolerances, negative tolerances, a
144/// tolerance request below the rounding level (absolute `0` with relative
145/// under `50 eps`), a zero panel budget, and any non-finite integrand
146/// value.
147pub fn integrate<F>(
148    mut f: F,
149    lower: f64,
150    upper: f64,
151    options: IntegrationOptions,
152) -> NumericResult<Integral>
153where
154    F: FnMut(f64) -> f64,
155{
156    finite(lower, "lower bound")?;
157    finite(upper, "upper bound")?;
158    let abs_tol = tolerance(options.absolute_tolerance, "absolute tolerance")?;
159    let rel_tol = tolerance(options.relative_tolerance, "relative tolerance")?;
160    if abs_tol == 0.0 && rel_tol < 50.0 * f64::EPSILON {
161        return Err(NumericError::InvalidArgument {
162            name: "tolerances",
163            reason: "must allow an error above rounding: absolute > 0 or relative >= 50 eps",
164        });
165    }
166    if options.max_panels == 0 {
167        return Err(NumericError::InvalidArgument {
168            name: "panel budget",
169            reason: "must be at least 1",
170        });
171    }
172    if lower == upper {
173        return Ok(Integral {
174            value: 0.0,
175            error_estimate: 0.0,
176            evaluations: 0,
177            panels: 0,
178            status: Status::Converged,
179        });
180    }
181    let (a, b, sign) = if lower < upper {
182        (lower, upper, 1.0)
183    } else {
184        (upper, lower, -1.0)
185    };
186    if !(b - a).is_finite() {
187        return Err(NumericError::InvalidArgument {
188            name: "interval",
189            reason: "is wider than f64 can represent",
190        });
191    }
192
193    let mut evaluations = 0usize;
194    let mut heap = BinaryHeap::new();
195    let mut settled: Vec<Panel> = Vec::new();
196    let first = kronrod(&mut f, a, b, &mut evaluations)?;
197    let mut total_value = first.value;
198    let mut total_error = first.error;
199    heap.push(first);
200    let target = |value: f64| abs_tol.max(rel_tol * value.abs());
201
202    let mut status = Status::Converged;
203    while total_error > target(total_value) {
204        if heap.len() + settled.len() >= options.max_panels {
205            status = Status::BudgetExhausted;
206            break;
207        }
208        let Some(worst) = heap.pop() else {
209            status = Status::RoundoffLimited;
210            break;
211        };
212        let mid = 0.5 * worst.a + 0.5 * worst.b;
213        // A panel too narrow to split further in f64 stays as it is.
214        let width = worst.b - worst.a;
215        if mid <= worst.a
216            || mid >= worst.b
217            || width <= 100.0 * f64::EPSILON * worst.a.abs().max(worst.b.abs())
218        {
219            settled.push(worst);
220            continue;
221        }
222        let left = kronrod(&mut f, worst.a, mid, &mut evaluations)?;
223        let right = kronrod(&mut f, mid, worst.b, &mut evaluations)?;
224        total_value += left.value + right.value - worst.value;
225        total_error += left.error + right.error - worst.error;
226        heap.push(left);
227        heap.push(right);
228        if heap.len() % 64 == 0 {
229            // Recompute the running sums to stop cancellation drift.
230            let (v, e) = sums(heap.iter().chain(settled.iter()));
231            total_value = v;
232            total_error = e;
233        }
234    }
235    let (value, error) = sums(heap.iter().chain(settled.iter()));
236    if status == Status::Converged && error > target(value) {
237        status = Status::RoundoffLimited;
238    }
239    Ok(Integral {
240        value: sign * value,
241        error_estimate: error,
242        evaluations,
243        panels: heap.len() + settled.len(),
244        status,
245    })
246}
247
248fn sums<'a>(panels: impl Iterator<Item = &'a Panel>) -> (f64, f64) {
249    let mut value = 0.0;
250    let mut compensation = 0.0;
251    let mut error = 0.0;
252    for p in panels {
253        // Kahan summation: many small panels near a singularity.
254        let y = p.value - compensation;
255        let t = value + y;
256        compensation = (t - value) - y;
257        value = t;
258        error += p.error;
259    }
260    (value, error)
261}
262
263fn kronrod<F: FnMut(f64) -> f64>(
264    f: &mut F,
265    a: f64,
266    b: f64,
267    evaluations: &mut usize,
268) -> NumericResult<Panel> {
269    let center = 0.5 * a + 0.5 * b;
270    let half = 0.5 * (b - a);
271    let mut eval = |x: f64| -> NumericResult<f64> {
272        *evaluations += 1;
273        let v = f(x);
274        if v.is_finite() {
275            Ok(v)
276        } else {
277            Err(NumericError::NonFiniteEvaluation {
278                name: "integrand",
279                at: Some(x),
280            })
281        }
282    };
283    let fc = eval(center)?;
284    let mut kronrod = WGK[7] * fc;
285    let mut gauss = WG[3] * fc;
286    let mut abs_sum = WGK[7] * fc.abs();
287    let mut values = [0.0f64; 15];
288    values[7] = fc;
289    for j in 0..7 {
290        let dx = half * XGK[j];
291        let f1 = eval(center - dx)?;
292        let f2 = eval(center + dx)?;
293        values[j] = f1;
294        values[14 - j] = f2;
295        kronrod += WGK[j] * (f1 + f2);
296        abs_sum += WGK[j] * (f1.abs() + f2.abs());
297        if j % 2 == 1 {
298            gauss += WG[j / 2] * (f1 + f2);
299        }
300    }
301    let mean = 0.5 * kronrod;
302    let mut asc = WGK[7] * (fc - mean).abs();
303    for j in 0..7 {
304        asc += WGK[j] * ((values[j] - mean).abs() + (values[14 - j] - mean).abs());
305    }
306    let result = kronrod * half;
307    let resabs = abs_sum * half.abs();
308    let resasc = asc * half.abs();
309    let mut error = ((kronrod - gauss) * half).abs();
310    if resasc != 0.0 && error != 0.0 {
311        error = resasc * (1.0f64).min((200.0 * error / resasc).powf(1.5));
312    }
313    if resabs > f64::MIN_POSITIVE / (50.0 * f64::EPSILON) {
314        error = error.max(50.0 * f64::EPSILON * resabs);
315    }
316    Ok(Panel {
317        a,
318        b,
319        value: result,
320        error,
321    })
322}
323
324/// An `n`-point Gauss-Legendre rule on `[-1, 1]`.
325///
326/// Exact for polynomials of degree up to `2n - 1`; it carries no error
327/// estimate, so use it only where the integrand's degree is known (or use
328/// [`integrate`]). Nodes come from Newton's method on `P_n` started at
329/// Tricomi's asymptotic estimates.
330#[derive(Debug, Clone, PartialEq)]
331pub struct GaussLegendre {
332    nodes: Vec<f64>,
333    weights: Vec<f64>,
334}
335
336impl GaussLegendre {
337    /// Most points supported; beyond it the three-term recurrence's
338    /// rounding dominates and an adaptive rule is the better tool.
339    pub const MAX_POINTS: usize = 512;
340
341    /// The `n`-point rule.
342    ///
343    /// # Errors
344    ///
345    /// Refuses `n == 0` and `n > MAX_POINTS`.
346    pub fn new(n: usize) -> NumericResult<Self> {
347        if n == 0 || n > Self::MAX_POINTS {
348            return Err(NumericError::InvalidArgument {
349                name: "point count",
350                reason: "must be between 1 and 512",
351            });
352        }
353        let mut nodes = vec![0.0; n];
354        let mut weights = vec![0.0; n];
355        let nf = n as f64;
356        for i in 0..n.div_ceil(2) {
357            let mut x = (std::f64::consts::PI * (i as f64 + 0.75) / (nf + 0.5)).cos();
358            let mut derivative = 0.0;
359            for _ in 0..100 {
360                let (p, dp) = legendre(n, x);
361                derivative = dp;
362                let step = p / dp;
363                x -= step;
364                if step.abs() <= 1e-16 * x.abs().max(1e-300) {
365                    break;
366                }
367            }
368            let (_, dp) = legendre(n, x);
369            if dp.is_finite() {
370                derivative = dp;
371            }
372            let w = 2.0 / ((1.0 - x * x) * derivative * derivative);
373            nodes[i] = -x;
374            nodes[n - 1 - i] = x;
375            weights[i] = w;
376            weights[n - 1 - i] = w;
377        }
378        if n % 2 == 1 {
379            nodes[n / 2] = 0.0;
380        }
381        Ok(Self { nodes, weights })
382    }
383
384    /// Nodes on `[-1, 1]`, ascending.
385    pub fn nodes(&self) -> &[f64] {
386        &self.nodes
387    }
388
389    /// Weights matching [`nodes`](Self::nodes); they sum to 2.
390    pub fn weights(&self) -> &[f64] {
391        &self.weights
392    }
393
394    /// Apply the rule to `f` on `[lower, upper]`.
395    ///
396    /// # Errors
397    ///
398    /// Refuses non-finite bounds and any non-finite value of `f`.
399    pub fn apply<F: FnMut(f64) -> f64>(
400        &self,
401        mut f: F,
402        lower: f64,
403        upper: f64,
404    ) -> NumericResult<f64> {
405        finite(lower, "lower bound")?;
406        finite(upper, "upper bound")?;
407        let center = 0.5 * lower + 0.5 * upper;
408        let half = 0.5 * upper - 0.5 * lower;
409        let mut sum = 0.0;
410        for (x, w) in self.nodes.iter().zip(&self.weights) {
411            let t = center + half * x;
412            let v = f(t);
413            if !v.is_finite() {
414                return Err(NumericError::NonFiniteEvaluation {
415                    name: "integrand",
416                    at: Some(t),
417                });
418            }
419            sum += w * v;
420        }
421        Ok(sum * half)
422    }
423}
424
425/// `P_n(x)` and `P_n'(x)` by the three-term recurrence.
426fn legendre(n: usize, x: f64) -> (f64, f64) {
427    let (mut p0, mut p1) = (1.0, x);
428    for k in 2..=n {
429        let kf = k as f64;
430        let p2 = ((2.0 * kf - 1.0) * x * p1 - (kf - 1.0) * p0) / kf;
431        p0 = p1;
432        p1 = p2;
433    }
434    if n == 0 {
435        return (1.0, 0.0);
436    }
437    let dp = n as f64 * (x * p1 - p0) / (x * x - 1.0);
438    (p1, dp)
439}
axiolid_numeric/quad.rs

axiolid_numeric/
quad.rs
