数值分析实验报告_清华大学__线性代数方程组的数值解法

线性代数方程组的数值解法 实验1.主元的选取与算法的稳定性 问题提出：Gauss消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss消去法作为数值算法的稳定性呢？Gauss消去法从理论算法到数值算法，其关键是主元的选择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。 实验内容：考虑线性方程组 编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss消去过程。 实验要求： （1）取矩阵 ，则方程有解 。取n=10计算矩阵的条件数。让程序自动选取主元，结果如何？ （2）现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。 （3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。 （4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。重复上述实验，观察记录并分析实验结果。 1.1程序清单 n=input('矩阵A的阶数:n='); A=6*diag(ones(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='); Ab for 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],:)=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 end end;    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); end x 1.2运行结果 (1)n=10，矩阵的条件数及自动选主元 Cond(A,1) =2.5575×103 Cond(A,2) =1.7276×103 Cond(A,inf) =2.5575×103 程序自动选择主元（列主元） a.输入数据 矩阵A的阶数:n=10 计算条件数使用p-范数,p=1 选主元方式(0:自动;1:手动)，r=0 b.计算结果 x=[1，1，1，1，1，1，1，1，1，1]T (2)n=10，手动选主元 a.每步消去过程总选取按模最小或按模尽可能小的元素作为主元 矩阵A的阶数:n=10 计算条件数使用p-范数,p=1 选主元方式(0:自动;1:手动)，r=1 第1步消元，请输入第1列所选元素所处的行数：1 第2步消元，请输入第2列所选元素所处的行数：2 …（实际选择时，第k步选择主元处于第k行） 最终计算得 x=[1., 1., 1., 1., 0., 1., 0., 1., 0., 1.]T b.每步消去过程总选取按模最大的元素作为主元 矩阵A的阶数:n=10 计算条件数使用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]T (3)n=20，手动选主元 a.每步消去过程总选取按模最小或按模尽可能小的元素作为主元 矩阵A的阶数:n=20 计算条件数使用p-范数,p=1 选主元方式(0:自动;1:手动)，r=1 第1步消元，请输入第1列所选元素所处的行数：1 第2步消元，请输入第2列所选元素所处的行数：2 …（实际选择时，第k步选择主元处于第k行） 最终计算得 x=[1.，1.，1.，1.，0.，1.，0.，1.，0.，1.，0.，1.，0.，1.，0.，1.，0.，1.，0.，1.]T b.每步消去过程总选取按模最大的元素作为主元 矩阵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，1]T (4)A分别为幻方矩阵，Hilbert矩阵，pascal矩阵和随机矩阵 将计算结果列于下：   11阶幻方矩阵 10阶Hilbert矩阵 选主元方式 模最大 模最小 模最大 模最小 x1 1. 3. 0. 0. x2 1. -1. 1. 1. x3 1. 1. 0. 0. x4 0. 0. 1. 1. x5 1. 1. 0. 0. x6 1. 1. 1. 1. x7 1. 0. 0. 0. x8 1. 1. 1. 1. x9 1. 1. 0. 0. x10 0. 1. 1. 1.             10阶Pascal方矩阵 10阶随机矩阵 选主元方式 模最大 模最小 模最大 模最小 x1 1. 1. 1. 0. x2 0. 1. 1. 0. x3 1. 1. 1. 0. x4 0. 1. 1. 1. x5 1. 1. 1. 0. x6 0. 1. 0. 1. x7 1. 1. 1. 0. x8 0. 1. 1. 1. x9 1. 1. 0. 0. x10 1. 1. 1. 1.           1.3简要分析 计算(1)表明：对于同一矩阵，不同范数定义的条件数是不同的；Gauss消去法在消去过程中选择模最大的主元能够得到比较精确的解。 计算(2)表明：通过比较每次选取模最大或模最小的元素的选主元方式，可以发现，在本题给定的问题中，选取模最大的元素作为主元比取模最小的元素作为主元时产生的结果更精确。因为这样做可以避免使用较小的数做为除数，以免发生结果数量级加大，使大数吃掉小数，产生舍入误差。 计算(3)表明：首先，n=20得到与计算(2),即n=10时一样的结论，即选取模最大的元素作为主元比取模最小的元素作为主元时产生的结果更精确；其次，与计算(2) （Cond(A10×10,1) =2.5575×103）比较，Cond(A20×20,1) =2.6214×106显著增大，且计算(3)的误差也远大于计算(2)，即，矩阵的条件数越大，产生的误差也越大。 计算(4)表明： Gauss消去法在消去过程中，主元的选择与算法的稳定性有密切的联系。一般来说，选取绝对值大的元素作为主元比绝对值小的元素作为主元时的计算结果更加精确。但这并不是绝对的，一些特殊的方阵，如Pascal方矩阵，则恰恰是选择模最小的元素作为主元时计算结果最精确（选模最小的元素只是一个表象，这种选主元方法优于其他选主元方法的本质是这种选择方法能使消去过程不产生浮点数，而全是整数运算，只有在回代过程中才有可能会产生浮点数）。在系数矩阵性质未知，或者说对于绝大多数的系数矩阵来说，选择模最大的元素作为主元是一种比较稳定和精确的方法。 实验2.病态的线性方程组的求解 问题提出：理论的分析表明，求解病态的线性方程组是困难的。实际情况是否如此，会出现怎样的现象呢？ 实验内容：考虑方程组Hx=b的求解，其中系数矩阵H为Hilbert矩阵， 这是一个著名的病态问题。通过首先给定解（例如取为各个分量均为1）再计算出右端b的办法给出确定的问题。 实验要求： （1）选择问题的维数为6，分别用Gauss消去法、J迭代法、GS迭代法和SOR迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比较，结论如何？ （2）逐步增大问题的维数，仍然用上述的方法来解它们，计算的结果如何？计算的结果说明了什么？ （3）讨论病态问题求解的算法 2.1程序清单 clear n=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=1.5; tol=1.0e-6; 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)<1.0e-20 disp('Singular matrix, stop!'); pause;  % 如果主元绝对值小于限值，则认为A奇异，终止 end 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; end; 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); end x_Gauss=x % Jacobi迭代法 BJ=D\(L+U); fJ=D\b; x(:,1)=1.2*ones(n,1); x(:,2)=BJ*x(:,1)+fJ; k=2; while (norm(x(:,k)-x(:,k-1))>=tol)&&(k=tol)&&(k=tol)&&(k工程

路基工程安全技术交底工程项目施工成本控制工程量增项单年度零星工程技术标正投影法基本原理