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

0、定义

已知函数 $f(x)$ 在区间 $[a,b]$ 上 $n+1$ 个互异节点， $a=x_{0}<x_{1}<\cdots < x_{n}=b$ 处的函数值为 $y_{i}=f(x_{i})$ ,若构造函数 $s(x)$ ,满足：

$s(x_{i})=y_{i}(i=0,1,\cdots ,n)$
在每个小区间 $[x_{i},x_{i+1}]$ 上是一个不超过三次的多项式
$s(x),s^{'}(x),s^{''}(x)$ 在 $[a,b]$ 上连续

则称 $s(x)$ 为 $f(x)$ 的三次样条插值函数。

根据定义知道规律为：

已知:

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

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

一、推导

插值和连续性：

$S_{i}(x_{i}) = y_{i}$

$S_{i}(x_{i+1}) = y_{i+1}$

其中 $i = 0, 1,\cdots , n-1$ .

微分连续性：

$S^{'}_{i}(x_{i+1})=S^{'}_{i+1}(x_{i+1})$

$S^{''}_{i}(x_{i+1}) = S^{''}_{i+1}(x_{i+1})$

其中 $i = 0, 1,\cdots , n-2$ .

样条曲线的微分式：

$S_{i}(x) = a_{i}+b_{i}(x-x_{i})+c_{i}(x-x_{i})^{2}+d_{i}(x-x_{i})^{3}$

$S^{'}_{i}(x) = b_{i}+2c_{i}(x-x_{i})+3d_{i}(x-x_{i})^{2}$

$S^{''}_{i} = 2c_{i}+6d_{i}(x-x_{i})$

将步长 $h_{i} = x_{i+1}-x_{i}$ 带入样条曲线的条件：

$S_{i}(x_{i}) = y_{i}$ $(i=0,1,\cdots ,n-1)$ $\Rightarrow$ $a_{i} = y_{i}$
$S_{i}(x_{i+1}) = y_{i+1}(i=0,1,\cdots ,n-1)\Rightarrow a_{i}+h_{i}b_{i}+h_{i}^{2}c_{i}+h_{i}^{3}d_{i}=y_{i+1}$
$S^{'}_{i}(x_{i+1})=S^{'}_{i+1}(x_{i+1})(i=0,1,\cdots ,n-2)\Rightarrow S^{'}_{i}(x_{i+1})=b_{i}+2c_{i}(x_{i+1}-x_{i})$ $+3d_{i}(x_{i+1}-x_{i})^{2}=b_{i}+2c_{i}h_{i}+3d_{i}h^{2}_{i}$ $S^{'}_{i+1}(x_{i+1})=b_{i+1}+2c_{i}(x_{i+1}-x_{i+1})+3d_{i}(x_{i+1}-x_{i+1})^{2}=b_{i+1}$
$\therefore$ 得

$b_{i}+2h_{i}c_{i}+3h^{2}_{i}d_{i}-b_{i-1}=0$
$S^{''}_{i}(x_{i+1}) = S^{''}_{i+1}(x_{i+1})$ ， $i = 0, 1,\cdots , n-2$ . $\Rightarrow 2c_{i}+6h_{i}d_{i}-2c_{i+1}=0$

设 $m_{i}=S^{''}_{i}(x_{i})=2c_{i}$ ，那么 $2c_{i}+6h_{i}d_{i}-2c_{i+1}=0$ 可改写为 $m_{i}+6h_{i}d_{i}-m_{i+1}=0\Rightarrow d_{i}=\frac{m_{i+1}-m_{i}}{6h_{i}}$
将 $c_{i},d_{i}$ 带入 $y_{i}+h_{i}b_{i}+h_{i}^{2}c_{i}+h_{i}^{3}d_{i}=y_{i+1} (a_{i}=y_{i})$ ，得 $b_{i}=\frac{y_{i+1}-y_{i}}{h_{i}} -\frac{h_{i}}{2}m_{i}-\frac{h_{i}}{6}(m_{i+1}-m_{i})$
将 $b_{i},c_{i},d_{i}$ 代入 $b_{i}+2h_{i}c_{i}+3h^{2}_{i}d_{i}-b_{i-1}=0 (i=0,1,\cdots ,n-2)$ ，得 $h_{i}m_{i}+2(h_{i}+h_{i+1})m_{i+1}+h_{i+1}m_{i+2}=6(\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{y_{i+1}-y_{i}}{h_{i}})$

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

$a_{i}=y_{i}$
$b_{i}=\frac{y_{i+1}-y_{i}}{h_{i}} -\frac{h_{i}}{2}m_{i}-\frac{h_{i}}{6}(m_{i+1}-m_{i})$
$c_{i}=\frac{m_{i}}{2}$
$d_{i}=\frac{m_{i+1}-m_{i}}{6h_{i}}$

二、边界条件建立方程

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

1、自由边界条件

两端没有任何限制，即 $S^{''}=0,m_{0}=0,m_{n}=0$ 。求解方程为：

$\begin{bmatrix} 1 &0 &0 & & \cdots &0 \\ h_{0} &2(h_{0}+h_{1}) &h_{1} &0 & & \\ 0 & h_{1} &2(h_{1}+h_{2}) &h_{2} &0 & \\ 0 & 0 &h_{2} &2(h_{2}+h_{3}) &h_{3} &\vdots \\ \vdots & & \ddots & \ddots & \ddots & \\ 0 & \cdots & &0 & 0&1 \end{bmatrix} \begin{bmatrix} m_{0}\\ m_{1}\\ m_{2}\\ m_{3}\\ \vdots \\ m_{n} \end{bmatrix} = 6\begin{bmatrix} 0\\ \frac{y_{2}-y_{1}}{h_{1}}-\frac{y_{1}-y_{0}}{h_{0}} \\ \frac{y_{3}-y_{2}}{h_{2}}-\frac{y_{2}-y_{1}}{h_{1}} \\ \frac{y_{4}-y_{3}}{h_{3}}-\frac{y_{3}-y_{2}}{h_{2}} \\ \vdots \\ \frac{y_{n}-y_{n-1}}{h_{n-1}}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}} \\ 0 \end{bmatrix}$

2、固定边界条件

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

$S^{'}_{0}(x_{0})=A$

$\Rightarrow b_{0}=A$

$\Rightarrow A=\frac{y_{1}-y_{0}}{h_{0}}-\frac{h_{0}}{2}m_{0}-\frac{h_{0}}{6}(m_{1}-m_{0})$

$\Rightarrow 2h_{0}m_{0}+h_{0}m_{1} = 6\begin{bmatrix} \frac{y_{1}-y_{0}}{h_{0}}-A \end{bmatrix}$

$S^{'}_{n-1}(x_{n})=B$

$\Rightarrow b_{n-1}=B$

$\Rightarrow h_{n-1}m_{n-1}+2h_{n-1}m_{n}=6\begin{bmatrix} B-\frac{y_{n}-y_{n-1}}{h_{n-1}} \end{bmatrix}$

则有

$\begin{bmatrix} 2h_{0}& h_{0} &0 &\cdots & \cdots &0 \\ h_{0} & 2(h_{0}+h_{1}) & h_{1}&0 & &\vdots \\ 0 & h_{1} &2(h_{1}+h_{2}) &h_{2} &0 & \vdots \\ \vdots &0 &\ddots & \ddots & \ddots &0 \\ 0& \cdots &0 & h_{n-2} &2(h_{n-2}+h_{n-1}) &h_{n-1} \\ 0& \cdots &\cdots &0 &h_{n-1} & 2h_{n-1} \end{bmatrix} \begin{bmatrix} m_{0}\\ m_{1}\\ m_{2}\\ m_{3}\\ \vdots \\ m_{n} \end{bmatrix} = 6\begin{bmatrix} 0\\ \frac{y_{2}-y_{1}}{h_{1}}-\frac{y_{1}-y_{0}}{h_{0}} \\ \frac{y_{3}-y_{2}}{h_{2}}-\frac{y_{2}-y_{1}}{h_{1}} \\ \frac{y_{4}-y_{3}}{h_{3}}-\frac{y_{3}-y_{2}}{h_{2}} \\ \vdots \\ \frac{y_{n}-y_{n-1}}{h_{n-1}}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}} \\ 0 \end{bmatrix}$

3、非节点边界条件

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

$S^{'''}_{0}(x_{1})=S^{'''}_{1}(x_{1})$

$S^{'''}_{n-2}(x_{n-1})=S^{'''}_{n-1}(x_{n-1})$

由 $S^{'''}(x)=6d_{i},d_{i}=\frac{m_{i+1}-m_{i}}{6h_{i}}$ . $\Rightarrow$ $h_{1}(m_{1}-m_{0})=h_{0}(m_{2}-m_[1]),h_{n-1}(m_{n-1}-m_{n-2})=h_{n-2}(m_{n}-m_{n-1})$

则有

$\begin{bmatrix} -h_{1}& h_{0}+h_{1} &-h_{0} &\cdots & \cdots &0 \\ h_{0} & 2(h_{0}+h_{1}) & h_{1}&0 & &\vdots \\ 0 & h_{1} &2(h_{1}+h_{2}) &h_{2} &0 & \vdots \\ \vdots &0 &\ddots & \ddots & \ddots &0 \\ 0& \cdots &0 & h_{n-2} &2(h_{n-2}+h_{n-1}) &h_{n-1} \\ 0& \cdots &\cdots &-h_{n-1} &h_{n-2}+h_{n-1} & -h_{n-2} \end{bmatrix}$

三、代码实现

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 = YH.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、算法流程

计算步长 $h_{i}=x_{i+1}-x_{i}$
根据边界条件填充矩阵 $H$
解矩阵方程 $M = Y / H$
由 $m_{i}$ 得到样条插值函数的系数

$a_{i}=y_{i}$
$b_{i}=\frac{y_{i+1}-y_{i}}{h_{i}} -\frac{h_{i}}{2}m_{i}-\frac{h_{i}}{6}(m_{i+1}-m_{i})$
$c_{i}=\frac{m_{i}}{2}$
$d_{i}=\frac{m_{i+1}-m_{i}}{6h_{i}}$

2、难点

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