插值与拟合 (一) ： 拉格朗日多项式插值 、Newton插值 、分段线性插值、Hermite插值 、样条插值、 B 样条函数插值、二维插值
插值：求过已知有限个数据点的近似函数。
拟合：已知有限个数据点，求近似函数，不要求过已知数据点，只要求在某种意义下它在这些点上的总偏差最小。
插值和拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二 者的数学方法上是完全不同的。而面对一个实际问题，究竟应该用插值还是拟合，有时容易确定，有时则并不明显。
下面介绍几种基本的、常用的插值：拉格朗日多项式插值、牛顿插值、分段线性插 值、Hermite 插值和三次样条插值。
目录
范德蒙特(Vandermonde)行列式 截断误差 / 插值余项
1.2 拉格朗日插值多项式 1.3 用 Matlab 作 Lagrange 插值
2.1 差商 : 定义与性质 2.2 Newton 插值公式
2.3 差分 ：向前差分、向后差分、中心差分 差分的两个性质
3.1 插值多项式的振荡 3.2 分段线性插值 3.3 用 Matlab 实现分段线性插值
4.1 Hermite 插值多项式 4.2 用 Matlab 实现 Hermite 插值
3 种类型的边界条件：完备/Lagrange 、自然边界条件、周期条件
5.4 三次样条插值在 Matlab 中的实现 例 1 机床加工
6.3 一维等距 B 样条函数插值 6.4 二维等距 B 样条函数插值
7.1 插值节点为网格节点 7.2 插值节点为散乱节点 习题
1 拉格朗日多项式插值
1.1 插值多项式
范德蒙特(Vandermonde)行列式
截断误差 / 插值余项
1.2 拉格朗日插值多项式
1.3 用 Matlab 作 Lagrange 插值
Matlab中没有现成的Lagrange插值函数，必须编写一个M文件实现Lagrange插值。 设n个节点数据以数组 x0 , y0 输入（注意 Matlat 的数组下标从 1 开始） ，m 个插值 点以数组 x输入，输出数组 y 为m 个插值。编写一个名为 lagrange.m 的 M 文件：
function y=lagrange(x0,y0,x);
n=length(x0);m=length(x);
for i=1:m
z=x(i);
s=0.0;
for k=1:n
p=1.0;
for j=1:n
if j~=k
p=p*(z-x0(j))/(x0(k)-x0(j));
end
end
s=p*y0(k)+s;
end
y(i)=s;
end
2 牛顿（Newton）插值
在导出 Newton 公式前，先介绍公式表示中所需要用到的差商、差分的概念及性质。
2.1 差商 : 定义与性质
2.2 Newton 插值公式
Newton 插值的优点
差商与导数的关系
2.3 差分 ：向前差分、向后差分、中心差分
当节点等距时，即相邻两个节点之差（称为步长）为常数，Newton 插值公式的形 式会更简单。此时关于节点间函数的平均变化率（差商）可用函数值之差（差分）来表 示。
差分的两个性质
（i）各阶差分均可表成函数值的线性组合，例如
（ii）各种差分之间可以互化。向后差分与中心差分化成向前差分的公式如下：
2.4 等距节点插值公式 、 Newton 向前插值公式
3 分段线性插值
3.1 插值多项式的振荡
高次插值多项式的这些缺陷，促使人们转而寻求简单的低次多项式插值。
3.2 分段线性插值
用 计算 x点的插值时，只用到 x左右的两个节点，计算量与节点个数n无关。 但n越大，分段越多，插值误差越小。实际上用函数表作插值计算时，分段线性插值就足够了，如数学、物理中用的特殊函数表，数理统计中用的概率分布表等。
3.3 用 Matlab 实现分段线性插值
用 Matlab 实现分段线性插值不需要编制函数程序，Matlab 中有现成的一维插值函 数 interp1。
y=interp1(x0,y0,x,'method') method 指定插值的方法，默认为线性插值。其值可为：
-
'nearest' 最近项插值
-
'linear' 线性插值
-
'spline' 逐段 3 次样条插值
-
'cubic' 保凹凸性 3 次插值
所有的插值方法要求 x0 是单调的。 当 x0 为等距时可以用快速插值法，使用快速插值法的格式为'*nearest'、'*linear'、 '*spline'、'*cubic'。
4 埃尔米特(Hermite)插值
4.1 Hermite 插值多项式
如果对插值函数，不仅要求它在节点处与函数同值，而且要求它与函数有相同的一 阶、二阶甚至更高阶的导数值，这就是 Hermite 插值问题。本节主要讨论在节点处插值 函数与函数的值及一阶导数值均相等的 Hermite 插值。
4.2 用 Matlab 实现 Hermite 插值
Matlab 中没有现成的 Hermite 插值函数，必须编写一个 M 文件实现插值。
function y=hermite(x0,y0,y1,x);
n=length(x0);m=length(x);
for k=1:m
yy=0.0;
for i=1:n
h=1.0;
a=0.0;
for j=1:n
if j~=i
h=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2;
a=1/(x0(i)-x0(j))+a;
end
end
yy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i));
end
y(k)=yy;
end
5 样条插值
许多工程技术中提出的计算问题对插值函数的光滑性有较高要求，如飞机的机翼外 形，内燃机的进、排气门的凸轮曲线，都要求曲线具有较高的光滑程度，不仅要连续， 而且要有连续的曲率，这就导致了样条插值的产生。
5.1 样条函数的概念
所谓样条（Spline）本来是工程设计中使用的一种绘图工具，它是富有弹性的细木 条或细金属条。绘图员利用它把一些已知点连接成一条光滑曲线（称为样条曲线），并使连接点处有连续的曲率。
内节点 、边界点、k 次样条函数空间
二次样条函数
三次样条函数
利用样条函数进行插值，即取插值函数为样条函数，称为样条插值。例如分段线性插值 是一次样条插值。下面我们介绍二次、三次样条插值。
5.2 二次样条函数插值
两类问题
证明这两类插值问题都是唯一可解的
5.3 三次样条函数插值
3 种类型的边界条件：完备/Lagrange 、自然边界条件、周期条件
5.4 三次样条插值在 Matlab 中的实现
在 Matlab 中数据点称之为断点。如果三次样条插值没有边界条件，最常用的方法， 就是采用非扭结（not-a-knot）条件。这个条件强迫第 1 个和第 2 个三次多项式的三阶 导数相等。对最后一个和倒数第 2 个三次多项式也做同样地处理。
Matlab 中三次样条插值也有现成的函数：
y=interp1(x0,y0,x,'spline')；
y=spline(x0,y0,x)；
pp=csape(x0,y0,conds)，y=ppval(pp,x)其中 x0,y0 是已知数据点，x 是插值点，y 是插值点的函数值。 对于三次样条插值，我们提倡使用函数 csape，csape 的返回值是 pp 形式，要求出插值点的函数值，必须调用函数 ppval。
pp=csape(x0,y0)：使用默认的边界条件，即 Lagrange 边界条件。
pp=csape(x0,y0,conds)中的 conds 指定插值的边界条件，其值可为：
- 'complete' 边界为一阶导数，即默认的边界条件
- 'not-a-knot' 非扭结条件
- 'periodic' 周期条件
- 'second' 边界为二阶导数，二阶导数的值[0, 0]。
- 'variational' 设置边界的二阶导数值为[0,0]。
对于一些特殊的边界条件，可以通过 conds 的一个 1× 2 矩阵来表示，conds 元素的 取值为 1，2。此时，使用命令
pp=csape(x0,y0_ext,conds) 其中 y0_ext=[left, y0, right]，这里 left 表示左边界的取值，right 表示右边界的取值。
conds(i)=j 的含义是给定端点i的 j 阶导数，即 conds 的第一个元素表示左边界的条 件，第二个元素表示右边界的条件；
conds=[2,1]表示左边界是二阶导数，右边界是一阶 导数，对应的值由 left 和 right 给出。
详细情况请使用帮助 help csape。
例 1 机床加工
解 编写以下程序：
clc,clear
x0=[0 3 5 7 9 11 12 13 14 15];
y0=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6];
x=0:0.1:15;
y1=lagrange(x0,y0,x); %调用前面编写的Lagrange插值函数
y2=interp1(x0,y0,x);
y3=interp1(x0,y0,x,'spline');
pp1=csape(x0,y0);
y4=ppval(pp1,x);
pp2=csape(x0,y0,'second');
y5=ppval(pp2,x);
fprintf('比较一下不同插值方法和边界条件的结果:\n')
fprintf('x y1 y2 y3 y4 y5\n')
xianshi=[x',y1',y2',y3',y4',y5'];
fprintf('%f\t%f\t%f\t%f\t%f\t%f\n',xianshi')
subplot(2,2,1), plot(x0,y0,'+',x,y1), title('Lagrange')
subplot(2,2,2), plot(x0,y0,'+',x,y2), title('Piecewise linear')
subplot(2,2,3), plot(x0,y0,'+',x,y3), title('Spline1')
subplot(2,2,4), plot(x0,y0,'+',x,y4), title('Spline2')
dyx0=ppval(fnder(pp1),x0(1)) %求x=0处的导数
ytemp=y3(131:151);
index=find(ytemp==min(ytemp));
xymin=[x(130+index),ytemp(index)] 计算结果略。 可以看出，拉格朗日插值的结果根本不能应用，分段线性插值的光滑性较差（特别 是在x =14 附近弯曲处），建议选用三次样条插值的结果。
6 B 样条函数插值方法
6.1 磨光函数
实际中的许多问题，往往是既要求近似函数（曲线或曲面）有足够的光滑性，又要 求与实际函数有相同的凹凸性，一般插值函数和样条函数都不具有这种性质。如果对于 一个特殊函数进行磨光处理生成磨光函数（多项式），则用磨光函数构造出样条函数作 为插值函数，既有足够的光滑性，而且也具有较好的保凹凸性，因此磨光函数在一维插 值（曲线）和二维插值（曲面）问题中有着广泛的应用。 由积分理论可知，对于可积函数通过积分会提高函数的光滑度，因此，我们可以利 用积分方法对函数进行磨光处理。
6.2 等距 B 样条函数
6.3 一维等距 B 样条函数插值
等距 B 样条函数与通常的样条有如下的关系：
6.4 二维等距 B 样条函数插值
7 二维插值
前面讲述的都是一维插值，即节点为一维变量，插值函数是一元函数（曲线）。若 节点是二维的，插值函数就是二元函数，即曲面。如在某区域测量了若干点（节点）的 高程（节点值），为了画出较精确的等高线图，就要先插入更多的点（插值点），计算这些点的高程（插值）。
7.1 插值节点为网格节点
Matlab 中有一些计算二维插值的程序。如
z=interp2(x0,y0,z0,x,y,'method') 如果是三次样条插值，可以使用命令
pp=csape({x0,y0},z0,conds,valconds)，z=fnval(pp,{x,y}) 解 编写程序如下：
clear,clc
x=100:100:500;
y=100:100:400;
z=[636 697 624 478 450
698 712 630 478 420
680 674 598 412 400
662 626 552 334 310];
pp=csape({x,y},z')
xi=100:10:500; yi=100:10:400
cz1=fnval(pp,{xi,yi})
cz2=interp2(x,y,z,xi,yi','spline')
[i,j]=find(cz1==max(max(cz1)))
x=xi(i),y=yi(j),zmax=cz1(i,j)
7.2 插值节点为散乱节点
对上述问题，Matlab 中提供了插值函数 griddata，其格式为：
ZI = GRIDDATA(X,Y,Z,XI,YI) 例 3 在某海域测得一些点（x,y）处的水深 z 由下表给出，在矩形区域（75,200） ×(-50,150) 内画出海底曲面的图形。
解 编写程序如下：
x=[129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.5];
y=[7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5];
z=-[4 8 6 8 6 8 8 9 9 8 8 9 4 9];
xi=75:1:200;
yi=-50:1:150;
zi=griddata(x,y,z,xi,yi','cubic')
subplot(1,2,1), plot(x,y,'*')
subplot(1,2,2), mesh(xi,yi,zi) 习题
上一篇： VSCode的好用的插件