高斯消元法解线性方程组.docx
《高斯消元法解线性方程组.docx》由会员分享，可在线阅读，更多相关《高斯消元法解线性方程组.docx（23页珍藏版）》请在冰豆网上搜索。
高斯消元法解线性方程组
实验一、高斯消元法解线性方程组
流程图：
KLi=input（'请输入你要解的方程的增广矩阵的列数:
'）;
e=input（'请输入你要求的非满秩矩阵的最小元绝对值:
'）;
fid=fopen（'Untitled.m','r'）;
A=fscanf（fid,'%f',[Li,inf]）;
A=A';
W=size（A）;
Ah=W
（1）;
Al=W
（2）;
forj=1:
Al-1;%找出矩阵的列主元并交换行
z=j;
fori=j:
Ah
ifabs（A（i,j））>=abs（A（z,j））
c=i;
i=z;
z=c;
else;
end
end
ifz==j
;
else
q=A（z,:
）;
A（z,:
）=A（j,:
）;
A（j,:
）=q;
end
ifabs（A（j,j））<=e
error（'方程系数矩阵为非满秩矩阵'）;
else;
end
end
forj=1:
Al-1
fori=j+1:
Ah
d=-A（i,j）/A（j,j）;
A（i,:
）=d*A（j,:
）+A（i,:
）;
end
end
W=size（A）;
Ah=W
（1）;
Al=W
（2）;
A1=A（1:
Ah,1:
Al-1）;
B=A（1:
Ah,Al）;
A=A1;
W=size（A）;
Ah=W
（1）;
Al=W
（2）;
X=zeros（Al,1）;
fori=Ah:
-1:
1
S=0;
forj=Al-1:
-1:
i
S=S+A（i,j+1）*X（j+1,1）;
end
X（i,1）=（B（i,1）-S）/A（i,i）;
end
X
运行结果：
附：
由P227第3题而产生的三角分解法
Li=input（'请输入你要解的方程的增广矩阵的列数:
'）;
e=input（'请输入你要求的非满秩矩阵的最小元绝对值:
'）;
fid=fopen（'Untitled.m','r'）;
A=fscanf（fid,'%f',[Li,inf]）;
A=A';
W=size（A）;
Ah=W
（1）;
Al=W
（2）;
Ae=A（1:
Ah,1:
Al-1）;
B=A（1:
Ah,Al）;
p=eye（Ah）;
forj=1:
Al-1;%找出矩阵的列主元并交换行
z=j;
fori=j:
Ah
ifabs（A（i,j））>=abs（A（z,j））
c=i;
i=z;
z=c;
else;
end
end
ifz==j
;
else
q=A（z,:
）;
A（z,:
）=A（j,:
）;
A（j,:
）=q;
pc=p（z,:
）;
p（z,:
）=p（j,:
）;
p（j,:
）=pc;
end
ifabs（A（j,j））<=e
error（'方程系数矩阵为非满秩矩阵'）;
else;
end
end
p;
forj=1:
Al-1
fori=j+1:
Ah
d=-A（i,j）/A（j,j）;
A（i,:
）=d*A（j,:
）+A（i,:
）;
end
end
U=A（1:
Ah,1:
Al-1）;
A=Ae;
L=（p*A）/U;
p
A
L
U
B=p*B
W=size（L）;%求解Ly=b
Ah=W
（1）;
Al=W
（2）;
Y=zeros（Al,1）;
fori=1:
Ah
S=0;
forj=1:
i-1
S=S+L（i,j）*Y（j,1）;
end
Y（i,1）=（B（i,1）-S）/L（i,i）;
end
Y
W=size（U）;%求解Ux=y
Ah=W
（1）;
Al=W
（2）;
X=zeros（Al,1）;
fori=Ah:
-1:
1
S=0;
forj=Al-1:
-1:
i
S=S+U（i,j+1）*X（j+1,1）;
end
X（i,1）=（Y（i,1）-S）/U（i,i）;
end
X
运行结果：
实验二、直接三角法解线性方程组
Li=input（'请输入你要解的方程的增广矩阵的列数:
'）;
e=input（'请输入你要求的非满秩矩阵的最小元绝对值:
'）;
fid=fopen（'Untitled.m','r'）;
A=fscanf（fid,'%f',[Li,inf]）;
A=A';
W=size（A）;
H=W
（1）;
AZ=A（1:
H,1:
H）;
B=A（1:
H,H+1）;
A=AZ;
forj=1:
H%求取矩阵U的第一行
ifabs（A（j,j））<=e
error（'方程系数矩阵为非满秩矩阵,请选择杨雁勇的下一个可以选择列主元的程序'）;
else;
end
end
forj=1:
H%求取矩阵U的第一行
A（1,j）=A（1,j）;
end
fori=2:
H%求取矩阵L的第1列
A（i,1）=A（i,1）/A（1,1）;
end
forr=2:
H%求取矩阵U的第2-H列和第2-H行
fori=r:
H
A（r,i）=A（r,i）-A（r,1:
r-1）*A（1:
r-1,i）;
end
fori=r+1:
H
A（i,r）=（A（i,r）-A（i,1:
r-1）*A（1:
r-1,r））/A（r,r）;
end
ifabs（A（r,r））<=e
error（'方程系数矩阵为非满秩矩阵,请选择杨雁勇的下一个可以选择列主元的程序'）;
else;
end
end
disp（'经过电算以后得到的A（由U与L中的元素置换而组成）'）
A
fori=1:
H%将A拆成L与U
forj=1:
H
ifi==j
L（i,j）=1;
elseifi>j
L（i,j）=A（i,j）;
else
L（i,j）=0;;
end
end
end
end
L
fori=1:
H%将A拆成L与U
forj=1:
H
ifi<=j
U（i,j）=A（i,j）;
else
U（i,j）=0;;
end
end
end
U
%回带求取解
W=size（L）;%求解Ly=b
Ah=W
（1）;
Al=W
（2）;
Y=zeros（Al,1）;
fori=1:
Ah
S=0;
forj=1:
i-1
S=S+L（i,j）*Y（j,1）;
end
Y（i,1）=（B（i,1）-S）/L（i,i）;
end
Y
W=size（U）;%求解Ux=y
Ah=W
（1）;
Al=W
（2）;
X=zeros（Al,1）;
fori=Ah:
-1:
1
S=0;
forj=Al-1:
-1:
i
S=S+U（i,j+1）*X（j+1,1）;
end
X（i,1）=（Y（i,1）-S）/U（i,i）;
end
X
运行结果：
实验三、列主元三角法解线性方程组
Li=input（'请输入你要解的方程的增广矩阵的列数:
'）;
e=input（'请输入你要求的非满秩矩阵的最小元绝对值:
'）;
fid=fopen（'Untitled.m','r'）;
A=fscanf（fid,'%f',[Li,inf]）;
A=A';
W=size（A）;
H=W
（1）;
AZ=A（1:
H,1:
H）;
B=A（1:
H,H+1）;
A=AZ;
K=1;
fori=1:
H
K=K*A（i,i）;
end
ifK==0
error（'方程中有对角元上元素=0'）;
else;
end
forj=1:
H%求取矩阵U的第一行
A（1,j）=A（1,j）;
end
fori=2:
H%求取矩阵L的第1列
A（i,1）=A（i,1）/A（1,1）;
end
forr=2:
H%求取矩阵U的第2-H列和第2-H行
%检验A（r,r）是否为一个特别小的数
ifabs（A（r,r））<=e
ss=zeros（r+1,1）;
s=[A（r,r）;ss];
fori=r+1:
H
S（i+1,1）=A（i,r）-A（i,1:
r-1）*A（1:
r-1,r）;
end
forj=r:
H%找出分块矩阵的列主元并交换行
z=j;
ifabs（S（r,1））>=abs（S（z,1））
c=j;
j=z;
z=c;
else;
end
end
ifz==r
;
else
q=A（z,r:
H）;
A（z,r:
H）=A（r,r:
H）;
A（r,r:
H）=q;
end
else;
end
fori=r:
H
A（r,i）=A（r,i）-A（r,1:
r-1）*A（1:
r-1,i）;
end
fori=r+1:
H
A（i,r）=（A（i,r）-A（i,1:
r-1）*A（1:
r-1,r））/A（r,r）;
end
ifabs（A（r,r））<=e
error（'该方程不是满秩矩阵'）;
else;
end
end
disp（'经过电算以后得到的A（由U与L中的元素置换而组成）'）
A
fori=1:
H%将A拆成L与U
forj=1:
H
ifi==j
L（i,j）=1;
elseifi>j
L（i,j）=A（i,j）;
else
L（i,j）=0;;
end
end
end
end
L
fori=1:
H%将A拆成L与U
forj=1:
H
ifi<=j
U（i,j）=A（i,j）;
else
U（i,j）=0;;
end
end
end
U
%回带求取解
W=size（L）;%求解Ly=b
Ah=W
（1）;
Al=W
（2）;
Y=zeros（Al,1）;
fori=1:
Ah
S=0;
forj=1:
i-1
S=S+L（i,j）*Y（j,1）;
end
Y（i,1）=（B（i,1）-S）/L（i,i）;
end
Y
W=size（U）;%求解Ux=y
Ah=W
（1）;
Al=W
（2）;
X=zeros（Al,1）;
fori=Ah:
-1:
1
S=0;
forj=Al-1:
-1:
i
S=S+U（i,j+1）*X（j+1,1）;
end
X（i,1）=（Y（i,1）-S）/U（i,i）;
end
X
补充：
SOR算法matlab实现
程序源代码：
e=input（'请输入你要求的精度（用无穷范数衡量）:
'）;
w=input（'请输入松弛因子:
'）;
fid=fopen（'Untitled.m','r'）;
A=fscanf（fid,'%f',[Li,inf]）;
A=A';
W=size（A）;
H=W
（1）;
AZ=A（1:
H,1:
H）;
B=A（1:
H,H+1）;
A=AZ;
fori=1:
H%将A拆成D与L与U
forj=1:
H
ifi==j
D（i,j）=A（i,j）;
else
D（i,j）=0;
end
end
end
D
fori=1:
H%将A拆成D与L与U
forj=1:
H
ifi>j
L（i,j）=-A（i,j）;
else
L（i,j）=0;
end
end
end
L
fori=1:
H%将A拆成L与U
forj=1:
H
ifiU（i,j）=-A（i,j）;
else
U（i,j）=0;;
end
end
end
U
X=zeros（H,1）;
forj=1:
1000
S=X;
fori=1:
H
a=X（i,1）;
X（i,1）=（B（i,1）-A（i,1:
i-1）*X（1:
i-1,1）-A（i,i+1:
H）*X（i+1:
H,1））/A（i,i）;
X（i,1）=a+w*（X（i,1）-a）;
end
disp（X'）
R=S-X;
c=abs（R（1,1））;
fori=1:
H
ifc<=abs（R（i,1））
c=abs（R（i,1））;
else
;
end
end
ifc<=e
break;
else;
end
end
运行结果：
实验体会：
通过实验我掌握了消元法解方程的一些基本算法以及用matlab实现矩阵的几种基本计算。
对MATLAB软件有了更深的了解。
同时，在实验中，在输入向量时错误也很多，比如将：
“（”写成“[”等,致使程序不能运行，无法得到预期的实验结果，有时候还会在结束循环体是落掉“end”等，通过此次练习，我基本上熟悉了Matlab编程的基技能，希望在以后的学习中还能有更多的机会运用它。
升级会员 