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

三次样条插值是一种常用的插值方法，它可以通过给定的一组数据点，生成一条光滑的曲线，以便于对数据进行分析和预测。以下是用C++实现三次样条插值的代码：

#include <iostream>
#include <vector>

using namespace std;

// 三次样条插值函数
void cubicSplineInterpolate(vector<double>& x, vector<double>& y, vector<double>& a)
{
    int n = x.size() - 1;
    vector<double> h(n + 1), alpha(n + 1), l(n + 1), mu(n + 1), z(n + 1), c(n + 1), b(n + 1), d(n + 1);

    // 计算步长
    for (int i = 0; i < n; ++i)
        h[i] = x[i + 1] - x[i];

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

    // 递推求解 l, mu, z
    l[0] = 1.0;
    mu[0] = z[0] = 0.0;

    for (int i = 1; i < n; ++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.0;
    z[n] = c[n] = 0.0;

    // 递推求解 c, b, d
    for (int j = n - 1; 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]);
    }

    // 计算系数 a
    for (int i = 0; i < n; ++i)
    {
        a[4 * i] = y[i];
        a[4 * i + 1] = b[i];
        a[4 * i + 2] = c[i];
        a[4 * i + 3] = d[i];
    }
}

int main()
{
    vector<double> x = { 0.0, 1.0, 2.0, 3.0, 4.0 };
    vector<double> y = { 0.0, 1.0, 4.0, 9.0, 16.0 };
    vector<double> a(4 * (x.size() - 1));

    cubicSplineInterpolate(x, y, a);

    cout << '三次样条插值结果：' << endl;
    for (int i = 0; i < a.size(); i += 4)
        cout << 'f(x) = ' << a[i] << ' + ' << a[i + 1] << '*(x-' << x[i / 4] << ') + ' << a[i + 2] << '*(x-' << x[i / 4] << ')^2 + ' << a[i + 3] << '*(x-' << x[i / 4] << ')^3' << endl;

    return 0;
}

在上述代码中，我们定义了一个 cubicSplineInterpolate 函数，它接受三个向量作为输入参数。其中，x 向量表示插值点的横坐标，y 向量表示插值点的纵坐标，a 向量表示输出的插值函数的系数。在函数中，我们首先计算步长 h 和 alpha 值，然后使用递推公式求解 l, mu 和 z，最后求解 c, b 和 d。最后，我们将系数 a 存储在 a 向量中，并输出结果。当然，这里的输入数据是预先给定的，实际使用时需要根据具体情况进行修改。