1、p=input(计算条件数使用p-范数,p=cond_A=cond(A,p)m,n=size(A);Ab=A b;r=input(选主元方式(0:自动;1:手动),r=Abfor i=1:n-1 switch r case(0) aii,ip=max(abs(Ab(i:n,i); ip=ip+i-1; case (1) ip=input(第,num2str(i),步消元,请输入第列所选元素所处的行数:); end; Ab(i ip,:)=Ab(ip i,: aii=Ab(i,i); for k=i+1:n Ab(k,i:n+1)=Ab(k,i:n+1)-(Ab(k,i)/aii)*Ab(i,i
2、:n+1); if r=1 Ab endend; x=zeros(n,1);x(n)=Ab(n,n+1)/Ab(n,n);for i=n-1:-1:1 x(i)=(Ab(i,n+1)-Ab(i,i+1:n)*x(i+1:n)/Ab(i,i);endx运行结果(1)n=10,矩阵的条件数及自动选主元Cond(A,1) =103Cond(A,2) = Cond(A,inf) =程序自动选择主元(列主元)a.输入数据n=10计算条件数使用p-范数,p=1手动),r=0b.计算结果x=1,1,1,1,1,1,1,1,1,1T(2)n=10,手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素
3、作为主元手动),r=1第1步消元,请输入第1列所选元素所处的行数:第2步消元,请输入第2列所选元素所处的行数:2(实际选择时,第k步选择主元处于第k行)最终计算得x=, , , , , , , , , Tb. 每步消去过程总选取按模最大的元素作为主元3(实际选择时,第k步选择主元处于第k+1行)(3)n=20,手动选主元n=20x=,Tx=1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1T(4)A分别为幻方矩阵,Hilbert矩阵,pascal矩阵和随机矩阵将计算结果列于下表:11阶幻方矩阵10阶Hilbert矩阵选主元方式模最大模最小x1x2x3x4x5x6x7
4、x8x9x1010阶Pascal方矩阵10阶随机矩阵简要分析计算(1)表明:对于同一矩阵,不同范数定义的条件数是不同的;Gauss消去法在消去过程中选择模最大的主元能够得到比较精确的解。计算(2)表明:通过比较每次选取模最大或模最小的元素的选主元方式,可以发现,在本题给定的问题中,选取模最大的元素作为主元比取模最小的元素作为主元时产生的结果更精确。因为这样做可以避免使用较小的数做为除数,以免发生结果数量级加大,使大数吃掉小数,产生舍入误差。计算(3)表明:首先,n=20得到与计算(2),即n=10时一样的结论,即选取模最大的元素作为主元比取模最小的元素作为主元时产生的结果更精确;其次,与计算(
5、2) (Cond(A1010,1) =103)比较,Cond(A2020,1) =106显着增大,且计算(3)的误差也远大于计算(2),即,矩阵的条件数越大,产生的误差也越大。计算(4)表明: Gauss消去法在消去过程中,主元的选择与算法的稳定性有密切的联系。一般来说,选取绝对值大的元素作为主元比绝对值小的元素作为主元时的计算结果更加精确。但这并不是绝对的,一些特殊的方阵,如Pascal方矩阵,则恰恰是选择模最小的元素作为主元时计算结果最精确(选模最小的元素只是一个表象,这种选主元方法优于其他选主元方法的本质是这种选择方法能使消去过程不产生浮点数,而全是整数运算,只有在回代过程中才有可能会产
6、生浮点数)。在系数矩阵性质未知,或者说对于绝大多数的系数矩阵来说,选择模最大的元素作为主元是一种比较稳定和精确的方法。实验2.病态的线性方程组的求解理论的分析表明,求解病态的线性方程组是困难的。实际情况是否如此,会出现怎样的现象呢?考虑方程组Hx=b的求解,其中系数矩阵H为Hilbert矩阵,这是一个着名的病态问题。通过首先给定解(例如取为各个分量均为1)再计算出右端b的办法给出确定的问题。(1)选择问题的维数为6,分别用Gauss消去法、J迭代法、GS迭代法和SOR迭代法求解方程组,其各自的结果如何将计算结果与问题的解比较,结论如何(2)逐步增大问题的维数,仍然用上述的方法来解它们,计算的结
7、果如何计算的结果说明了什么(3)讨论病态问题求解的算法clearn=6; % Hilbert矩阵的阶数H=hilb(n);cond_H=cond(H);b=H* ones(n,1); D=diag(diag(H);U=-triu(H,1);L=-tril(H,-1);w=;tol=;maxk=10000;Ab=H b;% Gauss消去法n-1 aii,ip=max(abs(Ab(i: ip=ip+i-1; if abs(aii)= tol)&(kmaxk) k=k+1; x(:,k)=BJ*x(:,k-1)+fJ;k_J=kx_J=x(:,k)% GS迭代法BG=(D-L)U;fG=(D-L
8、)b;,2)=BG*x(:,1)+fG;,k)=BG*x(:,k-1)+fG;k_G=kx_G=x(:% SOR迭代法lw=(D-w*L)(1-w)*D+w*U);fS=(D-w*L)b*w;,2)=lw*x(:,1)+fS;,k)=lw*x(:,k-1)+fS;k_SOR=kx_SOR=x(:运行结果及简要分析(1)6阶Hilbert矩阵计算方法GaussJ法GS法SOR法迭代矩阵的1-范数-迭代矩阵的谱半径迭代次数48966659951InfNaNGauss消元法得到了较为精确的解。以作为收敛的标准:J法迭代矩阵的谱半径大于1,迭代不收敛;GS法和SOR法迭代矩阵的谱半径都略小于1,迭代是
9、收敛的,但收敛速度非常慢,在很多次迭代之后仍与精确解有一定误差,GS收敛速度略快于SOR法。(2)n阶Hilbert矩阵计算从6阶到25阶的Hilbert矩阵,为观察收敛速度,以作为收敛的标准。计算结果如下(仅列出有代表性的6、7、10、11、12、13、21、25阶计算结果)6阶k75160774355932537X71034917601878X8X9X101133216061589X111231814471367X121330613071197X132125225862893X14X15X16X17X18X19X20X212523722532240X22X23X24X25得到的主要结论如下
10、:(1)Gauss消去法:Gauss消去法求得的解与精确解的误差随Hilbert矩阵阶数的增加而增加,Hilbert矩阵阶数不大于11时,误差较小(小于1%),对于要求不要高的工程问题,这样的误差可以接受。当阶数大于11时,误差迅速增加,当阶数为13时,误差已经超过100%,一般来说这样的近似解不可接受。当阶数为25时,误差达到38900%。即低阶Hilbert矩阵可用Gauss消元法求解。(2) Jacobi迭代方法:无论Hilbert矩阵为多少阶,Jacobi迭代矩阵的谱半径都大于1,迭代不稳定、不收敛。(3) GS迭代方法:GS迭代矩阵的谱半径略小于1,迭代收敛,但收敛速度非常缓慢。(4
11、) SOR迭代方法:取w=,SOR迭代矩阵的谱半径略小于1,迭代收敛,但收敛速度亦非常缓慢。在本实验中,多数情况SOR法较GS迭代更慢一些(收敛步数更大)。从上面的结果可以看出:Hilbert矩阵阶数较小时,可用Gauss消元法直接求解,解的精度比迭代法高;随着阶数增加,Gauss消元法的精度迅速下降,解变得不可靠,这时,有一些迭代法,如GS法和SOR法仍可继续求得比较精确的解。另外,在三种迭代法中,GS和SOR法相对Jacobi法更有优势,但这两种方法的迭代矩阵谱半径已经非常接近1(病态问题),收敛速度都很慢。(3)病态问题的求解求解病态问题,主要的方法是对原方程进行预处理,以降低系数矩阵的条件数。例如选择非奇异矩阵,使方程组化为等价方程组,原方程的解。原则上应使矩阵的条件数比有所改善。一般和可选择为三角形矩阵或对角矩阵。理论上最好选择对角阵和,满足:。对于1100阶的Hilbert矩阵,假设,从设为由Hilbert矩阵对角元素组成的对角阵,设是中,使得最小的一种。的结果如下图所示:虽然这里取的并不是最优的D矩阵,但在上图中,可观察到经过预处理后的DHD的条件数相比原Hilbert矩阵减小了,在一定程度上改善了原Hilbert矩阵的性质。1515