列主元高斯消去法和列主元三角分解法解线性方程041514Word格式文档下载.docx
《列主元高斯消去法和列主元三角分解法解线性方程041514Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《列主元高斯消去法和列主元三角分解法解线性方程041514Word格式文档下载.docx（14页珍藏版）》请在冰豆网上搜索。
1、列主元高斯消去法
设有线性方程组Ax=b,其中设A为非奇异矩阵。
方程组的增广矩阵为
第1步（kJ）:
首先在A的第一列中选取绝对值最大的元素切，作为第一步
的主兀素:
然后交换（A,b）的第1行与第I行元素，再进行消元计算。
设列主元素消去法已经完成第1步到第k<
L步的按列选主元，交换两行，消元
计算得到与原方程组等价的方程组A（k）x二b（k）
第k步计算如下：
对于k=l»
2,…，n-1
（1）按列选主元：
即确定t使
（2）如果tHk,则交换［A,b］第t行与第k行元素。
（3）消元计算
%J片加=——,（i=k+1,・••,〃）
Clkk
aj%+叫,（：
j=*+1,…/）
bj叫片（j=k+l,・・・M）
（4）回代求解
2、列主元三角龙卿•法
对方程组的增广矩卩轻=［人如过k-1步分解后，可变成如下形式:
M11
W12
…5"
%
•••
…％
>
1
W22
…"
心
li2k
♦••
皿j
…"
2”
■
•
91.2
m
Uk-\.k
3
in
…
♦♦•
…%
bk
♦
第k步分解，为了避免用绝对值很小的数5作除数，引进量
'
女-】
ukj一E—%（丿・=£
k+1,…M；
*=1,2,…/）
m-1
<
&
-1
/汝=（血一）/%（，=k+1,£
+2,・・•/;
=1,2,•••丿）
.wr-l
行与第k行元素互换，将（i,j）位置的新元素仍记为'
力或5,然后再做第k步分解，这时
%=sk（sk即交换前的耳），
lik=6Isk（j=k+l,k+2,・・・,〃）
且|Z;
J<
1（i=R+l,R+2,•••,“），
【列主元高斯消去法程序流程图】
%列主元法高斯消去法解线性方程Ax=b
%判断输入的方程组是否有误
【列主元高斯消去法Matlab主程序】
functionx=gaussl（A,b,c）if（length（A）"
=length（b））disp（'
输入方程有误！
）return;
end
disp（*原方程为AX=b：
'
）%显示方程组
A
b
dispC•）n=length（A）;
fork=l:
n-l%找列主元
[p,q]=max（abs（A（k:
n,k）））;
%找出第k列中的最大值，其下标为[p,q]
q=q+k-l;
%q在A（k:
n,k）中的行号转换为在A中的行号
讦abs（p）<
c
disp（'
列元素太小，det（A）QO'
）;
break;
%列主元所在行不是当前行，将当前行与列主元所在行交换（包括b）
else讦q>
ktempl=A（k,:
A（k,:
）=A（q/:
A（q/）=templ;
temp2=b（k,:
b（k,:
）=b（q,:
b（q,:
）=temp2;
%消元
fori=k+l:
n
m（i,k）=A（i,k）/A（k,k）;
%A（k,k）将A（i,k）i肖为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消元处理
dispC消元后所得到的上三角阵是J
A%显示消元后的系数矩阵
b（n）=b（n）/A（n,n）;
%回代求解・
fori=n-l:
-l:
l
b（i）=（b（i）-sum（A（izi+l:
n）*b（i+l:
n）））/A（izi）;
clearx;
AX=b的解x是J
x=b;
【调用函数解题】
CommandWindow
»
A=[034;
1-11：
212];
b=[l;
2;
3];
c二0・00001:
x=gaussl（Ajb,c）
原方程为AX二b:
034
1-11
212
b二
2
消元后所得到的上三角蹲是
A=
002
AS二b的斛垃是
X=
1.1667
-0.3333
0.5000
【列主元三角分解法程序流程图］
【列主元三角分解法Matlab主程序】①自己编的程序：
functionx=PLU（A,b,eps）if（length（A）"
）
return;
disp（■原方程为AX=b：
%定义函数列主元三角分解法函数
%显示方程组
Ab
disp（**）
n=length（A）;
A二[Ab];
%将人与b合并，得到增广矩阵
forr=l:
ifr==l
fori=l:
[cd]=max（abs（A（:
l）））;
%选取最大列向量，并做行交换
讦c<
=eps%最大值小于e,主元太小，程序结束
break;
else
d=d+l-l;
P=A（1,：
）；
A（l,:
）=A（d,:
A（d,:
）=p;
A（l,i）=A⑴i）;
A（l,2:
n）=A（lz2:
A（2:
n/l）=A（2:
n/l）/A（l,l）;
%求u（lzi）
elseu（r/r）=A（r/r）-A（r,l:
r-l）*A（l:
r-lzr）;
%按照方程求取u（r,i）
ifabs（u（r,r））<
=eps%如果u（r,r）小于e,则交换行P=A（r,:
A（r,:
）=A（r+l/:
A（r+l/）=p;
fori=r:
A（rzi）=A（r,i）-A（r,l:
r-l/i）;
%根据公式求解，并把结果存在矩阵A中
fori=r+l:
A（i/r）=（A（i/r）-A（i/l：
r-l,r））/A（r,r）;
%根据公式求解,并把结果存在矩阵A中end
y（l）=A（l,n+l）;
fori=2:
h=0;
i-l
h=h+A（i,k）*y（k）;
y（i）=A（i,n+l）-h;
%根据公式求解y（i）
x（n）=y（n）/A（n,n）;
fork=i+l:
h=h+A（i/k）*x（k）;
x（i）=（y（i）-h）/A（i,i）;
%根据公式求解x（i）
AX=b的解x是'
x=x]%输出方程的解
②可直接得到P儿U并解出方程解的的程序（查阅资料得子函数PLU1,其作用是将矩阵A分解成L乘以U的形式。
PLU2为调用PLU1解题的程序，是自己编的）
（I）.
function[I,u/P]=PLU1（A）%定义子函数，其功能为列主元三角分解系数矩阵A[m,n]=size（A）;
%判断系数矩阵是否为方阵
讦
errorf矩阵不是方阵J
return
ifdet（A）==O%判断系数矩阵能否被三角分解
error）'
矩阵不能被三角分解J
u=A;
p=eye（m）;
l=eye（m）;
%将系数矩阵三角分解，分别求出P,L,U
forj=i:
t（j）=u（j,i）;
t（j）=t（j）-u（j,k）*u（k,i）;
a=i;
b=abs（t（i））;
forj=i+l:
ifb<
abs（t（j））
b=abs（t（j））;
a=j;
讦a~=i
forj=l:
c=u（i,j）；
u（i,j）=u（aj）;
u（a,j）=c；
mc=P（iJ）；
P（iJ）=P（aJ）；
P（a,j）=c;
c=t（a）;
t（a）=t（0；
t（i）=c;
endu（i,i）=t（i）;
mu（j,i）=t（j）/t（i）;
endforj=i+l:
mfork=l:
i-lu（ij）=u（ij）-u（izk）*u（kj）;
l=tril（uz-l）+eye（m）;
u二triu（u’O）
（II）.
functionx=PLU2（A,b）
[Izu/P]=PLU1（A）
%定义列主元三角分解法的函数
%调用PLU分解系数矩阵A
m=length（A）;
v=b;
forq=l:
%由于A左乘p,故b也要左乘p
b（q）=sum（p（q/l:
m）*v（l:
m,l））;
end
b（l）=b（l）%求解方程Ly二b
l:
b（i）=（b（i）-sum（l（i,l:
i-l）
升级会员 