数值计算笔记之插值（四）三次样条插值

0、定义

已知函数在区间上个互异节点，处的函数值为,若构造函数,满足：

1. 
2. 在每个小区间上是一个不超过三次的多项式
3. 在上连续

则称为的三次样条插值函数。

根据定义知道规律为：

已知:

• n+1个数据点[xi, yi], i = 0, 1, …, n
•  每一分段都是三次多项式函数曲线
•  节点达到二阶连续
• 左右两端点处特性（自然边界，固定边界，非节点边界）

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

一、推导

插值和连续性：

  其中.

微分连续性：

其中.

样条曲线的微分式：

 

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

1.     
2. 
3.                                                                                                      

   得

   

4. ，.

• 设  ，那么  可改写为 
• 将  带入 ，得 
• 将  代入，得

至此，可以得到一个关于m为未知数的线性方程组，通过解这个方程组，确定其系数。

• 
• 
• 
• 

二、边界条件建立方程

边界条件有三种：*边界条件、固定边界条件、非节点边界条件。

1、*边界条件

两端没有任何限制，即。求解方程为：

2、固定边界条件

数据两端节点的微分值是已知的，令分别为A 和B。则有

则有

3、非节点边界条件

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

由  . 

则有

 

三、代码实现

1、matlab

在matlab中，spline函数可以实现三次样条插值。例：

clc;
clear;
//x = [ 0, 3, 5, 7, 9, 11, 12, 13, 14, 15 ];
//y = [0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6 ];
//x_ = 0:.1:15;
x = -pi:pi;
y = sin(x);
x_ = -pi:.1:pi;
p1 = spline(x,y,x_); % 分段三次样条插值
plot(x,y,'ko',x_,p1,'b.')
legend('插值节点','分段三次样条插值','location','southeast')

 可以看出，由点组成得曲线比较平滑。

2、C++

程序我有用到矩阵库Armadillo，没有安装的请点这个安装教程链接

#include<iostream>
#include<armadillo>
#include<vector>
using namespace std;
using namespace arma;

int main()
{
	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);
	const double min = x[0];
	const double max = x[size - 1];

	vector<double> xx, yy;              //插值点与计算的函数值
	for (double i = x[0]; i <= x[size - 1]; i = i + 0.1) {
		xx.push_back(i);
	}
	int size_xx = xx.size();

	vector<double> dx, 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);
	}

	mat H, Y, M;           // H * M = Y
	H.zeros(size, size);
	Y.zeros(size,1);
	//M.zeros(1, size);

	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 = solve(H,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();
	cout << "Hello World !" << endl;
}

结果用matlab显示为：

四、小结

1、算法流程

1. 计算步长  
2. 根据边界条件填充矩阵 
3. 解矩阵方程 
4. 由  得到样条插值函数的系数

• 
• 
• 
• 

2、难点

三次样条插值的整体推导，以及 H 矩阵的填充。