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

bacon_sci/interp/
mod.rs
