数值分析上机作业 .doc

数值实验数值实验综述：线性代数方程组的解法是一切科学计算的基础与核心问题。求解方法大致可分为直接法和迭代法两大类。直接法指在没有舍入误差的情况下经过有限次运算可求得方程组的精确解的方法，因此也称为精确法。当系数矩阵是方的、稠密的、无任何特殊结构的中小规模线性方程组时，Gauss消去法是目前最基本和常用的方法。如若系数矩阵具有某种特殊形式，则为了尽可能地减少计算量与存储量，需采用其他专门的方法来求解。Gauss消去等同于矩阵的三角分解，但它存在潜在的不稳定性，故需要选主元素。对正定对称矩阵，采用平方根方法无需选主元。方程组的性态与方程组的条件数有关，对于病态的方程组必须采用特殊的方法进行求解。实验一一、 实验名称：矩阵的LU分解.二、 实验目的：用不选主元的LU分解和列主元LU分解求解线性方程组 Ax=b, 并比较这两种方法.三、 实验内容与要求（1）用所熟悉的计算机语言将不选主元和列主元LU分解编成通用的子程序，然后用编写的程序求解下面的84阶方程组将计算结果与方程组的精确解进行比较，并就此谈谈你对Gauss消去法的看法。（2）写出追赶法求解三对角方程组的过程，并编写程序求该实验中的方程组（1）不选主元高斯消去法求解方程组function L,U=gauss1(A,b)n=length(A);for k=1:n-1 A(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);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.000000000000001.000000000000000.9999999999999981.000000000000000.9999999999999931.000000000000010.9999999999999721.000000000000060.9999999999998861.000000000000230.9999999999995451.000000000000910.9999999999981811.000000000003640.9999999999927251.000000000014550.9999999999708981.000000000058200.9999999998835921.000000000232820.9999999995343671.000000000931270.9999999981374691.000000003725060.9999999925498741.000000014900250.9999999701994971.000000059601010.9999998807979861.000000238404030.9999995231919461.000000953616110.9999980927677831.000003814464440.9999923710711301.000015257857740.9999694842845201.000061031430960.9998779371380811.000244125723840.9995117485523231.000976502895350.9980469942092911.003906011581410.9921879768371871.015624046325570.9687519073490881.062496185300920.8750076294018071.249984741181830.5000305176945351.99993896437811-0.9998779278249614.99975585192486-6.9995116889494716.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.124999994185857.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);for k=1:n-1 m,p=max(abs(A(:,k); A(k,p,:)=A(p,k,:); P(k,p,:)=P(p,k,:); if A(k,k)=0 A(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;主程序为：clearclc 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 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.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000000.9999999999999991.000000000000000.9999999999999951.000000000000010.9999999999999791.000000000000040.9999999999999171.000000000000170.9999999999996671.000000000000670.9999999999986661.000000000002670.9999999999946641.000000000010670.9999999999786571.000000000042690.9999999999146291.000000000170740.9999999996585261.000000000682930.9999999986342271.000000002731210.9999999945389091.000000010916850.9999999781876491.000000043539330.9999999132628241.000000172108410.9999996612469361.000000655651090.9999987761179611.000002098083500.999997202555339这个病态方程的解没有发散，用选主元的方法可以有效的提高解的精度。（2）追赶法clearclc a=ones(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:n L(i)=c(i)/U(i-1); U(i)=a(i)-L(i)*e(i-1);endfor i=2:n y(1)=b(1); y(i)=b(i)-L(i)*y(i-1);endfor i=n-1:-1:1 x(n)=y(n)/U(n); x(i)=(y(i)-e(i)*x(i+1)/U(i);endt3=tocx=x解为：11.000000000000001.000000000000001.000000000000000.9999999999999981.000000000000000.9999999999999931.000000000000010.9999999999999721.000000000000060.9999999999998861.000000000000230.9999999999995451.000000000000910.9999999999981811.000000000003640.9999999999927251.000000000014550.9999999999708981.000000000058200.9999999998835921.000000000232820.9999999995343671.000000000931270.9999999981374691.000000003725060.9999999925498741.000000014900250.9999999701994971.000000059601010.9999998807979861.000000238404030.9999995231919461.000000953616110.9999980927677831.000003814464440.9999923710711301.000015257857740.9999694842845201.000061031430960.9998779371380811.000244125723840.9995117485523231.000976502895350.9980469942092911.003906011581410.9921879768371871.015624046325570.9687519073490881.062496185300920.8750076294018071.249984741181830.5000305176945351.99993896437811-0.9998779278249614.99975585192486-6.9995116889494716.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.124999994185857.49999999-8355328.9999999716645128.9999999-33028126.999999965007744.9999998-125821439.000000234866688.999999-402628606.999999536838144.999998实验二一、 实验名称：实对称正定矩阵的的Cholesky分解.二、 实验目的：用平方根法和改进的平方根方法求解线性方程组 Ax=b.三、 实验内容与要求用所熟悉的计算机语言将Cholesky分解和改进的Cholesky分解编成通用的子程序，然后用编写的程序求解对称正定方程组Ax=b，其中（1） b随机的选取，系数矩阵为100阶矩阵（2） 系数矩阵为40阶Hilbert矩阵，即系数矩阵A的第i行第j列元素为，向量b的第i个分量为（3） 用实验一的程序求解这两个方程组，并比较所有的计算结果，然后评价各个方法的优劣。平方根方法求解function L=cholesky1(A) n=length(A);tic;for k=1:n A(k,k)=sqrt(A(k,k); A(k+1:n,k)=A(k+1:n,k)/A(k,k); for j=k+1:n A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k); endendt1=tocL=tril(A);改进的平方根方法function L,D=cholesky2(A)n=length(A);v=ones(n,1);for j=1:n for i=1:j v(i)=A(j,i)*A(i,i); end A(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阶矩阵clearclc a=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);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:n for j=1:n A(i,j)=1/(i+j-1); endendfor i=1:n b(i,1)=sum(A(i,1:i);endL=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（3）使用两种方法求解对称正定方程实验一的程序求解这两个方程组系数矩阵为100阶矩阵，b随机的选取clearclc a=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(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.08110999575538770.08110999575538770.03961271673513010.03961271673513010.08332236424819580.08332236424819580.05677250753957500.05677250753957500.04561976091128230.04561976091128230.06982084852344180.06982084852344180.07156896533172070.07156896533172070.09350340275652920.0935034027565292-0.0176913768174234-0.01769137681742340.08393274077464980.08393274077464980.04380256008395010.04380256008395010.09060812786984790.09060812786984790.04006636692640140.04006636692640140.03640827220458030.03640827220458030.07537429623801430.07537429623801430.01119637093722940.01119637093722940.04050493009573320.04050493009573320.08184861930182870.08184861930182870.04186136541798440.04186136541798440.07419894564851450.07419894564851450.06132736315090680.06132736315090680.05116771483781950.05116771483781950.01298252429737390.01298252429737390.06574156817441670.0657415681744167-0.00398198872207325-0.003981988722073250.05756113264893840.05756113264893840.05433044740427210.05433044740427210.06007895125568290.06007895125568290.07463189535612010.07463189535612010.08435421150843850.08435421150843850.06412921244310140.06412921244310140.04338274939644340.04338274939644340.08348978146786240.08348978146786240.05003249823912060.0500324982391206-0.00372439810062703-0.003724398100627030.004194421104410900.004194421104410900.08263975815507620.08263975815507620.03211871604449680.03211871604449680.08046959261205870.08046959261205870.008041032411180030.008041032411180030.04852516729707580.04852516729707580.05899863615683680.0589986361568368-0.00862814380102282-0.008628143801022820.05927381761595830.05927381761595830.03060338675858040.0306033867585804-0.00289622292871012-0.002896222928710120.04789142157058190.04789142157058190.01355199640021270.01355199640021270.009099010489365500.009099010489365500.01854164625207750.01854164625207750.01097869789753940.01097869789753940.01818628538741830.0181862853874183-0.00376937729910872-0.003769377299108720.06215989851481230.06215989851481230.01736830901086810.01736830901086810.04602386725693710.04602386725693710.06098969646510190.06098969646510190.03924220753637640.03924220753637640.04570424165372370.04570424165372370.03951643167749940.03951643167749940.004314606867324180.004314606867324180.04126977724732890.04126977724732890.07334491412740490.07334491412740490.07827923681943840.07827923681943840.01779012353994490.01779012353994490.01411386007380980.01411386007380980.04953263447327060.04953263447327060.05553936593168490.05553936593168490.03538553137263840.03538553137263840.007634271984816690.007634271984816690.09424726431143800.0942472643114380-0.00217379380602789-0.002173793806027890.009561880846566820.009561880846566820.01226441192208110.01226441192208110.009835121836620330.009835121836620330.05584481058813640.05584481058813640.05267541621732450.0526754162173245-0.00888920792018309-0.008889207920183090.08829455327037610.08829455327037610.05714505982467230.05714505982467230.06891653016117180.0689165301611718-0.00846870763880048-0.008468707638800480.07917504691965120.07917504691965120.07715880148052020.07715880148052020.08364205723635990.08364205723635990.07081893839685270.07081893839685270.06710737547898000.06710737547898000.04366629607837850.04366629607837850.009607082324810190.009607082324810190.03786534117938470.03786534117938470.01032900261718580.0103290026171858-0.00722411636327180-0.007224116363271800.09280170976048370.09280170976048370.01834872482798310.01834872482798310.02501710654607800.02501710654607800.02701404418659380.02701404418659380.03777873342415990.03777873342415990.06226680860065990.0622668086006599-0.0122484129646017-0.01224841296460170.08544550253839360.0854455025383936两种方法的解是相同的。系数矩阵为40阶Hilbert矩阵clearclcn=40;for i=1:n for j=1:n A(i,j)=1/(i+j-1); endendfor i=1:n b(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.583514253494510.4659068886998611682.710460948400.160278216644432-69474.43276579960.6821769420344461217187.329610031.30078576824103-11019964.26580121.7165929853842955145062.28769761.79424475284066-144412699.6643411.51265622719715105740127.2113960.924709225049974448345012.6047360.124921502308529-1367974445.83320-0.7751288988507831439837162.77148-1.66448050957300-205484678.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.039075533770914458860687.48536391.13210875052064-618484046.4644622.33870666132992706258577.5856473.52316498750377-224851400.4623004.63154192570550144185780.6689005.61558776795708-118116092.8653376.4343175430432116