sciforge_lib/maths/interpolation/
lagrange.rs1pub fn lagrange_interpolate(points: &[(f64, f64)], x: f64) -> f64 {
2 let n = points.len();
3 let mut result = 0.0;
4 for i in 0..n {
5 let mut basis = 1.0;
6 for j in 0..n {
7 if i != j {
8 basis *= (x - points[j].0) / (points[i].0 - points[j].0);
9 }
10 }
11 result += basis * points[i].1;
12 }
13 result
14}
15
16pub fn lagrange_barycentric(points: &[(f64, f64)], x: f64) -> f64 {
17 let n = points.len();
18 let mut weights = vec![1.0; n];
19 for i in 0..n {
20 for j in 0..n {
21 if i != j {
22 weights[i] /= points[i].0 - points[j].0;
23 }
24 }
25 }
26
27 let mut numer = 0.0;
28 let mut denom = 0.0;
29 for i in 0..n {
30 let diff = x - points[i].0;
31 if diff.abs() < 1e-30 {
32 return points[i].1;
33 }
34 let t = weights[i] / diff;
35 numer += t * points[i].1;
36 denom += t;
37 }
38 numer / denom
39}
40
41pub fn neville(points: &[(f64, f64)], x: f64) -> f64 {
42 let n = points.len();
43 let mut table: Vec<f64> = points.iter().map(|p| p.1).collect();
44 for k in 1..n {
45 for i in 0..n - k {
46 table[i] = ((x - points[i + k].0) * table[i] - (x - points[i].0) * table[i + 1])
47 / (points[i].0 - points[i + k].0);
48 }
49 }
50 table[0]
51}
52
53pub fn divided_differences(points: &[(f64, f64)]) -> Vec<f64> {
54 let n = points.len();
55 let mut dd: Vec<f64> = points.iter().map(|p| p.1).collect();
56 for k in 1..n {
57 for i in (k..n).rev() {
58 dd[i] = (dd[i] - dd[i - 1]) / (points[i].0 - points[i - k].0);
59 }
60 }
61 dd
62}
63
64pub fn newton_interpolate(points: &[(f64, f64)], x: f64) -> f64 {
65 let dd = divided_differences(points);
66 let mut result = dd[points.len() - 1];
67 for i in (0..points.len() - 1).rev() {
68 result = result * (x - points[i].0) + dd[i];
69 }
70 result
71}
72
73pub fn lagrange_derivative(points: &[(f64, f64)], x: f64) -> f64 {
74 let n = points.len();
75 let mut result = 0.0;
76 for i in 0..n {
77 let mut sum = 0.0;
78 for j in 0..n {
79 if j == i {
80 continue;
81 }
82 let mut prod = 1.0;
83 for k in 0..n {
84 if k == i || k == j {
85 continue;
86 }
87 prod *= (x - points[k].0) / (points[i].0 - points[k].0);
88 }
89 sum += prod / (points[i].0 - points[j].0);
90 }
91 result += points[i].1 * sum;
92 }
93 result
94}
95
96pub fn chebyshev_nodes(n: usize, a: f64, b: f64) -> Vec<f64> {
97 let pi = std::f64::consts::PI;
98 (0..n)
99 .map(|k| {
100 let t = ((2 * k + 1) as f64 * pi / (2.0 * n as f64)).cos();
101 0.5 * (a + b) + 0.5 * (b - a) * t
102 })
103 .collect()
104}
105
106pub fn chebyshev_interpolate(f: impl Fn(f64) -> f64, a: f64, b: f64, n: usize, x: f64) -> f64 {
107 let nodes = chebyshev_nodes(n, a, b);
108 let points: Vec<(f64, f64)> = nodes.iter().map(|&xi| (xi, f(xi))).collect();
109 lagrange_barycentric(&points, x)
110}
111
112pub fn rational_interpolate(points: &[(f64, f64)], x: f64) -> f64 {
113 let n = points.len();
114 let mut c: Vec<f64> = points.iter().map(|p| p.1).collect();
115 let mut d: Vec<f64> = c.clone();
116 let mut ns = 0;
117 let mut dif = (x - points[0].0).abs();
118 for (i, pi) in points.iter().enumerate().skip(1) {
119 let dift = (x - pi.0).abs();
120 if dift < dif {
121 ns = i;
122 dif = dift;
123 }
124 }
125 let mut y = points[ns].1;
126 ns -= 0;
127 for m in 1..n {
128 for i in 0..n - m {
129 let w = c[i + 1] - d[i];
130 let h = points[i + m].0 - x;
131 let t = (points[i].0 - x) * d[i] / h;
132 let dd = t - c[i + 1];
133 if dd.abs() < 1e-30 {
134 return y;
135 }
136 let dd = w / dd;
137 d[i] = c[i + 1] * dd;
138 c[i] = t * dd;
139 }
140 if 2 * (ns + 1) < n - m {
141 y += c[ns + 1];
142 } else {
143 y += d[ns];
144 ns = ns.saturating_sub(1);
145 }
146 }
147 y
148}
149
150pub fn newton_forward_difference(points: &[(f64, f64)], x: f64) -> f64 {
151 let n = points.len();
152 let h = if n > 1 {
153 points[1].0 - points[0].0
154 } else {
155 1.0
156 };
157 let s = (x - points[0].0) / h;
158 let mut table: Vec<Vec<f64>> = vec![points.iter().map(|p| p.1).collect()];
159 for k in 1..n {
160 let prev = &table[k - 1];
161 let diff: Vec<f64> = (0..prev.len() - 1).map(|i| prev[i + 1] - prev[i]).collect();
162 table.push(diff);
163 }
164 let mut result = table[0][0];
165 let mut s_prod = 1.0;
166 for (k, tk) in table.iter().enumerate().skip(1) {
167 s_prod *= (s - (k - 1) as f64) / k as f64;
168 if tk.is_empty() {
169 break;
170 }
171 result += s_prod * tk[0];
172 }
173 result
174}
175
176pub fn interpolation_error_bound(points: &[(f64, f64)], x: f64, max_deriv: f64) -> f64 {
177 let n = points.len();
178 let mut prod = 1.0;
179 for pi in points {
180 prod *= (x - pi.0).abs();
181 }
182 let mut factorial = 1.0;
183 for i in 1..=n {
184 factorial *= i as f64;
185 }
186 max_deriv * prod / factorial
187}
188
189pub fn lebesgue_constant(nodes: &[f64], eval_points: &[f64]) -> f64 {
190 let n = nodes.len();
191 let mut max_lambda = 0.0_f64;
192 for &x in eval_points {
193 let mut lambda = 0.0;
194 for i in 0..n {
195 let mut li = 1.0;
196 for j in 0..n {
197 if i != j {
198 li *= (x - nodes[j]) / (nodes[i] - nodes[j]);
199 }
200 }
201 lambda += li.abs();
202 }
203 max_lambda = max_lambda.max(lambda);
204 }
205 max_lambda
206}