number_diff/functions/
integration.rs

1use crate::{Elementary, Error, Func, Function, Round};
2
3/// types that implement the [Integrate](crate::Integrate) trait can safely be integrated within
4/// the domain ℝ.
5pub trait Integrate {
6    /// Method will return an instance of [Integral](crate::Integral) that can be taylored to
7    /// specific usecases.
8    ///
9    /// Example:
10    /// ```rust
11    /// let function = Function::from("sin(x)");
12    ///
13    /// let mut integral = function.integrate();
14
15    /// // specify bounds and precision for the integral
16    /// integral
17    ///     .set_lower_bound(0.)
18    ///     .set_upper_bound(PI / 2.)
19    ///     .set_precision(20000);
20    ///
21    /// // evaluate the integral
22    /// let value = integral.evaluate().unwrap();
23    /// // note that the value of the evaluated integral must be unwrapped if using the `integrate()`
24    /// // method because the method cannot guarantee that bounds have been set at the point of
25    /// // evaluating. The evaluate_integral() method which is implemented for any instance with the
26    /// // Integrate trait is safer and is guaranteed to yield a valid result.
27    ///
28    /// // round the value to 10 decimal places
29    /// value.round_to(10);
30    ///
31    /// assert_eq!(value, 1.);
32    /// ```
33    /// Also note that if the precision is not specified, it will default to 1000.
34    fn integrate(&self) -> Integral;
35
36    /// Method will return the value of the definite integral of the function evaluated from the
37    /// provided lower and upper bound.
38    ///
39    /// Example:
40    /// ```
41    /// let function = Function::from("cos(x)");
42    ///
43    /// // the evaluate_integral() method will automatically round the result to five decimal points.
44    /// // This is because higher precision cannot be guaranteed with using the standard precision set
45    /// // for the method. Provided that the function is defined for all values between the lower and
46    /// // upper bounds, the method will always return a valid result.
47    /// let value = function.evaluate_integral(0., PI);
48    ///
49    /// assert_eq!(value, 0.);
50    /// ```
51    fn evaluate_integral(&self, lower_bound: f64, upper_bound: f64) -> f64;
52}
53
54const STANDARD_PRECISION: usize = 1000;
55
56/// See [Integrate documentation](crate::Integrate) for usage and examples
57pub struct Integral {
58    function: Func,
59    lower_bound: Option<f64>,
60    upper_bound: Option<f64>,
61    precision: usize,
62}
63
64impl Integral {
65    pub fn vacant(function: Func) -> Self {
66        Self {
67            function,
68            lower_bound: None,
69            upper_bound: None,
70            precision: STANDARD_PRECISION,
71        }
72    }
73
74    pub fn set_lower_bound(&mut self, lower_bound: f64) -> &mut Self {
75        self.lower_bound = Some(lower_bound);
76        self
77    }
78
79    pub fn set_upper_bound(&mut self, upper_bound: f64) -> &mut Self {
80        self.upper_bound = Some(upper_bound);
81        self
82    }
83
84    pub fn set_precision(&mut self, precision: usize) -> &mut Self {
85        self.precision = precision;
86        self
87    }
88
89    pub fn evaluate(&self) -> Result<f64, Error> {
90        if let (Some(lower_bound), Some(upper_bound)) = (self.lower_bound, self.upper_bound) {
91            Ok(simpsons_rule(
92                &self.function,
93                lower_bound,
94                upper_bound,
95                self.precision,
96            ))
97        } else {
98            Err(Error::InternalError(String::from(
99                "Bounds of integration must be set in order to evaluate the integral",
100            )))
101        }
102    }
103}
104
105/// See [Integrate](crate::Integrate) for usage and examples.
106impl Integrate for Elementary {
107    fn integrate(&self) -> Integral {
108        Integral::vacant(self.clone().call())
109    }
110    /// Evaluating the integral gives a value of the integral with eight decimal places
111    fn evaluate_integral(&self, lower_bound: f64, upper_bound: f64) -> f64 {
112        unsafe {
113            let mut value = self
114                .integrate()
115                .set_lower_bound(lower_bound)
116                .set_upper_bound(upper_bound)
117                .evaluate()
118                .unwrap_unchecked(); // this unwrap will never fail because the upper and lower bounds
119                                     // will always be set
120            value.with_significant_figures(5)
121        }
122    }
123}
124
125/// See [Integrate](crate::Integrate) for usage and examples.
126impl Integrate for Function {
127    fn integrate(&self) -> Integral {
128        self.elementary().integrate()
129    }
130    fn evaluate_integral(&self, lower_bound: f64, upper_bound: f64) -> f64 {
131        self.elementary()
132            .evaluate_integral(lower_bound, upper_bound)
133    }
134}
135
136fn simpsons_rule(funciton: &Func, lower_bound: f64, upper_bound: f64, precision: usize) -> f64 {
137    // note that n must be an even number for Simpson's rule to work
138    let n = precision * 2;
139    let dx = (upper_bound - lower_bound) / n as f64;
140    let mut sum: f64 = (1..n)
141        .map(|x| {
142            if x % 2 == 0 {
143                2. * funciton(lower_bound + x as f64 * dx)
144            } else {
145                4. * funciton(lower_bound + x as f64 * dx)
146            }
147        })
148        .sum();
149    sum += funciton(lower_bound) + funciton(upper_bound);
150
151    sum * dx / 3.
152}
number_diff/functions/integration.rs

number_diff/functions/
integration.rs
