sequential_integration/engine/quadrature/simpson/
simpson_quadrature_double_integral.rs

1use fehler::throws;
2
3use super::{simpson_range::SimpsonRangeGenerator, utils as simpson_utils};
4use crate::{
5    engine::{
6        helper_equation_traits::EquationOfTwoVariable,
7        quadrature::{
8            FinalizeCalculation, GetQuadratureRange, GetStepSizeDoubleIntegral,
9            QuadratureDoubleIntegral,
10        },
11        range_generator::RangeGenerator,
12        Bounds, CalculationResult, CalculationStep,
13    },
14    errors::Error,
15};
16
17pub struct SimpsonQuadratureDoubleIntegral<E: Fn(f64, f64) -> f64> {
18    equation: E,
19    h: f64,
20    k: f64,
21}
22
23impl<E: Fn(f64, f64) -> f64> SimpsonQuadratureDoubleIntegral<E> {
24    #[throws]
25    pub fn new(equation: E, h: f64, k: f64) -> Self {
26        Self { equation, h, k }
27    }
28
29    #[throws]
30    fn calculate_simpson(&self, x_values: [f64; 3], y_values: [f64; 3]) -> f64 {
31        let mut f = vec![];
32
33        for x in x_values.iter() {
34            let mut f_y = vec![];
35            for y in y_values.iter() {
36                f_y.push((self.equation)(*x, *y));
37            }
38            f.push(f_y);
39        }
40
41        let result = f[0][0]
42            + f[2][0]
43            + f[0][2]
44            + f[2][2]
45            + 4. * (f[1][0] + f[0][1] + f[2][1] + f[1][2])
46            + 16. * f[1][1];
47        result
48    }
49
50    fn multiple_with_simpson_constant(value: f64, h: f64, k: f64) -> f64 {
51        h * k * value / 9.
52    }
53}
54
55impl<E: Fn(f64, f64) -> f64> EquationOfTwoVariable for SimpsonQuadratureDoubleIntegral<E> {
56    #[throws]
57    fn calculate(
58        &self,
59        x: CalculationStep,
60        bounds_x: Bounds,
61        y: CalculationStep,
62        bounds_y: Bounds,
63    ) -> CalculationResult {
64        let mut is_last_step = false;
65
66        let x = simpson_utils::SimpsonPoints::generate(x, bounds_x, self.h, &mut is_last_step);
67        let y = simpson_utils::SimpsonPoints::generate(y, bounds_y, self.k, &mut is_last_step);
68
69        let x_values = [x.v0, x.v1, x.v2];
70        let y_values = [y.v0, y.v1, y.v2];
71
72        let mut result = CalculationResult::new();
73        if is_last_step {
74            result.add_last(Self::multiple_with_simpson_constant(
75                self.calculate_simpson(x_values, y_values)?,
76                x.h,
77                y.h,
78            ));
79        } else {
80            result.add_common(self.calculate_simpson(x_values, y_values)?);
81        }
82
83        result
84    }
85}
86
87impl<E: Fn(f64, f64) -> f64> FinalizeCalculation for SimpsonQuadratureDoubleIntegral<E> {
88    #[throws]
89    fn finalize(&self, result: CalculationResult) -> f64 {
90        Self::multiple_with_simpson_constant(result.common, self.h, self.k) + result.last
91    }
92}
93
94impl<E: Fn(f64, f64) -> f64> GetStepSizeDoubleIntegral for SimpsonQuadratureDoubleIntegral<E> {
95    fn get_step_size(&self) -> (f64, f64) {
96        (self.h, self.k)
97    }
98}
99
100impl<E: Fn(f64, f64) -> f64> GetQuadratureRange for SimpsonQuadratureDoubleIntegral<E> {
101    #[throws]
102    fn get_range_generator(bounds: Bounds, h: f64) -> Option<Box<dyn RangeGenerator>> {
103        if let Some(range_generator) = SimpsonRangeGenerator::new(bounds, h)? {
104            Some(Box::new(range_generator) as Box<dyn RangeGenerator>)
105        } else {
106            None
107        }
108    }
109}
110
111impl<E: Fn(f64, f64) -> f64> QuadratureDoubleIntegral for SimpsonQuadratureDoubleIntegral<E> {}
sequential_integration/engine/quadrature/simpson/simpson_quadrature_double_integral.rs

sequential_integration/engine/quadrature/simpson/
simpson_quadrature_double_integral.rs
