gauss-quad 0.3.0

// Copyright 2019-2024 Dominique Dresen
// Copyright 2023-2026 Johanna Sörngård
// SPDX-License-Identifier: MIT OR Apache-2.0

//! Numerical integration using the Gauss-Chebyshev quadrature rule.
//!
//! This rule can integrate formulas on the form f(x) * (1 - x^2)^`a` on finite intervals, where `a` is either -1/2 or 1/2.

use crate::{
    __impl_node_weight_rule, Node, Weight,
    math::{cos, sin},
};

use alloc::boxed::Box;
use core::{f64::consts::PI, num::NonZeroUsize};

#[cfg(feature = "rayon")]
use rayon::iter::{
    IndexedParallelIterator, IntoParallelIterator, IntoParallelRefIterator, ParallelIterator,
};

/// A Gauss-Chebyshev quadrature scheme of the first kind.
///
/// Used to integrate functions of the form
/// f(x) / sqrt(1 - x^2) on finite intervals.
///
/// # Example
///
/// ```
/// # use gauss_quad::chebyshev::GaussChebyshevFirstKind;
/// # use approx::assert_relative_eq;
/// # use core::f64::consts::PI;
/// let rule = GaussChebyshevFirstKind::new(2.try_into().unwrap());
///
/// assert_relative_eq!(rule.integrate(0.0, 2.0, |x| x), PI);
/// ```
#[derive(Debug, Clone, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct GaussChebyshevFirstKind {
    node_weight_pairs: Box<[(Node, Weight)]>,
}

impl GaussChebyshevFirstKind {
    /// Create a new `GaussChebyshevFirstKind` rule that can integrate functions of the form f(x) / sqrt(1 - x^2).
    pub fn new(degree: NonZeroUsize) -> Self {
        let n = degree.get() as f64;

        Self {
            node_weight_pairs: (1..degree.get() + 1)
                .rev()
                .map(|i| (cos(PI * (2.0 * (i as f64) - 1.0) / (2.0 * n)), PI / n))
                .collect(),
        }
    }

    fn argument_transformation(x: f64, a: f64, b: f64) -> f64 {
        0.5 * ((b - a) * x + (b + a))
    }

    fn scale_factor(a: f64, b: f64) -> f64 {
        0.5 * (b - a)
    }

    #[cfg(feature = "rayon")]
    /// Same as [`new`](Self::new) but runs in parallel.
    pub fn par_new(degree: NonZeroUsize) -> Self {
        let n = degree.get() as f64;

        Self {
            node_weight_pairs: (1..degree.get() + 1)
                .into_par_iter()
                .rev()
                .map(|i| (cos(PI * (2.0 * (i as f64) - 1.0) / (2.0 * n)), PI / n))
                .collect(),
        }
    }

    /// Returns the value of the integral of the given `integrand` in the inverval \[`a`, `b`\].
    ///
    /// # Example
    ///
    /// ```
    /// # use gauss_quad::chebyshev::GaussChebyshevFirstKind;
    /// # use approx::assert_relative_eq;
    /// # use core::f64::consts::PI;
    /// let rule = GaussChebyshevFirstKind::new(2.try_into().unwrap());
    ///
    /// assert_relative_eq!(rule.integrate(-1.0, 1.0, |x| 1.5 * x * x - 0.5), PI / 4.0);
    /// ```
    pub fn integrate<F>(&self, a: f64, b: f64, mut integrand: F) -> f64
    where
        F: FnMut(f64) -> f64,
    {
        let result: f64 = self
            .node_weight_pairs
            .iter()
            .map(|(x, w)| integrand(Self::argument_transformation(*x, a, b)) * w)
            .sum();
        result * Self::scale_factor(a, b)
    }

    #[cfg(feature = "rayon")]
    /// Same as [`integrate`](Self::integrate) but runs in parallel.
    pub fn par_integrate<F>(&self, a: f64, b: f64, integrand: F) -> f64
    where
        F: Sync + Fn(f64) -> f64,
    {
        let result: f64 = self
            .node_weight_pairs
            .par_iter()
            .map(|(x, w)| integrand(Self::argument_transformation(*x, a, b)) * w)
            .sum();
        result * Self::scale_factor(a, b)
    }
}

__impl_node_weight_rule! {GaussChebyshevFirstKind, GaussChebyshevFirstKindNodes, GaussChebyshevFirstKindWeights, GaussChebyshevFirstKindIter, GaussChebyshevFirstKindIntoIter}

/// A Gauss-Chebyshev quadrature scheme of the second kind.
///
/// Used to integrate functions of the form
/// f(x) * sqrt(1 - x^2) on finite intervals.
///
/// # Example
///
/// ```
/// # use gauss_quad::chebyshev::GaussChebyshevSecondKind;
/// # use approx::assert_relative_eq;
/// # use core::f64::consts::PI;
/// let rule = GaussChebyshevSecondKind::new(2.try_into().unwrap());
///
/// assert_relative_eq!(rule.integrate(-1.0, 1.0, |x| x * x), PI / 8.0);
/// ```
#[derive(Debug, Clone, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct GaussChebyshevSecondKind {
    node_weight_pairs: Box<[(Node, Weight)]>,
}

impl GaussChebyshevSecondKind {
    /// Create a new `GaussChebyshev` rule that can integrate functions of the form f(x) * sqrt(1 - x^2).
    pub fn new(degree: NonZeroUsize) -> Self {
        let n = degree.get() as f64;

        Self {
            node_weight_pairs: (1..degree.get() + 1)
                .rev()
                .map(|i| {
                    let over_n_plus_1 = 1.0 / (n + 1.0);
                    let sin_val = sin(PI * i as f64 * over_n_plus_1);
                    (
                        cos(PI * i as f64 * over_n_plus_1),
                        PI * over_n_plus_1 * sin_val * sin_val,
                    )
                })
                .collect(),
        }
    }

    #[cfg(feature = "rayon")]
    /// Same as [`new`](Self::new) but runs in parallel.
    pub fn par_new(degree: NonZeroUsize) -> Self {
        let n = degree.get() as f64;

        Self {
            node_weight_pairs: (1..degree.get() + 1)
                .into_par_iter()
                .rev()
                .map(|i| {
                    let over_n_plus_1 = 1.0 / (n + 1.0);
                    let sin_val = sin(PI * i as f64 * over_n_plus_1);
                    (
                        cos(PI * i as f64 * over_n_plus_1),
                        PI * over_n_plus_1 * sin_val * sin_val,
                    )
                })
                .collect(),
        }
    }

    fn argument_transformation(x: f64, a: f64, b: f64) -> f64 {
        0.5 * ((b - a) * x + (b + a))
    }

    fn scale_factor(a: f64, b: f64) -> f64 {
        0.5 * (b - a)
    }

    /// Returns the value of the integral of the given `integrand` in the inverval \[`a`, `b`\].
    ///
    /// # Example
    ///
    /// ```
    /// # use gauss_quad::chebyshev::GaussChebyshevSecondKind;
    /// # use approx::assert_relative_eq;
    /// # use core::f64::consts::PI;
    /// let rule = GaussChebyshevSecondKind::new(2.try_into().unwrap());
    ///
    /// assert_relative_eq!(rule.integrate(-1.0, 1.0, |x| 1.5 * x * x - 0.5), -PI / 16.0);
    /// ```
    pub fn integrate<F>(&self, a: f64, b: f64, mut integrand: F) -> f64
    where
        F: FnMut(f64) -> f64,
    {
        let result: f64 = self
            .node_weight_pairs
            .iter()
            .map(|(x, w)| integrand(Self::argument_transformation(*x, a, b)) * w)
            .sum();
        result * Self::scale_factor(a, b)
    }

    #[cfg(feature = "rayon")]
    /// Same as [`integrate`](Self::integrate) but runs in parallel.
    pub fn par_integrate<F>(&self, a: f64, b: f64, integrand: F) -> f64
    where
        F: Fn(f64) -> f64 + Sync,
    {
        let result: f64 = self
            .node_weight_pairs
            .par_iter()
            .map(|(x, w)| integrand(Self::argument_transformation(*x, a, b)) * w)
            .sum();
        result * Self::scale_factor(a, b)
    }
}

__impl_node_weight_rule! {GaussChebyshevSecondKind, GaussChebyshevSecondKindNodes, GaussChebyshevSecondKindWeights, GaussChebyshevSecondKindIter, GaussChebyshevSecondKindIntoIter}

#[cfg(test)]
mod test {
    use crate::math::{cos, sin};
    use approx::assert_abs_diff_eq;

    use super::{GaussChebyshevFirstKind, GaussChebyshevSecondKind};

    use core::{f64::consts::PI, num::NonZeroUsize};

    #[test]
    fn check_sorted() {
        for deg in (2..100).step_by(10) {
            let rule1 = GaussChebyshevFirstKind::new(deg.try_into().unwrap());
            assert!(rule1.as_node_weight_pairs().is_sorted());
            let rule2 = GaussChebyshevSecondKind::new(deg.try_into().unwrap());
            assert!(rule2.as_node_weight_pairs().is_sorted());
        }
    }

    #[cfg(feature = "rayon")]
    #[test]
    fn check_par_sorted() {
        for deg in (2..100).step_by(10) {
            let rule1 = GaussChebyshevFirstKind::par_new(deg.try_into().unwrap());
            assert!(rule1.as_node_weight_pairs().is_sorted());
            let rule2 = GaussChebyshevSecondKind::par_new(deg.try_into().unwrap());
            assert!(rule2.as_node_weight_pairs().is_sorted());
        }
    }

    #[test]
    fn check_degree_1() {
        let rule = GaussChebyshevFirstKind::new(1.try_into().unwrap());

        assert_abs_diff_eq!(*rule.nodes().next().unwrap(), 0.0);
        assert_abs_diff_eq!(*rule.weights().next().unwrap(), PI);

        assert_abs_diff_eq!(rule.integrate(-1.0, 1.0, |x| x), 0.0);

        assert_abs_diff_eq!(rule.integrate(-1.0, 1.0, |x| x + 1.0), PI);

        // Calculated with wolfram alpha
        assert_abs_diff_eq!(rule.integrate(0.0, 1.0, |x| x), PI / 4.0);

        let rule = GaussChebyshevSecondKind::new(1.try_into().unwrap());

        assert_abs_diff_eq!(*rule.nodes().next().unwrap(), 0.0);
        assert_abs_diff_eq!(*rule.weights().next().unwrap(), PI / 2.0);

        assert_abs_diff_eq!(rule.integrate(-1.0, 1.0, |x| x), 0.0);

        // Calculated with wolfram alpha
        assert_abs_diff_eq!(rule.integrate(0.0, 1.0, |x| x), PI / 8.0);
    }

    #[test]
    fn check_chebyshev_1st_deg_5() {
        // Source: https://mathworld.wolfram.com/Chebyshev-GaussQuadrature.html
        let ans = [
            (-0.5 * (0.5 * (5.0 + f64::sqrt(5.0))).sqrt(), PI / 5.0),
            (-0.5 * (0.5 * (5.0 - f64::sqrt(5.0))).sqrt(), PI / 5.0),
            (0.0, PI / 5.0),
            (0.5 * (0.5 * (5.0 - f64::sqrt(5.0))).sqrt(), PI / 5.0),
            (0.5 * (0.5 * (5.0 + f64::sqrt(5.0))).sqrt(), PI / 5.0),
        ];

        let rule = GaussChebyshevFirstKind::new(5.try_into().unwrap());

        for ((x, w), (x_should, w_should)) in rule.into_iter().zip(ans.into_iter()) {
            assert_abs_diff_eq!(x, x_should);
            assert_abs_diff_eq!(w, w_should);
        }
    }

    #[test]
    fn check_chebyshev_2nd_deg_5() {
        // I couldn't find lists of nodes and weights to compare to. So this function computes
        // them itself with formulas from Wikipedia.

        let deg = NonZeroUsize::new(5).unwrap();
        let rule = GaussChebyshevSecondKind::new(deg);
        let deg = deg.get() as f64;

        for (i, (x, w)) in rule.into_iter().enumerate() {
            // Source: https://en.wikipedia.org/wiki/Chebyshev%E2%80%93Gauss_quadrature
            let ii = deg - i as f64;
            let x_should = cos(ii * PI / (deg + 1.0));
            let w_should = PI / (deg + 1.0) * sin(ii * PI / (deg + 1.0)).powi(2);

            assert_abs_diff_eq!(x, x_should);
            assert_abs_diff_eq!(w, w_should);
        }
    }

    #[test]
    fn check_integral_of_line() {
        let rule = GaussChebyshevFirstKind::new(2.try_into().unwrap());

        assert_abs_diff_eq!(rule.integrate(0.0, 2.0, |x| x), PI);
    }

    #[test]
    fn check_integral_of_legendre_2() {
        let deg = NonZeroUsize::new(2).unwrap();
        let rule1 = GaussChebyshevFirstKind::new(deg);
        let rule2 = GaussChebyshevSecondKind::new(deg);

        fn legendre_2(x: f64) -> f64 {
            1.5 * x * x - 0.5
        }

        assert_abs_diff_eq!(rule1.integrate(-1.0, 1.0, legendre_2), PI / 4.0);
        assert_abs_diff_eq!(rule2.integrate(-1.0, 1.0, legendre_2), -PI / 16.0);
    }

    #[test]
    fn check_integral_of_parabola() {
        let rule = GaussChebyshevSecondKind::new(2.try_into().unwrap());

        assert_abs_diff_eq!(rule.integrate(-1.0, 1.0, |x| x * x), PI / 8.0);
    }

    #[cfg(feature = "rayon")]
    #[test]
    fn test_par_integrate() {
        let rule1 = GaussChebyshevFirstKind::par_new(2.try_into().unwrap());
        let rule2 = GaussChebyshevSecondKind::par_new(2.try_into().unwrap());

        assert_abs_diff_eq!(rule1.par_integrate(0.0, 2.0, |x| x), PI);
        assert_abs_diff_eq!(rule2.par_integrate(-1.0, 1.0, |x| x * x), PI / 8.0);
    }
}