1、第七章 一维、二维插值方法实验,一维插值方法二维插值方法,l一维插值方法,1.引例1 已经测得在北纬32.3 海洋不同深度处的温度如下表：表7.1.1 根据这些数据，我们希望能合理地估计出其它深度（如500米、600米、1000米）处的水温。解决这个问题，可以通过构造一个与给定数据相适应的函数来解决，这是一个被称为插值的问题。,2插值问题的基本提法 对于给定的函数表,其中f (x) 在区间 a, b 上连续，x0，x1，xn为 a, b 上n 1个互不相同的点，要求在一个性质优良、便于计算的函数类 P(x) 中，选出一个使,P(xi ) = yi，i = 0, 1, , n (7.1.1),成
2、立的函数P(x) 作为 f (x) 的近似，这就是最基本的插值问题（见图7.1.1）。,为便于叙述，通常称区间 a, b 为插值区间，称点x0，x1，xn为插值节点，称函数类 P(x) 为插值函数类，称式 (7.1.1) 为插值条件，称函数 P(x) 为插值函数，称f (x) 为被插函数。求插值函数 P(x) 的方法称为插值法。,l一维数据插值的Matlab实现,2例： 7.1.1引例1的求解解法一 直接借助于Matlab命令，分别采用最临近插值方法、线性插值方法、三次插值方法和三次样条插值方法进行插值。Matlab程序如下：,x=466,714,950,1422,1634;y=7.04,4.
3、28,3.40,2.54,2.13;x1=linspace(400,1660,100);y1=interp1(x,y,x1,nearest);y2=interp1(x,y,x1,linear);y3=interp1(x,y,x1,cubic);y4=interp1(x,y,x1,spline);subplot(2,2,1),plot(x,y,o,x1,y1),title(最邻近插值曲线)subplot(2,2,2),plot(x,y,o,x1,y2),title(线性插值曲线)subplot(2,2,3),plot(x,y,o,x1,y3),title(三次插值曲线)subplot(2,2,4
4、),plot(x,y,o,x1,y4),title(三次样条插值曲线)x2=500:100:1600w1=interp1(x,y,x2,nearest)w2=interp1(x,y,x2,linear)w3=interp1(x,y,x2,cubic)w4=interp1(x,y,x2,spline),要想求得任一深度值xi处的水温，只需在上述程序运行完毕之后，继续运行如下命令即可：,xi=500:100:1600yi1=interp1(x,y,xi,nearest)yi2=interp1(x,y,xi,linear)yi3=interp1(x,y,xi,cubic)yi4=interp1(x,
5、y,xi,spline),下面是运行上述程序后得到的四种不同插值方法在500米到1600米不同深度处的水温值（见表7.2.1），以及四种不同的插值函数的曲线图（见图7.2.2）。从图中可以看出，三次插值和三次样条插值得到的插值曲线的光滑性更要好些。,图7.2.2 四种不同插值函数的曲线图,如果要想求得任一深度值xi处的水温，只需在上述程序运行完毕之后，继续运行如下命令即可：xi=500:100:1600yi1=interp1(x,y,xi,nearest)yi2=interp1(x,y,xi,linear)yi3=interp1(x,y,xi,cubic)yi4=interp1(x,y,xi,
6、spline),表7.2.1,l二维插值方法,二维规则插值二维散乱插值,l二维数据插值的Matlab实现,解 下面分别用最邻近插值、双线性插值、双三次插值和双三次样条插值，给出不同月份按纬度变化的气旋值（插值结果），并作出可视化图形如下。,最邻近插值的Matlab程序为： x=1:12; y=5:10:85; z=2.4, 18.7, 20.8, 22.1, 37.3, 48.2, 25.6, 5.3, 0.3 1.6, 21.4, 18.5, 20.1, 28.8, 36.6, 24.2, 5.3, 0 2.4, 16.2, 18.2, 20.5, 27.8, 35.5, 25.5, 5.4
7、, 0 3.2, 9.2, 16.6, 25.1, 37.2, 40, 24.6, 4.9, 0.3 1.0, 2.8, 12.9, 29.2, 40.3, 37.6, 21.1, 4.9, 0 0.5, 1.7, 10.1, 32.6, 41.7, 35.4, 22.2, 7.1, 0 0.4, 1.4, 8.3, 33.0, 46.2, 35, 20.2, 5.3, 0.1 0.2, 2.4, 11.2, 31.0, 39.9, 34.7, 21.2, 7.3, 0.2 0.5, 5.8, 12.5, 28.6, 25.9, 35.7, 22.6, 7, 0.3 0.8, 9.2, 21.
8、1, 32.0, 40.3, 39.5, 28.5, 8.6, 0 2.4, 10.3, 23.9, 28.1, 38.2, 40, 25.3, 6.3, 0.1 3.6, 16, 25.5, 25.6, 43.4, 41.9, 24.3, 6.6, 0.3; xi,yi=meshgrid(1:12,5:1:85); zi=interp2(x,y,z,xi,yi,nearest); figure surf(xi,yi,zi) xlabel(月份),ylabel(纬度),zlabel(气旋), axis(0 12 0 90 0 50) title(南半球气旋可视化图形),双线性插值、双三次插值、
9、双三次样条插值的Matlab程序为：分别将最邻近线性插值程序中的zi=interp2(x,y,z,xi,yi,nearest)改写为zi=interp2(x,y,z,xi,yi,linear)zi=interp2(x,y,z,xi,yi,cubic)zi=interp2(x,y,z,xi,yi,spline),例2 水道测量数据（AMCM86A题）在某海域测得一些点（x, y）处的水深z（单位：英尺）由表7.3.2给出，水深数据是在低潮时测得的。船的吃水深度为5英尺，问在矩形区域 (75，200) (50，150) 里的哪些地方船要避免进入。,解 （1）假设 由题目给出的信息是很少的，除了14
10、个位置的水深之外一无所知。显然，题目要求我们找出水深不到5英尺的区域。为了讨论方便，下面三个假设是合理的： 所给数据是精确的； 讨论区域的海底曲面是光滑的，更确切地说，可以认为曲面的一阶、二阶导数是连续的。因为我们可以认为讨论区域为浅水海域，由于长期的海水水流作用，形成的是以砾石或沙为主要组成部分的海底，不存在珊瑚礁、水底峡谷、山脊等不可意料的突变地形。 水深是一个按区域来划分的变量，在某个位置的水深与其周围区域的水深是相互依赖的，但这种依赖作用随距离的增大而减小。就我们讨论的问题来说，每一个给定数据点影响周围的每一个未知点，一个给定数据点离未知点越近，作用就越大。 （2）问题分析 根据假设，
11、海底曲面是连续光滑的，不存在珊瑚礁、水底峡谷、山脊等不可意料的突变地形，因而很自然的想法就是用某种光滑的拟合曲面去逼近已知的14个数据点或以14个已知的数据点为基础，利用二维插值补充一些点的水深，以求得水深不超过5米的区域。在此，我们采用二维插值方法，应用Matlab程序，作出矩形区域 (75，200) (50，150) 范围内的海底地形图、水深不超过5米的危险区域的平面图以及水深不超过5米的危险区域的海底地貌图，并求出水深不超过5米的危险海域范围。 （3）问题求解采用改进的Shepard方法，利用Matlab软件作出作出矩形区域 (75，200) (50，150) 范围内的海底地形图、水深不
12、超过5米的危险区域的平面图以及水深不超过5米的危险区域的海底地貌图（见图7.3.3），并求出水深不超过5米的危险海域范围为：115，2003，119。,（4）求解的Matlab程序如下：clear;x=129,140,108.5,88,185.5,195,105.5,157.5,107.5,77,81,162,162,117.5;y=7.5,141.5,28,147,22.5,137.5,85.5,-6.5,-81,3,56.5,-66.5,84,-38.5;subplot(2,2,1),plot(x,y,+),title(测量点分布图);z=-4,-8,-6,-8,-6,-8,-8,-9,-
13、9,-8,-8,-9,-4,-9;fnodes,minnq,rnw,rnq,ifail=e01sef(x,y,z);nx=100;px=linspace(75,200,nx);,ny=200;py=linspace(-50,150,ny);for i=1:ny for j=1:nx pf(i,j),ifail=e01sff(x,y,z,rnw,fnodes,px(j),py(i); endendsubplot(2,2,2),meshz(px,py,pf+5),title(75，200)x(-50，150) 范围内的海底地形图);a,b=find(pf=-5);amin=min(a);amax=
14、max(a);bmin=min(b);bmax=max(b);xmin=75+(200-75)/100)*bminxmax=75+(200-75)/100)*bmaxymin=-50+(150+50)/200)*aminymax=-50+(150+50)/200)*amax,for k=1:length(b) i0(k)=75+(200-75)/100)*b(k);endfor k=1:length(a) j0(k)=-50+(150+50)/200)*a(k);endsubplot(2,2,3),plot(i0,j0,+),title(水深不超过5米的危险区域的平面图);i1,j1=find
15、(pf-5);for k=1:length(i1) pf(i1(k),j1(k)=-5;endsubplot(2,2,4),meshc(px,py,pf),title(水深不超过5米的危险区域的海底地貌图); rotate3d,X=129 140 103.5 88 185.5 195 105.5 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; plot(X,Y,*); Z=-4 -8 -6 -8 -6 -8 -8 -9 -9 -8 -8 -9 -4-9; %a=linspace(75,200,100); %b=linspace(-50,150,200); x,y=meshgrid(75:.5:200,-50:.5:150); z=griddata(X,Y,Z,x,y,cubic); figure(2); meshz(x,y,z+5); figure(3); contour(x,y,z,-5,5);figure(4); meshz(x,y,z);rotate3d 同样可以作出海底地貌图，只是效果不是太好,
