下载文档原格式
《矩阵分析与数值分析》实验报告院系：姓名：学号：所在班号：任课老师：一．设错误！未找到引用源。
，分别编制从小到大和从大到小的顺序程序计算错误！未找到引用源。
并指出有效位数。
程序如下：function sum3j=input('请输入求和个数 "j":');A=0;B=0;double B;double A;for n=2:jm=n^2-1;t=1./m;A=A+t;enddisp('从小到大：')s=Afor n=j:-1:2m=n^2-1;t=1./m;B=B+t;enddisp('从大到小:')s=B运行结果：>> sum3请输入求和个数 "j":100从小到大：s =0.740049504950495从大到小:s =0.740049504950495>> sum3请输入求和个数 "j":10000从小到大：s =0.749900004999506从大到小:s =0.749900004999500>> sum3请输入求和个数 "j":1000000从小到大：s =0.749999000000522从大到小:s =0.749999000000500二、解线性方程组1．分别Jacobi 迭代法和Gauss-Seidel 迭代法求解线性方程组。
⎪⎪⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----000121001210012100124321x x x x 迭代法计算停止的条件为：6)()1(3110max -+≤≤<-k j k j j x x 。
解：（1）Jacobi 迭代法程序代码： function jacobi(A, b, N) clc;clear;A=[-2 1 0 0;1 -2 1 0;0 1 -2 1;0 0 1 -2]; b=[-1 0 0 0]'; N=100;n = size(A,1); D = diag(diag(A)); L = tril(-A,-1); U = triu(-A,1); Tj = inv(D)*(L+U); cj = inv(D)*b; tol = 1e-06; k = 1;format longx = zeros(n,1); while k <= Nx(:,k+1) = Tj*x(:,k) + cj;disp(k); disp('x = ');disp(x(:,k+1)); if norm(x(:,k+1)-x(:,k)) < toldisp('The procedure was successful')disp('Condition ||x^(k+1) - x^(k)|| < tol was met after k iterations') disp(k); disp('x = ');disp(x(:,k+1)); break endk = k+1; end结果输出The procedure was successfulCondition ||x^(k+1) - x^(k)|| < tol was met after k iterations 60 x =0.799998799067310.599998427958700.399998056850090.19999902842505（2）Gauss-Seidel迭代法程序代码：function gauss_seidel(A, b, N)clc;clear;A=[-2 1 0 0;1 -2 1 0;0 1 -2 1;0 0 1 -2];b=[-1 0 0 0]';N=100;n = size(A,1);D = diag(diag(A));L = tril(-A,-1);U = triu(-A,1);Tg = inv(D-L)*U;cg = inv(D-L)*b;tol = 1e-06;k = 1;x = zeros(n,1);while k <= Nx(:,k+1) = Tg*x(:,k) + cg;disp(k); disp('x = ');disp(x(:,k+1));if norm(x(:,k+1)-x(:,k)) < toldisp('The procedure was successful')disp('Condition ||x^(k+1) - x^(k)|| < tol was met after k iterations') disp(k); disp('x = ');disp(x(:,k+1));breakendk = k+1;end结果输出The procedure was successfulCondition ||x^(k+1) - x^(k)|| < tol was met after k iterations31x =0.799999213979350.599998971085610.399999167590770.199999583795392. 用Gauss列主元消去法、QR方法求解如下方程组：⎪⎪⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎭⎫ ⎝⎛---017435222331325212214321x x x x (1)Gauss 列主元消去法 程序代码：function x=Gaussmain(A,b) clc;clear; format longA=[1 2 1 2;2 5 3 -2;-2 -2 3 5;1 3 2 3]; b=[4 7 -1 0]'; N=length(A); x=zeros(N,1); y=zeros(N,1); c=0; d=0;A(:,N+1)=b; for k=1:N-1 for i=k:4if c<abs(A(i,k))d=i;c=abs(A(i,k)); end endy=A(k,:);A(k,:)=A(d,:); A(d,:)=y; for i=k+1:N c=A(i,k);for j=1:N+1A(i,j)=A(i,j)-A(k,j)*c/A(k,k); end end endb=A(:,N+1);x(N)=b(N)/A(N,N); for k=N-1:-1:1x(k)=(b(k)-A(k,k+1:N)*x(k+1:N))/A(k,k); end结果输出 ans =18.00000000000000 -9.571428571428576.00000000000000-0.42857142857143（2）QR方法：程序代码function QR(A,b)clc;clear;format longA=[1 2 1 2;2 5 3 -2;-2 -2 3 5;1 3 2 3];b=[4 7 -1 0]';[Q,R]=qr(A)y=Q'*b;x=R\y结果输出Q =-0.31622776601684 -0.04116934847963 -0.75164602800283 0.57735026918962 -0.63245553203368 -0.49403218175557 -0.15032920560056 -0.57735026918963 0.63245553203368 -0.74104827263336 -0.22549380840085 -0.00000000000000 -0.31622776601684 -0.45286283327594 0.60131682240226 0.57735026918963 R =-3.16227766016838 -6.00832755431992 -0.94868329805051 2.84604989415154 0 -2.42899156029822 -4.65213637819829 -4.15810419644272 0 0 -0.67648142520255 -0.52615221960200 0 0 0 4.04145188432738 x =17.99999999999989-9.571428571428515.99999999999997-0.42857142857143三、非线性方程的迭代解法1．用Newton迭代法求方程()06cos22x=-++=-xexf x的根，计算停止的条件为：6110-+<-kkxx；编程如下：function newton(f,df,x,a,a0)syms xf=input('please enter your equation:') a0=input('please enter you x(0):');df=diff(f)e=1e-6;a1=a0+1;N=0;while abs(a1-a0)>ea=a0-subs(f,a0)/subs(df,a0); a1=a0; a0=a; N=N+1; endfprintf('a=%0.6f',a) N运行结果： >> newtonplease enter your equation:exp(x)+2^(-x)+2*cos(x)-6 f =exp(x)+2^(-x)+2*cos(x)-6 please enter you x(0):2df =exp(x)-2^(-x)*log(2)-2*sin(x) a=1.829384 N =42．利用Newton 迭代法求多项式07951.2954.856.104.5x 234=+-+-x x x的所有实零点，注意重根的问题。
Hilbert 代数方程组的数值解法1 Hilbert 矩阵的条件数和矩阵的阶数的关系编制Matlab 程序：clear,clc format long %format short for n=1:14; A=hilb(n); %A=rand(n);condA(n)=cond(A,inf); endplot(log10(condA)) 得出图形图1：Hilbert 矩阵和阶次的关系由图可见a.Hilbert 矩阵的2-条件数的对数与阶次近似正比;b.阶次超过12后，计算困难，舍入误差会导致计算不准 结论：Hilbert 矩阵的2-条件数随阶次指数增长，关系大概是： Log 10(Cond(A))= 1.33791720780254(n-1)0246810121424681012141618nl o g 10(c o n d A )条件数的对数-阶次2.各种求解方法的对比2.1 直接法：Gauss 消去方法（程序见附录，下同）在双精度型变量下，即eps=2.220446049250313e-016，对于直接法Gauss 消去方法，求解Hilbert 矩阵方程组，在阶次n=13时，解得: X=[1 0.99997 1.0013 0.97881 1.1906 -0.035441 4.6171 -7.3967 14.089 -12.542 9.9162 -2.3817 1.5624]出现明显错误，可见Gauss 消去方法对病态问题比较敏感。
求解阶次不高。
2.2 Jacobi 迭代法很遗憾，jakobi 迭代法的迭代矩阵的谱半径随阶次线性上升（程序见附录2），普遍大于1，计算结果发送，无法迭代出满意的结果。
需放弃图2：Jacobi 迭代矩阵的谱半径2.3 Gauss-Seidel 迭代与SOR 迭代求解先分析SOR 迭代矩阵的谱半径和松弛因子w 的关系0510152025510152025阶次谱半径Jacobi 迭代收敛性分析图3：SOR 迭代矩阵的谱半径和松弛因子w 的关系可以发现SOR 迭代和G-S 迭代的谱半径都普遍接近于1，因此松弛因子的选择对计算结果的影响不大。
（Hilbert矩阵）病态线性方程组的求解理论分析表明，数值求解病态线性方程组很困难。
考虑求解如下的线性方程组的求解1Hx=b,期中H是Hilbert矩阵,H=（h j）四，h j=,।,j=1,2，…，nij-11,估计矩阵的2条件数和阶数的关系2 .对不同的n,取x=（1,1，…，1）w n,分别用Gauss消去，Jacobi迭代，Gauss-seidel迭代，SOR迭代和共轲梯度法求解，比较结果。
3 .结合计算结果，试讨论病态线性方程组的求解。
第1小题：condition.m。
蟆1小题程序t1=20;%>数n=20x1=1:t1;y1=1:t1;fori=1:t1H=hilb(i);y1(i)=log(cond(H));endplot(x1,y1);xlabel('阶数n');ylabel('2-条件数的对数(log(cond(H))');title('2-条件数的对数(10g(cond(H))与阶数n的关系图’);t2=200;%J>数n=200x2=1:t2;y2=1:t2;fori=1:t2H=hilb(i);y2(i)=log(cond(H));endplot(x2,y2);xlabel('阶数n');ylabel('2-条件数的对数(log(cond(H))');title('2-条件数的对数(10g(cond(H))与阶数n的关系图’);画出Hilbert矩阵2-条件数的对数和阶数的关系n=200时n=20时马KiMn的美方四*--wficEJ-帛=«』从图中可以看出，1)在n小于等于13之前，图像近似直线log(cond(H))-1.519n-1.8332)在n大于13之后，图像趋于平缓，并在一定范围内上下波动，同时随着n的增加稍有上升的趋势第2小题：soke.m%蹴2小题主程序N=4000;xGauss=zeros(N,1);xJacobi=zeros(N,1);xnJ=zeros(N,1);xGS=zeros(N,1);xnGS=zeros(N,1);xSOR=zeros(N,1);xnSOR=zeros(N,1);xCG=zeros(N,1);xnCG=zeros(N,1);forn=1:N;x=ones(n,1);t=1.1;%初始值偏差x0=t*x;%迭代初始值e=1.0e-8;%^定的误差A=hilb(n);b=A*x;max=100000000000;%可能最大的迭代次数w=0.5;%SORt代的松弛因子G=Gauss(A,b);[J,nJ]=Jacobi(A,b,x0,e,max);[GS,nGS]=G_S(A,b,x0,e,max);[S_R,nS_R]=SOR(A,b,x0,e,max,w);[C_G,nC_G]=CG(A,b,x0,e,max);normG=norm(G'-x);xGauss(n)=normG;normJ=norm(J-x);nJ;xJacobi(n)=normJ;xnJ(n)=nJ;normGS=norm(GS-x);nGS;xGS(n)=normGS;xnGS(n)=nGS;normS_R=norm(S_R-x);nS_R;xSOR(n)=normS_R;xnSOR(n)=nS_R;normC_G=norm(C_G-x);nC_G;xCG(n)=normC_G;xnCG(n)=nC_G;endGauss.m%Gauss消去法functionx=Gauss(A,b)n=length(b);l=zeros(n,n);x=zeros(1,n);脸肖去过程fori=1:n-1forj=i+1:nl(j,i)=A(j,i)/A(i,i);fork=i:nA(j,k)=A(j,k)-l(j,i)*A(i,k);endb(j)=b(j)-l(j,i)*b(i);endend%回代过程x(n)=b(n)/A(n,n);fori=n-1:-1:1c=A(i,:).*x;x(i)=(b(i)-sum(c(i+1:n)))/A(i,i);endJacobi.m%Jacobi迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function[x,n]=Jacobi(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=D\(L+U);f=D\b;x=B*x0+f;n=1;whilenorm(x-x0)>ex0=x;x=B*x0+f;n=n+1;ifn>mdisp('Jacobi迭代次数过多，迭代可能不收敛‘);break;endendG_S.m%Gauss-Seidel迭代,x0表示迭代初值，e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function[x,n]=G_S(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-L)\U;f=(D-L)\b;x=B*x0+f;n=1;whilenorm(x-x0)>ex0=x;x=B*x0+f;n=n+1;ifn>mdisp('Gauss-Seidel迭代次数过多，迭代可能不收敛’)；break;endendSOR.m%SOR超松弛迭代，x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数,w松弛因子function[x,n]=SOR(A,b,x0,e,m,w)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-w*L)\((1-w)*D+w*U);f=(D-w*L)\b*w;x=B*x0+f;n=1;whilenorm(x-x0)>ex0=x;x=B*x0+f;n=n+1;ifn>mdisp('SOR超松弛迭代次数过多，迭代可能不收敛’)；break;endendCG.m%C献朝梯度法，x0表示迭代初值，e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function[x,n]=CG(A,b,x0,e,m)r=b-A*x0;p=r;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=1;whilenorm(r1)>ebelta=(r1'*r1)/(r'*r);p=r1+belta*p;r=r1;x0=x;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=n+1;ifn>mdisp('CG共朝梯度法迭代次数过多，迭代可能不收敛’);break;endend第2小题选择不同的阶数n,取x=(l,l,…RI分别使用Gauss消去，Jacobi迭代，Gauss-seidel迭代,SOR迭代和共挽梯度法(CG)求解，比较结果。
病态线性代数方程组的求解理论的分析表明，求解病态的线性代数方程组是困难的。
考虑方程组Hx = b 的求解，其中H 为Hilbert 矩阵，n n ij h H ⨯=)(，11-+=j i h ij ，n j i ,...,2,1,=1. 估计Hilbert 矩阵2-条件数与阶数的关系；2. 选择问题的不同维数，分别用Gauss 消去法，Jacobi 迭代，GS 迭代和SOR 迭代求解，比较结果；3. 讨论病态问题求解的算法。
解： 1、取Hilbert 矩阵阶数最高分别为n=20和n=100。
采用Hilbert 矩阵的2-条件数作为实验的比较对象，画出的曲线如下图所示：lg(())n cond H n从图中可以看出，在n ≤13之前，图像接近直线，在n ＞13之后，图像趋于平缓，在一定的范围内上下波动。
为了比较图像的线性部分，作出一条直线与已知曲线进行比较。
比较直线的关系式为：833.1519.1))(lg(-=n H cond n ，结果下图所示。
nl g (c o n d (H n ))lg(cond(Hn))~n 关系图从图2中可以看出，当n 较小时，n H cond n ~))(lg(之间近似满足线性关系。
当n 继续增大到100时，n H cond n ~))(lg(关系图下图所示：从图中可以看出，图像的走势符合在n=20时的猜想，在n 大于一定的值之后，图像趋于平缓，且在一定范围内震荡，同时又有一定上升趋势，但上升速度很慢。
2、选择不同的阶数n ，设方程组的精确解为xz=(1,1,…,1)T进行计算，用四种方法解x_Guass1、x_Jacobi1、x_GS1、x_SOR1对比表如下nl g (c o n d (H n ))lg(cond(Hn))~n关系图nl g (c o n d (H n ))lg(cond(Hn))~n 关系图Gauss消去法求解：选择问题的阶数为3~8时，用Gauss消去法求得的解与精确解一致，当阶数为9~14时，解开始出现偏差，而且n越大，偏差越大。
数值分析第一次上机实习报告——线性方程组直接解法一、问题描述设 H n = [h ij ] ∈ R n ×n 是 Hilbert 矩阵, 即11ij h i j =+- 对n = 2,3,4，…13,(a) 取11n n x R ⨯⎛⎫ ⎪=∈ ⎪ ⎪⎝⎭，及n n b H x =，用Gauss 消去法和Cholesky 分解方法来求解n n H y b =，看看误差有多大.(b) 计算条件数：2()n cond H(c) 使用某种正则化方法改善(a)中的结果.二、方法描述1. Gauss 消去法Gauss 消去法一般用于系数矩阵稠密且没有任何特殊结构的线性方程组。
设H =[h ij ]，y = (y 1,y 2,…，y n )T . 首先对系数矩阵H n 进行LU 分解，对于k=1,2，…n,交替进行计算：1111),,1,,1(),1,2,,k kj kj kr rj r k ik ik ir rk r kk u h l u j k k n l a l u i k k n u -=-=⎧=-=+⎪⎪⎨⎪=-=++⎪⎩∑∑…… 由此可以得到下三角矩阵L=[l ij ]和上三角矩阵U=[u ij ]. 依次求解方程组Ly=b 和Ux=y ，111,1,2,,1(),,1,,1i i i ir r r n i i ir r r i ii y b l y i n x y u x i n n u -==+⎧=-=⎪⎪⎨⎪=-=-⎪⎩∑∑…… 即可确定最终解。
2. Cholesky 分解法对于系数矩阵对称正定的方程组n n H y b =，存在唯一的对角元素为正数的下三角矩阵L ，使得H=LL T 。
因此，首先对矩阵H n 进行Cholesky 分解，即1122111()1()j jj jj jk k j ij ij ik jk k jj l h l l h l l l -=-=⎧=-⎪⎪⎨⎪=-⎪⎩∑∑ 1,i j n =+… L 的元素求出之后，依次求解方程组Ly=b 和L T x=y ，即1111111(),2,3,i i i ik k k ii b y l y b l y i n l -=⎧=⎪⎪⎨⎪=-=⎪⎩∑… 11(),1,2,n n nn n i i ki k k i nn y x l x y l x i n n l =+⎧=⎪⎪⎨⎪=-=--⎪⎩∑…1 由此求得方程组n n H y b =的解。
实验一病态线性代数方程组的求解1.估计Hilbert矩阵2-条件数与阶数的关系运行tiaojianshu.m 输入m=10 可以得到如下表的结果2.选择不同维数，分别用Guass消去（LU分解），Jacobi迭代，GS 迭代，SOR迭代求解，比较结果。
说明：Hx=b，H矩阵可以由matlab直接给出，为了设定参考解，我们先设x为分量全1的向量，求出b,然后将H和b作为已知量，求x，与设定的参考解对比。
对于Jacobi迭代，GS迭代，SOR迭代，取迭代初值x0为0向量，迭代精度eps=1.0e-6，迭代次数<100000, SOR迭代中w=1.2和0.8分别计算。
a. n=5b. n=8c. n=10d. n=15取不同的n值，得到如下结果:对于Guass法，可以看出来，随着n的增大，求解结果误差变大，这是因为随着n增大，系数矩阵的条件数变大，微小的扰动就容易造成很大的误差。
最后得不到精确解。
对于Jacobi迭代，计算结果为Inf，说明是发散的。
对于GS迭代和SOR迭代，结果是收敛的，但是可以看出迭代次数比较多，并且对于不同维数GS和SOR收敛速度不一样，有时候GS快，有时SOR快。
对SOR取不同的w迭代速度也不一样，存在一个最优的松弛因子w。
并且可以知道，迭代次数多少跟初值x0也有关系。
3.讨论病态问题求解的算法。
通过上面的实验分析，可以看出，求解病态矩阵的时候要小心，否则可能得不到所要求的精确度。
可以采用高精度运算，用双倍多倍字长，使得由于误差放大而损失若干有效数字位之后，还能保留一些有效位。
另外可以通过对原方程作某些预处理，降低系数矩阵的条件数，因为cond(aA)=cond(A),所以不能通过将每一个方程乘上相同的常数来达到这个目标，可考虑将矩阵的每一行和每一列分别乘上不同的常数，亦即找到可逆的对角阵D1和D2将方程组化为D1AD2y=D1b，x=D2y这称为矩阵的平衡问题，但是这样计算量比原问题本身要多。
数值分析课程实验报告
题目：病态线性方程组的求解
理论分析表明，数值求解病态线性方程组很困难。
考虑求解如下的线性方程组的求解
Hx = b ，期中H 是Hilbert 矩阵，()ij n n H h ⨯=，1
1
ij h i j =+-，i ，j = 1,2，…，n
1. 估计矩阵的2条件数和阶数的关系
2. 对不同的n ，取(1,1,
,1)n
x =∈
，分别用Gauss 消去，Jacobi 迭代，Gauss-seidel 迭
代，SOR 迭代和共轭梯度法求解，比较结果。
3. 结合计算结果，试讨论病态线性方程组的求解。
解答过程
1.估计矩阵的2-条件数和阶数的关系
矩阵的2-条件数定义为：1
222
()Cond A A A
-=⨯，将Hilbert 矩阵带入有：
1222
()Cond H H H -=⨯
调用自编的Hilbert_Cond 函数对其进行计算，取阶数n = 50，可得从1阶到50阶的2-条件数，以五位有效数字输出，其中前10项见表1。
表1.前十阶Hilbert 矩阵的2-条件数
从表1可以看出，随着阶数每递增1，Hilbert 矩阵的2-条件数都至少增加一个数量级，但难以观察出明显的相依规律。
故考虑将这些数据点绘制在以n 为横轴、Cond （H ）2为纵轴的对数坐标系中（编程用Hilbert_Cond 函数同时完成了这个功能），生成结果如图1。
图1.不同阶数下Hilbert矩阵的2-条件数分布
由图可见，当维数较小时，在y-对数坐标系中Cond（H）2与n有良好的线性关系；但n超过10后，线性趋势开始波动，n超过14后更是几乎一直趋于平稳。
事实上，从n = 12开始，系统便已经开始提出警告：“Warning: Matrix is close to singular or badly scaled.Results may be inaccurate.”。
也就是说，当n较大时，H矩阵已经接近奇异，计算结果可能是不准确的。
通过查阅相关资料，我找到了造成这种现象的原因：在matlab中，用inv函数求条件数过大的矩阵的逆矩阵将是不可靠的。
而调用系统自带的专门对Hilbert矩阵求逆的invhilb（n）函数则不存在这个问题，生成结果如图2。
图2. 修正后的不同阶数下Hilbert矩阵的2-条件数分布
简便起见，取n 不大于10的前十项进行线性拟合，结果如图3。
C o n d (H )
n
图3.1~10阶Hilbert 矩阵2-条件数的线性拟合
由拟合结果知，Cond （H ）2与n 的关系为：
2lg () 1.65183 1.47863Cond H n
=-+
其线性相关系数r=0.99985，可见二者具有较好的线性关系。
对上式稍作变形得：
1.651831.478632()10n Cond H -+=
这就是对Hilbert 矩阵的-2条件数与其阶数n 的关系估计。
可见Hilbert 矩阵的2-条件数会随其阶数n 的增加呈指数增大趋势，因此当n 较大时Hilbert 矩阵将是严重病态的，甚至导致matlab 中inv 求逆运算失真。
2.对不同的n ，采用各种方法求解方程
调用自编的Calc 主函数（其中包括的Hilbert 函数以及b_函数可创建出对应阶数的H 矩阵以及向量b ，Gauss_Cal 函数、Jacobi_Cal 函数、Gauss_Seidel_Cal 函数、SOR_Cal 函数（该函数自动寻找最优松弛因子，然后以最优因子进行求解）以及CG_Cal 函数则可完成各
自方法的求解），分别取n = 2,5,10,20,50，对于迭代法设定终止计算精度为10
10ε-=，所得
计算结果以16位有效数字输出，分别见表2~表6。
表2.n=2的计算结果
从表2可看出，n=2时，四种迭代法都能够收敛，迭代次数最大为e+2量级（J法），最小仅要2次（CG法），并且五种解法都能给出非常精确的结果，最大误差为e-10量级（GS 法）。
表3.n=5的计算结果
从表3可以看出，n=5时，J法已经不收敛，其余解法依然收敛（但值得注意的是GS 法以及SOR法的谱半径以及非常接近1，达到了收敛的边缘），最大迭代次数达到e+6量级（GS法），最小需要7次（CG法），计算精度依然较高，最大误差为e-6量级。
SOR法比
GS法节省了较多的迭代次数。
表4.n=10的计算结果
从表4可以看出，n=10时，各种解法的收敛性与n=5时相同（GS法与SOR法的谱半径进一步趋近1），最大迭代次数达e+8量级（SOR法），最小需要8次（CG法），计算精度已经较低，最大误差达到e-3量级。
此时有一个异常现象：SOR法的迭代次数不但不比GS法少，反而超出好几倍。
但根据计算出的两种算法下的迭代矩阵谱半径，可以看出SOR 法为0.999999999871931，而GS法为0.999999999997045，在最优松弛因子之下的SOR法确实具有更小的迭代矩阵谱半径，因此应当更快收敛。
经检查，计算程序的编制应该没有错误，但计算结果确实与理论不符，希望老师指点迷津！
表5.n=20的计算结果
从表5可以看出，n=20时，各种解法的收敛性依然没变（但GS法以及SOR法的谱半径在matlab的十六位显示精度下已经等于1），最大迭代次数e+8量级（SOR法），最小需要10次（CG法），Gauss消元法的计算结果已失去意义（误差e+2量级），迭代法的计算误差均为e-3量级。
且此时SOR的迭代次数依然高于GS。
表6.n=50的计算结果
从表6可以看出，n=50时，计算结果的特点与n=20相似，不同之处在于Gauss消元法的误差进一步增大（e+3量级）。
3.综合讨论
根据上一部分的计算结果，可以总结出求解病态矩阵的两个特点：
1.收敛性差，收敛速度慢。
对于本例，当阶数n大于2以后J法就不收敛了，这是收敛性差的一个体现。
GS法和SOR法虽然一直保持收敛（事实上这是由于Hilbert矩阵是正定对称的，所以决定了GS法以及SOR法一定收敛），但随着n的增大它们的谱半径迅速趋近于1（n=20时matlab的16位显示精度已经无法显示出它们的谱半径与1的差别），根据理论我们知道，趋近于1的谱半径决定了计算的收敛速度将会非常慢，之前计算结果中的迭代次数N也很好地验证了这一点。
2.计算精度低，阶数较高时误差惊人。
根据理论，线性方程组系数矩阵的条件数的直观意义就是初始数据扰动的相对误差的放大倍数。
由于数值计算时的舍入误差等因素，初始数据的微小扰动难以避免，因此解病态方程组时，这种微小的初始相对误差将被病态矩阵巨大的条件数所放大，最后造成计算结果的巨大偏差。
对于本例，n=20时Gauss消元法的计算误差竟然达到e+2量级，使得计算结果完全失去意义，而其他几种迭代法的计算误差看似能够接受（e-3量级），但仔细分析上一节
的计算结果不难发现：从n=2到n=5再到n=10，几种迭代法的迭代次数迅速增大，且计算精度显著降低；而从n=10到n=20再到n=50，迭代次数和计算精度却趋于稳定，由于矩阵的条件数一直呈指数增大，因此这种计算结果的趋势明显是不合理的。
究其原因，我认为和一开始图1所示的2-条件数的计算失真是类似的。
当矩阵的条件数非常大时，inv指令的求逆结果不可靠，从而造成图1中2-条件数趋于平稳的那个阶段，而在GS迭代以及SOR迭代中也需要对矩阵进行求逆运算，因此很有可能又是因为inv对高条件数矩阵求逆的偏差造成最终计算结果的失真。
因此，我预计实际上当n=20时，GS法和SOR法的迭代次数以及计算误差要远远高于n=10的值，同理n=50时的结果又要远远大于n=20的值。
不难看出，造成上述两个特点的本质因素还是因为病态矩阵的高条件数，无论收敛性、收敛速度还是计算精度等，都与矩阵的条件数直接相关，因此改善病态矩阵的求解的根本方法还是在于降低矩阵的条件数。
因此我尝试着对Hilbert矩阵进行了一些预处理，并发现在某些预处理方法下条件数的阶数被有效减小，因此预处理是求解病态方程组的一个十分有效且必要的环节。
此外，若给定条件数并横向比较上一节中的不同求解方法，可以看出CG法的表现是最令人满意的，当矩阵阶数n不太大时（不超过10），CG法都能以非常少的迭代次数（只需要几次）以及较高的计算精度完成计算任务，相比之下，GS法以及SOR法的上千万次迭代次数显得非常吃力。
此外，事实上n较小时Gauss消元法的计算结果也比较理想，具有很快的计算速度和较高的计算精度。
我们知道，在传统意义上Gauss消元法的最大弊端在于n非常大时计算量大得令人难以接受（n三次方量级），而迭代法的单步计算量则为n 平方量级，但在求解病态问题时，n只能比较小，且迭代法的迭代次数却非常大，所以造成迭代法的计算速度反而远慢于Gauss消元法。
综上，求解病态方程组时首先应该对系数矩阵进行预处理，尽可能降低它的条件数；其次可以选择共轭梯度（CG）法对其进行求解。
如果阶数n不大，Gauss消元法也是比较有效的。
