rint/multi/
basic.rs

1use crate::multi::{Integrator, Rule};
2use crate::IntegralEstimate;
3use crate::Integrand;
4use crate::Limits;
5use crate::{InitialisationError, InitialisationErrorKind};
6
7/// A non-adaptive multi-dimensional integrator.
8///
9/// The [`Basic`] integrator applies a fully-symmetric integration [`Rule`] to approximate the
10/// integral of an $N$-dimensional function. It is non-adaptive: it runs exactly once on the input
11/// function. Thus, it is only suitable for the integration of smooth functions with no problematic
12/// regions in the integration region. If higher accuracy is required then the [`Adaptive`].
13///
14/// [`Adaptive`]: crate::multi::Adaptive
15///
16/// # Example
17///
18/// Here we present a calculation of [Catalan's constant] $G$ using the integral representation:
19/// $$
20/// G = \int_{0}^{1} \int_{0}^{1} \frac{1}{1 + x^{2} y^{2}} dy dx
21/// $$
22///
23/// [Catalan's constant]: <https://en.wikipedia.org/wiki/Catalan%27s_constant>
24///```rust
25/// use rint::{Integrand, Limits};
26/// use rint::multi::{Basic, Rule13};
27///
28/// const G: f64 = 0.915_965_594_177_219_015_054_603_514_932_384_110_774;
29/// const N: usize = 2;
30///
31/// struct Catalan;
32///
33/// impl Integrand for Catalan {
34///     type Point = [f64; N];
35///     type Scalar = f64;
36///
37///     fn evaluate(&self, coordinate: &[f64; N]) -> Self::Scalar {
38///         let x = coordinate[0];
39///         let y = coordinate[1];
40///
41///         1.0 / (1.0 + x.powi(2) * y.powi(2))
42///     }
43/// }
44///
45/// # use std::error::Error;
46/// # fn main() -> Result<(), Box<dyn Error>> {
47/// let catalan = Catalan;
48/// let limits = [Limits::new(0.0,1.0)?,Limits::new(0.0,1.0)?];
49/// let rule = Rule13::generate();
50/// let integral = Basic::new(&catalan, &rule, limits)?.integrate();
51///
52/// let result = integral.result();
53/// let error = integral.error();
54/// let abs_actual_error = (G - result).abs();
55/// let iters = integral.iterations();
56/// assert_eq!(iters, 1);
57/// assert!(abs_actual_error < error);
58/// # Ok(())
59/// # }
60///```
61#[derive(Debug, Clone, Copy, PartialEq, PartialOrd)]
62pub struct Basic<'a, I, R, const N: usize> {
63    function: &'a I,
64    rule: &'a R,
65    limits: [Limits; N],
66}
67
68impl<'a, I, const N: usize, const FINAL: usize, const TOTAL: usize>
69    Basic<'a, I, Rule<N, FINAL, TOTAL>, N>
70where
71    I: Integrand<Point = [f64; N]>,
72{
73    /// Create a new [`Basic`] multi-dimensional integrator.
74    ///
75    /// The user first defines a `function` which is something implementing the [`Integrand`] trait
76    /// and selects a fully-symmetric multi-dimensional integration [`Rule`], `rule`, to integrate
77    /// the function in the hypercube formed by the [`Limits`], `limits` in each of the $N$
78    /// integration directions.
79    ///
80    /// # Errors
81    ///
82    /// Will fail if $N < 2$ or $N > 15$. The routines probided in this module are developed
83    /// for dimensionalities between $2 \le N \le 15$.
84    pub fn new(
85        function: &'a I,
86        rule: &'a Rule<N, FINAL, TOTAL>,
87        limits: [Limits; N],
88    ) -> Result<Self, InitialisationError> {
89        if N < 2 || N > 15 {
90            return Err(InitialisationError::new(
91                InitialisationErrorKind::InvalidDimension(N),
92            ));
93        }
94        Ok(Self {
95            function,
96            rule,
97            limits,
98        })
99    }
100
101    /// Integrate the given function.
102    ///
103    /// Applies the user supplied integration rule to obtain an [`IntegralEstimate`], which is the
104    /// numerically evaluated estimate of the integral value and error, as well as the number of
105    /// function evaluations and integration routine iterations. Note: for the [`Basic`] integrator
106    /// the number of iterations is 1.
107    #[must_use]
108    pub fn integrate(&self) -> IntegralEstimate<I::Scalar> {
109        let integral = self.integrator().integrate();
110        IntegralEstimate::new()
111            .with_result(integral.result())
112            .with_error(integral.error())
113            .with_iterations(1)
114            .with_evaluations(self.rule.evaluations())
115    }
116
117    /// Create an [`Integrator`] from the integration data stored in a [`Basic`].
118    const fn integrator(&self) -> Integrator<'_, I, N, FINAL, TOTAL> {
119        Integrator::new(self.function, self.rule, self.limits)
120    }
121}
rint/multi/basic.rs

rint/multi/
basic.rs
