数值线性代数课程设计报告

燕山大学课程设计题目：Gauss消去法及其变形方法求解线性方程组学院〔系〕：年级专业：学号：学生姓名：指导教师：教师职称：燕山大学课程设计〔论文〕任务书院〔系〕：理学院基层教学单位：信息与计算科学系学号学生姓名专业〔班级〕设计题目Gauss消去法及其变形方法求解线性方程组设计技术参数Gauss消去法是一个古老的求解线性方程组的方法，但由于它的改良、变形得到的平方根法和改良的平方根法仍是在计算机上求解系数矩阵为中、低阶稠密矩阵的线性方程组常用的有效方法，本文对这几种方法进行了比拟。

用高斯消去法求解n阶线性方程组Ax=b的根本思想是在逐步消元的过程中，把方程组的系数矩阵化为上三角矩阵，从而将原方程组约化为容易求解的等价三角方程组，然后进行回代求解。

设计要求通过查阅资料，了解课题在实际中的应用；运用所学知识，建立问题的数学模型；利用MATLAB软件进行计算，并给出计算机程序；运用所学知识，对计算结果进行分析；要阅读不少于四篇的相关文献。工作量本课程设计要在一周内完成；报告的正文在4000字左右。工作计划第1-2天:阅读参考文献,收集资料.第3-5天:构建数学模型,编制计算机程序.第5-6天:上机运算,调试程序,计算结果分析.第7天:撰写报告.参考资料[1]张德丰.MATLAB数值分析与应用[M].北京：国防工业出版社.2007.1[2]李庆扬，王能超，易大义.数值分析[M].武汉：华中科技大学出版社.2006.7[3]萧树铁.数学实验[M].北京：高等教育出版社.1990[4]莱〔LayD.C.〕，刘深泉.数值线性代数及其应用[M].机械工业出版社.2005.8指导教师签字基层教学单位主任签字说明：此表一式四份，学生、指导教师、基层教学单位、系部各一份。年月日燕山大学课程设计评审意见表指导教师评语：成绩：指导教师：年月日辩论小组评语：成绩：评阅人：年月日课程设计总成绩：辩论小组成员签字：年月日目录摘要………………2一、问题的提出………………2二、数学模型的建立………3三、计算机程序………………3四、计算结果及分析………5五、参考文献…………………20六、心得体会…………………20

摘要： Gauss消去法是求解线性方程组的最根本的方法之一。实际计算中，线性方程组的系数矩阵常常具有对称正定性，即其各阶顺序主子式及全部特征值均大于零。矩阵的这一特性使它的三角分解也有更简单的形式，从而导出一些特殊的解法，如平方根法和改良的平方根法。关键词：线性方程组，Gauss消去法，平方根法，改良的平方根法，Hilbert矩阵。一、问题的提出1、先用你所熟悉的计算机语言将不选主元和列主元Gauss消去法编写成通用程序，然后用你编写的程序求解下面的84阶方程组最后，将你的计算结果与方程组的精确解进行比拟，并就此谈谈你对Gauss消去法的看法。2、先用你所熟悉的计算机语言将平方根法和改良的平方根法编写成通用程序，然后用你编写的程序求解对称正定方程组，其中〔1〕随机地选取，系数矩阵为100阶矩阵〔2〕系数矩阵为40阶Hilbert矩阵，即系数矩阵的第行第列元素为，向量的第个分量为。3、用1中的程序求解2中的两个方程组并比拟所有的计算结果，然后评论各个方法的优劣。。二、数学模型的建立１、Gauss消去法：用消去法解方程组的根本思想是用逐次消去未知数的方法把原来方程组Ax=b化为与其等价的三角方程组，而求解三角方程组就容易了。换句话说，上述过程就是用行的初等变换将原方程组系数矩阵化为简单形式，从而将求解原方程组的问题转化为求解简单方程组的问题。２、平方根法：平方根法就是利用对称正定矩阵的三角分解而得到的求解对称正定方程组的一种有效方法。平方根法递推公式可以证明对于对称正定矩阵A,可以唯一地分解成A=LLT,其中L是非奇异下三角形矩阵。３、改良的平方根法：平方根算法的计算量约为普通三角分解算法的一半，但由于要用到开平方，效率不是很高，所以，便有了为效率而存在的改良版平方根算法。对Ａ进行ＬＵ分解，即A=LU，Ａ对称正定，可得ｕｉｉ＞０，那么可得A=LDLT分解计算公式，Ｌ是对角元素为１的单位下三角矩阵，Ｄ为对角矩阵。三、计算机程序模型一：利用Gauss消去法对线性方程组Ax=b进行求解。用MATLAB建立m文件DelGauss.m,程序如下：functionx=DelGauss(a,b)[n,m]=size(a);nb=length(b);det=1;x=zeros(n,1);fork=1:n-1fori=k+1:nifa(k,k)==0returnendm=a(i,k)/a(k,k);forj=k+1:na(i,j)=a(i,j)-m*a(k,j);endb(i)=b(i)-m*b(k);enddet=det*a(k,k);enddet=det*a(n,n);fork=n:-1:1forj=k+1:nb(k)=b(k)-a(k,j)*x(j);endx(k)=b(k)/a(k,k);end模型二：利用平方根法对线性方程组Ax=b进行求解。用MATLAB建立m文件pingfg.m,程序如下：function[x]=pingfg(A,b)%乔累斯基分解[n,n]=size(A);L=zeros(n,n);%实际上不用为L申请空间，使用A即可L(1,1)=sqrt(A(1,1));fork=2:nL(k,1)=A(k,1)/L(1,1);endfork=2:n-1L(k,k)=sqrt(A(k,k)-sum(L(k,1:k-1).^2));fori=k+1:nL(i,k)=(A(i,k)-sum(L(i,1:k-1).*L(k,1:k-1)))/L(k,k);endendL(n,n)=sqrt(A(n,n)-sum(L(n,1:n-1).^2));%解下三角方程组Ly=by=zeros(n,1);fork=1:nj=1:k-1;y(k)=(b(k)-L(k,j)*y(j))/L(k,k);end%解上三角方程组L'x=yx=zeros(n,1);U=L';fork=n:-1:1j=k+1:n;x(k)=(y(k)-U(k,j)*x(j))/U(k,k);end模型三：利用改良的平方根法对线性方程组Ax=b进行求解。用MATLAB建立m文件ave.m,程序如下：function[x]=ave(A,b,n)%用改良平方根法求解Ax=bL=zeros(n,n);%L为n*n矩阵D=diag(n,0);%D为n*n的主对角矩阵S=L*D;fori=1:n%L的主对角元素均为1L(i,i)=1;endfori=1:nforj=1:n%验证A是否为对称正定矩阵if(eig(A)<=0)|(A(i,j)~=A(j,i))%A的特征值小于0或A非对称时，输出wrongdisp('wrong');break;endendendHilbert矩阵用MATLAB建立m文件Hil.m,程序如下：functionb=Hil()fork=1:40form=1:40s=0;t=s+1/(k+m-1);s=t;endb(k,1)=s;end四、计算结果及分析问题1：>>A=ones(1,83);>>B=6*ones(1,84);>>C=8*ones(1,83);>>D=diag(A,1)+diag(B,0)+diag(C,-1);>>b=[7;15*ones(82,1);14];>>DelGauss(D,b)ans=1.0e+008*0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.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方程组的精确解为x1=x2=…=x84=1.0000，与Gauss消去法求得的解差距很大，所得结果不够准确，计算简单但其消元过程有时不能进行到底而使求解出现解失真的情况。问题2：>>A=ones(1,99);B=10*ones(1,100);C=diag(A,1)+diag(B,0)+diag(A,-1);b=ones(100,1);pingfg(C,b)ans=0.09180.08250.08340.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08340.08250.0918>>A=ones(1,99);B=10*ones(1,100);C=diag(A,1)+diag(B,0)+diag(A,-1);b=ones(100,1);>>ave(C,b,100)ans=0.09180.08250.08340.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08340.08250.0918A=hilb(40);b=Hil();pingfg(A,b)ans=1.0e+007*0.0000-0.00000.0001-0.00190.0171-0.09390.3305-0.76851.1892-1.20770.7761-0.30260.06840.0135-0.0220-0.04500.0954-0.06750.0217-0.0032-0.00060.0013-0.00090.00020.0009-0.00150.0009-0.00060.0005-0.00020.0001-0.00010.00000.0000-0.0000-0.00000.0000-0.0000-0.00000.0000>>A=hilb(40);b=Hil();ave(A,b,40)ans=-0.00000.0000-0.00010.0010-0.00880.0459-0.14160.2546-0.25990.2041-0.28120.2680-0.10590.5017-0.7878-0.13990.36920.5230-0.23950.0245-0.4819-0.20150.6467-0.2634-0.02930.4175-0.43350.4124-0.4021-0.14970.15800.5615-0.2600-0.72070.34570.09950.27410.2579-0.79171.3333问题3：>>A=ones(1,99);B=10*ones(1,100);C=diag(A,1)+diag(B,0)+diag(A,-1);b=ones(100,1);>>DelGauss(C,b)ans=0.09180.08250.08340.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.08330.0833