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

1、计算方法实验报告 1【课题名称】用列主元高斯消去法和列主元三角分解法解线性方程【目的和意义】高斯消去法是一个古老的求解线性方程组的方法，但由它改进得到的选主元的高斯消去法则是目前计算机上常用的解低阶稠密矩阵方程组的有效方法。用高斯消去法解线性方程组的基本思想时用矩阵行的初等变换将系数矩阵A 约化为具有简单形式的矩阵(上三角矩阵、单位矩阵等)，而三角形方程组则可以直接回带求解用高斯消去法解线性方程组Ax b(其中A Rn xn)的计算量为：乘除法运算步骤为(n 1) n (2 n 5)6。相比之下，传统的克莱姆19法则则较为繁琐，如求解 20 阶线性方程组，克莱姆法则大约要5 10次乘法，而用高

2、斯消去法只需要 3060 次乘除法。在高斯消去法运算的过程中，如果出现abs(A(i,i)等于零或过小的情况，则会导致矩阵元素数量级严重增长和舍入误差的扩散，使得最后的计算结果不可靠，所以目前计算机上常用的解低阶稠密矩阵方程的快速有效的方法时列主元高斯消去法，从而使计算结果更加精确。2、列主元三角分解法高斯消去法的消去过程，实质上是将 A 分解为两个三角矩阵的乘积 A=LU，并求解 Ly=b 的过程。回带过程就是求解上三角方程组 Ux=y。所以在实际的运算中， 矩阵 L 和 U 可以直 接计算出，而不需要任何中间步骤， 从而在计算过程中将高斯消去法的步骤进行了进一步的n (n 1)(n 1)

3、n (2 n1)n (n 1)2 62(n 1) n (2 n 1 )n ( n1)n (n 1)6 22n3，加减运算步骤为MDASn ( n 1)2简略，大大提高了运算速度，这就是三角分解法采用选主元的方式与列主元高斯消去法一样，也是为了避免除数过小，从而保证了计算的精确度【计算公式】1、列主元高斯消去法设有线性方程组 Ax=b，其中设 A 为非奇异矩阵。方程组的增广矩阵为a11ai2LLa1nMbia21a22LLa2nMb2MMa,bMMai1MMM Man1an2LLannMbn第 1 步(k=1)：首先在 A 的第一列中选取绝对值最大的元素ai1，作为第一步的主元素al1max a

4、i101 i n然后交换(A，b)的第 1 行与第 I 行元素，再进行消元计算。设列主元素消去法已经完成第1 步到第 k-1 步的按列选主元，交换两行，消元计算得到与原方程组等价的方程组A(k)x=b(k)aa11a12La1ka(1)1nM bL(2)Ma22a2ka2nb2OMMMA,b A(k),b(k)akk)Lakn)M bkk)(k)ak 1,k(k)ak 1,nM(k)bk1MMM Maii第 k 步计算如下:对于 k=1 , 2，n-1Jk)atk(1) 按列选主元：即确定 t 使(2) 如果 t 殊，则交换A , b第 t 行与第(3 )消元计算aikaikmik,(i k

5、1,L , n) akkaijajmkakj,(i,jk 1L ,n)bbimkbk,(ik 1L ,n)消元乘数mik 满足:(4 )回代求解XnbnannXi(bij iaijXj)1-,(i n 1,n 2,-,1)max aik0k i nk 行元素。2、列主元三角分解法A A,bk 1UkjakjIkmUmj(j k,k 1,L , n;k 1,2,L , n)m 1k 1lik(aiklimUmk)/Ukk(i k 1,k 2 丄,n;k 1,2,L ,n)m 1k,k 1,L , n)，于是有Ukk=Sk。如果sunU12LU1,k 1U1kLu1jLU1ny1.1U22LU2,

6、k 1U2kLU2jLU2ny2MMMMMMML1,11k 1,2LUk 1,k 1Uk 1,kLUk 1,jLUk 1,nyk 1A気lk2Llk,k 1akkLakjLaknbkMMMMMMMli1li2Lli ,k 1aikLdLanbMMMMMMMIn1ln2Lln,k 13nkL3njLannbn对方程组的增广矩阵经过 k-1 步分解后，可变成如下形式:第 k 步分解，为了避免用绝对值很小的数Ukk作除数，引k 1limumkm 1maxS，则将矩阵的第t 行与第 k 行元素互换，将（i，j）位置的新元素仍记为l a -jj或jj，然后再做第 k 步分解，这UkkSk(Sk即交换前的

7、 St),likS /Sk(i k 1,k 2,L ,n)且 lik1 (i k 1,k2,L , n),【列主元高斯消去法程序流程图】辅入方程：系数矩陪九，力以及主兀九许昂、小值 c刘断輸丸的芳严程組是舌肖尿fc-1,25.n*l按别透期工云艮宜 T 标p, q=rti 3K 匚曲减赳匕:切）祁将巾的行号粧換苗左丸中的行岂，訊 q q+k- lBon (p) cAO.A心:工Atenipl;tEtnp2=Xk-；H【列主元高斯消去法Matlab 主程序】function x=gauss1(A,b,c)%列主元法高斯消去法解线性方程Ax=bif (le ngth(A)=le ngth(b)%判

8、断输入的方程组是否有误disp( 输 入 方 程 有 误 ！)return;enddisp(原方程为 AX=b :)Abdisp(- )%找列主元%找出第 k 列中的最大值，其下标为p,q%q 在 A(k:n,k)中的行号转换为在A 中的行if abs(p)ktemp 仁A(k,:);%列主元所在行不是当前行，将当前行与列%显示方程组n=length(A);for k=1: n-1p,q=max(abs(A(k: n,k); q=q+k_1;号A(k,:)=A(q,:);元所在行交换(包括 b)A(q,:)=temp1;temp2=b(k,:);b(k,:)=b(q,:);b(q,:)=tem

9、p2;endfor i=k+1: nm(i,k)=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(消元后所得到的上三角阵是)Ab(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;%消元%A(k,k) 将 A(i,k)消为 0 所乘系数%第 i 行消元处理%b 消元处理%显示消元后的系数矩阵2disp（AX=b 的解 x 是）x=b;【调

10、用函数解题】Command Window A=0 3 4：1 -1 1：2 1 21 b=l ;2;3; c=0.00301; 其=匚厘u孚$1 CAj c）原方程为:A =0341-I1212消元后師得到的上三角阵罡21203400 2AS：=b的鯛龙是1 16盯-0- 33330.5000【列主元三角分解法程序流程图】【列主元三角分解法Matlab 主程序】自己编的程序：fun ction x=PLU(A,b,eps)数if (le ngth(A)=length(b)%判断输入的方程组是否有误%定义函数列主元三角分解法函disp（输入方程有误！）return;enddisp(原方程为 AX

11、=b :)Abdisp(- )n=len gth(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,:)=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);elseu(r,r)=A(r,r)-A(r,1:r-1)*A(1:r-1,r);

12、if abs(u(r,r)v=epsP=A(r,:);A(r,:)=A(r+1,:);A(r+1,:)=p;elseend for 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,i)%如果 u (r,r)小于 e,则交换行%根据公式求解，并把结果存在矩阵 A%根据公式求解，并把结果存在矩阵endendendy(1)=A(1, n+1);for i=2: nh=0;for k=1:i-1h=

13、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)可直接得到 P,L,U 并解出方程解的的程序(查阅资料得子函数 PLU1，其作用 是将矩阵 A 分解成 L 乘以 U 的形式。PLU2 为调用 PLU1 解题的程序，是自己 编的)(I).fun ctio n l,u,p=PLU1(A

14、)%定义子函数，其功能为列主元三角分解系数矩阵 Am, n=size(A);%判断系数矩阵是否为方阵if m=nerror(矩阵不是方阵)returnendif det(A)=0%判断系数矩阵能否被三角分解error(矩阵不能被三角分解)endu=A;p=eye(m);l=eye(m);%将系数矩阵三角分解，分别求出PL,Ufor i=1:mfor j=i:mt(j)=u(j,i);for k=1:i-1t(j)=t(j)-u(j,k)*u(k,i);endenda=i;b=abs(t(i);for j=i+1:mif b x-PLU2 (Aj b)1u -2120340022;1 =K COOO000 1.0000Q0 5000-0. 5000K 0000p =b二312AX匸b的解*是L 1557一0.33330. 5000【编程疑难】这是第一