大连理工大学矩阵与数值分析上机作业代码

《矩阵与数值分析》上机作业实验报告陈晓学号：21204047机械工程学院教学班号：02任课教师：金光日习题1111给定n阶方程组Axb=，其中6186186186A⎛⎞⎜⎟⎜⎟⎜⎟=⎜⎟⎜⎟⎜⎟⎝⎠⋱⋱⋱，7151514b⎛⎞⎜⎟⎜⎟⎜⎟=⎜⎟⎜⎟⎜⎟⎝⎠⋮则方程组有解(1,1,,1)Tx=⋯。对10n=和84n=，分别用Gauss消去法和列主元消去法解方程组，并比较计算结果。1.11.11.11.1程序程序程序程序（1）高斯消元法n=10时，[A,b]=CreateMatrix(10)A=6100000000861000000008610000000086100000000861000000008610000000086100000000861000000008610000000086b=7151515151515151514x=GaussMethod(A,b)x=1.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000n=84时，[A,b]=CreateMatrix(84)x=GaussMethod(A,b)1.0e+008*0.00000.00000.00000.00000.00000.00000.00000.0000……-0.00010.0002-0.00030.0007-0.00130.0026-0.00520.0105-0.02090.0419-0.08360.1665-0.33030.6501-1.25822.3487-4.02635.3684（2）列主元消去法n=10时，[A,b]=CreateMatrix(10)x=MainElementMethod(A,b)x=1.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000n=84时，[A,b]=CreateMatrix(84)x=MainElementMethod(A,b)x=1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000……1.00001.00001.00001.00001.21.21.21.2计算结果分析计算结果分析计算结果分析计算结果分析对于高斯消元法，当n=10时，得到的计算结果为x=（1.0000，1.0000，1.0000，1.0000，……，1.0000）T，恰好为精确值；但当n=84时，得到的计算结果却为x=(0.0000,0.0000,0.0000,……，-0.0001,0.0002，……，2.3487，-4.0263，5.3684）T,完全偏离了精确值。对于列主元消去法，当n=10时，得到的计算结果为x=（1.0000，1.0000，1.0000，1.0000，……，1.0000）T，恰好为精确值；当n=10时，得到的计算结果为x=（1.0000，1.0000，1.0000，1.0000，……，1.0000）T，恰好为精确值。由于使用高斯消元法时可能会出现小主元作除数的情况，虽然在求解系数矩阵为10阶的方程时，由于步数较少，尚能求得精确解，但在求解求解系数矩阵为84阶的巨型方程时，由于求解步数太多，使得小主元作除数的情况累加，最终导致“大数除小数”的情况发生，因而所得解严重偏离精确解。而使用列主元消去法时，由于有效地避免了这一不利情况，因而无论是在求解系数矩阵为10阶还是84阶的方程时，都能得到精确解。1.31.31.31.3MMMM文件文件文件文件1.3.11.3.11.3.11.3.1CreateMatrixCreateMatrixCreateMatrixCreateMatrix.m.m.m.mfunction[A,b]=CreateMatrix(n)%用于存放习题1的题目信息，并建立构造题目中矩阵的函数%对矩阵A赋值A1=6*ones(1,n);A2=ones(1,n-1);A3=8*ones(1,n-1);A=diag(A1)+diag(A2,1)+diag(A3,-1);%对向量b赋值b=15*ones(n,1);b(1)=7;b(n)=14;1.3.21.3.21.3.21.3.2GaussMethod.mGaussMethod.mGaussMethod.mGaussMethod.mfunctionx=GaussMethod(A,b)%高斯消元法返回方程Ax=b的解%x是方程的解%A为方程的系数矩阵%b为方程的常数矩阵n=length(b);%定义单位下三角矩阵元胞数组，并将全部数组元素赋值为n阶单位阵l=cell(1,n-1);l(:)={eye(n)};U=A;fori=1:n-1forj=i+1:nl{i}(j,i)=-U(j,i)/U(i,i);endU=l{i}*U;end%计算L和UL=inv(l{1});fori=2:n-1L=L/l{i};endx=U\(L\b);1.3.31.3.31.3.31.3.3MainElementMethodMainElementMethodMainElementMethodMainElementMethod.m.m.m.mfunctionx=MainElementMethod(A,b)%列主元消去法返回方程Ax=b的解%x是方程的解%A为方程的系数矩阵%b为方程的常数矩阵n=length(b);%定义单位下三角矩阵元胞数组，并将全部数组元素赋值为n阶单位阵l=cell(1,n-1);l(:)={eye(n)};p=cell(1,n-1);p(:)={eye(n)};U=A;fori=1:n-1%找出列主元所在行，并交换首行，得到Pi[Max,row]=max(U(i:n,i));row=row+i-1;mid=p{i}(row,:);p{i}(row,:)=p{i}(i,:);p{i}(i,:)=mid;U=p{i}*U;forj=i+1:nl{i}(j,i)=-U(j,i)/Max;endU=l{i}*U;end%计算LL=eye(n);fori=1:n-2forj=i+1:n-1l{i}=p{j}*l{i}/p{j};endL=l{i}*L;endL=inv(L*l{n-1});%计算PP=eye(n);fori=1:n-1P=p{i}*P;end%计算xx=U\(L\P*b);习题2222设有方程组Axb=，其中2020AR×∈，30.50.250.530.50.250.50.50.250.530.50.250.53A−−⎛⎞⎜⎟−−⎜⎟⎜⎟−−⎜⎟=⎜⎟⎜⎟−−⎜⎟−−⎜⎟⎜⎟−−⎝⎠⋱⋱⋱⋱⋱⋱⋱⋱⋱⋱⋱⋱(a)选取不同的初始向量(0)x和不同的右端向量b，给定迭代误差要求，用Jacobi和Gauss-Seidel迭代法计算，观测得出的迭代向量序列是否收敛。若收敛，记录迭代次数，分析计算结果并得出你的结论。(b)选定初始向量初始向量(0)x和不同的右端向量b，如取(0)0,,(1,1,1)TxbAee===⋯。将A的主对角线元素成倍增长若干次，非主对角元素不变，每次用Jacobi法计算，要求迭代误差满足(1)()610kkxx+−∞−，比较收敛速度，分析现象并得出你的结论。2.12.12.12.1程序程序程序程序（1）取(0)x=0，b=（1,1，……，1,1）T,迭代误差取默认值105−，迭代次数上限取默认值50，使用Jacobi法进行迭代。test2b=ones(20,1)x0=zeros(20,1)[x,n]=JacobiMethod(A,b,x0)x=0.27660.23270.21590.22230.22270.22210.22220.22220.22220.22220.22220.22220.22220.22220.22210.22270.22230.21590.23270.2766n=17（2）取(0)x=0，b=（1,1，……，1,1）T,迭代误差取默认值105−，迭代次数上限取默认值50，使用Gauss-Seidel法进行迭代。test2b=ones(20,1)x0=zeros(20,1)[x,n]=Gauss_SeidelMethod(A,b,x0)x=0.27660.23270.21590.22230.22280.22210.22220.22220.22220.22220.22220.22220.22220.22220.22210.22280.22230.21590.23270.2766n=9（3）取(0)x=0，b=（1,0,1,0,……，1,0）T,迭代误差取默认值105−，迭代次数上限取默认值50，使用Jacobi法进行迭代。test2b=[1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0]'x0=zeros(20,1)[x,n]=JacobiMethod(A,b,x0)x=0.3238-0.09850.3116-0.08810.3110-0.08890.3111-0.08890.3111-0.08890.3111-0.08890.3111-0.08890.3111-0.08820.3105-0.09570.3313-0.0472n=16（4）取(0)x=0，b=（1,0,1,0,……，1,0）T,迭代误差取默认值105−，迭代次数上限取默认值50，使用Gauss-Seidel法进行迭代。test2b=[1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0]'x0=zeros(20,1)[x,n]=Gauss_SeidelMethod(A,b,x0)x=0.3238-0.09850.3115-0.08810.3110-0.08890.3111-0.08890.3111-0.08890.3111-0.08890.3111-0.08890.3111-0.08820.3104-0.09570.3313-0.0472n=9（5）取(0)x=（1,1，……，1,1）T，b=（1,1，……，1,1）T,迭代误差取默认值105−，迭代次数上限取默认值50，使用Jacobi法进行迭代。test2b=zeros(20,1)x0=zeros(20,1)[x,n]=JacobiMethod(A,b,x0)x=0.27660.23270.21590.22230.22280.22210.22220.22220.22220.22220.22220.22220.22220.22220.22210.22280.22230.21590.23270.2766n=19（6）取(0)x=（1,1，……，1,1）T，b=（1,1，……，1,1）T,迭代误差取默认值105−，迭代次数上限取默认值50，使用Gauss-Seidel法进行迭代。test2b=zeros(20,1)x0=zeros(20,1)[x,n]=Gauss_SeidelMethod(A,b,x0)x=0.27660.23270.21590.22230.22280.22210.22220.22220.22220.22220.22220.22220.22220.22220.22210.22280.22230.21590.23270.2766n=10（7）取(0)x=0，b=Ae,将A的主对角线元素分别增大1,2,3……50倍，迭代误差取106−，迭代次数上限取30，使用Jacobi法进行迭代，比较每次迭代次数。n=test2_2()n=Columns1through3021119887766666666555555555555555Columns31through50555555555555555555552.22.22.22.2计算结果分析计算结果分析计算结果分析计算结果分析由（1）~（6）中计算结果可见，一