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

use crate::algebra::abstr::Real;
use crate::analysis::integral::newton_cotes::ClosedFixedIntervalIterator;

#[cfg(feature = "serde")]
use serde::{Deserialize, Serialize};

/// Closed Newton-Cotes
///
/// <https://en.wikipedia.org/wiki/Newton-Cotes_formulas>
///
#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
#[derive(Clone, Debug)]
pub struct NewtonCotes<T>
where
    T: Real,
{
    weight: Vec<T>,
}

impl<T> NewtonCotes<T>
where
    T: Real,
{
    /// # Arguments
    ///
    /// n:
    /// 1 => Trapezoidal rule
    /// 2 => Simpson's rule
    /// 3 => Simpson's 3/8 rule
    /// 4 => Boole0s rule
    /// 5 =>
    ///
    /// # Panics
    ///
    /// Panics if n < 1 || n > 5
    ///
    /// # Examples
    /// ```
    /// # #[macro_use]
    /// # extern crate mathru;
    /// # fn main()
    /// # {
    /// use mathru::analysis::integral::newton_cotes::{NewtonCotes};
    ///
    /// let nc = NewtonCotes::new(1);
    /// let f = | x | {x};
    ///
    /// let integral = nc.integrate(f, 2.0, 4.0, 4);
    ///
    /// assert_relative_eq!(integral, 6.0)
    /// # }
    /// ```
    pub fn new(n: u8) -> NewtonCotes<T> {
        if !(1..=5).contains(&n) {
            panic!("'n' is out of bounds");
        }
        let weight = match n {
            1 => vec![T::from_f32(0.5), T::from_f32(0.5)],
            2 => vec![
                T::from_f64(1.0 / 6.0),
                T::from_f64(2.0 / 3.0),
                T::from_f64(1.0 / 6.0),
            ],
            3 => vec![
                T::from_f64(1.0 / 8.0),
                T::from_f64(3.0 / 8.0),
                T::from_f64(3.0 / 8.0),
                T::from_f64(1.0 / 8.0),
            ],
            4 => vec![
                T::from_f64(7.0 / 90.0),
                T::from_f64(16.0 / 45.0),
                T::from_f64(2.0 / 15.0),
                T::from_f64(16.0 / 45.0),
                T::from_f64(7.0 / 90.0),
            ],
            5 => vec![
                T::from_f64(19.0 / 288.0),
                T::from_f64(25.0 / 96.0),
                T::from_f64(25.0 / 144.0),
                T::from_f64(25.0 / 144.0),
                T::from_f64(25.0 / 96.0),
                T::from_f64(19.0 / 288.0),
            ],
            _ => panic!(""),
        };
        NewtonCotes { weight }
    }

    ///
    /// ```math
    /// \int_{a}^{b}f(x)\,dx
    /// ```
    ///
    /// # Arguments
    /// * 'a'
    /// * 'b'
    ///  num_subintervals: \[a,b\] into smaller subintervals,
    /// #
    pub fn integrate<F>(&self, f: F, a: T, b: T, num_subintervals: u32) -> T
    where
        F: Fn(T) -> T,
    {
        let mut sub_interval: ClosedFixedIntervalIterator<T> =
            ClosedFixedIntervalIterator::new(a, b, num_subintervals);

        let mut lower_bound: T = sub_interval.next().unwrap();
        let num_intervals = self.weight.len() as u32 - 1;
        let mut sum = T::zero();
        sub_interval.for_each(|x_i: T| {
            let interval: ClosedFixedIntervalIterator<T> =
                ClosedFixedIntervalIterator::new(lower_bound, x_i, num_intervals);
            let h = x_i - lower_bound;
            let sum_i = interval
                .zip(self.weight.iter())
                .map(|(x, w)| *w * f(x))
                // .sum()
                .fold(T::zero(), |a, b| a + b)
                * h;

            sum += sum_i;

            lower_bound = x_i;
        });

        sum
    }
}