thermolib/
algorithms.rs

1/// Reference:
2/// Richard Brent,
3/// Algorithms for Minimization Without Derivatives.
4pub fn brent_zero<F>(mut f: F, a: f64, b: f64) -> f64
5where
6    F: FnMut(f64) -> f64,
7{
8    let macheps = f64::EPSILON.sqrt();
9    let (mut a, mut fa) = (a, f(a));
10    let (mut b, mut fb) = (b, f(b));
11    if (fa * fb) > 0.0 {
12        return f64::NAN;
13    }
14    let (mut c, mut fc) = (a, fa);
15    let (mut d, mut e) = (b - a, b - a);
16    let (mut tol, mut m, mut p, mut q, mut r, mut s);
17    loop {
18        if fc.abs() < fb.abs() {
19            (a, fa) = (b, fb);
20            (b, fb) = (c, fc);
21            (c, fc) = (a, fa);
22        }
23        tol = macheps * (2.0 * b.abs() + 10.0);
24        m = 0.5 * (c - b);
25        if m.abs() <= tol || fb == 0.0 {
26            return b;
27        }
28        if e.abs() < tol || fa.abs() <= fb.abs() {
29            e = m;
30            d = e;
31        } else {
32            s = fb / fa;
33            if a == c {
34                p = 2.0 * m * s;
35                q = 1.0 - s;
36            } else {
37                q = fa / fc;
38                r = fb / fc;
39                p = s * (2.0 * m * q * (q - r) - (b - a) * (r - 1.0));
40                q = (q - 1.0) * (r - 1.0) * (s - 1.0);
41            }
42            if 0.0 < p {
43                q = -q;
44            } else {
45                p = -p;
46            }
47            s = e;
48            e = d;
49            if 2.0 * p < 3.0 * m * q - (tol * q).abs() && p < (0.5 * s * q).abs() {
50                d = p / q;
51            } else {
52                e = m;
53                d = e;
54            }
55        }
56        (a, fa) = (b, fb);
57        b += if tol < d.abs() {
58            d
59        } else if 0.0 < m {
60            tol
61        } else {
62            -tol
63        };
64        fb = f(b);
65        if (0.0 < fb && 0.0 < fc) || (fb.is_sign_negative() && fc.is_sign_negative()) {
66            (c, fc) = (a, fa);
67            e = b - a;
68            d = e;
69        }
70    }
71}
72/// Romberg Numerical Differentiation
73/// Reference:
74/// C. J. F. RIDDERS
75/// Accurate computation of F'(x) and F'(x) F"(x)
76/// Adv. Eng. Software, 1982, Vol. 4, No. 2, 75-76.
77const N_DIM: usize = 11;
78pub fn romberg_diff<F>(mut f: F, x: f64) -> f64
79where
80    F: FnMut(f64) -> f64,
81{
82    let tol = f64::EPSILON.sqrt() * 2.0;
83    let mut fx: [[f64; N_DIM]; N_DIM] = [[0.0; N_DIM]; N_DIM];
84    let mut h: f64 = if x > 1.0 { 0.01 * x } else { 0.1 * x };
85    fx[0][0] = (f(x + h) - f(x - h)) / 2.0 / h;
86    for i in 1..N_DIM {
87        h /= 2.0;
88        fx[0][i] = (f(x + h) - f(x - h)) / 2.0 / h;
89        for j in 1..i + 1 {
90            fx[j][i] = (fx[j - 1][i] * 4_f64.powi(j as i32) - fx[j - 1][i - 1])
91                / (4_f64.powi(j as i32) - 1.0);
92        }
93        if (fx[i][i] / fx[i - 1][i - 1] - fx[i - 1][i - 1] / fx[i][i]).abs() < tol {
94            return fx[i][i];
95        }
96    }
97    fx[N_DIM - 1][N_DIM - 1]
98}
99/// unit test
100#[cfg(test)]
101mod tests {
102    use super::*;
103    #[test]
104    fn test_algorithms() {
105        // test brent_zero
106        let x0 = brent_zero(|x: f64| x.sin() - x / 2.0, 1.0, 2.0);
107        let (x1, digits) = (1.8954942796e+00, 10_f64.powi(10));
108        assert_eq!(x1, (x0 * digits).round() / digits);
109        let x0 = brent_zero(|x: f64| 2.0 * x - (-x).exp(), 0.0, 1.0);
110        let (x1, digits) = (3.5173365631e-01, 10_f64.powi(11));
111        assert_eq!(x1, (x0 * digits).round() / digits);
112        let x0 = brent_zero(|x: f64| x * (-x).exp(), -1.0, 0.5);
113        let (x1, digits) = (-4.0352160429e-10, 10_f64.powi(20));
114        assert_eq!(x1, (x0 * digits).round() / digits);
115        let x0 = brent_zero(|x: f64| x.exp() - 1.0 / 100.0 / x / x, 0.0001, 20.0);
116        let (x1, digits) = (9.5344620258e-02, 10_f64.powi(12));
117        assert_eq!(x1, (x0 * digits).round() / digits);
118        let x0 = brent_zero(|x: f64| (x + 3.0) * (x - 1.0) * (x - 1.0), -5.0, 2.0);
119        let (x1, digits) = (-3.0000000000e+00, 10_f64.powi(10));
120        assert_eq!(x1, (x0 * digits).round() / digits);
121        // test romberg_diff
122        let (df0, digits) = (140.7377356, 10_f64.powi(7)); // real = 140.7377355
123        let df1 = romberg_diff(|x: f64| x.exp() / (x.sin() - x.powi(2)), 1.0);
124        assert_eq!(df0, (df1 * digits).round() / digits);
125        let (df0, digits) = (2.718281828, 10_f64.powi(9));
126        let df1 = romberg_diff(|x: f64| x.exp(), 1.0);
127        assert_eq!(df0, (df1 * digits).round() / digits);
128        let (df0, digits) = (22026.46579, 10_f64.powi(5));
129        let df1 = romberg_diff(|x: f64| x.exp(), 10.0);
130        assert_eq!(df0, (df1 * digits).round() / digits);
131    }
132}
thermolib/algorithms.rs

thermolib/
algorithms.rs
