插值和拟合(讲义).ppt

,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,插值与拟合在数学建模中的应用,数学学院 袁栩,1,一维插值,2,二维插值,注：,Hermite,插值（略）,Runge,现象,用,MATLAB,作插值计算,一维插值函数：,yi=interp1(x,，,y,，,xi,，,method),插值方法,被插值点,插值节点,xi,处的插值结果,nearest,最邻近插值；,linear,线性插值,(,缺省）；,spline,三次样条插值；,cubic,立方插值；,缺省时 分段线性插值,注意：所有的插值方法都要求,x,是单调，并且,xi,不能够超过,x,的范围,例,1.,用三次样条插值选取,11,个基点计算插值,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),此例，可以看出插值函数得到的函数图像与原函数很接近。,例,2,：从,1,点,12,点,的,11,小时内，每隔,1,小时测量一次温度，测得的温度的数值依次为：,5,，,8,，,9,，,15,，,25,，,29,，,31,，,30,，,22,，,25,，,27,，,24,试估计每隔,1/10,小时的温度值,hours=1:12;,temps=5 8 9 15 25 29 31 30 22 25 27 24;,h=1:0.1:12;,t=interp1(hours,temps,h,spline);%(,直接输出数据将是很多的,),plot(hours,temps,+,h,t,b,hours,temps,r:)%,作图,xlabel(Hour),ylabel(Degrees Celsius),x,y,机翼下轮廓线,例,3.,已知飞机下轮廓线上数据如下，求,x,每改变,0.1,时的,y,值,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:15;,y1=interp1(x0,y0,x);,y2=interp1(x0,y0,x,spline);,subplot(2,1,1),plot(x0,y0,k+,x,y1,r),grid,title(piecewise linear),subplot(2,1,2),plot(x0,y0,k+,x,y2,r),grid,title(spline),要求,x0,y0,单调；,x,，,y,可取,为矩阵，或,x,取行向量，,y,取为列向量，,x,y,的值分别不能超出,x0,y,0,的范围,z=interp2(x0,y0,z0,x,y,method),被插值点,插值方法,用,MATLAB,作网格节点数据的插值,插值节点,被插值点的函数值,nearest,最邻近插值；,linear,双线性插值（缺省值）；,cubic,双三次插值；,缺省时 双线性插值,.,例,4.,测得平板表面,3,5,网格点处的温度分别为：,82 81 80 82 84 79 63 61 65 81 84 84 82 85 86,试作出平板表面的温度分布曲面,z,=,f,(,x,y,),的图形,输入以下命令：,x=1:5;,y=1:3;,temps=82 81 80 82 84;79 63 61 65 81;84 84 82 85 86;,mesh(x,y,temps),1.,先在三维坐标画出原始数据，画出粗糙的温度分布曲线图,.,2,以平滑数据,在,x,、,y,方向上每隔,0.2,个单位的地方进行插值,.,再输入以下命令,:,xi=1:0.2:5;,yi=1:0.2:3;,zi=interp2(x,y,temps,xi,yi,cubic);,mesh(xi,yi,zi),画出插值后的温度分布曲面图,.,通过此例对最近邻点插值、双线性插值方法和双三次插值方法的插值效果进行比较,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);,插值函数,griddata,格式为,:,cz,=griddata,（,x,，,y,，,z,，,cx,，,cy,，,method,）,用,MATLAB,作散点数据的插值计算,要求,cx,取行向量，,cy,取为列向量,被插值点,插值方法,插值节点,被插值点的函数值,nearest,最邻近插值,linear,双线性插值,cubic,双三次插值,v4-MATLAB,提供的插值方法,缺省时,双线性插值,例,6.,在某海域测得一些点,(,x,y,),处的水深,z,由下表给出，船的吃水深度为,5,英尺，在矩形区域（,75,，,200,）,（,-50,，,150,）里的哪些地方船要避免进入,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;,z=-4-8-6-8-6-8-8-9-9-8-8-9-4-9;,cx=75:0.5:200;,cy=-50:0.5:150;,cz=griddata(x,y,z,cx,cy,cubic);,4.,作出水深小于,5,的海域范围,即,z,=5,的等高线,.,.,.,3,作海底曲面图,meshz(cx,cy,cz),rotate3d,xlabel(X),ylabel(Y),zlabel(Z),figure(2),contour(cx,cy,cz,-5-5);grid,hold on,plot(x,y,+),xlabel(X),ylabel(Y),用,MATLAB,解拟合问题,1.,线性最小二乘拟合,2.,非线性最小二乘拟合,用,MATLAB,作线性最小二乘拟合,1.,作多项式,f,(,x,)=,a,1,x,m,+,a,m,x,+,a,m+1,拟合,可利用已有程序,:,a=polyfit(x,y,m),2.,对超定方程组,可得最小二乘意义下的解,，用,3.,多项式在,x,处的值,y,可用以下命令计算：,y=polyval,（,a,，,x,）,输出拟合多项式系数,a=,a,1,a,m,a,m+1,(,数组,),）,输入同长度,的数组,x,，,y,拟合多项,式次数,即要求 出二次多项式,:,中 的,使得,:,例,7.,对下面一组数据作二次多项式拟合,1,）输入以下命令,：,x=0:0.1:1;,y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2;,R=(x.2)x ones(11,1);,A=Ry,解法,1,用解超定方程的方法,2,）计算结果,：,=-9.8108 20.1293 -0.0317,1,）输入以下命令：,x=0:0.1:1;,y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2;,A=polyfit(x,y,2),z=polyval(A,x);,plot(x,y,k+,x,z,r),%,作出数据点和拟合曲线的图形,2,）计算结果：,=-9.8108 20.1293 -0.0317,解法,2,用多项式拟合的命令,1.,lsqcurvefit,已知,数据点,：,xdata=,（,xdata,1,xdata,2,xdata,n,）,ydata=,（,ydata,1,ydata,2,ydata,n,）,用,MATLAB,作非线性最小二乘拟合,MATLAB,提供了两个求非线性最小二乘拟合的函数：,lsqcurvefit,和,lsqnonlin,两个命令都要先建立,M,文件,fun.m,，在其中定义函数,f,(,x,),，但两者定义,f,(,x,),的方式是不同的,可参考例题,.,用以求含参量,x,（向量）的向量值函数,F(x,xdata)=,（,F,（,x,xdata,1,）,F,（,x,，,xdata,n,）,T,中的参变量,x(,向量,),使得,输入格式为,:,(1),x=,lsqcurvefit,(fun,x0,xdata,ydata);,(2),x=lsqcurvefit(fun,x0,xdata,ydata,options);,(3),x=lsqcurvefit(fun,x0,xdata,ydata,options,grad);,(4),x,options=lsqcurvefit(fun,x0,xdata,ydata,);,(5),x,options,funval=lsqcurvefit(fun,x0,xdata,ydata,);,(6),x,options,funval,Jacob=lsqcurvefit(fun,x0,xdata,ydata,);,fun,是一个事先建立的定义函数,F(x,xdata),的,M,文件,自变量为,x,和,xdata,说明,:,x=lsqcurvefit(fun,x0,xdata,ydata,options);,迭代初值,已知数据点,选项见无,约束优化,lsqnonlin,用以求含参量,x,（向量）的向量值函数,f,(,x,),=(,f,1,(,x,),f,2,(,x,),f,n,(,x,),T,中的参量,x,，使得,最小,其中,f,i,（,x,）,=,f,（,x,，,xdata,i,，,ydata,i,）,=,F,(,x,xdata,i,)-,ydata,i,2.,lsqnonlin,已知数据点：,xdata=,（,xdata,1,，,xdata,2,，,，,xdata,n,）,ydata=,（,ydata,1,，,ydata,2,，,，,ydata,n,）,输入格式为：,1),x=,lsqnonlin,（,fun,，,x0,）,；,2),x=lsqnonlin,（,fun,，,x0,，,options,）,；,3),x=lsqnonlin,（,fun,，,x0,，,optionsgrad,）,；,4),x,，,options=lsqnonlin,（,fun,，,x0,，,）,；,5),x,，,options,，,funval=lsqnonlin,（,funx0,，,）,；,说明：,x=,lsqnonlin,（,fun,，,x0,，,options,）,；,fun,是一个事先建立的定义函数,f,(,x,),的,M,文件，,自变量为,x,迭代初值,选项见无,约束优化,例,8.,用下面一组数据拟合,中的参数,a,，,b,，,k,该问题即解最优化问题：,1,）编写,M,文件,curvefun1.m,function,f=curvefun1(,x,tdata,),f=x(1)+x(2)*exp(-0.02*x(3)*tdata),%,其中,x(1)=a;x(2)=b,；,x(3)=k;,2,）输入命令,tdata=100:100:1000,cdata=1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;,x0=0.2,0.05,0.05;,x=lsqcurvefit(curvefun1,x0,tdata,cdata),f=curvefun1(x,tdata),F,(,x,，,tdata,)=,，,x,=(,a,，,b,，,k,),解法,1,.,用命令,lsqcurvefit,3,）运算结果为,：,f=0.0043 0.0051 0.0056 0.0059 0.0061,0.0062 0.0062 0.0063 0.0063 0.0063,x=0.0063 -0.0034 0.2542,4,）结论,：,a,=0.0063,b,=-0.0034,k,=0.2542,解法,2,用命令,lsqnonlin,f,(,x,)=,F,(x,tdata,ctada,)=,x,=,（,a,，,b,，,k,）,1,）,编写,M,文件,curvefun2.m,function f=curvefun2(x),tdata=100:100:1000;,cdata=1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;,f=cdata-x(1)-x(2)*exp(-0.02*x(3)*tdata),2,）输入命令,:,x0=0.2,0.05,0.05;,x=lsqnonlin(,curvefun2,x0),f=curvefun2(x),函数,curvefun2,的自变量是,x,，,cdata,和,tdata,是已知参数，故应将,cdata tdata,的值写在,curvefun2.m,中,3,）运算结果为,f,=1.0e-003*(0.2322 -0.1243 -0.2495 -0.2413-0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792,x=0.0063 -0.0034 0.2542,可以看出,两个命令的计算结果是相同的,.,4,）结论,：,即拟合得,a,=,0.0063,b,=-0.0034,k,=0.2542,MATLAB,解应用问题实例,1.,电阻问题,2.,给药方案问题,*3.,水塔流量估计问题,电阻问题,温度,t,(,C)20.5 32.7 51.0 73.0 95.7,电阻,R,(,)765 826 873 942 1032,例,.,由数据,拟合,R,=,a,1,t,+,b,方法,1.,用命令,polyfit(x,y,m),t=20.5 32.5 51 73 95.7;,r=765 826 873 942 1032;,aa=polyfit(t,r,1);,a=aa(1),b=aa(2),y=polyval(aa,t);,plot(t,r,k+,t,y,r),得到,a,=3.3940,b,=702.4918,方法,2.,直接用,结果相同,t=20.5 32.5 51 73 95.7;,r=765 826 873 942 1032;,R=t ones(5,1);,aa=Rr;,a=aa(1),b=aa(2),y=polyval(aa,t);,plot(t,r,k+,t,y,r),一室模型,：将整个机体看作一个房室，称,中心室,，室内血药浓度是均匀的快速静脉注射后，浓度立即上升；然后迅速下降当浓度太低时，达不到预期的治疗效果；当浓度太高，又可能导致药物中毒或副作用太强临床上，每种药物有一个最小有效浓度,c,1,和一个最大有效浓度,c,2,设计给药方案时，要使血药浓度 保持在,c,1,c,2,之间本题设,c,1,=10,ug/ml,，,c,2,=25ug/ml.,拟 合 问 题 实 例,2,给药方案,一种新药用于临床之前，必须设计给药方案,.,药物进入机体后通过血液输送到全身，在这个过程中不断地被吸收、分布、代谢，最终排出体外，药物在血液中的浓度，即单位体积血液中的药物含量，称为,血药浓度,在实验方面,对某人用快速静脉注射方式一次注入该药物,300mg,后,在一定时刻,t,(h),采集血药,测得血药浓度,c,(ug/ml),如下表,:,t,(h)0.25 0.5 1 1.5 2 3 4 6 8,c,(,g/ml),19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01,要设计给药方案,必须知道给药后血药浓度随时间变化的规律从实验和理论两方面着手：,给药方案,1.,在快速静脉注射的给药方式下，研究血药浓度（单位体积血液中的药物含量）的变化规律,t,c,2,c,c,1,O,问题,2.,给定药物的最小有效浓度和最大治疗浓度，设计给药方案：每次注射剂量多大；间隔时间多长,分析,理论：用一室模型研究血药浓度变化规律,实验：对血药浓度数据作拟合，符合负指数变化规律,3.,血液容积,v,t,=0,注射剂量,d,血药浓度立即为,d,/,v,.,2.,药物排除速率与血药浓度成正比，比例系数,k,(0),模型假设,1.,机体看作一个房室，室内血药浓度均匀,一室模型,模型建立,在此，,d,=300mg,，,t,及,c,（,t,）在某些点处的值见前表，需经拟合求出参数,k,、,v.,用线性最小二乘拟合,c,(,t,),计算结果：,d=300;,t=0.25 0.5 1 1.5 2 3 4 6 8;,c=19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01;,y=log(c);,a=polyfit(t,y,1),k=-a(1),v=d/exp(a(2),程序：,用非线性最小二乘拟合,c,(,t,)-,用,lsq,curvefit,2.,主程序如下,clear,tdata=0.25 0.5 1 1.5 2 3 4 6 8;,cdata=19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01;,x0=10,0.5;,x=lsqcurvefit(,curvefun3,x0,tdata,cdata);,f=curvefun3(x,tdata),x,1.,用,M,文件,curvefun3.m,定义函数,function f=curvefun3(x,tdata),d=300,f=(x(1)d)*exp(-x(2)*tdata),%x(1)=v;x(2)=k,给药方案 设计,c,c,2,c,1,O,t,设每次注射剂量,D,间隔时间,血药浓度,c,(,t,),应,c,1,c,(,t,),c,2,初次剂量,D,0,应加大,给药方案记为：,2.,1.,计算结果：,给药方案：,c,1,=10,c,2,=25,k=0.2347,v=15.02,故可制定给药方案：,即,:,首次注射,375mg,，,其余每次注射,225mg,，,注射的间隔时间为,4h,估计水塔的流量,2.,解题思路,3.,算法设计与编程,1.,问题,某居民区有一供居民用水的圆柱形水塔，一般可以通过测量其水位来估计水的流量，但面临的困难是，当水塔水位下降到设定的最低水位时，水泵自动启动向水塔供水，到设定的最高水位时停止供水，这段时间无法测量水塔的水位和水泵的供水量通常水泵每天供水一两次，每次约两小时,.,水塔是一个高,12.2m,，直径,17.4m,的正圆柱按照设计，水塔水位降至约,8.2m,时，水泵自动启动，水位升到约,10.8m,时水泵停止工作,表,1,是某一天的水位测量记录，试估计任何时刻（包括水泵正供水时）从水塔流出的水流量，及一天的总用水量,流量估计的解题思路,拟合水位,时间函数,确定流量,时间函数,估计一天总用水量,拟合水位,时间函数,从测量记录看，一天有两个供水时段（以下称第,1,供水时段和第,2,供水时段），和,3,个水泵不工作时段（以下称第,1,时段,t,=0,到,t,=8.97,，第,2,次时段,t,=10.95,到,t,=20.84,和第,3,时段,t,=23,以后）对第,1,、,2,时段的测量数据直接分别作多项式拟合，得到水位函数为使拟合曲线比较光滑，多项式次数不要太高，一般在,36,由于第,3,时段只有,3,个测量记录，无法对这一时段的水位作出较好的拟合,确定流量,时间函数,对于第,1,、,2,时段只需将水位函数求导数即可，对于两个供水时段的流量，则用供水时段前后（水泵不工作时段）的流量拟合得到，并且将拟合得到的第,2,供水时段流量外推，将第,3,时段流量包含在第,2,供水时段内,一天总用水量的估计,总用水量等于两个水泵不工作时段和两个供水时段用水量之和，它们都可以由流量对时间的积分得到,算法设计与编程,1.,拟合第,1,、,2,时段的水位，并导出流量,2.,拟合供水时段的流量,3.,估计一天总用水量,4.,流量及总用水量的检验,1.,拟合第,1,时段的水位，并导出流量,设,t,，,h,为已输入的时刻和水位测量记录（水泵启动的,4,个时刻不输入），,第,1,时段,各时刻的流量可如下得：,1,）,c1=polyfit,（,t,（,1,：,10,），,h,（,1,：,10,），,3,）；,%,用,3,次多项式拟合第,1,时段水位，,c1,输出,3,次多项式的系数,2,）,a1=polyder,（,c1,）；,%a1,输出多项式（系数为,c1,）导数的系数,3,）,tp1=0,：,0.1,：,9,；,x1=-polyval,（,a1,，,tp1,）；,%x1,输出多项式（系数,a1,）在,tp1,点的函数值（取负后边为正值），即,tp1,时刻的流量,4,）,流量函数为：,拟合第,2,时段的水位，并导出流量,设,t,，,h,为已输入的时刻和水位测量记录（水泵启动的,4,个时刻不输入），,第,2,时段,各时刻的流量可如下得：,1,）,c2=polyfit(t(10.9:21),h(10.9:21),3),；,%,用,3,次多项式拟合第,2,时段水位，,c2,输出,3,次多项式的系数,2,）,a2=polyder(c2),；,%,a2,输出多项式（系数为,c2,）导数的系数,3,）,tp2=10.9:0.1:21;,x2=-polyval(a2,tp2);,%,x2,输出多项式（系数为,a2,）在,tp2,点的函数值（取负后边为正值），即,tp2,时刻的流量,4,）,流量函数为：,2.,拟合供水时段的流量,在第,1,供水时段（,t,=911,）之前（即第,1,时段）和之后（即第,2,时段）各取几点，其流量已经得到，用它们拟合第,1,供水时段的流量为使流量函数在,t,=9,和,t,=11,连续，我们简单地只取,4,个点，拟合,3,次多项式（即曲线必过这,4,个点），实现如下：,xx1=-polyval(a1,8 9);,%,取第,1,时段在,t,=8,9,的流量,xx2=-polyval(a2,11 12);,%,取第,2,时段在,t=11,12,的流量,xx12=xx1 xx2;,c12=polyfit(8 9 11 12,xx12,，,3);,%,拟合,3,次多项式,tp12=9,：,0.1,：,11,;,x12=polyval(c12,tp12);,%,x12,输出第,1,供水时段 各时刻的流量,拟合的流量函数为：,在第,2,供水时段之前取,t,=20,，,20.8,两点的流水量，在该时刻之后（第,3,时段）仅有,3,个水位记录，我们用差分得到流量，然后用这,4,个数值拟合第,2,供水时段的流量如下：,dt3=diff,（,t(22,：,24),）；,%,最后,3,个时刻的两两之差,dh3=diff,（,h(22,：,24),）；,%,最后,3,个水位的两两之差,dht3=-dh3./dt3,；,%t(22),和,t,(23),的流量,t3=20 20.8 t(22)t(23),；,xx3=-polyval(a2,，,t3(1,：,2),，,dht3),；,%,取,t3,各时刻的流量,c3=polyfit,（,t3,，,xx3,，,3,）；,%,拟合,3,次多项式,t3=20.8,：,0.1,：,24,；,x3=polyval,（,c3,，,tp3,）；,%x3,输出第,2,供水时段,（外推至,t,=24,）各时刻的流量,拟合的流量函数为：,3.,一天总用水量的估计,第,1,、,2,时段和第,1,、,2,供水时段流量的积分之和，就是一天总用水量虽然诸时段的流量已表为多项式函数，积分可以解析地算出，这里仍用数值积分计算如下：,y1=0.1*trapz(x1),；,%,第,1,时段用水量（仍按高度计），,0.1,为积分步长,y2=0.1*trapz(x2),；,%,第,2,时段用水量,y12=0.1*trapz(x12),；,%,第,1,供水时段用水量,y3=0.1*trapz(x3),；,%,第,2,供水时段用水量,y=(y1+y2+y12+y3)*237.8*0.01;,%,一天总用水量（）,计算结果：,y1,=146.2,y2,=266.8,y12,=47.4,y3,=77.3,，,y,=1250.4,4.,流量及总用水量的检验,计算出的,各时刻的流量,可用水位记录的数值微分来检验用水量,y1,可用第,1,时段水位测量记录中下降高度,968-822=146,来检验，类似地，,y2,用,1082-822=260,检验,供水时段流量,的一种,检验方法,如下：供水时段的用水量加上水位上升值,260,是该时段泵入的水量，除以时段长度得到水泵的功率（单位时间泵入的水量），而两个供水时段水泵的功率应大致相等第,1,、,2,时段水泵的功率可计算如下：,p1=(y12+260)/2,；,%,第,1,供水时段水泵功率（水量仍以高度计）,tp4=20.8,：,0.1,：,23,；,xp2=polyval,（,c3,，,tp4,）；,%xp2,输出第,2,供水时段各时刻流量,p2=(0.1*trapz(xp2)+260)/2.2,；,%,第,2,供水时段水泵的功率（水量仍以高度计）,计算结果,：,p1=154.5,，,p2=140.1,计算结果,流量函数为：,流量曲线见图,n,=(3,4),n,=(5,6),练习,1.,用给定的多项式，如,y,=,x,3,-6,x,2,+5,x,-3,，产生一组数据,(,x,i,y,i,，,i,=1,2,n,),再在,y,i,上添加随机干扰,(,可用,rand,产生,(0,1),均匀分布随机数,或用,rands,产生,N,(0,1),分布随机数,),，然后用,x,i,和添加了随机干扰的,y,i,作的,3,次多项式拟合，与原系数比较 如果作,2,或,4,次多项式拟合，结果如何？,实验作业,练习,2.,用电压,V,=10V,的电池给电容器充电，电容器上,t,时刻的电压为 ，其中,V,0,是电容器的初始电压,是充电常数,.,试由下面一组,t,，,v,数据确定,V,0,.,t(s),0.5,1,2,3,4,5,7,9,V(,伏）,6.36,6.48,7.26,8.22,8.66,8.99,9.43,9.63,实验作业,实验作业,