数值分析A实验报告

数值分析实验报告工程物理系二〇一六年一月十日

实验3.1（主元的选取与算法的稳定性）问题提出：Gauss消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss消去法作为数值算法的稳定性呢？Gauss消去法从理论算法到数值算法，其关键是主元的选择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。实验内容：考虑线性方程组编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss消去过程。实验要求：（1）取矩阵，则方程有解取n=10计算矩阵的条件数。分别用顺序GAUSS消元法、列主元GAUSS消元法、完全主元GAUSS消元法，结果如何？

（2）现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。（4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。重复上述实验，观察记录并分析实验结果。3.1.1程序清单formatlong;n=input('矩阵的阶数：n=');sp_M=input('矩阵的种类（1:Hilbert;2:随机矩阵；3:本题给出的矩阵；4：幻方矩阵）：sp_M=');switchsp_Mcase(1);A=hilb(n);case(2);A=round(8*rand(n));case(3);A=6*diag(ones(1,n),0)+8*diag(ones(1,n-1),-1)+diag(ones(1,n-1),1);case(4);A=magic(n);end;b=A*ones(n,1);p=input('计算条件数的p-范数，p=');cond_A=cond(A,p)Any1=zeros(1,n);Any20=zeros(n,1);Any21=zeros(n,1);Any12=eye(n);[m,n]=size(A);Ab=[Ab];Pro=input('计算方法（1：顺序高斯消元法；2,：列主元高斯消元法；3：完全主元高斯消元法；4：手动选主元法，Pro=');Abfori=1:n-1switchProcase(1);case(2);[aii,ip]=max(abs(Ab(i:n,i)));ip=ip+i-1;Any1=Ab(ip,:);Ab(ip,:)=Ab(i,:);Ab(i,:)=Any1;case(3);[Y,I]=max(max(abs(Ab(i:n,i:n))));%显示最大值列号I=I+i-1;[x1,r]=max(max(abs(Ab(i:n,i:n)')));%显示最大值行号r=r+i-1;Any2=Ab(:,I);Ab(:,I)=Ab(:,i);Ab(:,i)=Any2;%Ab阵I列与i列互换Any1=Ab(r,:);Ab(r,:)=Ab(i,:);Ab(i,:)=Any1;%Ab阵r行与i行互换Any21=Any12(:,I);Any12(:,I)=Any12(:,i);Any12(:,i)=Any21;%列交换跟踪case(4);ip=input(['第',num2str(i),'步消元，请输入第',num2str(i),'列所选元素所处行数：']);Any1=Ab(ip,:);Ab(ip,:)=Ab(i,:);Ab(i,:)=Any1;end;aii=Ab(i,i);fork=i+1:nif(aii~=0)Ab(k,i:n+1)=Ab(k,i:n+1)-(Ab(k,i)/aii)*Ab(i,i:n+1);elsebreak;end;end;Abend;x=zeros(n,1);x(n)=Ab(n,n+1)/Ab(n,n);fori=n-1:-1:1if(Pro==3)x(i)=(Ab(i,n+1)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);x=Any12^-1*x;elsex(i)=(Ab(i,n+1)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);end;endx3.1.2实验结果及分析（1）Cond(A,1)=2.557500000000000×103Cond(A,2)=1.727556024913821×103Cond(A,inf)=2.557500000000000×103未知数顺序高斯消元法列主元高斯消元法完全主元高斯消元法x11.00000000000000011x21.00000000000000011x31.00000000000000011x41.00000000000000111x50.99999999999999811x61.00000000000000411x70.99999999999999311x81.00000000000001211x90.99999999999997911x101.00000000000002811(2)手动选取主元，n=10未知数主元模最小或模尽可能小主元模最最大x11.0000000000000001x21.0000000000000001x31.0000000000000001x41.0000000000000011x50.9999999999999981x61.0000000000000041x70.9999999999999931x81.0000000000000121x90.9999999999999791x101.0000000000000281结果分析：由计算结果可知，主元取模最小结果没有主元取模最大结果好。由于最大主元与最小主元都没有过于小，都是相对好的一个数，所以结果差距很小。造成这样结果的原因是使用较小的数作为除数计算结果的误差会被放大。（3）n=20Cond(A,1)=2.621437500000000×106Cond(A,2)=1.789670565812014×106Cond(A,inf)=2.621437500000001×106未知数顺序高斯消元法列主元高斯消元法完全主元高斯消元法手动主元模最小或模尽可能小手动主元模最最大x11.000000000000000111.0000000000000001x21.000000000000000111.0000000000000001x31.000000000000000111.0000000000000001x41.000000000000001111.0000000000000011x50.999999999999998110.9999999999999981x61.000000000000004111.0000000000000041x70.999999999999993110.9999999999999931x81.000000000000014111.0000000000000141x90.999999999999972110.9999999999999721x101.000000000000057111.0000000000000571x110.999999999999886110.9999999999998861x121.000000000000227111.0000000000002271x130.999999999999547110.9999999999995471x141.000000000000902111.0000000000009021x150.999999999998209110.9999999999982091x161.000000000003524111.0000000000035241x170.999999999993179110.9999999999931791x181.000000000012732111.0000000000127321x190.999999999978173110.9999999999781731x201.000000000029102111.0000000000291021计算结果与二相比表明：1.选取模最大的元素作为主元会产生为精确的结果。2.在选取模最小元素作为主元时条件数越大误差越大。这是可以由摄动理论中事后误差估计是看出来的。（4）1、10阶Hilbert矩阵中各种方法解的情况Cond(A,1)=3.535336877175068×1013Cond(A,2)=1.602517099202631×1013Cond(A,inf)=3.535336877175066×1013未知数手动主元模最小或模尽可能小手动主元模最最大x10.9999999984611550.999999998758705x21.0000001311885301.000000106500618x30.9999972342961560.999997743217252x41.0000249388387821.000020435292621x50.9998818352853170.999902835789981x61.0003230379745661.000266409823714x70.9994724788979130.999563856087119x81.0005077323000321.000420698156409x90.9997343797724790.999779492550673x101.0000582339522601.0000484246488472、10阶随机矩阵中各种方法解的情况Cond(A,1)=70.699729258533296Cond(A,2)=1.581788216637637×102Cond(A,inf)=1.358379386423256×102未知数手动主元模最小或模尽可能小手动主元模最最大x10.9999999999999860.999999999999998x21.0000000000000381.000000000000005x30.9999999999999620.999999999999996x41.0000000000000281.000000000000007x50.9999999999999720.999999999999999x61.0000000000000071.000000000000000x71.0000000000000140.999999999999993x81.0000000000000060.999999999999998x90.9999999999999681.000000000000002x101.0000000000000301.0000000000000033、10阶幻方矩阵中各种方法解的情况Cond(A,1)=8.043128074178017×1017Cond(A,2)=7.670050592487186×1017Cond(A,inf)=8.129276476140978×1017未知数手动主元模最小或模尽可能小手动主元模最大x11.0000000122802001x20.9999999550964281x31.0000000962946981x40.9999999927228721x51.0000001239341051x60.9999999229185381x71.0000000309500471x80.9999999927228721x91.0000000007631181x101.0000000122802001一般来说，模最大元素作为主元比模最小的元素作为主元时的计算结果更精确。但一些方阵，如幻方矩阵，则是选择模最小的元素作为主元时计算结果最精确（选模最小的元素只是一个表象，这种选主元方法优于其他选主元方法的本质是这种选择方法能使消去过程不产生浮点数，而全是整数运算，只有在回代过程中才有可能会产生浮点数）。一般来说，需按模最大元素作为主元精度比较高。

实验3.3（病态的线性方程组的求解）问题提出：理论的分析表明，求解病态的线性方程组是困难的。实际情况是否如此，会出现怎样的现象呢？实验内容：考虑方程组Hx=b的求解，其中系数矩阵H为Hilbert矩阵，这是一个著名的病态问题。通过首先给定解（例如取为各个分量均为1）再计算出右端b的办法给出确定的问题。实验要求：（1）选择问题的维数为6，分别用Gauss消去法、J迭代法、GS迭代法和SOR迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比较，结论如何？（2）逐步增大问题的维数，仍然用上述的方法来解它们，计算的结果如何？计算的结果说明了什么？（3）讨论病态问题求解的算法3.3.1程序清单formatlong;n=input('矩阵的阶数：n=');Hibert_n=ones(n,n);forr=1:1:nforc=1:1:nHibert_n(r,c)=1/(r+c-1);end;end;b=Hibert_n*ones(n,1);Hibert_nb=[Hibert_nb];D=zeros(n,n);L=zeros(n,n);U=zeros(n,n);fori=2:1:nL(i,1:i-1)=-1.*Hibert_n(i,1:i-1);end;fori=1:1:n-1U(i,i+1:n)=-1.*Hibert_n(i,i+1:n);end;fori=1:1:nD(i,i)=Hibert_n(i,i);end;x=zeros(n,1);Medthod_number=input('计算方法（1：高斯消去法；2:Jacobi迭代；3：GS迭代；4：SOR迭代，Medthod_number=');Hibert_nbswitchMedthod_numbercase(1);fori=1:n-1Hibert_nbii=Hibert_nb(i,i);fork=i+1:nif(Hibert_nbii~=0)Hibert_nb(k,i:n+1)=Hibert_nb(k,i:n+1)-(Hibert_nb(k,i)/Hibert_nbii)*Hibert_nb(i,i:n+1);elsebreak;end;end;end;x=zeros(n,1);x(n)=Hibert_nb(n,n+1)/Hibert_nb(n,n);fori=n-1:-1:1x(i)=(Hibert_nb(i,n+1)-Hibert_nb(i,i+1:n)*x(i+1:n))/Hibert_nb(i,i);end;xcase(2);num_of_iter=0;norm_errorv=2;whilenorm_errorv>=10^-6xtemp=x;num_of_iter_temp=num_of_iter;x=D^-1*(L+U)*xtemp+D^-1*b;error_vector=x-xtemp;norm_errorv=norm(error_vector);num_of_iter=num_of_iter_temp+1;end;num_of_iterxcase(3);num_of_iter=0;norm_errorv=2;whilenorm_errorv>=10^-6xtemp=x;num_of_iter_temp=num_of_iter;x=(D-L)^-1*U*xtemp+(D-L)^-1*b;error_vector=x-xtemp;norm_errorv=norm(error_vector);num_of_iter=num_of_iter_temp+1;end;num_of_iterxcase(4);I=eye(n);B=I-D^-1*Hibert_n;spe_rB=max(abs(eig(B)))Wopt=1.5;num_of_iter=0;norm_errorv=2;Lw=(D-Wopt*L)^-1*((1-Wopt)*D+Wopt*U)whilenorm_errorv>=10^-6xtemp=x;num_of_iter_temp=num_of_iter;x=(D-Wopt*L)^-1*((1-Wopt)*D+Wopt*U)*xtemp+Wopt*(D-Wopt*L)^-1*b;error_vector=x-xtemp;norm_errorv=norm(error_vector);num_of_iter=num_of_iter_temp+1;end;Woptnum_of_iterxend;3.3.2运行结果及简要分析（1）n=6未知数GaussJGSSORx10.999999999999228Inf0.9999306341337850.999998116326860x21.000000000021937Inf1.0009194791104210.999536557238050x30.999999999851792NaN0.9981016368189121.005197221859359x41.000000000385370NaN0.9973780771021130.982933473818400x50.999999999574584NaN1.0089019779771601.021622489348798x61.000000000167680NaN0.9947049210396540.990659946952854迭代次数——4871740616769由结果知GAUSS法的结果最为精确。以error<=10^-6为收敛标准。GS法和SOR法收敛，但收敛的速度比较慢，SOR法略快于GS法。J法是发散的。（2）6阶未知数GaussJGSSORx10.999999999999228Inf0.9999306341337850.999998116326860x21.000000000021937Inf1.0009194791104210.999536557238050x30.999999999851792NaN0.9981016368189121.005197221859359x41.000000000385370NaN0.9973780771021130.982933473818400x50.999999999574584NaN1.0089019779771601.021622489348798x61.000000000167680NaN0.9947049210396540.990659946952854迭代次数——48717406167697阶未知数GaussJGSSORx10.999999999994453Inf0.9999993847606911.000077419518485x21.000000000221598NaN0.9994111921692410.997675627422646x30.999999997864368NaN1.0054159211198081.015125014152487x41.000000008303334NaN0.9869439523676270.965402133950665x50.999999984778069NaN1.0058756746177651.026302560980251x61.000000013152188NaN1.0106338440546081.004148021358175x70.999999995682024NaN0.9916438432417530.991184902946059迭代次数——4331174992368阶未知数GaussJGSSORx10.999999999966269-Inf1.0001010618346241.000110641227852x21.000000001809060NaN0.9974178826668250.997590430173341x30.999999976372675NaN1.0136096522484571.011909590979345x41.000000127868104NaN0.9794442972837040.981162475678398x50.999999655764117NaN0.9981715055106731.004262765594196x61.000000487042162NaN1.0140860365352361.005452790879320x70.999999653427125NaN1.0100862161756011.007936908023023x81.000000097774747NaN0.9869559021244590.991507742993205迭代次数——3968342379219阶未知数GaussJGSSORx10.999999999760212-Inf1.0001599663510921.000073440013547x21.000000016452158-Inf0.9969056403632590.998669263275490x30.999999722685967-Inf1.0126928847973041.004744337452915x41.000001973450782-Inf0.9875795363182760.998078248613115x50.999992779505688-Inf0.9915767143208420.991737673628083x61.000014713209384-Inf1.0049764504074231.002354905151437x70.999983130886114NaN1.0115911901907401.006475193208506x81.000010174801268NaN1.0059753549458001.004999791122969x90.999997489086183NaN0.9884243695012710.992808128443414迭代次数——368215713867510阶未知数GaussJGSSORx10.999999998754834Inf1.0001099892711001.000037307041751x21.000000106784973Inf0.9981823981941720.999664538362142x30.999997737861476Inf1.0056174964955880.998580945511472x41.000020479418515NaN0.9993251807758511.010945178465891x50.999902641847339NaN0.9915579706762920.985033582889975x61.000266907013337NaN0.9967960737850350.999312905811790x70.999563088470270NaN1.0052445079606501.001998084480765x81.000421401157600NaN1.0087341883658521.008630071839024x90.999779140796212NaN1.0038903039624241.003331412484509x101.000048498721852NaN0.9904417404459950.992405941672389迭代次数——347269513217411阶未知数GaussJGSSORx10.999999994751196-Inf1.0000678423903170.999999366078665x21.000000546746356-Inf0.9992080278903001.000670734390764x30.999985868343700-Inf1.0003082345229980.992687197515483x41.000157549468624NaN1.0069212532444701.022009340799829x50.999063537004343NaN0.9934019559184000.981143191862049x61.003286333127798NaN0.9923264105133340.996870387749579x70.992855789229368NaN0.9995569274820530.996749665623868x81.009726486881558NaN1.0065982879089851.008306329410864x90.991930155925812NaN1.0083329414393731.007566335246938x101.003729850349020NaN1.0028836723136931.003440352977447x110.999263885025643NaN0.9902968578591940.990482851656876迭代次数——330250402616312阶未知数GaussJGSSORx10.999999978015381-Inf1.0000259509399550.999954703394399x21.000002718658690-Inf1.0001743738870221.001785735264775x30.999916137983244-Inf0.9956523480373680.986627540764427x41.001124917423389-Inf1.0125908158347851.032053329265755x50.991858974446319NaN0.9961501109955520.979136034412255x61.035384713191155NaN0.9898767778078420.995292945294574x70.902314812813458NaN0.9947282800474210.991506439526719x81.175419475669440NaN1.0028964799507191.005967134258715x90.795767756284983NaN1.0084196081275921.008142963002572x101.148662570420705NaN1.0083434372718081.009201716248544x110.938525472294121NaN1.0018282708181941.001354831964326x121.011022484909254NaN0.9892065937798180.988879057012070迭代次数——316218012099613阶未知数GaussJGSSORx11.000000100845780-Inf0.9999797589992360.999898394014506x20.999984107801401-Inf1.0011764560567461.003102724057934x31.000614076745963-Inf0.9911882524521380.979985275134039x40.989763393555848-Inf1.0171284922885411.041822199898047x51.091975046430242-Inf0.9994937504545630.978339148728908x60.501003423820652-Inf0.9888375948008110.994598624657415x72.740851874285036-Inf0.9907570860497240.986468393979083x8-3.036058361310595-Inf0.9986697541967931.002570102424027x97.283884268701900-Inf1.0063066472730771.006486292294313x10-5.493529661418352-Inf1.0099777355518151.011620410991678x115.270874315267408-Inf1.0081295929919911.008018704221022x12-0.618191404610789NaN1.0005391576898431.000316599383520x131.268828909804817NaN0.9876890626527210.986641262097305迭代次数——304187201714914阶未知数GaussJGSSORx11.000000076434261-Inf0.9999254280475640.999821850368462x20.999989977372466-Inf1.0022742081668451.004719127001512x31.000309640223686-Inf0.9867031890927140.972637303241281x40.996220152145736-Inf1.0208247183007161.051064686323259x51.019948388320157-Inf1.0032765441908230.978784942485267x60.981737772434851-Inf0.9889493525394300.994964408389401x70.675197450814358-Inf0.9876706993474990.982038216196169x82.922896201795558-Inf0.9944160645575090.998668726968105x9-4.498100214715922NaN1.0029870007784921.003296483598423x1010.515964827807839NaN1.0092161216653711.011498111734201x11-9.428207690148314NaN1.0109778536801801.011357832701018x128.094762876614597NaN1.0075545030208791.007823563006471x13-1.740964386133816NaN0.9990553739301920.998296502569765x141.460245042732823NaN0.9860139144292810.984858091437531迭代次数——294161951466915阶未知数GaussJGSSORx10.999821850368462-Inf0.9998575571816410.999887622919398x21.004719127001512NaN1.0035255617953491.002415708045205x30.972637303241281NaN0.9821110358594410.988810342812522x41.051064686323259NaN1.0236479722255971.014030721059258x50.978784942485267NaN1.0074096037750861.002548533950694x60.994964408389401NaN0.9901678030952660.998016571956734x70.982038216196169NaN0.9855812613031880.991839105511616x80.998668726968105NaN0.9904621531699090.994102989105642x91.003296483598423NaN0.9990129280149990.997668442512087x101.011498111734201NaN1.0068192011970671.002777831762777x111.011357832701018NaN1.0113120959450541.006058415843252x121.007823563006471NaN1.0113368439412851.007066865308874x130.998296502569765NaN1.0066439110125291.004799747956600x140.984858091437531NaN0.9974992219972080.999304390919186x150.999821850368462NaN0.9844255314424750.990589637285929迭代次数——286144985007516阶未知数GaussJGSSORx11.000000038309608-Inf0.9998154791644970.999912245138576x20.999992356721026-Inf1.0037679427425541.001674705174921x31.000358126460529NaN0.9838633778807720.993777266155035x40.992915973395254NaN1.0162269207347921.002826246938333x51.074614485124896NaN1.0101717951463951.009452652207939x60.529206540133373NaN0.9964068272273280.999442700448315x72.896111260937623NaN0.9891272502238620.995130123070590x8-4.015764564087196NaN0.9895886176802570.992776438860727x99.670665346779543NaN0.9947899305784160.995497169367632x10-8.180072412591157NaN1.0013820248397070.999350896127904x115.338569007974306NaN1.0068625014596781.003550249866065x123.212675695583656NaN1.0097038543794491.006125996340453x13-3.852719521183188NaN1.0091571637045981.006509290332898x144.216873886870227NaN1.0050100426552411.004108083739434x15-0.004602922381234NaN0.9973850721265090.998887412656845x161.121176658583945NaN0.9865937795648230.990897211484367迭代次数——2782633351324矩阵阶数678910111213141516GS1740611749379212157126951250402180118720161951449826333SOR16769923683423867532174261632099617149146695007551324结果分析：1、从6阶一直算到16阶，6到12阶GAUSS法精度最高，从12阶开始误差急剧增大，从13阶以后GAUSS法结果误差很大。2、J法一直不收敛，原因是迭代矩阵谱半径大于1。3、GS法和SOR（w=1.5）法一直收敛，但收敛速度较慢。4、SOR法在大多数情况下收敛速度比GS法慢。（3）病态问题求解对原方程进行适当变换，取非奇异矩阵P,Q使得（PAQ）y=Pb,原方程的解x=Qy。用此种方法的目的是使矩阵PAQ的条件数cond(PAQ)<=cond(A),即使得条件数有所改善，对于1~100阶Hilbert矩阵一般选择对角阵D1=D2=D,设DM为Hilbert矩阵的对角阵。设使得cond(DHD)最小的一种。这样条件数的大小就会被改善。

4.1非线性方程组的解法问题提出：非线性方程组的求解方法很多，基本的思想是线性化。不同的方法效果如何，要靠计算的实践来分析、比较。

实验内容：考虑算法

(1)牛顿法

(2)拟牛顿法

分别编写它们的matlab程序。

实验要求：

(1)用上述方法，分别计算两个例子。在达到精度相同的前提下，比较迭代次数、浮点运算次数和CPU时间等。(2)取其他初值结果又如何？反复选取不同的初值比较其结果。（3）总结归纳你的实验结果，试说明各种方法适用的问题4.1.1程序清单*f1.m文件functiony=f1(x)y(1)=12*x(1)-x(2)^2-4*x(3)-7;y(2)=x(1)^2+10*x(2)-x(3)-11;y(3)=x(2)^3+10*x(3)-8;*ff1.m文件functiony=ff1(x)y(1,:)=[12,-2*x(2),-4];y(2,:)=[2*x(1),10,-1];y(3,:)=[0,3*x(2)^2,10];*f2.m文件functiony=f2(x)y(1)=3*x(1)-cos(x(2)*x(3))-0.5;y(2)=x(1)^2-81*(x(2)+0.1)^2+sin(x(3))+1.06;y(3)=exp(-1*x(2)*x(1))+20*x(3)+1/3*(10*pi-3);*ff2.m文件functiony=ff2%f2的雅克比矩阵y(1,:)=[3,sin(x(2)*x(3)*x(3)),sin(x(2)*x(3))*x(2)];y(2,:)=[2*x(1),162*(x(2)+0.1),cos(x(3))];y(3,:)=[exp(-1*x(2)*x(1))*(-1*x(2)),exp(-1*x(2)*x(1))*(-1*x(1)),20];以函数（1）为例牛顿法x10=[1,1,1]';x20=[0,0,0]';n=15;i=1;increx=[1,1,1]';x(:,1)=x10;incx=increx;tic; while(norm(incx)>10^(-6))&(i<n)incx=-ff1(x(:,1))\f1(x(:,1))';x(:,1)=x(:,1)+incx;i=i+1;end;iterations=i-1t=tocx(:,1)拟牛顿法x10=[1,1,1]';x20=[0,0,0]';n=15;i=1;increx=[1,1,1]';x(:,1)=x10;incx=increx;A=ff1(x(:,1));F=f1(x(:,1));tic;while(norm(incx)>10^(-6))&(i<n)xt=x(:,1);yt=f1(x(:,1));x(:,1)=x(:,1)-A\f1(x(:,1))';incx=x(:,1)-xt;incy=f1(x(:,1))'-yt';A=A+(incy-A*incx)*incx'/(incx'*incx);i=i+1;end;itera