1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
//! Core quadrature rule types.
use num_traits::Float;
#[cfg(not(feature = "std"))]
use alloc::vec::Vec;
/// A precomputed quadrature rule on the reference interval \[-1, 1\].
///
/// Stores nodes and weights that can be reused across many integrations.
/// Rule construction allocates; integration with a precomputed rule does not.
#[derive(Debug, Clone)]
pub struct QuadratureRule<F> {
/// Quadrature nodes on \[-1, 1\].
pub nodes: Vec<F>,
/// Corresponding quadrature weights.
pub weights: Vec<F>,
}
impl<F: Float> QuadratureRule<F> {
/// Returns the number of points in this rule.
#[inline]
#[must_use]
pub fn order(&self) -> usize {
self.nodes.len()
}
/// Integrate `f` over \[a, b\] using this rule.
///
/// Applies an affine transform from \[-1, 1\] to \[a, b\]:
/// x = (b - a) / 2 * t + (a + b) / 2
///
/// The result is scaled by (b - a) / 2.
///
/// # Panics
///
/// Panics if the conversion of `2.0` to `F` via `F::from` fails, which
/// should not occur for standard floating-point types.
#[inline]
pub fn integrate<G>(&self, a: F, b: F, f: G) -> F
where
G: Fn(F) -> F,
{
let two = F::from(2.0).unwrap();
let half_width = (b - a) / two;
let midpoint = (a + b) / two;
let mut sum = F::zero();
for (node, weight) in self.nodes.iter().zip(self.weights.iter()) {
let x = half_width * *node + midpoint;
sum = sum + *weight * f(x);
}
sum * half_width
}
/// Integrate `f` over \[a, b\] using a composite rule with `n_panels` equal subintervals.
///
/// Each panel applies this quadrature rule independently, then sums the results.
/// Total function evaluations = `self.order() * n_panels`.
///
/// # Panics
///
/// Panics if `n_panels` or the panel indices cannot be converted to `F`
/// via `F::from`, which should not occur for standard floating-point types
/// with reasonable panel counts.
pub fn integrate_composite<G>(&self, a: F, b: F, n_panels: usize, f: G) -> F
where
G: Fn(F) -> F,
{
if n_panels == 0 {
return F::zero();
}
let n = F::from(n_panels).unwrap();
let panel_width = (b - a) / n;
let mut total = F::zero();
for i in 0..n_panels {
let panel_a = a + F::from(i).unwrap() * panel_width;
let panel_b = panel_a + panel_width;
total = total + self.integrate(panel_a, panel_b, &f);
}
total
}
/// Parallel composite integration over \[a, b\] with `n_panels` subintervals.
///
/// Equivalent to [`integrate_composite`](Self::integrate_composite) up to
/// floating-point summation order (rayon's tree reduction sums panels in a
/// different order, so results can differ in the last ULPs), but evaluates
/// panels in parallel using rayon. Requires `F: Send + Sync` and
/// `f: Fn(F) -> F + Sync`.
///
/// # Panics
///
/// Panics if `n_panels` or the panel indices cannot be converted to `F`
/// via `F::from`, which should not occur for standard floating-point types
/// with reasonable panel counts.
#[cfg(feature = "parallel")]
pub fn integrate_composite_par<G>(&self, a: F, b: F, n_panels: usize, f: G) -> F
where
F: Send + Sync,
G: Fn(F) -> F + Sync,
{
use rayon::prelude::*;
if n_panels == 0 {
return F::zero();
}
let n = F::from(n_panels).unwrap();
let panel_width = (b - a) / n;
(0..n_panels)
.into_par_iter()
.map(|i| {
let panel_a = a + F::from(i).unwrap() * panel_width;
let panel_b = panel_a + panel_width;
self.integrate(panel_a, panel_b, &f)
})
.reduce(F::zero, |a, b| a + b)
}
}
#[cfg(test)]
mod tests {
use crate::GaussLegendre;
#[test]
fn zero_panels_returns_zero() {
// n_panels == 0 must short-circuit instead of computing (b - a) / 0.
let gl = GaussLegendre::new(5).unwrap();
let r = gl.rule().integrate_composite(0.0, 1.0, 0, |x: f64| x);
assert_eq!(r, 0.0);
}
}