大学数学实验五_线性代数方程组的数值解法

1 大学数学实验五 线性代数方程组的数值解法 实验报告实验报告 【实验目的】 1、学会用 MATLAB 软件数值求解线性代数方程组，对迭代法的收敛性和解 的稳定性作初步分析。 2、通过实例学习用线性代数方程组解决简化的实际问题。 【实验内容】 3 已知方程组已知方程组 Ax=b，其中，其中 ，定义为 ，定义为 试通过迭代法求解此方程组， 认识迭代法收敛的含义以及迭代初值和方程组系数矩阵性质对试通过迭代法求解此方程组， 认识迭代法收敛的含义以及迭代初值和方程组系数矩阵性质对 收敛速度的影响。实验要求：收敛速度的影响。实验要求： （1） 选取不同的初始向量选取不同的初始向量 x(0)和不同的方程组的右端项向量和不同的方程组的右端项向量 b，给定迭代误差要求，用雅，给定迭代误差要求，用雅 可比迭代法和高斯可比迭代法和高斯-赛德尔迭代法计算， 观测得到的迭代向量序列是否均收敛？若收赛德尔迭代法计算， 观测得到的迭代向量序列是否均收敛？若收 敛，记录迭代次数，分析计算结果并得出结论。敛，记录迭代次数，分析计算结果并得出结论。 （2） 取定右端向量取定右端向量 b 和初始向量和初始向量 x(0)，将，将 A 的主对角线元素成倍增长若干次，的主对角线元素成倍增长若干次，非主对角非主对角 线元素不变，每次用雅可比迭代法计算，要求迭代误差满足线元素不变，每次用雅可比迭代法计算，要求迭代误差满足| x(k+1)- x(k)|e k=k+1; xj=Bj*xj+fj; end function xg,k=GaussSeidel(A,X0,b,e) D=diag(diag(A); n=length(A); L=-(tril(A)-D); U=-(triu(A)-D); fg=(D-L)b; Bg=(D-L)U; xg=X0; k=0; while norm(A*xg-b)/norm(b)e k=k+1; xg=Bg*xg+fg; end 3 题目要求选取不同的初始向量和右端向量，并分别用雅可比迭代法和高斯-赛德尔迭代 法求解，输出最终解和迭代次数，给定的迭代误差 e 为 1e-7。初始向量分别取零向量、1 向 量和1, 2, 3, , 19, 20，右端向量分别取 1 向量和1, 2, 3, , 19, 20。 解向量X为： B =ones(20,1)时，解向量 X=0.4816, 0.5734, 0.6328, 0.6521, 0.6609, 0.6643, 0.6657, 0.6663, 0.6665, 0.6666, 0.6666, 0.6665, 0.6663, 0.6657, 0.6643, 0.6609, 0.6521, 0.6328, 0.5734, 0.4816 B =1:20时，解向量 X=0.7247, 1.3444, 2.0071, 2.6690, 3.3344, 4.0004, 4.6668, 5.3333, 5.9998, 6.6661, 7.3320, 7.9967, 8.6585, 9.3133, 9.9501, 10.5455, 11.0251, 11.2817, 10.6973, 9.3897 雅可比迭代法雅可比迭代法 迭代次数 k X0=zeros(20,1) X0=ones(20,1) X0=1:20 B =ones(20,1) 23 22 27 B =1:20 23 23 22 高斯高斯-赛德尔迭代法赛德尔迭代法 迭代次数 k X0=zeros(20,1) X0=ones(20,1) X0=1:20 B =ones(20,1) 14 14 17 B =1:20 15 15 14 对雅可比迭代法和高斯-赛德尔迭代法的函数略作修改，以便输出图形，观察迭代向量 序列是否收敛。 修改后的修改后的 Jacobi 函数函数 A=Amatrix(20,3,-1/2,-1/4); %根据需要删除注释符号 X0=zeros(20,1); %初始向量（零向量） %X0=ones(20,1); %初始向量（1向量） B=ones(20,1); %右端向量（1向量） %X0=1:20; %初始向量（1, 2, 3, , 19, 20） %B=1:20; %右端向量（1, 2, 3, , 19, 20） X,k= Jacobi(A,X0,B,1e-7); %X,k= GaussSeidel (A,X0,B,1e-7); X k function xj,k,P=Jacobi(A,X0,b,e) D=diag(diag(A); n=length(A); L=-(tril(A)-D); U=-(triu(A)-D); fj=Db; Bj=D(L+U); xj=X0; k=0; P=X0; while norm(A*xj-b)/norm(b)e k=k+1; xj=Bj*xj+fj; P=P,xj; end 4 修改后的 Jacobi 函数多输出了矩阵 P， 矩阵 P 可视为一个行向量， 其每个元素均为迭代 k 次后得到的 xk。这样以 k 为横轴，解向量为纵轴，可输出图形观察 xk是否收敛。函数 GaussSeidel 也需作同样修改，修改后的函数在此不再赘述。 输出图形输出图形直观地观察迭代向量序列是否收敛直观地观察迭代向量序列是否收敛 雅可比迭代法雅可比迭代法 B =ones(20,1) B =1:20 X0=zeros(20,1) X0=ones(20,1) X0=1:20 A=Amatrix(20,3,-1/2,-1/4); %根据需要删除注释符号 X0=zeros(20,1); %初始向量（零向量） %X0=ones(20,1); %初始向量（1向量） B=ones(20,1); %右端向量（1向量） %X0=1:20; %初始向量（1, 2, 3, , 19, 20） %B=1:20; %右端向量（1, 2, 3, , 19, 20） X,k,P= Jacobi(A,X0,B,1e-7); % X,k,P= GaussSeidel(A,X0,B,1e-7); hold on; grid on; plot(0:k,P) 5 高斯高斯-赛德尔迭代法赛德尔迭代法 B =ones(20,1) B =1:20 X0=zeros(20,1) X0=ones(20,1) X0=1:20 由于给定的矩阵 A 是严格对角占优的，所以雅可比迭代和高斯-赛德尔迭代均收敛。从 上两表中的图片也可看出，得到的迭代向量序列都是收敛的。 结果讨论：结果讨论： 1、比较雅可比迭代法和高斯、比较雅可比迭代法和高斯-赛德尔迭代法，发现同样条件下，雅可比方法的迭代次数多于赛德尔迭代法，发现同样条件下，雅可比方法的迭代次数多于 高斯高斯-赛德尔方法，说明后者具有更好的收敛速度。赛德尔方法，说明后者具有更好的收敛速度。 2、给定的初始向量和最终的解向量越接近，迭代收敛越快，迭代次数越少。给定的初始向量和最终的解向量越接近，迭代收敛越快，迭代次数越少。 取定右端向量 b 为1:20， 初始向量为 1 向量， 改变矩阵 A 的主对角线元素分别为 3、 6、 9、 81、300，分别用雅可比迭代法和高斯-赛德尔迭代法计算，记录下迭代次数，并绘出解向量-迭 代次数的 -k 曲线。 6 输出结果输出结果 迭代次数迭代次数 k 雅可比迭代法雅可比迭代法 高斯高斯-赛德尔迭代法赛德尔迭代法 n=3 k=16 k=11 n=6 k=8 k=6 n=9 k=6 k=5 n=81 k=4 k=3 A=Amatrix(20,3,-1/2,-1/4); %修改数字3即可改变A的主对角线元素 X0=ones(20,1); B=1:20; X,k,P= Jacobi(A,X0,B,1e-5); %可改用GaussSeidel函数 hold on; grid on; plot(0:k,P) k 7 n=300 k=3 k=3 结果讨论：结果讨论： 1、 随着矩阵随着矩阵 A 主对角线元素的增大，主对角线元素的增大，两种迭代方法的收敛速度均有加快，迭代次数越来越两种迭代方法的收敛速度均有加快，迭代次数越来越 小。小。 2、 当矩阵当矩阵 A 主对角线元素增大到一定程度后，无论它再怎样增大，收敛速度不会再主对角线元素增大到一定程度后，无论它再怎样增大，收敛速度不会再明显明显加加 快，迭代次数不会再减小。快，迭代次数不会再减小。 3、 相同条件下，雅可比迭代法不如高斯相同条件下，雅可比迭代法不如高斯-赛德尔迭代法收敛得快，但当矩阵赛德尔迭代法收敛得快，但当矩阵 A 主对角线元主对角线元 素增大到一定程度后，两种方法的迭代次数相同。素增大到一定程度后，两种方法的迭代次数相同。 4 一年生植物的繁殖问题一年生植物的繁殖问题 若一棵植物秋季产种的平均数为若一棵植物秋季产种的平均数为 c，种子能够活过冬天的比例为，种子能够活过冬天的比例为 b，1 岁的种子能在春岁的种子能在春 季发芽的比例为季发芽的比例为 a1， 2 岁的种子能在春天发芽的比例为岁的种子能在春天发芽的比例为 a2。 假定种子最多可以活动两个冬天。 假定种子最多可以活动两个冬天。 如果在繁衍过程中， 我们知道了某一年植物的数量， 并希望若干年后它的数量达到给定的规如果在繁衍过程中， 我们知道了某一年植物的数量， 并希望若干年后它的数量达到给定的规 模，就需要在原模型的基础上改进，求模，就需要在原模型的基础上改进，求第二年（及第二年（及以后诸年以后诸年）这种植物的数量。这种植物的数量。假设假设 c=10， a1=0.5，a2=0.25，b=0.20，第一年有，第一年有 50 棵植物，且第棵植物，且第 50 年后有年后有 600 棵植物，试用追赶法、棵植物，试用追赶法、 稀疏系数矩阵和满矩阵求解。若稀疏系数矩阵和满矩阵求解。若 b 有有 10%误差，估计对结果的影响。误差，估计对结果的影响。 模型：模型： 已知某年该植物的数量为 x0，记第 k 年的植物数量为 xk，那么有 xk + pxk-1 + qxk-2 = 0 (k = 2, 3, , n) 其中 p = -a1bc，q = -a2b(1-a1)bc。若要求 n 年后数量达到 xn，则 Ax = b 其中 ， ， 8 用稀疏系数矩阵求解。用稀疏系数矩阵求解。 输出的植物数量-年数曲线如下。 输出的植物的数量输出的植物的数量随时间变化的数据随时间变化的数据。 第第 1 年年 50.0000 第第 11 年后年后 97.3970 第第 22 年后年后 162.6504 第第 1 年后年后 61.5030 第第 12 年后年后 102.0451 第第 23 年后年后 170.4124 第第 2 年后年后 64.0030 第第 13 年后年后 106.9149 第第 24 年后年后 178.5450 第第 3 年后年后 67.0782 第第 14 年后年后 112.0172 第第 25 年后年后 187.0656 第第 4 年后年后 70.2783 第第 15 年后年后 117.3629 第第 26 年后年后 195.9928 第第 5 年后年后 73.6322 第第 16 年后年后 122.9638 第第 27 年后年后 205.3461 第第 6 年后年后 77.1462 第第 17 年后年后 128.8319 第第 28 年后年后 215.1458 第第 7 年后年后 80.8278 第第 18 年后年后 134.9801 第第 29 年后年后 225.4131 第第 8 年后年后 84.6851 第第 19 年后年后 141.4217 第第 30 年后年后 236.1703 第第 9 年后年后 88.7265 第第 20 年后年后 148.1707 第第 31 年后年后 247.4410 第第10年后年后 92.9607 第第 21 年后年后 155.2418 第第 32 年后年后 259.2495 c=10;a1=0.5;a2=0.25;b=0.20; p=-a1*b*c;q=-a2*b*(1-a1)*b*c; A1=p*eye(49,49); A2=zeros(48,1),eye(48,48);zeros(1,49); A3=zeros(1,48);q*eye(48,48),zeros(49,1); A=A1+A2+A3; B=zeros(49,1);B(1,1)=-q*50;B(49,1)=-600; A1=sparse(A);B1=sparse(B);X=A1B1; X1=50;X;600 hold on; grid on; plot(0:50,X1) 9 第第33年后年后 271.6216 第第 39 年后年后 359.2874 第第 45 年后年后 475.2474 第第34年后年后 284.5840 第第 40 年后年后 376.4335 第第 46 年后年后 497.9275 第第35年后年后 298.1651 第第 41 年后年后 394.3979 第第 47 年后年后 521.6898 第第36年后年后 312.3943 第第 42 年后年后 413.2196 第第 48 年后年后 546.5862 第第37年后年后 327.3026 第第 43 年后年后 432.9395 第第 49 年后年后 572.6707 第第38年后年后 342.9223 第第 44 年后年后 453.6005 第第 50 年后年后 600.0000 用满矩阵求解。用满矩阵求解。 输出的结果输出的结果： X 为零向量，表明用稀疏矩阵和用满矩阵算出来的结果一致。 t1= 4.3665 10-4s，t2= 0.5849s，表明 MATLAB 软件用满矩阵计算所用的时间比用稀疏矩 阵计算所用的时间长许多。 用追赶法求解。用追赶法求解。 clear all; c=10;a1=0.5;a2=0.25;b=0.20; p=-a1*b*c;q=-a2*b*(1-a1)*b*c; A1=p*eye(49,49); A2=zeros(48,1),eye(48,48);zeros(1,49); A3=zeros(1,48);q*eye(48,48),zeros(49,1); A=A1+A2+A3; B=zeros(49,1);B(1,1)=-q*50;B(49,1)=-600; A1=sparse(A);A2=full(A1);B1=sparse(B);B2=full(B1); tic;X1=A1B1;t1=toc; tic;X2=A2B2;t2=toc; X=X1-X2; X t1 t2 c=10;a1=0.5;a2=0.25;b=0.20; p=-a1*b*c;q=-a2*b*(1-a1)*b*c; A1=p*eye(49,49); A2=zeros(48,1),eye(48,48);zeros(1,49); A3=zeros(1,48);q*eye(48,48),zeros(49,1); A=A1+A2+A3; B=zeros(49,1);B(1,1)=-q*50;B(49,1)=-600; L,U=lu(A); %矩阵A的LU分解 y(1,1)=-q*50; 10 输出的结果：输出的结果： Delta 为零向量，表明用追赶法和用直接法算出来的结果一致。 估计估计 b 的扰动对结果的影响。的扰动对结果的影响。 利用 cond(A)语句求矩阵 A 的条件数，得到的结果是 27.3659 1，表明线性方程组 Ax = b 是病态的，可以预计 b 的扰动对解的影响很大。 下面改变 b 的值，输出植物数量-年数曲线，直观地观察 b 的扰动对结果的影响。程序 用的是第 8 页上的程序（改变 b 的值即可） ，在此不再赘述。 （b=0.2 90%0.2 100%的情况） for i=2:49 y(i,1)=B(i,1)-L(i,i-1)*y(i-1,1); end x(49,1)=y(49,1)/U(49,49);i=48; while (i=1) x(i,1)=(y(i,1)-U(i,i+1)*x(i+1,1)/U(i,i);i=i-1; end x1=50;x;600; X=AB; % 直接法 X1=50;X;600; Delta=X1-x1 % 比较直接法和追赶法求得的解 11 （b=0.2 100%0.2 110%的情况） 可见，可见，b 的扰动对结果的影响很大，的扰动对结果的影响很大，故在测定故在测定 b 的值时要尽量准确，并可采用多次估计的值时要尽量准确，并可采用多次估计 取平均值的方法减小误差。取平均值的方法减小误差。 5 投入产出分析投入产出分析 设国民经济由农业、 制造业、 服务业三个部门构成。 已知某年它们之间的投入产出关系、设国民经济由农业、 制造业、 服务业三个部门构成。 已知某年它们之间的投入产出关系、 外部需求、初始投入等如下表所示。外部需求、初始投入等如下表所示。 产出产出 投入投入 农业农业 制造业制造业 服务业服务业 外部需求外部需求 总产出总产出 农业农业 15 20 30 35 100 制造业制造业 30 10 45 115 200 服务业服务业 20 60 0 70 150 初始投入初始投入 35 110 75 总投入总投入 100 200 150 单位：亿元单位：亿元 假定每个部门的产出与它的投入成正比， 由上表可确定这三个部门的投入产出表， 如下表所假定每个部门的产出与它的投入成正比， 由上表可确定这三个部门的投入产出表， 如下表所 示。示。 产出产出 投入投入 农业农业 制造业制造业 服务业服务业 农业农业 0.15 0.1 0.2 制造业制造业 0.3 0.05 0.3 服务业服务业 0.2 0.3 0 12 （1）如果今年对农业、制造业、服务业的外部需求分别为）如果今年对农业、制造业、服务业的外部需求分别为 50,150,100 亿元，问这三个部门亿元，问这三个部门 的总产出分别为多少？的总产出分别为多少？ （2）如果三个部门的外部需求分别增加）如果三个部门的外部需求分别增加 1 个单位，问它们的总产出分别个单位，问它们的总产出分别增加多少？增加多少？ 模型：模型： 设有 n 个部门， 记某段时间内第 i 个部门的总产出为 xi， 其中对第 j 个部门的投入为 xij， 外部需求为 di，则有 投入系数记为 aij，有 故有 记投入系数矩阵 A=(aij)n n， 产出向量 x=(x1, x2, x3, xn)T， 需求向量 d=(d1, d2, d3, dn)T，则有 x = Ax + d 即 (I A)x = d 通过求解线性方程组解答第一小问。通过求解线性方程组解答第一小问。 输出的结果为 x = 139.2801, 267.6056, 208.1377T。 结论：结论： 农业、制造业、服务业部门的总产出分别为农业、制造业、服务业部门的总产出分别为 139.2801，267.5056，208.1377 亿元。亿元。 通过求解线性方程组解答第二小问。通过求解线性方程组解答第二小问。 已知 (I A)x = d 故有 (I A) x = d A=0.15,0.1,0.2;0.3,0.05,0.3;0.2,0.3,0; I=eye(3,3); d=50;150;100; x=(I-A)d; x 13 假设农业的外部需求增加 1 个单位。 输出的结果为 x = 1.3459, 0.5634, 0.4382T。 结论：结论： 农业的外部需求增加农业的外部需求增加 1 个单位时，农业部门的产出增加个单位时，农业部门的产出增加 1.3459 亿元，制造业部门的产亿元，制造业部门的产 出增加出增加 0.5634 亿元，服务业部门的产出增加亿元，服务业部门的产出增加 0.4382 亿元。亿元。 假设制造业的外部需求增加 1 个单位。 （程序与本页上方的程序类似