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
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
use operation::extra_ops::PowOps;
#[allow(unused_imports)]
use structure::matrix::*;
#[allow(unused_imports)]
use structure::polynomial::*;
#[allow(unused_imports)]
use structure::vector::*;
#[allow(unused_imports)]
use util::non_macro::*;
pub fn cubic_spline(node_x: Vector, node_y: Vector) -> Vec<Polynomial> {
let n = node_x.len() - 1;
assert_eq!(n, node_y.len() - 1);
let mut h = vec![0f64; n];
let mut b = vec![0f64; n];
let mut v = vec![0f64; n];
let mut u = vec![0f64; n];
for i in 0..n {
if i == 0 {
h[i] = node_x[i + 1] - node_x[i];
b[i] = (node_y[i + 1] - node_y[i]) / h[i];
} else {
h[i] = node_x[i + 1] - node_x[i];
b[i] = (node_y[i + 1] - node_y[i]) / h[i];
v[i] = 2. * (h[i] + h[i - 1]);
u[i] = 6. * (b[i] - b[i - 1]);
}
}
let mut m = matrix(vec![0f64; (n - 1) * (n - 1)], n - 1, n - 1, Col);
for i in 0..n - 2 {
m[(i, i)] = v[i + 1];
m[(i + 1, i)] = h[i + 1];
m[(i, i + 1)] = h[i + 1];
}
m[(n - 2, n - 2)] = v[n - 1];
let z_inner = m.inv().unwrap() * Vec::from(&u[1..]).to_matrix();
let mut z = vec![0f64];
z.extend(&z_inner.data);
z.push(0f64);
let mut s: Vec<Polynomial> = Vec::new();
for i in 0..n {
let t_i = node_x[i];
let t_i1 = node_x[i + 1];
let z_i = z[i];
let z_i1 = z[i + 1];
let h_i = h[i];
let y_i = node_y[i];
let y_i1 = node_y[i + 1];
let temp1 = poly(vec![1f64, -t_i]);
let temp2 = poly(vec![1f64, -t_i1]);
let term1 = temp1.powi(3) * (z_i1 / (6f64 * h_i));
let term2 = temp2.powi(3) * (-z_i / (6f64 * h_i));
let term3 = temp1 * (y_i1 / h_i - z_i1 * h_i / 6.);
let term4 = temp2 * (-y_i / h_i + h_i * z_i / 6.0);
s.push(term1 + term2 + term3 + term4);
}
return s;
}