用C++实现的三次样条插值

 在写代码之前需要先了解什么是三次样条插值：

什么是三次样条插值？

所谓三次样条插值对于一个区间（a,b）将区间分成x0 = a < x1 ......xn-1 < b = xn 的n-1个区间，我们需要通过已知的n+1个点来模拟一个未知的函数，在三次样条插值中我们采用分段的方法来做这件事情。

三次样条插值得到的分段函数保证一下条件成立，而这些条件也是用来求解每一段样条插值的条件：

1 模拟出来的函数在已知点的函数值等于f的函数值

2模拟出来的分段函数是二阶连续的也就是说导数和二阶导数在分段的交界点是相等的（连续性）

3需要知道在a和b点的二阶导数的情况，或者二阶导数在这n+1个点的变化规律（凹凸性）

下面直接转载其内在的规律

已知:

a. n+1个数据点[xi, yi], i = 0, 1, …, n

b. 每一分段都是三次多项式函数曲线

c. 节点达到二阶连续

d. 左右两端点处特性（自然边界，固定边界，非节点边界）

根据定点，求出每段样条曲线方程中的系数，即可得到每段曲线的具体表达式。

 

插值和连续性：

, 其中 i = 0, 1, …, n-1

微分连续性：

 , 其中 i = 0, 1, …, n-2

样条曲线的微分式：

 

将步长 带入样条曲线的条件：

a. 由 (i = 0, 1, …, n-1)推出

b. 由 (i = 0, 1, …, n-1)推出

c. 由  (i = 0, 1, …, n-2)推出

由此可得：

d. 由  (i = 0, 1, …, n-2)推出

 

设 ，则

a.  可写为：

 ，推出

b. 将ci, di带入  可得：

c. 将bi, ci, di带入 (i = 0, 1, …, n-2)可得：

这样我们可以构造一个以m为未知数的线性方程组，而且在端点条件已知的情况下我们是知道其中几个未知数的值的

端点条件

由i的取值范围可知，共有n-1个公式， 但却有n+1个未知量m 。要想求解该方程组，还需另外两个式子。所以需要对两端点x0和xn的微分加些限制。 选择不是唯一的，3种比较常用的限制如下。

a. *边界(Natural)

首尾两端没有受到任何让它们弯曲的力，即 。具体表示为 和 

则要求解的方程组可写为：

 

 

b. 固定边界(Clamped)

首尾两端点的微分值是被指定的，这里分别定为A和B。则可以推出

将上述两个公式带入方程组，新的方程组左侧为

c. 非节点边界(Not-A-Knot)

指定样条曲线的三次微分匹配，即

根据 和 ，则上述条件变为

新的方程组系数矩阵可写为：

 

 

右下图可以看出不同的端点边界对样条曲线的影响：

 

1.3 算法总结

假定有n+1个数据节点

a. 计算步长 (i = 0, 1, …, n-1)

b. 将数据节点和指定的首位端点条件带入矩阵方程

c. 解矩阵方程，求得二次微分值。该矩阵为三对角矩阵，具体求法参见我的上篇文章：三对角矩阵的求解。

d. 计算样条曲线的系数：

其中i = 0, 1, …, n-1

e. 在每个子区间 中，创建方程

准备工作：

环境：Visul Studio 2019 、安装好C++矩阵计算工具eigen,安装方法点这里

 我写的是自然边界条件的三次样条插值代码

#include<iostream>
#include<vector>
#include "Eigen/Dense"
#include<cmath>
#include <fstream>

using namespace std;
using namespace Eigen;

int main() {
	const double min = 0;
	const double max = 15;
	double x[] = { 0, 3, 5, 7, 9, 11, 12, 13, 14, 15 };
	double y[] = { 0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6 };
	int size = sizeof(x) / sizeof(double);

	vector<double> xx, yy;
	double step = 0.1;
	double value = x[0];
	int num = (max - min) / step;
	for (double i = 0; i <= num; i++) {
		xx.push_back(value);
		value = value + step;
	}
	int size_xx = xx.size();

	vector<double> dx;
    vector<double> dy;
	for (int i = 0; i < size - 1; i++) {
		double temp_dx = x[i + 1] - x[i];
		dx.push_back(temp_dx);
		double temp_dy = y[i + 1] - y[i];
		dy.push_back(temp_dy);
	}

	MatrixXd H = MatrixXd::Random(size, size);
	for (int i = 0; i < size; i++) {
		for (int j = 0; j < size; j++) {
			H(i, j) = 0;
		}
	}
	VectorXd Y(size);
	for (int i = 0; i < size; i++) {
		Y(i) = 0;
	}
	VectorXd M(size);
	for (int i = 0; i < size; i++) {
		M(i) = 0;
	}

	H(0, 0) = 1;
	H(size - 1, size - 1) = 1;
	for (int i = 1; i < size - 1; i++) {
		H(i, i - 1) = dx[i - 1];
		H(i, i) = 2 * (dx[i - 1] + dx[i]);
		H(i, i + 1) = dx[i];
		Y(i) = 3 * (dy[i] / dx[i] - dy[i - 1] / dx[i - 1]);
	}

	M = H.colPivHouseholderQr().solve(Y);

	vector<double> ai, bi, ci, di;
	for (int i = 0; i < size - 1; i++) {
		ai.push_back(y[i]);
		di.push_back((M(i + 1) - M(i)) / (3 * dx[i]));
		bi.push_back(dy[i] / dx[i] - dx[i] * (2 * M(i) + M(i + 1)) / 3);
		ci.push_back(M(i));
	}

	vector<int> x_, xx_;
	for (int i = 0; i < size; i++) {
		int temp = x[i] / 0.1;
		x_.push_back(temp);
	}
	for (int i = 0; i < size_xx; i++) {
		int temp = xx[i] / 0.1;
		xx_.push_back(temp);
	}

	for (int i = 0; i < size_xx; i++) {
		int k = 0;
		for (int j = 0; j < size - 1; j++) {
			if (xx_[i] >= x_[j] && xx_[i] < x_[j + 1]) {
				k = j;
				break;
			}
			else if (xx[i] == x[size - 1]) {
				k = size - 1;
			}
		}
		//yy(i) = y[i] + bi(k) * (xx[i] - x[k]) + 1 / 2.0 * M(i) * pow((xx[i] - x[k]) , 2) + di(k) * pow((xx[i] - x[k]),3);
		double temp = ai[k] + bi[k] * (xx[i] - x[k]) + M(k) * pow((xx[i] - x[k]), 2) + di[k] * pow((xx[i] - x[k]), 3);
		yy.push_back(temp);
	}

	std::ofstream output;
	output.open("Spline.txt");
	for (unsigned i = 0; i < size_xx; i++) {
		output << xx[i] << '\t' << yy[i] << std::endl;
	}
	output.close();
}