numerical_differentiation/
lib.rs

1use cached::proc_macro::cached;
2use nalgebra::{DMatrix, DVector, RowDVector};
3use rayon::iter::{IndexedParallelIterator, IntoParallelIterator, ParallelIterator};
4
5#[cached]
6fn factorial(n: u128) -> u128 {
7    match n {
8        0 => 1,
9        n => n * factorial(n - 1),
10    }
11}
12
13/// https://en.wikipedia.org/wiki/Finite_difference_coefficient#Arbitrary_stencil_points
14#[cached]
15fn get_finite_difference_coefficients(stencils: Vec<i128>, d: u128) -> Vec<f64> {
16    let n = stencils.len();
17
18    let stencils = RowDVector::from_iterator(n, stencils.iter().map(|&i| i as f64));
19
20    let mut stencil_matrix = DMatrix::from_element(n, n, 0f64);
21    let mut pow = RowDVector::from_element(n, 1.);
22
23    stencil_matrix.set_row(0, &pow);
24    for mut row in stencil_matrix.row_iter_mut().skip(1) {
25        pow = pow.component_mul(&stencils);
26        row.copy_from(&pow);
27    }
28
29    stencil_matrix.try_inverse_mut();
30
31    let mut factorial_column = DVector::from_element(n, 0.);
32    factorial_column[d as usize] = factorial(d) as f64;
33
34    let coefficients = stencil_matrix * factorial_column;
35
36    coefficients.data.as_vec().to_owned()
37}
38
39pub fn get_epsilon(x: f64, d: u128) -> f64 {
40    x * f64::EPSILON.powf(1. / (d as f64 + 1.))
41}
42
43pub fn differentiate_with_stencils<F: Fn(f64) -> f64 + Sync>(
44    x: f64,
45    f: F,
46    d: u128,
47    stencils: Vec<i128>,
48) -> f64 {
49    assert!((d as usize) < stencils.len());
50
51    let epsilon = get_epsilon(x, d);
52    let coefficients = get_finite_difference_coefficients(stencils.clone(), d);
53    let numerator: f64 = stencils
54        .into_par_iter()
55        .zip(coefficients.into_par_iter())
56        .map(|(stencil, coefficient)| coefficient * f(x + stencil as f64 * epsilon))
57        .sum();
58    let denominator = epsilon.powi(d as i32);
59
60    numerator / denominator
61}
62
63#[cached]
64fn get_stencils(d: u128) -> Vec<i128> {
65    let start = -(d as i128) / 2;
66    let end = (d as i128 + 1) / 2;
67
68    (start..=end).collect()
69}
70
71pub fn differentiate<F: Fn(f64) -> f64 + Sync>(x: f64, f: F, d: u128) -> f64 {
72    let stencils = get_stencils(d);
73
74    differentiate_with_stencils(x, f, d, stencils)
75}