三次样条插值 C++实现
时间: 2023-07-13 08:31:00 浏览: 105
以下是一个简单的C++实现三次样条插值的示例代码:
```c++
#include <iostream>
#include <vector>
#include <cmath>
using namespace std;
vector<double> spline(const vector<double>& x, const vector<double>& y)
{
int n = x.size();
vector<double> a(n), b(n), c(n), d(n), h(n), alpha(n), l(n), mu(n), z(n), c_new(n);
for(int i = 1; i < n; i++)
{
h[i] = x[i] - x[i-1];
alpha[i] = 3 * (y[i] - y[i-1]) / h[i] - 3 * (y[i-1] - y[i-2]) / h[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;
c[n-1] = 0;
for(int j = n-2; j >= 0; j--)
{
c[j] = z[j] - mu[j] * c[j+1];
b[j] = (y[j+1] - y[j]) / h[j] - h[j] * (c[j+1] + 2 * c[j]) / 3;
d[j] = (c[j+1] - c[j]) / (3 * h[j]);
a[j] = y[j];
}
return vector<double>{a, b, c, d};
}
double spline_eval(double x, const vector<double>& a, const vector<double>& b, const vector<double>& c, const vector<double>& d, const vector<double>& xs)
{
int n = xs.size();
int k = 0;
for(int i = 1; i < n; i++)
{
if(x <= xs[i])
{
k = i-1;
break;
}
}
double dx = x - xs[k];
return a[k] + b[k] * dx + c[k] * dx * dx + d[k] * dx * dx * dx;
}
int main()
{
vector<double> x{0, 1, 2, 3, 4};
vector<double> y{0, 1, 4, 9, 16};
auto params = spline(x, y);
auto a = params[0];
auto b = params[1];
auto c = params[2];
auto d = params[3];
for(double x_val = 0; x_val <= 4; x_val += 0.1)
{
cout << "f(" << x_val << ") = " << spline_eval(x_val, a, b, c, d, x) << endl;
}
return 0;
}
```
在上面的代码中,spline函数的输入参数为x和y,表示需要插值的数据点的横坐标和纵坐标。该函数返回一个包含四个元素的vector,分别表示插值函数的四个系数。spline_eval函数的输入参数为x、a、b、c、d和xs,分别表示需要求解的横坐标、插值函数的四个系数和插值区间的横坐标。该函数返回对应于x的插值函数的值。
在示例代码中,我们使用了一个简单的二次函数f(x) = x^2进行插值,并输出了在插值区间内的所有取值。