数值分析上机作业2.doc

1、数值实验数值实验综述：线性代数方程组的解法是一切科学计算的基础与核心问题。求解方法大致可分为直接法和迭代法两大类。直接法指在没有舍入误差的情况下经过有限次运算可求得方程组的精确解的方法，因此也称为精确法。当系数矩阵是方的、稠密的、无任何特殊结构的中小规模线性方程组时，Gauss 消去法是目前最基本和常用的方法。如若系数矩阵具有某种特殊形式，则为了尽可能地减少计算量与存储量，需采用其他专门的方法来求解。Gauss 消去等同于矩阵的三角分解，但它存在潜在的不稳定性，故需要选主元素。对正定对称矩阵，采用平方根方法无需选主元。方程组的性态与方程组的条件数有关，对于病态的方程组必须采用特殊的方法进行求解

2、。实验一一、 实验名称：矩阵的 LU 分解.二、 实验目的：用不选主元的 LU 分解和列主元 LU 分解求解线性方程组 Ax=b, 并比较这两种方法.三、 实验内容与要求（1）用所熟悉的计算机语言将不选主元和列主元 LU 分解编成通用的子程序，然后用编写的程序求解下面的 84 阶方程组将计算结果与方程组的精确解进行比较，并就此谈谈你对 Gauss 消去法的看法。（2）写出追赶法求解三对角方程组的过程，并编写程序求该实验中的方程组（1）不选主元高斯消去法求解方程组function L,U=gauss1(A,b)n=length(A);for k=1:n-1A(k+1:n,k)=A(k+1:n,k

3、)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);endL=tril(A(:,1:n),-1)+eye(n);U=triu(A);主程序为：Clear； clc；a=ones(84,1)*6;b=ones(83,1);c=ones(83,1)*8;A=diag(a)+diag(b,1)+diag(c,-1);d=ones(82,1)*15;b=7;d;14;L U=gauss1(A,b);n=length(A);y=ones(n,1);y=Lb;x=ones(n,1);x=Uy解为：x=11.000000000000001

4、.000000000000001.000000000000000.9999999999999981.000000000000000.9999999999999931.000000000000010.9999999999999721.000000000000060.9999999999998861.000000000000230.9999999999995451.000000000000910.9999999999981811.000000000003640.9999999999927251.000000000014550.9999999999708981.000000000058200.999

5、9999998835921.000000000232820.9999999995343671.000000000931270.9999999981374691.000000003725060.9999999925498741.000000014900250.9999999701994971.000000059601010.9999998807979861.000000238404030.9999995231919461.000000953616110.9999980927677831.000003814464440.9999923710711301.000015257857740.999969

6、4842845201.000061031430960.9998779371380811.000244125723840.9995117485523231.000976502895350.9980469942092911.003906011581410.9921879768371871.015624046325570.9687519073490881.062496185300920.8750076294018071.249984741181830.5000305176945351.99993896437811-0.9998779278249614.99975585192486-6.9995116

7、889494716.9990233182979-30.998046398191864.9960918427676-126.992179871071256.984344484284-510.9686279371361024.93701174855-2046.873046994204096.74218797682-8190.4687519073216383.8750076293-32764.500030517565531.0001220701-131055.000488281262097.001953124-524127.0078124981048001.03125000-2094975.1249

8、99994185857.49999999-8355328.9999999716645128.9999999-33028126.999999965007744.9999998-125821439.000000234866688.999999-402628606.999999536838144.999998可见，这是一个病态方程，从 56 个跟开始发散。可见|A|=7.48288838313422e+50 不等于零，所以是非奇异矩阵，这主要是一个小主元问题引起的，需要更高精度的计算机来计算。选主元的高斯消去法function L,U,P=gauss2(A,b)n,l=size(A);P=eye(n

9、);for k=1:n-1m,p=max(abs(A(:,k);A(k,p,:)=A(p,k,:);P(k,p,:)=P(p,k,:);if A(k,k)=0A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);endendL=tril(A(:,1:n),-1)+eye(n);U=triu(A);P;主程序为：clearclca=ones(84,1)*6;b=ones(83,1);c=ones(83,1)*8;A=diag(a)+diag(b,1)+diag(c,-1);d=ones(8

10、2,1)*15;b=7;d;14;L U P=gauss2(A,b);n=length(A);y=ones(n,1);y=L(P*b);x=ones(n,1);x=Uy解为：1.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000

11、000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001

12、.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.00000000000

13、0000.9999999999999991.000000000000000.9999999999999951.000000000000010.9999999999999791.000000000000040.9999999999999171.000000000000170.9999999999996671.000000000000670.9999999999986661.000000000002670.9999999999946641.000000000010670.9999999999786571.000000000042690.9999999999146291.00000000017074

14、0.9999999996585261.000000000682930.9999999986342271.000000002731210.9999999945389091.000000010916850.9999999781876491.000000043539330.9999999132628241.000000172108410.9999996612469361.000000655651090.9999987761179611.000002098083500.999997202555339这个病态方程的解没有发散，用选主元的方法可以有效的提高解的精度。（2）追赶法clearclca=ones

15、(84,1)*6;e=ones(84,1);c=ones(84,1)*8;d=ones(82,1)*15;b=7;d;14;n,l=size(b);tic;U(1)=a(1);for i=2:nL(i)=c(i)/U(i-1);U(i)=a(i)-L(i)*e(i-1);endfor i=2:ny(1)=b(1);y(i)=b(i)-L(i)*y(i-1);endfor i=n-1:-1:1x(n)=y(n)/U(n);x(i)=(y(i)-e(i)*x(i+1)/U(i);endt3=tocx=x解为：11.000000000000001.000000000000001.0000000000

16、00000.9999999999999981.000000000000000.9999999999999931.000000000000010.9999999999999721.000000000000060.9999999999998861.000000000000230.9999999999995451.000000000000910.9999999999981811.000000000003640.9999999999927251.000000000014550.9999999999708981.000000000058200.9999999998835921.0000000002328

17、20.9999999995343671.000000000931270.9999999981374691.000000003725060.9999999925498741.000000014900250.9999999701994971.000000059601010.9999998807979861.000000238404030.9999995231919461.000000953616110.9999980927677831.000003814464440.9999923710711301.000015257857740.9999694842845201.000061031430960.

18、9998779371380811.000244125723840.9995117485523231.000976502895350.9980469942092911.003906011581410.9921879768371871.015624046325570.9687519073490881.062496185300920.8750076294018071.249984741181830.5000305176945351.99993896437811-0.9998779278249614.99975585192486-6.9995116889494716.9990233182979-30.

19、998046398191864.9960918427676-126.992179871071256.984344484284-510.9686279371361024.93701174855-2046.873046994204096.74218797682-8190.4687519073216383.8750076293-32764.500030517565531.0001220701-131055.000488281262097.001953124-524127.0078124981048001.03125000-2094975.124999994185857.49999999-835532

20、8.9999999716645128.9999999-33028126.999999965007744.9999998-125821439.000000234866688.999999-402628606.999999536838144.999998实验二一、 实验名称：实对称正定矩阵的 的 Cholesky 分解.A二、 实验目的：用平方根法和改进的平方根方法求解线性方程组 Ax=b.三、 实验内容与要求用所熟悉的计算机语言将 Cholesky 分解和改进的 Cholesky 分解编成通用的子程序，然后用编写的程序求解对称正定方程组 Ax=b，其中（1） b 随机的选取，系数矩阵为 100

21、阶矩阵（2） 系数矩阵为 40 阶 Hilbert 矩阵，即系数矩阵 A 的第 i 行第 j 列元素为，向量 b 的第 i 个分量为（3） 用实验一的程序求解这两个方程组，并比较所有的计算结果，然后评价各个方法的优劣。平方根方法求解function L=cholesky1(A)n=length(A);tic;for k=1:nA(k,k)=sqrt(A(k,k);A(k+1:n,k)=A(k+1:n,k)/A(k,k);for j=k+1:nA(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);endendt1=tocL=tril(A);改进的平方根方法function L,D=

22、cholesky2(A)n=length(A);v=ones(n,1);for j=1:nfor i=1:jv(i)=A(j,i)*A(i,i);endA(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1,1);A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1,1)/A(j,j);endD=diag(diag(A);L=tril(A(:,1:n),-1)+eye(n);使用两种方法求解对称正定方程（1）系数矩阵为 100 阶矩阵clearclca=ones(100,1)*10;b=ones(99,1);c=ones(99,1);A=diag

23、(a)+diag(b,1)+diag(c,-1);n=length(A);b=rand(n,1);L=cholesky1(A);y=ones(n,1);y=Lb;x=ones(n,1);x=Ly% L D=cholesky2(A);% z=ones(n,1);% z=Lb;% y=ones(n,1);% y=Dz;% x=ones(n,1);% x=Ly（2）系数矩阵为 40 阶 Hilbert 矩阵clearclcn=40;for i=1:nfor j=1:nA(i,j)=1/(i+j-1);endendfor i=1:nb(i,1)=sum(A(i,1:i);endL=cholesky1(

24、A);y=ones(n,1);y=Lb;x=ones(n,1);x=Ly% L D=cholesky2(A);% z=ones(n,1);% z=Lb;% y=ones(n,1);% y=Dz;% x=ones(n,1);% x=Ly（3）使用两种方法求解对称正定方程实验一的程序求解这两个方程组系数矩阵为 100 阶矩阵，b 随机的选取clearclca=ones(100,1)*10;b=ones(99,1);c=ones(99,1);A=diag(a)+diag(b,1)+diag(c,-1);n=length(A);b=rand(n,1);L1,U1=gauss1(A,b);y1=ones

25、(n,1);y1=L1b;x1=ones(n,1);x1=U1y1;L2,U2,P=gauss2(A,b);y2=ones(n,1);y2=L2(P*b);x2=ones(n,1);x2=U2y2;x=x1,x2解为：0.0811099957553877 0.08110999575538770.0396127167351301 0.03961271673513010.0833223642481958 0.08332236424819580.0567725075395750 0.05677250753957500.0456197609112823 0.04561976091128230.0698

26、208485234418 0.06982084852344180.0715689653317207 0.07156896533172070.0935034027565292 0.0935034027565292-0.0176913768174234 -0.01769137681742340.0839327407746498 0.08393274077464980.0438025600839501 0.04380256008395010.0906081278698479 0.09060812786984790.0400663669264014 0.04006636692640140.036408

27、2722045803 0.03640827220458030.0753742962380143 0.07537429623801430.0111963709372294 0.01119637093722940.0405049300957332 0.04050493009573320.0818486193018287 0.08184861930182870.0418613654179844 0.04186136541798440.0741989456485145 0.07419894564851450.0613273631509068 0.06132736315090680.0511677148

28、378195 0.05116771483781950.0129825242973739 0.01298252429737390.0657415681744167 0.0657415681744167-0.00398198872207325 -0.003981988722073250.0575611326489384 0.05756113264893840.0543304474042721 0.05433044740427210.0600789512556829 0.06007895125568290.0746318953561201 0.07463189535612010.0843542115

29、084385 0.08435421150843850.0641292124431014 0.06412921244310140.0433827493964434 0.04338274939644340.0834897814678624 0.08348978146786240.0500324982391206 0.0500324982391206-0.00372439810062703 -0.003724398100627030.00419442110441090 0.004194421104410900.0826397581550762 0.08263975815507620.03211871

30、60444968 0.03211871604449680.0804695926120587 0.08046959261205870.00804103241118003 0.008041032411180030.0485251672970758 0.04852516729707580.0589986361568368 0.0589986361568368-0.00862814380102282 -0.008628143801022820.0592738176159583 0.05927381761595830.0306033867585804 0.0306033867585804-0.00289

31、622292871012 -0.002896222928710120.0478914215705819 0.04789142157058190.0135519964002127 0.01355199640021270.00909901048936550 0.009099010489365500.0185416462520775 0.01854164625207750.0109786978975394 0.01097869789753940.0181862853874183 0.0181862853874183-0.00376937729910872 -0.003769377299108720.

32、0621598985148123 0.06215989851481230.0173683090108681 0.01736830901086810.0460238672569371 0.04602386725693710.0609896964651019 0.06098969646510190.0392422075363764 0.03924220753637640.0457042416537237 0.04570424165372370.0395164316774994 0.03951643167749940.00431460686732418 0.004314606867324180.04

33、12697772473289 0.04126977724732890.0733449141274049 0.07334491412740490.0782792368194384 0.07827923681943840.0177901235399449 0.01779012353994490.0141138600738098 0.01411386007380980.0495326344732706 0.04953263447327060.0555393659316849 0.05553936593168490.0353855313726384 0.03538553137263840.007634

34、27198481669 0.007634271984816690.0942472643114380 0.0942472643114380-0.00217379380602789 -0.002173793806027890.00956188084656682 0.009561880846566820.0122644119220811 0.01226441192208110.00983512183662033 0.009835121836620330.0558448105881364 0.05584481058813640.0526754162173245 0.0526754162173245-0

35、.00888920792018309 -0.008889207920183090.0882945532703761 0.08829455327037610.0571450598246723 0.05714505982467230.0689165301611718 0.0689165301611718-0.00846870763880048 -0.008468707638800480.0791750469196512 0.07917504691965120.0771588014805202 0.07715880148052020.0836420572363599 0.08364205723635

36、990.0708189383968527 0.07081893839685270.0671073754789800 0.06710737547898000.0436662960783785 0.04366629607837850.00960708232481019 0.009607082324810190.0378653411793847 0.03786534117938470.0103290026171858 0.0103290026171858-0.00722411636327180 -0.007224116363271800.0928017097604837 0.092801709760

37、48370.0183487248279831 0.01834872482798310.0250171065460780 0.02501710654607800.0270140441865938 0.02701404418659380.0377787334241599 0.03777873342415990.0622668086006599 0.0622668086006599-0.0122484129646017 -0.01224841296460170.0854455025383936 0.0854455025383936两种方法的解是相同的。系数矩阵为 40 阶 Hilbert 矩阵cle

38、arclcn=40;for i=1:nfor j=1:nA(i,j)=1/(i+j-1);endendfor i=1:nb(i,1)=sum(A(i,1:i);endL1,U1=gauss1(A,b);y1=ones(n,1);y1=L1b;x1=ones(n,1);x1=U1y1;L2,U2,P=gauss2(A,b);y2=ones(n,1);y2=L2(P*b);x2=ones(n,1);x2=U2y2;x=x1,x2解为：-9.58351425349451 0.4659068886998611682.71046094840 0.160278216644432-69474.4327657

39、996 0.6821769420344461217187.32961003 1.30078576824103-11019964.2658012 1.7165929853842955145062.2876976 1.79424475284066-144412699.664341 1.51265622719715105740127.211396 0.924709225049974448345012.604736 0.124921502308529-1367974445.83320 -0.7751288988507831439837162.77148 -1.66448050957300-205484

40、678.547075 -2.44466823987476-785518769.385150 -3.03692156421836581622696.264047 -3.38570949019895508452075.087937 -3.45953515110102-2193142960.61393 -3.249747700036192748141030.24443 -2.76800810135934-2133814708.44930 -2.042924791003932337988306.63981 -1.11626769843022-1514395741.54634 -0.0390755337

41、70914458860687.4853639 1.13210875052064-618484046.464462 2.33870666132992706258577.585647 3.52316498750377-224851400.462300 4.63154192570550144185780.668900 5.61558776795708-118116092.865337 6.434317543043211667484676.49910 7.05509432679855-1737458598.78151 7.45425748969768-567709464.670292 7.617340

42、87629818537496130.366378 7.538932652310172116002843.20996 7.22223205645186-2468219728.20071 6.67835920730599-911417410.448053 5.925473005759592305367098.68411 4.98774952660395-1535255766.31145 3.894269522052141542851149.18856 2.67785911234374-468128587.158900 1.37392270489563-13694657.5457891 0.0193

43、019005578695-747502864.596306 -2949.23188216782461675146.163823 2970.90076474062两种方法的解完全不同，用选主元的方法的解更加收敛。实验三实验名称：直接法的时间复杂性试验。实验目的：分别用三种不同方法求解线性方程组 Ax=b，不同工作量得出不同时间。实验内容与要求：生成方程组 中矩阵 和右端项 ，分别用 ， 和三角分bAxbxAb()*invb解法计算，并分别记录所花费的 CPU 时间，进行分析比较。实验要求：（1）取 ，随机生成 的一条主对角线元素和二条次对角线元素，使 为严30nAA格对角占优的三对角阵和 ；其中

44、三条对角线元素分别用三个一维数组存储；b（2）用 Matlab 语言自编 M 文件分别用 ， 和追赶法计算这三xb()*invAb对角方程组；并分别记录所花费的 CPU 时间；(3) 分析结果，得出你的结论。 xAbfunction x1,t1=cpu1(A,b)tic;x1=Ab;t1=toc ()*xinvAbfunction x2,t2=cpu2(A,b)tic;x2=inv(A)*b;t2=toc追赶法function x3,t3=cpu3(A,b)n=length(A);% A=A,b;tic;for k=1:n-1A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1

45、:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);endL=tril(A(:,1:n),-1)+eye(n);U=triu(A);y=ones(n,1);y=Lb;x3=ones(n,1);x3=Uy;t3=toc取 n=300，随机生成 A 的一条主对角线元素和二条次对角线元素，使 A 为严格对角占优的三对角阵和 b；其中三条对角线元素分别用三个一维数组存储；clearclca=rand(300,1)*10;b=rand(299,1);c=rand(299,1);A=diag(a)+diag(b,1)+diag(c,-1);b=rand(300,1);x1,t1=cpu1(A,b);x2,t2=cpu2(A,b);x3,t3=cpu3(A,b);可以得出以下结论：t1 =0.0431t2 =0.2868t3 =0.0890可见 x=Ab 高斯消去法用时最短，然后是追赶法，最慢的是 x=inv(A)*b。