C++ 实现三次样条插值算法

三次样条插值是一种通过已知的数据点构建出一条光滑曲线的方法。在 C++ 中，可以使用以下步骤实现三次样条插值：

1. 构建三对角矩阵：

首先，需要根据已知数据点构建出一个三对角矩阵，其中对角线元素为 2，次对角线元素为 1，边界处的元素需要特殊处理。具体实现如下：

void build_tridiagonal_matrix(vector<double>& h, vector<double>& a, vector<double>& b, vector<double>& c, vector<double>& d) {
    int n = h.size() - 1;
    a.resize(n + 1, 0);
    b.resize(n + 1, 0);
    c.resize(n + 1, 0);
    d.resize(n + 1, 0);
    for (int i = 1; i < n; i++) {
        a[i] = h[i];
        b[i] = 2 * (h[i] + h[i + 1]);
        c[i] = h[i + 1];
        d[i] = 3 * ((d[i + 1] - d[i]) / h[i + 1] - (d[i] - d[i - 1]) / h[i]);
    }
    // 边界处的元素需要特殊处理
    b[0] = 1;
    c[0] = 0;
    d[0] = 0;
    a[n] = 0;
    b[n] = 1;
    d[n] = 0;
}

其中，h 为已知数据点之间的距离，a、b、c、d 分别为三对角矩阵的对角线元素、次对角线元素和常数项。

2. 解三对角线性方程组：

使用追赶法或高斯消元法求解三对角线性方程组，得到每个数据点处的二阶导数。具体实现如下：

void solve_tridiagonal_matrix(vector<double>& a, vector<double>& b, vector<double>& c, vector<double>& d, vector<double>& m) {
    int n = a.size() - 1;
    for (int i = 1; i < n; i++) {
        double w = a[i] / b[i - 1];
        b[i] = b[i] - w * c[i - 1];
        d[i] = d[i] - w * d[i - 1];
    }
    m[n] = d[n] / b[n];
    for (int i = n - 1; i >= 0; i--) {
        m[i] = (d[i] - c[i] * m[i + 1]) / b[i];
    }
}

其中，m 为每个数据点处的二阶导数。

3. 计算插值函数：

根据已知数据点和每个数据点处的二阶导数，计算插值函数。具体实现如下：

double cubic_spline_interpolation(double x, vector<double>& x_data, vector<double>& y_data, vector<double>& m) {
    int n = x_data.size() - 1;
    int k = 0;
    while (k < n && x > x_data[k + 1]) {
        k++;
    }
    double h = x_data[k + 1] - x_data[k];
    double a = (m[k + 1] - m[k]) / (6 * h);
    double b = m[k] / 2;
    double c = (y_data[k + 1] - y_data[k]) / h - h * (2 * m[k] + m[k + 1]) / 6;
    double d = y_data[k];
    return a * pow(x - x_data[k], 3) + b * pow(x - x_data[k], 2) + c * (x - x_data[k]) + d;
}

其中，x 为待插值的自变量，x_data 和 y_data 为已知数据点的横纵坐标，m 为每个数据点处的二阶导数。

完整代码如下：