class CubicSplineInterpolation {
private:
std::vector x; // 存储插值点横坐标
std::vector y; // 存储插值点纵坐标
std::vector h; // 存储每个小区间的步长
std::vector a; // 存储三次样条函数的系数a
std::vector b; // 存储三次样条函数的系数b
std::vector c; // 存储三次样条函数的系数c
std::vector d; // 存储三次样条函数的系数d
public:
// 构造函数，传入插值点横坐标和纵坐标
CubicSplineInterpolation(const std::vector& _x, const std::vector& _y) {
if (_x.size() != _y.size())
throw std::invalid_argument("The sizes of x and y must be the same.");
x = _x;
y = _y;
int n = x.size();
h.resize(n - 1);
a.resize(n);
b.resize(n);
c.resize(n - 1);
d.resize(n - 1);
for (int i = 0; i < n - 1; i++)
h[i] = x[i + 1] - x[i];
// 利用Thomas算法求解三次样条函数的系数
std::vector alpha(n - 1);
std::vector beta(n - 1);
alpha[0] = beta[0] = 0;
for (int i = 1; i < n - 1; i++) {
alpha[i] = h[i - 1] / (h[i - 1] + h[i]);
beta[i] = (3 * (y[i + 1] - y[i]) / h[i] - 3 * (y[i] - y[i - 1]) / h[i - 1]) / (h[i - 1] + h[i]);
}
alpha[n - 1] = beta[n - 1] = 0;
std::vector c_prime(n);
std::vector l(n);
std::vector mu(n);
std::vector z(n);
l[0] = 1;
mu[0] = z[0] = 0;
for (int i = 1; i < n - 1; i++) {
l[i] = 2 - alpha[i] * mu[i - 1];
mu[i] = alpha[i] / l[i];
z[i] = (beta[i] - alpha[i] * z[i - 1]) / l[i];
}
l[n - 1] = 1;
z[n - 1] = c_prime[n - 1] = 0;
for (int i = n - 2; i >= 0; i--) {
c_prime[i] = z[i] - mu[i] * c_prime[i + 1];
c[i] = (3 * (y[i + 1] - y[i]) / h[i] - 2 * c_prime[i] - c_prime[i + 1]) / h[i];
d[i] = (c_prime[i + 1] + c_prime[i] - 2 * (y[i + 1] - y[i]) / (h[i] * h[i])) / (h[i] * h[i]);
a[i] = y[i];
b[i] = (y[i + 1] - y[i]) / h[i] - h[i] * (2 * c_prime[i] + c_prime[i + 1]) / 3;
}
a[n - 1] = y[n - 1];
b[n - 1] = c[n - 1] = d[n - 1] = 0;
}
// 根据横坐标x计算插值函数的纵坐标
double operator()(const double x_val) {
int n = x.size();
int i = 0;
while (i < n - 1 && x[i + 1] <= x_val)
i++;
double t = (x_val - x[i]) / h[i];
double y_val = a[i] + b[i] * t + c[i] * std::pow(t, 2) + d[i] * std::pow(t, 3);
return y_val;
}
}
