matlab插值与拟合(命令与示例).doc

matlab插值与拟合(命令与示例)

———————————————————————————————— 作者：
———————————————————————————————— 日期：

21

个人收集整理 勿做商业用途
目录
【一维插值】interp1 1
yi = interp1(x,y,xi,method) 1
例1 1
例2 2
【二维插值】interp2 3
ZI = interp2(X，Y,Z,XI,YI,method） 3
插值方式比较示例 3
例3 3
例4 3
【三角测量和分散数据插值】 3
【数据拟合】 3
例5 3
例6 3

【一维插值】interp1
yi = interp1(x,y，xi,method)
例1
在1—12的11小时内，每隔1小时测量一次温度，测得的温度依次为：5，8,9,15，25，29，31，30,22，25，27,24。试估计每隔1/10小时的温度值。

建立M文件temp。m
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,’kp’，h,t，’b'）;

例2
已知飞机下轮廓线上数据如下,求x每改变0。1时的y值.

建立M文件plane.m
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=interp1(x0,y0,x,’nearest'）;
y2=interp1(x0,y0，x）;
y3=interp1（x0,y0，x，’spline')；
plot（x0,y0，'kp’，x，y1,’r')

plot(x0，y0，’kp'，x,y2,’r'）

plot（x0，y0，’kp',x，y3，’r')

【二维插值】interp2
ZI = interp2（X,Y，Z，XI,YI,method）
插值方式比较示例
l 用较大间隔产生peaks函数数据点
[x,y] = meshgrid(—3：1:3）;
z = peaks(x，y);
surf（x,y，z）

l 产生一个较好的网格
[xi,yi］ = meshgrid(-3：0。25：3);
l 利用最近邻方式插值
zi1 = interp2（x,y，z，xi，yi，’nearest’）;surf（xi,yi,zi1)

l 双线性插值方式
zi2 = interp2（x,y，z,xi,yi，'bilinear')；surf（xi，yi,zi2）

l 双立方插值方式
zi3 = interp2(x，y,z,xi，yi，'bicubic’）；surf(xi，yi，zi3）

l 不同插值方式构造的等高线图对比
contour(xi，yi,zi1）

contour（xi,yi,zi2）

contour（xi,yi,zi3）

例3
测得平板表面3*5网格点处的温度分别为:
82 81 80 82 84
79 63 61 65 81
84 84 82 85 86
试作出平板表面的温度分布曲面z=f(x,y)的图形.
建立M文件wendu.m
xi=1：0。2：5;
yi=1：0.2：3;
zi=interp2(x,y,temps,xi'，yi，'cubic'）；
mesh（xi,yi,zi)；

例4
某山区测得一些地点的高度如下表所示，平面区域为，试作出该山区的地貌图和等高线图。比较几种插值方法。

建立M文件moutain。m
x=0：400：5600；
y=0：400：4800;
z=[370 470 550 600 670 690 670 620 580 450 400 300 100 150 250;.。。
510 620 730 800 850 870 850 780 720 650 500 200 300 350 320；。。。
650 760 880 970 1020 1050 1020 830 900 700 300 500 550 480 350;。。.
740 880 1080 1130 1250 1280 1230 1040 900 500 700 780 750 650 550;。.。
830 980 1180 1320 1450 1420 1400 1300 700 900 850 840 380 780 750；。..
880 1060 1230 1390 1500 1500 1400 900 1100 1060 950 870 900 930 950;。..
910 1090 1270 1500 1200 1100 1350 1450 1200 1150 1010 880 1000 1050 1100;。.。
950 1190 1370 1500 1200 1100 1550 1600 1550 1380 1070 900 1050 1150 1200;...
1430 1430 1460 1500 1550 1600 1550 1600 1600 1600 1550 1500 1500 1550 1550；.。。
1420 1430 1450 1480 1500 1550 1510 1430 1300 1200 980 850 750 550 500;..。
1380 1410 1430 1450 1470 1320 1280 1200 1080 940 780 620 460 370 350;...
1370 1390 1410 1430 1440 1140 1110 1050 950 820 690 540 380 300 210;..。
1350 1370 1390 1400 1410 960 940 880 800 690 570 430 290 210 150］;

figure（1)；
meshz（x，y,z)

xi=0：50：5600；
yi=0：50:4800

figure（2)
z1i=interp2（x，y，z,xi，yi'，’nearest')；
surfc（xi，yi，z1i)

figure（3）
z2i=interp2(x，y，z,xi，yi');
surfc（xi，yi，z2i)

figure（4）
z3i=interp2（x，y,z,xi，yi'，'cubic’)；
surfc(xi,yi,z3i）

figure（5)
subplot（1,2，1）,contour（xi，yi，z2i，10）;
subplot(1,2,2),contour（xi，yi，z3i，10）；

【三角测量和分散数据插值】
凸包（Convex Hulls）
load seamount
plot（x，y，’。’,'markersize'，10）
k = convhull(x,y);
hold on， plot（x（k)，y（k),’—r')， hold off
grid on

德洛涅三角（Delaunay Triangulation）
load seamount
plot(x，y,’。','markersize'，12）
xlabel(’Longitude’)， ylabel（’Latitude’）
grid on

tri = delaunay（x,y)；
hold on, triplot(tri,x,y）， hold off

figure
hidden on
trimesh（tri，x,y,z）
grid on
xlabel（’Longitude’）; ylabel(’Latitude'); zlabel（'Depth in Feet'）

figure
［xi,yi］ = meshgrid(210。8:。01:211。8,-48.5：.01：—47.9）;
zi = griddata（x,y,z,xi,yi，’cubic’);
［c,h］ = contour(xi，yi,zi，'b—'）;
clabel（c，h)
xlabel（'Longitude’）， ylabel（'Latitude’)

火龙尼图形（Voronoi Diagrams）
load seamount
voronoi（x,y)
grid on
xlabel（'Longitude')， ylabel('Latitude')

【数据拟合】
例5
对下面一组数据作二次多项式拟合

x=［0.1 0.2 0。4:.1:1］;
y=［1.978 3。28 6。16 7。34 7.66 9.58 9.48 9。30 11。2］；
A=polyfit（x,y,2）；
z=polyval(A,x)；
plot（x，y，’k+'，x，z,’r'）

例6
用下面一组数据拟合中的参数a,b，k.
方法1:用lsqcurvefit
建立M文件curvefun1.m
function f=curvefun1(x，tdata)
f=x(1)+x（2)＊exp(-0.02＊x(3）*tdata)
％其中 x(1）=a； x（2）=b；x(3)=k；

输入命令：
tdata=100：100：1000；
cdata=1e-03*［4。54，4。99,5.35,5。65，5.90,6.10,6。26，6。39,6。50，6。59］；
x0=［0.2,0。05，0.05];
x=lsqcurvefit （'curvefun1’,x0,tdata，cdata）

x =
0.0063 —0.0034 0.2542

方法2：用lsqnonlin
建立M文件curvefun2.m
function f=curvefun2(x)
tdata=100:100：1000；
cdata=1e-03*［4.54,4。99，5。35，5。65,5.90,6.10,6。26，6。39，6.50,6.59］;
f=x（1)+x(2）*exp（-0。02＊x(3)*tdata)- cdata
输入命令：
x0=［0。2，0.05，0.05］；
x=lsqnonlin（’curvefun2’，x0)
x =
0.0063 -0。0034 0。2542