解线性方程组的直接法的Matlab实现1.docVIP

学校编码：* 分类号 密级 学号：* UDC 本科毕业论文(设计) 解线性方程组的直接法的Matlab实现 学生姓名：*所属院部：*专 业：*指导教师：*2014年 5 月 14 日摘要：给出用MATLAB解线性方程组的各种方法,用MATLAB直接操作,不用编程,便可立即求出线性方程组的解.方法直观、简便、速度快,具有较强的实用性,另外提供了Jacobi迭代法程序.关键字：线性方程组 数值解 程序设计 MATLAB Jacobi迭代法 数据结构1 引言在自然科学与社会科学的研究中，常常需要求解线性代数方程组，如实验数据的曲线、曲面的拟合和用差分法或有限元法解偏微分方程等都要用到线性代数方程组的求解。由于从不同的问题导出的线性代数方程组的系数矩阵不同，比如：矩阵阶数的大小、矩阵中的非零元稠密情况等，粗略地，系数矩阵可以分为低阶稠密矩阵和大型稀疏矩阵。关于线性代数程组的求解，主要分为直接法和迭代法两大类，在理论上，用直接法可以通过有限步的计算得到精确解，而迭代法是通过逐次迭代逼近来求得近似解，实际上，由于舍入误差的影响，由直接法得到的解也不精确，因此，在某些需要高精度解的问题中，常常把由直接法得到的解再运用迭代法迭代若干步，以提高解的精度，一般地说，对于低阶稠密的线性代数方程组以及大型带形方程组的求解，采用直接法比较有效，而对于大型稀疏（非带形）方程组则用迭代法求解比较有利。当然，采用直接法，还是迭代法，还是直接法与迭代法交替运用，要根据具体情况确定。本章主要讨论一些最基本的直接法，并在此基础上讨论它的变形情况，大量的科技与工程实际问题，常常归结为解线性代数方程组，有关线性方程组解的存在性和唯一性在“线性代数”理论中已经作过详细介绍，本章的主要任务是讨论系数行列式不为零的n阶非齐次线性方程组Ax=b的两类主要求方法：直接法（精确法）和迭代法。对它的解法我们最熟悉的就是主元消去法,但它只是适用于A是低价稠密的矩阵,对于由工程技术中产生的大型稀疏矩阵方程组(即A的阶数n很大,但零元素较多,例如求某些偏微分方程数值解所产生的线性方程组,n104),还需利用迭代法求解。如在计算机内存和运算两方面,都可以根据A中有大量零元素的特点采用迭代法。本文将介绍两种常见的迭代:Jacobi迭代法和Gauss-Seidel迭代,并用迭代法在数学软件Matlab上实现线性方程组的解。1迭代法的基本思想迭代法是按照某种规则构造一个向量序列x(k),使其极限向量x*是Ax=b的精确解。因此,对迭代法来说一般有下面几个问题:(1)如何构造迭代序列?(2)构造的迭代法序列是否收敛?在什么情况下收敛?(3)如果收敛,收敛的速度如何?我们应该给予量的刻划,用以比较各种迭代法收敛的快慢。2 相关知识线性方程组的概念及分类线性方程组的一般形式为a11x1+a12x2+a1nxn=b1a21x1+a22x2+a2nxn=b2am1x1+am2x2+amnxn=bn(1)若记X=x1x2(x n)T,b=b1 b2(bn)TA=a11 a12a1na21 a22a2nam1 am2a mn则线性方程组(1)记为AX=b.(2)若b的元素不全为零,则称方程组(1)为非齐次线性方程组;若b的元素全为零,即b1=b2=bn=0,则AX=0.(3)并称方程组(3)为齐次线性方程组,也称作方程组(2)的导出方程组,称(A b)=a11 a12a1nb1a21 a22a2nb2am1 am2amnb n为线性方程组(1)的增广矩阵,记作A.若在方程组(1)中,当mn, 3、算法用高斯消元法解线性方程组bAX的MATLAB程序输入的量：系数矩阵A和常系数向量b；输出的量：系数矩阵A和增广矩阵B的秩RA,RB, 方程组中未知量的个数n和有关方程组解X及其解的信息.function RA,RB,n,X=gaus(A,b) B=A b; n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica0, disp(请注意：因为RA=RB，所以此方程组无解.) returnend if RA=RB if RA=n disp(请注意：因为RA=RB=n，所以此方程组有唯一解.) X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1（1）LU分解法lu 分解法解线性方程组function x=luxiaoyuan(A,b)m,n=size(A);l u=lu(A);s=inv(l)*A,b;x=ones(m,1);for i=m:-1:1 h=s(i,m+1); for j=m:-1:1; if j=i h=h-x(j)*s(i,j); end end x(i)=h/s(i,i)end （2）高斯消元法高斯消元法的基本思想：Ax=b其对应的增广矩阵为为(A，b)（1）对换(A，b)某两行的顺序（2）(A，b)中的某行乘以一个不为零的数。（3）把(A，b)某一行乘以一个常数后加到另一行其增广矩阵变为(A，b)变换后的方程组为Ax=b与原方程组等价（同解）例：用高斯消去法找出下列方程組的解原方程组增广矩阵如果增广矩阵第一行第一列不为零，将第二，三行第一列化为零。然后第三行的第二列化为零，则其解便可求得Gauss消元法的matlab程序：% % % % -function x = gaussElim(A,b)% This subroutine will perform Gaussian elmination% on the matrix that you pass to it.% i.e., given A and b it can be used to find x,% Ax = b% A must be a squre matrixN = length(b);% Perform Gaussian Eliminationif A(1,1)=0 n_temp=find(A(:,1)=0); c=A(n_temp(1),:); A(n_temp(1),:)=A(1,:); A(1,:)=c;endfor j=2:N%the jth row m = A(j:N,j-1)/A(j-1,j-1);%column vector A(j:N,:) = A(j:N,:) - m*A(j-1,:); b(j:N) = b(j:N) - m*b(j-1); if A(j,j)=0 &j=N n_temp=find(A(j:end,j)=0); c=A(j-1+n_temp(1),:); A(j-1+n_temp(1),:)=A(j,:); A(j,:)=c; end End% Perform back substitutionif A(N,N)=0 error(the coefficient matrix of the linear equations is sigular);endx = zeros(N,1);x(N) = b(N)/A(N,N);for j=N-1:-1:1, x(j) = (b(j)-A(j,j+1:N)*x(j+1:N)/A(j,j); end% % % % -列主元消元法MATLAB程序：% % % % -function x = gaussElim_pivot_column(A,b)% This subroutine will perform Gaussian elmination with partial pivoting% on the matrix that you pass to it.% i.e., given A and b it can be used to find x,% Ax = b% A must be a squre matrixN = length(b);max_Aj=max(abs(A(:,1);n_temp=find(abs(A(:,1)=max_Aj);c=A(n_temp(1),:);A(n_temp(1),:)=A(1,:);A(1,:)=c;d=b(n_temp(1);b(n_temp(1)=b(1)b(1)=d;% Perform Gaussian Eliminationfor j=2:N%the jth row m = A(j:N,j-1)/A(j-1,j-1);%column vector A(j:N,:) = A(j:N,:) - m*A(j-1,:); b(j:N) = b(j:N) - m*b(j-1); if j=N max_Aj=max(abs(A(j:end,j); n_temp=find(abs(A(j:end,j)=max_Aj); c=A(j-1+n_temp(1),:); A(j-1+n_temp(1),:)=A(j,:); A(j,:)=c d=b(j-1+n_temp(1); b(j-1+n_temp(1)=b(j) b(j)=d; end end% Perform back substitutionif A(N,N)=0 error(the coefficient matrix of the linear equations is sigular);endx = zeros(N,1);x(N) = b(N)/A(N,N); for j=N-1:-1:1, x(j) = (b(j)-A(j,j+1:N)*x(j+1:N)/A(j,j);end% % % % -如果进一步地对增广矩阵进行初等变换，使上三角阵变换为对角矩阵，则称为Gauss-Jordan消元法。Gauss-Jordan消元法：用初等变换化系数矩阵为对角矩阵的方法称为Gauss-Jordan消元法。直接分解法（3） 高斯消元法对系数矩阵A消元的过程相当于将A分解为一个上三角阵和一个下三角阵的过程。如果直接分解A得到L和U，A=LU。这时方程Ax=b化为LUx=b。令Ux=y，则Ly=b解出来y，再由Ux=y解出来x。这就是直接分解法。LU decomposition is a matrix decomposition which writes a matrix as the product of a lower triangular matrix and an upper triangular matrix. 如果直接分解A为A=LU。，当A为单位下三角阵（对角线元素全部为1），U是上三角阵时，称为Dolittle分解；当L是下三角阵，U是单位上三角阵时，称为Courant分解。Let A be a square matrix. An LU decomposition is a decomposition of the formwhere L and U are lower and upper triangular matrices (of the same size), respectively. This means that L has only zeros above the diagonal and U has only zeros below the diagonal. For a 3*3 matrix, this becomes:An LDU decomposition is a decomposition of the formwhere D is a diagonal matrix and L and U are unit triangular matrices, meaning that all the entries on the diagonals of L and U are one.Matlab中的lu函数可直接进行矩阵的LU分解。LU LU factorization. L,U = LU(A) stores an upper triangular matrix in U and a psychologically lower triangular matrix (i.e. a product of lowertriangular and permutation matrices) in L, so that A = L*U. A can be（4）JOR雅可比超松弛迭代法求线性方程组Ax=b的解function x,n=JOR(A,b,x0,w,eps,M) if nargin=4 eps= 1.0e-6; M = 10000; elseif nargin =5 M = 10000; end if(w=2) %收敛条件要求error; return; end D=diag(diag(A); %求A的对角矩阵 B=w*inv(D); %迭代过程 x=x0; n=0; %迭代次数 tol=1; %迭代过程 while tol=eps x=x0-B*(A*x0-b); n = n+1; tol = norm(x-x0); x0 = x; if(n=M) disp(Warning: 迭代次数太多，可能不收敛！); return; end end （5）twostep两步迭代法求线性方程组Ax=b的解 function x,n=twostep(A,b,x0,eps,varargin) if nargin=3 eps= 1.0e-6; M = 200; elseif nargin=eps x0 =x; x12=B1*x0+f1; x =B2*x12+f2; n=n+1; if(n=M) disp(Warning: 迭代次数太多，可能不收敛！); return; end end （6）fastdown最速下降法求线性方程组Ax=b的解 function x,n= fastdown(A,b,x0,eps) if(nargin = 3) eps = 1.0e-6; end x=x0; n=0; tol=1; while(toleps) % （7）conjgrad共轭梯度法求线性方程组Ax=b的解 function x,n= conjgrad (A,b,x0) r1 = b-A*x0; p = r1; n = 0; for i=1:rank(A) %以下过程可参考算法流程 if(dot(p,A*p) 1.0e-50) %循环结束条件 break; end alpha = dot(r1,r1)/dot(p,A*p); x = x0+ alpha*p; r2 = r1- alpha*A*p; if(r2 1.0e-50) %循环结束条件 break; end belta = dot(r2,r2)/dot(r1,r1); 4结果高斯的计算步骤很多要进行一系列矩阵的计算，而且必须满足矩阵A的行、列相等才有唯一解，不然就会有不同的结果，而且计算时间很长，速度很慢，而追赶法计算比较快，不用进行矩阵的计算就可以得出结果。参考文献（1） 周晓阳。数学实现与Matlab，武汉：华中科技大学出版，2002（2） 王锁林，高杰，孙信得，Matlab混合编程与工程。北京：清华大学出版社，2008,5:120（3） 黄明游，刘帆，刘涛数值计算方法。北京：科学出