graplot/
polynomial.rs

1use crate::{Plot, PlotArg, XEnd};
2
3#[derive(Debug, Clone)]
4/// Use the Polynomial struct or polynomial() function to create a polynomial function that runs through all given points.
5/// # Example
6/// ```
7/// use graplot::{x, Plot, Polynomial};
8///
9/// let poly = Polynomial::new(&[2., 3., 1.], &[2., 3., 2.]);
10/// let plot = Plot::new((poly, x(10.)));
11/// plot.show();
12/// ```
13pub struct Polynomial {
14    coefficients: Vec<f64>,
15}
16
17impl Polynomial {
18    pub fn new(xs: &[f64], ys: &[f64]) -> Polynomial {
19        polynomial(xs, ys)
20    }
21}
22
23pub fn polynomial(xs: &[f64], ys: &[f64]) -> Polynomial {
24    let degree = xs.len() - 1;
25
26    let mut coeffs = Vec::<f64>::with_capacity(xs.len() * xs.len());
27
28    for x in xs {
29        for pow in (0..=degree).rev() {
30            coeffs.push(x.powi(pow as i32));
31        }
32    }
33    Polynomial {
34        coefficients: solve_lu(xs.len(), &coeffs, ys),
35    }
36}
37
38impl PlotArg for Polynomial {
39    fn as_plot(&self) -> crate::Plot {
40        let mut xs = [0.; 200];
41
42        let mut add = -100f64;
43        for x in &mut xs {
44            *x = add / 100.;
45            add += 1.;
46        }
47
48        let mut ys = [0.; 200];
49
50        let degree = self.coefficients.len() - 1;
51
52        for (i, y) in ys.iter_mut().enumerate() {
53            for (pow, coefficient) in self.coefficients.iter().enumerate() {
54                *y += coefficient * xs[i].powi((degree - pow) as i32);
55            }
56        }
57        Plot {
58            xs: vec![xs.to_vec()],
59            ys: vec![ys.to_vec()],
60            line_desc: vec![Default::default()],
61            axis_desc: Default::default(),
62            desc: Default::default(),
63        }
64    }
65}
66
67impl PlotArg for (Polynomial, XEnd) {
68    fn as_plot(&self) -> crate::Plot {
69        let mut xs = [0.; 200];
70
71        let mut add = -100f64;
72        for x in &mut xs {
73            *x = (add / 100.) * self.1 .0;
74            add += 1.;
75        }
76
77        let mut ys = [0.; 200];
78
79        let degree = self.0.coefficients.len() - 1;
80
81        for (i, y) in ys.iter_mut().enumerate() {
82            for (pow, coefficient) in self.0.coefficients.iter().enumerate() {
83                *y += coefficient * xs[i].powi((degree - pow) as i32);
84            }
85        }
86        Plot {
87            xs: vec![xs.to_vec()],
88            ys: vec![ys.to_vec()],
89            line_desc: vec![Default::default()],
90            axis_desc: Default::default(),
91            desc: Default::default(),
92        }
93    }
94}
95
96impl PlotArg for (Polynomial, &str) {
97    fn as_plot(&self) -> Plot {
98        let mut plot = self.0.as_plot();
99        plot.line_desc = vec![self.1.into()];
100        plot
101    }
102}
103
104impl PlotArg for (Polynomial, XEnd, &str) {
105    fn as_plot(&self) -> Plot {
106        let mut plot = (self.0.clone(), self.1).as_plot();
107        plot.line_desc = vec![self.2.into()];
108        plot
109    }
110}
111
112pub fn solve_lu(n: usize, lhs: &[f64], rhs: &[f64]) -> Vec<f64> {
113    let mut lu = vec![0f64; n * n];
114    let mut sum;
115    for i in 0..n {
116        for j in i..n {
117            sum = 0.;
118            for k in 0..i {
119                sum += lu[i * n + k] * lu[k * n + j];
120            }
121            lu[i * n + j] = lhs[i * n + j] - sum;
122        }
123        for j in (i + 1)..n {
124            sum = 0.;
125            for k in 0..i {
126                sum += lu[j * n + k] * lu[k * n + i];
127            }
128            lu[j * n + i] = (1. / lu[i * n + i]) * (lhs[j * n + i] - sum)
129        }
130    }
131
132    let mut y = vec![0.; n];
133    for i in 0..n {
134        sum = 0.;
135        for k in 0..i {
136            sum += lu[i * n + k] * y[k];
137        }
138        y[i] = rhs[i] - sum;
139    }
140
141    let mut x = vec![0.; n];
142    for i in (0..n).rev() {
143        sum = 0.;
144        for k in (i + 1)..n {
145            sum += lu[i * n + k] * x[k];
146        }
147        x[i] = (1. / lu[i * n + i]) * (y[i] - sum);
148    }
149    x
150}
graplot/polynomial.rs

graplot/
polynomial.rs
