《在科学计算中的应用5》由会员分享，可在线阅读，更多相关《在科学计算中的应用5（90页珍藏版）》请在金锄头文库上搜索。
1、第五章 多项式、插值与数据拟合多项式MATLAB命令插值Lagrange插值Hermite插值Runge现象和分段插值分段插值样条插值的MATLAB表示数据拟合多项式拟合函数线性组合的曲线拟合方法最小二乘曲线拟合5.1 关于多项式MATLAB命令一个多项式的幂级数形式可表示为：也可表为嵌套形式或因子形式 N阶多项式n个根，其中包含重根和复根。若多项式所有系数均为实数，则全部复根都将以共轭对的形式出现 幂系数：在MATLAB里，多项式用行向量表示，其元素为多项式的系数，并从左至右按降幂排列。 例： 被表示为 p=2 1 4 5 poly2sym(p) ans = 2*x3+x2+4*x+5Roo
2、ts: 多项式的零点可用命令roots求的。 例： r=roots(p) 得到 r = 0.2500 + 1.5612i 0.2500 - 1.5612i -1.0000 所有零点由一个列向量给出。Poly: 由零点可得原始多项式的各系数，但可能相差一个常数倍。 例： poly(r)ans = 1.0000 0.5000 2.0000 2.5000 注意：若存在重根，这种转换可能会降低精度。 例： r=roots(1 -6 15 -20 15 -6 1) r = 1.0042 + 0.0025i 1.0042 - 0.0025i 1.0000 + 0.0049i 1.0000 - 0.0049
3、i 0.9958 + 0.0024i 0.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.8125 76.0000polyfit:给定n+1个点将可以唯一确定一个n阶多项式。利用命令polyfit可容易确定多项式的系数。 例： x=1.1,
4、2.3,3.9,5.1; y=3.887,4.276,4.651,2.117; a=polyfit(x,y,length(x)-1)a = -0.2015 1.4385 -2.7477 5.4370 poly2sym(a) ans = -403/2000*x3+2877/2000*x2-27477/10000*x+5437/1000 多项式为Polyfit的第三个参数是多项式的阶数。多项式积分： 功能：求多项式积分 调用格式：py=poly_itg(p) p:被积多项式的系数 py：求积后多项式的系数 poly_itg.m function py=poly_itg(p) n=length(p)
5、; 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.mfunction p3=poly_add(p1,p2)n1=length(p1);n2=length(p2
6、);if n1=n2 p3=p1+p2;endif n1n2 p3=p1+zeros(1,n1-n2),p2;endif n1 a=2,-5,6,-1,9; b=3,-90,-18; c=conv(a,b)c = 6 -195 432 -453 9 -792 -162 q,r=deconv(c,b)q = 2 -5 6 -1 9r = 0 0 0 0 0 0 0 poly2sym(c) ans = 6*x6-195*x5+432*x4-453*x3+9*x2-792*x-162 5.2 插值5.2.1 Lagrange插值方法介绍 对给定的n个插值点 及对应的函数值 ，利用构造的n-1次Lag
7、range插值多项式，则对插值区间内任意x的函数值y可通过下式求的：MATLAB实现function y=lagrange(x0,y0,x)ii=1:length(x0); y=zeros(size(x);for i=ii ij=find(ii=i); y1=1; for j=1:length(ij), y1=y1.*(x-x0(ij(j); end y=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.5108
8、26,-0.356675,-0.223144; lagrange(x,y,0.54,0.55,0.78)ans = -0.6161 -0.5978 -0.2484 （ 精确解-0.616143)5.2.2 Hermite插值方法介绍 不少实际问题不但要求在节点上函数值相等，而且要求导数值也相等，甚至要求高阶导数值也相等，满足这一要求的插值多项式就是Hermite插值多项式。下面只讨论函数值与一阶导数值个数相等且已知的情况。 已知n个插值点 及对应的函数值 和一阶导数值 。则对插值区间内任意x的函数值y的Hermite插值公式：MATLAB实现% hermite.mfunction y=herm
9、ite(x0,y0,y1,x)n=length(x0); m=length(x);for k=1:m yy=0.0; for i=1:n h=1.0; a=0.0; for j=1:n if j=i h=h*(x(k)-x0(j)/(x0(i)-x0(j)2; a=1/(x0(i)-x0(j)+a; end end yy=yy+h*(x0(i)-x(k)*(2*a*y0(i)-y1(i)+y0(i); end y(k)=yy;end算例:对给定数据,试构造Hermite多项式求出sin0.34的近似值。 x0=0.3,0.32,0.35; y0=0.29552,0.31457,0.34290;
10、 y1=0.95534,0.94924,0.93937; format long; y=hermite(x0,y0,y1,0.34)y = 0.33348889007407 sin(0.34) 与精确值比较ans = 0.33348709214081 x=0.3:0.005:0.35;y=hermite(x0,y0,y1,x); plot(x,y) y2=sin(x); hold on plot(x,y2,-r)5.2.3 Runge现象问题的提出：根据区间a,b上给出的节点做插值多项式p(x)的近似值，一般总认为p(x)的次数越高则逼近f(x)的精度就越好，但事实并非如此。反例： 在区间-5
11、,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)插值曲线 hold on plot(x0,y1,-b)原曲线为解决Rung问题，引入分段插值。算法分析：所谓分段插值就是通过插值点用折线或低次曲线连接起来逼近原曲线。MATLAB实现 可调用内部函数。命令1 interp1 功能 ： 一维数据插值（表格查找）。该命
12、令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。 格式1 yi = interp1(x,Y,xi) %返回插值向量yi，每一元素对应于参量xi，同时由向量x与Y的内插值决定。参量x指定数据Y的点。若Y为一矩阵，则按Y的每列计算。5.2.4 分段插值格式2 yi = interp1(Y,xi) %假定x=1:N，其中N为向量Y的长度，或者为矩阵Y的行数。格式3 yi = interp1(x,Y,xi,method) %用指定的算法计算插值： nearest：最近邻点插值，直接完成计算； linear：线性插值（缺省方式），直接完成计算； splin
13、e：三次样条函数插值。 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,1
14、23.203,131.669,. 150.697,179.323,203.212,226.505,249.633,256.344,267.893; p1995 = interp1(year,product,1995)p1995 = 252.9885 x = 1900:1:2010; y = interp1(year,product,x,cubic); plot(year,product,o,x,y)例:已知的数据点来自函数已知的数据点来自函数根据生成的数据进行插值处理，得出较平滑的曲线根据生成的数据进行插值处理，得出较平滑的曲线直接生成数据。直接生成数据。 x=0:.12:1; y=(x.2-
15、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(y
16、0-y4)ans = 0.0177 0.0086 0.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,-)命令2 interp2 功能 二维数据内插值 格式1 ZI = interp2(X,Y,Z,XI,YI) %返回矩阵ZI，其元素包含对应于
17、参量XI与YI（可以是向量、或同型矩阵）的元素。参量X与Y必须是单调的，且相同的划分格式，就像由命令meshgrid生成的一样。若Xi与Yi中有在X与Y范围之外的点，则相应地返回NaN。格式2 ZI = interp2(Z,XI,YI) %缺省地，X=1:n、Y=1:m，其中m,n=size(Z)。再按第一种情形进行计算。格式3 ZI = interp2(X,Y,Z,XI,YI,method) %用指定的算法method计算二维插值： linear：双线性插值算法（缺省算法）； nearest：最临近插值； spline：三次样条插值； cubic：双三次插值。算例： years=1950:1
18、0:1990; service=10:10:30; wage = 150.697 199.592 187.625 179.323 195.072 250.287 203.212 179.092 322.767 226.505 153.706 426.730 249.633 120.281 598.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,-
19、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
20、=(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
21、，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)误差分析
22、 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, plot
23、3(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(x
24、0,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); hold on, plot(x,y,x); line(x0,y0) 误差的
25、二维等高线图命令3 interp3 三维网格生成用meshgrid( )函数,调用格式: x,y,z=meshgrid(x1,y1,z1) 其中x1,y1,z1为这三维所需要的分割形式，应以向量形式给出，返回x,y,z为网格的数据生成，均为三维数组。 griddata3( ) 三维非网格形式的插值拟合命令4 interpn n维网格生成用ndgrid( )函数,调用格式: x1,x2,xn=ndgridv1,v2,vn griddatan( ) n维非网格形式的插值拟合interp3 ( )、 interpn( )调用格式同interp2( )函数一致；函数一致；griddata3( )、 g
26、riddatan( )调用格式同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
27、; max(err(:)ans = 0.04195.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,:); hold on,
28、sp = form: pp breaks: 0 0.4000 1 2 3.1416 coefs: 4x4 double pieces: 4 order: 4 dim: 1 ezplot(sin(t),0,pi); plot(x0,y0,o) sp.coefsans = -0.1627 0.0076 0.9965 0 -0.1627 -0.1876 0.9245 0.3894 0.0244 -0.4804 0.5238 0.8415 0.0244 -0.4071 -0.3637 0.9093在(0.4000, 1)区间内，插值多项式可以表示为：例点，用三次样条插值的方法对这些数据进行拟合 x=0
29、:.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 = Columns 1 through 7 0 0.1200 24.7396 -19.3588 4.5151 0 0.4800 0.1200 0.2400 24.7396 -10.4526 0.9377 0.3058 0.6000 0.2400 0.3600 4.507
30、1 -1.5463 -0.5022 0.3105 0.7200 0.3600 0.4800 1.9139 0.0762 -0.6786 0.2358 0.8400 Columns 8 through 12 0.6000 -0.2404 0.7652 -0.5776 0.1588 0.7200 -0.4774 0.6787 -0.4043 0.1001 0.8400 -0.4559 0.5068 -0.2621 0.0605 0.9600 -0.4559 0.3427 -0.1601 0.0356格式 S=csapi(x1,x2,xn,z)处理多个自变量的网格数据三次样条插值类： x0=-3:.
31、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 = 0 2 4 5 8 12 12.8 17.2 19.9 20; y = exp(x).*sin(x); xx = 0:.25:20;
32、yy = spline(x,y,xx); plot(x,y,o,xx,yy)格式5.2.6 B样条函数及其MATLAB表示 S=spapi(k,x,y)例例 x0=0,0.4,1,2,pi; y0=sin(x0); ezplot(sin(t),0,pi); hold on sp1=csapi(x0,y0); fnplt(sp1,-); % 三次分段多项式样条插值 sp2=spapi(5,x0,y0); fnplt(sp2,:) % 5 次 B 样条插值y=sin(t)和 x=0:.12:1; y=(x.2-3*x+5).*exp(-5*x).*sin(x); ezplot(x2-3*x+5)*
33、exp(-5*x)*sin(x),0,1), hold on sp1=csapi(x,y); fnplt(sp1,-); sp2=spapi(5,x,y); fnplt(sp2,:)5.2.7 基于样条插值的数值微积分运算基于样条插值的数值微分运算格式： Sd=fnder(S,k)该函数可以求取S的k阶导数。格式： Sd=fnder(S,k1,kn)可以求取多变量函数的偏导数例： syms x; f=(x2-3*x+5)*exp(-5*x)*sin(x); ezplot(diff(f),0,1), hold on x=0:.12:1; y=(x.2-3*x+5).*exp(-5*x).*sin
34、(x); sp1=csapi(x,y);建立三次样条函数 dsp1=fnder(sp1,1); fnplt(dsp1,-)绘制样条图 sp2=spapi(5,x,y);5阶次B样条 dsp2=fnder(sp2,1); fnplt(dsp2,:); axis(0,1,-0.8,5)例：拟合曲面 x0=-3:.3:3; y0=-2:.2:2; x,y=ndgrid(x0,y0); z=(x.2-2*x).*exp(-x.2-y.2-x.*y); sp=spapi(5,5,x0,y0,z); B样条dspxy=fnder(sp,1,1); fnplt(dspxy)%生成样条图理论方法： syms
35、x y; z=(x2-2*x)*exp(-x2-y2-x*y); ezsurf(diff(diff(z,x),y),-3 3,-2 2)对符号变量表达式做三维表面图基于样条插值的数值积分运算格式： f=fnint(S)其中S为样条函数。例：考虑 中较稀疏的样本点,用样条积分的方式求出定积分及积分函数。 x=0,0.4,1 2,pi; y=sin(x); sp1=csapi(x,y); a=fnint(sp1,1); 建立三次样条函数并积分 xx=fnval(a,0,pi); xx(2)-xx(1)ans = 2.0191 sp2=spapi(5,x,y); b=fnint(sp2,1); xx
36、=fnval(b,0,pi); xx(2)-xx(1)ans = 1.9999绘制曲线 ezplot(-cos(t)+2,0,pi); hold on不定积分可上下平移 fnplt(a,-); fnplt(b,:)5.3 数据拟合用插值的方法对一函数进行近似用插值的方法对一函数进行近似,要求所得到要求所得到的插值多项式经过已知插值节点的插值多项式经过已知插值节点;在在n比较大的比较大的情况下情况下,插值多项式往往是高次多项式插值多项式往往是高次多项式,这也就这也就容易出现振荡现象（龙格现象），即虽然在插容易出现振荡现象（龙格现象），即虽然在插值节点上没有误差值节点上没有误差,但在插值节点之外插
37、值误但在插值节点之外插值误差变得很大差变得很大,从从“整体整体”上看上看,插值逼近效果将插值逼近效果将变得变得“很差很差”。所所谓谓数数据据拟拟合合是是求求一一个个简简单单的的函函数数,例例如如是是一一个个低低次次多多项项式式,不不要要求求通通过过已已知知的的这这些些点点,而而是是要要求求在在整整体体上上“尽尽量量好好”的的逼逼近近原原函函数数。这这时时,在在每每个个已已知知点点上上就就会会有有误误差差,数数据据拟拟合合就就是是从从整体上使误差整体上使误差,尽量的小一些。尽量的小一些。5.3.1 多项式拟合n次多项式：曲线与数据点 的残差为：残差的平方和为：为使其最小化，可令R关于 的偏导数为
38、零，即：或 或矩阵形式：多项式拟合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,
39、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+.420
40、37e-6Taylor幂级数展开： syms x; 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);
41、 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幂级数展开效果将更差。 syms x; 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./
42、(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.2200 2.3397 -0.6797 0
43、.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,
44、y进行对数变换： x1=log(x); y1=log(y); plot(x1,y1) A=x1, ones(size(x1); c=Ay1c = -1.2339 -0.2630 exp(c(2)ans = 0.7687 x=0:0.1:1; y=(x.2-3*x+5).*exp(-5*x).*sin(x); n=8; A=; for i=1:n+1, A(:,i)=x.(n+1-i); end c=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
45、.9869*x+.42037e-6例例5.3.3 最小二乘曲线拟合格式： a, jm=lsqcurvefit(Fun,a0,x,y)例 x=0:.1:10; y=0.12*exp(-0.213*x)+0.54*exp(-0.17*x).*sin(1.23*x); f=inline(a(1)*exp(-a(2)*x)+a(3)*exp(-a(4)*x).*sin(a(5)*x),a,x); xx,res=lsqcurvefit(f,1,1,1,1,1,x,y); xx,resOptimization terminated successfully: Relative function value
46、 changing by less than OPTIONS.TolFunans = 0.1197 0.2125 0.5404 0.1702 1.2300res = 7.1637e-007修改最优化选项： ff=optimset; ff.TolFun=1e-20; ff.TolX=1e-15; % 修改精度限制 xx,res=lsqcurvefit(f,1,1,1,1,1,x,y,ff); xx,res 变量界Optimization terminated successfully: Relative function value changing by less than OPTIONS.T
47、olFunans = 0.1200 0.2130 0.5400 0.1700 1.2300res = 9.5035e-021绘制曲线： x1=0:0.01:10; y1=f(xx,x1); plot(x1,y1,x,y,o)例例 x=0.1:0.1:1; y=2.3201,2.6470,2.9707,3.2885,3.6008,3.9090,4.2147,4.5191,4.8232,5.1275;function y=c8f3(a,x)y=a(1)*x+a(2)*x.2.*exp(-a(3)*x)+a(4); a=lsqcurvefit(c8f3,1;2;2;3,x,y); aMaximum number of function evaluations exceeded; increase options.MaxFunEvalsans = 2.4575 2.4557 1.4437 2.0720绘制曲线： y1=c8f3(a,x); plot(x,y,x,y1,o)
