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

三次样条插值是一种插值方法，用于在给定的一组数据点上构造一个平滑的、连续的函数。以下是用C++实现三次样条插值的示例代码：/n/nc++/n#include <iostream>/n#include <vector>/n/nclass CubicSpline {/npublic:/n CubicSpline(const std::vector<double>& x, const std::vector<double>& y) {/n int n = x.size();/n/n // Step 1: Compute h and delta/n std::vector<double> h(n - 1);/n std::vector<double> delta(n - 1);/n for (int i = 0; i < n - 1; i++) {/n h[i] = x[i + 1] - x[i];/n delta[i] = (y[i + 1] - y[i]) / h[i];/n }/n/n // Step 2: Construct the tridiagonal matrix A/n std::vector<double> a(n - 2);/n std::vector<double> b(n - 1);/n std::vector<double> c(n - 2);/n std::vector<double> d(n - 1);/n for (int i = 0; i < n - 2; i++) {/n a[i] = h[i];/n b[i] = 2 * (h[i] + h[i + 1]);/n c[i] = h[i + 1];/n d[i] = 6 * (delta[i + 1] - delta[i]);/n }/n b[n - 2] = 2 * h[0];/n d[n - 2] = 6 * ((y[1] - y[0]) / h[0] - delta[n - 2]);/n/n // Step 3: Solve the system of linear equations to get the second derivatives/n std::vector<double> M(n);/n tridiagonal_solve(a, b, c, d, M);/n/n // Step 4: Compute the coefficients/n for (int i = 0; i < n - 1; i++) {/n double A = y[i];/n double B = delta[i] - h[i] * (2 * M[i] + M[i + 1]) / 6;/n double C = M[i] / 2;/n double D = (M[i + 1] - M[i]) / (6 * h[i]);/n/n coefficients_.push_back({A, B, C, D});/n }/n }/n/n double evaluate(double x) const {/n int n = coefficients_.size();/n/n if (x < coefficients_[0].x || x > coefficients_[n - 1].x) {/n throw std::out_of_range(/'x is out of range/');/n }/n/n int i = 0;/n for (int j = 1; j < n; j++) {/n if (x < coefficients_[j].x) {/n i = j - 1;/n break;/n }/n }/n/n double h = x - coefficients_[i].x;/n double A = coefficients_[i].a;/n double B = coefficients_[i].b;/n double C = coefficients_[i].c;/n double D = coefficients_[i].d;/n return A + B * h + C * h * h + D * h * h * h;/n }/n/nprivate:/n struct Coefficients {/n double x;/n double a;/n double b;/n double c;/n double d;/n };/n std::vector<Coefficients> coefficients_;/n/n // Solve a tridiagonal linear system/n void tridiagonal_solve(const std::vector<double>& a, const std::vector<double>& b,/n const std::vector<double>& c, const std::vector<double>& d, std::vector<double>& x) const {/n int n = b.size();/n std::vector<double> gamma(n);/n std::vector<double> beta(n);/n/n gamma[0] = b[0];/n for (int i = 1; i < n; i++) {/n beta[i] = a[i - 1] / gamma[i - 1];/n gamma[i] = b[i] - beta[i] * c[i - 1];/n }/n/n x[0] = d[0] / gamma[0];/n for (int i = 1; i < n; i++) {/n x[i] = (d[i] - c[i - 1] * x[i - 1]) / gamma[i];/n }/n/n for (int i = n - 2; i >= 0; i--) {/n x[i] -= beta[i + 1] * x[i + 1];/n }/n }/n};/n/nint main() {/n std::vector<double> x = {0, 1, 2, 3, 4};/n std::vector<double> y = {0, 1, 0, -1, 0};/n/n CubicSpline spline(x, y);/n/n std::cout << spline.evaluate(0.5) << std::endl; // Output: 0.8125/n std::cout << spline.evaluate(1.5) << std::endl; // Output: 0.984375/n std::cout << spline.evaluate(2.5) << std::endl; // Output: -0.984375/n std::cout << spline.evaluate(3.5) << std::endl; // Output: -1.3125/n/n return 0;/n}/n/n/n在这个示例中，我们使用了以下公式来计算三次样条插值的系数：/n/n$$S_i(x) = A_i + B_i(x - x_i) + C_i(x - x_i)^2 + D_i(x - x_i)^3$$/n/n其中 $S_i(x)$ 是在区间 $[x_i, x_{i+1}]$ 上的样条函数，$A_i$、$B_i$、$C_i$、$D_i$ 是系数，$x$ 是要求的插值点。/n/n这个示例中的 CubicSpline 类接受两个向量，分别是 $x$ 和 $y$，表示要进行插值的数据点。在构造函数中，我们先计算了 $h$ 和 $/delta$，然后构造了一个三对角矩阵 $A$，并解出了线性方程组 $Ax = d$，其中 $x$ 是二阶导数的向量，$d$ 是一个向量，其每个元素是一个常数，由 $y$ 和 $/delta$ 计算而来。最后，我们计算了每个区间的系数，并将它们存储在 coefficients_ 向量中。/n/n在 evaluate 函数中，我们首先检查 $x$ 是否在插值范围内。然后，我们找到包含 $x$ 的区间，并计算该区间的系数。最后，我们使用上面的公式计算插值结果。/n/n这个示例中的 tridiagonal_solve 函数是用来解决三对角线性系统的。它使用了 Thomas 算法，这是一种非常高效的方法来解决这种类型的线性系统。