数值分析线性方程组实验报告包含代码

线性代数方程组的直接解法实验目的：线性方程组求解的直接法编程实现，实验内容：线性方程组求解的高斯消去法算法实现线性方程组求解的主元素消去法算法实现线性方程组求解的LU分解得方法算法实现线性方程组求解追赶法算法实现实验比较：高斯消去、主元素消去、LU分解都用实例这个进行比较。知识理论解线性方程组的方法大致分为直接法和迭代法。直接法是指假设计算过程中不产生舍入误差，经过有限次的运算可求的方程组精确解的方法。方程（2-1）记为：AX=b；一、高斯（Gauss）消元法(1).Gauss消元法是最基本的一种方法。先逐次消去变量，将方程组化成同解的上三角方程组。 消元过程：先逐次消去变量，将方程组化成同解的上三角方程组基本思想 回代过程：按方程的相反顺序求解三角形方程组，得到原方程组的解。(2)Gauss消去法的求解思路为：若先让第一个方程组保持不变，利用它消去其余方程组中的,使之变成一个关于变元的n-1阶方程组。按照（1）中的思路继续运算得到更为低阶的方程组。经过n-1步的消元后，得到一个三角方程。利用求解公式回代得到线性方程组的解。消元过程：第一次消元，（i=1,2,3……，n）.将（2-1）中第i个方程减去第一个方程乘以（i=1,2,3……，n），完成第一次消元，（2-2）其中：；简记为：其中：按上述方法完成n-1次消元后得到。同解的三角方程组：简记为：回代过程：按逆序逐步回代得到方程的解。（3）算法：（5）Matlab程序代码：function[RA,RB,n,X]=gaus(A,b)B=[A,b];n=length(b);RA=rank(A);RB=rank(B);zhicha=RB-RA;if(zhicha~=0)disp('因为RA≠RB，所以此方程组无解');return;endifRA==RBifRA==ndisp('因为RA=RB=n,所以此方程有唯一解');X=zeros(n,1);forp=1:n-1fork=p+1:nm=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);forq=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('因为RA=RB<n,所以此方程有无穷解');endend运行检测：>>A=[111;12-33;-183-1];b=[6;15;-15];[RA,RB,n,X]=gaus(A,b)因为RA=RB=n,所以此方程有唯一解>>RA=3RB=3n=3X=1.00002.00003.0000二、主元素法1.列主元素法消元(1)基本思想：在每次消元前，在要消去未知数的系数中找到绝对值最大的系数作为主元，通过对换行将其换到对角线上，然后进行消元.(2)消元过程：与Gauss很类似，每次对对角线换成最大的值，后面过程与Gauss基本相同。如此经过n-1步,（2-1）的增广矩阵[A,b]被化为上三角矩阵；回代过程：同Gauss算法一样回代求解。(3)算法：det<-1对于k=1,2,3,…n-1;|aik,k|=max|aik|(k≤i≤n)(i)如果aik=0,则停止计算（det(A)=0）(ii)如果aik=k,则zhuan（i）换行akj<-->aik,j(j=k,k+1,…n) bk<-->bik后边就是Gauss消元(4).Matlab程序functionX=liezhu(A,b)B=[A,b];n=length(b);RA=rank(A);RB=rank(B);ifRA==RBifRA==ndisp('因为系数矩阵与增广矩阵的秩相同且为n，此方程有唯一解')X=zeros(n,1);forp=1:n-1[y,j]=max(abs(B(p:n,p)));C=B(p,:);B(p,:)=B(j+p-1,:);B(j+p-1,:)=C;endforp=1:n-1fork=p+1:nm=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);forq=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('因为系数矩阵与增广矩阵的秩相同且小于n，方程组有无穷解')endelsedisp('因为系数矩阵与增广矩阵的秩不相同，所以此方程组无解')endend测试结果：>>clearallclcA=[111;12-33;-183-1];b=[615-15]';liezhu(A,b)>>因为系数矩阵与增广矩阵的秩相同且为n，此方程有唯一解ans=1.00002.00003.0000例2测试结果>>clearallclcA=[0.501.13.1;2.04.50.36;5.00.966.5];b=[6.00.0200.96]';liezhu(A,b)>>ans=-2.60001.00002.0000 三、直接三角分解法1.高斯消元的矩阵表示第一次消元经过K次消元后得到：消元过程：则（1）LU分解：因为A的各元素已知，可求出L和U的各元素。看A的第k行，有，再看A的第k列，有。这样就可以计算计算出U的第k行和L的第k列的元素。从k=1到k=n交替计算U和L，就能逐步计算出U和L的全部元素。分两步求解LUx=b，先解方程组Ly=b，因为L是下三角矩阵，求解只要逐步向前回代。第二步是解方程组Ux=y，因为U是上三角矩阵，用逐次向后回代的方法求出x。3.Matlab程序 function[L,U,Y,X]=Doolittle(A,b)n=length(b);L=zeros(n,n);U=zeros(n,n);%先将A进行LU分解fori=1:n%把第一列第一行求出forj=i:nifi==1A(i,j)=A(i,j);ifj>iA(j,i)=A(j,i)/A(i,i);endendend%下面计算U的第i行值forj=i:nifi~=1m=0;fork=1:i-1m=m+A(i,k)*A(k,j);endA(i,j)=A(i,j)-m;endend%下面计算L的第i列值forj=i:nifi~=1&&j>il=0;fork=1:i-1l=l+A(j,k)*A(k,i);endA(j,i)=(A(j,i)-l)/A(i,i);endendend%将LU的值分别加入fori=1:nforj=i:nU(i,j)=A(i,j);endforj=1:iifi==jL(i,j)=1;elseL(i,j)=A(i,j);endendend%LY=B，UX=Y，进行求解Ab=[L,b]Y=zeros(n,1);Y(1)=Ab(1,n+1)/Ab(1,1);fori=2:nforj=1:i-1Y(i)=Y(i)+Ab(i,j)*Y(j);endY(i)=(Ab(i,n+1)-Y(i))/Ab(i,i);endAb=[U,Y];X=zeros(n,1);X(n)=Ab(n,n+1)/Ab(n,n);fori=n-1:-1:1forj=i+1:nX(i)=X(i)+Ab(i,j)*X(j);endX(i)=(Ab(i,n+1)-X(i))/Ab(i,i);endend运行检测：>>A=[111;12-33;-183-1];b=[6;15;-15];[L,U,Y,X]=Doolittle(A,b)运行结果：>>Ab=1.0000006.000012.00001.0000015.0000-18.0000-1.40001.0000-15.0000L=1.00000012.00001.00000-18.0000-1.40001.0000U=1.00001.00001.00000-15.0000-9.0000004.4000Y=6.0000-57.000013.2000X=1.00002.00003.0000四、解三角方程组的追赶法1.设系数矩阵为三对角矩阵则方程组Ax=f称为三对角方程组。2.设矩阵A非奇异，A有Crout分解A=LU，其中L为下三角矩阵，U为单位上三角矩阵，记先依次求出L，U中的元素后，令Ux=y，先求解下三角方程组Ly=f得出y，再求解上三角方程组Ux=y。求解三对角方程组的2追赶法将矩阵三角分解的计算与求解两个三角方程组的计算放在一起，使算法吏为紧凑。其计算公式为：3.算法：1.输入a=(a2,…，an),b=(b1,…，bn),c=(c1,…，cn-1),d=(d1,…，dn),n2.对i=2,3，…，n ai/bi-1=>ai(li) bi-ci-1ai=>bi(ui) di-aidi-1=>di(yi)3.dn/bn=>dn（xn）4.对i=n-1,n-2…，1 (di-cidi+1)/bi=>di(xi)5.输出d=x,停机4.Matlab程序：clc;clear;%如果需要验证其他矩阵的解，只需更改A、B和n的值即可A=[2,1,0,0;1,3,1,0;0,1,4,-1;0,0,2,3]B=[3;5;4;5]n=4;%获得A中不为0的数存到相应数组中fori=1:nforj=1:nifi==jb(i)=A(i,j);endifi==j-1c(j-1)=A(i,j);endifj==i-1a(i)=A(i,j);endendend%将变换后的值存入相应矩阵中u(1)=b(1);fori=2:n;l(i)=a(i)/u(i-1);u(i)=b(i)-l(i)*c(i-1);endfori=1:nforj=1:nifi==jL(i,j)=1;U(i,j)=u(i);endifi==j-1U(i,j)=c(i);endifj==i-1L(i,j)=l(i);endendend%LY=B，UX=Y，进行求解disp('增广矩阵');AB=[L,B]B=zeros(n,1);B(1)=AB(1,n+1)/AB(1,1);fori=2:nforj=1:i-1B(i)=B(i)+AB(i,j)*B(j);endB(i)=(AB(i,n+1)-B(i))/AB(i,i);enddisp('解方程LY=B得');Y=BAB=[U,B];B=zeros(n,1);B(n)=AB(n,n+1)/AB(n,n);fori=n-1:-1:1forj=i+1:nB(i)=B(i)+AB(i,j)*B(j);endB(i)=(AB(i,n+1)-B(i))/AB(i,i);enddisp('方程的解为：');X=B运行结果：A=21001310014