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

- 1- 解线性方程组的直接法的Matlab 实现 姓名 ********** 摘要:给出用MATLAB 解线性方程组的各种方法,用MATLAB 直接操作,不用编程,便可立即求出线性方程组的解.方法直观、简便、速度快,具有较强的实用性,另外提供了Jacobi 迭代法程序. 1 引言 线性方程组Ax=b 是我们在科学和工程计算中经常出现的数学模型,大量的科技与工程实际问题,常常归结为解线性代数方程组,有关线性方程组解的存在性和唯一性在“线性代数”理论中已经作过详细介绍,本章的主要任务是讨论系数行列式不为零的n 阶非齐次线性方程组Ax=b 的两类主要求方法:直接法(精确法)和迭代法。对它的解法我们最熟悉的就是主元消去法,但它只是适用于A 是低价稠密的矩阵,对于由工程技术中产生的大型稀疏矩阵方程组(即A 的阶数n 很大,但零元素较多,例如求某些偏微分方程数值解所产生的线性方程组,n ≥104),还需利用迭代法求解。如在计算机内存和运算两方面,都可以根据A 中有大量零元素的特点采用迭代法。本文将介绍两种常见的迭代:Jacobi 迭代法和Gauss-Seidel 迭代,并用迭代法在数学软件Matlab 上实现线性方程组的解。1迭代法的基本思想迭代法是按照某种规则构造一个向量序列{x(k)},使其极限向量x*是Ax=b 的精确解。因此,对迭代法来说一般有下面几个问题:(1)如何构造迭代序列?(2)构造的迭代法序列是否收敛?在什么情况下收敛?(3)如果收敛,收敛的速度如何?我们应该给予量的刻划,用以比较各种迭代法收敛的快慢。 2 相关知识 线性方程组的概念及分类线性方程组的一般形式为a11x1+a12x2+…+a1nxn=b1a21x1+a22x2+…+a2nxn=b2am1x1+am2x2+…+amnxn=b{n(1)若记X=x1x2(…x n)T,b=b1 b2(…bn)TA=a11 a12…a1na21 a22…a2n am1 am2…a mn 则线性方程组(1)记为AX=b.(2)若b 的元素不全为零,则称方程组 (1)为非齐次线性方程组;若b 的元素全为零,即b1=b2=…=bn=0,则AX=0.(3) 并称方程组(3)为齐次线性方程组,也称作方程组(2)的导出方程组,称(A b)=a11 a12…a1n …b1a21 a22…a2n …b2 am1 am2…amn …b 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 zhica>0, 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’) 变换后的方程组为 A‘x=b’ 与原方程组等价(同解) 因此,高斯消去法实际上就是通过上述三种变换,把(A,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 matrix N = length(b); % Perform Gaussian Elimination if A(1,1)==0 n_temp=find(A(:,1)~=0); c=A(n_temp(1),:); A(n_temp(1),:)=A(1,:); A(1,:)=c; end for 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 substitution if A(N,N)==0 error('the coefficient matrix of the linear equations is sigular'); end x = 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 % % % % -------------------------------------------- 在消元法中,舍入误差的增长是十分重要的。如果abs(A(k,k))很小,在运算中用它作为分母会导致舍入误差的扩散。 如果调换方程组的次序,使得运算中做分母的绝对值尽量大,以减少舍入误差的影响。 如果在高斯消元法的每一步中,选取绝对值最大的元素,将其调到主对角线位置再做消元,则称为列主元消元法。 例如第k步,则将第k列绝对值最大的元素调到A(k,k)的位置,再消元。 列主元消元法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 matrix N = 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 Elimination for 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 substitution if A(N,N)==0 error('the coefficient matrix of the linear equations is sigular'); end x = 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消元法。 由于将上三角阵化为对角矩阵的工作量比直接从上三角阵回待求解的工作量。因此,用回待求解方程组比Gauss-Jordan继续消元更适用。故本文不对Gauss-Jordan消元做进一步考虑。 (3)直接分解法 (4)高斯消元法对系数矩阵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 form where 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 form where 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 lower triangular and permutation matrices) in L, so that A = L*U. A can be rectangular. 注意:用Matlab中的[L,U]=lu(A)命令时,分解出来的U是上对角矩阵没错,可是L却不是下对角矩阵。 其原因在于: 定理:当一个方阵A的1~(n-1)阶顺序主子式都不为0时,A存在唯一的doolittle分解。 对于满足条件的这个方阵,采用直接三角分解法可以得到这个唯一的L、U矩阵。但事实上,m atlab 的lu函数并不是采用直接三角分解法,因为直接三角分解法存在着引入很大舍入误差的可能(直接分解法过程中,当对角元为零或很小时)。 matlab的lu函数采用的是列主元三角分解法,又称为PLU分解,实现LU=PA。该分解法要求的前提条件是A为非奇异矩阵。 4结果 高斯的计算步骤很多要进行一系列矩阵的计算,而且必须满足矩阵A的行、列相等才有唯一解,不然就会有不同的结果,而且计算时间很长,速度很慢,而追赶法计算比较快,不用进行矩阵的计算就可以得出结果。 致谢 在完成之际,我要特别感谢我的指导老师***老师的热情关怀和悉心指导.本论文在xxx导师的悉心指导下完成的。导师渊博的专业知识、严谨的治学态度,精益求精的工作作风,诲人不倦的高尚师德,严于律己、宽以待人的崇高风范,朴实无法、平易近人的人格魅力对本人影响深远。不仅使本人树立了远大的学习目标、掌握了基本的研究方法,还使本人明白了许多为人处事的道理。本次论文从选题到完成,每一步都是在导师的悉心指导下完成的,倾注了导师大量的心血。在此,谨向导师表示崇高的敬意和衷心的感谢!在写论文的过程中,遇到了很多的问题,在老师的耐心指导下,问题都得以解决。所以在此,再次对老师道一声:老师,谢谢您! 时光匆匆如流水,转眼便是大学毕业时节,春梦秋云,聚散真容易。离校日期已日趋渐进,毕业论文的完成也随之进入了尾声。从开始进入课题到论文的顺利完成,一直都离不开老师、同学、朋友给我热情的帮助,在这里请接受我诚挚的谢意!在此我向xxxx学校xxxxxxx专业的所有老师表示衷心的感谢,谢谢你们三年的辛勤栽培,谢谢你们在教学的同时更多的是传授我们做人的道理,谢谢三年里面你们孜孜不倦的教诲! 声明人(签名): 指导教师(签名): 年月日