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

三次样条插值是一种常用的数据插值方法，可以通过一些已知数据点，生成一个连续的三次函数曲线，该曲线在每个数据点处的函数值与实际数据点的值相等，且在相邻数据点之间具有一定的平滑性。

以下是C++代码实现三次样条插值的函数：

#include <iostream>
#include <vector>

using namespace std;

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

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

    // 计算alpha
    for (int i = 2; 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[1] = 1.0; mu[1] = z[1] = 0.0;
    for (int i = 2; 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];
    }

    // 计算c, b, d
    l[n] = 1.0; z[n] = c[n] = d[n] = 0.0;
    for (int i = n - 1; i >= 1; i--)
    {
        c[i] = z[i] - mu[i] * c[i + 1];
        b[i] = (y[i + 1] - y[i]) / h[i] - h[i] * (c[i + 1] + 2.0 * c[i]) / 3.0;
        d[i] = (c[i + 1] - c[i]) / (3.0 * h[i]);
        a[i] = y[i];
    }
}

int main()
{
    vector<double> x{ 1, 2, 3, 4, 5 };
    vector<double> y{ 2.0, 3.1, 3.9, 4.2, 4.1 };
    vector<double> a(x.size()), b(x.size()), c(x.size()), d(x.size());

    // 计算样条插值
    spline(x, y, a, b, c, d);

    // 输出结果
    for (int i = 1; i < x.size(); i++)
    {
        cout << "S" << i << "(x) = ";
        if (c[i] != 0)
            cout << c[i] << "(x - " << x[i] << ")^3 + ";
        if (b[i] != 0)
            cout << b[i] << "(x - " << x[i] << ")^2 + ";
        if (a[i] != 0)
            cout << a[i] << "(x - " << x[i] << ") + ";
        cout << d[i] << endl;
    }

    return 0;
}

在实现中，我们使用了向量来存储已知数据点的x坐标和y坐标，以及计算出的a、b、c、d系数。spline函数的输入参数是x和y向量，输出参数是a、b、c、d向量。

spline函数的实现过程如下：

计算向量h，其中h[i] = x[i] - x[i - 1]
计算向量alpha，其中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[1] = 1.0，mu[1] = z[1] = 0.0，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]，其中i从2到n-1
计算向量c、b、d、a，其中l[n] = 1.0，z[n] = c[n] = d[n] = 0.0，c[i] = z[i] - mu[i] * c[i + 1]，b[i] = (y[i + 1] - y[i]) / h[i] - h[i] * (c[i + 1] + 2.0 * c[i]) / 3.0，d[i] = (c[i + 1] - c[i]) / (3.0 * h[i])，a[i] = y[i]，其中i从n-1到1

最后，我们可以将得到的a、b、c、d系数代入样条插值函数的公式，得到每个数据点之间的三次函数表达式。在上面的代码中，我们输出了每个数据点之间的函数表达式。