实验名称：列主元消去法解方程组
1引言
我们知道，高斯消去法是一个古老的解线性方程组的方法。
而在用高斯消去法解Ax=b 时，其中设A为非奇异矩阵，可能出现a kk)0的情况，这时必须进行带行交换的高斯消去法。
但在实际计算中即使a kk)0但其绝对值很小时，用a；：〉作除数，会导致中间结果矩阵A(k)元素数量级严重增长和舍入误差的扩散，使得最后的结果不可靠。
因此，小主元可能导致计算的失败，我们应该避免采用绝对值很小的主元素。
为此，我们在高斯消去法的每一步应该在系数矩阵或消元后的低阶矩阵中选取绝对值最大的元素作为主元素，保持乘数1，以便减少计算过程中舍入误差对计算解的影响。
一种方式是完全主元消去法，这种消去法是在每次选主元时，选择a(二k max a j k)0
kJ k k in」
k j n
为主元素。
这种方法是解低阶稠密矩阵方程组的有效方法，但这种方法在选取主元时要花费一定的计算机时间。
实际计算中我们常采用部分选主元的的消去法。
列主元消去法即在每次
选主元时，仅依次按列选取绝对值最大的元素作为主元素，且仅交换两行，再进行消元计算
2实验目的和要求
运用matlab编写一个.m文件，要求用列主元消去法求解方程组(实现PA=L)
111111
2 11111
3 2 1111
4 3 2 1 1 1
5 4 3 2 1 1
6 5 4 3 2 1
7 6 5 4 3 2 要求输出以下内容:
(1) 计算解x；
x.
117 1x28 1x310 1x413 1x517
1
x6
22 128 x 7
(2) L,U ；
(3) 整形数组IP (i ) (i=1,2，…，n-1 )(记录主行信息)
3算法原理与流程图
(1) 算法原理
设有线性方程组Ax=b，其中设A为非奇异矩阵。
方程组的增广矩阵为
an a12L L a1n M bi
a21a22L L a2n M b2
［代口
M
M M M
a i1
M M M M
a n1a n2L L a nn M
b n
第1步(k=1):首先在A的第一列中选取绝对值最大的元素a^，作为第一步的主元素3in max
耳！0，然后交换(A，b)的第1行与第i1行元素，再进行消元计算。
1 1 i n
设列主元素消去法已经完成第1步到第k-1步的按列选主元，交换两行，消元计算得到与原方程组等价的方程组A(k)x b(k)
J1) c ⑴
an ai2
L a(k)a1(n)M
⑵
L (2)⑵
M
a22a2k a2n b2()
O M M M [A,b] [A(k),b(k)]a kk)L a kk)M.(k)bk
a k k)u a k k)1,n M(k)
d 1
M M M M
a nk)L a(k)叫
n
M
第k步计算如下：
对于k=1，2,…，n-1
(1)按列选主元：即确定ik使
时max无0
(2) 如果a ik,k 0，则A为非奇异矩阵，停止计算
(3) 如果ik 丰k ,则交换[A , b]第ik 行与第k 行元素。
(4) 消元计算
(5) 回代求解
(3)
4程序代码及注释
for k=1:n-1
a ik
E ik
电，(i
a kk
1L ,n)
a j
b E k a
kj ,
mb.
(i ，j (i
k 1,L ,n)
k 1,L ,n)
消元乘数m ik 满足: E ik
a ik a kk
1,(i k 1,L ,n)
0 / a nn
(当 a
nn
b
(b i
j b n
n a ij b j ) / a ii ,
I
i 1
0) n 1,K ,2,1
计算解(X 1,L ,X n )在常数项 b(n)内得至叽
(2)
流程图见图1
IP(k)=k;
if i<j U(i,j)=A(i,j);
elseif i==j
L(i,j)=1； U(i,j)=A(i,j);
else
L(i,j)=A(i,j);
end
end
end
end
输入n, A , b,£
det=1.0
k=1,2,…，n-1
输出det( A) 0
按列选主元a k ,k max k i n a ik0
c0= a ik ,k
血
停机
(I
(j
c0|<
Y
ik=k
换行
a kj a ik ,j (j
k, L , n)
b k b ik,det det
消元计算
k+1,L
a ij
k+1,L
det
回代求解
b n
det
(b
a nn
a nn(当a nn0)
nn nn
n
a^ b j)/ a i ,I 1,2, L n
j i 1
* det
输岀计算解及行列式值b
(I)( I=1,2，…，n)及det
,n)
,n)
停机
,a ik
* a kj
,b i
a kk det
a ik / a kk
6
附图1
4.0000
5.0000
6.0000 0
5算例分析
1、测试示例 >> A=[ 1 2 3 4 5 6];
>> b=[3 7 11];
>> [x L U IP P]=gauss(A,b) Error! Please in put aga in! >> A=[ 1 2 3 0 0 0 4 5 6]; >> b=[3 7 11];
>> [x L U IP P]=gauss(A,b)
The equati ons have no unique soluti on! x = NaN -Inf Inf
1.0000
0.2500 1.0000
0 1.0000
0 0.7500 1.5000
IP =
3 3
P =
1 0 0
0 0 1
0 1 0
2、计算过程
(1)首先输入系数矩阵A和矩阵b
>> A=[
1 1 1 1 1 1 1
2 1 1 1 1 1 1
3 2 1 1 1 1 1
4 3 2 1 1 1 1
5 4 3 2 1 1 1
6 5 4 3 2 1 1
7 6 5 4 3 2 1];
>> b=[7 8 10 13 17 22 28];
(2)输出结果
>> [x L U IP P]=gauss(A,b)
x =
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
1.000000 0000
0.2857 1.000000 000
0.42860.8000 1.000000 00
0.57140.60000.7500 1.0000000
0.71430.40000.50000.6667 1.000000
0.1429-0.2000-0.2500-0.3333-0.5000 1.00000
0.85710.20000.25000.33330.5000 --1.0000 1.0000
7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000
0 -0.7143 -0.4286 -0.1429 0.1429 0.4286 0.7143
0 -0.8000 -0.6000 -0.4000 -0.2000 0.0000
0 -0.7500 -0.5000 -0.2500 0.0000
0 -0.6667 -0.3333 -0.0000
0 0.5000 1.0000
0 1.0000
IP =
6讨论与结论
1、时间复杂度：
>> tic;[x L U IP P]=gauss(A,b);toc
Elapsed time is 0.000856 sec on ds.
2、程序优化
初次编程时，没有考虑到给一个变量赋初值的情况。
虽然在MATLAB^变量不赋初值是完全允许的，但是由于一个变量中含有多个元素时，每次改变该数组的长度，便会增加计算机时间。
另外，给程序加上一定的判断条件及报错信息，一定程度上有程序优化的作用。
因此，本程序中的以下程序段都起到了程序优化的作用。
参考文献
[1] 易大义，沈云宝，李有法•计算方法(第2版)，浙江大学出版社.p.29-53.
[2] 张琨高思超毕靖编著MATLAB2010从入门到精通电子工业出版社
0 0 0 1 0 0 0
0 0 0 0 0 1 0。