MATLAB实现：数学插值

关于插值问题，在MATLAB软件中有很多现成的函数可以调用，熟悉这些函数，可以游刃有余地解决许多实际问题。下面介绍几种常用的插值函数。

1. 一维插值

一维插值函数interp1，语法为

y=interpl(x0,y0,x,'method')

其中：method指定插值的方法，缺省时默认为线性插值。当x0为等距时可以用快速插值法，格式为'*nearest'、'*linear'、'*spline'、'*cubic'。

'nearest'最近项插值

'linear'线性插值

'spline'三次样条插值

'cubic'立方插值

所有的插值方法要求x0是单调的，并且不能超过x的范围。

2. 三次样条插值

在MATLAB软件中数据点称为断点。如果三次样条插值没有边界条件，最常用的方法就是采用非扭结（not-a-knot）条件。这个条件强迫第一个和第二个三次多项式的三阶导数相等。对最后一个和倒数第二个三次多项式也做同样的处理。

MATLAB软件中三次样条插值有如下几种函数：

其中：x0,y0是已知数据点；x是插值点；y是插值点的函数值。

对于三次样条插值，提倡使用函数csape。csape的返回值是pp形式，要求插值点的函数值，必须调用函数ppval。

'complete' 边界为一阶导数，一阶导数的值在valconds参数中给出，若忽略valconds参数，则按默认情况处理。

'not-a-knot' 非扭结条件

'periodic' 周期条件

'second' 边界为二阶导数，二阶导数的值在valconds参数中给出，若忽略 valconds参数，二阶导数的默认值为[0,0]。

'variational' 设置边界的二阶导数值为[0,0]。

对于一些特殊的边界条件，可以通过conds的一个1×2矩阵来表示，conds元素的取值为0、1、2。

conds(i)=j的含义是给定端点i的j阶导数，即conds的第一个元素表示左边界的条件，第二个元素表示右边界的条件，conds=[2,1]表示左边界是二阶导数，右边界是一阶导数，对应的值由valconds给出。详细情况请使用“帮助”文档doc csape。

例1 用interp1对正弦函数进行分段线性插值。

解 在MATLAB命令窗口中输入以下命令：

＞＞x=0: 2*pi;

＞＞y=sin(x);

＞＞x0=0: 0.5: 2*pi

＞＞y0=interp1(x, y, x0);

＞＞plot(x, y, 's', x0, y0)

注：例1中用默认的分段线性插值对已知的7个正弦函数的数据点进行插值，用plot画出插值结果。从图5.1.1中可以看出，分段线性就是联结两个邻近的已知点的线性函数插值计算该区间内插值点上的函数值。

图5.1.1

例2 根据表5.1.1的数据点求出其拉格朗日插值多项式，并计算当x=1.6时y的值。

表5.1.1　数据点

解 在MATLAB命令窗口中输入以下命令：

基本函数可以从工具箱中调用，复杂一点的函数则需要自己编写，如这里的Language函数。详见本章附录1。

例3 根据表5.1.2的数据点求出其Hermite插值多项式，并计算当x=1.144时y的值。

表5.1.2　数据点

解 在MATLAB命令窗口中输入以下命令：(www.xing528.com)

注意，这里的Hermite函数也需要自己编写（详见本章附录2）。

例4 从1点到12点的11小时内，每隔1小时测量一次温度，测得的温度的数值依次为5，8，9，15，25，29，31，30，22，25，27，24。试估计每隔1/10小时的温度值。

解 在MATLAB命令窗口中输入以下命令：

＞＞hours=1: 12;

＞＞temps=[5 8 9 15 25 29 31 30 22 25 27 24];

＞＞h=1: 0.1: 12;

＞＞t=interp1(hours, temps, h, 'spline');

＞＞plot(hours, temps, '+', h, t, hours, temps, 'r: ')

＞＞xlabel('Hour'), ylabel('Degrees Celsius')

计算结果略（见图5.1.2）。

图5.1.2

3. 二维插值

前面讲述的都是一维插值，即节点为一维变量，插值函数是一元函数（曲线）。若节点是二维的，插值函数就是二元函数，即曲面。如在某区域测量了若干节点的节点值，为了画出较精确的曲面图，就要先插入更多的点，计算这些点的插值。

（1）插值节点为网格节点。

已知m×n个节点: (xi ,y j,z ij )(i=1,2,…,m; j=1,2,…,n)，且x1＜…＜x m;y1…＜yn 。求点(x,y)处的插值z。

MATLAB中有一些计算二维插值的命令。如：

z=interp2( x0 ,y0,z0 ,x,y,'method')

其中：x0，y0分别为m维和n维向量，表示节点；z0为n×m矩阵，表示节点值；x，y为一维数组，表示插值点，x与y应是方向不同的向量，即一个是行向量，另一个是列向量；z为矩阵，它的行数为y的维数，列数为x的维数，表示得到的插值；′Method′的用法同上面的一维插值。

如果是三次样条插值，可以使用如下命令：

其中：x0，y0分别为m维和n维向量；z0为m×n矩阵；z为矩阵，它的行数为x的维数，列数为y的维数，表示得到的插值。具体使用方法同一维插值。

例5 测得平板表面3×5网格点处的温度分别为

解 在MATLAB命令窗口中输入以下命令：

（2）插值节点为散乱节点。

已知n个节点(xi ,yi ,zi )，i=1,2,…n ，求点（x,y）处的插值z。

对上述问题，可以构造一个二元函数z=f (x, y)，通过全部已知节点，用f(x,y)计算插值。MATLAB中提供了插值函数griddata，其格式为

其中：x，y，z均为n维向量，指明所给数据点的横坐标、纵坐标和竖坐标；向量XI，YI是给定的网格点的横坐标和纵坐标；返回值ZI为网格( XI,YI )处的函数值。XI与YI应是方向不同的向量，即一个是行向量，另一个是列向量。

例6 在某海域测得一些点（x, y）处的水深z由表5.1.3给出，船的吃水深度为5英尺，在矩形区域( 75,200 )×(-50,150)的哪些地方，船要避免进入。

表5.1.3　数据点

%程序一：插值并做出海底曲面图（见图5.1.3）

图5.1.3　海底曲线图

z1(z1＞=5)=nan; %将水深大于5的置为nan

meshc(x1, y1, z1)

%程序二：插值并做出水深小于5的海域范围。