conspire/math/special/
mod.rs

1#[cfg(test)]
2mod test;
3
4use crate::{ABS_TOL, REL_TOL};
5use std::f64::consts::E;
6
7const NINE_FIFTHS: f64 = 9.0 / 5.0;
8
9/// Returns the inverse Langevin function.
10///
11/// ```math
12/// x = \mathcal{L}^{-1}(y)
13/// ```
14/// The first few terms of the Maclaurin series are used when $`|y|<3e^{-3}`$.
15/// ```math
16/// x \sim 3y + \frac{9}{5}y^3 + \mathrm{ord}(y^5)
17/// ```
18/// Two iterations of Newton's method are used to improve upon an initial guess given by [inverse_langevin_approximate()] otherwise.
19/// The resulting maximum relative error is below $`1e^{-12}`$.
20pub fn inverse_langevin(y: f64) -> f64 {
21    let y_abs = y.abs();
22    if y_abs >= 1.0 {
23        panic!()
24    } else if y_abs <= 3e-3 {
25        3.0 * y + NINE_FIFTHS * y.powi(3)
26    } else {
27        let mut x = inverse_langevin_approximate(y_abs);
28        for _ in 0..2 {
29            x += (y_abs - langevin(x)) / langevin_derivative(x);
30        }
31        if y < 0.0 { -x } else { x }
32    }
33}
34
35/// Returns an approximation of the inverse Langevin function.[^cite]
36///
37/// [^cite]: R. Jedynak, [Math. Mech. Solids **24**, 1992 (2019)](https://doi.org/10.1177/1081286518811395).
38///
39/// ```math
40/// \mathcal{L}^{-1}(y) \approx \frac{2.14234 y^3 - 4.22785 y^2 + 3y}{(1 - y)(0.71716 y^3 - 0.41103 y^2 - 0.39165 y + 1)}
41/// ```
42/// This approximation has a maximum relative error of $`8.2e^{-4}`$.
43pub fn inverse_langevin_approximate(y: f64) -> f64 {
44    (2.14234 * y.powi(3) - 4.22785 * y.powi(2) + 3.0 * y)
45        / (1.0 - y)
46        / (0.71716 * y.powi(3) - 0.41103 * y.powi(2) - 0.39165 * y + 1.0)
47}
48
49/// Returns the Lambert W function.
50///
51/// ```math
52/// y = W_0(x)
53/// ```
54pub fn lambert_w(x: f64) -> f64 {
55    if x == -1.0 / E {
56        -1.0
57    } else if x == 0.0 {
58        0.0
59    } else if x == E {
60        1.0
61    } else if x < -1.0 / E {
62        panic!()
63    } else {
64        let mut w = if x < 0.0 {
65            (E * x * (1.0 + (1.0 + E * x).sqrt()).ln()) / (1.0 + E * x + (1.0 + E * x).sqrt())
66        } else if x < E {
67            x / E
68        } else {
69            x.ln() - x.ln().ln()
70        };
71        let mut error = w * w.exp() - x;
72        while error.abs() >= ABS_TOL && (error / x).abs() >= REL_TOL {
73            w *= (1.0 + (x / w).ln()) / (1.0 + w);
74            error = w * w.exp() - x;
75        }
76        w
77    }
78}
79
80/// Returns the Langevin function.
81///
82/// ```math
83/// \mathcal{L}(x) = \coth(x) - x^{-1}
84/// ```
85pub fn langevin(x: f64) -> f64 {
86    if x == 0.0 {
87        0.0
88    } else {
89        1.0 / x.tanh() - 1.0 / x
90    }
91}
92
93/// Returns derivative of the Langevin function.
94///
95/// ```math
96/// \mathcal{L}'(x) = x^{-2} - \sinh^{-2}(x)
97/// ```
98pub fn langevin_derivative(x: f64) -> f64 {
99    1.0 / x.powi(2) - 1.0 / x.sinh().powi(2)
100}
conspire/math/special/mod.rs

conspire/math/special/
mod.rs
