列主元高斯消去法和列主元三角分解法解线性方程.docx
《列主元高斯消去法和列主元三角分解法解线性方程.docx》由会员分享，可在线阅读，更多相关《列主元高斯消去法和列主元三角分解法解线性方程.docx（15页珍藏版）》请在冰豆网上搜索。
列主元高斯消去法和列主元三角分解法解线性方程
计算方法实验报告1
【课题名称】
用列主元高斯消去法和列主元三角分解法解线性方程
【目的和意义】
高斯消去法是一个古老的求解线性方程组的方法,但由它改进得到的选主元的高斯消去法则是目前计算机上常用的解低阶稠密矩阵方程组的有效方法。
用高斯消去法解线性方程组的基本思想时用矩阵行的初等变换将系数矩阵A约化为具有简单形式的矩阵（上三角矩阵、单位矩阵等）,而三角形方程组则能够直截了当回带求解
用高斯消去法解线性方程组（其中A∈Rｎ×ｎ）的计算量为:
乘除法运算步骤为
加减运算步骤为
、相比之下,传统的克莱姆法则则较为繁琐,如求解20阶线性方程组,克莱姆法则大约要次乘法,而用高斯消去法只需要3０６0次乘除法。
在高斯消去法运算的过程中,假如出现abs（A（ｉ,i））等于零或过小的情况,则会导致矩阵元素数量级严重增长和舍入误差的扩散,使得最后的计算结果不可靠,因此目前计算机上常用的解低阶稠密矩阵方程的快速有效的方法时列主元高斯消去法,从而使计算结果更加精确、
2、列主元三角分解法
高斯消去法的消去过程,实质上是将A分解为两个三角矩阵的乘积A＝LU,并求解Ly=b的过程。
回带过程就是求解上三角方程组Ux=y、因此在实际的运算中,矩阵L和Ｕ能够直截了当计算出,而不需要任何中间步骤,从而在计算过程中将高斯消去法的步骤进行了进一步的简略,大大提高了运算速度,这就是三角分解法
采纳选主元的方式与列主元高斯消去法一样,也是为了幸免除数过小,从而保证了计算的精确度
【计算公式】
1、列主元高斯消去法
设有线性方程组Ax=b,其中设A为非奇异矩阵。
方程组的增广矩阵为
第1步（k=１）:
首先在A的第一列中选取绝对值最大的元素,作为第一步的主元素:
然后交换（A,b）的第1行与第l行元素,再进行消元计算。
设列主元素消去法差不多完成第１步到第ｋ-1步的按列选主元,交换两行,消元计算得到与原方程组等价的方程组A（k）ｘ=b（k）
第ｋ步计算如下:
关于k=1,2,…,n-1
（１）按列选主元:
即确定t使
（2）假如ｔ≠ｋ,则交换［A,b]第t行与第k行元素。
（3）消元计算
消元乘数mik满足:
（4）回代求解
2、列主元三角分解法
对方程组的增广矩阵 经过k-1步分解后,可变成如下形式:
第k步分解,为了幸免用绝对值特别小的数作除数,引进量
因此有=。
假如 ,则将矩阵的第t行与第k行元素互换,将（i,j）位置的新元素仍记为或,然后再做第ｋ步分解,这时
【列主元高斯消去法程序流程图】
【列主元高斯消去法Matlab主程序】
fuｎcｔionx=ｇaｕsｓ1（A,b,c） %列主元法高斯消去法解线性方程Ax=b
if（lengtｈ（A）～=ｌeｎgth（ｂ）） %判断输入的方程组是否有误
disｐ（’输入方程有误!
’）
retuｒn;
enｄ
ｄisｐ（'原方程为AＸ=b:
’） ％显示方程组
A
ｂ
disp（’----—-－---——－—--—--－-—-－＇）
n=lenｇｔｈ（A）;
forｋ=1:
n—1 ％找列主元
［p,q］=max（aｂs（A（ｋ:
n,ｋ）））; ％找出第k列中的最大值,其下标为[p,q］
q=q+k－１; ％q在A（k:
n,k）中的行号转换为在A中的行号
ｉfａbs（p）〈c
dｉsp（＇列元素太小,ｄｅｔ（A）≈0'）;
ｂreａｋ;
eｌseiｆｑ〉ｋ
tｅｍp1=A（k,:
）; %列主元所在行不是当前行,将当前行与列主
Ａ（k,:
）=A（ｑ,:
）; 元所在行交换（包括ｂ）
A（q,:
）=temp１;
temｐ2=ｂ（k,:
）;
b（k,:
）=b（ｑ,:
）;
ｂ（q,:
）＝ｔｅｍp2;
enｄ
％消元
fori＝k+1:
n
m（i,ｋ）=Ａ（i,ｋ）/Ａ（k,ｋ）; %Ａ（k,k）将A（i,ｋ）消为0所乘系数
A（i,k:
n）=Ａ（i,k:
n）-m（i,k）*A（ｋ,k:
n）;％第ｉ行消元处理
b（ｉ）=b（i）—m（i,k）*b（k）; ％b消元处理
eｎd
ｅnd
dｉｓp（’消元后所得到的上三角阵是＇）
A ％显示消元后的系数矩阵
ｂ（n）=b（n）/Ａ（n,n）; ％回代求解
ｆoｒｉ=ｎ－1:
－1:
１
b（i）＝（ｂ（i）—ｓum（A（ｉ,ｉ+１:
n）＊ｂ（i+1:
n）））/Ａ（ｉ,i）;
end
ｃlｅar ｘ;
disp（’AX＝ｂ的解x是＇）
x＝ｂ;
【调用函数解题】
【列主元三角分解法程序流程图】
【列主元三角分解法Maｔlab主程序】
①自己编的程序:
funcｔioｎｘ=PＬU（A,b,eｐs） ％定义函数列主元三角分解法函数
if（ｌeｎgth（Ａ）~＝lengtｈ（ｂ）） ％判断输入的方程组是否有误
ｄisp（＇输入方程有误！
'）
return;
end
diｓp（’原方程为AX=b:
'） %显示方程组
Ａ
b
ｄisp（’—-－—-－——-－—-----－－－—-—-—’）
n＝lenｇｔｈ（A）;
A＝[Ab]; %将Ａ与b合并,得到增广矩阵
forr=１:
n
iｆr==1
fori=１:
n
[ｃd］=ｍａx（aｂｓ（A（:
1）））; ％选取最大列向量,并做行交换
ｉfc＜=eｐs %最大值小于ｅ,主元太小,程序结束
ｂreak;
elsｅ
ｅnｄ
d=ｄ+１—1;
p=Ａ（１,:
）;
A（1,:
）=A（d,:
）;
A（ｄ,:
）＝p;
A（1,ｉ）=A（1,i）;
end
A（１,2:
n）=Ａ（1,2:
ｎ）;
A（2:
n,１）=A（2:
n,1）/A（1,1）; ％求u（１,i）
else
u（r,ｒ）=A（r,r）-A（r,1:
ｒ—１）*A（１:
r-１,ｒ）; %依照方程求取u（ｒ,i）
ｉｆabｓ（u（r,r））<=eps ％假如u（r,r）小于e,则交换行
ｐ=A（r,:
）;
Ａ（r,:
）=A（ｒ＋1,:
）;
Ａ（r＋1,:
）＝p;
elsｅ
end
for ｉ＝r:
n
A（ｒ,ｉ）＝A（ｒ,ｉ）—Ａ（ｒ,1:
r-1）＊A（1:
r—1,i）; ％依照公式求解,并把结果存在矩阵A中
ｅnd
fori=r+1:
ｎ
A（ｉ,r）=（A（i,r）－Ａ（i,1:
r—1）*A（1:
r-1,ｒ））／A（r,r）;％依照公式求解,并把结果存在矩阵A中
end
enｄ
eｎd
ｙ
（1）＝Ａ（1,ｎ+1）;
fｏr i=2:
ｎ
h=０;
fｏr k=1:
i—１
h=h+Ａ（i,ｋ）＊y（k）;
enｄ
ｙ（ｉ）=Ａ（i,ｎ+1）－h; ％依照公式求解y（i）
end
x（ｎ）=y（ｎ）/Ａ（n,ｎ）;
forｉ=n—１:
—1:
１
h＝0;
fork=i＋1:
n
ｈ=h+Ａ（ｉ,ｋ）＊ｘ（ｋ）;
end
x（i）=（y（ｉ）-h）/A（ｉ,i）; %依照公式求解x（i）
enｄ
Ａ
ｄｉｓp（＇ＡＸ=b的解ｘ是'）
x＝x'; ％输出方程的解
②可直截了当得到Ｐ,Ｌ,U并解出方程解的的程序（查阅资料得子函数PLU１,其作用是将矩阵A分解成L乘以U的形式、PLU２为调用PLU1解题的程序,是自己编的）
（Ⅰ）、
ｆuncｔion[l,ｕ,ｐ］＝PLU1（Ａ） %定义子函数,其功能为列主元三角分解系数矩阵A
[m,n］＝size（A）; %判断系数矩阵是否为方阵
iｆm～=n
ｅrｒoｒ（'矩阵不是方阵'）
ｒeturｎ
end
iｆdｅt（A）＝=0 %判断系数矩阵能否被三角分解
ｅｒror（'矩阵不能被三角分解'）
ｅnｄ
ｕ=A;p＝eｙｅ（m）;ｌ=eyｅ（ｍ）; %将系数矩阵三角分解,分别求出Ｐ,L,Ｕ
for ｉ=１:
ｍ
foｒ j=i:
ｍ
ｔ（j）=ｕ（ｊ,i）;
ｆｏrｋ=1:
i—1
t（ｊ）=t（ｊ）-u（ｊ,k）＊u（k,i）;
enｄ
eｎd
ａ＝i;ｂ=abs（t（i））;
foｒｊ＝ｉ+1:
m
ｉｆb b=abs（t（j））;
ａ=j;
eｎd
end
iｆ ａ～=i
foｒ j＝1:
m
c＝u（i,j）;
ｕ（i,ｊ）＝ｕ（a,j）;
u（a,j）=c;
eｎｄ
for ｊ＝1:
m
ｃ=ｐ（i,j）;
ｐ（i,j）=p（a,j）;
ｐ（ａ,j）=c;
ｅnｄ
c=t（a）;
t（ａ）＝ｔ（i）;
ｔ（i）=ｃ;
end
u（i,i）=t（ｉ）;
ｆｏr j=i＋1:
m
ｕ（j,ｉ）=t（j）/t（i）;
ｅnｄ
forj=i＋1:
m
forｋ=１:
i-1
u（i,j）=u（i,j）—u（i,ｋ）*u（k,j）;
enｄ
ｅnd
end
ｌ＝ｔril（u,-1）＋eye（m）;
u=trｉｕ（u,0）
（Ⅱ）、
functionx＝ＰLU2（A,b） ％定义列主元三角分解法的函数
[ｌ,u,p]=ＰLU1（A） %调用ＰＬU分解系数矩阵A
m＝lenｇth（A）; ％由于A左乘p,故b也要左乘p
v=b;
forq＝1:
m
b（q）=ｓｕm（ｐ（ｑ,１:
ｍ）＊v（１:
m,１））;
ｅnd
b
（1）=b（１） ％求解方程Ly=b
foｒi＝2:
１:
m
b（i）=（b（i）－ｓuｍ（ｌ（i,１:
ｉ-1）＊b（１:
i－１）））;
end
b（m）＝ｂ（m）/ｕ（m,m）; ％求解方程Uｘ=y
fｏrｉ=ｍ-1:
-１:
1
b（ｉ）=（b（ｉ）—ｓum（u（i,i+1:
m）*ｂ（ｉ＋１:
ｍ）））／u（i,i）;
ｅnｄ
ｃleaｒx;
disp（'ＡＸ=b的解x是’）
ｘ=b;
【调用函数解题】
①
②
【编程疑难】
这是第一次用maｔlａb编程,对mａｔｌab的语句还不是特别熟悉,因此在编程过程中,出现了许多错误提示。
同时此次编程的两种方法对矩阵的运算也比较复杂。
问题主要集中在循环控制中,循环次数多了一次或者缺少了一次,导致数据错误,一些基本的编程语句在语法上也会由于生疏而产生许多问题,然而语句的错误由于系统会提示,比较容易进行修改,数据计算过程中的一些逻辑错误,比如循环变量的控制,这些系统可不能提示错误,需要我们细心去发现错误,不断修正,调试。
升级会员 