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

计算方法实验报告1【课题名称】用列主元高斯消去法和列主元三角分解法解线性方程【目的和意义】高斯消去法是一个古老的求解线性方程组的方法，但由它改进得到的选主元的高斯消去法则是目前计算机上常用的解低阶稠密矩阵方程组的有效方法。

用高斯消去法解线性方程组的基本思想时用矩阵行的初等变换将系数矩阵A 约化为具有简单形式的矩阵（上三角矩阵、单位矩阵等），而三角形方程组则可以直接回带求解 用高斯消去法解线性方程组b Ax =（其中A ∈Rn ×n ）的计算量为：乘除法运算步骤为32(1)(1)(21)(1)(1)262233n n n n n n n n n n nMD n ----+=+++=+-，加减运算步骤为(1)(21)(1)(1)(1)(25)6226n n n n n n n n n n AS -----+=++=。

相比之下，传统的克莱姆法则则较为繁琐，如求解20阶线性方程组，克莱姆法则大约要19510⨯次乘法，而用高斯消去法只需要3060次乘除法。

在高斯消去法运算的过程中，如果出现abs(A(i,i))等于零或过小的情况，则会导致矩阵元素数量级严重增长和舍入误差的扩散，使得最后的计算结果不可靠，所以目前计算机上常用的解低阶稠密矩阵方程的快速有效的方法时列主元高斯消去法，从而使计算结果更加精确。

2、列主元三角分解法高斯消去法的消去过程，实质上是将A 分解为两个三角矩阵的乘积A=LU ，并求解Ly=b 的过程。

回带过程就是求解上三角方程组Ux=y 。

所以在实际的运算中，矩阵L 和U 可以直接计算出，而不需要任何中间步骤，从而在计算过程中将高斯消去法的步骤进行了进一步的简略，大大提高了运算速度，这就是三角分解法 采用选主元的方式与列主元高斯消去法一样，也是为了避免除数过小，从而保证了计算的精确度【计算公式】1、 列主元高斯消去法设有线性方程组Ax=b ，其中设A 为非奇异矩阵。

方程组的增广矩阵为第1步（k=1）：首先在A 的第一列中选取绝对值最大的元素1l a ，作为第一步的主元素:111211212222112[,]n n n l n nn n a a a a b a a a b a a a b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦a b然后交换（A ，b ）的第1行与第l 行元素，再进行消元计算。

设列主元素消去法已经完成第1步到第k-1步的按列选主元，交换两行，消元计算得到与原方程组等价的方程组 A(k)x=b(k)第k 步计算如下：对于k=1，2，…，n-1（1）按列选主元：即确定t 使 （2）如果t ≠k ，则交换[A ，b]第t 行与第k 行元素。

（3）消元计算消元乘数mik 满足:（4）回代求解2、 列主元三角分解法 对方程组的增广矩阵 经过k-1步分解后，可变成如下形式：111max 0l i i n a a ≤≤=≠(1)(1)(1)(1)(1)1112111(2)(2)(2)(2)22222()(()1)()()()()()1,1()(,)()[,][,] k k k k nk k nk n k k k k k kk kn k k k k n k k k n nn a a a a b a a a b a a b a b b a a a +++⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥→=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦A b A b ()()max 0k k tk ik k i na a ≤≤=≠,(1,,)ik ik ik kka a m i k n a ←=-=+, (,1,,), (1,,)ij ij ik kji i ik k a a m a i j k n b b m b i k n ←+=+⎧⎨←+=+⎩⎪⎪⎩⎪⎪⎨⎧--=-←←∑+=)1,,2,1(,)(1n n i a x a b x a b x ii n i j j ij i i nnn n [,]A A b =11121,11111222,122221,11,1,1,211,11,2121,112,112,1k k k k k k k j n k k j n k k k i i i k n n kk kj kn k ik ij in i nknjk k k j k n n nnk k n a a a b A a u u u u u u y l l l l l l ll l l l u u u u u y u u u u y a a b a a b l a -------------⎡→⎣⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎦第k 步分解，为了避免用绝对值很小的数kku 作除数，引进量1111 (,1,,;1,2,,) ()/ (1,2,,;1,2,,)k kj kj km mj m k ik ik im mk kkm u a l u j k k n k n l a l u u i k k n k n -=-=⎧=-=+=⎪⎪⎨⎪=-=++=⎪⎩∑∑11(,1,,)k i ik im mk m s a l u i k k n -==-=+∑，于是有kk u =ks 。

如果 ，则将矩阵的第t 行与第k 行元素互换，将（i ，j ）位置的新元素仍记为jjl 或jja ，然后再做第k 步分解，这时【列主元高斯消去法程序流程图】max t i k i ns s ≤≤= ()/ 1,2,,)1 (1,2,,),kk k k t iki k ik u s s s l s s i k k n l i k k n ===++≤=++即交换前的，（且【列主元高斯消去法Matlab主程序】function x=gauss1(A,b,c) %列主元法高斯消去法解线性方程Ax=bif (length(A)~=length(b)) %判断输入的方程组是否有误disp('输入方程有误！')return;enddisp('原方程为AX=b：') %显示方程组Abdisp('------------------------')n=length(A);for k=1:n-1 %找列主元[p,q]=max(abs(A(k:n,k))); %找出第k列中的最大值，其下标为[p,q] q=q+k-1; %q在A(k:n,k)中的行号转换为在A中的行号if abs(p)<cdisp('列元素太小，det(A)≈0');break;elseif q>ktemp1=A(k,:); %列主元所在行不是当前行，将当前行与列主A(k,:)=A(q,:); 元所在行交换(包括b)A(q,:)=temp1;temp2=b(k,:);b(k,:)=b(q,:);b(q,:)=temp2;end%消元for i=k+1:nm(i,k)=A(i,k)/A(k,k); %A(k,k)将A(i,k)消为0所乘系数A(i,k:n)=A(i,k:n)-m(i,k)*A(k,k:n); %第i行消元处理b(i)=b(i)-m(i,k)*b(k); %b消元处理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);endclear x;disp('AX=b的解x是') x=b;【调用函数解题】【列主元三角分解法Matlab主程序】clear;augm=input('Please input the augumental matric:');[n,t]=size(augm);for i=1:n %d控制列三角主元的层数b=augm(i,:);p=i; %4到10行是选列主元并交换for w=i+1:nif b(1,1)<=augm(w,i)b=augm(w,:);p=w;endendc=augm(i,:); augm(i,:)=b;augm(p,:)=c;if p~=i %只有发生了换行才将这种效果输出augm,endfor j=i:t %首先变换与augm（i,i）同一行的元素，其列指标从i到t s=0;for k=1:i-1s=s+augm(i,k)*augm(k,j);endaugm(i,j)=(augm(i,j)-s);endfor I=i+1:n%再变换与augm（i,i）同一列的元素，其行指标I从i+1到n,列数为i s=0;for k=1:i-1 %下三角部分列数与第i层一致s=s+augm(I,k)*augm(k,i);endaugm(I,i)=(augm(I,i)-s)/augm(i,i);endaugm,endx(n)=augm(n,t)/augm(n,n); %回代for i=n-1:-1:1s=0;for j=n:-1:i+1s=s+x(j)*augm(i,j);endx(i)=(augm(i,t)-s)/augm(i,i),end【调用函数解题】输入[A b]可以清晰地看到augm（A）的变化过程：结果输出：【列主元三角分解法程序流程图】【编程疑难】这是第一次用matlab编程，对matlab的语句还不是非常熟悉，因此在编程过程中，出现了许多错误提示。

并且此次编程的两种方法对矩阵的运算也比较复杂。

问题主要集中在循环控制中，循环次数多了一次或者缺少了一次，导致数据错误，一些基本的编程语句在语法上也会由于生疏而产生许多问题，但是语句的错误由于系统会提示，比较容易进行修改，数据计算过程中的一些逻辑错误，比如循环变量的控制，这些系统不会提示错误，需要我们细心去发现错误，不断修正，调试。