曲线曲面的插值与拟合方法次课PPT教案

1、曲线曲面的插值与拟合方法次课曲线曲面的插值与拟合方法次课 (2.3457)=？第1页/共43页第2页/共43页106*40原始图像：原始图像：左边：左边：最近邻插值最近邻插值放大放大450%右边：右边：双三次插值双三次插值放大放大450%第3页/共43页第4页/共43页分段线性插值：分段线性插值：用直线用直线(线性线性)连接数据点列上相邻的两点。连接数据点列上相邻的两点。比如比如在两点在两点xi-1,xi上线性插值函数为上线性插值函数为11111( ),0,1,2,.,iiiiiiiiiixxxxq xyyxxxxxxxin拉格朗日插值：拉格朗日插值：用用n次拉格朗日插值多项式次拉格朗日插值多

2、项式0()()nnkkkLxylx连接数据点列上相邻的连接数据点列上相邻的n+1个点。个点。Pszjs71第5页/共43页2200011221220010202110120122021()()()()()() ()()() ()() ()() ()() ()() ()kkkLxylxylxy lxylxxxxxLxyxxxxxxxxyxxxxxxxxyxxxx第6页/共43页第7页/共43页第8页/共43页q1q2x0 x1x2第9页/共43页interps syntaxOne-dimensional r times longer data interpolation y = interp(y

3、,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);第10页/共43页interp1s syntaxOne-dimensional data interpolation yi = interp1(x,y,xi,method)nearest Nea

4、rest neighbor interpolationlinear Linear interpolation (default)spline Cubic spline interpolationcubic 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;

5、yi=interp1(x,y,xi,spline),xi=0:1/3600:24; yi=interp1(x,y,xi,spline); plot(xi,yi,b-);第11页/共43页interp2s syntaxTwo-dimensional data interpolation ZI = interp2(X,Y,Z,XI,YI,method)nearest Nearest neighbor interpolationlinear Bilinear interpolation (default)spline Cubic spline interpolationcubic Bicubuc i

6、nterpolation第12页/共43页动画展示：三维空间中的曲面等距格点动画展示：三维空间中的曲面等距格点第13页/共43页题例题例 粗糙山顶曲面的平滑处理粗糙山顶曲面的平滑处理(等距情形等距情形)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,yi

7、i,zii) %绘制平滑处理后的山顶曲面绘制平滑处理后的山顶曲面hold on;xx,yy=meshgrid(x,y);plot3(xx,yy,z+0.1,ob);第14页/共43页题例题例 粗糙山顶曲面的平滑处理粗糙山顶曲面的平滑处理(等距情形等距情形)第15页/共43页griddatas syntaxData 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 linea

8、r interpolationcubic Triangle-based cubic (default) nearest Nearest neighbor v4 MATLAB 4 griddata methodMATLAB二维插值函数二维插值函数griddata，可以将平面或曲面上的可以将平面或曲面上的散乱点散乱点插值为插值为规则网格规则网格第16页/共43页题例题例 粗糙山顶曲面的平滑处理粗糙山顶曲面的平滑处理( (散乱情形散乱情形) )rand(seed,0)x = rand(100,1)*4-2; y = rand(100,1)*4-2;z = x.*exp(-x.2-y.2);plot3

9、(x,y,z,o);hold onti = -2:.25:2; XI,YI = meshgrid(ti,ti);ZI = griddata(x,y,z,XI,YI);mesh(XI,YI,ZI);第17页/共43页题例题例 墨西哥草帽的平滑处理墨西哥草帽的平滑处理( (散乱情形散乱情形) ) 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 onxlin = linspace(min(x),max(x),33

10、);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; 第18页/共43页第19页/共43页第20页/共43页通用程序通用程序tbp69.m可近似计算时间段内的用水量可近似计算时间段内的用水量格式为：格式为：tbp69(ts,tf) 其中其中ts为起点时间，为起点时间，tf为终点时间为终点时间第21页/共43页第22页/共43页第23页/共43页第24页/共43页第25页/共43页t=1:16;c=4 6.4 8 8.

11、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)第26页/共43页经验函数形式：经验函数形式：已经拟定为多项式函数：已经拟定为多项式函数：y= at2 +bt+ c剩下的工作是确定拟合原则：剩下的工作是确定拟合原则：可选的法则很多，其中最常用的是最小二乘法则可选的法则很多，其中最常用的是最小二乘法则(method of Least Squares)，即，即各点偏差平方和最小各点偏差平方和最小高斯和勒让德关于最小二乘法的发明权第27页/共43页第28页/共43页第29页/共43页第30页/共43

12、页由质量守恒考察由质量守恒考察 t,t+t 时间段时间段B向向A中渗透物质：中渗透物质：VA*CA(t+t)-VA*CA(t) = SKCB(t)-CA(t)t 推出推出dCA(t)/dt = SK/VA*CB(t)-CA(t) 两边除以两边除以t, 令令t0又由质量守恒考察整个容器中物质总量始终不变：又由质量守恒考察整个容器中物质总量始终不变：VA*CA(t)+VB*CB(t) = VA*aA+VB*aB 推出推出CA(t) = aA+VB/VA*aB-VB/VA*CB(t) 代入上式代入上式2 推出推出dCB(t)/dt = -SK(1/VA+1/VB)CB(t)+SK(aA/VB+aB/

13、VA)CB(0)=aB 初值条件初值条件 此带初值微分方程可由此带初值微分方程可由dsolve求解求解在上式中，已知的包括在上式中，已知的包括VA，VB，S以及一组以及一组t和和CB(t)值值未知的包括未知的包括aA，aB，K，下面通过数据拟合确定渗透率，下面通过数据拟合确定渗透率K第31页/共43页在上式中，代入已知值在上式中，代入已知值VA=VB=1000cm3，S=10cm2令令a=(aA*VA+aB*VB)/(VA+VB)，b=VA(aB-aA)/(VA+VB)简化之后的表达式为：简化之后的表达式为：CB(t)=a+b*exp(-0.02*k*t)编写被调编写被调M文件文件 tbp79

14、.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*ex

15、p(-0.002*t)第32页/共43页第33页/共43页同名升级同名升级换名升级换名升级第34页/共43页第35页/共43页下面还是以引例的采样数据为例，进行演示：下面还是以引例的采样数据为例，进行演示：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导入数据导入数据绘制散点绘制散点图图第36页/共43页进行拟合进行拟合拟合方法拟合方法结果和误差分析结果和误差分析 这里可供选择的这里可供选择的拟合类型拟合类型和和可选参数可选参数比较多比较多,包括多项式函数包括多项式函

16、数,指数函数指数函数,幂函数等幂函数等,如何确定如何确定最优的方案最优的方案？第37页/共43页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 meas

17、ures 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. 第38页/共43页分析数据分析数据 这里可以计算拟合函数在采样点的一阶二阶导函数值，并绘制相应的函数图像这里可以计算拟合函数在采样点的一阶二阶导