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 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152
use crate::{Elementary, Error, Func, Function, Round};
/// types that implement the [Integrate](crate::Integrate) trait can safely be integrated within
/// the domain ℝ.
pub trait Integrate {
/// Method will return an instance of [Integral](crate::Integral) that can be taylored to
/// specific usecases.
///
/// Example:
/// ```rust
/// let function = Function::from("sin(x)");
///
/// let mut integral = function.integrate();
/// // specify bounds and precision for the integral
/// integral
/// .set_lower_bound(0.)
/// .set_upper_bound(PI / 2.)
/// .set_precision(20000);
///
/// // evaluate the integral
/// let value = integral.evaluate().unwrap();
/// // note that the value of the evaluated integral must be unwrapped if using the `integrate()`
/// // method because the method cannot guarantee that bounds have been set at the point of
/// // evaluating. The evaluate_integral() method which is implemented for any instance with the
/// // Integrate trait is safer and is guaranteed to yield a valid result.
///
/// // round the value to 10 decimal places
/// value.round_to(10);
///
/// assert_eq!(value, 1.);
/// ```
/// Also note that if the precision is not specified, it will default to 1000.
fn integrate(&self) -> Integral;
/// Method will return the value of the definite integral of the function evaluated from the
/// provided lower and upper bound.
///
/// Example:
/// ```
/// let function = Function::from("cos(x)");
///
/// // the evaluate_integral() method will automatically round the result to five decimal points.
/// // This is because higher precision cannot be guaranteed with using the standard precision set
/// // for the method. Provided that the function is defined for all values between the lower and
/// // upper bounds, the method will always return a valid result.
/// let value = function.evaluate_integral(0., PI);
///
/// assert_eq!(value, 0.);
/// ```
fn evaluate_integral(&self, lower_bound: f64, upper_bound: f64) -> f64;
}
const STANDARD_PRECISION: usize = 1000;
/// See [Integrate documentation](crate::Integrate) for usage and examples
pub struct Integral {
function: Func,
lower_bound: Option<f64>,
upper_bound: Option<f64>,
precision: usize,
}
impl Integral {
pub fn vacant(function: Func) -> Self {
Self {
function,
lower_bound: None,
upper_bound: None,
precision: STANDARD_PRECISION,
}
}
pub fn set_lower_bound(&mut self, lower_bound: f64) -> &mut Self {
self.lower_bound = Some(lower_bound);
self
}
pub fn set_upper_bound(&mut self, upper_bound: f64) -> &mut Self {
self.upper_bound = Some(upper_bound);
self
}
pub fn set_precision(&mut self, precision: usize) -> &mut Self {
self.precision = precision;
self
}
pub fn evaluate(&self) -> Result<f64, Error> {
if let (Some(lower_bound), Some(upper_bound)) = (self.lower_bound, self.upper_bound) {
Ok(simpsons_rule(
&self.function,
lower_bound,
upper_bound,
self.precision,
))
} else {
Err(Error::InternalError(String::from(
"Bounds of integration must be set in order to evaluate the integral",
)))
}
}
}
/// See [Integrate](crate::Integrate) for usage and examples.
impl Integrate for Elementary {
fn integrate(&self) -> Integral {
Integral::vacant(self.clone().call())
}
/// Evaluating the integral gives a value of the integral with eight decimal places
fn evaluate_integral(&self, lower_bound: f64, upper_bound: f64) -> f64 {
unsafe {
let mut value = self
.integrate()
.set_lower_bound(lower_bound)
.set_upper_bound(upper_bound)
.evaluate()
.unwrap_unchecked(); // this unwrap will never fail because the upper and lower bounds
// will always be set
value.with_significant_figures(5)
}
}
}
/// See [Integrate](crate::Integrate) for usage and examples.
impl Integrate for Function {
fn integrate(&self) -> Integral {
self.elementary().integrate()
}
fn evaluate_integral(&self, lower_bound: f64, upper_bound: f64) -> f64 {
self.elementary()
.evaluate_integral(lower_bound, upper_bound)
}
}
fn simpsons_rule(funciton: &Func, lower_bound: f64, upper_bound: f64, precision: usize) -> f64 {
// note that n must be an even number for Simpson's rule to work
let n = precision * 2;
let dx = (upper_bound - lower_bound) / n as f64;
let mut sum: f64 = (1..n)
.map(|x| {
if x % 2 == 0 {
2. * funciton(lower_bound + x as f64 * dx)
} else {
4. * funciton(lower_bound + x as f64 * dx)
}
})
.sum();
sum += funciton(lower_bound) + funciton(upper_bound);
sum * dx / 3.
}