線性最小二乘法擬合

《MATLAB程序設(shè)計(jì)實(shí)踐》課程考核一.實(shí)現(xiàn)以下科學(xué)計(jì)算算法，并舉一例應(yīng)用之。(參考書籍《精通MATLAB科學(xué)計(jì)算》，王正林等著，電子工業(yè)出版社，2009年)"線性最小二乘法擬合"解答：一，最小二乘法擬合；方法介紹在科學(xué)實(shí)驗(yàn)的統(tǒng)計(jì)方法研究中，往往要從一組實(shí)驗(yàn)數(shù)據(jù)(xi，yi)中尋找出自變量乂和因變量y之間的關(guān)系y=f(x)。由于觀測數(shù)據(jù)往往不準(zhǔn)確，因此并不要求y=f(x)經(jīng)過所有的點(diǎn)(xi,yi),而只要求在給定的點(diǎn)xi上誤差8=f(xi)-yi按照某種標(biāo)準(zhǔn)達(dá)到最小，通常采用歐氏范數(shù)。八2作為誤差量度的標(biāo)準(zhǔn)。這就是所謂的最小二乘法。MATLAB實(shí)現(xiàn)在MATLAB中實(shí)現(xiàn)最小二乘法擬合通?？梢圆捎萌缦峦緩剑豪胮olyfit函數(shù)進(jìn)行多項(xiàng)式擬合。polyfit函數(shù)源程序: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);流程圖：［例題］設(shè)y二span{1,x,x八2},用最小二乘法擬合如表所示數(shù)據(jù)。X0.51.01.52.02.53.0y1.752.453.814.808.008.60用polyfit功能函數(shù)進(jìn)行擬合。開始流程■/輸入數(shù)/據(jù)調(diào)用polyfit函數(shù)擬合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('實(shí)驗(yàn)值'，'擬合值')end在matlab命令窗口中輸入：task1()運(yùn)行結(jié)果: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編程解決以下科學(xué)計(jì)算問題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解答：此電路中設(shè)節(jié)點(diǎn)電壓為變量，共有ua,ub,uc,ud四個(gè)。將各支路電流均用這些電壓來表示.,U—UI=sb'i￡+R2,U-U《=十'3對b點(diǎn)和c點(diǎn)列出節(jié)點(diǎn)電流方程：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時(shí)，這兩個(gè)方程中只有兩個(gè)未知數(shù)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)，可以推出類似的矩陣表達(dá)式，只是輸入輸出量不同流程■:源程序代碼：代入矩陣方程代入矩陣方程源程序代碼：代入矩陣方程代入矩陣方程r1=2;r2=4;r3=4;r4=4;r5=2;r6=12;r7=12;%為元件賦值%將各系數(shù)矩陣賦值?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]%列出系數(shù)矩陣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]%列出系數(shù)矩陣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)運(yùn)行結(jié)果：輸入U(xiǎn)s=10時(shí)，Ubc=2.2222,Ude=0.7407,I7=0.3704輸入U(xiǎn)de=4時(shí)，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.求汁算演教差漸表?！暌浴蹬c出牛顏孝項(xiàng)式*：（3）在給定值H由上求中頓?頃式M（H）1*Z,3,1的值*:富⑷比較（3》中的結(jié)果與實(shí)際兩數(shù)值.差商表：源代碼：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運(yùn)行結(jié)果：下三角狀的牛頓差商表如下：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

>>寫出牛頓多項(xiàng)式Nk(x),k=1,2,3,4.在給定值x處求牛頓多項(xiàng)式Nk(x),k=1,2,3,4的值。比較(3)中的結(jié)果與實(shí)際函數(shù)值。(2),(3),(4)編為一個(gè)程序,(2)(3)以結(jié)果表示出來,(4)以圖形表示出來流程圖：k*nk*n源程序代碼：%輸入?yún)?shù)中x,y為元素個(gè)數(shù)相等的向量，七為待估計(jì)的點(diǎn)，可以為數(shù)字或向量。%輸出參數(shù)中P2為所求得的牛頓插值多項(xiàng)式，z為利用多項(xiàng)式所得的t的函數(shù)值。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),'次多項(xiàng)式'))vpa(pol