《六Matlab插值(精品).ppt》由会员分享,可在线阅读,更多相关《六Matlab插值(精品).ppt(59页珍藏版)》请在淘文阁 - 分享文档赚钱的网站上搜索。
1、MATLAB插值插值2/10/2023天 津 科 技 大 学 数 学 系 谢中华办公室:泰达4号楼211室 电话:60600830E-mail:MATLAB统计分析与应用:统计分析与应用:40个案例分析个案例分析http:/http:/ 一维插值一维插值 二维插值二维插值 高维插值高维插值MATLAB插值插值2/10/2023第一节第一节 一维插值一维插值MATLAB插值插值2/10/2023 一维插值问题的一维插值问题的数学描述数学描述为为:已知某一函数已知某一函数 y=g(x)(g(x)的解析表达式可能十分复杂的解析表达式可能十分复杂,也可以是未知的)在区间也可以是未知的)在区间a,b上上
2、n+1个互异点个互异点 xj 处的函数值处的函数值 y j,j=0,1,n,还知道还知道g(x)在在a,b上有若干阶导数,如何上有若干阶导数,如何求出求出g(x)在在a,b上任一点上任一点x的近似值的近似值。一、一、数学描述及基本思想数学描述及基本思想数学描述数学描述MATLAB插值插值2/10/2023 一维插值方法的一维插值方法的基本思想基本思想是:根据是:根据g(x)在区间在区间a,b上上n+1个互异点个互异点 x j(称为节点)的函数值(称为节点)的函数值 y j,j=0,1,n,求,求一个足够光滑、简单便于计算的函数一个足够光滑、简单便于计算的函数f(x)(称为插值函数)(称为插值函
3、数)作为作为g(x)的近似表达式,使得的近似表达式,使得f(x j)=y j,j=0,1,n。(1)然后计算然后计算f(x)在区间在区间a,b(称为插值区间称为插值区间)上点上点x(称为插值(称为插值点)的值作为原函数点)的值作为原函数g(x)(称为被插函数)在此点的近似(称为被插函数)在此点的近似值。求插值函数值。求插值函数f(x)的方法称为插值方法,(的方法称为插值方法,(1)式称为插)式称为插值条件值条件。代数多项式比较简单,常用多项式作为插值函数代数多项式比较简单,常用多项式作为插值函数。2基本思想基本思想MATLAB插值插值2/10/2023 假设假设f(x)是一个满足插值条件(是一
4、个满足插值条件(1)的次数不超过)的次数不超过n的代数的代数多项式,即多项式,即 f(x)=a0+a1x+a n x n为满足(为满足(1)的插值函数的插值函数,则,则f(x)的的n+1个待定系数个待定系数 a 0,a 1,a n 满足满足二、多项式插值二、多项式插值MATLAB插值插值2/10/2023缺点:当缺点:当 n 比较大时,方程组很难解。比较大时,方程组很难解。由于节点互异,行列式不等于由于节点互异,行列式不等于0 0,方程组存在唯一解,方程组存在唯一解以上方程组的系数行列式为范德蒙行列式以上方程组的系数行列式为范德蒙行列式MATLAB插值插值2/10/2023三、三、Lagran
5、ge(拉格朗日)(拉格朗日)插值插值1拉格朗日线性插值拉格朗日线性插值插值多项式:插值多项式:直线的两点式表达式直线的两点式表达式分别称为节点分别称为节点x0和和x1的一次的一次插值插值基函数基函数插值函数为基函数的线性组合,组合系数插值函数为基函数的线性组合,组合系数就是对应节点上的函数值就是对应节点上的函数值在在 上,过点上,过点MATLAB插值插值2/10/20232.拉格朗日二次插值拉格朗日二次插值二次插值多项式:二次插值多项式:对于对于 ,可假定:,可假定:基函数满足:基函数满足:基函数如何确定?基函数如何确定?从而从而过点过点MATLAB插值插值2/10/2023r 构造一组插值基
6、函数(构造一组插值基函数(n次次多项式多项式)缺点:缺点:当当 n 比较大时,插值多项式比较大时,插值多项式Ln(x)的收敛性与稳定性变的收敛性与稳定性变 差,逼近效果不理想。差,逼近效果不理想。li(x)是是 n次次多项式多项式,且满足且满足r 由插值基函数由插值基函数li(x)构造构造 n次拉格朗日插值次拉格朗日插值多项式多项式如下:如下:3.拉格朗日拉格朗日n 次插值次插值MATLAB插值插值2/10/2023r 拉格朗日拉格朗日n次插值的误差估计次插值的误差估计拉格朗日插值产生的截断误差为拉格朗日插值产生的截断误差为Rn(x),则,则MATLAB插值插值2/10/2023 连接相邻数据
7、点连接相邻数据点(x j-1,y j-1)、(x j,y j)得到得到n条线段,它条线段,它们组成一条折线。把区间们组成一条折线。把区间a,b上这上这n条折线段表示的函数称条折线段表示的函数称为被插函数为被插函数g(x)关于这关于这n+1个数据点的分段线性插值函数,记个数据点的分段线性插值函数,记作作I(x),则,则缺点:缺点:分段线性插值函数在节点处的一阶导数一般不存在,分段线性插值函数在节点处的一阶导数一般不存在,光滑性不高光滑性不高 I(xi)=g(xi),且,且I(x)在在a,b上连续。上连续。四、四、分段线性分段线性插值插值MATLAB插值插值2/10/2023 若构造插值基函数若构
8、造插值基函数则则I(x)可写成如下形式:可写成如下形式:MATLAB插值插值2/10/2023五、三次样条插值五、三次样条插值 是一个细的、可弯曲的木制或塑料条,在飞机或轮是一个细的、可弯曲的木制或塑料条,在飞机或轮船等的设计制造过程中为描绘出光滑的外形曲线船等的设计制造过程中为描绘出光滑的外形曲线(放样放样)所用的工具所用的工具 从物理上讲,样条满足插值点的约束,同时使势能从物理上讲,样条满足插值点的约束,同时使势能达到最小达到最小 三次样条本质上是一段一段的三次多项式拼合而成三次样条本质上是一段一段的三次多项式拼合而成的曲线,在拼接处的曲线,在拼接处,不仅函数是连续的不仅函数是连续的,且一
9、阶和二阶导且一阶和二阶导数也是连续的数也是连续的,19461946年年,Schoenberg,Schoenberg将样条引入数学将样条引入数学,即即所谓的所谓的样条函数样条函数1.什么是样条?什么是样条?MATLAB插值插值2/10/20232.三次样条插值函数三次样条插值函数MATLAB插值插值2/10/20233.三次样条插值原理三次样条插值原理在在n n个小区间构造个小区间构造S(xS(x),共有,共有n n个三次多项式,需确定个三次多项式,需确定4n4n个参数个参数在所有节点上在所有节点上n+1个方程个方程在除端点外的节点上在除端点外的节点上3(n-1)个方程个方程这样就得到这样就得到
10、4n-2个方程,为保证待定参数的唯一性,还差两个个方程,为保证待定参数的唯一性,还差两个方程。为此,常用的方法是对边界节点除函数值外附加要求,方程。为此,常用的方法是对边界节点除函数值外附加要求,这就是所谓的边界条件。根据实际问题的不同,三次样条插值这就是所谓的边界条件。根据实际问题的不同,三次样条插值常用到下列常用到下列三类边界条件三类边界条件。MATLAB插值插值2/10/2023r m边界条件:边界条件:即给定端点处的一阶导数值。即给定端点处的一阶导数值。r M边界条件:边界条件:即给定端点处的二阶导数值。即给定端点处的二阶导数值。r 周期边界条件周期边界条件:当当y=g(x)是以是以
11、b-a=x0 -xn 为周期的周期函数时,要求为周期的周期函数时,要求S(x)也是周期函数,故端点处要满足也是周期函数,故端点处要满足S(a+0)=S(b-0),S(a+0)=S(b-0),此条件称为周期条件此条件称为周期条件。MATLAB插值插值2/10/2023r 拉格朗日插值拉格朗日插值Matlab程序程序1:function y=lagrange(x0,y0,x)ii=1:length(x0);y=zeros(size(x);for i=ii ij=find(ii=i);y1=1;for j=1:length(ij)y1=y1.*(x-x0(ij(j);end y=y+y1*y0(i)
12、/prod(x0(i)-x0(ij);end六、一维插值的六、一维插值的Matlab实现实现1.拉格朗日插值拉格朗日插值对于拉格朗日多项式插值法,对于拉格朗日多项式插值法,Matlab中没有专用指令,可根中没有专用指令,可根据其算法编制如下程序据其算法编制如下程序MATLAB插值插值2/10/2023r 拉格朗日插值拉格朗日插值Matlab程序程序2:function y=lagrange(x0,y0,x)n=numel(x0);m=numel(x);x0=x0(:);y0=y0(:);x=x(:);y=zeros(n,m);for i=1:n y(i,:)=y0(i)*prod(repmat
13、(x,n-1,1)-.repmat(x0(1:n=i),1,m)/prod(x0(1:n=i)-x0(i);endy=sum(y);%替代写法:替代写法:%fun=(i)y0(i)*prod(repmat(x,n-1,1)-.%repmat(x0(1:n=i),1,m)/prod(x0(1:n=i)-x0(i);%y=sum(cell2mat(arrayfun(fun,1:n,UniformOutput,0);MATLAB插值插值2/10/2023例例3-1 考虑一个著名的例子考虑一个著名的例子,f(x)=1/(1+25x2),-1=x=1,假,假设已知其中一些点的坐标,则可以采用下面的命令进
14、行设已知其中一些点的坐标,则可以采用下面的命令进行lagrange插值插值 x0=-1+2*0:10/10;y0=1./(1+25*x0.2);x=-1:.01:1;y=lagrange(x0,y0,x);ya=1./(1+25*x.2);plot(x,ya,x,y,:)legend(原始图像原始图像,Lagrange插值插值)拟合效果不好拟合效果不好r Runge(龙格)现象(龙格)现象MATLAB插值插值2/10/2023r Matlab命令:命令:用用Matlab实现分段线性插值不需要编制函数程序,实现分段线性插值不需要编制函数程序,Matlab中有中有现成的一维插值函数现成的一维插值函
15、数interp1。y=interp1(x0,y0,x,method)method指定插值的方法,默认为线性插值。其值可为:指定插值的方法,默认为线性插值。其值可为:nearest 最近点等值方式插值最近点等值方式插值 linear 线性插值线性插值 spline 三次样条插值三次样条插值 cubic(新版本改为(新版本改为pchip)三次三次Hermite插值插值或立方插值。或立方插值。所有的插值方法要求所有的插值方法要求x0是单调的。是单调的。样条插值还可以这样:样条插值还可以这样:y=spline(x0,y0,x)2.分段线性插值与样条插值分段线性插值与样条插值注意:注意:向量向量x0 是
16、单调的。若是单调的。若y0为矩阵,则对为矩阵,则对y0的每一列进行插值,的每一列进行插值,x 可以为向量。可以为向量。MATLAB插值插值2/10/2023例例3-2 对例对例3-1中数据进行分段线性插值和样条插值中数据进行分段线性插值和样条插值 x0=-1+2*0:10/10;y0=1./(1+25*x0.2);x=-1:.01:1;y=interp1(x0,y0,x);ya=1./(1+25*x.2);subplot(2,1,1)plot(x,ya,x,y,r:)title(分段线性插值分段线性插值)legend(原始图像原始图像,分段线性插值分段线性插值)y=interp1(x0,y0,
17、x,spline);subplot(2,1,2)plot(x,ya,x,y,r:)title(样条插值样条插值)legend(原始图像原始图像,样条插值样条插值)MATLAB插值插值2/10/2023MATLAB插值插值2/10/2023例例3-3 已知某转子流量计在已知某转子流量计在100-1000 mL/min流量范围内,流量范围内,刻度值与校正值有如下关系(见表刻度值与校正值有如下关系(见表3-1)。试用线性插值法)。试用线性插值法计算流量计的刻度值为计算流量计的刻度值为780时,实际流量为多少?时,实际流量为多少?刻度值校正值刻度值校正值100105.3600605.8200207.2
18、700707.4300308.1800806.7400406.9900908.0500507.51000107.9表表3-1 例例3-3数据数据MATLAB插值插值2/10/2023解:解:Matlab计算程序如下:计算程序如下:X=100,200,300,400,500,600,700,800,900,1000;Y=105.3,207.2,308.1,406.9,507.5,605.8,707.4,806.7,908.0,107.9;Xk=780;Yk=interp1(X,Y,Xk)执行结果:执行结果:Yk=786.8400这里:这里:X和和Y分别表示样本点的刻度值和校正值;分别表示样本点的
19、刻度值和校正值;Xk和和Yk分别表示插值点的刻度值和校正值。分别表示插值点的刻度值和校正值。MATLAB插值插值2/10/2023r Matlab命令:命令:y=interp1(x0,y0,x,spline);y=spline(x0,y0,x);pp=csape(x0,y0,conds);pp=csape(x0,y0,conds,valconds);y=ppval(pp,x)。pp=csapi(x0,y0);y=fnval(pp,x)或或 y=csapi(x0,y0,x);%不带边界条件的三次样条插值。不带边界条件的三次样条插值。其中其中x0,y0是已知数据点,是已知数据点,x是插值点,是插值
20、点,y是插值点的函数值。是插值点的函数值。3.三次样条插值三次样条插值多项式函数多项式函数(polyfun)样条工具箱(样条工具箱(spline)MATLAB插值插值2/10/2023注注:对于三次样条插值,我们提倡使用函数对于三次样条插值,我们提倡使用函数csape,csape的返回的返回值是值是pp形式,要求插值点的函数值,必须调用函数形式,要求插值点的函数值,必须调用函数ppval。pp=csape(x0,y0);%使用默认的边界条件,即使用默认的边界条件,即Lagrange边界条件。边界条件。pp=csape(x0,y0,conds,valconds)其中由其中由conds指定插值的边
21、界条件,其值可为:指定插值的边界条件,其值可为:complete%边界为一阶导数边界为一阶导数 一阶导数的值在一阶导数的值在valconds参数中给出,参数中给出,若忽略若忽略valconds参数,则按缺省情况处理。参数,则按缺省情况处理。not-a-knot%非扭结条件非扭结条件 periodic%周期条件周期条件 second%边界为二阶导数边界为二阶导数,二阶导数的值在二阶导数的值在valconds参数中给出,若参数中给出,若 忽略忽略valconds参数,二阶导数的缺省值为参数,二阶导数的缺省值为0,0。MATLAB插值插值2/10/2023 注注:对于一些特殊的边界条件,其调用详细情
22、况请使用帮助对于一些特殊的边界条件,其调用详细情况请使用帮助 help csape例例3-4 机床加工待加工零件的外形根据工艺要求由一组数据机床加工待加工零件的外形根据工艺要求由一组数据(x,y)给出(在平面情况下),用程控铣床加工时每一刀只能沿给出(在平面情况下),用程控铣床加工时每一刀只能沿x方向和方向和y 方向走非常小的一步,这就需要从已知数据得到加工方向走非常小的一步,这就需要从已知数据得到加工所要求的步长很小的所要求的步长很小的(x,y)坐标。坐标。x035791112131415y01.21.72.02.12.01.81.21.01.6表中给出的表中给出的x,y数据位于机翼断面的下
23、轮廓线上,假设需要得到数据位于机翼断面的下轮廓线上,假设需要得到x坐标每改变坐标每改变0.1时的时的y坐标。试完成加工所需数据,画出曲线,坐标。试完成加工所需数据,画出曲线,并求出并求出x=0处的曲线斜率和处的曲线斜率和 范围内范围内y 的最小值。的最小值。要求用三次样条插值方法。要求用三次样条插值方法。要求用三次样条插值方法。要求用三次样条插值方法。表表3-2 例例3-4数据数据MATLAB插值插值2/10/2023 Matlab程序程序: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:
24、15;y1=interp1(x0,y0,x,spline);pp1=csape(x0,y0);%样条插值工具箱函数样条插值工具箱函数y2=ppval(pp1,x);%计算计算x对应的对应的y值值pp2=csape(x0,y0,second);y3=ppval(pp2,x);xydata=x,y1,y2,y3subplot(1,2,1)plot(x0,y0,+,x,y1)title(Spline1)subplot(1,2,2)plot(x0,y0,+,x,y2)title(Spline2)dx=diff(x);dy=diff(y2);dy_dx=dy./dx;dy_dx0=dy_dx(1)yte
25、mp=y2(13=x&x=15);ymin=min(ytemp);xmin=x(y2=ymin);xymin_1315=xmin,ymin 部分部分结果结果dy_dx0=0.4972xymin_1315=13.8000 0.9851MATLAB插值插值2/10/2023MATLAB插值插值2/10/2023例例3-5 假设已知的数据点来自函数。假设已知的数据点来自函数。x=0:.12:1;y=(x.2-3*x+5).*exp(-5*x).*sin(x);%原始数据原始数据xdat=0:.02:1;ydat=(xdat.2-3*xdat+5).*exp(-5*xdat).*sin(xdat);%
26、插值点坐标插值点坐标y1=interp1(x,y,xdat);%线性插值结果线性插值结果 y2=interp1(x,y,xdat,pchip);%三次三次Hermite插值结果插值结果y3=interp1(x,y,xdat,spline);%三次样条插值结果三次样条插值结果y4=interp1(x,y,xdat,nearest);%最近点等值方式插值结果最近点等值方式插值结果plot(xdat,y1,y2,y3,y4,:,x,y,o,xdat,ydat)legend(线性插值线性插值,三次三次Hermite插值插值,三次样条插值三次样条插值,最近点等值方式插值最近点等值方式插值,原始数据原始数
27、据,插值点插值点)产生数据的命令为:产生数据的命令为:x=0:.12:1;y=(x.2-3*x+5).*exp(-5*x).*sin(x);试根据生成的数据进行插值处理,得出较平滑的曲线。试根据生成的数据进行插值处理,得出较平滑的曲线。Matlab程序程序:MATLAB插值插值2/10/2023MATLAB插值插值2/10/2023例例3-63-6 为了绘制乙醇水二元溶液的气液平衡曲线,查到为了绘制乙醇水二元溶液的气液平衡曲线,查到2424点点实验数据如表实验数据如表4-44-4所示。表中所示。表中x x和和y y分别表示液相和气相中乙醇摩分别表示液相和气相中乙醇摩尔分数。但表中缺少低浓度点的
28、数据,且数据分布不均匀。试尔分数。但表中缺少低浓度点的数据,且数据分布不均匀。试用三次样条插值法对以下插值点插值并画出乙醇水二元溶液用三次样条插值法对以下插值点插值并画出乙醇水二元溶液的气液平衡曲线。的气液平衡曲线。表表3-3 3-3 乙醇水二元溶液的气液平衡曲线实验值乙醇水二元溶液的气液平衡曲线实验值MATLAB插值插值2/10/2023x=0,0.018,0.054,0.078,0.111,0.147,0.172,0.203,0.210,0.241,.0.297,0.320,0.381,0.391,0.423,0.425,0.497,0.507,0.567,0.685,.0.767,0.8
29、38,0.882,0.894;y=0,0.163,0.328,0.386,0.452,0.491,0.512,0.524,0.520,0.546,.0.573,0.577,0.600,0.608,0.621,0.623,0.653,0.656,0.684,0.743,.0.796,0.846,0.884,0.894;xk=0.010,0.018,0.020,0.050,0.078,0.100,0.118,0.200,0.300,0.400,.0.425,0.500,0.600,0.700,0.800,0.838,0.866,0.899;yk=spline(x,y,xk);%yk=csapi(x
30、,y,xk);%pp=csape(x,y,second);%yk=ppval(pp,xk);plot(x,y,*,xk,yk,r);grid;legend(样本点,三次样条插值,Location,NorthWest)Matlab程序一程序一:这里:这里:x和和y分别表示样本点的液相和气相中乙醇分别表示样本点的液相和气相中乙醇摩尔分数,摩尔分数,xk和和yk分别表示插值点的液相和气相分别表示插值点的液相和气相中乙醇摩尔分数。中乙醇摩尔分数。MATLAB插值插值2/10/2023配合界面演示配合界面演示MATLAB插值插值2/10/2023MATLAB插值插值2/10/2023第二节第二节 二维插
31、值二维插值MATLAB插值插值2/10/2023 已知二元函数已知二元函数g(x,y)在某矩形区域在某矩形区域R=a,b c,d上互异上互异节点节点(x i,y j)的函数值的函数值 z ij,如何求出在,如何求出在R 上任一点上任一点(x,y)处处的函数值的函数值g(x,y)的近似值。的近似值。根据已知数据点根据已知数据点(x i,y j,z ij),构造一个具有一定光滑性、,构造一个具有一定光滑性、简单易于计算的函数简单易于计算的函数z=f(x,y)(已知曲面)作为(已知曲面)作为z=g(x,y)(未知曲面)的近似,使曲面(未知曲面)的近似,使曲面z=f(x,y)通过已知数据点,即通过已知
32、数据点,即 f(x i,y j)=z ij。然后计算点。然后计算点(x,y)的函数值的函数值f(x,y)作为作为g(x,y)的近似值。的近似值。一、数学描述及基本思想一、数学描述及基本思想数学描述数学描述2基本思想基本思想MATLAB插值插值2/10/2023二、网格节点插值二、网格节点插值MATLAB插值插值2/10/2023 已知已知 m n 个节点个节点 其中其中互不相同,不妨设互不相同,不妨设 构造一个二元函数构造一个二元函数通过全部已知节点通过全部已知节点,即即再用再用计算插值,即计算插值,即MATLAB插值插值2/10/2023q 最邻近点插值最邻近点插值网格节点插值法主要有以下几
33、种形式:网格节点插值法主要有以下几种形式:q 分片线性插值分片线性插值q 分片双线性插值分片双线性插值q 分片双三次样条插值分片双三次样条插值MATLAB插值插值2/10/2023三、散乱节点插值三、散乱节点插值MATLAB插值插值2/10/2023已知已知n n个节点个节点其中其中互不相同,互不相同,构造一个二元函数构造一个二元函数通过全部已知节点通过全部已知节点,即即再用再用计算插值,即计算插值,即MATLAB插值插值2/10/2023四、二维插值的四、二维插值的Matlab实现实现1.网格节点插值网格节点插值要求要求x0,y0单调;单调;x,y可取可取为矩阵,或为矩阵,或x 取行向量,取
34、行向量,y 取为列向取为列向量,量,x,y 的值分别不能超出的值分别不能超出x0,y0的范围。的范围。(1)z =interp2(x0,y0,z0,x,y,method)被插值点被插值点插值方法插值方法插值插值节点节点被插值点的被插值点的函数值函数值nearestnearest 最邻近插值最邻近插值linearlinear 双线性插值双线性插值splinespline 三次样条三次样条cubiccubic 双三次插值双三次插值缺省时缺省时 双线性插值双线性插值MATLAB插值插值2/10/2023例例3-7 假设由二元函数假设由二元函数可以计算出一些较稀疏的网格数据,试根据这些数据对整个函可以
35、计算出一些较稀疏的网格数据,试根据这些数据对整个函数曲面进行各种插值拟合,比较拟合结果,并进行误差分析。数曲面进行各种插值拟合,比较拟合结果,并进行误差分析。原始数据与原始面图原始数据与原始面图x,y=meshgrid(-3:.6:3,-2:.4:2);z=(x.2-2*x).*exp(-x.2-y.2-x.*y);%原始数据原始数据figure(1)surf(x,y,z)%原始数据面图原始数据面图axis(-3,3,-2,2,-0.7,1.5)%设置坐标轴范围设置坐标轴范围box ontitle(原始数据面图原始数据面图)xdat,ydat=meshgrid(-3:.2:3,-2:.2:2)
36、;%插值点坐标插值点坐标MATLAB插值插值2/10/2023线性插值面图线性插值面图zdat1=interp2(x,y,z,xdat,ydat);%线性插值结果线性插值结果 figure(2)surf(xdat,ydat,zdat1)%线性插值面图线性插值面图axis(-3,3,-2,2,-0.7,1.5)%设置坐标轴范围设置坐标轴范围box ontitle(线性插值面图线性插值面图)立方插值面图立方插值面图zdat2=interp2(x,y,z,xdat,ydat,cubic);%立方插值算法结果立方插值算法结果figure(3)surf(xdat,ydat,zdat2)%立方插值算法面图
37、立方插值算法面图axis(-3,3,-2,2,-0.7,1.5)%设置坐标轴范围设置坐标轴范围box ontitle(立方插值算法面图立方插值算法面图)MATLAB插值插值2/10/2023样条插值面图样条插值面图zdat3=interp2(x,y,z,xdat,ydat,spline);%样条插值结果样条插值结果figure(4)surf(xdat,ydat,zdat3)%样条插值面图样条插值面图axis(-3,3,-2,2,-0.7,1.5)%设置坐标轴范围设置坐标轴范围box ontitle(样条插值面图样条插值面图)样条插值误差面图样条插值误差面图zdat=(xdat.2-2*xdat
38、).*exp(-xdat.2-ydat.2-xdat.*ydat);figure(6)surf(xdat,ydat,abs(zdat-zdat3)%样条插值对应的误差面图样条插值对应的误差面图axis(-3,3,-2,2,0,0.025)%设置坐标轴范围设置坐标轴范围box ontitle(样条插值算法对应的误差面图样条插值算法对应的误差面图)插值点处的函数值插值点处的函数值MATLAB插值插值2/10/2023MATLAB插值插值2/10/2023MATLAB插值插值2/10/2023 z =ppval(pp,x,y)(2)pp =csape(x0,y0,z0,conds,valconds)
39、其中其中x0,y0,z0为已知的数据,为已知的数据,x0为为m维向量,维向量,y0为为n维向量,维向量,z0为为mn维矩阵。维矩阵。x,y为一维数组表示插值点(为一维数组表示插值点(x为行向量,为行向量,y为为列向量。列向量。若不考虑边界条件,以上两条语句可用下面命令代替若不考虑边界条件,以上两条语句可用下面命令代替z=csapi(x0,y0,z0,x,y);被插值点被插值点插值边界条件及取插值边界条件及取值,同一维情形值,同一维情形插值节点插值节点被插值点被插值点的函数值的函数值系数矩阵系数矩阵MATLAB插值插值2/10/2023例例3-8 在一丘陵地带测量高程,在一丘陵地带测量高程,x和
40、和y方向每隔方向每隔100米测一个点,米测一个点,得高程如下表,试拟合一曲面,确定合适的模型,并由此找出得高程如下表,试拟合一曲面,确定合适的模型,并由此找出最高点和该点的高程。最高点和该点的高程。高程x100200300400500y100636697624478450200698712630478420300680674598412400400662626552334310表表3-4 3-4 某地高程数据某地高程数据MATLAB插值插值2/10/2023原始数据面图原始数据面图x=100:100:500;y=100:100:400;z=636 697 624 478 450698 712
41、630 478 420680 674 598 412 400662 626 552 334 310;%原始数据原始数据figure(1)surf(x,y,z)title(原始数据面图原始数据面图)MATLAB插值插值2/10/2023插值数据面图插值数据面图xdat,ydat=ndgrid(100:10:500,100:10:400);%插值点网格数据插值点网格数据%pp=csape(x,y,z);%zdat=fnval(pp,100:10:500,100:10:400);%经插值运算得到的经插值运算得到的插值点处函数值插值点处函数值zdat=csapi(x,y,z,100:10:500,10
42、0:10:400);%经插值运算得到经插值运算得到的插值点处函数值的插值点处函数值figure(2)surf(xdat,ydat,zdat)%插值面图插值面图title(插值数据面图插值数据面图)MATLAB插值插值2/10/2023最高点的高程求解最高点的高程求解x_temp=100:500;%以以1为步长产生数据为步长产生数据y_temp=100:400;z_temp=csapi(x,y,z,100:500,100:400);z_max=max(z_temp(:);x_max,y_max=find(z_temp=z_max);xyz_max=x_temp(x_max),y_temp(y_m
43、ax),z_max%最高点和该点的高程最高点和该点的高程MATLAB插值插值2/10/2023z =griddata(x0,y0,z0,x,y,method)2.散乱节点插值散乱节点插值被插被插值点值点插值插值方法方法插值插值节点节点被插值点被插值点的函数值的函数值nearest 最邻近插值最邻近插值linear 双线性插值双线性插值cubic 双三次插值双三次插值v4 Matlab提供的插值方法提供的插值方法缺省时缺省时 双线性插值双线性插值MATLAB插值插值2/10/2023其中其中x0,y0,z0为已知的数据,这里并不要求是网格型的,可以为已知的数据,这里并不要求是网格型的,可以是任意
44、分布的,均由向量给出。是任意分布的,均由向量给出。x,y为插值点,可以是单个为插值点,可以是单个点,可以是向量或网格型矩阵,得出的应该维数和点,可以是向量或网格型矩阵,得出的应该维数和x,y一致,一致,表示插值的结果。表示插值的结果。method 表示插值算法,可以为表示插值算法,可以为linear,cubic,nearest,v4其中其中v4是是Matlab4.0版本中提供的插值算法,公认效果较好,版本中提供的插值算法,公认效果较好,但没有一个正式的名称,所以用但没有一个正式的名称,所以用v4表明该算法。表明该算法。MATLAB插值插值2/10/2023例例3-9 已经有了已经有了 x,y,
45、z 的坐标值,如的坐标值,如 x=1 2 3 4;y=2 3 4 5;z=5 6 7 4。由这些坐标点通过插值方法拟合出三维曲面由这些坐标点通过插值方法拟合出三维曲面。Matlab程序程序x=1 2 3 4;y=2 3 4 5;z=5 6 7 4;x1,y1=meshgrid(1:0.1:4,2:0.1:5);z1=griddata(x,y,z,x1,y1,v4);surf(x1,y1,z1)MATLAB插值插值2/10/2023第三节第三节 高维插值高维插值MATLAB插值插值2/10/2023 Matlab命令命令:1.z=interp3(x0,y0,z0,x1,y1,method)2.z=interpn(x1,x2,xn,z0,xd1,xd2,xdn,method)3.z=csapi(x1,x2,xn,z0,xd1,xd2,xdn);4.v=griddata3(x0,y0,z0,v0,x,y,z,method);5.v=griddata3(X,Y,XI method);具体用法请参考具体用法请参考Matlab帮助帮助二、高维插值的二、高维插值的Matlab实现实现一、数学描述及插值原理一、数学描述及插值原理内容略内容略
限制150内