use_calculus/
integral.rs

1use crate::CalculusError;
2
3/// A finite interval for definite integration.
4#[derive(Debug, Clone, Copy, PartialEq)]
5pub struct IntegrationInterval {
6    start: f64,
7    end: f64,
8}
9
10impl IntegrationInterval {
11    /// Creates an interval from `start` to `end`.
12    #[must_use]
13    pub const fn new(start: f64, end: f64) -> Self {
14        Self { start, end }
15    }
16
17    /// Creates an interval from finite bounds.
18    ///
19    /// # Errors
20    ///
21    /// Returns [`CalculusError::NonFiniteBound`] when `start` or `end` is not
22    /// finite.
23    pub fn try_new(start: f64, end: f64) -> Result<Self, CalculusError> {
24        CalculusError::validate_bound("start", start)?;
25        CalculusError::validate_bound("end", end)?;
26
27        Ok(Self::new(start, end))
28    }
29
30    /// Validates that the stored bounds are finite.
31    ///
32    /// # Errors
33    ///
34    /// Returns the same error variants as [`Self::try_new`].
35    pub fn validate(self) -> Result<Self, CalculusError> {
36        Self::try_new(self.start, self.end)
37    }
38
39    /// Returns the start bound.
40    #[must_use]
41    pub const fn start(&self) -> f64 {
42        self.start
43    }
44
45    /// Returns the end bound.
46    #[must_use]
47    pub const fn end(&self) -> f64 {
48        self.end
49    }
50
51    /// Returns the signed interval width, `end - start`.
52    #[must_use]
53    pub const fn width(&self) -> f64 {
54        self.end - self.start
55    }
56}
57
58/// Numerical integration configuration.
59#[derive(Debug, Clone, Copy, PartialEq, Eq)]
60pub struct Integrator {
61    subintervals: usize,
62}
63
64impl Integrator {
65    /// Creates an integrator from a subinterval count.
66    #[must_use]
67    pub const fn new(subintervals: usize) -> Self {
68        Self { subintervals }
69    }
70
71    /// Creates an integrator from a positive subinterval count.
72    ///
73    /// # Errors
74    ///
75    /// Returns [`CalculusError::ZeroSubintervals`] when `subintervals == 0`.
76    pub fn try_new(subintervals: usize) -> Result<Self, CalculusError> {
77        CalculusError::validate_subintervals(subintervals)?;
78        Ok(Self::new(subintervals))
79    }
80
81    /// Validates that the stored subinterval count is positive.
82    ///
83    /// # Errors
84    ///
85    /// Returns the same error variants as [`Self::try_new`].
86    pub fn validate(self) -> Result<Self, CalculusError> {
87        Self::try_new(self.subintervals)
88    }
89
90    /// Returns the stored subinterval count.
91    #[must_use]
92    pub const fn subintervals(&self) -> usize {
93        self.subintervals
94    }
95
96    /// Approximates a definite integral with the trapezoidal rule.
97    ///
98    /// # Errors
99    ///
100    /// Returns [`CalculusError`] when the interval or subinterval count is
101    /// invalid, or when sampled evaluations are not finite.
102    pub fn trapezoidal<F>(
103        self,
104        function: F,
105        interval: IntegrationInterval,
106    ) -> Result<f64, CalculusError>
107    where
108        F: FnMut(f64) -> f64,
109    {
110        trapezoidal_integral(function, interval, self.subintervals)
111    }
112
113    /// Approximates a definite integral with Simpson's rule.
114    ///
115    /// # Errors
116    ///
117    /// Returns [`CalculusError`] when the interval or subinterval count is
118    /// invalid, when an odd number of subintervals is used, or when sampled
119    /// evaluations are not finite.
120    pub fn simpson<F>(
121        self,
122        function: F,
123        interval: IntegrationInterval,
124    ) -> Result<f64, CalculusError>
125    where
126        F: FnMut(f64) -> f64,
127    {
128        simpson_integral(function, interval, self.subintervals)
129    }
130}
131
132/// Approximates a definite integral with the trapezoidal rule.
133///
134/// # Errors
135///
136/// Returns [`CalculusError`] when the interval or subinterval count is
137/// invalid, or when sampled evaluations are not finite.
138#[must_use = "integral estimates should be used or handled"]
139#[allow(clippy::cast_precision_loss)]
140pub fn trapezoidal_integral<F>(
141    mut function: F,
142    interval: IntegrationInterval,
143    subintervals: usize,
144) -> Result<f64, CalculusError>
145where
146    F: FnMut(f64) -> f64,
147{
148    let interval = interval.validate()?;
149    let subintervals = CalculusError::validate_subintervals(subintervals)?;
150    let step = interval.width() / subintervals as f64;
151    let start = interval.start();
152    let end = interval.end();
153    let first = evaluate(&mut function, start)?;
154    let last = evaluate(&mut function, end)?;
155    let mut sum = 0.5 * (first + last);
156
157    for index in 1..subintervals {
158        let sample = step.mul_add(index as f64, start);
159        sum += evaluate(&mut function, sample)?;
160    }
161
162    Ok(sum * step)
163}
164
165/// Approximates a definite integral with Simpson's rule.
166///
167/// # Errors
168///
169/// Returns [`CalculusError`] when the interval or subinterval count is
170/// invalid, when an odd number of subintervals is used, or when sampled
171/// evaluations are not finite.
172///
173/// # Examples
174///
175/// ```
176/// use use_calculus::{IntegrationInterval, simpson_integral};
177///
178/// let interval = IntegrationInterval::try_new(0.0, 1.0)?;
179/// let area = simpson_integral(|x| x * x, interval, 64)?;
180///
181/// assert!((area - (1.0 / 3.0)).abs() < 1.0e-6);
182/// # Ok::<(), use_calculus::CalculusError>(())
183/// ```
184#[must_use = "integral estimates should be used or handled"]
185#[allow(clippy::cast_precision_loss)]
186pub fn simpson_integral<F>(
187    mut function: F,
188    interval: IntegrationInterval,
189    subintervals: usize,
190) -> Result<f64, CalculusError>
191where
192    F: FnMut(f64) -> f64,
193{
194    let interval = interval.validate()?;
195    let subintervals = CalculusError::validate_even_subintervals(subintervals)?;
196    let step = interval.width() / subintervals as f64;
197    let start = interval.start();
198    let end = interval.end();
199    let first = evaluate(&mut function, start)?;
200    let last = evaluate(&mut function, end)?;
201    let mut odd_sum = 0.0;
202    let mut even_sum = 0.0;
203
204    for index in 1..subintervals {
205        let sample = step.mul_add(index as f64, start);
206        let value = evaluate(&mut function, sample)?;
207
208        if index.is_multiple_of(2_usize) {
209            even_sum += value;
210        } else {
211            odd_sum += value;
212        }
213    }
214
215    let edge_sum = 4.0_f64.mul_add(odd_sum, first + last);
216    let total = 2.0_f64.mul_add(even_sum, edge_sum);
217
218    Ok((step / 3.0) * total)
219}
220
221fn evaluate<F>(function: &mut F, input: f64) -> Result<f64, CalculusError>
222where
223    F: FnMut(f64) -> f64,
224{
225    let input = CalculusError::validate_point("sample", input)?;
226    let value = function(input);
227
228    CalculusError::validate_evaluation(input, value)
229}
230
231#[cfg(test)]
232mod tests {
233    use super::{
234        CalculusError, IntegrationInterval, Integrator, simpson_integral, trapezoidal_integral,
235    };
236
237    fn assert_close(left: f64, right: f64, tolerance: f64) {
238        assert!(
239            (left - right).abs() <= tolerance,
240            "expected {left} to be within {tolerance} of {right}"
241        );
242    }
243
244    #[test]
245    fn validates_interval_bounds() {
246        assert!(matches!(
247            IntegrationInterval::try_new(f64::NAN, 1.0),
248            Err(CalculusError::NonFiniteBound { bound: "start", .. })
249        ));
250    }
251
252    #[test]
253    fn validates_subinterval_counts() {
254        assert!(matches!(
255            Integrator::try_new(0),
256            Err(CalculusError::ZeroSubintervals)
257        ));
258        assert!(matches!(
259            simpson_integral(|x| x, IntegrationInterval::new(0.0, 1.0), 3),
260            Err(CalculusError::OddSubintervalCount(3))
261        ));
262    }
263
264    #[test]
265    fn computes_trapezoidal_integrals() -> Result<(), CalculusError> {
266        let interval = IntegrationInterval::try_new(0.0, 1.0)?;
267        let area = trapezoidal_integral(|x| x * x, interval, 2_048)?;
268
269        assert_close(area, 1.0 / 3.0, 1.0e-6);
270        Ok(())
271    }
272
273    #[test]
274    fn computes_simpson_integrals() -> Result<(), CalculusError> {
275        let interval = IntegrationInterval::try_new(-1.0, 1.0)?;
276        let area = simpson_integral(|x| x * x, interval, 64)?;
277
278        assert_close(area, 2.0 / 3.0, 1.0e-10);
279        Ok(())
280    }
281
282    #[test]
283    fn rejects_non_finite_evaluations() {
284        assert!(matches!(
285            trapezoidal_integral(|_| f64::NAN, IntegrationInterval::new(0.0, 1.0), 8),
286            Err(CalculusError::NonFiniteEvaluation { .. })
287        ));
288    }
289}
use_calculus/integral.rs

use_calculus/
integral.rs
