c++ 三次样条插值 自然边界条件

三次样条插值中的自然边界条件是指在插值区间的两端点处，二阶导数值为0。这样做的目的是为了使得插值函数的两端点更加平滑，避免出现不必要的振荡。

下面是使用自然边界条件进行三次样条插值的 C++ 代码实现：

#include <iostream>
#include <vector>

using namespace std;

void cubic_spline_interpolation(vector<double>& x, vector<double>& y, vector<double>& M, double x_new, double& y_new) {
    int n = x.size();
    int i = 0;
    while (i < n && x[i] < x_new) {
        i++;
    }
    double h = x[i] - x[i-1];
    double a = (M[i] - M[i-1]) / (6 * h);
    double b = M[i-1] / 2;
    double c = (y[i] - y[i-1]) / h - (2 * h * M[i-1] + h * M[i]) / 6;
    double d = y[i-1];
    y_new = a * pow(x_new - x[i-1], 3) + b * pow(x_new - x[i-1], 2) + c * (x_new - x[i-1]) + d;
}

int main() {
    // example data
    vector<double> x = {0, 1, 2, 3, 4};
    vector<double> y = {0, 1, 4, 9, 16};
    
    int n = x.size();
    vector<double> h(n-1);
    for (int i = 0; i < n-1; i++) {
        h[i] = x[i+1] - x[i];
    }
    
    // solve tridiagonal linear system to get M
    vector<double> a(n-2);
    vector<double> b(n-2);
    vector<double> c(n-2);
    vector<double> d(n-2);
    for (int i = 0; i < n-2; i++) {
        a[i] = h[i];
        b[i] = 2 * (h[i] + h[i+1]);
        c[i] = h[i+1];
        d[i] = 6 * ((y[i+2] - y[i+1]) / h[i+1] - (y[i+1] - y[i]) / h[i]);
    }
    vector<double> M(n);
    M[0] = 0;
    M[n-1] = 0;
    for (int i = 0; i < n-2; i++) {
        double ratio = a[i] / b[i];
        b[i+1] -= ratio * c[i];
        d[i+1] -= ratio * d[i];
    }
    M[n-2] = d[n-3] / b[n-3];
    for (int i = n-4; i >= 0; i--) {
        M[i+1] = (d[i] - c[i] * M[i+2]) / b[i];
    }
    
    // interpolate at x_new = 2.5
    double x_new = 2.5;
    double y_new;
    cubic_spline_interpolation(x, y, M, x_new, y_new);
    cout << "Interpolated value at x = " << x_new << " is y = " << y_new << endl;
    
    return 0;
}

其中，解三对角线性方程组的部分可以使用 Thomas 算法来进行优化，这里为了方便直接使用了高斯消元法。