实验指导书_函数逼近

湖北民族学院理学院《数值分析》实验指导书陈以平编写1实验指导书_函数逼近与曲线拟合曲线拟合：由一组实验数据},,2,1),,{(niyxii，选择一个较简单的函数)(xf（如多项式），在一定准则下，最接近这组数据.函数逼近：已知一个较为复杂的连续函数],[),(baxxy，要求选择一个较简单的函数)(xf，在一定准则下最接近)(xy.掌握内容：切比雪夫零点插值、最佳平方逼近、曲线拟合的最小二乘法的Matlab实现一、切比雪夫零点插值1．算法原理设1[,],()nnfCabLx是以,0,1,,22kkbaabxtkn为插值节点的Lagrange插值多项式，其中21cos,0,1,,2(1)kktknn是Chebyshev多项式1()nTx的零点，此时插值误差最小，误差（余项）估计公式为：1(1)21()()()()max()2(1)!nnnnnaxbbaRxfxLxfxn2．实例例1（P64例4）求()xfxe在[0,1]上的四次Lagrange插值多项式4()Lx,插值节点用5()Tx的零点，并估计误差401max().xxeLx首先编写Lagrange插值主程序,并保存为lagrange.mfunction[L,C,l,L1]=lagrange(X,Y)%输入的量:n+1个节点(xi,yi)的横坐标向量X，纵坐标向量Y；%输出的量:n次拉格朗日插值多项式L及其系数向量C，基函数l及其系数矩阵L1m=length(X);L=ones(m,m);fork=1:mV=1;fori=1:mifk~=iV=conv(V,poly(X(i)))/(X(k)-X(i));endendL1(k,:)=V;l(k,:)=poly2sym(V);endC=Y*L1;L=Y*l;l=vpa(l,4);L=vpa(L,4);然后编写如下程序,并保存为lc_P64eg4.m%P64eg4求Lagrange-Chebyshev多项式的程序lc_P64eg4.m%X为插值节点湖北民族学院理学院《数值分析》实验指导书陈以平编写2%L为Lagrange-Chebyshev多项式的表达式%C为按降幂排列的Lagrange-Chebyshev多项式的系数%R为误差限clear;formatshortgs=1;fork=0:4X(k+1)=(1+cos((2*k+1)*pi/10))/2;Y(k+1)=exp(X(k+1));s=s*(k+1);endX,[L,C,l,L1]=lagrange(X,Y);L,C,R=2.71828/(s*2^9)x=linspace(X(1),X(5),50);y=polyval(C,x);y1=exp(x);plot(X,Y,'r*',x,y,'r-')figure,plot(X,Y,'r*',x,y1,'b-')figure,plot(X,Y,'r*',x,y,'r-',x,y1,'b-')在Matlab窗口中运行程序lc_P64eg4.mlc_P64eg4X=0.975530.793890.50.206110.024472L=.9988*x+.1403*x^3+.6942e-1*x^4+.5098*x^2+1.000C=0.0694160.140280.509780.998761R=4.4243e-005图1.()xfxe在[0,1]上的四次Lagrange-Chebyshev插值多项式湖北民族学院理学院《数值分析》实验指导书陈以平编写3例2(P65例5)设211()fxx,在[-5,5]上利用11()Tx的零点作插值点,构造10次拉格朗日插值多项式10()Lx,并与等距节点的10()Lx近似()fx作比较.首先编写如下程序,并保存为lc_P65eg5.m%P64eg4求Lagrange-Chebyshev多项式的程序lc_P65eg5.m%X为Chebyshev零点%L为Lagrange-Chebyshev多项式的表达式%C为按降幂排列的Lagrange-Chebyshev多项式的系数%L1为等距节点多项式的表达式%C1为按降幂排列的等距节点多项式的系数clear;formatshortfork=0:10X(k+1)=(10*cos((2*k+1)*pi/22))/2;Y(k+1)=1/(1+X(k+1)^2);endX,[L,C,l,L1]=lagrange(X,Y);L,C,x=linspace(X(1),X(11),100);y=polyval(C,x);y1=1./(1+x.^2);X1=linspace(-5,5,11);Y1=1./(1+X1.^2);[L1,C1,l1,L11]=lagrange(X1,Y1);L1,C1,X2=-5:0.1:5;Y2=polyval(C1,X2);plot(X,Y,'r*',x,y,'r-',x,y1,'b-',X1,Y1,'O',X2,Y2,'k-')gtext('切比雪夫插值点r*')gtext('切比雪夫插值曲线r-')gtext('被插函数曲线b-')gtext('等距节点插值点o')gtext('等距节点插值多项式k-')在Matlab窗口中运行程序lc_P65eg5.m,得如下图形结果(这里,数值结果略,请自己运行观察)图2.211()fxx在[-5,5]上的插值多项式10()Lx与10()Lx湖北民族学院理学院《数值分析》实验指导书陈以平编写4二、最佳平方多项式逼近1．求定义在区间[a,b]上的已知函数最佳平方逼近多项式的算法(1)设已知函数()fx的最佳平方逼近多项式为01()nnpxaaxax,由最佳平方逼近的定义有:010012(,,)(,,,,)niFaaaina其中20101(,,)()()bnnnaFaaafxaaxaxdx(2)求多项式()px系数的法方程:CaD,其中1211222111212()(),()(bbbbbnnaaaaabbbbbnnaaaaabbbbbnnnnnaaaaabbbbnnnnnaaaadxxdxxdxxdxfxdxxdxxdxxdxxdxxfxdxCDxdxxdxxdxxdxxfxdxxdxxdxxdxxdxxfx)badx2．求定义在区间[a,b]上的已知函数最佳平方逼近多项式的主程序function[coff,d]=ZJPF(func,n,a,b)%func:已知函数%n:最佳平方逼近多项式的最高次数%a:逼近区间的左端点%b:逼近区间的右端点%coff:按升幂排列的逼近多项式的系数%d:法方程的右端项C=zeros(n+1,n+1);var=findsym(sym(func));func=func/var;fori=1:n+1C(1,i)=(power(b,i)-power(a,i))/i;%算法中的C矩阵的第一行func=func*var;d(i,1)=int(sym(func),var,a,b);%算法中的D向量的第一行endfori=2:n+1C(i,1:n)=C(i-1,2:n+1);f1=power(b,n+i);f2=power(a,n+i);C(i,n+1)=(f1-f2)/(n+i);%形成C矩阵endcoff=C\d;%求解逼近多项式的系数湖北民族学院理学院《数值分析》实验指导书陈以平编写53.例题(教材P68例6)设211()fxx,求[0,1]上的一次最佳平方逼近多项式.首先编写被逼近函数的M文件,并保存为funcp68eg6.mfunctionf=funcp68eg6()symsxf=sqrt(1+x^2);然后编写如下程序,并保存为p68_eg6.m%p68_eg6.mclear;[coff,d]=ZJPF(funcp68eg6,1,0,1);coff=vpa(coff)x=findsym(sym(funcp68eg6));func=funcp68eg6*funcp68eg6;PFWC=vpa(int(sym(func),x,0,1)-coff'*d,4)运行结果p68_eg6coff=[.9343200492928959528618412015402][.4269470508068461683106156461091]PFWC=.12e-2三、曲线拟合的最小二乘法1.多项式拟合及其MATLAB程序面对一组数据,,,2,1),,(niyxii用线性最小二乘法作曲线拟合时，如果选取一组函数02(),(),,()mxxx为)(,,,,12nmxxxm，则拟合曲线为多项式11mmmaxaxay.一般m=2,3,不宜过高.对于指数曲线，拟合前需作变量代换，化为系数参数的线性函数.2.用MATLAB作线性最小二乘拟合的多项式拟合用MATLAB作线性最小二乘拟合的多项式拟合有现成程序.调用格式为:a=polyfit(x,y,m)其中输入参数x,y为要拟合的数据,是长度自定义的数组,m为拟合多项式的次数,输出参数a为拟合多项式:11mmmaxaxay的系数向量a),,,(11mmaaa.(注意:按降幂排列)湖北民族学院理学院《数值分析》实验指导书陈以平编写63.多项式在x处的值y的计算:y=polyval(a,x)例给出一组数据点),(iiyx列入下表中，试用线性最小二乘法求拟合曲线，并用估计其误差，作出拟合曲线.xi-2.9-1.9-1.1-0.800.11.52.73.6yi53.9433.6820.8816.928.798.984.179.1219.88解（1）首先根据给出的数据点),(iiyx，用下列MATLAB程序画出散点图.在MATLAB工作窗口输入程序x=[-2.9-1.9-1.1-0.800.11.52.73.6];y=[53.9433.6820.8816.928.798.984.179.1219.88];plot(x,y,'r*'),legend('数据点(xi,yi)')xlabel('x'),ylabel('y'),title('例7.4.1的数据点(xi,yi)的散点图')运行后屏幕显示数据的散点图，见图3.（2）因为数据的散点图3的变化趋势与二次多项式很接近，所以选取一组函数2,,1xx，令3221)(axaxaxf,其中ka是待定系数)3,2,1(k.（3）用作线性最小二乘拟合的多项式拟合的MATLAB程序求待定系数ka)3,2,1(k.输入程序a=polyfit(x,y,2)运行后输出按拟合多项式的系数a=2.8302-7.37219.1382故拟合多项式为2138.91372.72830.2)(2xxxf.（4）编写下面的MATLAB程序估计其误差，并做出拟合曲线和数据的图形.输入程序xi=[-2.9-1.9-1.1-0.800.11.52.73.6];y=[53.9433.6820.8816.928.798.984.179.1219.88];n=length(xi);f=2.8302.*xi.^2-7.3721.*xi+9.1382x=-2.9:0.001:3.6;F=2.8302.*x.^2-7.3721.*x+8.79;fy=abs(f-y);fy2=fy.^2;Er=sqrt((sum(fy2)/n)),plot(xi,y,'r*',x,F,'b-'),legend('数据点(xi,yi)','拟合曲线y=f(x)')xlabel('x'),ylabel('y'),运行后屏幕显示数据),(iiyx与拟合函数f的均方根误差E2及其数据点(xi,yi)和拟合曲线y=f(x)的图形，见图4.湖北民族学院理学院《数值分析》实验指导书陈以平编写7图3图4