线性方程组的直接解法迭代解法

1广东金融学院实验报告课程名称：数值分析实验编号及实验名称线性方程组的直接解法/迭代解法系别应用数学系姓名学号班级实验地点实验日期实验时数2指导教师同组其他成员无成绩2实验目的及要求实验目的：题一：通过数值实验，从中体会解线性方程组选主元的必要性和LU分解法的优点，以及方程组系数矩阵和右端向量的微小变化对解向量的影响。比较各种直接接法在解线性方程组中的效果；题二：认识各种迭代法收敛的含义、影响各迭代法收敛速度的因素。实验要求：题一：（1）在MATLAB中编写程序用列主元高斯消去法和LU分解求解上述方程组，输出Ax=b中矩阵A及向量b和A=LU分解中的L及U，detA及解向量x.（2）将方程组中的2.099999改为2.1，5.900001改为5.9，用列主元高斯消去法求解变换后的方程组，输出解向量x及detA，并与（1）中的结果比较。（3）用MATLAB的内部函数inv求出系数矩阵的逆矩阵，再输入命令x=inv(A)*b,即可求出方程组的解。请与列主元高斯消去法和LU分解法求出的解进行比较，体会选主元的方法具有良好的数值稳定性。用MATLAB的内部函数det求出系数行列式的值，并与（1）、（2）中输出的系数行列式的值进行比较。（4）比较以上各种直接解法在解线性方程组中的效果。题二：（1）选取不同的初始向量X及右端向量b，给定迭代误差要求，用Jacobi迭代法和Gauss-Seidel迭代法求解，观察得到的序列是否收敛？若收敛，记录迭代次数，分析计算结果并得出你的结论。列出算法清单。（2）用SOR迭代法求上述方程组的解，松弛系数取的不同的三个值，在时停止迭代，记录迭代次数，分析计算结果与松弛系数的关系并得出你的结论。（3）用MATLAB的内部函数inv求出系数矩阵的逆矩阵，再输入命令x=inv(A)*b,即可求出上述各个方程组的解，并与上述三种方法求出的解进行比较。请将比较结果列入下表。方程组的解1x2x3x迭代次数误差精确解Jacibi解法Gause-seidel解法SOR解法)0(215)1()(10kkXX3实验环境及相关情况（包含使用软件、实验设备、主要仪器及材料等）1.Win72.Matlab7.0实验内容及步骤（包含简要的实验步骤流程）实验内容：题一：解下列线性方程组题二：研究解线性方程组迭代法的收敛性、收敛速度以及SOR方法中最佳松弛因子的选取问题，用迭代法求解，其中512184114A，15217b实验结果（包括程序或图表、结论陈述、数据记录及分析等，可附页）题一：直接解法解线性方程组(1)列主元高斯消去法与LU分解求解列主元高斯消去法:编写matalab程序（见附录gaosi.m），输出矩阵𝐴=[10.000−7.0000.0002.50000.0001−5.000−1.50.0000.0000.0000.0006.0000022.3000.0005.080]向量𝑏=[818.3005.0800]解向量：𝑋=(0，−1，1，1)𝑇其中系数行列式的值det(A)=762.00009LU分解求解：编写matalab程序（见附录zhjLU.m和LU.m），执行输出:15900001.582012151526099999.23107104321xxxxbAxbAx4𝐿=[1.00000.0000−3.0001.000000.00000.00000.00000.00000.5000−25000000.2000−24000001.00000.00000.96001.0000]𝑈=[10.0000−7.00000.0000−0.0000010.00001.00000.00002.30000.00000.00000.00000.00001500000057500000.00005.0800]在matlab命令窗口输入L*U，可以得到A=L*U,即分解结果正确。解向量:𝑋=(0，−1，1，1)𝑇其中系数行列式的值det(A)=-762.00009(2)将方程组中的2.099999改为2.1，5.900001改为5.9,使用列主元法求解后与（1）比较结果修改程序（见附录gaosi2.m），输出解向量：𝑋=(0，−1，1，1)𝑇其中系数行列式的值det(A)=762，与(1)中结果对比：方法解向量|det(A)|列主元高斯消去法𝑋=(0，−1，1，1)𝑇762.00009微小变化后列主元法𝑋=(0，−1，1，1)𝑇762从上表可以看出，微小变化后运用列主元高斯消去法输出解向量和|det(A)|几乎不变，可以初步认为列主元高斯消去算法的稳定性较好,且说运用该方法计算这类方程有效。同时也可以说明这个方程是良态的。(3)内部函数inv直接求解直接使用inv函数编写matlab程序(如附录inv.m)，输出向量解：𝑋=(0，−1，1，1)𝑇其中det(A)=-762.00009。结合(1)、(2)，可列出全部结果表如下：方法解向量|det(A)|列主元高斯消去法𝑋=(0，−1，1，1)𝑇762.00009LU分解求解𝑋=(0，−1，1，1)𝑇762.00009微小变化后列主元法𝑋=(0，−1，1，1)𝑇762Inv内部函数法𝑋=(0，−1，1，1)𝑇762.00009由上表可看出，列主元高斯消去法、LU分解求解和Inv内部函数法求得的解向量都一致，都可以较精确的求解原线性方程，在原方程组系数发生微小变化时，列主元法仍能稳定计算出方程组的解，说明其具有一定的稳定性。系数行列式的值（|det(A)|）几乎都相同。(4)比较以上各种直接解法在解线性方程组中的效果使用matlab内部函数计算得出的解向量x=inv(A)*b，即为该方程组的解。与列主元法和LU分解法得出的结果进行比较可知，他们都是等同的。经过系数微小变化后列主元法仍能稳定的输出解向量，说明其具有良好的数值稳定性。而LU法具有较简介的优点，也可以准确输出方程的解。5------------------------------------------------------------题一程序附录-----------------------------------------------(1)%gaosi.m文件%高斯消元法clcclearformatlonggA=[10-7018;-32.099999625.900001;5-15-15;21021;]%10最大，第一次不用选主元A(2,:)=A(2,:)-A(1,:)*A(2,1)/A(1,1)A(3,:)=A(3,:)-A(1,:)*A(3,1)/A(1,1)A(4,:)=A(4,:)-A(1,:)*A(4,1)/A(1,1)tempo=A(3,:),A(3,:)=A(2,:),A(2,:)=tempo%第二次选主元A(3,:)=A(3,:)-A(2,:)*A(3,2)/A(2,2)A(4,:)=A(4,:)-A(2,:)*A(4,2)/A(2,2)%6最大,不同进行第三次选主元A(4,:)=A(4,:)-A(3,:)*A(4,3)/A(3,3)x(4)=A(4,5)/A(4,4)x(3)=(A(3,5)-A(3,4)*x(4))/A(3,3)x(2)=(A(2,5)-A(2,4)*x(4)-A(2,3)*x(3))/A(2,2)x(1)=(A(1,5)-A(1,4)*x(4)-A(1,3)*x(3)-A(1,2)*x(2))/A(1,1)deta=A(1,1).*A(2,2).*A(3,3).*A(4,4)%det(A)%zhjLU文件%实现对矩阵A的LU分解的函数function[L,U]=zhjLU(A)[n,n]=size(A);L=zeros(n,n);U=zeros(n,n);fori=1:nL(i,i)=1;endfork=1:nforj=k:nU(k,j)=A(k,j)-sum(L(k,1:k-1).*U(1:k-1,j)');endfori=k+1:nL(i,k)=(A(i,k)-sum(L(i,1:k-1).*U(1:k-1,k)'))/U(k,k);endend%LU.m文件%LU分解与求解clcclear6A=[10-701;-32.09999962;5-15-1;2102;];[L,U]=zhjLU(A)%求LU[n,m]=size(A);b=[85.90000151]';y(1)=b(1);fori=2:nforj=1:i-1b(i)=b(i)-L(i,j)*y(j);endy(i)=b(i);endx(n)=y(n)/U(n,n);fori=(n-1):-1:1forj=n:-1:i+1y(i)=y(i)-U(i,j)*x(j);endx(i)=y(i)/U(i,i);enddet=U(1,1)*U(2,2)*U(3,3)*U(4,4)------(2)%gaosi2.m文件%高斯消元法clcclearformatlonggA=[10-7018;-32.1625.9;5-15-15;21021;]%10最大，第一次不用选主元A(2,:)=A(2,:)-A(1,:)*A(2,1)/A(1,1)A(3,:)=A(3,:)-A(1,:)*A(3,1)/A(1,1)A(4,:)=A(4,:)-A(1,:)*A(4,1)/A(1,1)tempo=A(3,:),A(3,:)=A(2,:),A(2,:)=tempo%第二次选主元A(3,:)=A(3,:)-A(2,:)*A(3,2)/A(2,2)A(4,:)=A(4,:)-A(2,:)*A(4,2)/A(2,2)%6最大,不同进行第三次选主元A(4,:)=A(4,:)-A(3,:)*A(4,3)/A(3,3)x(4)=A(4,5)/A(4,4)x(3)=(A(3,5)-A(3,4)*x(4))/A(3,3)x(2)=(A(2,5)-A(2,4)*x(4)-A(2,3)*x(3))/A(2,2)x(1)=(A(1,5)-A(1,4)*x(4)-A(1,3)*x(3)-A(1,2)*x(2))/A(1,1)deta=A(1,1).*A(2,2).*A(3,3).*A(4,4)%det(A)------(3)%inv.m文件%inv函数解法clc7clearformatlonga=[10-701;-32.09999962;5-15-1;2102;]b=[85.90000151]'x=inv(a)*bdet(a)-----------------------------------题一分割线----------------------------------题二：迭代解法解线性方程组（1）选取不同的初始向量X及右端向量b，给定迭代误差要求，用Jacobi迭代法和Gauss-Seidel迭代法求解，观察得到的序列是否收敛？若收敛，记录迭代次数，分析计算结果并得出你的结论。列出算法清单。Jacobi迭代法：Step1:判断收敛编辑程序如下：-----雅可比迭代收敛判断--------clcclearA=[4-11;4-81;-215;];B=[7-2115]';n=length(A);D=zeros(n);L=zeros(n);U=zeros(n);fori=1:nforj=1:nifi==jD(i,j)=A(i,j)elseifijL(i,j)=-A(i,j);elseU(i,j)=-A(i,j);endendendJ=inv(D)*(L+U)PR=max(abs(eig(J)))---------------------------结果：得雅可比迭代矩阵谱半径()=(𝐿𝑈)=0.3351,即对于矩阵A雅可比迭代法收敛。Step2:迭代次数)0(8编写程序如下：------雅可比迭代求解-------clcclearA=[4-