bacon_sci_1/interp/
mod.rs

1/* This file is part of bacon.
2 * Copyright (c) Wyatt Campbell.
3 *
4 * See repository LICENSE for information.
5 */
6
7use super::polynomial::Polynomial;
8use nalgebra::ComplexField;
9
10mod spline;
11pub use spline::*;
12
13/// Create a Lagrange interpolating polynomial.
14///
15/// Create an nth degree polynomial matching the n points (xs[i], ys[i])
16/// using Neville's iterated method for Lagrange polynomials. The result will
17/// match no derivatives.
18///
19/// # Examples
20/// ```
21/// use bacon_sci::interp::lagrange;
22/// use bacon_sci::polynomial::Polynomial;
23/// fn example() {
24///     let xs: Vec<_> = (0..10).map(|i| i as f64).collect();
25///     let ys: Vec<_> = xs.iter().map(|x| x.cos()).collect();
26///     let poly = lagrange(&xs, &ys, 1e-6).unwrap();
27///     for x in xs {
28///         assert!((x.cos() - poly.evaluate(x)).abs() < 0.00001);
29///     }
30/// }
31/// ```
32pub fn lagrange<N: ComplexField>(
33    xs: &[N],
34    ys: &[N],
35    tol: N::RealField,
36) -> Result<Polynomial<N>, String> {
37    if xs.len() != ys.len() {
38        return Err("lagrange: slices have mismatched dimension".to_owned());
39    }
40
41    let mut qs = vec![Polynomial::with_tolerance(tol)?; xs.len() * xs.len()];
42    for (ind, y) in ys.iter().enumerate() {
43        qs[ind] = polynomial![*y];
44    }
45
46    for i in 1..xs.len() {
47        let mut poly_2 = polynomial![N::one(), -xs[i]];
48        poly_2.set_tolerance(tol)?;
49        for j in 1..=i {
50            let mut poly_1 = polynomial![N::one(), -xs[i - j]];
51            poly_1.set_tolerance(tol)?;
52            let idenom = N::one() / (xs[i] - xs[i - j]);
53            let numer =
54                &poly_1 * &qs[i + xs.len() * (j - 1)] - &poly_2 * &qs[(i - 1) + xs.len() * (j - 1)];
55            qs[i + xs.len() * j] = numer * idenom;
56        }
57    }
58
59    for i in 0..=qs[xs.len() * xs.len() - 1].order() {
60        if qs[xs.len() * xs.len() - 1].get_coefficient(i).abs() < tol {
61            qs[xs.len() * xs.len() - 1].purge_coefficient(i);
62        }
63    }
64
65    qs[xs.len() * xs.len() - 1].purge_leading();
66    Ok(qs[xs.len() * xs.len() - 1].clone())
67}
68
69pub fn hermite<N: ComplexField>(
70    xs: &[N],
71    ys: &[N],
72    derivs: &[N],
73    tol: N::RealField,
74) -> Result<Polynomial<N>, String> {
75    if xs.len() != ys.len() {
76        return Err("hermite: slices have mismatched dimension".to_owned());
77    }
78    if xs.len() != derivs.len() {
79        return Err("hermite: derivatives have mismatched dimension".to_owned());
80    }
81
82    let mut zs = vec![N::zero(); 2 * xs.len()];
83    let mut qs = vec![N::zero(); 4 * xs.len() * xs.len()];
84
85    for i in 0..xs.len() {
86        zs[2 * i] = xs[i];
87        zs[2 * i + 1] = xs[i];
88        qs[2 * i] = ys[i];
89        qs[2 * i + 1] = ys[i];
90        qs[2 * i + 1 + (2 * xs.len())] = derivs[i];
91
92        if i != 0 {
93            qs[2 * i + (2 * xs.len())] = (qs[2 * i] - qs[2 * i - 1]) / (zs[2 * i] - zs[2 * i - 1]);
94        }
95    }
96
97    for i in 2..2 * xs.len() {
98        for j in 2..=i {
99            qs[i + j * (2 * xs.len())] = (qs[i + (j - 1) * (2 * xs.len())]
100                - qs[i - 1 + (j - 1) * (2 * xs.len())])
101                / (zs[i] - zs[i - j]);
102        }
103    }
104
105    let mut hermite = polynomial![N::zero()];
106    for i in (1..2 * xs.len()).rev() {
107        hermite += qs[i + i * (2 * xs.len())];
108        hermite *= polynomial![N::one(), -xs[(i - 1) / 2]]
109    }
110    hermite += qs[0];
111
112    for i in 0..=hermite.order() {
113        if hermite.get_coefficient(i).abs() < tol {
114            hermite.purge_coefficient(i);
115        }
116    }
117
118    hermite.purge_leading();
119    Ok(hermite)
120}
bacon_sci_1/interp/mod.rs

bacon_sci_1/interp/
mod.rs
