计算方法实验报告1 【课题名称】 用列主元高斯消去法和列主元三角分解法解线性方程
【目的和意义】 高斯消去法是一个古老的求解线性方程组的方法，但由它改进得到的选主元的高斯消去法则是目前计算机上常用的解低阶稠密矩阵方程组的有效方法。 用高斯消去法解线性方程组的基本思想时用矩阵行的初等变换将系数矩阵A约化为具有简单形式的矩阵（上三角矩阵、单位矩阵等），而三角形方程组则可以直接回带求解 用高斯消去法解线性方程组bAx（其中A∈Rn×n）的计算量为：乘除法运算步骤为32(1)(1)(21)(1)(1)262233nnnnnnnnnnnMDn
，加减运算步骤为
(1)(21)(1)(1)(1)(25)6226nnnnnnnnnnAS。相比之下，传统的克莱姆
法则则较为繁琐，如求解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的第一列中选取绝对值最大的元素1la，作为第一步的主元素: 然后交换（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）消元计算
111211212222112[,]nn
nlnnnna
aaabaaab
aaab
ab
111max0liinaa
(1)(1)(1)(1)(1)
1112111(2)(2)(2)(2)22222()(()1)()()()()()1,1()(,)()[,][,] kkkknkkn
kn
kkkkk
kkknk
kk
knkkknnnaaaabaaabaababbaaa
AbAb
()()max0kktkikkinaa
,(1,,)ikikikkkaamikna, (,1,,), (1,,)ijijikkjiiikkaamaijknbbmbikn
消元乘数mik满足: （4）回代求解
2、 列主元三角分解法 对方程组的增广矩阵 经过k-1步分解后，可变成如下形式：
第k步分解，为了避免用绝对值很小的数kku作除数，引进量 1111 (,1,,;1,2,,) ()/ (1,2,,;1,2,,)kkjkjkmmjmkikikimmkkkmualujkknknlaluuikknkn
11(,1,,)kiikimmkmsaluikkn
，于是有kku=ks。如果 ，则将矩阵的第t
行与第k行元素互换，将（i，j）位置的新元素仍记为jjl或jja，然后再做第k步分解，这时
【列主元高斯消去法程序流程图】
)1,,2,1(,)(1nniaxabxab
x
iinijjijii
nnnn
[,]AAb11121,11111222,122221,11,1,1,211,11,2121,112,112,1kkkkkkkjnkkjnkkkiiiknnkkkjknkikijininknjkkkjknnnnkknaaabAauuuuuuyllllllllllluuuuuyuuuuyaabaabla
maxtikinss ()/ 1,2,,)1 (1,2,,),kkkktikikikussslssikknlikkn即交换前的，（且 【列主元高斯消去法Matlab主程序】 function x=gauss1(A,b,c) %列主元法高斯消去法解线性方程Ax=b if (length(A)~=length(b)) %判断输入的方程组是否有误 disp('输入方程有误！') return; end
disp('原方程为AX=b：') %显示方程组 A b disp('------------------------') 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) disp('列元素太小，det(A)≈0'); break;
elseif q>k temp1=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:n m(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消元处理 end end disp('消元后所得到的上三角阵是') A %显示消元后的系数矩阵 b(n)=b(n)/A(n,n); %回代求解 for i=n-1:-1:1 b(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主程序】 ①自己编的程序： function x=PLU(A,b,eps) %定义函数列主元三角分解法函数 if (length(A)~=length(b)) %判断输入的方程组是否有误 disp('输入方程有误！') return; end disp('原方程为AX=b：') %显示方程组 A b disp('------------------------') n=length(A); A=[A b]; %将A与b合并，得到增广矩阵 for r=1:n if r==1 for i=1:n [c d]=max(abs(A(:,1))); %选取最大列向量，并做行交换 if c<=eps %最大值小于e，主元太小，程序结束 break; else end d=d+1-1; p=A(1,:); A(1,:)=A(d,:); A(d,:)=p; A(1,i)=A(1,i); end A(1,2:n)=A(1,2:n); A(2:n,1)=A(2:n,1)/A(1,1); %求u(1,i) else u(r,r)=A(r,r)-A(r,1:r-1)*A(1:r-1,r); %按照方程求取u(r,i) if abs(u(r,r))<=eps %如果u（r,r）小于e，则交换行