1、姚文 200631000091 06 信计第 1 页 共 12 页数学实验报告一、实验目的（1）复习用 Gauss_Seidel，Jacobi ，超松弛迭代法求解线性方程组。（2）比较用以上各种方法求解问题所要求满足的条件二、实验问题构造 3 阶 Hilbert 矩阵 A，并使其精确解 x 为 3 维全一向量（1 1 1） ，进而构造出右端向量 b。对于方程组 Ax=b，从理论上分析 Jacobi 迭代法和 Gauss_Seidel 迭代法是否收敛，并求解各自的谱半径和收敛速度。逐步增加 n 的值，重复求解上述问题，并总结计算结果。三、问题的分析对于三阶 Hilbert 矩阵，以下分析其用 G
2、auss_Seidel 迭代法和用 Jacobi 迭代法所得的迭代矩阵的谱半径，并由此分析其用相应方法求解的可行性。（1）在 MATLAB 的 command window 中输入以下命令： format rat %有理数格式 a=hilb(3)a =1 1/2 1/3 1/2 1/3 1/4 1/3 1/4 1/5 D=diag(diag(a),0)D =1 0 0 0 1/3 0 0 0 1/5 %取矩阵 a 的对角元素 L=-tril(a,-1)L =0 0 0 -1/2 0 0 -1/3 -1/4 0 姚文 200631000091 06 信计第 2 页 共 12 页 %取矩阵 a 的
3、下三角矩阵的负值 U=-triu(a,1)U =0 -1/2 -1/3 0 0 -1/4 0 0 0 %取矩阵 a 的上三角矩阵的负值 B=inv(D-L)*UB =0 -1/2 -1/3 0 3/4 -1/4 0 -5/48 125/144 %矩阵 B 为矩阵 a 用 Gauss_seidel 法求解时得到的的迭代矩阵 pr=max(abs(eig(B)pr =2101/2142 %迭代矩阵 B 的谱半径由上述结果可见，用 Gauss_Seidel 法得到的迭代矩阵谱半径小于一，用此方法求解，结果收敛。（2）在 MATLAB 的 command window 中输入以下命令： format
4、rat a=hilb(3); D=diag(diag(a),0)D =1 0 0 0 1/3 0 0 0 1/5 %取矩阵 a 的对角元素 B=inv(D)*(D-a)B =0 -1/2 -1/3 姚文 200631000091 06 信计第 3 页 共 12 页-3/2 0 -3/4 -5/3 -5/4 0 %Jacobi 方法的迭代矩阵 pr=max(abs(eig(B)pr =1449/841 %迭代矩阵的谱半径由上述结果可见，用 Jacobi 法得到的迭代矩阵谱半径大于一，用此方法求解，结果将不收敛。下面编写程序，对于输入的正整数 w，构造维数为 w 的 Hilbert 矩阵，并以全一
5、向量为精确解构造右端向量，然后分别用 Gauss_Seidel 方法和 Jacobi 方法求解方程组，并与精确解比较误差。（a）用 Gauss_Seidel 迭代法分析该问题，用 MATLAB 编辑程序如下：%Gauss_Seidel 迭代法的程序functionres,k,pr,wc=Gauss_Seidel(err,M)%输出结果为方程组的解 res,迭代次数 k 和系数矩阵的谱半径 pr%err 为误差控制，M 为限制的最大迭代次数format longw=input(please input a number:); %输入一个正整数 w，为 hilbert 矩阵的维数A=hilb(w)
6、; %根据输入的正整数构造 hilbert 矩阵 Am=size(A); %求 A 的维数n=m(1);xjq=ones(n,1); %方程组的精确解x0=zeros(n,1);b=A*xjq; %方程组的右端向量D=diag(diag(A),0);L=-tril(A,-1);U=-triu(A,1);B=inv(D-L)*U;g=inv(D-L)*b;pr=max(abs(eig(B); %求迭代矩阵 B 的谱半径x1=B*x0+g;k=0;姚文 200631000091 06 信计第 4 页 共 12 页while(norm(x1-x0),2)err)k=k+1;x0=x1;x1=B*x0
7、+g;if(kM)disp(迭代次数过大！),return;endend k;res=x1;wc=norm(res-xjq,2); %wc 为所求结果和精确解之差的二范数（b）用 jacobi 迭代法分析该问题，用 MATLAB 编辑程序如下：%jacobi 迭代法的程序functionres,k,pr,wc=Jacobi(err,M)%输入：err 为误差控制，M 为最大迭代次数%输出：方程组的解 res，迭代次数 k 和系数矩阵的谱半径 pr，精确解与所得解的误差 wcformat longw=input(please input a number:); %输入一个正整数 w，为 hilb
8、ert 矩阵的维数A=hilb(w); %根据输入的正整数构造 hilbert 矩阵 Am=size(A); %求 A 的维数n=m(1);xjq=ones(n,1); %方程组的精确解x0=zeros(n,1);b=A*xjq; %方程组的右端向量D=diag(diag(A),0); %D 为以 A 的对角元素为其对角元素的对角阵B=inv(D)*(D-A); %B 为迭代矩阵g=inv(D)*b; %g 为右端向量pr=max(abs(eig(B); %求迭代矩阵 B 的谱半径x1=B*x0+g; %用 x0 和 x1 存放相邻两次迭代结果k=0; %用 k 记录迭代次数while(nor
9、m(x1-x0),2)err) %用二范数循环控制k=k+1;x0=x1;x1=B*x0+g; if(kM)disp(迭代次数过大！),return;end姚文 200631000091 06 信计第 5 页 共 12 页endk;res=x1;wc=norm(res-xjq,2); %wc 为所求结果和精确解之差的二范数下面程序使用超松弛迭代法，求解以上问题：%超松弛方法function res,wc=SOR(err,w)%输入：误差控制 wc，松弛因子 w%输出：方程组的解 res，精确解与所得解之间的误差%w 为松弛因子，w=1 为 Gauss-Seidel 迭代法,w1 称为超松弛法f
10、ormat longl=input(please input a number:); %输入一个正整数 l，为 hilbert 矩阵的维数A=hilb(l); %根据输入的正整数构造 hilbert 矩阵 Am=size(A); %求 A 的维数n=m(1);xjq=ones(n,1); %方程组的精确解x0=zeros(n,1);b=A*xjq; %方程组的右端向量n=length(b);x=zeros(n,1);for k=1:1000xk=x;error=0;for i=1:nsum=0;for j=1:nsum=sum+A(i,j)*x(j);endx(i)=xk(i)+w*(b(i)
11、-sum)/A(i,i);error=error+abs(x(i)-xk(i);endif error/norm(x,1) err=0.000001; M=100000; res,k,pr,wc=Gauss_Seidel(err,M)得到结果：please input a number:输入 3 之后得到：res =0.999992820695231.000036640590620.99996616476968k =480pr =0.98085893099526wc =5.038747968330873e-005(2)若输入命令： res,k,pr,wc=Jacobi(err,M)得：plea
12、se input a number:3res =InfInfNaNk =1304姚文 200631000091 06 信计第 7 页 共 12 页pr =1.72294966962992wc =NaN（3）下面分别取 n=6，8，10，得到的结果如下：（3.1）n=6： res,k,pr,wc=Gauss_Seidel(err,M)please input a number:6res =0.999930634133761.000919479110870.998101636818390.997378077097901.008901977986930.99470492103412k =17405p
13、r =0.99999829925262wc =0.01089089478727 res,k,pr,wc=Jacobi(err,M)please input a number:6姚文 200631000091 06 信计第 8 页 共 12 页res =InfInfNaNNaNNaNNaNk =486pr =4.30853103479331wc =NaN（3.2）n=8： res,k,pr,wc=Gauss_Seidel(err,M)please input a number:8res =1.000101061834730.997417882666001.013609652247160.9794
14、44297292960.998171505508771.014086036515591.010086216191120.98695590212338k =8341姚文 200631000091 06 信计第 9 页 共 12 页pr =0.99999999762484wc =0.03298601453316 res,k,pr,wc=Jacobi(err,M)please input a number:8res =-InfNaNNaNNaNNaNNaNNaNNaNk =395pr =6.04213434207628wc =NaN（3.3）n=10： res,k,pr,wc=Gauss_Seide
15、l(err,M)please input a number:10res =1.00010998927114姚文 200631000091 06 信计第 10 页 共 12 页0.998182398194581.005617496487180.999325180802160.991557970661490.996796073759101.005244507984471.008734188342541.003890304015690.99044174041455k =26950pr =0.99999999999705wc =0.01808718126312 res,k,pr,wc=Jacobi(e
16、rr,M)please input a number:10res =InfInfInfNaNNaNNaNNaNNaNNaNNaNk =姚文 200631000091 06 信计第 11 页 共 12 页346pr =7.77981513192998wc =NaN（4）下面使用超松弛方法求解上述问题：在 MATLAB 的 command window 中输入以下命令： err=0.000001; w=1; res,wc=SOR(err,w)please input a number:3res =0.998595589658141.007167605506770.99338117701964wc
17、=0.00985676187976事实上，当 w=1 时，超松弛方法和 Gauss_Seidel 得到的结果是一样的。通过选取不通的 w，并比较 wc 的大小，可得当 w=1.736 时，sor 方法所得结果最精确。 w=1.736; res,wc=SOR(err,w)please input a number:3res =0.999926971809941.000001999288241.00015674197053wc =1.729310817139307e-004姚文 200631000091 06 信计第 12 页 共 12 页五、问题求解结果的分析和结论由（四）中的计算结果可见：对于
18、同一个问题，利用不同的方法，会得到完全不同的结果。对于本实验中的 Hilbert 矩阵，当 n=3 时，用 Gauss_Seidel 法得到的迭代矩阵的谱范数为 0.980858930995261，因而此方法不收敛；精确解和所得解之间的误差wc=NaN（无穷大） 。并且可以看出，随着 Hilbert 矩阵的维数 n 的增大，上述两种方法所得的迭代矩阵的谱半径也随之增大，所得解和精确解之间的误差也越来越大。六、实验的总结和体会Gauss_Seidel 方法和 Jacobi 方法求解线性方程组，要求所得迭代矩阵的谱半径小于一，否则所得结果将不会收敛。SOR 方法一般比 Gauss_Seidel 方法和 Jacobi 方法更精确，但是松弛因子事先不知，需要自己确定，一般先取 w=1，然后在其左右变动，最终确定一个 w 值，使得结果比较精确。通过本实验，是我们更“亲身”体会到了各种不同的求解线性方程组的方法的要求和其优劣。
