prospicio_math/
integrate.rs

1//! Numerical integration.
2
3/// Nodes in `(0, 1)` and weights of the 8-point Gauss–Legendre rule on
4/// `[-1, 1]`; the rule is symmetric, so each node `x` stands for `±x`.
5const GL8_X: [f64; 4] = [
6    0.183_434_642_495_649_8,
7    0.525_532_409_916_329,
8    0.796_666_477_413_626_7,
9    0.960_289_856_497_536_3,
10];
11const GL8_W: [f64; 4] = [
12    0.362_683_783_378_362,
13    0.313_706_645_877_887_3,
14    0.222_381_034_453_374_5,
15    0.101_228_536_290_376_3,
16];
17
18/// `∫_lo^hi f` by the 8-point Gauss–Legendre rule, exact for polynomials
19/// of degree up to 15. `f` may fail; the first error is returned.
20///
21/// One panel, no error estimate: callers split the range into panels
22/// where `f` is smooth (between thresholds, for instance) and refine
23/// until the result stops moving.
24///
25/// ```
26/// use prospicio_math::integrate::gauss_legendre;
27///
28/// let v = gauss_legendre(|x| Ok::<_, ()>(x.powi(15)), 0.0, 1.0).unwrap();
29/// assert!((v - 1.0 / 16.0).abs() < 1e-15);
30/// ```
31pub fn gauss_legendre<E>(
32    mut f: impl FnMut(f64) -> Result<f64, E>,
33    lo: f64,
34    hi: f64,
35) -> Result<f64, E> {
36    let (mid, half) = (0.5 * (lo + hi), 0.5 * (hi - lo));
37    let mut sum = 0.0;
38    for (x, w) in GL8_X.iter().zip(GL8_W) {
39        sum += w * (f(mid - half * x)? + f(mid + half * x)?);
40    }
41    Ok(sum * half)
42}
43
44#[cfg(test)]
45mod tests {
46    use super::*;
47
48    #[test]
49    fn integrates_smooth_functions() {
50        let v = gauss_legendre(|x| Ok::<_, ()>(x.exp()), 0.0, 1.0).unwrap();
51        assert!((v - (1f64.exp() - 1.0)).abs() < 1e-15);
52    }
53
54    #[test]
55    fn passes_errors_through() {
56        let v = gauss_legendre(|x| if x > 0.5 { Err("bad") } else { Ok(x) }, 0.0, 1.0);
57        assert_eq!(v, Err("bad"));
58    }
59}
prospicio_math/integrate.rs

prospicio_math/
integrate.rs
