高斯消元法解线性方程组.docx
《高斯消元法解线性方程组.docx》由会员分享,可在线阅读,更多相关《高斯消元法解线性方程组.docx(23页珍藏版)》请在冰豆网上搜索。
![高斯消元法解线性方程组.docx](https://file1.bdocx.com/fileroot1/2023-1/10/efabd123-b930-499f-a257-c1380e004aa0/efabd123-b930-499f-a257-c1380e004aa01.gif)
高斯消元法解线性方程组
实验一、高斯消元法解线性方程组
流程图:
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编程的基技能,希望在以后的学习中还能有更多的机会运用它。