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

1、实验 5：线性代数方程组的数值解法习题3：已知方程组,其中，定义为： 试通过迭代法求解此方程组，认识迭代法收敛的含义以及迭代初值和方程组系数矩阵性质对收敛速度的影响。实验要求：（1） 选取不同的初始向量x0和不同的方程组右端向量b，给定迭代误差要求，用雅可比迭代法和高斯-赛德尔迭代法计算，观测得到的迭代向量序列是否均收敛？若收敛，记录迭代次数，分析计算结果并得出结论；（2） 取定右端向量b和初始向量x0，将A的主对角线元素成倍的增长若干次，非主对角元素不变，每次用雅可比迭代法计算，要求迭代误差满足,比较收敛速度，分析现象并得出结论。1、 程序设计（可直接粘贴运行）1） Jacobi迭代法fun

2、ction y=jacobi(a,b,x0,e,m) %定义jacobi函数，其中：a,b为线性方程组中的矩阵和右端向量；x0为初始值；%e和m分别为人为设定的精度和预计迭代次数；运行结果y为迭代的结果和所有中间值组成的%矩阵y=0;%对y初始化d=diag(diag(a);%按雅可比迭代标准形形式取主对角元素作为矩阵Du=-triu(a,1);%取上三角矩阵ul=-tril(a,-1);%取下三角矩阵lbj=d-1*(l+u);fj=d-1*b;x=x0,zeros(20,m-1);%初始化x,其中x1=x0,即初始值 for k=1:m%人为规定迭代次数，防止不收敛迭代导致死循环x(:,k

3、+1)=bj*x(:,k)+fj;%jacobi迭代 if norm(x(:,k+1)-x(:,k),inf)<e %判断迭代后是否满足迭代中止条件： y=x(:,1:k+1);%赋给y所有中间值和迭代结果 sizej=k;%若去掉;号，则输出迭代次数 break %并结束迭代 end%若不成立，继续迭代end %以下部分为验证迭代公式收敛的方法，仅需运行一次即可，因为收敛性完全由A矩阵决定,而A%在本题是固定不变的；通过判断中B的谱半径或范数大小（B在jacobi迭代法中%为矩阵bj），即可得知收敛性:e=eig(bj)%输出全部特征值，若，则收敛n1=norm(bj,1);%计算1-

4、范数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:mx(:,k+1)=bgs*x(:,k)+fgs; if norm(x(:,k+

5、1)-x(:,k),inf)<ey=x(:,1:k+1); sizeg=k; break endende=eig(bgs)n1=norm(bgs,1);n2=norm(bgs);nn=norm(bgs,inf);min(n1 n2 nn) 3） 操作函数：%构造矩阵An=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);%还原为满矩阵%通过给定不同的初始向量

6、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);bj=d-1*(l+u);fj=d-1*b;x=x0,zeros(20,m-1); n1=norm(bj,1);%计算范

7、数n2=norm(bj);nn=norm(bj,inf);q=min(n1 n2 nn);y(1)=q;%输出结果1：范数的最小值，判断收敛速度的方法 for k=1:mx(:,k+1)=bj*x(:,k)+fj; if norm(x(:,k+1)-x(:,k),inf)<e y(2)=k;%输出结果2：迭代次数 break endend %改变A矩阵主对角元素的值，比较jacobi迭代的收敛速度，即迭代误差满足%时的迭代次数b=(1:20)'%取定bx0=20*ones(20,1);%取定x0e=10-5;%取定精度r=20; %设置主对角元素扩大的最大倍数 y=0;m=50;

8、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);endy%输出对角线扩大倍数最小范数迭代次数2、 运行结果分析1） 不同初值（x0）和右端项（b）条件下的解的情况表一：收敛性判断1-范数2-范数-范数Jacobi0.01630.01670.01630.0167Gauss0.00080.00840.00840.0084可以看到，矩阵A无论是谱半径或是任意范数的值都小于1，可知在A不变的情况下，Ja

9、cobi和Gauss法必然收敛。2） b取不同的值，x0=20*ones(20,1), e=10-5, m=50 条件下的情况对比迭代次数B=1:20B=10:10:200B=20:-1:1B=20*ones(20,1)B=2000*ones(20,1)Jacobi2426242330Gauss1617161520根据1）分析的结果，可以证明无论b取任何值，采用两种方法迭代均收敛，但b的值的变化会影响迭代的次数；且Gauss迭代法总是比Jacobi迭代法收敛速度更快。下表列出的是B=1:20情况下部分结算结果，可以很明显的看到两种迭代法的收敛速度不同:JacobiK=1K=2K=3K=4K=5

10、K=6 K=22K=23K=24标准值X15.33332.75001.53941.08740.88580.79820.72470.72470.72470.7247X29.00004.33332.66441.92481.60821.46521.34441.34441.34441.3444X311.00005.80563.71992.78012.36262.17172.00722.00722.00722.0072GaussK=1K=2K=3K=4K=5K=6K=14K=15K=16标准值X15.33332.05401.13540.85390.76570.73770.72470.72470.7247

11、0.7247X26.55562.94321.84591.50331.39491.36051.34441.34441.34441.3444X37.53703.73862.55532.18152.06272.02492.00722.00722.00722.0072由于精度问题，在迭代的最后几次中从显示的数位已经不能看出标准值与计算值得差别，但是若采用long显示设定，就可以看到更多位小数的显示，其结果符合最初设定的精度e，数据繁琐，略。3） x0取不同的值，b=20*ones(20,1), e=10-5, m=50 条件下的情况对比迭代次数X0=1:20X0=10:10:200X0=20:-1:1

12、X0=20*ones(20,1)X0=2000*ones(20,1)Jacobi2227222331Gauss1518151520根据1）分析的结果，可以证明无论X0取任何值，采用两种方法迭代均收敛，但x0的值的变化会影响迭代的次数；且Gauss迭代法总是比Jacobi迭代法收敛速度更快。下表列出的是x0=1:20情况下部分结算结果，可以很明显的看到两种迭代法的收敛速度不同:JacobiK=1K=2K=3K=4K=5K=6 K=20K=21K=22标准值X17.25008.62509.22289.45369.55419.59749.63279.63279.63279.6327X27.66679

13、.958310.81411.18311.34111.41111.468311.468311.468311.4683X38.166710.75711.81612.28212.48712.57912.656012.656012.656012.6560GaussK=1K=2K=3K=4K=5K=6K=13K=14K=15标准值X17.25008.93529.42679.57209.61509.62769.63279.63279.63279.6327X28.708310.65411.22811.39811.44811.46311.468311.468311.468311.4683X39.805611.

14、81312.40812.58412.63612.65012.656012.656012.656012.65604） b=(1:20)'x0=20*ones(20,1);e=10-5; m=50;（固定各个参数）r=20;（改变a的主对角元素，依次扩大1倍、2倍20倍）运行结果：主对角元素扩大倍数三个范数的最小值q迭代次数1.00000.489321.00002.00000.244712.00003.00000.16319.00004.00000.12238.00005.00000.09798.00006.00000.08167.00007.00000.06997.00008.00000

15、.06127.00009.00000.05446.000010.00000.04896.000011.00000.04456.000012.00000.04086.000013.00000.03766.000014.00000.03506.000015.00000.03266.000016.00000.03066.000017.00000.02886.000018.00000.02726.000019.00000.02585.000020.00000.02455.0000结果分析：根据判断迭代是否收敛的原理，若A是严格对角占优的，即，则Jacobi迭代和Gauss迭代均收敛。由此推测，A的对角

16、占优程度可能是影响收敛速度的关键因素。由上述试验表面，对A的对角元素数值的扩大的确可以明显的改善Jacobi迭代的收敛速度，与预测的情况相符。主对角占优程度可以影响收敛速度，但是这个结论因为数学水平有限没有得到严格的推导。但关于A的范数和收敛速度的关系在书上有比较完整的证明，即有公式：（1）q为A矩阵的三种范数中最小的值，x*为原线性方程组的解，可知q越小，序列收敛越快。从上表中可以看出，迭代次数的减小(收敛速度的增加)，是和q的减小同步发生的，公式（1）得到验证。 还可以看到，随着扩大倍数的增大，迭代速度的增长率在不断的下降。可以直观的理解为,主对角元素的相对优势在扩大一定倍数之后已经相当明

17、显，即使再增大若干倍其对于非对角元素仍然是具有压倒性优势的，优势地位不会因此而改变很多，根据前面的假设，收敛速度与对角占优程度有关，速度的改变也将减缓。习题8 种群的繁殖与稳定收获：种群的数量因繁殖而增加，因自然死亡而减少，对于人工饲养的种群（比如家畜）而言，为了保证稳定的收获，各个年龄的种群数量应维持不变，种群因雌性个体的繁殖而改变，为方便起见以下种群数量均指其中的雌性。种群年龄记作k=1，2，n，当年年龄k的种群数量记作xk，繁殖率记作bk（每个雌性个体1年繁殖的数量），自然存活率记作sk（sk=1-dk，dk为1年的死亡率），收获量记作hk，则来年年龄k的种群数量应为，要求各个年龄的种群

18、数量每年维持不变就是要使。(1)若已知bk，sk，给定收获量hk，建立求各个年龄的稳定种群数量xk的模型（用矩阵、向量表示）。(2)设n=5，b1=b2=b5=0，b3=5，b4=3，s1=s4=0.4，s2=s3=0.6，如果要求h1h5为500，400，200，200，100，求(3)要使h1h5均为500， 如何达到？1.模型的建立根据题目给出的模型和各个参量，建立线性方程组。（1）参数b,s,h的意义如题目所述，向量，表示年龄为k的种群在第i个统计时间段中的数量。题目中还给出了使各年龄的种群数量每年维持不变的要求，即在k取一定值时，不随i的变化而变化。方程（1）变为如下形式：（2）得到

19、了方程（2）的形式，就可以通过直接或间接（迭代）方法求解线性方程组。2.程序设计1）构造矩阵n=5;b=zeros(1,n);b(3)=5;b(4)=3;s=0.4 0.6 0.6 0.4;h=0 500 400 200 100'ss=diag(s);m=b;ss,zeros(n-1,1)-eye(5);x1=100*ones(5,1); %初值x1，计算guass迭代使用到x=gauss1(m,h,x1); %试通过gauss迭代法得到方程组的解xx=mh; %xx为matlab通过内定方法得到的解，可以认为是精确的 2）高斯函数: 试采用高斯迭代法求解方程，内部设定精度和迭代次数，

20、输出函迭代的结果和中间值，并输出全部特征值和范数的最小值function x=gauss1(a,b,x0) d=diag(diag(a);u=-triu(a,1);l=-tril(a,-1);bgs=(d-l)-1*u;fgs=(d-l)-1*b;e=10.-3; %设定精度m=20; %迭代次数x=x0;for k=1:mx(:,k+1)=bgs*x(:,k)+fgs;if norm(x(:,k+1)-x(:,k),inf)<e;x(:,1:k+1);breakendend e=eig(bgs) %输出全部特征值n1=norm(bgs,1); n2=norm(bgs);nn=norm(

21、bgs,inf);q=min(n1 n2 nn) %输出范数的最小值3、 运行结果分析1） 问题1的方程已在“模型的建立”中给出2） 将条件n=5，b1=b2=b5=0，b3=5，b4=3，s1=s4=0.4，s2=s3=0.6，h=500,400,200 ,100,100带入方程（2）,用左除程序和gauss迭代分别求解x1x5运行结果：e = 0 0 0 0 1.6320q = 6.4974x=1.0e+007 *colume 0 1 234 16 17 18 19 20xx = 1.0e+003 * 8.4810 2.8924 1.3354 0.6013 0.1405由特征值的和范数大于

22、1的结果可以分析出，采用gauss迭代法不能收敛，具体运算结果x发散现象很明显，迭代四次就会出现1e5数量级的值，验证了上面的推断。方程的正确结果为xx，即稳定状态下，种群处于1、2、3、4、5年龄段的个体数目分别为：8481、2892、1335、601、141，数量维持不变。3） 改变h1h5均为500，再次求解x,并讨论结果由于矩阵m（2）方程中左端矩阵，没有改变，可知方程仍然不可以用gauss或者jacobi迭代法进行求解。直接由左除运算得到：xx = 1.0e+004 * 1.1772 0.4209 0.2025 0.0715 -0.0214保持自然存活率和繁殖率不变，在对每个年龄段的

23、人工收获量都为500的情况下，达到稳定时，种群处于1、2、3、4、5年龄段的个体数目分别为：11772、4209、2025、715、-214，数量维持不变。可以从数据中看到，年龄段5的个体数目出现了负值，这在现实中是不可能出现的。分析：问题1：由方程（2）的建立过程可以看出，影响种群某年龄段的主要因素是上一年龄段的存活率和人工收获量，两者对于该年龄段种群数量产生相反的作用，而特殊的，对于年龄段1，即种群中刚刚出的个体数量，其仅仅受到各个种群的繁殖率的影响。出现解为负值的情况在数学上是由方程本身的性质决定的，是正常现象，但联系本题涉及的实际情况却反映出了较大的问题：它说明此种人工繁殖和收获的方法

24、是不能实现的。问题2：不收敛性，由于方程（2）在题目要求的初始条件下其范数和特征根均大于1，即采用迭代法不能收敛于方程的解。这种情况在现实中意味着：若人工给定初始种群数量时，采用不同于方程（2）的解所规定的数量，在给定的繁殖率、存活率、收获量的情况下，这种人工繁殖的方法不可能趋向于稳定。解决1:关于问题1：解的不符合现实的情况，可以通过改变参数（s、b、h）来解决。下面针对年龄段5出现负值的问题提出一种最简单的解决办法：方案：在现有人工饲养条件和种群质量（对应着参数b、s）下，减少对年龄段5的个体的收集量，以保证稳定解的存在。y=0;for k=200:10:500 h(5)=k;%减小h(5

25、),即对年龄段5的人工收集量 k0=(k-190)/10;xx=mh%输出此时方程的解y(1,k0)=xx(5);%将新的新得到h(5)的值赋予yendkk=500:-10:200;plot(kk,y(1,:),grid;%作图运行结果： 对年龄段5的收获量取不同值时方程的部分解： 1.0e+004 *X(1)X(2)X(3)X(4)X(5)H(5)=2001.17720.42090.20250.07150.00863001.17720.42090.20250.0715-0.00144001.17720.42090.20250.0715-0.01145001.17720.42090.20250

26、.0715-0.0214由于年龄段5的数量x（5）仅仅可能对年龄段1通过繁殖率发生影响，而b5为0，所以改变年龄段5的情况不会对其他年龄段书数量造成任何影响。仅减小小对年龄段5的收集量，就可以使方程存在合理的稳定解；由图形还可以得到，当h(5)大约为286时，方程的解x(5)为零，此值为对年龄段5的最大收获量。可以看出，本法可以维持其他年龄段的收集两不变，且不需增加成本改善人工饲养条件和种群质量（改变参数b、s），比较经济。 解决2：解的的不收敛性，实际上比起1更加严重，因为对方程(2)的迭代不能自动收敛，对应着实际种群数量不能自发趋向稳定，这对于人工饲养造成了非常大的麻烦，因为必须保证初始的各个年龄段个体的数量分布就满足方程2解的形式 ，且由于现实中的情况是不可能严格满足方程（2）的，总会有无法控制的个体数量波动，这又要求对种群数量的随时监测和调整，