MATLAB插值程序.docx
《MATLAB插值程序.docx》由会员分享，可在线阅读，更多相关《MATLAB插值程序.docx（10页珍藏版）》请在冰豆网上搜索。
MATLAB插值程序
插值程序
例1：
采用拉格朗日多项式插值：
选取不同插值节点个数n+1，其中n为插值多项式的次数，当n分别取2,4,6,8,10时，绘出插值结果图形.
1、Large1:
functiony=lagr1（x0,y0,x）
n=length（x0）;m=length（x）;
fori=1:
m
z=x（i）;
s=0.0;
fork=1:
n
p=1.0;
forj=1:
n
ifj~=k
p=p*（z-x0（j））/（x0（k）-x0（j））;
end
end
s=p*y0（k）+s;
end
y（i）=s;
end
2、lch:
m=101;
x=-5:
10/（m-1）:
5;
y=1./（1+x.^2）;z=0*x;
plot（x,z,'r',x,y,'LineWidth',1.5）,
gtext（'y=1/（1+x^2）'）,pause
n=3;
x0=-5:
10/（n-1）:
5;
y0=1./（1+x0.^2）;
y1=lagr1（x0,y0,x）;
holdon,plot（x,y1,'b'）,gtext（'n=2'）,pause,
holdoff
n=5;
x0=-5:
10/（n-1）:
5;
y0=1./（1+x0.^2）;
y2=lagr1（x0,y0,x）;
holdon,plot（x,y2,'b:
'）,gtext（'n=4'）,pause,
holdoff
n=7;
x0=-5:
10/（n-1）:
5;
y0=1./（1+x0.^2）;
y3=lagr1（x0,y0,x）;holdon,
plot（x,y3,'r'）,gtext（'n=6'）,pause,
holdoff
n=9;
x0=-5:
10/（n-1）:
5;
y0=1./（1+x0.^2）;
y4=lagr1（x0,y0,x）;holdon,
plot（x,y4,'r:
'）,gtext（'n=8'）,pause,
holdoff
n=11;
x0=-5:
10/（n-1）:
5;
y0=1./（1+x0.^2）;
y5=lagr1（x0,y0,x）;holdon,
plot（x,y5,'m'）,gtext（'n=10'）
例2：
用分段线性插值法求插值,并观察插值误差.
（1）xch11:
x=linspace（-6,6,100）;
y=1./（x.^2+1）;
x1=linspace（-6,6,5）;
y1=1./（x1.^2+1）;
plot（x,y,x1,y1,x1,y1,'o','LineWidth',1.5）,
gtext（'n=4'）,
（2）xch12
x=linspace（-6,6,100）;
y=1./（x.^2+1）;
x1=linspace（-6,6,11）;
y1=1./（x1.^2+1）;
plot（x,y,x1,y1,x1,y1,'o','LineWidth',1.5）,
gtext（'n=10'）,
（3）xch13
x=linspace（-6,6,100）;
y=1./（x.^2+1）;
x1=linspace（-6,6,21）;
y1=1./（x1.^2+1）;
plot（x,y,x1,y1,x1,y1,'o','LineWidth',1.5）,
gtext（'n=20'）,
（4）xch14
x=linspace（-6,6,100）;
y=1./（x.^2+1）;
x1=linspace（-6,6,41）;
y1=1./（x1.^2+1）;
plot（x,y,x1,y1,x1,y1,'o','LineWidth',1.5）,
gtext（'n=40'）,
例2:
用三次样条插值选取11个基点计算插值
ych：
x0=linspace（-5,5,11）;
y0=1./（1+x0.^2）;
x=linspace（-5,5,100）;
y=interp1（x0,y0,x,'spline'）;
x1=linspace（-5,5,100）;
y1=1./（1+x1.^2）;
plot（x1,y1,'k',x0,y0,'+',x,y,'r'）;
例3：
在1-12的11小时内，每隔1小时测量一次温度，测得的温度依次为：
5，8，9，15，25，29，31，30，22，25，27，24。
试估计每隔1/10小时的温度值。
temp:
hours=1:
12;
temps=[589152529313022252724];
h=1:
0.1:
12;
t=interp1（hours,temps,h,'spline'）;
plot（hours,temps,'+',h,t,hours,temps,'r:
'）
xlabel（'Hour'）,ylabel（'DegreesCelsius'）
例4：
已知飞机下轮廓线上数据如下，求x每改变0.1时的y值。
plane:
x0=[035791112131415];
y0=[01.21.72.02.12.01.81.21.01.6];
x=0:
0.1:
15;
y1=lagr1（x0,y0,x）;
y2=interp1（x0,y0,x）;
y3=interp1（x0,y0,x,'spline'）;
subplot（3,1,1）
plot（x0,y0,'k+',x,y1,'r'）
grid
title（'lagrange'）
subplot（3,1,2）
plot（x0,y0,'k+',x,y2,'r'）
grid
title（'piecewiselinear'）
subplot（3,1,3）
plot（x0,y0,'k+',x,y3,'r'）
grid
title（'spline'）
例5：
测得平板表面3*5网格点处的温度分别为：
828180828479636165818484828586试作出平板表面的温度分布曲面z=f（x,y）的图形。
wendu:
x=1:
5;
y=1:
3;
temps=[8281808284;7963616581;8484828586];
mesh（x,y,temps）
pause
xi=1:
0.2:
5;
yi=1:
0.2:
3;
zi=interp2（x,y,temps,xi',yi,'cubic'）;
figure
（2）
mesh（xi,yi,zi）
例6：
山区地貌：
在某山区测得一些地点的高程如下表。
平面区域为
1200<=x<=4000,1200<=y<=3600）
试作出该山区的地貌图和等高线图，并对几种插值方法进行比较。
X
Y
1200
1600
2000
2400
2800
3200
3600
4000
1200
1130
1250
1280
1230
1040
900
500
700
1600
1320
1450
1420
1400
1300
700
900
850
2000
1390
1500
1500
1400
900
1100
1060
950
2400
1500
1200
1100
1350
1450
1200
1150
1010
2800
1500
1200
1100
1550
1600
1550
1380
1070
3200
1500
1550
1600
1550
1600
1600
1600
1550
3600
1480
1500
1550
1510
1430
1300
1200
980
通过此例对最近邻点插值、双线性插值方法和双三次插值方法的插值效果进行比较。
moutain:
x=0:
400:
5600;
y=0:
400:
4800;
z=[370470550600670690670620580450400300100150250;...
510620730800850870850780720650500200300350320;...
650760880970102010501020830900700300500550480350;...
740880108011301250128012301040900500700780750650550;...
830980118013201450142014001300700900850840380780750;...
88010601230139015001500140090011001060950870900930950;...
9101090127015001200110013501450120011501010880100010501100;...
9501190137015001200110015501600155013801070900105011501200;...
143014301460150015501600155016001600160015501500150015501550;...
1420143014501480150015501510143013001200980850750550500;...
138014101430145014701320128012001080940780620460370350;...
13701390141014301440114011101050950820690540380300210;...
13501370139014001410960940880800690570430290210150];
figure
（1）;
meshz（x,y,z）
xlabel（'X'）,ylabel（'Y'）,zlabel（'Z'）
xi=0:
50:
5600;
yi=0:
50:
4800;
figure
（2）
z1i=interp2（x,y,z,xi,yi','nearest'）;
surfc（xi,yi,z1i）
xlabel（'X'）,ylabel（'Y'）,zlabel（'Z'）
figure（3）
z2i=interp2（x,y,z,xi,yi'）;
surfc（xi,yi,z2i）
xlabel（'X'）,ylabel（'Y'）,zlabel（'Z'）
figure（4）
z3i=interp2（x,y,z,xi,yi','cubic'）;
surfc（xi,yi,z3i）
xlabel（'X'）,ylabel（'Y'）,zlabel（'Z'）
figure（5）
subplot（1,3,1）,contour（xi,yi,z1i,10,'r'）;
subplot（1,3,2）,contour（xi,yi,z2i,10,'r'）;
subplot（1,3,3）,contour（xi,yi,z3i,10,'r'）;
实验作业
山区地貌：
在某山区测得一些地点的高程如下表：
（平面区域:
1200<=x<=4000,1200<=y<=3600），
试作出该山区的地貌图和等高线图，并对几种插值方法进行比较。
升级会员 