mathhook_core/calculus/integrals/numerical/
gaussian.rs

1//! Gaussian quadrature integration
2//!
3//! Implements Gauss-Legendre quadrature for numerical integration.
4//! Uses precomputed nodes and weights for high accuracy.
5use super::{IntegrationConfig, IntegrationResult, NumericalIntegrator};
6use crate::error::MathError;
7/// Gaussian quadrature integrator using Gauss-Legendre nodes
8pub struct GaussianQuadrature {
9    order: usize,
10}
11impl GaussianQuadrature {
12    /// Create a new Gaussian quadrature integrator
13    ///
14    /// # Arguments
15    ///
16    /// * `order` - Number of quadrature points (2, 3, 4, or 5)
17    ///
18    /// # Examples
19    ///
20    /// ```rust
21    /// use mathhook_core::calculus::integrals::numerical::GaussianQuadrature;
22    ///
23    /// let integrator = GaussianQuadrature::new(5);
24    /// ```
25    pub fn new(order: usize) -> Self {
26        assert!(
27            (2..=5).contains(&order),
28            "Gaussian quadrature order must be between 2 and 5"
29        );
30        Self { order }
31    }
32    /// Get Gauss-Legendre nodes and weights for the specified order
33    fn get_nodes_and_weights(&self) -> (&[f64], &[f64]) {
34        match self.order {
35            2 => (&GAUSS_NODES_2, &GAUSS_WEIGHTS_2),
36            3 => (&GAUSS_NODES_3, &GAUSS_WEIGHTS_3),
37            4 => (&GAUSS_NODES_4, &GAUSS_WEIGHTS_4),
38            5 => (&GAUSS_NODES_5, &GAUSS_WEIGHTS_5),
39            _ => panic!("Invalid Gaussian quadrature order"),
40        }
41    }
42    /// Transform integration from [a, b] to [-1, 1]
43    fn transform_point(&self, x: f64, a: f64, b: f64) -> f64 {
44        0.5 * ((b - a) * x + (b + a))
45    }
46    /// Compute scaled Jacobian for interval transformation
47    fn jacobian(&self, a: f64, b: f64) -> f64 {
48        0.5 * (b - a)
49    }
50}
51impl NumericalIntegrator for GaussianQuadrature {
52    fn integrate<F>(
53        &self,
54        f: F,
55        a: f64,
56        b: f64,
57        _config: &IntegrationConfig,
58    ) -> Result<IntegrationResult, MathError>
59    where
60        F: Fn(f64) -> f64,
61    {
62        if a >= b {
63            return Err(MathError::InvalidInterval { lower: a, upper: b });
64        }
65        let (nodes, weights) = self.get_nodes_and_weights();
66        let jac = self.jacobian(a, b);
67        let mut sum = 0.0;
68        for (node, weight) in nodes.iter().zip(weights.iter()) {
69            let x = self.transform_point(*node, a, b);
70            sum += weight * f(x);
71        }
72        let value = jac * sum;
73        let error_estimate = self.estimate_error(&f, a, b, value);
74        Ok(IntegrationResult {
75            value,
76            error_estimate,
77            iterations: 1,
78            subdivisions: 1,
79        })
80    }
81}
82impl GaussianQuadrature {
83    /// Estimate error using Richardson extrapolation
84    fn estimate_error<F>(&self, f: &F, a: f64, b: f64, full_value: f64) -> f64
85    where
86        F: Fn(f64) -> f64,
87    {
88        let mid = (a + b) / 2.0;
89        let (nodes, weights) = self.get_nodes_and_weights();
90        let jac_left = self.jacobian(a, mid);
91        let mut sum_left = 0.0;
92        for (node, weight) in nodes.iter().zip(weights.iter()) {
93            let x = self.transform_point(*node, a, mid);
94            sum_left += weight * f(x);
95        }
96        let value_left = jac_left * sum_left;
97        let jac_right = self.jacobian(mid, b);
98        let mut sum_right = 0.0;
99        for (node, weight) in nodes.iter().zip(weights.iter()) {
100            let x = self.transform_point(*node, mid, b);
101            sum_right += weight * f(x);
102        }
103        let value_right = jac_right * sum_right;
104        let split_value = value_left + value_right;
105        (split_value - full_value).abs() / 15.0
106    }
107}
108const GAUSS_NODES_2: [f64; 2] = [-0.5773502691896257, 0.5773502691896257];
109const GAUSS_WEIGHTS_2: [f64; 2] = [1.0, 1.0];
110const GAUSS_NODES_3: [f64; 3] = [-0.7745966692414834, 0.0, 0.7745966692414834];
111const GAUSS_WEIGHTS_3: [f64; 3] = [0.5555555555555556, 0.8888888888888888, 0.5555555555555556];
112const GAUSS_NODES_4: [f64; 4] = [
113    -0.8611363115940526,
114    -0.3399810435848563,
115    0.3399810435848563,
116    0.8611363115940526,
117];
118const GAUSS_WEIGHTS_4: [f64; 4] = [
119    0.3478548451374538,
120    0.6521451548625461,
121    0.6521451548625461,
122    0.3478548451374538,
123];
124const GAUSS_NODES_5: [f64; 5] = [
125    -0.906_179_845_938_664,
126    -0.538_469_310_105_683,
127    0.0,
128    0.538_469_310_105_683,
129    0.906_179_845_938_664,
130];
131const GAUSS_WEIGHTS_5: [f64; 5] = [
132    0.2369268850561891,
133    0.4786286704993665,
134    0.5688888888888889,
135    0.4786286704993665,
136    0.2369268850561891,
137];
138#[cfg(test)]
139mod tests {
140    use super::*;
141    #[test]
142    fn test_gaussian_quadrature_polynomial() {
143        let integrator = GaussianQuadrature::new(3);
144        let config = IntegrationConfig::default();
145        let result = integrator.integrate(|x| x * x, 0.0, 1.0, &config).unwrap();
146        assert!((result.value - 1.0 / 3.0).abs() < 1e-10);
147    }
148    #[test]
149    fn test_gaussian_quadrature_sine() {
150        let integrator = GaussianQuadrature::new(5);
151        let config = IntegrationConfig::default();
152        let result = integrator
153            .integrate(|x| x.sin(), 0.0, std::f64::consts::PI, &config)
154            .unwrap();
155        assert!((result.value - 2.0).abs() < 1e-6);
156    }
157    #[test]
158    fn test_gaussian_quadrature_exponential() {
159        let integrator = GaussianQuadrature::new(5);
160        let config = IntegrationConfig::default();
161        let result = integrator
162            .integrate(|x| x.exp(), 0.0, 1.0, &config)
163            .unwrap();
164        let expected = std::f64::consts::E - 1.0;
165        assert!((result.value - expected).abs() < 1e-10);
166    }
167    #[test]
168    fn test_gaussian_error_estimate() {
169        let integrator = GaussianQuadrature::new(3);
170        let config = IntegrationConfig::default();
171        let result = integrator
172            .integrate(|x| x * x * x, 0.0, 1.0, &config)
173            .unwrap();
174        assert!(result.error_estimate < 1e-10);
175    }
176    #[test]
177    fn test_gaussian_invalid_interval() {
178        let integrator = GaussianQuadrature::new(3);
179        let config = IntegrationConfig::default();
180        let result = integrator.integrate(|x| x, 1.0, 0.0, &config);
181        assert!(result.is_err());
182    }
183}
mathhook_core/calculus/integrals/numerical/gaussian.rs

mathhook_core/calculus/integrals/numerical/
gaussian.rs
