matlab多项式、插值与数据拟合.ppt

第五章多项式、插值与数据拟合,多项式MATLAB命令插值Lagrange插值Hermite插值Runge现象和分段插值分段插值样条插值的MATLAB表示数据拟合多项式拟合函数线性组合的曲线拟合方法最小二乘曲线拟合B样条函数及其MATLAB表示,5.1关于多项式MATLAB命令,一个多项式的幂级数形式可表示为：也可表为嵌套形式或因子形式N阶多项式n个根，其中包含重根和复根。若多项式所有系数均为实数，则全部复根都将以共轭对的形式出现,幂系数：在MATLAB里，多项式用行向量表示，其元素为多项式的系数，并从左至右按降幂排列。例：被表示为p=2145poly2sym(p)ans=2*x3+x2+4*x+5Roots:多项式的零点可用命令roots求的。例：r=roots(p)得到r=0.2500+1.5612i0.2500-1.5612i-1.0000所有零点由一个列向量给出。,Poly:由零点可得原始多项式的各系数，但可能相差一个常数倍。例：poly(r)ans=1.00000.50002.00002.5000注意：若存在重根，这种转换可能会降低精度。例：r=roots(1-615-2015-61)r=1.0042+0.0025i1.0042-0.0025i1.0000+0.0049i1.0000-0.0049i0.9958+0.0024i0.9958-0.0024i舍入误差的影响，与计算精度有关。,polyval:可用命令polyval计算多项式的值。例：计算y(2.5)c=3,-7,2,1,1；xi=2.5;yi=polyval(c,xi)yi=23.8125如果xi是含有多个横坐标值的数组，则yi也为与xi长度相同的向量。c=3,-7,2,1,1;xi=2.5,3;yi=polyval(c,xi)yi=23.812576.0000,polyfit:给定n+1个点将可以唯一确定一个n阶多项式。利用命令polyfit可容易确定多项式的系数。例：x=1.1,2.3,3.9,5.1;y=3.887,4.276,4.651,2.117;a=polyfit(x,y,length(x)-1)a=-0.20151.4385-2.74775.4370poly2sym(a)ans=-403/2000*x3+2877/2000*x2-27477/10000*x+5437/1000多项式为Polyfit的第三个参数是多项式的阶数。,多项式积分：功能：求多项式积分调用格式：py=poly_itg(p)p:被积多项式的系数py：求积后多项式的系数poly_itg.mfunctionpy=poly_itg(p)n=length(p);py=p.*n:-1:1.(-1),0不包括最后一项积分常数,多项式微分：Polyder:求多项式一阶导数的系数。调用格式为:b=polyder(c)c为多项式y的系数，b是微分后的系数，其值为：,两个多项式的和与差：命令poly_add:求两个多项式的和,其调用格式为:c=poly_add(a,b)多项式a减去b，可表示为:c=poly_add(a,-b),功能：两个多项式相加调用格式：b=poly_add(p1,p2)b:求和后的系数数组poly_add.mfunctionp3=poly_add(p1,p2)n1=length(p1);n2=length(p2);ifn1=n2p3=p1+p2;endifn1n2p3=p1+zeros(1,n1-n2),p2;endifn1a=2,-5,6,-1,9;b=3,-90,-18;c=conv(a,b)c=6-195432-4539-792-162q,r=deconv(c,b)q=2-56-19r=0000000poly2sym(c)ans=6*x6-195*x5+432*x4-453*x3+9*x2-792*x-162,5.2插值5.2.1Lagrange插值,方法介绍对给定的n个插值点及对应的函数值，利用构造的n-1次Lagrange插值多项式，则对插值区间内任意x的函数值y可通过下式求的：MATLAB实现,functiony=lagrange(x0,y0,x)ii=1:length(x0);y=zeros(size(x);fori=iiij=find(ii=i);y1=1;forj=1:length(ij),y1=y1.*(x-x0(ij(j);endy=y+y1*y0(i)/prod(x0(i)-x0(ij);end算例：给出f(x)=ln(x)的数值表，用Lagrange计算ln(0.54)的近似值。x=0.4:0.1:0.8；y=-0.916291,-0.693147,-0.510826,-0.356675,-0.223144;lagrange(x,y,0.54,0.55,0.78)ans=-0.6161-0.5978-0.2484（精确解-0.616143),5.2.2Hermite插值,方法介绍不少实际问题不但要求在节点上函数值相等，而且要求导数值也相等，甚至要求高阶导数值也相等，满足这一要求的插值多项式就是Hermite插值多项式。下面只讨论函数值与一阶导数值个数相等且已知的情况。已知n个插值点及对应的函数值和一阶导数值。则对插值区间内任意x的函数值y的Hermite插值公式：,MATLAB实现%hermite.mfunctiony=hermite(x0,y0,y1,x)n=length(x0);m=length(x);fork=1:myy=0.0;fori=1:nh=1.0;a=0.0;forj=1:nifj=ih=h*(x(k)-x0(j)/(x0(i)-x0(j)2;a=1/(x0(i)-x0(j)+a;endendyy=yy+h*(x0(i)-x(k)*(2*a*y0(i)-y1(i)+y0(i);endy(k)=yy;end,算例:对给定数据,试构造Hermite多项式求出sin0.34的近似值。x0=0.3,0.32,0.35;y0=0.29552,0.31457,0.34290;y1=0.95534,0.94924,0.93937;y=hermite(x0,y0,y1,0.34)y=0.3335sin(0.34)与精确值比较ans=0.3335,x=0.3:0.005:0.35;y=hermite(x0,y0,y1,x);plot(x,y)y2=sin(x);holdonplot(x,y2,-r),5.2.3Runge现象,问题的提出：根据区间a,b上给出的节点做插值多项式p(x)的近似值，一般总认为p(x)的次数越高则逼近f(x)的精度就越好，但事实并非如此。反例：在区间-5,5上的各阶导数存在，但在此区间上取n个节点所构成的Lagrange插值多项式在全区间内并非都收敛。取n=10,用Lagrange插值法进行插值计算。,x=-5:1:5;y=1./(1+x.2);x0=-5:0.1:5;y0=lagrange(x,y,x0);y1=1./(1+x0.2);绘制图形plot(x0,y0,-r)插值曲线holdonplot(x0,y1,-b)原曲线为解决Rung问题，引入分段插值。,算法分析：所谓分段插值就是通过插值点用折线或低次曲线连接起来逼近原曲线。MATLAB实现可调用内部函数。命令1interp1功能：一维数据插值（表格查找）。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。格式1yi=interp1(x,Y,xi)%返回插值向量yi，每一元素对应于参量xi，同时由向量x与Y的内插值决定。参量x指定数据Y的点。若Y为一矩阵，则按Y的每列计算。算例对于t，beta、alpha分别有两组数据与之对应，用分段线性插值法计算当t=321,440,571时beta、alpha的值。,5.2.4分段插值,temp=300,400,500,600;beta=1000*3.33,2.50,2.00,1.67;alpha=10000*0.2128,0.3605,0.5324,0.7190;ti=321,400,571;propty=interp1(temp,beta,alpha,ti);propty=interp1(temp,beta,alpha,ti,linear);ti,proptyans=1.0e+003*0.32103.15572.43820.40002.50003.60500.57101.76576.6489,格式2yi=interp1(Y,xi)%假定x=1:N，其中N为向量Y的长度，或者为矩阵Y的行数。格式3yi=interp1(x,Y,xi,method)%用指定的算法计算插值：nearest：最近邻点插值，直接完成计算；linear：线性插值（缺省方式），直接完成计算；spline：三次样条函数插值。cubic：分段三次Hermite插值。其它，如v5cubic。对于超出x范围的xi的分量，使用方法nearest、linear、v5cubic的插值算法，相应地将返回NaN。对其他的方法，interp1将对超出的分量执行外插值算法。yi=interp1(x,Y,xi,method,extrap)yi=interp1(x,Y,xi,method,extrapval)%确定超出x范围的xi中的分量的外插值extrapval，其值通常取NaN或0。,算例year=1900:10:2010;product=75.995,91.972,105.711,123.203,131.669,.150.697,179.323,203.212,226.505,249.633,256.344,267.893;p1995=interp1(year,product,1995)p1995=252.9885x=1900:1:2010;y=interp1(year,product,x,cubic);plot(year,product,o,x,y),例:已知的数据点来自函数根据生成的数据进行插值处理，得出较平滑的曲线直接生成数据。x=0:.12:1;y=(x.2-3*x+5).*exp(-5*x).*sin(x);plot(x,y,x,y,o),调用interp1()函数：x1=0:.02:1;y0=(x1.2-3*x1+5).*exp(-5*x1).*sin(x1);y1=interp1(x,y,x1);y2=interp1(x,y,x1,cubic);y3=interp1(x,y,x1,spline);y4=interp1(x,y,x1,nearest);plot(x1,y1,y2,y3,y4,:,x,y,o,x1,y0)误差分析max(abs(y0(1:49)-y2(1:49),max(abs(y0-y3),max(abs(y0-y4)ans=0.01770.00860.1598,x0=-1+2*0:10/10;y0=1./(1+25*x0.2);x=-1:.01:1;y=lagrange(x0,y0,x);%Lagrange插值ya=1./(1+25*x.2);plot(x,ya,x,y,:),例,y1=interp1(x0,y0,x,cubic);y2=interp1(x0,y0,x,spline);plot(x,ya,x,y1,:,x,y2,-),命令2interp2功能二维数据内插值（表格查找）格式1ZI=interp2(X,Y,Z,XI,YI)%返回矩阵ZI，其元素包含对应于参量XI与YI（可以是向量、或同型矩阵）的元素。参量X与Y必须是单调的，且相同的划分格式，就像由命令meshgrid生成的一样。若Xi与Yi中有在X与Y范围之外的点，则相应地返回NaN。格式2ZI=interp2(Z,XI,YI)%缺省地，X=1:n、Y=1:m，其中m,n=size(Z)。再按第一种情形进行计算。格式3ZI=interp2(X,Y,Z,XI,YI,method)%用指定的算法method计算二维插值：linear：双线性插值算法（缺省算法）；nearest：最临近插值；spline：三次样条插值；cubic：双三次插值。,算例：years=1950:10:1990;service=10:10:30;wage=150.697199.592187.625179.323195.072250.287203.212179.092322.767226.505153.706426.730249.633120.281598.243;w=interp2(service,years,wage,15,1975)w=190.6288,例x,y=meshgrid(-3:.6:3,-2:.4:2);z=(x.2-2*x).*exp(-x.2-y.2-x.*y);surf(x,y,z),axis(-3,3,-2,2,-0.7,1.5),选较密的插值点，用默认的线性插值算法进行插值x1,y1=meshgrid(-3:.2:3,-2:.2:2);z1=interp2(x,y,z,x1,y1);surf(x1,y1,z1),axis(-3,3,-2,2,-0.7,1.5),立方和样条插值：z1=interp2(x,y,z,x1,y1,cubic);z2=interp2(x,y,z,x1,y1,spline);surf(x1,y1,z1),axis(-3,3,-2,2,-0.7,1.5)figure;surf(x1,y1,z2),axis(-3,3,-2,2,-0.7,1.5),算法误差的比较z=(x1.2-2*x1).*exp(-x1.2-y1.2-x1.*y1);surf(x1,y1,abs(z-z1),axis(-3,3,-2,2,0,0.08)figure;surf(x1,y1,abs(z-z2),axis(-3,3,-2,2,0,0.025),二维一般分布数据的插值,功能:可对非网格数据进行插值格式:z=griddata(x0,y0,z0,x,y,method)v4:MATLAB4.0提供的插值算法，公认效果较好；linear：双线性插值算法（缺省算法）；nearest：最临近插值；spline：三次样条插值；cubic：双三次插值。例：在x为-3，3，y为2，2矩形区域随机选择一组坐标，用v4与cubic插值法进行处理，并对误差进行比较。,x=-3+6*rand(200,1);y=-2+4*rand(200,1);z=(x.2-2*x).*exp(-x.2-y.2-x.*y);x1,y1=meshgrid(-3:.2:3,-2:.2:2);z1=griddata(x,y,z,x1,y1,cubic);surf(x1,y1,z1),axis(-3,3,-2,2,-0.7,1.5)z2=griddata(x,y,z,x1,y1,v4);figure;surf(x1,y1,z2),axis(-3,3,-2,2,-0.7,1.5),误差分析z0=(x1.2-2*x1).*exp(-x1.2-y1.2-x1.*y1);surf(x1,y1,abs(z0-z1),axis(-3,3,-2,2,0,0.15)figure;surf(x1,y1,abs(z0-z2),axis(-3,3,-2,2,0,0.15),例：在x为3，3，y为2，2矩形区域随机选择一组坐标中，对分布不均匀数据，进行插值分析。x=-3+6*rand(200,1);y=-2+4*rand(200,1);z=(x.2-2*x).*exp(-x.2-y.2-x.*y);%生成已知数据plot(x,y,x)%样本点的二维分布figure,plot3(x,y,z,x),axis(-3,3,-2,2,-0.7,1.5)，grid,去除在(-1,-1/2)点为圆心，以0.5为半径的圆内的点。x=-3+6*rand(200,1);y=-2+4*rand(200,1);%重新生成样本点z=(x.2-2*x).*exp(-x.2-y.2-x.*y);ii=find(x+1).2+(y+0.5).20.52);%找出满足条件的点坐标x=x(ii);y=y(ii);z=z(ii);plot(x,y,x)t=0:.1:2*pi,2*pi;x0=-1+0.5*cos(t);y0=-0.5+0.5*sin(t);line(x0,y0)%在图形上叠印该圆，可见，圆内样本点均已剔除,用新样本点拟合出曲面x1,y1=meshgrid(-3:.2:3,-2:.2:2);z1=griddata(x,y,z,x1,y1,v4);surf(x1,y1,z1),axis(-3,3,-2,2,-0.7,1.5),误差分析z0=(x1.2-2*x1).*exp(-x1.2-y1.2-x1.*y1);surf(x1,y1,abs(z0-z1),axis(-3,3,-2,2,0,0.1)contour(x1,y1,abs(z0-z1),30);holdon,plot(x,y,x);line(x0,y0)误差的二维等高线图,命令3interp3三维网格生成用meshgrid()函数,调用格式:x,y,z=meshgrid(x1,y1,z1)其中x1,y1,z1为这三维所需要的分割形式，应以向量形式给出，返回x,y,z为网格的数据生成，均为三维数组。griddata3()三维非网格形式的插值拟合命令4interpnn维网格生成用ndgrid()函数,调用格式:x1,x2,xn=ndgridv1,v2,vngriddatan()n维非网格形式的插值拟合interp3()、interpn()调用格式同interp2()函数一致；griddata3()、griddatan()调用格式同griddata()函数一致。,例：通过函数生成一些网格型样本点，试根据样本点进行拟合，并给出拟合误差。x,y,z=meshgrid(-1:0.2:1);x0,y0,z0=meshgrid(-1:0.05:1);V=exp(x.2.*z+y.2.*x+z.2.*y).*cos(x.2.*y.*z+z.2.*y.*x);V0=exp(x0.2.*z0+y0.2.*x0+z0.2.*y0).*cos(x0.2.*y0.*z0+z0.2.*y0.*x0);V1=interp3(x,y,z,V,x0,y0,z0,spline);err=V1-V0;max(err(:)ans=0.0419,5.2.5样条插值的MATLAB表示,定义一个三次样条函数类：S=csapi(x,y)其中x=x1,x2,.,xn,y=y1,y2,yn为样本点。S返回样条函数对象的插值结果，包括子区间点、各区间点三次多项式系数等。可用fnplt()绘制出插值结果，其调用格式：fnplt(S)对给定的向量xp，可用fnval()函数计算，其调用格式：yp=fnval(S,xp)其中得出的yp是xp上各点的插值结果。,例：x0=0,0.4,1,2,pi;y0=sin(x0);sp=csapi(x0,y0),fnplt(sp,:);holdon,sp=form:ppbreaks:00.4000123.1416coefs:4x4doublepieces:4order:4dim:1ezplot(sin(t),0,pi);plot(x0,y0,o)sp.coefsans=-0.16270.00760.99650-0.1627-0.18760.92450.38940.0244-0.48040.52380.84150.0244-0.4071-0.36370.9093,在(0.4000,1)区间内，插值多项式可以表示为：,例,点，用三次样条插值的方法对这些数据进行拟合,x=0:.12:1;y=(x.2-3*x+5).*exp(-5*x).*sin(x);sp=csapi(x,y);fnplt(sp)c=sp.breaks(1:4)sp.breaks(2:5)sp.coefs(1:4,:),.sp.breaks(5:8)sp.breaks(6:9)sp.coefs(5:8,:)c=Columns1through700.120024.7396-19.35884.515100.48000.12000.240024.7396-10.45260.93770.30580.60000.24000.36004.5071-1.5463-0.50220.31050.72000.36000.48001.91390.0762-0.67860.23580.8400,Columns8through120.6000-0.24040.7652-0.57760.15880.7200-0.47740.6787-0.40430.10010.8400-0.45590.5068-0.26210.06050.9600-0.45590.3427-0.16010.0356,格式S=csapi(x1,x2,xn,z),处理多个自变量的网格数据三次样条插值类：,x0=-3:.6:3;y0=-2:.4:2;x,y=ndgrid(x0,y0);%注意这里只能用ndgrid,否则生成的z矩阵顺序有问题z=(x.2-2*x).*exp(-x.2-y.2-x.*y);sp=csapi(x0,y0,z);fnplt(sp);,例,函数spline功能三次样条数据插值格式yy=spline(x,y,xx)例:对离散分布在y=exp(x)sin(x)函数曲线上的数据点进行样条插值计算：x=024581212.817.219.920;y=exp(x).*sin(x);xx=0:.25:20;yy=spline(x,y,xx);plot(x,y,o,xx,yy),5.3数据拟合,用插值的方法对一函数进行近似,要求所得到的插值多项式经过已知插值节点;在n比较大的情况下,插值多项式往往是高次多项式,这也就容易出现振荡现象（龙格现象），即虽然在插值节点上没有误差,但在插值节点之外插值误差变得很大,从“整体”上看,插值逼近效果将变得“很差”。所谓数据拟合是求一个简单的函数,例如是一个低次多项式,不要求通过已知的这些点,而是要求在整体上“尽量好”的逼近原函数。这时,在每个已知点上就会有误差,数据拟合就是从整体上使误差,尽量的小一些。,5.3.1多项式拟合,n次多项式：曲线与数据点的残差为：残差的平方和为：为使其最小化，可令R关于的偏导数为零，即：,或或矩阵形式：,多项式拟合MATLAB命令：polyfit格式：p=polyfit(x,y,n),x0=0:.1:1;y0=(x0.2-3*x0+5).*exp(-5*x0).*sin(x0);p3=polyfit(x0,y0,3);vpa(poly2sym(p3),10)%可以如下显示多项式ans=2.839962923*x3-4.789842696*x2+1.943211631*x+.5975248921e-1,例,绘制拟合曲线：x=0:.01:1;ya=(x.2-3*x+5).*exp(-5*x).*sin(x);y1=polyval(p3,x);plot(x,y1,x,ya,x0,y0,o),就不同的次数进行拟合：p4=polyfit(x0,y0,4);y2=polyval(p4,x);p5=polyfit(x0,y0,5);y3=polyval(p5,x);p8=polyfit(x0,y0,8);y4=polyval(p8,x);plot(x,ya,x0,y0,o,x,y2,x,y3,x,y4),拟合最高次数为8的多项式：vpa(poly2sym(p8),5)ans=-8.2586*x8+43.566*x7-101.98*x6+140.22*x5-125.29*x4+74.450*x3-27.672*x2+4.9869*x+.42037e-6Taylor幂级数展开：symsx;y=(x2-3*x+5)*exp(-5*x)*sin(x);vpa(taylor(y,9),5)ans=5.*x-28.*x2+77.667*x3-142.*x4+192.17*x5-204.96*x6+179.13*x7-131.67*x8多项式表示数据模型是不唯一的，即是两个多项式函数完全不同。在某一区域内其曲线可能特别近似。,多项式拟合的效果并不一定总是很精确的。x0=-1+2*0:10/10;y0=1./(1+25*x0.2);x=-1:.01:1;ya=1./(1+25*x.2);p3=polyfit(x0,y0,3);y1=polyval(p3,x);p5=polyfit(x0,y0,5);y2=polyval(p5,x);p8=polyfit(x0,y0,8);y3=polyval(p8,x);p10=polyfit(x0,y0,10);y4=polyval(p10,x);plot(x,ya,x,y1,x,y2,-.,x,y3,-,x,y4,:),例,用Taylor幂级数展开效果将更差。symsx;y=1/(1+25*x2);p=taylor(y,x,10)p=1-25*x2+625*x4-15625*x6+390625*x8多项式拟合效果x1=-1:0.01:1;ya=1./(1+25*x1.2);y1=subs(p,x,x1);plot(x1,ya,-,x1,y1),5.3.2函数线性组合的曲线拟合方法,该方程的最小二乘解为:,其中,例,x=0,0.2,0.4,0.7,0.9,0.92,0.99,1.2,1.4,1.48,1.5;y=2.88;2.2576;1.9683;1.9258;2.0862;2.109;2.1979;2.5409;2.9627;3.155;3.2052;A=ones(size(x)，exp(-3*x),cos(-2*x).*exp(-4*x)，x.2;c=Ay;c1=cc1=1.22002.3397-0.67970.8700,图形显示x0=0:0.01:1.5;A1=ones(size(x0)exp(-3*x0),cos(-2*x0).*exp(-4*x0)x0.2;y1=A1*c;plot(x0,y1,x,y,x),数据分析x=1.1052,1.2214,1.3499,1.4918,1.6487,1.8221,2.0138,.2.2255,2.4596,2.7183,3.6693;y=0.6795,0.6006,0.5309,0.4693,0.4148,0.3666,0.3241,.0.2864,0.2532,0.2238,0.1546;plot(x,y,x,y,*),例,分别对x,y进行对数变换：x1=log(x);y1=log(y);plot(x1,y1),A=x1,ones(size(x1);c=Ay1c=-1.2339-0.2630exp(c(2)ans=0.7687,x=0:0.1:1;y=(x.2-3*x+5).*exp(-5*x).*sin(x);n=8;A=;fori=1:n+1,A(:,i)=x.(n+1-i);endc=Ay;vpa(poly2sym(c),5)ans=-8.2586*x8+43.566*x7-101.98*x6+140.22*x5-125.29*x4+74.450*x3-27.672*x2+4.9869*x+.42037e-6,例,5.3.3最小二乘曲线拟合,格式：a,jm=lsq