第四讲曲线曲面的插值与拟合方法

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,第四讲 插值与拟合之插值(上),内容：,插值是离散函数逼近的重要方法，利用它,可通过函数在有限个点处的取值状况，估,算出函数在其他点处的近似值,目的：,学习插值的基本思想和方法，掌握,Matlab,的一维/二维等距和非等距插值函数,要求：,掌握,M,atlab,插值函数，处理插值应用问题,了解拉格朗日和分段线性插值的基本思想,了解三次样条插值的提法和思路,掌握插值函数,interp,interp1,interp2,griddata,掌握水塔用水量的计算,(,-,体积,-,流速,-,积分,),关于插值与拟合的,区别,面对工程实践和科学计算中的采集得到数据,(,xi,yi),，,我们总是试图去揭示,x,与,y,之间的关系,，即用近似的,y=f(x),来表示，那么我们通常可以采用两种方法：,插值与拟合,插值与拟合的,区别,在于,插值试图去通过,已知点,了解,未知,点,处的函数值；而拟合则在于在,整体上用某种,已知函数,去拟合数,据点列所在,未知函数,的性态。,关键区别,在于插值要求,必须经过,已知点列，拟合只求,尽量靠近,不必经过！拟合将在本讲下介绍,引例1 函数查表问题：,已知标准正态分布函数表，求,表中没有的值,(2.34)=0.99036 (2.35)=0.99061,求,(2.3457) (2.35-2.3457)/(2.35-2.34)*,(2.34),+,(2.3457-2.34)/(2.35-2.34)*,(2.3,5),引例2 地图绘制问题：,假如我们在地图边界获取了一些边界点的坐标，连接这些边界点形成闭合曲线，可以用来近似表示真实边界线，如何,更准确地逼近,真实边界线？,函数查表与地图边界线绘制,(2.3457)=？,如何更准确地,逼近,真实边界线？,插值在数码图像放大中的应用,引例,3,图像插值放大：,数码相机运用插值的方法可以,创造出,比传感器实际像素更多的图像，这种处理称为“数码变焦”。,106*40,原始图像：,左边：,最近邻插值,放大,450%,右边：,双三次插值,放大,450%,插值在图像三维重建中的应用,Surface recostruction from scattered points cloud,分段线性插值和拉格朗日插值,分段线性插值：,用直线(线性)连接数据点列上相邻的两点。,比如在两点,x,i-1,x,i,上线性插值函数为,拉格朗日插值：,用,n,次拉格朗日插值多项式,连接数据点列上相邻的,n+1,个点。,Pszjs71,拉格朗日插值基函数的构造,比如 在三个点,x,0,x,1,x,2,上,lagrange,插值函数为,(,线性插值,是拉格朗日插值,最简单,的情形),分段三次埃尔米特插值条件数,分段三次埃尔米特插值：,线性插值,在每一小段上(两点之间)，用到,2个,条件,q(x,i,)=y,i,，,所以确定了一个,线性,插值函数；,三次埃尔米特插值,在每一小段上，用到,4个,条件,q(x,i,)=y,i,， q(x,i,)=y,i,所以确定一个,3次多项式,插值函数。,分段插值主要是为了避免高次插值可能出现的,大幅度振荡,现象，在实际应用中通常采用分段低次插值来,提高近似程度,，比如可用分段线性插值或分段三次埃尔米特插值来逼近已知函数，但它们的,总体光滑性较差,，为了克服这一缺点，三次样条插值成为比较理想的工具。,三次样条(spline)插值的概念,样条的概念出自工程设计和机械加工(飞机、船舶外形曲线设计)中的绘图工具(曲线尺)，简单说就是具有,连续二阶导数的三次插值多项式函数,。,三次样条(spline)插值的条件数,首先从段数,n=2,分析：我们知道在每一小段的三次多项式有4个系数，所以如下图，总共需要有4*,2,=,8,个方程来确定； 由,q(x,i,)=y,i,可以确定2*,2,=,4,个方程，又由内部节点,q,1,(x,i,)= q,2,(x,i,),和,q,1,(x,i,)= q,2,(x,i,),可以确定2*(,2,-1)=,2,个方程，看来剩下的8-(4+2)=2个方程只有靠,外部给定,(边界条件)了,q,1,q,2,x,0,x,1,x,2,一维曲线等距插值函数interp,interps syntax,One-dimensional,r times,longer,data interpolation,y = interp(y,r,),题例,在原始数据点中增倍插值,x=0:0.001:1;,y=sin(2*pi*30*x)+sin(2*pi*60*x);,yi=interp(y,4);,subplot(1,2,1);,stem(y(1:30);,title(Original Points);,subplot(1,2,2);,stem(yi(1:120);,title(Interpolated Points);,一维曲线等距插值函数interp1,interp1s syntax,One-dimensional data interpolation,yi = interp1(,x,y,xi,method,),nearest Nearest neighbor interpolation,linear Linear interpolation (default),spline, Cubic spline interpolation,cubic Piecewise cubic Hermite interpolation,题例,在一天,24,小时内,从零点开始每间隔,2,小时测得的环境温度，推测在,15,点,6,分的的温度,x=0:2:24; y=12,9,9,10,18,24,28,27,25,20,18,15,13;,plot(x,y,-ro); hold on;,xi=15.1;,yi=interp1(x,y,xi,spline),xi=0:1/3600:24; yi=interp1(x,y,xi,spline); plot(xi,yi,b-);,二维曲面等距插值函数interp2,interp,2,s syntax,Two-dimensional data interpolation,ZI = interp2(,X,Y,Z,XI,YI,method,),nearest Nearest neighbor interpolation,linear Bilinear interpolation (default),spline, Cubic spline interpolation,cubic Bicubuc interpolation,二维曲面等距插值函数interp2,动画展示：三维空间中的曲面等距格点,二维曲面等距插值函数interp2,题例,粗糙山顶曲面的平滑处理,(,等距情形,),load mountain.mat,%,载入山顶地形数据,mesh(x,y,z),%,绘制原始山顶地形图,xi=linspace(0,5,50);,yi=linspace(0,6,80);,xii,yii=meshgrid(xi,yi);,zii=interp2(x,y,z,xii,yii,spline);,%,三次样条插值,figure;,surf(xii,yii,zii) %,绘制平滑处理后的山顶曲面,hold on;,xx,yy=meshgrid(x,y);,plot3(xx,yy,z+0.1,ob);,二维曲面等距插值函数interp2,题例,粗糙山顶曲面的平滑处理,(,等距情形,),二维曲面散乱插值函数griddata,griddatas syntax,Data interpolation,for scattered points,ZI = griddata(x,y,z,XI,YI),XI,YI,ZI = griddata(x,y,z,xi,yi),. = griddata(.,method),linear Triangle-based linear interpolation,cubic, Triangle-based cubic (default),nearest Nearest neighbor,v4 MATLAB 4 griddata method,MATLAB,二维插值函数,griddata,，,可以将平面或曲面上的,散乱点,插值为,规则网格,二维曲面散乱插值函数griddata,题例,粗糙山顶曲面的平滑处理,(,散乱情形,),rand(seed,0),x = rand(100,1)*4-2;,y = rand(100,1)*4-2;,z = x.*exp(-x.2-y.2);,plot3(x,y,z,o);,hold on,ti = -2:.25:2;,XI,YI = meshgrid(ti,ti);,ZI = griddata(x,y,z,XI,YI);,mesh(XI,YI,ZI);,二维曲面散乱插值函数griddata,题例,墨西哥草帽的平滑处理,(,散乱情形,),x = rand(100,1)*16 - 8;,y = rand(100,1)*16 - 8;,r = sqrt(x.2 + y.2) + eps;,z = sin(r)./r;,plot3(x,y,z,.,MarkerSize,15),hold on,xlin = linspace(min(x),max(x),33);,ylin = linspace(min(y),max(y),33);,X,Y = meshgrid(xlin,ylin);,Z = griddata(x,y,z,X,Y,cubic);,mesh(X,Y,Z); axis tight;,南半球气旋变化的可视图形,山区地貌的可视化图形,水塔用水量估计通用程序,通用程序,tbp69.m,可近似计算时间段内的用水量,格式为：,tbp69(ts,tf),其中,ts,为起点时间，,tf,为终点时间,实验一：水塔用水量估计,Thats all3Q!,第四讲 插值与拟合之拟合(下),内容：,拟合是离散函数逼近的重要方法，利用它,可通过函数在有限个点处的取值状况，,拟,合出近似替代函数，,进而估算出函数在其,他点处的近似值。,目的：,学习拟合的基本思想和方法，掌握,Matlab,的多项式,/,一般拟合函数,/,曲线拟合工具箱,要求：,掌握,Matlab,拟合函数，处理拟合应用问题,了解基于最小二乘法则拟合的基本思想,掌握拟合函数,polyfit lsqcurvefit curvefit,掌握,cftool,曲线拟合工具箱,(,多目标函数多法则,),关于数据拟合的两个要素.,在工程实践和科学计算中，用某种,经验函数,解析式,y=f(x),来近似刻画,采集数据(,x,y),之间的关系的方法就叫拟合，所谓“拟合”有,“最贴近”,之意,。,与插值不同，拟合的主要目,标是要离散点,尽量靠近,拟合函,数。一般过程是，我们首先根,据采样点的散点分布图，大致,推测,x,与,y,之间的,经验函数,形式,(比如多项式、指数函数等)，,然后依据,某种法则,(比如最常用的最小二乘法则)，确定出的经验函数解析式中的待定参数。其中,经验函数,和,拟合法则,是拟合的两个关键要素！,引例,1,化合物浓度随时间变化的规律：,与插值面临的问题相似，我们被要求去,求解或预测,表格中没有的因变量对应值，与插值的解决思路不同，我们试图获得,比较完备,的解决方案：设计并求出离散数据点的近似替代函数，有了近似函数解析式，就可以进一步代值计算或作图分析。,为了揭示浓度,y,与时间,t,之间呈现的函数规律，我们首先作出散点图，帮助,分析和设计,经验函数,化合物浓度随时间变化的规律,化合物浓度随时间变化的规律,如图, 化合物浓度,y,随时间,t,大致呈抛物线状（二次函数）变化，这种分析和判断来自,已有经验,.,t=1:16;,c=4 6.4 8 8.4 9.28 9.5 9.7 9.86 10 10.2 10.32 10.42 10.5 10.55 10.58 10.6;,plot(t,c,-ro),化合物浓度随时间变化的规律,经验函数形式：,已经拟定为多项式函数：,y= at,2,+bt+ c,剩下的工作是确定拟合原则：,可选的法则很多，其中最常用的是最小二乘法则,(,method of,L,east,S,q,uares,),，即,各点偏差平方和最小,高斯和勒让德关于最小二乘法的发明权,化合物浓度随时间变化的规律,对经验函数形式确定的补充说明：,拟合函数解析式选用,什么形式,(用多项式，还是用幂函数?)，主要取决于,采样点的分布,无疑，那么如何求出这些含有待定参数的解析式呢？,把各点偏差的平方和最小作为一个目标函数，实际上考虑为极值问题，极值点导数为零。,具体计算时，我们在把经验函数用一系列,拟合基函数,线性表出同时，在,j,个采样点对待定参数,Cj,求偏导(=0),获取,j,个方程，进而解出,Cj，,具体参见数值计算,SZJSp9091,在本例中，已经拟定拟合的目标函数为多项式函数：,y= at,2,+bt+ c ，,所以只要解出三个待定参数,a,b,c，,问题即获解决,基于最小二乘的多项式拟合函数,基于最小二乘的多项式拟合函数,polyfit,：,Polynomial curve fitting .Syntax：,p = polyfit(x,y,n),其中,n,是拟合多项式的阶数，,不能超过,（散点数据对数,-1,）,下面回到化合物浓度随时间变化的引例：,t=1:16; c=4 6.4 8 8.4 9.28 9.5 9.7 9.86 10 10.2 10.32 10.42 10.5 10.55 10.58 10.6; plot(t,c,ko); hold on;,%,作散点图,p2=polyfit(t,c,2);,y2=poly2str(p2,t),%,作多次拟合比较,p5=polyfit(t,c,5);,y5=poly2sym(p5,t),f=inline(y5),ti=0:.001:20; plot(ti,polyval(p2,ti),b-,ti,f(ti),r-);,disp(,化合物在刻度,11.2,的浓度,近似值,为,num2str(f(11.2),disp(,化合物在刻度,17.8,的浓度,预测值,为,num2str(f(17.8),stem(11.2 17.8,f(11.2) f(17.8),r); xlabel(,时间,t); ylabel(,化合物浓度,c); title(,化合物浓度随时间变化的规律,),引例,2,确定医用薄膜渗透率的数学模型：,某种医用薄膜允许一种物质分子,从高浓度溶液,V,B,穿过薄膜向低浓度,溶液,V,A,中扩散。,通过单位面积膜,S,分子扩散的速度与膜两侧溶液的浓,度差成正比，,比例系数K表示薄膜,被该物质分子穿透的能力，称为,渗透率,，,定时测量薄膜,V,B,侧的溶液浓度值,C,B,，以此确定K的值VA=VB=1000cm,3,S=10cm,2,容器的B部分溶液浓度,C,B,的测试结果如下表,：,( CB,单位为,mg/,cm,3,),确定医用薄膜渗透率的数学模型,确定医用薄膜渗透率的数学模型,由质量守恒考察, t,t+,t,时间段,B,向,A,中渗透物质：,VA*CA(t+t)-VA*CA(t) = SKCB(t)-CA(t)t,推出,d,CA(t)/dt = SK/VA*CB(t)-CA(t),两边除以,t,令,t,0,又由质量守恒考察整个容器中物质总量始终不变：,V,A*CA(t)+,V,B*CB(t) =,V,A*aA+,V,B*aB,推出,CA(t) = aA+VB/VA*aB-VB/VA*CB(t),代入上式,2,推出,d,CB(t)/dt = -SK(1/VA+1/VB)CB(t)+SK(aA/VB+aB/VA),CB(0)=aB,初值条件 此带初值微分方程可由,dsolve,求解,在上式中，已知的包括,VA,，,VB,，,S,以及一组,t,和,CB(t),值,未知的包括,aA,，,aB,，,K,，下面通过数据拟合确定渗透率,K,确定医用薄膜渗透率的数学模型,在上式中，代入已知值,VA=VB=1000cm,3,，,S=10cm,2,令,a=(aA*VA+aB*VB)/(VA+VB),，,b=VA(aB-aA)/(VA+VB),简化之后的表达式为：,CB(t)=a+b*exp(-0.02*k*t),编写被调,M,文件,tbp79.m,function CB=tbp79(x,t),CB=x(1)+x(2)*exp(-0.02*x(3)*t);,编写主调,M,文件,fittbp79.m,（片段）,x=curvefit(tbp79,x0,t,CB),%curvefit拟合及图像,x=lsqcurvefit(tbp79,x0,t,CB),%lsqcurvefit拟合及图像,求解结果：,a=x(1)=0.0070; b=x(2)=-0.0030;,k=x(3)=0.1012,进一步求解：,aA=0.01,；,aB=0.004,最终数学模型：,CB(t)=0.007-0.003*exp(-0.002*t),基于最小二乘的一般拟合函数,基于最小二乘的一般拟合函数,lsqcurvefit,：,Solve nonlinear curve-fitting (data-fitting) problems in the least-squares sense.Syntax：,x = lsqcurvefit(fun,x0,xdata,ydata),x,resnorm = lsqcurvefit(.),范例,:,function F=myfun(x,xdata),F=x(1)*xdata.2+x(2)*sin(xdata)+x(3)*xdata.3;,%,下面是主调函数,fitmyfun.m,的部分,xdata = 3.6 7.7 9.3 4.1 8.6 2.8 1.3 7.9 10.0 5.4;,ydata = 16.5 150.6 263.1 24.7 208.5 9.9 2.7 163.9 325.0 54.3;,x0 = 10, 10, 10;% Starting guess,x,resnorm = lsqcurvefit(myfun,x0,xdata,ydata),MATLAB工具箱的版本更新.,同名升级,换名升级,CFTool曲线拟合工具箱简介,基于,MATLAB,的曲线拟合问题，已经提供独立的,toolbox,供调用，该,toolbox,采用,GUI,界面，功能强大，下面简单介绍如何使用该,Toolbox,解决一般曲线拟合问题。,在,command window,中键入指令,cftool,即可启动曲线拟合工具箱。在该集成环境里面，可以实现,多种经验函数，多种法则,的曲线拟合，实时绘制图像并进行误差分析。,需要注意,的是：在进入,Curve Fitting Toolbox,环境进行曲线拟合之前，需要,预先,在,workspace,输入或载入供拟合的,数据源,CFTool-选择Data导入数据,下面还是以引例的采样数据为例，进行演示：,t=1:16;y=4 6.4 8 8.4 9.28 9.5 9.7 9.86 10 10.2 10.32 10.42 10.5 10.55 10.58 10.6;,cftool,导入数据,绘制散点图,CFTool-选择Fitting拟合数据,进行拟合,拟合方法,结果和误差分析,这里可供选择的,拟合类型,和,可选参数,比较多,包括多项式函数,指数函数,幂函数等,如何确定,最优的方案,？,CFTool-拟合效果评价指标,SSE,- The sum of squares due to error. This statistic measures the deviation of the responses from the fitted values of the responses.,A value closer to 0 indicates a better fit.,R-square,- The coefficient of multiple determination. This statistic measures how successful the fit is in explaining the variation of the data.,A value closer to 1 indicates a better fit.,Adjusted R-square,- The degree of freedom adjusted R-square. A value closer to 1 indicates a better fit. It is generally the best indicator of the fit quality when you add additional coefficients to your model.,RMSE,- The root mean squared error.,A value closer to 0 indicates a better fit.,CFTool-选择Analysis分析数据,分析数据,这里可以计算拟合函数在采样点的一阶二阶导函数值，并绘制相应的函数图像,CFTool-拟合薄膜渗透数学模型,导入数据：,load bmstl.mat,cftool,拟合数据：,syms a b c d t,f = a*exp(b*t) + c*exp(d*t),cbvst=,subs,(f,a,b,c,d,0.008839,-0.0001336,-0.004818,-0.001422);,f =,inline,(cbvst),figure,t=100:1:1000;,plot(t,f(t),b-),CFTool-拟合薄膜渗透数学模型,分析数据：,实验一：,P87 Malthus,人口模型实验二：,P87,旧车价格预测,Thats all3Q!,