1、-本页仅作为文档封面，使用时请直接删除即可-内页可以根据需求调整合适字体及大小-数值分析实验报告-清华大学-线性代数方程组的数值解法(总15页)线性代数方程组的数值解法实验1. 主元的选取与算法的稳定性问题提出：Gauss消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss消去法作为数值算法的稳定性呢？Gauss消去法从理论算法到数值算法，其关键是主元的选择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。实验内容：考虑线性方程组 编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss消去过程。实验要
2、求：（1）取矩阵，则方程有解。取n=10计算矩阵的条件数。让程序自动选取主元，结果如何？（2）现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。（4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。重复上述实验，观察记录并分析实验结果。程序清单n=input(矩阵A的阶数:n=);A=6*diag(ones
3、(1,n)+diag(ones(1,n-1),1)+8*diag(ones(1,n-1),-1);b=A*ones(n,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),步消元，请输入第,num2str(i),列所选元素所处的行数：); end; Ab(i ip,:)=
4、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:n+1); end; 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) = 103Cond(A,inf) =103程序自动选择主元（列主元）a.输入数据矩
5、阵A的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动)，r=0b.计算结果x=1，1，1，1，1，1，1，1，1，1T(2)n=10，手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元矩阵A的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动)，r=1第1步消元，请输入第1列所选元素所处的行数：1第2步消元，请输入第2列所选元素所处的行数：2（实际选择时，第k步选择主元处于第k行）最终计算得x=, , , , , , , , , Tb. 每步消去过程总选取按模最大的元素作为主元矩阵A的阶数:n=10计算条件数使用p-范数,
6、p=1选主元方式(0:自动;1:手动)，r=1第1步消元，请输入第1列所选元素所处的行数：2第2步消元，请输入第2列所选元素所处的行数：3（实际选择时，第k步选择主元处于第k+1行）最终计算得x=1，1，1，1，1，1，1，1，1，1T(3)n=20，手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元矩阵A的阶数:n=20计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动)，r=1第1步消元，请输入第1列所选元素所处的行数：1第2步消元，请输入第2列所选元素所处的行数：2（实际选择时，第k步选择主元处于第k行）最终计算得x=，Tb. 每步消去过程总选取按模最大的元
7、素作为主元矩阵A的阶数:n=20计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动)，r=1第1步消元，请输入第1列所选元素所处的行数：2第2步消元，请输入第2列所选元素所处的行数：3（实际选择时，第k步选择主元处于第k+1行）最终计算得x=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矩阵选主元方式模最大模最小模最大模最小x1x2x3x4x5x6x7x8x9x1010阶Pascal方矩阵10阶随机矩阵选主元方式模最大模最小模最
8、大模最小x1x2x3x4x5x6x7x8x9x10简要分析计算(1)表明：对于同一矩阵，不同范数定义的条件数是不同的；Gauss消去法在消去过程中选择模最大的主元能够得到比较精确的解。计算(2)表明：通过比较每次选取模最大或模最小的元素的选主元方式，可以发现，在本题给定的问题中，选取模最大的元素作为主元比取模最小的元素作为主元时产生的结果更精确。因为这样做可以避免使用较小的数做为除数，以免发生结果数量级加大，使大数吃掉小数，产生舍入误差。计算(3)表明：首先，n=20得到与计算(2),即n=10时一样的结论，即选取模最大的元素作为主元比取模最小的元素作为主元时产生的结果更精确；其次，与计算(2
9、) （Cond(A1010,1) =103）比较，Cond(A2020,1) =106显着增大，且计算(3)的误差也远大于计算(2)，即，矩阵的条件数越大，产生的误差也越大。计算(4)表明： Gauss消去法在消去过程中，主元的选择与算法的稳定性有密切的联系。一般来说，选取绝对值大的元素作为主元比绝对值小的元素作为主元时的计算结果更加精确。但这并不是绝对的，一些特殊的方阵，如Pascal方矩阵，则恰恰是选择模最小的元素作为主元时计算结果最精确（选模最小的元素只是一个表象，这种选主元方法优于其他选主元方法的本质是这种选择方法能使消去过程不产生浮点数，而全是整数运算，只有在回代过程中才有可能会产生
10、浮点数）。在系数矩阵性质未知，或者说对于绝大多数的系数矩阵来说，选择模最大的元素作为主元是一种比较稳定和精确的方法。实验2.病态的线性方程组的求解问题提出：理论的分析表明，求解病态的线性方程组是困难的。实际情况是否如此，会出现怎样的现象呢？实验内容：考虑方程组Hx=b的求解，其中系数矩阵H为Hilbert矩阵， 这是一个着名的病态问题。通过首先给定解（例如取为各个分量均为1）再计算出右端b的办法给出确定的问题。实验要求：（1）选择问题的维数为6，分别用Gauss消去法、J迭代法、GS迭代法和SOR迭代法求解方程组，其各自的结果如何将计算结果与问题的解比较，结论如何（2）逐步增大问题的维数，仍然
11、用上述的方法来解它们，计算的结果如何计算的结果说明了什么（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消去法for i=1:n-1 aii,ip=max(abs(Ab(i:n,i); ip=ip+i-1; Ab(i ip,:)=Ab(ip i,:); aii=Ab(i,i); if abs(aii)= tol)&(k= tol)
12、&(k= tol)&(kmaxk) k=k+1; x(:,k)=lw*x(:,k-1)+fS;endk_SOR=kx_SOR=x(:,k)运行结果及简要分析(1)6阶Hilbert矩阵计算方法GaussJ法GS法SOR法迭代矩阵的1-范数-迭代矩阵的谱半径-迭代次数-48966659951x1Infx2Infx3NaNx4NaNx5NaNx6NaNGauss消元法得到了较为精确的解。以作为收敛的标准：J法迭代矩阵的谱半径大于1，迭代不收敛；GS法和SOR法迭代矩阵的谱半径都略小于1，迭代是收敛的，但收敛速度非常慢，在很多次迭代之后仍与精确解有一定误差，GS收敛速度略快于SOR法。(2)n阶Hi
13、lbert矩阵计算从6阶到25阶的Hilbert矩阵，为观察收敛速度，以作为收敛的标准。计算结果如下（仅列出有代表性的6、7、10、11、12、13、21、25阶计算结果）计算方法GaussJ法GS法SOR法6阶-k-489751607x1Infx2Infx3NaNx4NaNx5NaNx6NaN7阶-k-4355932537x1Infx2NaNx3NaNx4NaNx5NaNx6NaNX7NaN10阶-k-34917601878x1Infx2NaNx3NaNx4NaNx5NaNx6NaNX7NaNX8NaNX9NaNX10NaN11阶-k-33216061589x1Infx2NaNx3NaNx4
14、NaNx5NaNx6NaNX7NaNX8NaNX9NaNX10NaNX11NaN12阶-k-31814471367x1Infx2NaNx3NaNx4NaNx5NaNx6NaNX7NaNX8NaNX9NaNX10NaNX11NaNX12NaN13阶-k-30613071197x1Infx2Infx3Infx4NaNx5NaNx6NaNX7NaNX8NaNX9NaNX10NaNX11NaNX12NaNX13NaN21阶-k-25225862893x1NaNx2NaNx3NaNx4NaNx5NaNx6NaNX7NaNX8NaNX9NaNX10NaNX11NaNX12NaNX13NaNX14NaNX1
15、5NaNX16NaNX17NaNX18NaNX19NaNX20NaNX21NaN25阶-k-23722532240x1NaNx2NaNx3NaNx4NaNx5NaNx6NaNX7NaNX8NaNX9NaNX10NaNX11NaNX12NaNX13NaNX14NaNX15NaNX16NaNX17NaNX18NaNX19NaNX20NaNX21NaNX22NaNX23NaNX24NaNX25NaN得到的主要结论如下：(1)Gauss消去法：Gauss消去法求得的解与精确解的误差随Hilbert矩阵阶数的增加而增加，Hilbert矩阵阶数不大于11时，误差较小（小于1%），对于要求不要高的工程问题，
16、这样的误差可以接受。当阶数大于11时，误差迅速增加，当阶数为13时，误差已经超过100%，一般来说这样的近似解不可接受。当阶数为25时，误差达到38900%。即低阶Hilbert矩阵可用Gauss消元法求解。(2) Jacobi迭代方法：无论Hilbert矩阵为多少阶，Jacobi迭代矩阵的谱半径都大于1，迭代不稳定、不收敛。(3) GS迭代方法：GS迭代矩阵的谱半径略小于1，迭代收敛，但收敛速度非常缓慢。(4) SOR迭代方法：取w=，SOR迭代矩阵的谱半径略小于1，迭代收敛，但收敛速度亦非常缓慢。在本实验中，多数情况SOR法较GS迭代更慢一些（收敛步数更大）。从上面的结果可以看出：Hilb
17、ert矩阵阶数较小时，可用Gauss消元法直接求解，解的精度比迭代法高；随着阶数增加，Gauss消元法的精度迅速下降，解变得不可靠，这时，有一些迭代法，如GS法和SOR法仍可继续求得比较精确的解。另外，在三种迭代法中，GS和SOR法相对Jacobi法更有优势，但这两种方法的迭代矩阵谱半径已经非常接近1（病态问题），收敛速度都很慢。(3)病态问题的求解求解病态问题，主要的方法是对原方程进行预处理，以降低系数矩阵的条件数。例如选择非奇异矩阵，使方程组化为等价方程组，原方程的解。原则上应使矩阵的条件数比有所改善。一般和可选择为三角形矩阵或对角矩阵。理论上最好选择对角阵和，满足：。对于1100阶的Hilbert矩阵，假设，从设为由Hilbert矩阵对角元素组成的对角阵，设是中，使得最小的一种。的结果如下图所示：虽然这里取的并不是最优的D矩阵，但在上图中，可观察到经过预处理后的DHD的条件数相比原Hilbert矩阵减小了，在一定程度上改善了原Hilbert矩阵的性质。1515
