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

三次样条插值是一种常用于数据拟合的方法，它可以通过已知的数据点，拟合出一条光滑的曲线来，以便于对未知数据的估计和预测。下面是用C++实现三次样条插值的代码：

#include <iostream>
#include <vector>
#include <cmath>

using namespace std;

vector<double> spline(vector<double> x, vector<double> y)
{
    int n = x.size();
    vector<double> h(n - 1), alpha(n - 1), l(n), mu(n - 1), z(n), c(n), b(n - 1), d(n - 1);

    for (int i = 0; i < n - 1; i++) {
        h[i] = x[i + 1] - x[i];
    }

    for (int i = 1; i < n - 1; i++) {
        alpha[i] = 3.0 / h[i] * (y[i + 1] - y[i]) - 3.0 / h[i - 1] * (y[i] - y[i - 1]);
    }

    l[0] = 1;
    mu[0] = 0;
    z[0] = 0;

    for (int i = 1; i < n - 1; i++) {
        l[i] = 2.0 * (x[i + 1] - x[i - 1]) - h[i - 1] * mu[i - 1];
        mu[i] = h[i] / l[i];
        z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / l[i];
    }

    l[n - 1] = 1;
    z[n - 1] = 0;
    c[n - 1] = 0;

    for (int j = n - 2; j >= 0; j--) {
        c[j] = z[j] - mu[j] * c[j + 1];
        b[j] = (y[j + 1] - y[j]) / h[j] - h[j] * (c[j + 1] + 2.0 * c[j]) / 3.0;
        d[j] = (c[j + 1] - c[j]) / (3.0 * h[j]);
    }

    vector<double> coeffs(n * 4);

    for (int i = 0; i < n - 1; i++) {
        coeffs[i * 4] = y[i];
        coeffs[i * 4 + 1] = b[i];
        coeffs[i * 4 + 2] = c[i];
        coeffs[i * 4 + 3] = d[i];
    }

    coeffs[(n - 1) * 4] = y[n - 1];

    return coeffs;
}

double evaluate(vector<double> coeffs, double x)
{
    int n = coeffs.size() / 4;
    int i = 0;

    while (i < n - 1 && x > coeffs[i * 4]) {
        i++;
    }

    i--;

    double dx = x - coeffs[i * 4];
    double result = coeffs[i * 4] + dx * (coeffs[i * 4 + 1] + dx * (coeffs[i * 4 + 2] + dx * coeffs[i * 4 + 3]));

    return result;
}

int main()
{
    vector<double> x = { 0, 1, 2, 3, 4 };
    vector<double> y = { 0, 1, 0, -1, 0 };
    vector<double> coeffs = spline(x, y);

    for (double i = 0; i <= 4; i += 0.1) {
        cout << 'x = ' << i << ', y = ' << evaluate(coeffs, i) << endl;
    }

    return 0;
}

该代码中，spline函数用于计算三次样条插值的系数，evaluate函数用于根据插值系数和输入的x值，计算对应的y值。在主函数中，我们给出了一些简单的测试数据，然后用evaluate函数计算插值曲线对应的y值，并将其输出到控制台中。

注：该代码中的插值曲线是分段的，每一段都是一个三次多项式。如果需要得到一条光滑的整体曲线，可以将相邻的多项式进行拼接。