列主元高斯消去法和列主元三角分解法解线性方程

1、实用文案计算方法实验报告 1【课题名称】用列主元高斯消去法和列主元三角分解法解线性方程【目的和意义】高斯消去法是一个古老的求解线性方程组的方法， 但由它改进得到的选主元的高斯消去法则是目前计算机上常用的解低阶稠密矩阵方程组的有效方法。用高斯消去法解线性方程组的基本思想时用矩阵行的初等变换将系数矩阵 A 约化为具 有简单形式的矩阵(上三角矩阵、单位矩阵等) ，而三角形方程组则可以直接回带求解n (n 1) (n 1) n (2 n1)n (n 1)262(n 1)n(2 n 1) n( n1)n (n 1)622n3n2n3 ，加减运算步骤为MDASn (n 1)2用高斯消去法解线性方程组 Ax

2、 b (其中 ARnn)的计算量为：乘除法运算步骤为(n 1)n (2 n 5)6 。相比之下，传统的克莱姆19法则则较为繁琐，如求解 20 阶线性方程组，克莱姆法则大约要 5 10 次乘法，而用高斯 消去法只需要 3060 次乘除法。在高斯消去法运算的过程中，如果出现 abs(A(i,i) 等于零或过小的情况，则会导致矩阵 元素数量级严重增长和舍入误差的扩散， 使得最后的计算结果不可靠， 所以目前计算机上常 用的解低阶稠密矩阵方程的快速有效的方法时列主元高斯消去法， 从而使计算结果更加精确。 2、列主元三角分解法高斯消去法的消去过程， 实质上是将 A 分解为两个三角矩阵的乘积 A=LU ，并

3、求解 Ly=b 的过程。 回带过程就是求解上三角方程组 Ux=y 。所以在实际的运算中， 矩阵 L 和 U 可以直 接计算出， 而不需要任何中间步骤， 从而在计算过程中将高斯消去法的步骤进行了进一步的标准文档实用文案简略，大大提高了运算速度，这就是三角分解法采用选主元的方式与列主元高斯消去法一样， 也是为了避免除数过小， 从而保证了计算的精确度【计算公式】1、 列主元高斯消去法设有线性方程组 Ax=b ，其中设 A 为非奇异矩阵。方程组的增广矩阵为a11a12LLa1nMb1a21a22LLa2nMb2MMa,bMMal1MMMMan1an2LLannMbn第 1 步( k=1 )：首先在 A

4、 的第一列中选取绝对值最大的元素 al1 ，作为第一步的主元素al1 max ai1 01 i n然后交换( A，b)的第 1 行与第 l 行元素，再进行消元计算。设列主元素消去法已经完成第 1 步到第 k-1 步的按列选主元，交换两行，消元计算得到与原方程组等价的方程组 A(k)x=b(k)(1) (1) a11a12La(1)1ka(1)1nMb1(1)(2)L(2)(2)M(2)a22a2ka2nb2OMMM标准文档 (k) (k)A,bA(k),b(k) ak(kk) Lak(kn)Mbk(k)(k) ak 1,k(k) ak 1,nM(k) bk 1MMMM实用文案第k 步计算如下：

5、对于 k=1 ，2 ，n-1(k)(k)atkmax kinaik01 ）按列选主元：即确定 t 使2）如果 t k ，则交换 A，b第t 行与第 k 行元素。3 ）消元计算aikaikmik,(i k 1,L ,n) akkaijaijmik akj ,(i, jk 1,L ,n)bibi mikbk ,(ik 1,L ,n)消元乘数4 ）回代求解xnnbnann(biaijxj )ji1xi 标准文档 ij i 1 ,(i n 1,n 2, ,1)aii实用文案2、 列主元三角分解法 A A,b对方程组的增广矩阵 经过 k-1 步分解后，可变成如下形式：第 k 步分解，为了避免用绝对值很小

6、的数ukk作除数，引进量u11u12Lu1,k 1u1kLu1jLu1ny1l21u22Lu2,k 1u2kLu2jLu2ny2MMMMMMMlk 1,1lk 1,2Luk 1,k 1uk 1,kLuk 1, jLuk 1,nyk 1Alk1lk2Llk,k 1akkLakjLaknbkMMMMMMMli1li2Lli ,k 1aikLaijLainbiMMMMMMMln1ln2Lln,k 1ankLanjLannbnk1ukj akjlkmumj( j k,k 1,L ,n;k 1,2,L ,n)m1k1lik (aiklimumk) /ukk (i k 1,k 2,L ,n;k 1,2,L

7、 ,n)m1siaikk1limumkm1k,k 1,L ,n) u s，于是有 ukk = sk 。如果stmk iaxn si ，则将矩阵的第t 行与第 k 行元素互换，将（i，j）位置的新元素仍记为lajj 或 jj ，然后再做第 k 步分解，这标准文档实用文案ukk sk(sk即交换前的 st )，lik si / sk( i k 1,k 2,L ,n)且 lik 1 (i k 1,k 2,L ,n),列主元高斯消去法程序流程图】列主元高斯消去法 Matlab 主程序】function x=gauss1(A,b,c)% 列主元法高斯消去法解线性方程 Ax=b标准文档实用文案if (le

8、ngth(A)=length(b)disp( 输入方程有误！ )return;enddisp( 原方程为 AX=b ：)Abdisp()n=length(A);for k=1:n-1p,q=max(abs(A(k:n,k);q=q+k-1;号if abs(p)k标准文档实用文案元所在行交换 ( 包括 b)A(k,:)=A(q,:);A(q,:)=temp1;temp2=b(k,:);b(k,:)=b(q,:);b(q,:)=temp2; end%消元%A(k,k) 将 A(i,k) 消为 0 所乘系数% 第 i 行消元处理%b 消元处理%显示消元后的系数矩阵for i=k+1:n m(i,k)

9、=A(i,k)/A(k,k);A(i,k:n)=A(i,k:n)-m(i,k)*A(k,k:n); b(i)=b(i)-m(i,k)*b(k);endenddisp( 消元后所得到的上三角阵是 )A b(n)=b(n)/A(n,n); % 回代求解for i=n-1:-1:1b(i)=(b(i)-sum(A(i,i+1:n)*b(i+1:n)/A(i,i);end clear x;标准文档实用文案disp（AX=b 的解 x 是 ）x=b;调用函数解题】标准文档实用文案列主元三角分解法程序流程图】标准文档实用文案【 列主元三角分解法 Matlab 主程序】%定义函数列主元三角分解法函自己编的程

10、序：function x=PLU(A,b,eps)%判断输入的方程组是否有误数if (length(A)=length(b)标准文档实用文案disp( 输入方程有误！ )return;enddisp( 原方程为 AX=b ：)Abdisp()n=length(A);A=A b;for r=1:nif r=1for i=1:nc d=max(abs(A(:,1);if c=eps结束break;elseendd=d+1-1;p=A(1,:);A(1,:)=A(d,:);%显示方程组%将 A 与 b 合并，得到增广矩阵%选取最大列向量，并做行交换% 最大值小于 e，主元太小，程序标准文档A(d,:

11、)=p;实用文案A(1,i)=A(1,i);endA(1,2:n)=A(1,2:n);A(2:n,1)=A(2:n,1)/A(1,1);else u(r,r)=A(r,r)-A(r,1:r-1)*A(1:r-1,r);if abs(u(r,r)=eps p=A(r,:); A(r,:)=A(r+1,:); A(r+1,:)=p;elseendfor i=r:nA(r,i)=A(r,i)-A(r,1:r-1)*A(1:r-1,i);中endfor i=r+1:nA(i,r)=(A(i,r)-A(i,1:r-1)*A(1:r-1,r)/A(r,r);A中% 求 u(1,i)%按照方程求取 u(r,

12、i)%如果 u (r,r)小于 e，则交换行%根据公式求解，并把结果存在矩阵 A% 根据公式求解， 并把结果存在矩阵endendend标准文档实用文案y(1)=A(1,n+1);for i=2:nh=0;for k=1:i-1h=h+A(i,k)*y(k);endy(i)=A(i,n+1)-h;endx(n)=y(n)/A(n,n);for i=n-1:-1:1h=0;for k=i+1:nh=h+A(i,k)*x(k);endx(i)=(y(i)-h)/A(i,i);endAdisp(AX=b 的解 x 是)x=x;%根据公式求解 y(i)% 根据公式求解 x(i)%输出方程的解标准文档实用

13、文案可直接得到 P,L,U 并解出方程解的的程序(查阅资料得子函数 PLU1 ，其作用 是将矩阵 A分解成 L乘以 U的形式。PLU2为调用 PLU1 解题的程序，是自己 编的)().function l,u,p=PLU1(A)%定义子函数，其功能为列主元三角分解系数矩阵Am,n=size(A);%判断系数矩阵是否为方阵if m=nerror( 矩阵不是方阵 )returnendif det(A)=0% 判断系数矩阵能否被三角分解error( 矩阵不能被三角分解 ) endu=A;p=eye(m);l=eye(m);% 将系数矩阵三角分解，分别求出P,L,Ufor i=1:mfor j=i:m

14、t(j)=u(j,i);for k=1:i-1标准文档实用文案t(j)=t(j)-u(j,k)*u(k,i);endenda=i;b=abs(t(i);for j=i+1:mif babs(t(j) b=abs(t(j); a=j;endendif a=ifor j=1:m c=u(i,j); u(i,j)=u(a,j); u(a,j)=c;endfor j=1:m c=p(i,j); p(i,j)=p(a,j); p(a,j)=c;endc=t(a);标准文档实用文案t(a)=t(i);t(i)=c;endu(i,i)=t(i);for j=i+1:mu(j,i)=t(j)/t(i);end

15、 for j=i+1:mfor k=1:i-1u(i,j)=u(i,j)-u(i,k)*u(k,j);endendendl=tril(u,-1)+eye(m);u=triu(u,0)%定义列主元三角分解法的().function x=PLU2(A,b)% 调用 PLU 分解系数矩函数l,u,p=PLU1(A)阵A标准文档实用文案m=length(A);%由于 A 左乘 p ，故 b 也要左乘pv=b;for q=1:mb(q)=sum(p(q,1:m)*v(1:m,1);endb(1)=b(1)%求解方程Ly=bfor i=2:1:mb(i)=(b(i)-sum(l(i,1:i-1)*b(1:i-1);%求解方程end b(m)=b(m)/u(m,m);Ux=yfor i=m-1:-1:1b(i)=(b(i)-sum(u(i,i+1:m)*b(i+1:m)/u(i,i);endclear x;disp(AX=b 的解 x 是)标准文档实用文案x=b;【调用函数解题】标准文档实用文案编程疑难】标准文