数学实验 5：线性代数方程组的数值解法

2、 运行结果分析 1） 不同初值（x0）和右端项（b）条件下的解的情况
表一：收敛性判断
1-范数 Jacobi 0.0163 0.0167 2-范数 0.0163 -范数 0.0167
Gauss 0.0008 0.0084 0.0084 0.0084 可以看到，矩阵A无论是谱半径或是任意范数的值都小于1，可知在 A不变的情况下，Jacobi和Gauss法必然收敛。 2） b取不同的值，x0=20*ones(20,1), e=10^-5, m=50 条 件下的情况对比 迭代 B= B= B= B=20*ones(20,1) B=2000*ones(20,1) 次数 [1:20]’ [10:10:200]’ [20:-1:1]’ Jacobi 24 26 24 23 30 Gauss 16 17 16 15 20 根据1）分析的结果，可以证明无论b取任何值，采用两种方法迭代 均收敛，但b的值的变化会影响迭代的次数；且Gauss迭代法总是比 Jacobi迭代法收敛速度更快。 下表列出的是B=[1:20]’情况下部分结算结果，可以很明显的看到两 种迭代法的收敛速度不同: Jacobi K=1 K=2 K=3 K=4 K=5 K=6 … K=22 K=23 X1 X2 X3 5.3333 9.0000 2.7500 1.5394 1.0874 0.8858 0.7982 … 0.7247 0.7247 4.3333 2.6644 1.9248 1.6082 1.4652 … 1.3444 1.3444 K=2 K=3 K=4 K=5 K=6 … K=14 K=15
由于精度问题，在迭代的最后几次中从显示的数位已经不能看出 标准值与计算值得差别，但是若采用long显示设定，就可以看到更多 位小数的显示，其结果符合最初设定的精度e，数据繁琐，略。 3） x0取不同的值，b=20*ones(20,1), e=10^-5, m=50 条 件下的情况对比 迭代 X0= X0= X0= X0=20*ones(20,1) X0=2000*ones(20,1) 次数 [1:20]’ [10:10:200]’ [20:-1:1]’ Jacobi 22 27 22 23 31 Gauss 15 18 15 15 20 根据1）分析的结果，可以证明无论X0取任何值，采用两种方法迭代 均收敛，但x0的值的变化会影响迭代的次数；且Gauss迭代法总是比 Jacobi迭代法收敛速度更快。 下表列出的是x0=[1:20]’情况下部分结算结果，可以很明显的看到两 种迭代法的收敛速 度不同: Jacobi K=1 K=2 K=3 K=4 K=5 K=6 … K=20 K=21 X1 X2 X3 X1 X2 X3 7.2500 8.6250 9.2228 9.4536 9.5541 9.5974 … 9.6327 9.6327 7.6667 9.9583 10.814 11.183 11.341 11.411 … 11.4683 11.4683 8.1667 10.757 11.816 12.282 12.487 12.579 … 12.6560 12.6560 K=Байду номын сангаас K=3 K=4 K=5 K=6 … K=13 K=14 9.6327 7.2500 8.9352 9.4267 9.5720 9.6150 9.6276 … 9.6327
break end end e=eig(bgs) n1=norm(bgs,1); n2=norm(bgs); nn=norm(bgs,inf); min([n1 n2 nn]) 3） 操作函数： %构造矩阵A n=20; a1=sparse(1:n,1:n,3,n,n); 造，比较方便 a2=sparse(1:n-1,2:n,-0.5,n,n); a3=sparse(1:n-2,3:n,-0.25,n,n); a=a1+a2+a3+a2'+a3'; a=full(a);
11.0000 5.8056 3.7199 2.7801 2.3626 2.1717 … 2.0072 2.0072
Gauss K=1 X1 X2 X3 5.3333 6.5556 7.5370
2.0540 1.1354 0.8539 0.7657 0.7377 … 0.7247 0.7247 2.9432 1.8459 1.5033 1.3949 1.3605 … 1.3444 1.3444 3.7386 2.5553 2.1815 2.0627 2.0249 … 2.0072 2.0072
%按稀疏矩阵的输入法构
%还原为满矩阵
%通过给定不同的初始向量x0或者右端项b，以及规定不同的误差要 求，进行jacobi和gauss %迭代，得到的结果y1、y2位两种迭代的次数，同时输出迭代结果，便 于分析 b= x0= e= m= y1=jacobi(a,b,x0,e,m); y2=gauss(a,b,x0,e,m); 4） 改变矩阵A： 先对jacobi函数作一定修正，方便分析，命名为jacobi2,如下: function y=jacobi2(a,b,x0,e,m) d=diag(diag(a)); u=-triu(a,1); l=-tril(a,-1);
结果分析： 根据判断迭代是否收敛的原理，若A是严格对角占优的，即，则 Jacobi迭代和Gauss迭代均收敛。由此推测，A的对角占优程度可能是影 响收敛速度的关键因素。由上述试验表面，对A的对角元素数值的扩大 的确可以明显的改善Jacobi迭代的收敛速度，与预测的情况相符。 主对角占优程度可以影响收敛速度，但是这个结论因为数学水平有
小值q 0.4893 0.2447 0.1631 0.1223 0.0979 0.0816 0.0699 0.0612 0.0544 0.0489 0.0445 0.0408 0.0376 0.0350 0.0326 0.0306 0.0288 0.0272 0.0258 0.0245 21.0000 12.0000 9.0000 8.0000 8.0000 7.0000 7.0000 7.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 5.0000 5.0000
实验 5：线性代数方程组的数值解法
习题3：
已知方程组,其中，定义为： 试通过迭代法求解此方程组，认识迭代法收敛的含义以及迭代初值和方 程组系数矩阵性质对收敛速度的影响。实验要求： （1） 选取不同的初始向量x0和不同的方程组右端向量b，给定迭 代误差要求，用雅可比迭代法和高斯-赛德尔迭代法计算， 观测得到的迭代向量序列是否均收敛？若收敛，记录迭代 次数，分析计算结果并得出结论； （2） 取定右端向量b和初始向量x0，将A的主对角线元素成倍的 增长若干次，非主对角元素不变，每次用雅可比迭代法计 算，要求迭代误差满足,比较收敛速度，分析现象并得出结 论。 1、 程序设计（可直接粘贴运行） 1） Jacobi迭代法 function y=jacobi(a,b,x0,e,m) %定义jacobi函数，其中：a,b为线性方程组中的矩阵和右端向量；x0 为初始值； %e和m分别为人为设定的精度和预计迭代次数；运行结果y为迭代的结 果和所有中间值组成的 %矩阵 y=0; %对y初始化 d=diag(diag(a)); %按雅可比迭代标准形 形式取主对角元素作为矩阵D u=-triu(a,1); %取上三角矩阵u l=-tril(a,-1); %取下三角矩阵l bj=d^-1*(l+u); fj=d^-1*b; x=[x0,zeros(20,m-1)]; %初始化x,其中x1=x0,即 初始值 for k=1:m %人为规定迭代次 数，防止不收敛迭代导致死循环 x(:,k+1)=bj*x(:,k)+fj; %jacobi迭代 if norm(x(:,k+1)-x(:,k),inf)<e
%判断迭代后是否满足迭代中止条件： y=x(:,1:k+1); 和迭代结果 sizej=k ; 输出迭代次数 break end 迭代 end
%赋给y所有中间值 %若去掉;号，则 %并结束迭代 %若不成立，继续
%以下部分为验证迭代公式收敛的方法，仅需运行一次即可，因为收敛 性完全由A矩阵决定,而A %在本题是固定不变的；通过判断中B的谱半径或范数大小（B在jacobi 迭代法中%为矩阵bj），即可得知收敛性: e=eig(bj) %输出全部特征值，若，则收敛 n1=norm(bj,1); %计算1-范数 n2=norm(bj); %计算2-范数 nn=norm(bj,inf)； %计算-范数 q=min([n1 n2 nn]) %由于谱半径不超过人以一种范数，所以只要范数的最小值q<1,也可判 断迭代法收敛 2） Gauss迭代法：与Jacobi程序结构相同，不再注释 function y=gauss(a,b,x0,e,m) y=0; d=diag(diag(a)); u=-triu(a,1); l=-tril(a,-1); x=[x0,zeros(20,m-1)]; bgs=(d-l)^-1*u; fgs=(d-l)^-1*b; for k=1:m x(:,k+1)=bgs*x(:,k)+fgs; if norm(x(:,k+1)-x(:,k),inf)<e y=x(:,1:k+1); sizeg=k;
限没有得到严格的推导。但关于A的范数和收敛速度的关系在书上有比 较完整的证明，即有公式： （1） q为A矩阵的三种范数中最小的值，x*为原线性方程组的解，可知q越 小，序列收敛越快。从上表中可以看出，迭代次数的减小(收敛速度的 增加)，是和q的减小同步发生的，公式（1）得到验证。 还可以看到，随着扩大倍数的增大，迭代速度的增长率在不断的下 降。可以直观的理解为,主对角元素的相对优势在扩大一定倍数之后已 经相当明显，即使再增大若干倍其对于非对角元素仍然是具有压倒性优 势的，优势地位不会因此而改变很多，根据前面的假设，收敛速度与对 角占优程度有关，速度的改变也将减缓。
bj=d^-1*(l+u); fj=d^-1*b; x=[x0,zeros(20,m-1)]; n1=norm(bj,1); %计算范数 n2=norm(bj); nn=norm(bj,inf); q=min([n1 n2 nn]); y(1)=q; %输出结果1：范数的最 小值，判断收敛速度的方法 for k=1:m x(:,k+1)=bj*x(:,k)+fj; if norm(x(:,k+1)-x(:,k),inf)<e y(2)=k; %输出结果2：迭 代次数 break end end %改变A矩阵主对角元素的值，比较jacobi迭代的收敛速度，即迭代误 差满足%时的迭代次数 b=(1:20)'; %取定b x0=20*ones(20,1); %取定x0 e=10^-5; %取定精度 r=20; %设置主对角元素 扩大的最大倍数 y=0; m=50; for k=1:r; a1=k*sparse(1:n,1:n,3,n,n); %只需更改A的主对角元 素 a=a1+a2+a3+a2'+a3'; %重新构造A a=full(a); y(k,1:3)=[k,jacobi2(a,b,x0,e,m)]; end y %输出对角线扩大倍 数\最小范数\迭代次数
Gauss K=1
8.7083 10.654 11.228 11.398 11.448 11.463 … 11.4683 11.4683 9.8056 11.813 12.408 12.584 12.636 12.650 … 12.6560 12.6560 4） b=(1:20)';x0=20*ones(20,1);e=10^-5; m=50; （固定各个参数） r=20;（改变a的主对角元素，依次扩大1倍、2倍……20倍） 运行结果： 主对角元素扩 三个范数的最 迭代次数
大倍数 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.0000 11.0000 12.0000 13.0000 14.0000 15.0000 16.0000 17.0000 18.0000 19.0000 20.0000