Crate gauss_quad

Source
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, e.g.:

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;

// initialize the quadrature rule
let degree = 10;
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);

§Setting up a quadrature rule

Using a quadrature rule takes two steps:

  1. Initialization
  2. 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);

§Specific quadrature rules

Different rules may take different parameters.

For example, the GaussLaguerre rule requires both a degree and an alpha parameter.

GaussLaguerre is also defined as an improper integral over the domain [0, ∞). This means no domain bounds are needed in the integrate call.

// initialize the quadrature rule
let degree = 2;
let alpha = 0.5;
let quad = GaussLaguerre::new(degree, alpha)?;

// use the rule to integrate a function
let integral = quad.integrate(|x| x * x);

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

§Errors

Quadrature rules are only defined for a certain set of input values. For example, every rule is only defined for degrees where degree > 1.

let degree = 1;
assert!(GaussLaguerre::new(degree, 0.1).is_err());

Specific rules may have other requirements. GaussJacobi for example, requires alpha and beta parameters larger than -1.0.

let degree = 10;
let alpha = 0.1;
let beta = -1.1;

assert!(GaussJacobi::new(degree, alpha, beta).is_err())

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 Fn(f64) -> f64 trait, i.e. functions of one f64 parameter.


// initialize the quadrature rule
let degree = 2;
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);

§Double integrals

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

let rule = GaussLegendre::new(3)?;

// 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.

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).

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.

Structs§

Type Aliases§

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