样条插值算法的C++ 实现

样条插值算法是一种常用的数值分析方法，用于对一组数据进行插值。下面是一个简单的C++实现：

#include <vector>

// 二次样条插值算法
double spline_interpolation(double x, const std::vector<double>& X, const std::vector<double>& Y, const std::vector<double>& M)
{
    int n = X.size() - 1;
    int i = 0;
    while (i < n && X[i+1] < x) {
        i++;
    }

    double h = X[i+1] - X[i];
    double a = (X[i+1] - x) / h;
    double b = (x - X[i]) / h;
    double y = a * Y[i] + b * Y[i+1] + ((a*a*a - a)*M[i] + (b*b*b - b)*M[i+1]) * (h*h) / 6.0;

    return y;
}

// 计算样条插值的一阶导数
std::vector<double> spline_first_derivative(const std::vector<double>& X, const std::vector<double>& Y)
{
    int n = X.size();
    std::vector<double> M(n);

    // 初始化
    M[0] = 0;
    M[n-1] = 0;

    std::vector<double> A(n-2);
    std::vector<double> B(n-2);
    std::vector<double> C(n-2);
    std::vector<double> F(n-2);

    for (int i = 0; i < n-2; i++) {
        A[i] = X[i+1] - X[i];
        C[i] = X[i+2] - X[i+1];
        B[i] = 2.0 * (A[i] + C[i]);
        F[i] = 6.0 * ((Y[i+2] - Y[i+1]) / C[i] - (Y[i+1] - Y[i]) / A[i]);
    }

    // 三对角矩阵求解
    for (int i = 1; i < n-2; i++) {
        B[i] -= A[i-1] * C[i-1] / B[i-1];
        F[i] -= A[i-1] * F[i-1] / B[i-1];
    }

    M[n-2] = F[n-3] / B[n-3];
    for (int i = n-4; i >= 0; i--) {
        M[i+1] = (F[i] - C[i] * M[i+2]) / B[i];
    }

    return M;
}

// 计算样条插值的二阶导数
std::vector<double> spline_second_derivative(const std::vector<double>& X, const std::vector<double>& Y)
{
    int n = X.size();
    std::vector<double> M(n);

    // 初始化
    M[0] = 0;
    M[n-1] = 0;

    std::vector<double> A(n-2);
    std::vector<double> B(n-2);
    std::vector<double> C(n-2);
    std::vector<double> F(n-2);

    for (int i = 0; i < n-2; i++) {
        A[i] = X[i+1] - X[i];
        C[i] = X[i+2] - X[i+1];
        B[i] = 2.0 * (A[i] + C[i]);
        F[i] = 6.0 * ((Y[i+2] - Y[i+1]) / C[i] - (Y[i+1] - Y[i]) / A[i]);
    }

    // 三对角矩阵求解
    for (int i = 1; i < n-2; i++) {
        B[i] -= A[i-1] * C[i-1] / B[i-1];
        F[i] -= A[i-1] * F[i-1] / B[i-1];
    }

    M[n-2] = F[n-3] / B[n-3];
    for (int i = n-4; i >= 0; i--) {
        M[i+1] = (F[i] - C[i] * M[i+2]) / B[i];
    }

    return M;
}

使用方法：

1. 先调用`spline_first_derivative`函数计算一阶导数`M`；
2. 再调用`spline_second_derivative`函数计算二阶导数`M`；
3. 最后在需要插值的点上调用`spline_interpolation`函数计算插值结果。

注意，这里实现的是二次样条插值算法，如果需要更高阶的插值可以参考其他资料进行实现。