gauss_quad/lib.rs
1//! # gauss-quad
2//!
3//! **gauss-quad** is a [Gaussian quadrature](https://en.wikipedia.org/wiki/Gaussian_quadrature) library for numerical integration.
4//!
5//! ## Quadrature rules
6//!
7//! **gauss-quad** implements the following quadrature rules:
8//! * [Gauss-Legendre](https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_quadrature)
9//! * [Gauss-Jacobi](https://en.wikipedia.org/wiki/Gauss%E2%80%93Jacobi_quadrature)
10//! * [Gauss-Laguerre](https://en.wikipedia.org/wiki/Gauss%E2%80%93Laguerre_quadrature) (generalized)
11//! * [Gauss-Hermite](https://en.wikipedia.org/wiki/Gauss%E2%80%93Hermite_quadrature)
12//! * [Gauss-Chebyshev](https://en.wikipedia.org/wiki/Chebyshev%E2%80%93Gauss_quadrature)
13//! * [Midpoint](https://en.wikipedia.org/wiki/Riemann_sum#Midpoint_rule)
14//! * [Simpson](https://en.wikipedia.org/wiki/Simpson%27s_rule)
15//!
16//! ## Using **gauss-quad**
17//!
18//! To use any of the quadrature rules in your project, first initialize the rule with
19//! a specified degree and then you can use it for integration, e.g.:
20//! ```
21//! use gauss_quad::GaussLegendre;
22//! # use gauss_quad::legendre::GaussLegendreError;
23//! // This macro is used in these docs to compare floats.
24//! // The assertion succeeds if the two sides are within floating point error,
25//! // or an optional epsilon.
26//! use approx::assert_abs_diff_eq;
27//!
28//! // initialize the quadrature rule
29//! let degree = 10;
30//! let quad = GaussLegendre::new(degree)?;
31//!
32//! // use the rule to integrate a function
33//! let left_bound = 0.0;
34//! let right_bound = 1.0;
35//! let integral = quad.integrate(left_bound, right_bound, |x| x * x);
36//! assert_abs_diff_eq!(integral, 1.0 / 3.0);
37//! # Ok::<(), GaussLegendreError>(())
38//! ```
39//!
40//! Select the degree, n, such that 2n-1 is the largest degree of polynomial that you want to integrate with the rule.
41//!
42//! ## Setting up a quadrature rule
43//!
44//! Using a quadrature rule takes two steps:
45//! 1. Initialization
46//! 2. Integration
47//!
48//! First, rules must be initialized using some specific input parameters.
49//!
50//! Then, you can integrate functions using those rules:
51//! ```
52//! # use gauss_quad::*;
53//! # let degree = 5;
54//! # let alpha = 1.2;
55//! # let beta = 1.2;
56//! # let a = 0.0;
57//! # let b = 1.0;
58//! # let c = -10.;
59//! # let d = 100.;
60//! let gauss_legendre = GaussLegendre::new(degree)?;
61//! // Integrate on the domain [a, b]
62//! let x_cubed = gauss_legendre.integrate(a, b, |x| x * x * x);
63//!
64//! let gauss_jacobi = GaussJacobi::new(degree, alpha, beta)?;
65//! // Integrate on the domain [c, d]
66//! let double_x = gauss_jacobi.integrate(c, d, |x| 2.0 * x);
67//!
68//! let gauss_laguerre = GaussLaguerre::new(degree, alpha)?;
69//! // no explicit domain, Gauss-Laguerre integration is done on the domain [0, ∞).
70//! let piecewise = gauss_laguerre.integrate(|x| if x > 0.0 && x < 2.0 { x } else { 0.0 });
71//!
72//! let gauss_hermite = GaussHermite::new(degree)?;
73//! // again, no explicit domain since Gauss-Hermite integration is done over the domain (-∞, ∞).
74//! let golden_polynomial = gauss_hermite.integrate(|x| x * x - x - 1.0);
75//! # Ok::<(), Box<dyn std::error::Error>>(())
76//! ```
77//!
78//! ## Specific quadrature rules
79//!
80//! Different rules may take different parameters.
81//!
82//! For example, the `GaussLaguerre` rule requires both a `degree` and an `alpha`
83//! parameter.
84//!
85//! `GaussLaguerre` is also defined as an improper integral over the domain [0, ∞).
86//! This means no domain bounds are needed in the `integrate` call.
87//! ```
88//! # use gauss_quad::laguerre::{GaussLaguerre, GaussLaguerreError};
89//! # use approx::assert_abs_diff_eq;
90//! # use core::f64::consts::PI;
91//! // initialize the quadrature rule
92//! let degree = 2;
93//! let alpha = 0.5;
94//! let quad = GaussLaguerre::new(degree, alpha)?;
95//!
96//! // use the rule to integrate a function
97//! let integral = quad.integrate(|x| x * x);
98//!
99//! assert_abs_diff_eq!(integral, 15.0 * PI.sqrt() / 8.0, epsilon = 1e-14);
100//! # Ok::<(), GaussLaguerreError>(())
101//! ```
102//!
103//! ## Errors
104//!
105//! Quadrature rules are only defined for a certain set of input values.
106//! For example, every rule is only defined for degrees where `degree > 1`.
107//! ```
108//! # use gauss_quad::GaussLaguerre;
109//! let degree = 1;
110//! assert!(GaussLaguerre::new(degree, 0.1).is_err());
111//! ```
112//!
113//! Specific rules may have other requirements.
114//! `GaussJacobi` for example, requires alpha and beta parameters larger than -1.0.
115//! ```
116//! # use gauss_quad::jacobi::GaussJacobi;
117//! let degree = 10;
118//! let alpha = 0.1;
119//! let beta = -1.1;
120//!
121//! assert!(GaussJacobi::new(degree, alpha, beta).is_err())
122//! ```
123//! Make sure to read the specific quadrature rule's documentation before using it.
124//!
125//! ## Passing functions to quadrature rules
126//!
127//! The `integrate` method takes any integrand that implements the [`FnMut(f64) -> f64`](FnMut) trait, i.e. functions of
128//! one `f64` parameter.
129//!
130//! ```
131//! # use gauss_quad::legendre::{GaussLegendre, GaussLegendreError};
132//! # use approx::assert_abs_diff_eq;
133//!
134//! // initialize the quadrature rule
135//! let degree = 2;
136//! let quad = GaussLegendre::new(degree)?;
137//!
138//! // use the rule to integrate a function
139//! let left_bound = 0.0;
140//! let right_bound = 1.0;
141//!
142//! let integral = quad.integrate(left_bound, right_bound, |x| x * x);
143//!
144//! assert_abs_diff_eq!(integral, 1.0 / 3.0);
145//! # Ok::<(), GaussLegendreError>(())
146//! ```
147//!
148//! ## Double integrals
149//!
150//! It is possible to use this crate to do double and higher integrals:
151//! ```
152//! # use gauss_quad::legendre::{GaussLegendre, GaussLegendreError};
153//! # use approx::assert_abs_diff_eq;
154//! let rule = GaussLegendre::new(3)?;
155//!
156//! // integrate x^2*y over the triangle in the xy-plane where x ϵ [0, 1] and y ϵ [0, x]:
157//! let double_int = rule.integrate(0.0, 1.0, |x| rule.integrate(0.0, x, |y| x * x * y));
158//!
159//! assert_abs_diff_eq!(double_int, 0.1);
160//! # Ok::<(), GaussLegendreError>(())
161//! ```
162//! However, the time complexity of the integration then scales with the number of nodes to
163//! the power of the depth of the integral, e.g. O(n³) for triple integrals.
164//!
165//! ## Feature flags
166//!
167//! `serde`: implements the [`Serialize`](serde::Serialize) and [`Deserialize`](serde::Deserialize) traits from
168//! the [`serde`](https://crates.io/crates/serde) crate for the quadrature rule structs.
169//!
170//! `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).
171
172// Only enable the nighlty `doc_auto_cfg` feature when
173// the `docsrs` configuration attribute is defined.
174// The config in Cargo.toml means that this is defined when we are on docs.rs (which uses the nightly compiler)
175// or if the environment variable RUSTDOCFLAGS is set as `RUSTDOCFLAGS="--cfg docsrs"`.
176// This lets us get a tag on docs.rs that says "Available on crate feature ... only".
177#![cfg_attr(docsrs, feature(doc_auto_cfg))]
178
179use nalgebra::{Dyn, Matrix, VecStorage};
180
181type DMatrixf64 = Matrix<f64, Dyn, Dyn, VecStorage<f64, Dyn, Dyn>>;
182
183pub mod chebyshev;
184mod data_api;
185mod gamma;
186pub mod hermite;
187pub mod jacobi;
188pub mod laguerre;
189pub mod legendre;
190pub mod midpoint;
191pub mod simpson;
192
193#[doc(inline)]
194pub use chebyshev::{GaussChebyshevFirstKind, GaussChebyshevSecondKind};
195#[doc(inline)]
196pub use data_api::{Node, Weight};
197#[doc(inline)]
198pub use hermite::GaussHermite;
199#[doc(inline)]
200pub use jacobi::GaussJacobi;
201#[doc(inline)]
202pub use laguerre::GaussLaguerre;
203#[doc(inline)]
204pub use legendre::GaussLegendre;
205#[doc(inline)]
206pub use midpoint::Midpoint;
207#[doc(inline)]
208pub use simpson::Simpson;