gauss_quad/simpson/
mod.rs

1//! Numerical integration using a Simpson's rule.
2//!
3//! A popular quadrature rule (also known as Kepler's barrel rule). It can be derived
4//! in the simplest case by replacing the integrand with a parabola that has the same
5//! function values at the end points a & b, as well as the Simpson m=(a+b)/2, which
6//! results in the integral formula
7//! S(f) = (b-a)/6 * [ f(a) + 4f(m) + f(b) ]
8//!
9//! Dividing the interval \[a, b\] into N neighboring intervals of length h = (b-a)/N and
10//! applying the Simpson rule to each subinterval, the integral is given by
11//!
12//! S(f) = h/6 * [ f(a) + f(b) + 2*Sum_{k=1..N-1} f(x_k) + 4*Sum_{k=1..N} f( (x_{k-1} + x_k)/2 )]
13//!
14//! with x_k = a + k*h.
15//!
16//! ```
17//! use gauss_quad::simpson::Simpson;
18//! # use gauss_quad::simpson::SimpsonError;
19//! use approx::assert_abs_diff_eq;
20//!
21//! use core::f64::consts::PI;
22//!
23//! let eps = 0.001;
24//!
25//! let n = 10;
26//! let quad = Simpson::new(n)?;
27//!
28//! // integrate some functions
29//! let integrate_euler = quad.integrate(0.0, 1.0, |x| x.exp());
30//! assert_abs_diff_eq!(integrate_euler, 1.0_f64.exp() - 1.0, epsilon = eps);
31//!
32//! let integrate_sin = quad.integrate(-PI, PI, |x| x.sin());
33//! assert_abs_diff_eq!(integrate_sin, 0.0, epsilon = eps);
34//! # Ok::<(), SimpsonError>(())
35//! ```
36
37#[cfg(feature = "rayon")]
38use rayon::prelude::{IndexedParallelIterator, IntoParallelRefIterator, ParallelIterator};
39
40use crate::{Node, __impl_node_rule};
41
42use std::backtrace::Backtrace;
43
44/// A Simpson's rule.
45///
46/// ```
47/// # use gauss_quad::simpson::{Simpson, SimpsonError};
48/// // initialize a Simpson rule with 100 subintervals
49/// let quad: Simpson = Simpson::new(100)?;
50///
51/// // numerically integrate a function from -1.0 to 1.0 using the Simpson rule
52/// let approx = quad.integrate(-1.0, 1.0, |x| x * x);
53/// # Ok::<(), SimpsonError>(())
54/// ```
55#[derive(Debug, Clone, PartialEq)]
56#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
57pub struct Simpson {
58    /// The dimensionless Simpsons nodes.
59    nodes: Vec<Node>,
60}
61
62impl Simpson {
63    /// Initialize a new Simpson rule with `degree` being the number of intervals.
64    ///
65    /// # Errors
66    ///
67    /// Returns an error if given a degree of zero.
68    pub fn new(degree: usize) -> Result<Self, SimpsonError> {
69        if degree >= 1 {
70            Ok(Self {
71                nodes: (0..degree).map(|d| d as f64).collect(),
72            })
73        } else {
74            Err(SimpsonError::new())
75        }
76    }
77
78    /// Integrate over the domain [a, b].
79    pub fn integrate<F>(&self, a: f64, b: f64, integrand: F) -> f64
80    where
81        F: Fn(f64) -> f64,
82    {
83        let n = self.nodes.len() as f64;
84
85        let h = (b - a) / n;
86
87        // first sum over the interval edges. Skips first index to sum 1..n-1
88        let sum_over_interval_edges: f64 = self
89            .nodes
90            .iter()
91            .skip(1)
92            .map(|&node| integrand(a + node * h))
93            .sum();
94
95        // sum over the midpoints f( (x_{k-1} + x_k)/2 ), as node N is not included,
96        // add it in the final result
97        let sum_over_midpoints: f64 = self
98            .nodes
99            .iter()
100            .skip(1)
101            .map(|&node| integrand(a + (2.0 * node - 1.0) * h / 2.0))
102            .sum();
103
104        h / 6.0
105            * (2.0 * sum_over_interval_edges
106                + 4.0 * sum_over_midpoints
107                + 4.0 * integrand(a + (2.0 * n - 1.0) * h / 2.0)
108                + integrand(a)
109                + integrand(b))
110    }
111
112    #[cfg(feature = "rayon")]
113    /// Same as [`integrate`](Simpson::integrate) but runs in parallel.
114    pub fn par_integrate<F>(&self, a: f64, b: f64, integrand: F) -> f64
115    where
116        F: Fn(f64) -> f64 + Sync,
117    {
118        let n = self.nodes.len() as f64;
119
120        let h = (b - a) / n;
121
122        let (sum_over_interval_edges, sum_over_midpoints): (f64, f64) = rayon::join(
123            || {
124                self.nodes
125                    .par_iter()
126                    .skip(1)
127                    .map(|&node| integrand(a + node * h))
128                    .sum::<f64>()
129            },
130            || {
131                self.nodes
132                    .par_iter()
133                    .skip(1)
134                    .map(|&node| integrand(a + (2.0 * node - 1.0) * h / 2.0))
135                    .sum::<f64>()
136            },
137        );
138
139        h / 6.0
140            * (2.0 * sum_over_interval_edges
141                + 4.0 * sum_over_midpoints
142                + 4.0 * integrand(a + (2.0 * n - 1.0) * h / 2.0)
143                + integrand(a)
144                + integrand(b))
145    }
146}
147
148__impl_node_rule! {Simpson, SimpsonIter, SimpsonIntoIter}
149
150/// The error returned by [`Simpson::new`] if given a degree of 0.
151#[derive(Debug)]
152pub struct SimpsonError(Backtrace);
153
154use core::fmt;
155impl fmt::Display for SimpsonError {
156    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
157        write!(f, "the degree of the Simpson rule must be at least 1.")
158    }
159}
160
161impl SimpsonError {
162    /// Calls [`Backtrace::capture`] and wraps the result in a `SimpsonError` struct.
163    fn new() -> Self {
164        Self(Backtrace::capture())
165    }
166
167    /// Returns a [`Backtrace`] to where the error was created.
168    ///
169    /// This backtrace is captured with [`Backtrace::capture`], see it for more information about how to make it display information when printed.
170    #[inline]
171    pub fn backtrace(&self) -> &Backtrace {
172        &self.0
173    }
174}
175
176impl std::error::Error for SimpsonError {}
177
178#[cfg(test)]
179mod tests {
180    use super::*;
181
182    #[test]
183    fn check_simpson_integration() {
184        let quad = Simpson::new(2).unwrap();
185        let integral = quad.integrate(0.0, 1.0, |x| x * x);
186        approx::assert_abs_diff_eq!(integral, 1.0 / 3.0, epsilon = 0.0001);
187    }
188
189    #[cfg(feature = "rayon")]
190    #[test]
191    fn par_check_simpson_integration() {
192        let quad = Simpson::new(2).unwrap();
193        let integral = quad.par_integrate(0.0, 1.0, |x| x * x);
194        approx::assert_abs_diff_eq!(integral, 1.0 / 3.0, epsilon = 0.0001);
195    }
196
197    #[test]
198    fn check_simpson_error() {
199        let simpson_rule = Simpson::new(0);
200        assert!(simpson_rule.is_err());
201        assert_eq!(
202            format!("{}", simpson_rule.err().unwrap()),
203            "the degree of the Simpson rule must be at least 1."
204        );
205    }
206
207    #[test]
208    fn check_derives() {
209        let quad = Simpson::new(10).unwrap();
210        let quad_clone = quad.clone();
211        assert_eq!(quad, quad_clone);
212        let other_quad = Simpson::new(3).unwrap();
213        assert_ne!(quad, other_quad);
214    }
215}
gauss_quad/simpson/mod.rs

gauss_quad/simpson/
mod.rs
