cheby 0.3.0

Unit-safe Chebyshev approximation and spectral numerics for Rust.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69

//! Coefficient fitting from samples and functions.
//!
//! The fixed-size fitter samples at roots of `T_N` and applies the direct
//! discrete-cosine coefficient formula. With `T_j(x_k) = cos(j θ_k)` evaluated
//! by Chebyshev recurrence in `j`, the work is `O(N²)` multiplies-and-adds
//! and only `O(N)` calls to `cos` (used to build the node grid once).

use crate::core::{nodes, ChebyScalar, ChebySeries, ChebySeriesOn, Domain, NodeKind};

/// Compute Chebyshev coefficients from values sampled at roots of `T_N`.
///
/// Uses the Chebyshev three-term recurrence to generate `T_j(x_k)` for each
/// sample `k` in a single pass over `j`, avoiding `N²` `cos()` calls.
#[inline]
pub fn fit_coeffs<T: ChebyScalar, const N: usize>(values: &[T; N]) -> [T; N] {
    let mut coeffs = [T::zero(); N];
    if N == 0 {
        return coeffs;
    }
    let xs = nodes::<N>(NodeKind::Roots);
    for (k, &v) in values.iter().enumerate() {
        let x = xs[k];
        // T_0 = 1
        coeffs[0] = coeffs[0] + v;
        if N > 1 {
            // T_1 = x
            coeffs[1] = coeffs[1] + v * x;
            let two_x = 2.0 * x;
            let mut tkm1 = 1.0_f64;
            let mut tk = x;
            for coeff in coeffs.iter_mut().take(N).skip(2) {
                let tkp1 = two_x * tk - tkm1;
                *coeff = *coeff + v * tkp1;
                tkm1 = tk;
                tk = tkp1;
            }
        }
    }
    let nf = N as f64;
    coeffs[0] = coeffs[0] / nf;
    let scale = 2.0 / nf;
    for coeff in coeffs.iter_mut().skip(1) {
        *coeff = *coeff * scale;
    }
    coeffs
}

/// Fit a Chebyshev series on normalized `[-1, 1]`.
#[inline]
pub fn fit_from_fn<T: ChebyScalar, const N: usize>(f: impl Fn(f64) -> T) -> ChebySeries<T, N> {
    let xs = nodes::<N>(NodeKind::Roots);
    let values = core::array::from_fn(|k| f(xs[k]));
    ChebySeries::new(fit_coeffs(&values))
}

/// Fit a Chebyshev series on a physical domain.
#[inline]
pub fn fit_from_fn_on<T, X, const N: usize>(
    domain: Domain<X>,
    f: impl Fn(X) -> T,
) -> ChebySeriesOn<T, X, N>
where
    T: ChebyScalar,
    X: crate::core::ChebyTime,
{
    let xs = crate::core::nodes::nodes_mapped::<X, N>(domain, NodeKind::Roots);
    let values = core::array::from_fn(|k| f(xs[k]));
    ChebySeriesOn::new(domain, ChebySeries::new(fit_coeffs(&values)))
}