用c++实现三次样条插值

三次样条插值是一种常用的数值计算方法，可以用于在给定的一组数据点上构造一条平滑的曲线。以下是用C++实现三次样条插值的代码：

#include <iostream>
#include <vector>

using namespace std;

struct Spline {
    vector<double> a, b, c, d, x;
    int n;

    void build(vector<double> &y) {
        n = y.size();
        a = y;
        b.resize(n);
        c.resize(n);
        d.resize(n);
        x.resize(n);

        // step 1: compute x
        for (int i = 0; i < n; i++) {
            x[i] = i;
        }

        // step 2: solve tridiagonal system
        vector<double> h(n), alpha(n), l(n), mu(n), z(n);
        for (int i = 1; i < n; i++) {
            h[i] = x[i] - x[i - 1];
            alpha[i] = 3.0 / h[i] * (a[i + 1] - a[i]) - 3.0 / h[i - 1] * (a[i] - a[i - 1]);
        }
        l[0] = 1;
        mu[0] = z[0] = 0;
        for (int i = 1; i < n - 1; i++) {
            l[i] = 2 * (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] = (a[j + 1] - a[j]) / h[j] - h[j] * (c[j + 1] + 2 * c[j]) / 3;
            d[j] = (c[j + 1] - c[j]) / (3 * h[j]);
        }
    }

    double evaluate(double t) {
        int i = 0, j = n - 1;
        while (i < j) {
            int k = (i + j) / 2;
            if (x[k] < t) {
                i = k + 1;
            } else {
                j = k;
            }
        }
        double dx = t - x[i];
        return a[i] + b[i] * dx + c[i] * dx * dx + d[i] * dx * dx * dx;
    }
};

int main() {
    vector<double> x = {0, 1, 2, 3, 4, 5};
    vector<double> y = {0, 1, 4, 9, 16, 25};
    Spline spline;
    spline.build(y);
    for (int i = 0; i < x.size(); i++) {
        cout << "x = " << x[i] << ", y = " << spline.evaluate(x[i]) << endl;
    }
    return 0;
}

在上述代码中，我们首先定义了一个结构体Spline，用于表示三次样条插值。其中，a、b、c、d分别表示三次样条函数的系数，x表示数据点的横坐标，n表示数据点的个数。

在函数build中，我们首先计算出数据点的横坐标x，然后根据三次样条插值的求解方法，解出三对角线矩阵的系数。具体而言，我们计算出向量h、alpha、l、mu、z，然后通过追赶法求解三对角线矩阵的系数。

在函数evaluate中，我们首先根据给定的横坐标t，找到其对应的区间[x[i], x[i+1])，然后根据三次样条函数的表达式计算出纵坐标。

最后，在main函数中，我们定义了一个数据点的横坐标x和纵坐标y，然后调用spline.build(y)求解三次样条插值，最后，通过spline.evaluate(x[i])计算出对应的纵坐标。