C++中实现三次样条插值

### C++ 实现三次样条插值算法 三次样条插值是一种用于平滑数据点之间关系的方法，在数值分析中有广泛应用。下面是一个完整的C++程序来实现这一功能。 #### 定义结构体存储节点信息 为了方便操作，定义一个`Node`结构用来保存各个节点的位置及其对应的函数值： ```cpp struct Node { double x; // 节点横坐标 double y; // 节点纵坐标 }; ``` #### 计算二阶导数数组 根据给定的边界条件计算每个内部节点处的一维二次微分系数\[b_i\], \[c_i\] 和 \[d_i\][^1]: ```cpp void calculateSecondDerivatives(const std::vector<Node>& nodes, std::vector<double>& d2ydx2) { int n = nodes.size(); // 初始化矩阵A和向量B Matrix A(n - 2); Vector B(n - 2); for (int i = 0; i < n - 2; ++i) { A[i][i] = 4; if (i != 0) { A[i][i-1] = 1; A[i-1][i] = 1; } B[i] = 3 * ((nodes[i+2].y - nodes[i+1].y)/(nodes[i+2].x - nodes[i+1].x) - (nodes[i+1].y - nodes[i].y )/(nodes[i+1].x - nodes[i].x)); } // 解线性方程组 Ay=B 得到中间结点处的二阶导数 solveTridiagonalSystem(A, B, &d2ydx2); // 设置端点处的自然边界条件（即两端点处的二阶导数为零） d2ydx2.insert(d2ydx2.begin(), 0.); d2ydx2.push_back(0.); } ``` 这里使用了追赶法求解三对角系统的线性代数子过程`solveTridiagonalSystem()`，具体细节可以参见其他资料[^2]。 #### 插值函数构建 有了上述准备之后就可以编写实际执行插值运算的部分了: ```cpp double interpolate(double targetX, const std::vector<Node>& nodes, const std::vector<double>& d2ydx2) { auto it = lower_bound(nodes.begin(), nodes.end(), Node{targetX, 0}, [](const Node& a, const Node& b){return a.x<b.x;}); --it; size_t idx = distance(nodes.begin(), it); double h = nodes[idx + 1].x - nodes[idx].x; double t = (targetX - nodes[idx].x)/h; return ( pow((1-t), 3)*nodes[idx].y + pow(t, 3)*nodes[idx + 1].y + h*h*(pow((1-t), 3)*(t+.5)+pow(t, 3)*(-t+.5))*d2ydx2[idx]/6 + h*h*(-pow((1-t), 3)*t+pow(t, 3)*(1-.5*t))*d2ydx2[idx + 1]/6 ); } ``` 这段代码实现了基于输入的目标位置\(x\)找到最近一对相邻控制点，并利用这些点的信息以及预先计算好的二阶导数来进行最终的结果预测。 #### 主函数调用示例 最后给出如何组织以上各部分完成整个流程的一个简单例子： ```cpp #include <iostream> #include <algorithm> #include <cmath> // ... 上述所有辅助函数声明 ... int main() { // 输入样本点集 std::vector<Node> samplePoints = {{0., 0.}, {1., .8}, {2., .9}, {3., .1}}; // 排序并去重 sort(samplePoints.begin(), samplePoints.end(), [](const Node& lhs, const Node& rhs){ return fabs(lhs.x-rhs.x)>1e-7 && lhs.x<rhs.x; }); unique(samplePoints.begin(), samplePoints.end(), [](const Node& lhs, const Node& rhs){ return fabs(lhs.x-rhs.x)<1e-7; }); // 准备工作：获取二阶导数序列 vector<double> secondDervs; calculateSecondDerivatives(samplePoints, secondDervs); // 测试查询若干个新位置上的估计值 cout << "Interpolated values:" << endl; for(int k=0;k<=10;++k){ double queryPt=(samplePoints.front().x*k+samplePoints.back().x*(10-k))/10.; printf("%.2f -> %.4fn",queryPt ,interpolate(queryPt,samplePoints,secondDervs)); } return EXIT_SUCCESS; } ``` 此段代码展示了怎样读取一组离散的数据点作为训练集合，接着通过排序去除重复项确保后续处理的有效性；随后调用了之前提到过的几个关键组件完成了从原始观测到任意指定坐标的转换任务。

好的，我可以帮你实现分段三次样条插值。以下是一个简单的 C++ 代码示例： ```c++ #include <iostream> #include <vector> using namespace std; // 定义一个结构体用于存储插值系数 struct SplineCoefficients { double a, b, c, d; }; // 定义一个函数用于计算分段三次样条插值的系数 vector<SplineCoefficients> calculateSplineCoefficients(const vector<double>& x, const vector<double>& y) { int n = x.size(); vector<double> h(n - 1), alpha(n - 1), l(n), mu(n - 1), z(n); vector<SplineCoefficients> coefficients(n - 1); for (int i = 0; i < n - 1; i++) { h[i] = x[i + 1] - x[i]; } for (int i = 1; i < n - 1; i++) { alpha[i] = 3.0 / h[i] * (y[i + 1] - y[i]) - 3.0 / h[i - 1] * (y[i] - y[i - 1]); } l[0] = 1; mu[0] = 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; coefficients[n - 2].c = 0; for (int i = n - 2; i >= 0; i--) { coefficients[i].c = z[i] - mu[i] * coefficients[i + 1].c; coefficients[i].b = (y[i + 1] - y[i]) / h[i] - h[i] * (coefficients[i + 1].c + 2 * coefficients[i].c) / 3.0; coefficients[i].d = (coefficients[i + 1].c - coefficients[i].c) / (3.0 * h[i]); coefficients[i].a = y[i]; } return coefficients; } // 定义一个函数用于计算分段三次样条插值的值 double evaluateSpline(const SplineCoefficients& coefficients, double x) { double dx = x - coefficients.a; return coefficients.a + coefficients.b * dx + coefficients.c * dx * dx + coefficients.d * dx * dx * dx; } int main() { vector<double> x = {1, 2, 4, 5}; vector<double> y = {3, 5, 4, 6}; vector<SplineCoefficients> coefficients = calculateSplineCoefficients(x, y); for (int i = 0; i < coefficients.size(); i++) { cout << "f(x) = " << coefficients[i].a << " + " << coefficients[i].b << "(x - " << coefficients[i].a << ") + " << coefficients[i].c << "(x - " << coefficients[i].a << ")^2 + " << coefficients[i].d << "(x - " << coefficients[i].a << ")^3" << endl; } cout << endl; for (double i = 1; i <= 5; i += 0.1) { for (int j = 0; j < coefficients.size(); j++) { if (i >= x[j] && i <= x[j + 1]) { cout << "f(" << i << ") = " << evaluateSpline(coefficients[j], i) << endl; } } } return 0; } ``` 这个代码实现了一个简单的分段三次样条插值，并且在主函数中示范了如何使用它来计算插值函数在某个点的值。你可以根据你的具体需求修改代码。