计算插值的软件很多,这里我们只介绍如何用MATLAB做一维插值和高维插值.
4.4.1 一维插值
MATLAB中的插值函数为interp1(),其调用格式为
yi=interp1(x,y,xi, 'method')
其中x,y为观测数据点,xi为插值(自变量)向量,yi为xi的插值结果(函数值).'method '表示采用的插值方法.MATLAB提供的插值方法有几种: 'nearest' 最邻近插值; 'linear '线性插值; 'spline' 三次样条插值; 'cubic ' 立方插值.缺省时表示线性插值.
注意:所有的插值方法都要求x是单调的,并且xi不能够超过x的范围.
例4.2 在一天24小时内,从零点开始每间隔2小时测得的环境温度为(摄氏度)
12,9,9,10,18,24,28,27,25,20,18,15,13
推测在每一秒时的温度.并利用不同的插值方法描绘温度曲线
键入：
x=0:2:24;
y=[12 9 9 10 18 24 28 27 25 20 18 15 13];
xi=0:1/3600:24;
yi=interp1(x,y,xi,'nearest');
hold on
plot(xi,yi,'r');
yi=interp1(x,y,xi,'linear');
plot(xi,yi,'g');
yi=interp1(x,y,xi,'spline');
plot(xi,yi,'b');
yi=interp1(x,y,xi,'cubic');
plot(xi,yi,'y');
还有其他的插值函数,如 interplq ,interpft, spline, intep2, interp3, interpN.
4.4.2 高维插值
N维插值函数 interpN( )
其中N可以为2,3,…,如N=2为二维插值,调用格式为
zi=interp2(x,y,z,Xi,Yi,’method’)
其中x,y为横纵坐标上的坐标点,{(x,y)}=mashgrid(x,y)生成平面网格点,z为观测到的在网格点上的二元函数值.{(x,y,z)}构成空间插值节点.引入两个向量xi,yi.xi为横坐标上的插值点,yi为纵坐标上的插值点.便可给出
[Xi,Yi]=meshgrid(xi,yi)
zi为新的或者是加细了的网格点上产生的插值结果(函数值).
'method' 表示采用的插值方法.`nearest` 最邻近插值,`linear`线性插值,`cubic`双三次插值.缺省时表示线性插值.所有的插值方法都要求x和y是单调的网格,x和y可以是等距的也可以是不等距的.
例如,产生一个山顶函数peaks曲面.
1)产生peaks的粗糙近似山顶曲面
[x,y,z]=peaks(10);
hold on
mesh(x,y,z)
2)通过插值作出更加精细的山顶曲面
xi=-3:.1:3;
yi=xi;
[Xi,Yi]=meshgrid(xi,yi);
Zi=interp2(x,y,z,Xi,Yi,`cubic`);
mesh(Xi,Yi,Zi)
表4.1 是气象学家测量得到的气象资料,它们分别表示在南半球地区按不同纬度、不同月份的平均气旋数字.根据这些数据,绘制出气旋分布曲面图形.
解 下面用二维(interp2)三次(cubic)插值方法,可以得到不同月份按纬度变化的气旋值(插值结果),然后再作出其可视化图形,如图4.1MATLAB程序如下:
x=1:12;
y=5:10:85;
z1=[2.4 1.6 2.4 3.2 1.0 0.5 0.4 0.2 0.5 0.8 2.4 3.6];
z2=[18.7 21.4 16.2 9.2 2.8 1.7 1.4 2.4 5.8 9.2 10.3 16];
z3=[20.8 18.5 18.2 16.5 12.9 10.1 8.3 11.2 12.5 21.1 23.9 25.5];
z4=[22.1 20.1 20.5 25.1 29.2 32.6 33.0 31.0 28.6 32.0 28.1 25.6];
z5=[37.3 28.8 27.8 37.2 40.3 41.7 46.2 39.9 35.9 40.3 38.2 43.4];
z6=[48.2 36.6 35.5 40 37.6 35.4 35 34.7 35.7 39.5 40 41.9];
z7=[25.6 24.2 25.5 24.6 21.1 22.2 20.2 21.2 22.6 28.5 25.3 24.3];
z8=[5.3 5.3 5.4 4.9 4.9 7.1 5.3 7.3 7 8.6 6.3 6.6 ];
z9=[0.3 0 0 0.3 0 0 0.1 0.2 0.3 0 0.1 0.3];
z=[z1;z2;z3;z4;z5;z6;z7;z8;z9];
[xi,yi]=meshgrid(1:12,5:1:85);
zi=interp2(x,y,z,xi,yi,'cudic');
mesh(xi,yi,zi)
xlabel('月份'),ylabel('纬度'),zlabel('气旋'),
title('南半球气旋可视化图形')
还有两个二维插值函数e01sef和e01sff,它们分别被用于求散点数据的插值函数和插值函数值,通常是两者配合使用,其调用格式为：
[fnodes,a,rnw,b,c]=e01sef(x,y,z);
[pf(i,j),ifail]=e01sff(x,y,z,rnw,fnodes,px(j),py(i);
其中x,y,z为插值节点,均为向量；px(j),py(i)为被插值点;pf(I,j)为被插值.
表4.1 南半球地区按不同纬度 不同月份的平均气旋数据
0-10 10-20 20-30 30-40 40-50 50-60 60-70 70-80 80-90
1月 2.4 18.7 20.8 22.1 37.3 48.2 25.6 5.3 0.3
2月 1.6 21.4 18.5 20.1 28.8 36.6 24.2 5.3 0
3月 2.4 16.2 18.2 20.5 27.8 35.5 25.5 5.4 0
4月 3.2 9.2 16.6 25.1 37.2 40 24.6 4.9 0.3
5月 1.0 2.8 12.9 29.2 40.3 37.6 21.1 4.9 0
6月 0.5 1.7 10.1 32.6 41.7 35.4 22.2 7.1 0
7月 0.4 1.4 8.3 33.0 46.2 35 20.2 5.3 0.1
8月 0.2 2.4 11.2 31.0 39.9 34.7 21.2 7.3 0.2
9月 0.5 5.8 12.5 28.6 35.9 35.7 22.6 7 0.3
10月 0.8 9.2 21.1 32.0 40.3 39.5 28.5 8.6 0
11月 2.4 10.3 23.9 28.1 38.2 40 25.3 6.3 0.1
12月 3.6 16 25.5 25.6 43.4 41.9 24.3 6.6 0.3
点(px(j),py(j)处的插值结果,其他输出参数涉及插值算法,可以不用了解.有兴趣的读者可以查看MATLAB的帮助文件.e01sef的输出fnodes和rnw为确定插值函数的参数,它们是e01sff需要的输入参数,因此两函数需配合使用.这两个函数将在§4.6中使用
