C++ 实现三次样条插值：曲线拟合算法

三次样条插值是一种常用的曲线拟合方法，其基本思想是将给定的离散数据点通过三次函数拟合，使得曲线在每个数据点处一阶、二阶导数连续，从而获得更加平滑的曲线。下面给出 C++ 实现的代码。

假设给定的数据点为(x[i], y[i])，其中 i=0,1,2,…,n-1，要求求出拟合曲线在区间[x[i], x[i+1]]上的系数 a[i], b[i], c[i], d[i]，使得在该区间内的拟合函数为：

S(x) = a[i] + bi + ci^2 + di^3, x∈[x[i], x[i+1]]

其中，a[i], b[i], c[i], d[i] 为待求系数。

求解三次样条插值的系数

首先，需要求解出每个区间内的系数。对于第 i 个区间[x[i], x[i+1]]，设 S(x) 在该区间内的一阶、二阶导数分别为 S'(x) 和 S''(x)，并令 h[i]=x[i+1]-x[i]，则有：

S(x[i]) = y[i]

S(x[i+1]) = y[i+1]

S'(x[i]) = S'(x[i+1])

S''(x[i]) = S''(x[i+1])

S''(x) = 2c[i] + 6di, x∈[x[i], x[i+1]]

由于 S(x) 是三次函数，因此可以表示为：

S(x) = a[i] + bi + ci^2 + di^3

代入上述四个条件，得到如下方程组：

a[i] = y[i]

b[i] = S'(x[i])

c[i] = (3(y[i+1]-y[i])/h[i]^2 - 2S'(x[i])/h[i] - S'(x[i+1])/h[i]) / h[i]

d[i] = (2(y[i]-y[i+1])/h[i]^3 + (S'(x[i])+S'(x[i+1]))/h[i]^2 - 2c[i]/h[i]) / h[i]

对于 i=1,2,…,n-2，还需要添加两个附加条件，使得插值函数在首尾两个数据点处一阶导数值等于给定数据点的一阶导数值。具体地，设 f'(x[0]) 和 f'(x[n-1]) 分别表示给定数据点的一阶导数值，然后添加以下两个条件：

S'(x[0]) = f'(x[0])

S'(x[n-1]) = f'(x[n-1])

这样就可以得到一个 (n+1)×(n+1) 的线性方程组，其中未知数为 a[0], b[0], c[0], d[0], a[1], b[1], c[1], d[1], …, a[n-1], b[n-1], c[n-1], d[n-1]。用高斯消元法或 LU 分解法求解该方程组即可得到系数。

使用三次样条插值函数求解

得到系数之后，就可以使用插值函数求解了。对于任意给定的 x，首先找到其所在的区间[i, i+1]，然后代入上述的插值函数，即可得到插值结果。

以下是完整的 C++ 代码实现：

#include <iostream>
#include <vector>

using namespace std;

// 三次样条插值类
class CubicSplineInterpolator {
public:
    // 构造函数，输入数据点
    CubicSplineInterpolator(const vector<double>& x, const vector<double>& y) : x(x), y(y) {
        n = x.size();
        // 初始化系数矩阵
        A.resize(n + 1, vector<double>(n + 1, 0.0));
        b.resize(n + 1, 0.0);
        // 计算系数
        calculateCoefficients();
    }

    // 插值函数，输入 x 值，返回插值结果
    double interpolate(double x) {
        // 找到 x 所在的区间
        int i = findInterval(x);
        // 计算插值结果
        return a[i] + b[i] * (x - x[i]) + c[i] * pow(x - x[i], 2) + d[i] * pow(x - x[i], 3);
    }

private:
    // 计算系数
    void calculateCoefficients() {
        // 计算 h[i]
        h.resize(n - 1); 
        for (int i = 0; i < n - 1; ++i) {
            h[i] = x[i + 1] - x[i];
        }
        // 构建系数矩阵 A
        A[0][0] = 1.0;
        A[n][n] = 1.0;
        for (int i = 1; i < n; ++i) {
            A[i][i - 1] = 1.0 / h[i - 1];
            A[i][i] = -2.0 / h[i - 1] - 2.0 / h[i];
            A[i][i + 1] = 1.0 / h[i];
        }
        // 构建向量 b
        b[0] = 0.0; // 首尾两点的一阶导数设为 0
        b[n] = 0.0;
        for (int i = 1; i < n; ++i) {
            b[i] = 3.0 * (y[i + 1] - y[i]) / h[i] - 3.0 * (y[i] - y[i - 1]) / h[i - 1];
        }
        // 使用高斯消元法求解线性方程组
        solveLinearEquation(A, b);
        // 计算系数
        a.resize(n); 
        b.resize(n); 
        c.resize(n); 
        d.resize(n); 
        for (int i = 0; i < n - 1; ++i) {
            a[i] = y[i];
            b[i] = b[i];
            c[i] = (3.0 * (y[i + 1] - y[i]) / pow(h[i], 2) - 2.0 * b[i] / h[i] - b[i + 1] / h[i]) / h[i];
            d[i] = (2.0 * (y[i] - y[i + 1]) / pow(h[i], 3) + (b[i] + b[i + 1]) / pow(h[i], 2) - 2.0 * c[i] / h[i]) / h[i];
        }
    }

    // 找到 x 所在的区间
    int findInterval(double x) {
        for (int i = 0; i < n - 1; ++i) {
            if (x >= x[i] && x <= x[i + 1]) {
                return i;
            }
        }
        return -1; // 未找到区间
    }

    // 使用高斯消元法求解线性方程组
    void solveLinearEquation(vector<vector<double>>& A, vector<double>& b) {
        int n = A.size();
        // 消元过程
        for (int i = 0; i < n - 1; ++i) {
            // 将主元归一
            double pivot = A[i][i];
            for (int j = i; j < n; ++j) {
                A[i][j] /= pivot;
            }
            b[i] /= pivot;
            // 消去其他行中的主元
            for (int k = i + 1; k < n; ++k) {
                double factor = A[k][i];
                for (int j = i; j < n; ++j) {
                    A[k][j] -= factor * A[i][j];
                }
                b[k] -= factor * b[i];
            }
        }
        // 回代过程
        for (int i = n - 1; i >= 0; --i) {
            // 计算当前未知数的值
            b[i] -= accumulate(b.begin() + i + 1, b.end(), 0.0);
            // 消去该未知数的影响
            for (int j = i - 1; j >= 0; --j) {
                b[j] -= A[j][i] * b[i];
            }
        }
    }

private:
    vector<double> x; // x 坐标
    vector<double> y; // y 坐标
    int n; // 数据点数
    vector<vector<double>> A; // 系数矩阵
    vector<double> b; // 右端向量
    vector<double> h; // 每个区间的长度
    vector<double> a; // 三次样条插值系数 a
    vector<double> b; // 三次样条插值系数 b
    vector<double> c; // 三次样条插值系数 c
    vector<double> d; // 三次样条插值系数 d
};

int main() {
    // 输入数据点
    vector<double> x = {0.0, 1.0, 2.0, 3.0};
    vector<double> y = {1.0, 2.0, 3.0, 4.0};
    // 创建三次样条插值对象
    CubicSplineInterpolator interpolator(x, y);
    // 插值结果
    double interpolatedValue = interpolator.interpolate(1.5);
    cout << "插值结果: " << interpolatedValue << endl;
    return 0;
}

该代码中，首先定义了一个 CubicSplineInterpolator 类，用于实现三次样条插值算法。该类包含了构造函数、插值函数以及一些私有函数，用于计算系数、寻找区间以及求解线性方程组。

在 main 函数中，首先创建了一个 CubicSplineInterpolator 对象，并传入数据点。然后使用 interpolate 函数对指定的 x 值进行插值，并输出插值结果。

该代码可以作为参考，您可以根据自己的需要进行修改和扩展。