线性最小二乘法拟合

《MATLAB程序设计实践》课程考核一.实现以下科学计算算法，并举一例应用之。(参考书籍《精通MATLAB科学计算》，王正林等著，电子工业出版社，2009年)"线性最小二乘法拟合"解答：一，最小二乘法拟合；方法介绍在科学实验的统计方法研究中，往往要从一组实验数据(xi，yi)中寻找出自变量乂和因变量y之间的关系y=f(x)。由于观测数据往往不准确，因此并不要求y=f(x)经过所有的点(xi,yi),而只要求在给定的点xi上误差8=f(xi)-yi按照某种标准达到最小，通常采用欧氏范数。八2作为误差量度的标准。这就是所谓的最小二乘法。MATLAB实现在MATLAB中实现最小二乘法拟合通常可以采用如下途径：利用polyfit函数进行多项式拟合。polyfit函数源程序:function[p,S,mu]=polyfit(x,y,n)%POLYFITFitpolynomialtodata.%P=POLYFIT(X,Y,N)findsthecoefficientsofapolynomialP(X)of%degreeNthatfitsthedataYbestinaleast-squaressense.Pisa%rowvectoroflengthN+1containingthepolynomialcoefficientsin%descendingpowers,P(1)*X八N+P(2)*X八(N-1)+...+P(N)*X+P(N+1).%%[P,S]=POLYFIT(X,Y,N)returnsthepolynomialcoefficientsPanda%structureSforusewithPOLYVALtoobtainerrorestimatesfor%predictions.Scontainsfieldsforthetriangularfactor(R)fromaQR%decompositionoftheVandermondematrixofX,thedegreesoffreedom

(df),andthenormoftheresiduals(normr).IfthedataYarerandom,%anestimateofthecovariancematrixofPis(Rinv*Rinv')*normr八2/df,%whereRinvistheinverseofR.%%[P,S,MU]=POLYFIT(X,Y,N)findsthecoefficientsofapolynomialin%XHAT=(X-MU(1))/MU(2)whereMU(1)=MEAN(X)andMU(2)=STD(X).This%centeringandscalingtransformationimprovesthenumericalproperties%ofboththepolynomialandthefittingalgorithm.%%WarningmessagesresultifNis>=length(X),ifXhasrepeated,or%nearlyrepeated,points,orifXmightneedcenteringandscaling.%%ClasssupportforinputsX,Y:%float:double,single%%SeealsoPOLY,POLYVAL,ROOTS.%Copyright1984-2004TheMathWorks,Inc.%$Revision:5.17.4.5$$Date:2004/07/0517:01:37$%Theregressionproblemisformulatedinmatrixformatas:%%%%y=%%%%%y=%%%%y=V*por32[xxx1][p3p2p1p0]%wherethevectorpcontainsthecoefficientstobefound.Fora%7thorderpolynomial,matrixVwouldbe:%V=[x.八7x.八6x.八5x.八4x.八3x.八2xones(size(x))];if~isequal(size(x),size(y))error('MATLAB:polyfit:XYSizeMismatch',...'XandYvectorsmustbethesamesize.')endx=x(:);y=y(：);ifnargout>2mu=[mean(x);std(x)];x=(x-mu(1))/mu(2);end%ConstructVandermondematrix.V(:,n+1)=ones(length(x),1,class(x));forj=n:-1:1V(:,j)=x.*V(:,j+1);end%Solveleastsquaresproblem.[Q,R]=qr(V,0);ws=warning('off','all');p=R\(Q'*y);%Sameasp=V\y;warning(ws);ifsize(R,2)>size(R,1)warning('MATLAB:polyfit:PolyNotUnique',...'Polynomialisnotunique;degree>=numberofdatapoints.')elseifcondest(R)>1.0e10ifnargout>2warning('MATLAB:polyfit:RepeatedPoints',...'Polynomialisbadlyconditioned.Removerepeateddatapoints.')elsewarning('MATLAB:polyfit:RepeatedPointsOrRescale',...['Polynomialisbadlyconditioned.Removerepeateddatapoints\nortrycenteringandscalingasdescribedinHELPPOLYFIT.'])endendr=y-V*p;p=p.';%Polynomialcoefficientsarerowvectorsbyconvention.%Sisastructurecontainingthreeelements:thetriangularfactorfroma%QRdecompositionoftheVandermondematrix,thedegreesoffreedomand%thenormoftheresiduals.S.R=R;S.df=length(y)-(n+1);S.normr=norm(r);流程图：［例题］设y二span{1,x,x八2},用最小二乘法拟合如表所示数据。X0.51.01.52.02.53.0y1.752.453.814.808.008.60用polyfit功能函数进行拟合。开始流程■/输入数/据调用polyfit函数拟合I终止建立m文件源代码：'、functiontask1()x=0.5:0.5:3.0;y=[1.752.453.814.808.008.60];a=polyfit(x,y,2)x1=0.5:0.05:3.0y1=a(3)+a(2)*x1+a(1)*x1.“2plot(x,y,'*')holdonplot(x1,y1,'-r')legend('实验值'，'拟合值')end在matlab命令窗口中输入：task1()运行结果:task1()a=0.49001.25010.8560x1=Columns1through90.50000.55000.60000.65000.70000.75000.80000.85000.9000Columns10through180.95001.00001.05000.95001.00001.05001.10001.15001.20001.25001.30001.3500Columns19through271.40001.45001.50001.40001.45001.50001.55001.60001.65001.70001.75001.8000Columns28through361.85001.90001.95001.85001.90001.95002.00002.05002.10002.15002.20002.2500Columns37through452.30002.35002.40002.30002.35002.40002.45002.50002.55002.60002.65002.7000Columns46through512.75002.80002.75002.80002.85002.90002.95003.0000y1=Columns1through91.60361.69181.78251.60361.69181.78251.87561.97122.06922.16972.27262.3780Columns10through182.48592.59612.70892.82412.94173.06183.18433.30933.4367Columns19through273.56663.69893.83373.56663.69893.83373.97094.11064.25284.39734.54444.6939Columns28through364.84585.00025.15704.84585.00025.15705.31635.47805.64225.80885.97796.1494Columns37through456.32346.49996.67876.32346.49996.67876.86017.04397.23017.41887.60997.8035Columns46through517.99958.19808.39898.60238.80819.0164编程解决以下科学计算问题1.按图8-1-1所示的梯形直流电路，问2nr~~卜—1A404n4nr，+mlJ1211LJF1-rIS8-1-J例H-1-1的电路图（1）如Us=10V，求Ubc，I7，Ude（2）如Ude=4V，求Ubc，I7，Us解答：此电路中设节点电压为变量，共有ua,ub,uc,ud四个。将各支路电流均用这些电压来表示.,U—UI=sb'i£+R2,U-U《=十'3对b点和c点列出节点电流方程：I1=I3+I7，I3=I4+I6，将上述各值代入化简，得梏+二+上）顷31URb3、345/—UUI7=黄'I4=I5=E7451U=二RcR+R1、R+R12(1—++—"RR+RR)「U]aa「U]「U[s_「a]Ubc」=11aL2112a22」Ubc」=L0J=13aL23」对问题（1）给出Us时，这两个方程中只有两个未知数Ub和Uc，可以写成矩阵方程其中：a11=1/(r1+r2)+1/r3+1/r7;a12=-1/r3;a13=1/(r1+r2);a21=1/r3;a22=-1/r3-1/(r4+r5)-1/r6;a23=0并用矩阵左除解出Ub和Uc再由RUU=U-U,U=5U,I=b,bcbcdeR—+Rc7~R~即可得到问题(1)的解。对问(2)，可以推出类似的矩阵表达式，只是输入输出量不同流程■:源程序代码：代入矩阵方程代入矩阵方程源程序代码：代入矩阵方程代入矩阵方程r1=2;r2=4;r3=4;r4=4;r5=2;r6=12;r7=12;%为元件赋值%将各系数矩阵赋值„a11=1/(r1+r2)+1/r3+1/r7;a12=-1/r3;a13=1/(r1+r2);a21=1/r3;a22=-1/r3-1/(r4+r5)-1/r6;a23=0;Us=input('Us=');%输入解(1)的已知条件A1=[a11,a12;a21,a22]%列出系数矩阵A1u=A1\[a13*Us;0]%u=[ub;uc]Ubc=u(1)-u(2),Ude=u(2)*r5/(r4+r5),I7=u(1)/r7%解(1)Ude=input('Ude='),%输入解(2)的已知条件A2=[A1,[a13;a23];0,r5/(r4+r5),0]%列出系数矩阵A2u=A2\[0;0;Ude],%u=[ub;uc;us]Ubc=u(1)-u(2),I7=u(1)/r7,Us=(r1+r2)*(u(1)*a11+u(2)*a12)%解(2)运行结果：输入Us=10时，Ubc=2.2222,Ude=0.7407,I7=0.3704输入Ude=4时，Ubc=12.0000,I7=2.0000,Us=54.0000如下：task2Us=10A1=0.5000-0.25000.2500-0.5000u=4.44442.2222Ubc=2.2222Ude=0.7407I7=0.3704Ude=4Ude=4A2=0.5000-0.25000.16670.2500-0.5000000.33330u=24.000012.0000Ubc=12.0000I7=2.0000Us=54.0000>>2.求汁算演教差渐表。£以〉与出牛颜孝项式*：（3）在给定值H由上求中顿•顷式M（H）1*Z,3,1的值*:富⑷比较（3》中的结果与实际两数值.差商表：源代码：x=[4.05.06.07.08.0];y=[2.0002.36072.449492.645752.82843];n=length(x);newton=[x',y'];forj=2:nfori=n:-1:1ifi>=jy(i)=(y(i)-y(i-1))./(x(i)-x(i-j+1));elsey(i)=0;endendnewton=[newton,y'];enddisp('下三角状的牛顿差商表如下：')newton运行结果：下三角状的牛顿差商表如下：newton=4.00002.000000005.00002.36070.36070006.00002.44950.0888-0.1360007.00002.64580.19630.05370.063208.00002.82840.1827-0.0068-0.0202-0.0209

>>写出牛顿多项式Nk(x),k=1,2,3,4.在给定值x处求牛顿多项式Nk(x),k=1,2,3,4的值。比较(3)中的结果与实际函数值。(2),(3),(4)编为一个程序,(2)(3)以结果表示出来,(4)以图形表示出来流程图：k*nk*n源程序代码：%输入参数中x,y为元素个数相等的向量，七为待估计的点，可以为数字或向量。%输出参数中P2为所求得的牛顿插值多项式，z为利用多项式所得的t的函数值。x=[45678];y=[22.36072.449492.645752.82843];t=[4.5,7.5];n=length(x);chaS(1)=y(1);fori=2:nx1=x;y1=y;x1(i+1:n)=[];y1(i+1:n)=[];n1=length(x1);s1=0;forj=1:n1t1=1;fork=1:n1ifk==jcontinue;elset1=t1*(x1(j)-x1(k));endends1=s1+y1(j)/t1;endchaS(i)=s1;endb(1,:)=[zeros(1,n-1)chaS(1)];cl=cell(1,n-1);fori=2:nu1=1;forj=1:i-1u1=conv(u1,[1-x(j)]);cl(i-1}=u1;endcl{i-1}=chaS(i)*cl{i-1};b(i,:)=[zeros(1,n-i),cl(i-1}];endp2=b(1,:);forq=2:np2=p2+b(q,:);disp(strcat('牛顿',num2str(q-1),'次多项式'))vpa(pol