C++ 实现三次样条插值 - 数据拟合与函数逼近

三次样条插值是一种广泛应用于数据拟合和函数逼近的方法，它可以在数据点之间插入一条平滑的曲线，从而更好地描述数据的特征。

以下是用 C++ 实现三次样条插值的示例代码：

#include <iostream>
#include <vector>

using namespace std;

vector<double> calculateCoefficients(vector<double> x, vector<double> y) {
    int n = x.size();

    // Step 1: Calculate delta-x and delta-y arrays
    vector<double> dx(n - 1), dy(n - 1);
    for (int i = 0; i < n - 1; i++) {
        dx[i] = x[i + 1] - x[i];
        dy[i] = y[i + 1] - y[i];
    }

    // Step 2: Calculate the tridiagonal matrix and right-hand side vector
    vector<double> A(n);
    vector<double> B(n);
    vector<double> C(n);
    vector<double> D(n);
    for (int i = 1; i < n - 1; i++) {
        A[i] = dx[i - 1];
        B[i] = 2 * (dx[i - 1] + dx[i]);
        C[i] = dx[i];
        D[i] = 3 * (dy[i] / dx[i] - dy[i - 1] / dx[i - 1]);
    }

    // Step 3: Apply boundary conditions
    A[0] = 0;
    B[0] = 1;
    C[0] = 0;
    D[0] = 0;

    A[n - 1] = 0;
    B[n - 1] = 1;
    C[n - 1] = 0;
    D[n - 1] = 0;

    // Step 4: Solve the tridiagonal system using Thomas Algorithm
    for (int i = 1; i < n; i++) {
        double m = A[i] / B[i - 1];
        B[i] = B[i] - m * C[i - 1];
        D[i] = D[i] - m * D[i - 1];
    }

    vector<double> M(n);
    M[n - 1] = D[n - 1] / B[n - 1];

    for (int i = n - 2; i >= 0; i--) {
        M[i] = (D[i] - C[i] * M[i + 1]) / B[i];
    }

    // Step 5: Calculate the cubic polynomial coefficients
    vector<double> a(n - 1), b(n - 1), c(n - 1), d(n - 1);
    for (int i = 0; i < n - 1; i++) {
        a[i] = y[i];
        b[i] = dy[i] / dx[i] - dx[i] * M[i] / 3 - dx[i] * M[i + 1] / 6;
        c[i] = M[i] / 2;
        d[i] = (M[i + 1] - M[i]) / (6 * dx[i]);
    }

    vector<double> coefficients;
    coefficients.insert(coefficients.end(), a.begin(), a.end());
    coefficients.insert(coefficients.end(), b.begin(), b.end());
    coefficients.insert(coefficients.end(), c.begin(), c.end());
    coefficients.insert(coefficients.end(), d.begin(), d.end());

    return coefficients;
}

double evaluateSpline(double x, vector<double> coefficients, vector<double> xData) {
    int n = xData.size();
    int i = 0;

    // Find the interval in which x resides
    while (i < n && x > xData[i]) {
        i++;
    }

    i--;

    // Evaluate the cubic polynomial in that interval
    double dx = x - xData[i];
    double a = coefficients[i];
    double b = coefficients[n + i];
    double c = coefficients[2 * n + i];
    double d = coefficients[3 * n + i];

    return a + b * dx + c * dx * dx + d * dx * dx * dx;
}

int main() {
    vector<double> xData = {0, 1, 2, 3}; // x-coordinates of data points
    vector<double> yData = {0, 1, 4, 9}; // y-coordinates of data points

    vector<double> coefficients = calculateCoefficients(xData, yData);

    // Evaluate the spline at some test points
    double xTest1 = 0.5;
    double yTest1 = evaluateSpline(xTest1, coefficients, xData);
    cout << "Spline value at " << xTest1 << " is " << yTest1 << endl;

    double xTest2 = 2.5;
    double yTest2 = evaluateSpline(xTest2, coefficients, xData);
    cout << "Spline value at " << xTest2 << " is " << yTest2 << endl;

    return 0;
}

该代码首先定义了一个calculateCoefficients函数，该函数接受两个向量x和y，表示数据点的x坐标和y坐标。它首先计算出dx和dy两个向量，表示相邻数据点之间的差值。然后，它构造了一个三对角矩阵，并计算了右手边向量。接下来，应用边界条件，然后使用Thomas算法解决三对角系统。最后，它计算了三次多项式的系数。

接下来，它定义了一个evaluateSpline函数，该函数接受一个x值、系数向量和数据点的x向量。该函数首先找到x所在的区间，然后在该区间内对三次多项式进行求值。

最后，该代码定义了一个main函数，该函数创建了一些测试点，并使用calculateCoefficients函数计算三次样条插值系数。然后，它在测试点上评估样条，并输出结果。

注意，此代码仅用于说明目的，可能需要进行更多的错误检查和边界条件处理。