Struct AdaptiveSingularity 

pub struct AdaptiveSingularity<I> { /* private fields */ }

An adaptive Gauss-Kronrod integrator with extrapolation for integrable singularities.

Overview

The AdaptiveSingularity integrator applies a Gauss-Kronrod integration Rule to approximate the integral of a one-dimensional function. Unlike the Basic routine, the routine implemented by AdaptiveSingularity is adaptive. After the initial integration, each iteration of the routine picks the previous integration area which has the largest error estimate and bisects this region, updating the estimate to the integal and the total approximated error. Adaptive routines concentrate new subintervals around the region of highest error. If this region contains an integrable singularity, then the adaptive routine of Adaptive may fail to obtain a suitable estimate. To overcome this, AdaptiveSingularity uses an additional Wynn epsilon-algorithm extrapolation proceedure to manage these integrable singularities and provide a suitable estimate. The combination of adaptive iteration and extrapolation makes this routine suitable as a general purpose integrator, and will often out-perform the Adaptive routine (though this should be benchmarked on a case-by-case basis). It is valid for integrals defined on both finite and (semi-)infinite integration intervals. A choice of Gauss-Kronrod integration rule is not required for the AdaptiveSingularity integrator. Instead, for general functions on finite intervals the 21-point rule Rule::gk21() is used, while for (semi-)infinite intervals (see below) the lower 15-point rule Rule::gk15() is used.

The adaptive routine will return the first approximation, result, to the integral which has an absolute error smaller than the tolerance tol encoded through the Tolerance enum, where

• Tolerance::Absolute specifies absolute tolerance and returns final estimate when error <= tol,
• Tolerance::Relative specifies a relative error and returns final estimate when error <= tol * abs(result),
• Tolerance::Either to return a result as soon as either the relative or absolute error bound has been satisfied.

The routine will end when either one of the tolerance conditions have been satisfied or an error has occurred, see Error for more details.

Finite, infinite and semi-infinite integration regions

As well as being able to calculate integrals over integrable singularities within finite integration limits , AdaptiveSingularity can also be used for integrations over infinite and semi-infinite integration limits. Each of these has a dedicated constructor:

• AdaptiveSingularity::finite: finite interval $x \in (b,a)$
• AdaptiveSingularity::infinite: infinite interval $x \in (-\infty,+\infty)$
• AdaptiveSingularity::semi_infinite_lower: infinite lower limit $x \in (-\infty,a)$
• AdaptiveSingularity::semi_infinite_upper: infinite upper limit, $x \in (b,+\infty)$

The integrations over (semi-)infinite intervals are managed by transforming to a finite interval. This can introduce integrable singularities even to smooth functions, which is why the additional extrapolation provided by AdaptiveSingularity should be used for these.

Example: finite integration region

Here we present a calculation of Catalan’s constant $G$ using the integral representation: $$ G = \int_{0}^{1} \frac{\mathrm{arcsinh} x}{\sqrt{1 - x^{2}}} d x $$ which has a finite integration region $x \in (0,1)$ but difficult integral behaviour as the integration variable $x \to 1$. We use AdaptiveSingularity::finite to approximate $G$.

 use rint::{Integrand, Limits, Tolerance};
 use rint::quadrature::AdaptiveSingularity;

 const G: f64 = 0.915_965_594_177_219_015_054_603_514_932_384_110_774;
 const TOL: f64 = 1.0e-12;

 struct Catalan;

 impl Integrand for Catalan {
     type Point = f64;
     type Scalar = f64;

     fn evaluate(&self, x: &Self::Point) -> Self::Scalar {
         x.asinh() / (1.0 - x.powi(2)).sqrt()
     }
 }

 let catalan = Catalan;
 let limits = Limits::new(0.0,1.0)?;
 let tolerance = Tolerance::Relative(TOL);
 let integral = AdaptiveSingularity::finite(catalan, limits, tolerance, 1000)?
     .integrate()?;  

 let result = integral.result();
 let error = integral.error();
 let abs_actual_error = (G - result).abs();
 let tol = TOL * result.abs();
 let iters = integral.iterations();
 assert_eq!(iters, 10);
 assert!(abs_actual_error < error);
 assert!(error < tol);

Example: infinite integration region

Here we present a calculation of Euler’s constant $\gamma$ (std::f64::consts::EULER_GAMMA) using the integral representation: $$ \gamma = -\int_{-\infty}^{+\infty} x e^{x - e^{x}} dx $$ which has an infinite integration region $x\in(-\infty,+\infty)$. We use AdaptiveSingularity::infinite to approximate $\gamma$.

 use rint::{Integrand, Tolerance};
 use rint::quadrature::AdaptiveSingularity;

 const GAMMA: f64 = std::f64::consts::EULER_GAMMA;
 const TOL: f64 = 1.0e-12;

 struct Gamma;

 impl Integrand for Gamma {
     type Point = f64;
     type Scalar = f64;

     fn evaluate(&self, x: &Self::Point) -> Self::Scalar {
         let exponent = x - x.exp();
         -x * exponent.exp()
     }
 }

 let gamma = Gamma;
 let tolerance = Tolerance::Relative(TOL);
 let integral = AdaptiveSingularity::infinite(gamma, tolerance, 1000)?
     .integrate()?;

 let result = integral.result();
 let error = integral.error();
 let abs_actual_error = (GAMMA - result).abs();
 let tol = TOL * result.abs();
 let iters = integral.iterations();
 assert_eq!(iters, 11);
 assert!(abs_actual_error < error);
 assert!(error < tol);

Example: semi-infinite integration region

Here we present a calculation of Catalan’s constant $G$ using the integral representation: $$ G = \frac{\pi}{2} \int_{1}^{+\infty} \frac{(x^{4} - 6 x^{2} + 1) \ln \ln x}{(1+x^{2})^{3}} d x $$ which has a semi-infinite integration region $x \in (1,+\infty)$ We use AdaptiveSingularity::semi_infinite_upper to approximate $G$.

 use rint::{Integrand, Limits, Tolerance};
 use rint::quadrature::AdaptiveSingularity;

 const G: f64 = 0.915_965_594_177_219_015_054_603_514_932_384_110_774;
 const TOL: f64 = 1.0e-12;

 struct Catalan;

 impl Integrand for Catalan {
     type Point = f64;
     type Scalar = f64;

     fn evaluate(&self, x: &Self::Point) -> Self::Scalar {
         let FRAC_PI_2 = std::f64::consts::FRAC_PI_2;
         let polynomial = x.powi(4) - 6.0 * x.powi(2) + 1.0;
         let denominator = (1.0 + x.powi(2)).powi(3);
         let lnlnx = x.ln().ln();
         
         FRAC_PI_2 * polynomial * lnlnx / denominator
     }
 }

 let catalan = Catalan;
 let lower = 1.0;
 let tolerance = Tolerance::Relative(TOL);
 let integral = AdaptiveSingularity::semi_infinite_upper(catalan, lower, tolerance, 1000)?
     .integrate()?;

 let result = integral.result();
 let error = integral.error();
 let abs_actual_error = (G - result).abs();
 let tol = TOL * result.abs();
 let iters = integral.iterations();
 assert_eq!(iters, 32);
 assert!(abs_actual_error < error);
 assert!(error < tol);

Implementations

impl<I> AdaptiveSingularity<I>

where I: Integrand<Point = f64>,

pub fn finite( function: I, limits: Limits, tolerance: Tolerance, max_iterations: usize, ) -> Result<Self, InitialisationError>

Generate a new AdaptiveSingularity integrator for functions with possible integrable singularities on finite intervals.

Initialise an adaptive Gauss-Kronrod integrator with extrapolation for approximating integrals of the form, $$ \int_{b}^{a} f(x) dx. $$ with finite integration bounds $x \in (b,a)$.

Arguments:

• function: A user supplied function to be integrated which is something implementing the Integrand trait.
• limits: The interval over which the function should be integrated, Limits.
• tolerance: The tolerance requested by the user. Can be either an absolute tolerance or relative tolerance. Determines the exit condition of the integration routine, see Tolerance.
• max_iterations: The maximum number of iterations that the routine should use to try to satisfy the requested tolerance.

Errors

Function can return an error if it receives bad user input. This is primarily related to using invalid values for the tolerance. The returned error is an InitialisationError. See Tolerance and InitialisationError for more details.

Example

Here we present a calculation of [Catalan’s constant] $G$ using the integral representation: $$ G = \int_{0}^{1} \frac{\mathrm{arcsinh} x}{\sqrt{1 - x^{2}}} d x $$ which has a finite integration region $x \in (0,1)$ but difficult integral behaviour as the integration variable $x \to 1$. We use AdaptiveSingularity::finite to approximate $G$.

 use rint::{Integrand, Limits, Tolerance};
 use rint::quadrature::AdaptiveSingularity;

 const G: f64 = 0.915_965_594_177_219_015_054_603_514_932_384_110_774;
 const TOL: f64 = 1.0e-12;

 struct Catalan;

 impl Integrand for Catalan {
     type Point = f64;
     type Scalar = f64;

     fn evaluate(&self, x: &Self::Point) -> Self::Scalar {
         x.asinh() / (1.0 - x.powi(2)).sqrt()
     }
 }

 let catalan = Catalan;
 let limits = Limits::new(0.0,1.0)?;
 let tolerance = Tolerance::Relative(TOL);
 let integral = AdaptiveSingularity::finite(catalan, limits, tolerance, 1000)?
     .integrate()?;

 let result = integral.result();
 let error = integral.error();
 let abs_actual_error = (G - result).abs();
 let tol = TOL * result.abs();
 let iters = integral.iterations();
 assert_eq!(iters, 10);
 assert!(abs_actual_error < error);
 assert!(error < tol);

pub fn integrate( &self, ) -> Result<IntegralEstimate<I::Scalar>, IntegrationError<I::Scalar>>

Integrate the function and return a IntegralEstimate integration result.

Applies adaptive Gauss-Kronrod integration with extrapolation via the Wynn epsilon-algorithm to the user supplied function implementing the Integrand trait to generate an IntegralEstimate upon successful completion.

Errors

The error type IntegrationError will return both the error IntegrationErrorKind and the IntegralEstimate obtained before an error was encountered. The integration routine has several ways of failing:

• The user supplied Tolerance could not be satisfied within the maximum number of iterations, see IntegrationErrorKind::MaximumIterationsReached.

• A roundoff error was detected. Can occur when the calculated numerical error from an internal integration is smaller than the estimated roundoff, but larger than the tolerance requested by the user, or when too many successive iterations do not reasonably improve the integral value and error estimate, see IntegrationErrorKind::RoundoffErrorDetected.

• Bisection of the highest error region into two subregions results in subregions with integraion limits that are too small, see IntegrationErrorKind::BadIntegrandBehaviour.

• The integral does not converge. Occurs if at least 5 bisections with extrapolation have failed to improve the integral value and error estimates, see IntegrationErrorKind::DoesNotConverge.

• The integral is divergent or slowly convergent. Occurs if the ratio of the extrapolation improved integral estimate and the non-extrapolated estimate is too large or small, or if the calculated error is larger than the calculated value of the integral, see IntegrationErrorKind::DivergentOrSlowlyConverging.

• An error is encountered when initialising the integration workspace. This is an internal error, which should not occur downstream , see IntegrationErrorKind::UninitialisedWorkspace.

pub const fn limits(&self) -> Limits

Return the integration Limits.

impl<I> AdaptiveSingularity<InfiniteInterval<I>>

where I: Integrand<Point = f64>,

pub fn infinite( function: I, tolerance: Tolerance, max_iterations: usize, ) -> Result<Self, InitialisationError>

Generate a new AdaptiveSingularity integrator for functions with possible integrable singularities on infinite intervals.

Initialise an adaptive Gauss-Kronrod integrator with extrapolation for approximating integrals of the form, $$ \int_{-\infty}^{+\infty} f(x) dx. $$ with finite integration bounds $x \in (-\infty,+\infty)$.

Arguments:

• function: A user supplied function to be integrated which is something implementing the Integrand trait.
• tolerance: The tolerance requested by the user. Can be either an absolute tolerance or relative tolerance. Determines the exit condition of the integration routine, see Tolerance.
• max_iterations: The maximum number of iterations that the routine should use to try to satisfy the requested tolerance.

Errors

Function can return an error if it receives bad user input. This is primarily related to using invalid values for the tolerance. The returned error is an InitialisationError. See Tolerance and InitialisationError for more details.

Example

Here we present a calculation of [Euler’s constant] $\gamma$ (std::f64::consts::EULER_GAMMA) using the integral representation: $$ \gamma = -\int_{-\infty}^{+\infty} x e^{x - e^{x}} dx $$ which has an infinite integration region $x\in(-\infty,+\infty)$. We use AdaptiveSingularity::infinite to approximate $\gamma$.

 use rint::{Integrand, Tolerance};
 use rint::quadrature::AdaptiveSingularity;

 const GAMMA: f64 = std::f64::consts::EULER_GAMMA;
 const TOL: f64 = 1.0e-12;

 struct Gamma;

 impl Integrand for Gamma {
     type Point = f64;
     type Scalar = f64;

     fn evaluate(&self, x: &Self::Point) -> Self::Scalar {
         let exponent = x - x.exp();
         -x * exponent.exp()
     }
 }

 let gamma = Gamma;
 let tolerance = Tolerance::Relative(TOL);
 let integral = AdaptiveSingularity::infinite(gamma, tolerance, 1000)?
     .integrate()?;

 let result = integral.result();
 let error = integral.error();
 let abs_actual_error = (GAMMA - result).abs();
 let tol = TOL * result.abs();
 let iters = integral.iterations();
 assert_eq!(iters, 11);
 assert!(abs_actual_error < error);
 assert!(error < tol);

impl<I> AdaptiveSingularity<SemiInfiniteIntervalPositive<I>>

where I: Integrand<Point = f64>,

pub fn semi_infinite_upper( function: I, lower: f64, tolerance: Tolerance, max_iterations: usize, ) -> Result<Self, InitialisationError>

Generate a new AdaptiveSingularity integrator for functions with possible integrable singularities on semi-infinite intervals with an upper infinite limit.

Initialise an adaptive Gauss-Kronrod integrator with extrapolation for approximating integrals of the form, $$ \int_{b}^{+\infty} f(x) dx. $$ with finite integration bounds $x \in (b,+\infty)$.

Arguments:

• function: A user supplied function to be integrated which is something implementing the Integrand trait.
• lower: the lower integration bound.
• tolerance: The tolerance requested by the user. Can be either an absolute tolerance or relative tolerance. Determines the exit condition of the integration routine, see Tolerance.
• max_iterations: The maximum number of iterations that the routine should use to try to satisfy the requested tolerance.

Errors

Function can return an error if it receives bad user input. This is primarily related to using invalid values for the tolerance. The returned error is an InitialisationError. See Tolerance and InitialisationError for more details.

Example

Here we present a calculation of [Catalan’s constant] $G$ using the integral representation: $$ G = \frac{\pi}{2} \int_{1}^{+\infty} \frac{(x^{4} - 6 x^{2} + 1) \ln \ln x}{(1+x^{2})^{3}} d x $$ which has a semi-infinite integration region $x \in (1,+\infty)$ We use AdaptiveSingularity::semi_infinite_upper to approximate $G$.

 use rint::{Integrand, Limits, Tolerance};
 use rint::quadrature::AdaptiveSingularity;

 const G: f64 = 0.915_965_594_177_219_015_054_603_514_932_384_110_774;
 const TOL: f64 = 1.0e-12;

 struct Catalan;

 impl Integrand for Catalan {
     type Point = f64;
     type Scalar = f64;

     fn evaluate(&self, x: &Self::Point) -> Self::Scalar {
         let FRAC_PI_2 = std::f64::consts::FRAC_PI_2;
         let polynomial = x.powi(4) - 6.0 * x.powi(2) + 1.0;
         let denominator = (1.0 + x.powi(2)).powi(3);
         let lnlnx = x.ln().ln();
         
         FRAC_PI_2 * polynomial * lnlnx / denominator
     }
 }

 let catalan = Catalan;
 let lower = 1.0;
 let tolerance = Tolerance::Relative(TOL);
 let integral = AdaptiveSingularity::semi_infinite_upper(catalan, lower, tolerance, 1000)?
     .integrate()?;

 let result = integral.result();
 let error = integral.error();
 let abs_actual_error = (G - result).abs();
 let tol = TOL * result.abs();
 let iters = integral.iterations();
 assert_eq!(iters, 32);
 assert!(abs_actual_error < error);
 assert!(error < tol);

impl<I> AdaptiveSingularity<SemiInfiniteIntervalNegative<I>>

where I: Integrand<Point = f64>,

pub fn semi_infinite_lower( function: I, upper: f64, tolerance: Tolerance, max_iterations: usize, ) -> Result<Self, InitialisationError>

Generate a new AdaptiveSingularity integrator for functions with possible integrable singularities on semi-infinite intervals with an upper infinite limit.

Initialise an adaptive Gauss-Kronrod integrator with extrapolation for approximating integrals of the form, $$ \int_{-\infty}^{a} f(x) dx. $$ with finite integration bounds $x \in (-\infty,a)$.

Arguments:

• function: A user supplied function to be integrated which is something implementing the Integrand trait.
• upper: the upper integration bound.
• tolerance: The tolerance requested by the user. Can be either an absolute tolerance or relative tolerance. Determines the exit condition of the integration routine, see Tolerance.
• max_iterations: The maximum number of iterations that the routine should use to try to satisfy the requested tolerance.

Errors

Function can return an error if it receives bad user input. This is primarily related to using invalid values for the tolerance. The returned error is an InitialisationError. See Tolerance and InitialisationError for more details.

Example

Here we present a calculation of [Catalan’s constant] $G$ using the integral representation: $$ G = \frac{\pi}{2} \int^{-1}_{-\infty} \frac{(x^{4} - 6 x^{2} + 1) \ln \ln (-x)}{(1+x^{2})^{3}} d x $$ which has a semi-infinite integration region $x \in (1,+\infty)$. We use AdaptiveSingularity::semi_infinite_upper to approximate $G$.

 use rint::{Integrand, Limits, Tolerance};
 use rint::quadrature::AdaptiveSingularity;

 const G: f64 = 0.915_965_594_177_219_015_054_603_514_932_384_110_774;
 const TOL: f64 = 1.0e-12;

 struct Catalan;

 impl Integrand for Catalan {
     type Point = f64;
     type Scalar = f64;

     fn evaluate(&self, x: &Self::Point) -> Self::Scalar {
         let FRAC_PI_2 = std::f64::consts::FRAC_PI_2;
         let polynomial = x.powi(4) - 6.0 * x.powi(2) + 1.0;
         let denominator = (1.0 + x.powi(2)).powi(3);
         let lnlnx = (-x).ln().ln();
         
         FRAC_PI_2 * polynomial * lnlnx / denominator
     }
 }
 let catalan = Catalan;
 let upper = -1.0;
 let tolerance = Tolerance::Relative(TOL);
 let integral = AdaptiveSingularity::semi_infinite_lower(catalan, upper, tolerance, 1000)?
     .integrate()?;

 let result = integral.result();
 let error = integral.error();
 let abs_actual_error = (G - result).abs();
 let tol = TOL * result.abs();
 let iters = integral.iterations();
 assert_eq!(iters, 32);
 assert!(abs_actual_error < error);
 assert!(error < tol);

Trait Implementations

impl<I: Clone> Clone for AdaptiveSingularity<I>

fn clone(&self) -> AdaptiveSingularity<I>

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl<I: Debug> Debug for AdaptiveSingularity<I>

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl<I: PartialEq> PartialEq for AdaptiveSingularity<I>

fn eq(&self, other: &AdaptiveSingularity<I>) -> bool

Tests for self and other values to be equal, and is used by ==.

fn ne(&self, other: &Rhs) -> bool

Tests for !=. The default implementation is almost always sufficient, and should not be overridden without very good reason.

impl<I: PartialOrd> PartialOrd for AdaptiveSingularity<I>

fn partial_cmp(&self, other: &AdaptiveSingularity<I>) -> Option<Ordering>

This method returns an ordering between self and other values if one exists.

fn lt(&self, other: &Rhs) -> bool

Tests less than (for self and other) and is used by the < operator.

fn le(&self, other: &Rhs) -> bool

Tests less than or equal to (for self and other) and is used by the <= operator.

fn gt(&self, other: &Rhs) -> bool

Tests greater than (for self and other) and is used by the > operator.

fn ge(&self, other: &Rhs) -> bool

Tests greater than or equal to (for self and other) and is used by the >= operator.

impl<I: Copy> Copy for AdaptiveSingularity<I>

impl<I> StructuralPartialEq for AdaptiveSingularity<I>