数值计算04-插值与拟合.ppt

1,数值计算,第四章插值与拟合,我们经常会遇到大量的数据需要处理，而处理数据的关键就在于这些算法，例如数据拟合、参数估计、插值等数据处理算法。此类问题在MATLAB中有很多现成的函数可以调用，熟悉MATLAB，这些方法都能游刃有余的用好。,在实际中，常常要处理由实验或测量所得到的一些离散数据。,例：从1点12点的11小时内，每隔1小时测量一次温度，测得的温度的数值依次为：5，8，9，15，25，29，31，30，22，25，27，24试估计每隔1/10小时的温度值,%interp_1hours=1:12;temps=589152529313022252724;plot(hours,temps,rs),在实际中，常常要处理由实验或测量所得到的一些离散数据。,例：从1点12点的11小时内，每隔1小时测量一次温度，测得的温度的数值依次为：5，8，9，15，25，29，31，30，22，25，27，24试估计每隔1/10小时的温度值,在实际中，常常要处理由实验或测量所得到的一些离散数据。,例：从1点12点的11小时内，每隔1小时测量一次温度，测得的温度的数值依次为：5，8，9，15，25，29，31，30，22，25，27，24试估计每隔1/10小时的温度值,在实际中，常常要处理由实验或测量所得到的一些离散数据。,例：从1点12点的11小时内，每隔1小时测量一次温度，测得的温度的数值依次为：5，8，9，15，25，29，31，30，22，25，27，24试估计每隔1/10小时的温度值,插值与拟合方法就是要通过这些数据去确定某一类已知函数的参数或寻求某个近似函数，使所得到的近似函数与已知数据有较高的拟合精度。如果要求这个近似函数（曲线或曲面）经过所已知的所有数据点，则称此类问题为插值问题。（不需要函数表达式）,如果不要求近似函数通过所有数据点，而是要求它能较好地反映数据变化规律的近似函数的方法称为数据拟合。（必须有函数表达式）近似函数不一定（曲线或曲面）通过所有的数据点。,9,第四章插值与拟合,4.1插值的含义4.2插值的Matlab实现4.3拟合问题,10,4.1插值的含义,一维插值:已知函数在n+1个节点的值,求任一点的函数值,11,4.1插值的含义,一维插值:已知函数在n+1个节点的值,求任一点的函数值,二维插值:已知函数在个节点的值,求任一点的函数值,一维插值的定义,2.2插值的Matlab实现,基本格式：yc=interp1(x,y,cx,method)%x,y分别表示已知数据点的横、纵坐标向量，x必须单调；%cx为需要插值的横坐标数据%method为插值方法，有nearest最近邻点插值linear线性插值（默认）spline三次样条插值cubic三次插值注：Lagrange,Newton插值法需自编程序。,线性插值主要参数,15,例2.9：对函数进行11个点的三次样条插值：x=-5,-4,-3,-2,-1,0,1,2,3,4,5;y=1/26,1/17,0.1,0.2,0.5,1,0.5,0.2,0.1,1/17,1/26;xi=-5:0.01:5;yi=interp1(x,y,xi);plot(xi,yi,x,y,o)holdonf=inline(1/(x2+1);fplot(f,-5,5,r),例2.9：对函数进行11个点的三次样条插值：x=-5,-4,-3,-2,-1,0,1,2,3,4,5;y=1/26,1/17,0.1,0.2,0.5,1,0.5,0.2,0.1,1/17,1/26;xi=-5:0.01:5;yi=interp1(x,y,xi,spline);plot(xi,yi,x,y,o)holdonf=inline(1/(x2+1);fplot(f,-5,5,r),例2.9：对函数进行11个点的三次样条插值：x=-5,-4,-3,-2,-1,0,1,2,3,4,5;y=1/26,1/17,0.1,0.2,0.5,1,0.5,0.2,0.1,1/17,1/26;xi=-5:0.01:5;yi=interp1(x,y,xi,cubic);plot(xi,yi,x,y,o)holdonf=inline(1/(x2+1);fplot(f,-5,5,r),二维插值的定义,第一种（网格节点）：,已知mn个节点,第二种（散乱节点）：,注意：最邻近插值一般不连续。具有连续性的最简单的插值是分片线性插值。,最邻近插值,将四个插值点（矩形的四个顶点）处的函数值依次简记为：,分片线性插值,f(xi,yj)=f1，f(xi+1,yj)=f2，f(xi+1,yj+1)=f3，f(xi,yj+1)=f4,插值函数为：,第二片(上三角形区域)：(x,y)满足,插值函数为：,注意：(x,y)当然应该是在插值节点所形成的矩形区域内。显然，分片线性插值函数是连续的；,分两片的函数表达式如下：,第一片(下三角形区域)：(x,y)满足,双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形式如下：,其中有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数。,双线性插值,要求x0,y0单调；x，y可取为矩阵，或x取行向量，y取为列向量，x,y的值分别不能超出x0,y0的范围。,z=interp2(x0,y0,z0,x,y,method),用MATLAB作网格节点数据的插值,nearest最邻近插值linear双线性插值cubic双三次插值缺省时,双线性插值,例2.10：测得平板表面3*5网格点处的温度分别为：828180828479636165818484828586试作出平板表面的温度分布曲面z=f(x,y)的图形。,输入以下命令：x=1:5;y=1:3;temps=8281808284;7963616581;8484828586;mesh(x,y,temps),1.先在三维坐标画出原始数据，画出粗糙的温度分布曲图.,%输入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi,yi,cubic);mesh(xi,yi,zi),2以平滑数据,在x、y方向上每隔0.2个单位的地方进行插值,画出插值后的温度分布曲面图.,32,已知某处山区地形选点测量坐标数据为：x=00.511.522.533.544.55y=00.511.522.533.544.555.56海拔高度数据为：z=8990878592919693908782929698999591898684828496989592908885848381858081828995969392898686828587989996978885828382858994959392918684888892939495898786838192929697989693958482818485858182808081859093958486819899989796958487808185828384879095868880828184858683828180828788899899979698949287,例2.11,画出该山区地形地貌图.,33,插值前插值后,Tomountain.m,最近邻点插值,双线性插值,双三次插值,插值函数griddata格式为:,cz=griddata（x，y，z，cx，cy，method）,用MATLAB作散点数据的插值计算,要求cx取行向量，cy取为列向量,nearest最邻近插值linear双线性插值cubic双三次插值v4-MATLAB提供的插值方法缺省时,双线性插值,例2.13:在某海域测得一些点(x,y)处的水深z由下表给出，船的吃水深度为5英尺，在矩形区域（75，200）（-50，150）里的哪些地方船要避免进入,4.作出水深小于5的海域范围,即z=5的等高线.,2.在矩形区域(75,200)(-50,150)进行插值。,1.输入插值基点数据,3.作海底曲面图,步骤：,%程序一：插值并作海底曲面图x=129.0140.0103.588.0185.5195.0105.5157.5107.577.081.0162.0162.0117.5;y=7.5141.523.0147.022.5137.585.5-6.5-813.056.5-66.584.0-33.5;z=48686889988949;x1=75:1:200;y1=-50:1:150;x1,y1=meshgrid(x1,y1);z1=griddata(x,y,z,x1,y1,v4);meshc(x1,y1,z1),海底曲面图,%程序二：插值并作出水深小于5的海域范围。x1=75:1:200;y1=-50:1:150;x1,y1=meshgrid(x1,y1);z1=griddata(x,y,z,x1,y1,v4);%插值z1(z1=5)=nan;%将水深大于5的置为nan，这样绘图就不会显示出来meshc(x1,y1,z1),水深小于5的海域范围,在实际应用中，经常遇到下列数据处理问题：,已知函数在m个点上的数据表，寻求其近似函数。,2.3曲线拟合问题,所谓”曲线拟合”,是指根据给定的数据表,寻找一个简单的表达式来”拟合”该组数据,此处的”拟合”的含义为:不要求该表达式对应的近似曲线完全通过所有的数据点,只要求该近似曲线能够反映数据的基本变化趋势.,已知函数在m个点上的数据表，寻求其近似函数。,设的近似函数为,其中是某函数族中的已知线性无关函数。,(1)最小二乘多项式拟合,称为残差向量,寻求一组常数，要求,的2-范数达到最小。,则得到最小二乘问题：,上述问题的解也称为方程组的最小二乘解。,当时称之为超定（或矛盾）方程组。,例2.14：考察某种纤维的强度与其拉伸倍数的关系.下表是实际测定的24个纤维样品的强度与相应的拉伸倍数的数据记录:,可以看出,纤维强度随拉伸倍数增加而增加,并且24个点大致分布在一条直线附近,该直线称为这一问题的数学模型。,因此可认为强度与拉伸倍数之间的主要关系是线性关系,怎样确定a,b,使得直线能较好地反映所给数据的基本“变化趋势”?,采用最小二乘的思想,令,问题转化为求参数使达到最小值。,这种求线性函数y=a+bx的过程称为线性拟合。,(2)非线性拟合,某些非线性拟合问题可转化为线性拟合问题,线性化处理：,令,则,由线性拟合方法可得到和，从而得到和。,如何选取拟合函数？,1.通过机理分析建立数学模型来确定f(x)；,2.将数据(xi,yi)i=1,n作图，通过直观判断确定f(x)：,用MATLAB解拟合问题,1、线性最小二乘拟合,2、非线性最小二乘拟合,用MATLAB作线性最小二乘拟合,1.作多项式f(x)=a1xm+amx+am+1拟合,可利用已有程序:,a=polyfit(x,y,m),2.对超定方程组,3.多项式在x处的值y可用以下命令计算：y=polyval（a，x）,例2.15对下面一组数据作二次多项式拟合,1）输入以下命令：x=0:0.1:1;y=-04471.9783.286.167.087.347.669.569.489.3011.2;R=(x.2)xones(10,1)；A=Ry,解法1用解超定方程的方法,2）计算结果：=-9.810820.1293-0.0317,1）输入以下命令：x=0:0.1:1;y=-0.4471.9783.286.167.087.347.669.569.489.3011.2;A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,k+,x,z,r)%作出数据点和拟合曲线的图形,2）计算结果：=-9.810820.1293-0.0317,解法2用多项式拟合的命令,1.lsqcurvefit已知数据点：xdata=（xdata1，xdata2，xdatan），ydata=（ydata1，ydata2，ydatan）,用MATLAB作非线性最小二乘拟合,Matlab的提供了两个求非线性最小二乘拟合的函数：lsqcurvefit和lsqnonlin。两个命令都要先建立M-文件fun.m，在其中定义函数f(x)，但两者定义f(x)的方式是不同的.,lsqcurvefit用以求含参量x（向量）的向量值函数F(x,xdata)=（F（x，xdata1），F（x，xdatan）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,);,说明：x=lsqcurvefit(fun,x0,xdata,ydata,options);,lsqnonlin用以求含参量x（向量）的向量值函数f(x)=(f1(x),f2(x),fn(x)T中的参量x，使得最小。其中fi（x）=f（x，xdatai，ydatai）=F(x,xdatai)-ydatai,2.lsqnonlin,已知数据点：xdata=（xdata1，xdata2，xdatan）ydata=（ydata1，ydata2，ydatan）,输入格式为：1）x=lsqnonlin（fun，x0）；2）x=lsqnonlin（fun，x0，options）；3）x=lsqnonlin（fun，x0，options，grad）；4）x，options=lsqnonlin（fun，x0，）；5）x，options，funval=lsqnonlin（fun，x0，）；,说明：x=lsqnonlin（fun，x0，options）；,例2.16用下面一组数据拟合中的参数a，b，k,该问题即解最优化问题：,1）编写M-文件curvefun1.mfunctionf=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:1000cdata=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.00430.00510.00560.00590.00610.00620.00620.00630.00630.0063x=0.0063-0.00340.2542,4）结论：a=0.0063,b=-0.0034,k=0.2542,解法2用命令lsqnonlinf(x)=F(x,tdata,ctada)=x=（a，b，k）,1）编写M-文件curvefun2.mfunctionf=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=x(1)+x(2)*exp(-0.02*x(3)*tdata)-cdata,2）输入命令:x0=0.2,0.05,0.05;x=lsqnonlin(curvefun2,x0)f=curvefun2(x),函数curvefun2的自变量是x，cdata和tdata是已知参数，故应将cdatatdata的值写在curvefun2.m中,3）运算结果为f=1.0e-003*(0.2322-0.1243-0.2495-0.2413-0.1668-0.07240.02410.11590.20300.2792x=0.0063-0.00340.2542,可以看出,两个命令的计算结果是相同的.,4）结论：即拟合得a=0.0063b=-0.0034k=0.2542,如何预报人口的增长人口的增长是当前世界上引起普遍关注的问题，并且我们会发现在不同的刊物预报同一时间的人口数字不相同，这显然是由于用了不同的人口模型计算的结果。我国是世界第一人口大国，基本上地球每九个人中就有一个中国人。有效地控制我国人口的增长是使我过全面进入小康社会、到21世纪中叶建成富强民主文明的社会主义国家的需要。而有效控制人口增长的前提是要认识人口数量的变化规律，建立人口模型，作出较准确的预报。,例:如何预报人口的增长,例2.17：1949年1994年我国人口数据资料如下：年份xi1949195419591964196919741979198419891994人口数yi5.46.06.77.08.19.19.810.311.311.8建模分析我国人口增长的规律,预报1999年我国人口数。,模型一：假设人口随时间线性地增加,模型：参数估计观测值的模型：拟合的精度:误差平方和。可以算出：a=-283.2320b=0.1480模型：y=1.93+0.146x,则可看成是线性方程,用polyfit命令计算得：,模型二：指数增长模型,可变为,Y,A,=,+,BX,a=2.33,b=0.0179,则所求模型为:,%程序nihe_2.m如下：x=1949195419591964196919741979198419891994;y=5.46.06.77.08.19.19.810.311.311.8;a=polyfit(x,y,1);x1=1949:10:1994;y1=a(2)+a(1)*x1;b=polyfit(x,log(y),1);y2=exp(b(2)*exp(b(1)*x1);plot(x,y,*)holdonplot(x1,y1,-r)holdonplot(x1,y2,-k)legend(原曲线,模型一曲线,模型二曲线),77,例3：体重约70kg的某人在短时间内喝下2瓶啤酒后，隔一定时间测量他的血液中酒精含量（mg/100ml），得到数据如表所示。试用所给数据用函数进行拟合，并求出未知常数。,78,分析：,由于拟合函数形式已经确定，且不是多项式函数，但若取对数可得,这样对参数是线性的。,因此可考虑先对数据h进行对数变换，再调用lsqcurvefit()函数进行拟合。,79,clearall;closeall;t=0.250.50.7511.522.533.544.55678910111213141516;h=3068758282776868585150413835282518151210774;h1=log(h);%对数变换f=inline(a(1)+a(2).*log(t)+a(3).*t,a,t);x,r=lsqcurvefit(f,1,0.5,-0.5,t,h1)%求参数lna,b,c的拟合值,x=4.48340.4709-0.2663r=0.4097,输出结果：,80,近年来我国的电信事业发展迅速，现已成世界第一电信大国。据统计某市的在过去近年中通信工具的拥有量如下表（单位为万台）：,应用案例1：通信工具的发展趋势,问题的提出：,使用matlab语句plot(x0,y0,r.,Markersize,10),81,应用案例1：通信工具的发展趋势,82,应用案例1：通信工具的发展趋势,O,x,lny,246810,.,.,.,.,.,.,.,.,.,1.00.80.60.40.2,83,clearall;closeall;x=123456789;y=1.782.242.743.744.455.316.928.8510.97;plot(x,y,r.,Markersize,10)holdonx1=1:9;b=polyfit(x,log(y),1);y2=exp(b(2)*exp(b(1)*x1);plot(x1,y2,-k),在实际问题中，有时散点图作出后未必是多项式的图