列主元高斯消去法和列主元三角分解法解线性方程Word文档格式.docx
《列主元高斯消去法和列主元三角分解法解线性方程Word文档格式.docx》由会员分享，可在线阅读，更多相关《列主元高斯消去法和列主元三角分解法解线性方程Word文档格式.docx（18页珍藏版）》请在冰豆网上搜索。
（4）回代求解
1、
列主元三角分解法
对方程组的增广矩阵经过k-1步分解后，可变成如下形式：
第k步分解，为了避免用绝对值很小的数
作除数，引进量
，于是有
=
如果，则将矩阵的第t行与第k行元素互换，将（i，j）位置的新元素仍记为
或
，然后再做第k步分解，这时
【列主元高斯消去法程序流程图】
【列主元高斯消去法Matlab主程序】
functionx=gauss1（A,b,c）%列主元法高斯消去法解线性方程Ax=b
if（length（A）~=length（b））%判断输入的方程组是否有误
disp（'
输入方程有误！
'
）
return;
end
disp（'
原方程为AX=b：
）%显示方程组
A
b
------------------------'
n=length（A）;
fork=1:
n-1%找列主元
[p,q]=max（abs（A（k:
n,k）））;
%找出第k列中的最大值，其下标为[p,q]
q=q+k-1;
%q在A（k:
n,k）中的行号转换为在A中的行号
ifabs（p）<
c
列元素太小，det（A）≈0'
）;
break;
elseifq>
k
temp1=A（k,:
%列主元所在行不是当前行，将当前行与列主
A（k,:
）=A（q,:
元所在行交换（包括b）
A（q,:
）=temp1;
temp2=b（k,:
b（k,:
）=b（q,:
b（q,:
）=temp2;
%消元
fori=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
消元后所得到的上三角阵是'
）
A%显示消元后的系数矩阵
b（n）=b（n）/A（n,n）;
%回代求解
fori=n-1:
-1:
1
b（i）=（b（i）-sum（A（i,i+1:
n）*b（i+1:
n）））/A（i,i）;
clearx;
AX=b的解x是'
x=b;
【调用函数解题】
【列主元三角分解法程序流程图】
【列主元三角分解法Matlab主程序】
①自己编的程序：
functionx=PLU（A,b,eps）%定义函数列主元三角分解法函数
if（length（A）~=length（b））%判断输入的方程组是否有误
）%显示方程组
A=[Ab];
%将A与b合并，得到增广矩阵
forr=1:
ifr==1
fori=1:
[cd]=max（abs（A（:
1）））;
%选取最大列向量，并做行交换
ifc<
=eps%最大值小于e，主元太小，程序结束
else
d=d+1-1;
p=A（1,:
A（1,:
）=A（d,:
A（d,:
）=p;
A（1,i）=A（1,i）;
A（1,2:
n）=A（1,2:
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）
ifabs（u（r,r））<
=eps%如果u（r,r）小于e，则交换行
p=A（r,:
A（r,:
）=A（r+1,:
A（r+1,:
fori=r:
A（r,i）=A（r,i）-A（r,1:
r-1,i）;
%根据公式求解，并把结果存在矩阵A中
fori=r+1:
A（i,r）=（A（i,r）-A（i,1:
r-1,r））/A（r,r）;
%根据公式求解，并把结果存在矩阵A中
y
（1）=A（1,n+1）;
fori=2:
h=0;
fork=1:
i-1
h=h+A（i,k）*y（k）;
y（i）=A（i,n+1）-h;
%根据公式求解y（i）
x（n）=y（n）/A（n,n）;
fori=n-1:
fork=i+1:
h=h+A（i,k）*x（k）;
x（i）=（y（i）-h）/A（i,i）;
%根据公式求解x（i）
A
x=x'
;
%输出方程的解
②可直接得到P,L,U并解出方程解的的程序（查阅资料得子函数PLU1，其作用是将矩阵A分解成L乘以U的形式。
PLU2为调用PLU1解题的程序，是自己编的）
（Ⅰ）.
function[l,u,p]=PLU1（A）%定义子函数，其功能为列主元三角分解系数矩阵A
[m,n]=size（A）;
%判断系数矩阵是否为方阵
ifm~=n
error（'
矩阵不是方阵'
return
ifdet（A）==0%判断系数矩阵能否被三角分解
矩阵不能被三角分解'
u=A;
p=eye（m）;
l=eye（m）;
%将系数矩阵三角分解，分别求出P,L,U
m
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+1:
ifb<
abs（t（j））
b=abs（t（j））;
a=j;
ifa~=i
forj=1:
c=u（i,j）;
u（i,j）=u（a,j）;
u（a,j）=c;
c=p（i,j）;
p（i,j）=p（a,j）;
p（a,j）=c;
end
c=t（a）;
t（a）=t（i）;
t（i）=c;
u（i,i）=t（i）;
u（j,i）=t（j）/t（i）;
u（i,j）=u（i,j）-u（i,k）*u（k,j）;
l=tril（u,-1）+eye（m）;
u=triu（u,0）
（Ⅱ）.
functionx=PLU2（A,b）%定义列主元三角分解法的函数
[l,u,p]=PLU1（A）%调用PLU分解系数矩阵A
m=length（A）;
%由于A左乘p，故b也要左乘p
v=b;
forq=1:
b（q）=sum（p（q,1:
m）*v（1:
m,1））;
b
（1）=b
（1）%求解方程Ly=b
fori=2:
1:
b（i）=（b（i）-sum（l（i,1:
i-1）*b（1:
i-1）））;
b（m）=b（m）/u（m,m）;
%求解方程Ux=y
fori=m-1:
b（i）=（b（i）-sum（u（i,i+1:
m）*b（i+1:
m）））/u（i,i）;
①
②
【编程疑难】
这是第一次用matlab编程，对matlab的语句还不是非常熟悉，因此在编程过程中，出现了许多错误提示。
并且此次编程的两种方法对矩阵的运算也比较复杂。
问题主要集中在循环控制中，循环次数多了一次或者缺少了一次，导致数据错误，一些基本的编程语句在语法上也会由于生疏而产生许多问题，但是语句的错误由于系统会提示，比较容易进行修改，数据计算过程中的一些逻辑错误，比如循环变量的控制，这些系统不会提示错误，需要我们细心去发现错误，不断修正，调试。
升级会员 