2024美赛数学建模常用数学建模模型之——插值法

插值：求过已知有限个数据点的近似函数。

拟合：已知有限个数据点，求近似函数，不要求过已知数据点，只要求在某种意义 下它在这些点上的总偏差最小。

插值和拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。而面对一个实际问题，究竟应该用插值还是拟合，有时容易确定，有时则并不明显。

1 插值方法

下面介绍几种基本的、常用的插值：拉格朗日多项式插值、牛顿插值、分段线性插值、Hermite 插值和三次样条插值。

1.1 拉格朗日多项式插值

1.1.1 插值多项式

用多项式作为研究插值的工具，称为代数插值。其基本问题是：已知函数 f (x)在 区间[a,b]上n +1个不同点 x0 , x1 ,L, xn 处的函数值 yi = f (xi) (i = 0,1,L,n) ，求一个 至多n 次多项式

1.1.2 拉格朗日插值多项式

实际上比较方便的作法不是解方程（3）求待定系数，而是先构造一组基函数

1.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

1.2 牛顿（Newton）插值

在导出 Newton 公式前，先介绍公式表示中所需要用到的差商、差分的概念及性质。

1.2.1 差商

1.2.2 Newton 插值公式

1.2.3 差分

当节点等距时，即相邻两个节点之差（称为步长）为常数，Newton 插值公式的形式会更简单。此时关于节点间函数的平均变化率（差商）可用函数值之差（差分）来表示。

1.2.4 等距节点插值公式

1.3 分段线性插值

1.3.1 插值多项式的振荡

1.3.2 分段线性插值

1.3.3 用 Matlab 实现分段线性插值

1.4 埃尔米特(Hermite)插值

1.4.1 Hermite 插值多项式

如果对插值函数，不仅要求它在节点处与函数同值，而且要求它与函数有相同的一 阶、二阶甚至更高阶的导数值，这就是 Hermite 插值问题。本节主要讨论在节点处插值函数与函数的值及一阶导数值均相等的 Hermite 插值。

1.4.2 用 Matlab 实现 Hermite 插值

Matlab 中没有现成的 Hermite 插值函数，必须编写一个 M 文件实现插值。

设n 个节点的数据以数组 x0 （已知点的横坐标）， y0 （函数值）， y1（导数值）

输入（注意 Matlat 的数组下标从 1 开始），m 个插值点以数组 x 输入，输出数组 y 为m

个插值。编写一个名为 hermite.m 的 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

1.5 样条插值

许多工程技术中提出的计算问题对插值函数的光滑性有较高要求，如飞机的机翼外 形，内燃机的进、排气门的凸轮曲线，都要求曲线具有较高的光滑程度，不仅要连续， 而且要有连续的曲率，这就导致了样条插值的产生。

1.5.1 样条函数的概念

所谓样条（Spline）本来是工程设计中使用的一种绘图工具，它是富有弹性的细木条或细金属条。绘图员利用它把一些已知点连接成一条光滑曲线（称为样条曲线），并使连接点处有连续的曲率。

数学上将具有一定光滑性的分段多项式称为样条函数。具体地说，给定区间[a,b] 的一个分划

1.5.2 二次样条函数插值

1.5.3 三次样条函数插值

1.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。

例 1 机床加工

待加工零件的外形根据工艺要求由一组数据(x, y) 给出（在平面情况下），用程控

铣床加工时每一刀只能沿 x 方向和 y 方向走非常小的一步，这就需要从已知数据得到加

工所要求的步长很小的(x, y) 坐标。

表 1 中给出的 x, y 数据位于机翼断面的下轮廓线上，假设需要得到 x 坐标每改变

0.1 时的 y 坐标。试完成加工所需数据，画出曲线，并求出 x = 0 处的曲线斜率和

13 ≤ x ≤ 15 范围内 y 的最小值。

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 附近弯曲处），建议选用三次样条插值的结果。

1.6 B 样条函数插值方法

1.6.1 磨光函数

实际中的许多问题，往往是既要求近似函数（曲线或曲面）有足够的光滑性，又要 求与实际函数有相同的凹凸性，一般插值函数和样条函数都不具有这种性质。如果对于 一个特殊函数进行磨光处理生成磨光函数（多项式），则用磨光函数构造出样条函数作 为插值函数，既有足够的光滑性，而且也具有较好的保凹凸性，因此磨光函数在一维插 值（曲线）和二维插值（曲面）问题中有着广泛的应用。 由积分理论可知，对于可积函数通过积分会提高函数的光滑度，因此，我们可以利 用积分方法对函数进行磨光处理。

1.6.2 等距 B 样条函数

对于任意的函数 f (x) ，定义其步长为 1 的中心差分算子δ 如下：

1.6.3 一维等距 B 样条函数插值

等距 B 样条函数与通常的样条有如下的关系：

1.6.4 二维等距 B 样条函数插值

1.7 二维插值

前面讲述的都是一维插值，即节点为一维变量，插值函数是一元函数（曲线）。若

节点是二维的，插值函数就是二元函数，即曲面。如在某区域测量了若干点（节点）的

高程（节点值），为了画出较精确的等高线图，就要先插入更多的点（插值点），计算这

些点的高程（插值）。

1.7.1 插值节点为网格节点

解 编写程序如下：

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)

1.7.2 插值节点为散乱节点

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)

参考：http://t.csdnimg.cn/ETkcm