第五章 多项式、插值与数据拟合 多项式MATLAB命令 插值 数据拟合 Lagrange插值 Hermite插值 Runge现象和分段插值 第五章 多项式、插值与数据拟合 多项式MATLAB命令 插值 Lagrange插值 Hermite插值 Runge现象和分段插值 分段插值 样条插值的MATLAB表示 数据拟合 多项式拟合 函数线性组合的曲线拟合方法 最小二乘曲线拟合 B样条函数及其MATLAB表示

5.1 关于多项式MATLAB命令 一个多项式的幂级数形式可表示为： 也可表为嵌套形式 或因子形式 N阶多项式n个根，其中包含重根和复根。若多项式所有系数均为实数，则全部复根都将以共轭对的形式出现

幂系数：在MATLAB里，多项式用行向量表示，其元素为多项式的系数，并从左至右按降幂排列。 例： 被表示为 >> p=[2 1 4 5] >> poly2sym(p) ans = 2*x^3+x^2+4*x+5 Roots: 多项式的零点可用命令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.0049i 0.9958 + 0.0024i 0.9958 - 0.0024i 舍入误差的影响，与计算精度有关。

polyval: 可用命令polyval计算多项式的值。 >> 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) 23.8125 76.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.2015 1.4385 -2.7477 5.4370 >> poly2sym(a) ans = -403/2000*x^3+2877/2000*x^2-27477/10000*x+5437/1000 多项式为 Polyfit的第三个参数是多项式的阶数。

多项式积分： 功能：求多项式积分 调用格式：py=poly_itg(p) p:被积多项式的系数 py：求积后多项式的系数 poly_itg.m function py=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)

function p3=poly_add(p1,p2) n1=length(p1); n2=length(p2); 功能：两个多项式相加 调用格式：b=poly_add(p1,p2) b:求和后的系数数组 poly_add.m function p3=poly_add(p1,p2) n1=length(p1); n2=length(p2); if n1==n2 p3=p1+p2;end if n1>n2 p3=p1+[zeros(1,n1-n2),p2];end if n1<n2 p3=[zeros(1,n2-n1),p1]+p2;end

m阶多项式与n阶多项式的乘积是d=m+n阶的多项式: 计算 系数的MATLAB命令是:c=conv(a,b) 多项式 除多项式 的除法满足: 其中 是商, 是除法的余数。多项式 和 可由命令deconv算出。 例：[q, r]=deconv(a,b)

例 >> 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 9 r = 0 0 0 0 0 0 0 >> poly2sym(c) ans = 6*x^6-195*x^5+432*x^4-453*x^3+9*x^2-792*x-162

5.2 插值 5.2.1 Lagrange插值 方法介绍 对给定的n个插值点 及对应的函数值 ，利用构造的n-1次Lagrange插值多项式，则对插值区间内任意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.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.2 Hermite插值 方法介绍 不少实际问题不但要求在节点上函数值相等，而且要求导数值也相等，甚至要求高阶导数值也相等，满足这一要求的插值多项式就是Hermite插值多项式。下面只讨论函数值与一阶导数值个数相等且已知的情况。 已知n个插值点 及对应的函数值 和一阶导数值 。则对插值区间内任意x的函数值y的Hermite插值公式：

MATLAB实现 % hermite.m function y=hermite(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 yy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i)); y(k)=yy;

算例:对给定数据,试构造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.3335 >> sin(0.34) ％与精确值比较 ans =

>> 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,5]上的各阶导数存在，但在此区间上取n个节点所构成的Lagrange插值多项式在全区间内并非都收敛。 取n=10,用Lagrange插值法进行插值计算。

为解决Rung问题，引入分段插值。 >> 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问题，引入分段插值。

5.2.4 分段插值 算法分析：所谓分段插值就是通过插值点用折线或低次曲线连接起来逼近原曲线。 MATLAB实现 可调用内部函数。 命令1 interp1 功能 ： 一维数据插值（表格查找）。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。 格式1 yi = interp1(x,Y,xi) %返回插值向量yi，每一元素对应于参量xi，同时由向量x与Y的内插值决定。参量x指定数据Y的点。若Y为一矩阵，则按Y的每列计算。 算例 对于t，beta 、alpha分别有两组数据与之对应，用分段线性插值法计算当t=321, 440, 571时beta 、alpha的值。

>> 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,propty] ans = 1.0e+003 * 0.3210 3.1557 2.4382 0.4000 2.5000 3.6050 0.5710 1.7657 6.6489

格式2 yi = interp1(Y,xi) %假定x=1:N，其中N为向量Y的长度，或者为矩阵Y的行数。 格式3 yi = 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.9885 >> x = 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.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) 功能 二维数据内插值（表格查找） 格式1 ZI = interp2(X,Y,Z,XI,YI) %返回矩阵ZI，其元素包含对应于参量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:10: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,-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’) 二维一般分布数据的插值 功能:可对非网格数据进行插值 格式: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，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).^2>0.5^2); % 找出满足条件的点坐标 >> 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); hold on, plot(x,y,‘x’); line(x0,y0) ％误差的二维等高线图

命令3 interp3 三维网格生成用meshgrid( )函数,调用格式: [x,y,z]=meshgrid(x1,y1,z1) 其中x1,y1,z1为这三维所需要的分割形式，应以向量形式给出，返回x,y,z为网格的数据生成，均为三维数组。 griddata3( ) 三维非网格形式的插值拟合 命令4 interpn n维网格生成用ndgrid( )函数,调用格式: [x1,x2,…,xn]=ndgrid[v1,v2,…,vn] griddatan( ) 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,':'); hold on, 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.coefs ans = -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:.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.5071 -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:.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; >>yy = spline(x,y,xx); >>plot(x,y,'o',xx,yy)

5.2.6 基于样条插值的数值微积分运算 基于样条插值的数值微分运算 格式： Sd=fnder(S,k) 该函数可以求取S的k阶导数。 Sd=fnder(S,[k1,…,kn]) 可以求取多变量函数的偏导数

例： >> syms x; f=(x^2-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(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 x y; z=(x^2-2*x)*exp(-x^2-y^2-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=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比较大的情况下,插值多项式往往是高次多项式,这也就容易出现振荡现象（龙格现象），即虽然在插值节点上没有误差,但在插值节点之外插值误差变得很大,从“整体”上看,插值逼近效果将变得“很差”。 所谓数据拟合是求一个简单的函数,例如是一个低次多项式,不要求通过已知的这些点,而是要求在整体上“尽量好”的逼近原函数。这时,在每个已知点上就会有误差,数据拟合就是从整体上使误差,尽量的小一些。

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*x^3-4.789842696*x^2+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); >> plot(x,ya,x0,y0,'o',x,y2,x,y3,x,y4)

拟合最高次数为8的多项式： >> vpa(poly2sym(p8),5) ans = -8.2586*x^8+43.566*x^7-101.98*x^6+140.22*x^5-125.29*x^4+74.450*x^3-27.672*x^2+4.9869*x+.42037e-6 Taylor幂级数展开： >> syms x; y=(x^2-3*x+5)*exp(-5*x)*sin(x); >> vpa(taylor(y,9),5) 5.*x-28.*x^2+77.667*x^3-142.*x^4+192.17*x^5-204.96*x^6+179.13*x^7-131.67*x^8 多项式表示数据模型是不唯一的，即是两个多项式函数完全不同。在某一区域内其曲线可能特别近似。

例 多项式拟合的效果并不一定总是很精确的。 >> 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幂级数展开效果将更差。 多项式拟合效果 >> syms x; y=1/(1+25*x^2); p=taylor(y,x,10) p = 1-25*x^2+625*x^4-15625*x^6+390625*x^8 多项式拟合效果 >> 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=A\y; c1=c' c1 = 1.2200 2.3397 -0.6797 0.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=[A\y1']‘ -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=A\y; vpa(poly2sym(c),5) ans = -8.2586*x^8+43.566*x^7-101.98*x^6+140.22*x^5-125.29*x^4+74.450*x^3-27.672*x^2+4.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',res Optimization terminated successfully: Relative function value changing by less than OPTIONS.TolFun ans = 0.1197 0.2125 0.5404 0.1702 1.2300 res = 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.TolFun ans = 0.1200 0.2130 0.5400 0.1700 1.2300 res = 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); a' Maximum number of function evaluations exceeded; increase options.MaxFunEvals ans = 2.4575 2.4557 1.4437 2.0720

绘制曲线： >> y1=c8f3(a,x); plot(x,y,x,y1,’o’)

5.3.4 B样条函数及其MATLAB表示 格式 S=spapi(k,x,y)

例 y=sin(t)和 >> 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 样条插值

>> x=0:.12:1; y=(x.^2-3*x+5).*exp(-5*x).*sin(x); >> ezplot('(x^2-3*x+5)*exp(-5*x)*sin(x)',[0,1]), hold on >> sp1=csapi(x,y); fnplt(sp1,'--'); sp2=spapi(5,x,y); fnplt(sp2,':')