Crate gauss_quad

Expand description

§gauss-quad

gauss-quad is a Gaussian quadrature library for numerical integration.

§Quadrature rules

gauss-quad implements the following quadrature rules:

§Using gauss-quad

To use any of the quadrature rules in your project, first initialize the rule with a specified degree and then you can use it for integration. The degree is a non-zero positive integer that determines the number of nodes used in the quadrature rule, and as such the number of points at which the integrand is evaluated.

use gauss_quad::GaussLegendre;
// This macro is used in these docs to compare floats.
// The assertion succeeds if the two sides are within floating point error,
// or an optional epsilon.
use approx::assert_abs_diff_eq;
use core::num::NonZeroUsize;

// initialize the quadrature rule
let degree = NonZeroUsize::new(10).unwrap();
let quad = GaussLegendre::new(degree);

// Use the rule to integrate a function
let left_bound = 0.0;
let right_bound = 1.0;
let integral = quad.integrate(left_bound, right_bound, |x| x * x);
assert_abs_diff_eq!(integral, 1.0 / 3.0);

Select the degree, n, such that 2n-1 is the largest degree of polynomial that you want to integrate with the rule.

§Setting up a quadrature rule

Using a quadrature rule takes two steps:

Initialization
Integration

First, rules must be initialized using some specific input parameters.

Then, you can integrate functions using those rules:

let gauss_legendre = GaussLegendre::new(degree);
// Integrate on the domain [a, b]
let x_cubed = gauss_legendre.integrate(a, b, |x| x * x * x);

let gauss_jacobi = GaussJacobi::new(degree, alpha, beta);
// Integrate on the domain [c, d]
let double_x = gauss_jacobi.integrate(c, d, |x| 2.0 * x);

let gauss_laguerre = GaussLaguerre::new(degree, alpha);
// No explicit domain, Gauss-Laguerre integration is done on the domain [0, ∞).
let piecewise = gauss_laguerre.integrate(|x| if x > 0.0 && x < 2.0 { x } else { 0.0 });

let gauss_hermite = GaussHermite::new(degree);
// Again, no explicit domain since Gauss-Hermite integration is done over the domain (-∞, ∞).
let golden_polynomial = gauss_hermite.integrate(|x| x * x - x - 1.0);

§Different rules have different requirements

Quadrature rules are only defined for a certain set of input values. They take parameters of types that guarantee valid input.

GaussJacobi for example, requires two parameters named “alpha” and “beta”, and they must be larger than -1.0. Therefore they are of type FiniteAboveNegOneF64, which can only be created from finite non-NAN values above -1.0:

use gauss_quad::{GaussJacobi, FiniteAboveNegOneF64};
use core::num::NonZeroUsize;

let degree = NonZeroUsize::new(10).unwrap();
let alpha = FiniteAboveNegOneF64::new(0.1).unwrap();

// Trying to create a `beta` below -1.0 results in `None`.
let beta = FiniteAboveNegOneF64::new(-1.1);
assert!(beta.is_none());

// This is valid:
let beta = FiniteAboveNegOneF64::new(-0.5).unwrap();
let quad = GaussJacobi::new(degree, alpha, beta);

GaussLaguerre is used to evaluate an improper integral over the domain [0, ∞). This means no domain bounds are needed in the integrate call.

// Initialize the quadrature rule
let degree = NonZeroUsize::new(2).unwrap();
let alpha = FiniteAboveNegOneF64::new(0.5).unwrap();
let quad = GaussLaguerre::new(degree, alpha);

// Use the rule to integrate a function, note the lack of domain bounds
let integral = quad.integrate(|x| x * x);

assert_abs_diff_eq!(integral, 15.0 * PI.sqrt() / 8.0, epsilon = 1e-14);

Make sure to read the specific quadrature rule’s documentation before using it.

§Passing functions to quadrature rules

The integrate method takes any integrand that implements the FnMut(f64) -> f64 trait, i.e. functions of one f64 parameter that can modify internal state.


// Initialize the quadrature rule
let degree = NonZeroUsize::new(2).unwrap();
let quad = GaussLegendre::new(degree);

// Use the rule to integrate a closure
let left_bound = 0.0;
let right_bound = 1.0;

let integral = quad.integrate(left_bound, right_bound, |x| x * x);

assert_abs_diff_eq!(integral, 1.0 / 3.0);

// You can also pass a function pointer
fn x_cubed(x: f64) -> f64 {
   x * x * x
}

let integral_x_cubed = quad.integrate(left_bound, right_bound, x_cubed);
assert_abs_diff_eq!(integral_x_cubed, 1.0 / 4.0);

The parallel par_integrate methods instead take a function that can’t modify any state, Fn(f64) -> f64, and is Sync to enable parallel evaluation.

§Double integrals

It is possible to use this crate to do double and higher integrals:

let rule = GaussLegendre::new(3.try_into().unwrap());

// integrate x^2*y over the triangle in the xy-plane where x ϵ [0, 1] and y ϵ [0, x]:
let double_int = rule.integrate(0.0, 1.0, |x| rule.integrate(0.0, x, |y| x * x * y));

assert_abs_diff_eq!(double_int, 0.1);

However, the time complexity of the integration then scales with the number of nodes to the power of the depth of the integral, e.g. O(n³) for triple integrals.

§Feature flags

serde: implements the Serialize and Deserialize traits from the serde crate for the quadrature rule structs.

rkyv: implements the Serialize, Deserialize, and Archive traits from the rkyv crate.

rayon: enables a parallel version of the integrate function on the quadrature rule structs. Can speed up integration if evaluating the integrand is expensive (takes ≫100 µs). The rayon crate depends on the standard library, so this also enables the std feature.

zerocopy: imlements the KnownLayout trait from the zerocopy crate for the quadrature rule structs.

One of the below features must be enabled:

libm (enabled by default): depends on the libm crate and uses it as the math backend. Does nothing if the std feature is enabled.

std: links the standard library and uses it as the math backend. When this feature is disabled the crate is no_std compatible.

Modules§

chebyshev: Numerical integration using the Gauss-Chebyshev quadrature rule.
hermite: Numerical integration using the Gauss-Hermite quadrature rule.
jacobi: Numerical integration using the Gauss-Jacobi quadrature rule.
laguerre: Numerical integration using the generalized Gauss-Laguerre quadrature rule.
legendre: Numerical integration using the Gauss-Legendre quadrature rule.
midpoint: Numerical integration using the midpoint rule.
simpson: Numerical integration using a Simpson’s rule.
trapezoid: Numerical integration using the trapezoid rule.

Structs§

FiniteAboveNegOneF64: A wrapper around an f64 that ensures the value is greater than -1.0, finite, and not NAN.
GaussChebyshevFirstKind: A Gauss-Chebyshev quadrature scheme of the first kind.
GaussChebyshevSecondKind: A Gauss-Chebyshev quadrature scheme of the second kind.
GaussHermite: A Gauss-Hermite quadrature scheme.
GaussJacobi: A Gauss-Jacobi quadrature scheme.
GaussLaguerre: A Gauss-Laguerre quadrature scheme.
GaussLegendre: A Gauss-Legendre quadrature scheme.
InfNanNegOneOrLessError: The error that that is returned when trying to convert a f64 value that is less than or equal to -1.0 into a FiniteAboveNegOneF64 with the TryFrom trait.
Midpoint: A midpoint rule.
Simpson: A Simpson’s rule.
Trapezoid: A trapezoid rule.

Enums§

ParseFiniteAboveNegOneF64Error: An error that can occur when parsing a &str into a FiniteAboveNegOneF64.

Type Aliases§

Node: A node in a quadrature rule.
Weight: A weight in a quadrature rule.

Crate gauss_quad

Crate gauss_quad Copy item path

§gauss-quad

§Quadrature rules

§Using gauss-quad

§Setting up a quadrature rule

§Different rules have different requirements

§Passing functions to quadrature rules

§Double integrals

§Feature flags

Modules§

Structs§

Enums§

Type Aliases§

Crate gauss_quad