列主元消去法解方程组实验报告.docx
《列主元消去法解方程组实验报告.docx》由会员分享，可在线阅读，更多相关《列主元消去法解方程组实验报告.docx（11页珍藏版）》请在冰豆网上搜索。
列主元消去法解方程组实验报告
实验名称：
列主元消去法解方程组
1引言
我们知道，高斯消去法是一个古老的解线性方程组的方法。
而在用高斯消去法解Ax=b时，其中设A为非奇异矩阵，可能出现
的情况，这时必须进展带行交换的高斯消去法。
但在实际计算中即使
但其绝对值很小时，用
作除数，会导致中间结果矩阵
元素数量级严重增长和舍入误差的扩散，使得最后的结果不可靠。
因此，小主元可能导致计算的失败，我们应该防止采用绝对值很小的主元素。
为此，我们在高斯消去法的每一步应该在系数矩阵或消元后的低阶矩阵中选取绝对值最大的元素作为主元素，保持乘数
，以便减少计算过程中舍入误差对计算解的影响。
一种方式是完全主元消去法，这种消去法是在每次选主元时，选择
为主元素。
这种方法是解低阶稠密矩阵方程组的有效方法，但这种方法在选取主元时要花费一定的计算机时间。
实际计算中我们常采用局部选主元的的消去法。
列主元消去法即在每次选主元时，仅依次按列选取绝对值最大的元素作为主元素，且仅交换两行，再进展消元计算。
2实验目的和要求
运用matlab编写一个.m文件，要求用列主元消去法求解方程组〔实现PA=LU〕：
要求输出以下容：
（1）计算解x；
（2）L,U；
（3）整形数组IP〔i〕〔i=1,2，…，n-1〕〔记录主行信息〕
3算法原理与流程图
（1）算法原理
设有线性方程组Ax=b，其中设A为非奇异矩阵。
方程组的增广矩阵为
第1步〔k=1〕：
首先在A的第一列中选取绝对值最大的元素
，作为第一步的主元素:
，然后交换〔A，b〕的第1行与第i1行元素，再进展消元计算。
设列主元素消去法已经完成第1步到第k-1步的按列选主元，交换两行，消元计算得到与原方程组等价的方程组
第k步计算如下：
对于k=1，2，…，n-1
〔1〕按列选主元：
即确定ik使
〔2〕如果
，那么A为非奇异矩阵，停顿计算。
〔3〕如果ik≠k，那么交换[A，b]第ik行与第k行元素。
〔4〕消元计算
消元乘数
满足:
〔5〕回代求解
计算解
在常数项b（n）得到。
（2）流程图见图1
（3）
4程序代码及注释
%列主元消去法解方程组Ax=b，实现PA=LU
function[x,L,U,IP,P]=gauss（A,b）
%x为方程组的解，IP用来记录行信息
%每次选列主元时，将A的第k行与第IP（k）行进展交换
n=length（b）;[p,q]=size（A）;
%当输入的系数矩阵不为方阵，或方阵维数与b不符时，报错
ifp~=q||p~=n
fprintf（'Error!
Pleaseinputagain!
'）;
end
%为提高运行速度，给IP，P，L，U赋初值
IP=zeros（1,n-1）;
L=zeros（n,n）;
U=zeros（n,n）;
P=eye（n）;
x=zeros（1,n）;
det=1.0;
%按列选主元，并进展行交换，记录行信息
fork=1:
n-1
IP（k）=k;
form=k+1:
n
ifabs（A（m,k））>abs（A（k,k））
IP（k）=m;
end
end
I=eye（n）;
ifIP（k）~=k
fori=1:
n
p（i）=I（k,i）;
I（k,i）=I（IP（k）,i）;
I（IP（k）,i）=p（i）;
end
A=I*A;
b=I*b';
b=b';
end
%进展消元计算
fori=k+1:
n
A（i,k）=A（i,k）/A（k,k）;
b（i）=b（i）-A（i,k）*b（k）;
forj=k+1:
n
A（i,j）=A（i,j）-A（i,k）*A（k,j）;
end
det=det*A（k,k）;
P=I*P;
end
%回代求解
x（n）=b（n）/A（n,n）;
fori=n-1:
-1:
1
sum=0.0;
forj=i+1:
n
sum=sum+A（i,j）*x（j）;
end
x（i）=（b（i）-sum）/A（i,i）;
end
det=det*A（n,n）;
ifdet==0
fprintf（'Theequationshavenouniquesolution!
'）;
end
%输出PA=LU中的L,U的信息
fori=1:
n
forj=1:
n
ifiU（i,j）=A（i,j）;
elseifi==j
L（i,j）=1;
U（i,j）=A（i,j）;
else
L（i,j）=A（i,j）;
end
end
end
end
附图1
5算例分析
1、测试例如
>>A=[
12
34
56];
>>b=[3711];
>>[xLUIPP]=gauss（A,b）
Error!
Pleaseinputagain!
>>A=[
123
000
456];
>>b=[3711];
>>[xLUIPP]=gauss（A,b）
Theequationshavenouniquesolution!
x=
NaN-InfInf
L=
1.000000
0.25001.00000
001.0000
U=
4.00005.00006.0000
00.75001.5000
000
IP=
33
P=
100
001
010
2、计算过程
〔1〕首先输入系数矩阵A和矩阵b
>>A=[
1111111
2111111
3211111
4321111
5432111
6543211
7654321];
>>b=[781013172228];
〔2〕输出结果
>>[xLUIPP]=gauss（A,b）
x=
1.00001.00001.00001.00001.00001.00001.0000
L=
1.0000000000
0.28571.000000000
0.42860.80001.00000000
0.57140.60000.75001.0000000
0.71430.40000.50000.66671.000000
0.1429-0.2000-0.2500-0.3333-0.50001.00000
0.85710.20000.25000.33330.5000-1.00001.0000
U=
7.00006.00005.00004.00003.00002.00001.0000
0-0.7143-0.4286-0.14290.14290.42860.7143
00-0.8000-0.6000-0.4000-0.20000.0000
000-0.7500-0.5000-0.25000.0000
0000-0.6667-0.3333-0.0000
000000.50001.0000
0000001.0000
IP=
723457
P=
1000000
0100000
0010000
0001000
0000100
0000001
0000010
6讨论与结论
1、时间复杂度：
>>tic;[xLUIPP]=gauss（A,b）;toc
Elapsedtimeis0.000856seconds.
2、程序优化
初次编程时，没有考虑到给一个变量赋初值的情况。
虽然在MATLAB中变量不赋初值是完全允许的，但是由于一个变量中含有多个元素时，每次改变该数组的长度，便会增加计算机时间。
另外，给程序加上一定的判断条件及报错信息，一定程度上有程序优化的作用。
因此，本程序中的以下程序段都起到了程序优化的作用。
%当输入的系数矩阵不为方阵，或方阵维数与b不符时，报错
ifp~=q||p~=n
fprintf（'Error!
Pleaseinputagain!
'）;
end
%为提高运行速度，给IP，P，L，U赋初值
IP=zeros（1,n-1）;
L=zeros（n,n）;
U=zeros（n,n）;
P=eye（n）;
x=zeros（1,n）;
det=1.0;
参考文献
[1]易大义，云宝，有法.计算方法（第2版），大学.p.29-53.
[2]琨高思超毕靖编著MATLAB2010从入门到精通电子工业
升级会员 