基于MATLAB的数值分析(6)

第六章数值逼近6.1关于多项式的MATLAB命令一、多项式表达方式的约定多项式用行向量表示用比较习惯的方式显示多项式：pp=poly2str(p,'x')1121()...nnnnPxaxaxaxa121[...]nnpaaaa【例】多项式可表示为p=[2145]pp=poly2str(p,'x')32()245Pxxxx二、多项式运算函数r=roots(p)：求多项式的零点p=poly(r):以r为零点的多项式p=poly(A):A的特征多项式PA=polyval(p,S):按数组运算规则，计算多项式的值其中S，PA为矩阵PM=polyvalm(p,S):按矩阵运算规则，计算多项式的值,其中S，PM为矩阵p=conv(p1,p2)：多项式的乘积[q,r]=deconv(p1,p2)：多项式的除法，p1/p2p1(x)=p2(x)q(x)+r(x)【例】由给定根向量求多项式系数向量。R=[-0.5,-0.3+0.4*i,-0.3-0.4*i];P=poly(R)PPR=poly2str(P,'x')P=1.00001.10000.55000.1250PPR=x^3+1.1x^2+0.55x+0.125【例】求多项式的零点。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.0024i665432(1)615201561yxxxxxxx注：尽管利用MATLAB使得从系数转换到零点或从零点转换到系数都非常容易，但是使用时一定要注意计算的精度。如果存在重根，这种转换可能会降低精度。对于数值计算，计算重根是最困难的问题之一。【例】求3阶方阵A的特征多项式。A=[111213;141516;171819];PA=poly(A)PPA=poly2str(PA,'x')PA=1.0000-45.0000-18.0000-0.0000PPA=x^3-45x^2-18x-2.8387e-015【例】求的“商”及“余”多项式。p1=conv([1,0,2],conv([1,4],[1,1]));p2=[1011];[q,r]=deconv(p1,p2);cq='商多项式为';cr='余多项式为';disp([cq,poly2str(q,'x')])disp([cr,poly2str(r,'x')])商多项式为x+5余多项式为5x^2+4x+31)1)(4)(2(32sssss【例】两种多项式求值指令的差别。S=pascal(4)P=poly(S);PP=poly2str(P,'s')PA=polyval(P,S)PM=polyvalm(P,S)S=1111123413610141020PP=s^4-29s^3+72s^2-29s+1PM=1.0e-010*0.00160.00330.00900.02050.00450.01010.02860.06970.00950.02100.06530.15960.01630.03870.12260.3019PA=1.0e+004*0.00160.00160.00160.00160.00160.0015-0.0140-0.05630.0016-0.0140-0.2549-1.20890.0016-0.0563-1.2089-4.3779可以用命令polyval或polyvalm计算多项式的值p=[2145];xi=2.5;yi=polyval(p,xi)yi=52.5000yi=polyvalm(p,xi)yi=52.5000【例】求多项式的积分poly_itg.mfunctionpy=poly_itg(p)n=length(p)py=[p.*[n:-1:1].^(-1),0]【例】将Chebyshev多项式展开为幂级数形式Cheby_pw.mfunctionpn=Cheby_pw(n)pbb=[1];ifn==0,pn=pbb;break;endpb=[10];ifn==1,pn=pb;break;endfori=2:npn=2*[pb,0]-[0,0,pbb];pbb=pb;pb=pn;end6.2Lageange插值多项式求由n+1维向量x,y给定数据的n次插值多项式pp=polyfit(x,y,n)【例】观察Lagrange插值基函数的形状x=[12345678];y=[01000000];p=polyfit(x,y,7);x0=1:0.2:8;y0=polyval(p,x0);plot(x,y,'ro',x0,y0)Lagr1.mfunctiony=lagr1(x0,y0,x)n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;endy=lagr1(x0,y0,x):Lagrange插值。给出n个插值节点，计算m个插值点。chzh00.mn=11;m=51;x=-1:2/(m-1):1;y=1./(1+25*x.^2);z=0*x;x0=-1:2/(n-1):1;y0=1./(1+25*x0.^2);y1=lagr1(x0,y0,x);N1=6;x00=-1:2/(n1-1):1;y00=1./(1+25*x00.^2);y2=lagr1(x00,y00,x);z1=-0.4:0.2:1.6;x1=0*z1;plot(x,z,'k',x1,z1,'k',x,y,'r:',x,y1,'b',x,y2,'c')gtext('n=5'),gtext('n=10'),gtext('y=1./(1+25*x.^2)')21(),[1,1]125fxxx观察Lagrange多项式的振荡现象，即Runge现象。观察函数：6.3插值多项式的误差及Chebyshev点的插值4sin()sin()xyxexxLx用0上的5个等距点对函数进行插值估计。误差函数定义为例：functionwuch1m=101;n=5;x=0:pi/(m-1):pi;y=sin(x);z=0*x;x0=0:pi/(n-1):pi;y0=sin(x0);y1=lagr1(x0,y0,x);er1=(y-y1)*100;z1=0:0.2:1;x1=0*z1;plot(x,z,'k',x1,z1,'k',x,er1,'b',x,y,'r')gtext('e(x)*100'),gtext('sin(x)')可以观察到误差曲线是震荡的，在接近端点的区间上最大，这种误差特性在等距多项式插值中非常典型，实际上，误差分布形状上的变化还取决于被插函数的性质以及插值区间的大小|b-a|.(1)011()()()()()()()()()()(1)!nnnnexfxLxLxfxxxxxxxxLxn分析Lagrange插值的误差:其中(1)(1)()()()()max()nnabexLxfLxf2.[0,],5wuchm：等距节点3.[0,],52wuchm：等距节点4.[0,],7wuchm：等距节点5.[0,],7wuchm：个Chebyshev点降低误差的三条途径functionwuch0(n)m=101;x=-1:2/(m-1):1;y=1./(1+x.^2);z=0*x;x0=-1:2/(n-1):1;y0=1./(1+x0.^2);y1=lagr1(x0,y0,x);er1=(y-y1)*100;x00=sort(roots(Cheby_pw(n)));y00=1./(1+x00.^2);y2=lagr1(x00,y00,x);er2=(y-y2)*100;z1=-0.1:0.2:0.2;x1=0*z1;plot(x,z,'k',x1,z1,'k',x,er1,'b',x,er2,'r')gtext('Lagr.'),gtext('Cheby.')n=11总结：（1）建议尽可能在小区间上使用多项式插值。（2）只能在一定范围内依靠增加插值点个数提高插值精度，如果插值点个数过多往往会适得其反。6.4分段插值interp1:一维分段插值interp2:二维分片插值interp3:三维分块插值spline:三次样条插值y0=interp1(x,y,x0,'插值方法选项')plot(x,y,x0,y0):绘插值曲线‘插值方法选项’有以下四种（1）'nearest':最接近点插值（2）'linear':分段线性插值（3）'cubic':分段三次插值（4）'spline':分段三次样条插值演示程序intplt.m一、一维分段插值二、三次样条插值(1)y0=interp1(x,y,x0,'spline')y0=interp2(x,y,x0,'spline')（2）y0=spline(x,y,x0)(一维)（3）pp=spline(x,y);（结构数组）y0=fnval(pp,x0);poly(x0,y0),pp.coefs（系数）drp=fnder(pp,n),drv=fnval(drp,a)求导数itp=fnint(pp),itv=fnval(itp,[a,b])*[-1,1]求积分fnplt(pp),fnplt(drp,n),fnplt(itp)绘图（4）pp=csape(x,y,'边界条件选项');（5）pp=csapi(x,y);演示程序sptest.m【例】根据连续时间函数的采样数据，利用spline重构该连续函数，并检查重构误差。t=-5:0.5:5;w=exp(-abs(t));N0=length(t);tt=linspace(t(1),t(end),10*N0);ww=spline(t,w,tt);error=max(abs(ww-exp(-abs(tt))))plot(tt,ww,'b');holdonstem(t,w,'filled','r');holdofferror=0.0840-50500.10.20.30.40.50.60.70.80.91【例】对于函数，很容易求得。本例将借此演示样条函数求数值不定积分、导函数的能力。xysinxxxdxxS0cos1sin)(xycos（1）不定积分样条函数、导数样条函数的求取和精度分析。x=(0:0.1:1)*2*pi;y=sin(x);pp=spline(x,y);int_pp=fnint(pp);der_pp=fnder(pp);%xx=(0:0.01:1)*2*pi;err_yy=max(abs(ppval(pp,xx)-sin(xx)))err_int=max(abs(ppval(int_pp,xx)-(1-cos(xx))))err_der=max(abs(ppval(der_pp,xx)-cos(xx)))err_yy=0.0026err_int=0.0010err_der=0.0253（2）不定积分样条函数、导数样条函数的使用%计算y(x)在区间[1,2]上的定积分DefiniteIntegral.bySpline=ppval(int_pp,[1,2])*[-1;1];DefiniteIntegral.byTheory=(1-cos(2))-(1-cos(1));%计算dy(3)/dxDerivative.bySpline=fnval(der_pp,3);Derivative.byTheory=cos(3);Derivative.byDiference=(sin(3.01)-sin(3))/0.01;DefiniteIntegral,DerivativeDefiniteIntegral=bySpline:0.9563byTheory:0.9564Derivative=bySp