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

1、解线性方程组的直接法的matlab实现姓名*摘要：给出用matlab解线性方程组的各种方法,用matlab直接操作,不用编程,便可立即求出线性方程组的解.方法直观、简便、速度快,具有较强的实用性,另外提供了jacobi迭代法程序.关键字：线性方程组 数值解 程序设计 matlab jacobi迭代法 数据结构1 引言线性方程组ax=b是我们在科学和工程计算中经常出现的数学模型,大量的科技与工程实际问题，常常归结为解线性代数方程组，有关线性方程组解的存在性和唯一性在“线性代数”理论中已经作过详细介绍，本章的主要任务是讨论系数行列式不为零的n阶非齐次线性方程组ax=b的两类主要求方法：直接法（精确

2、法）和迭代法。对它的解法我们最熟悉的就是主元消去法,但它只是适用于a是低价稠密的矩阵,对于由工程技术中产生的大型稀疏矩阵方程组(即a的阶数n很大,但零元素较多,例如求某些偏微分方程数值解所产生的线性方程组,n104),还需利用迭代法求解。如在计算机内存和运算两方面,都可以根据a中有大量零元素的特点采用迭代法。本文将介绍两种常见的迭代:jacobi迭代法和gauss-seidel迭代,并用迭代法在数学软件matlab上实现线性方程组的解。1迭代法的基本思想迭代法是按照某种规则构造一个向量序列x(k),使其极限向量x*是ax=b的精确解。因此,对迭代法来说一般有下面几个问题:(1)如何构造迭代序列

3、?(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)并称方程

4、组(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(

5、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); fo

6、r 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与原方程组等价（同解）因此，高斯消去法实际上就是通过上述三种变换，把（a，b）化为等价的上三角形式。例：用高斯消去法找出下列方程組的解原方程组增广矩阵跟著以上的方法來運算，這個矩陣可以轉

7、變為以下的樣子：如果增广矩阵第一行第一列不为零，将第二，三行第一列化为零。然后第三行的第二列化为零，则其解便可求得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);% perfor

8、m 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)

9、,:); 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% % % % -在消元法中，舍入误差的增长是十分重要的。如果abs（a（k，k

10、）很小，在运算中用它作为分母会导致舍入误差的扩散。如果调换方程组的次序，使得运算中做分母的绝对值尽量大，以减少舍入误差的影响。如果在高斯消元法的每一步中，选取绝对值最大的元素，将其调到主对角线位置再做消元，则称为列主元消元法。例如第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

11、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

12、 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

13、 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消元法：用初等变换化系数矩阵为对角矩阵的方法称为gau

14、ss-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

15、 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 s

16、ize), 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 th

17、e 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.注意：用matl

18、ab中的l,u=lu(a)命令时，分解出来的u是上对角矩阵没错，可是l却不是下对角矩阵。其原因在于：定理：当一个方阵a的1（n-1）阶顺序主子式都不为0时，a存在唯一的doolittle分解。对于满足条件的这个方阵，采用直接三角分解法可以得到这个唯一的l、u矩阵。但事实上，matlab的lu函数并不是采用直接三角分解法，因为直接三角分解法存在着引入很大舍入误差的可能（直接分解法过程中，当对角元为零或很小时）。matlab的lu函数采用的是列主元三角分解法，又称为plu分解，实现lu=pa。该分解法要求的前提条件是a为非奇异矩阵。4结果高斯的计算步骤很多要进行一系列矩阵的计算，而且必须满足矩阵a的行、列相等才有唯一解，不然就会有不同的结果，而且计算时间很长，速度很慢，而追赶法计算比较快，不用进行矩阵的计算就可以得出结果。致 谢在论文完成之际，我要特别感谢我的指导老师*老师的热情关怀和悉心指导本论文在xxx导师的悉心指导下完成的。导师渊博的专业知识、严谨的治学态度，精益求精的工作作风，诲人不倦的高尚师德，严于律己、宽以待人的崇高风范，朴实无法、平易近人的人格魅力对本人影响深远。不仅使本人树立了远大的学习目标、掌握了基本的研究方法，还使本人明白了许多为人处事的道理。本次论文从选题到完成，每一步都是在导师的悉心指导下完成的，倾注了导师大量的心