数值分析线性方程组的实验报告包含代码解析.docx
《数值分析线性方程组的实验报告包含代码解析.docx》由会员分享，可在线阅读，更多相关《数值分析线性方程组的实验报告包含代码解析.docx（20页珍藏版）》请在冰点文库上搜索。
数值分析线性方程组的实验报告包含代码解析
线性代数方程组的直接解法
实验目的：
线性方程组求解的直接法编程实现，
实验内容：
线性方程组求解的高斯消去法算法实现
线性方程组求解的主元素消去法算法实现
线性方程组求解的LU分解得方法算法实现
线性方程组求解追赶法算法实现
实验比较：
高斯消去、主元素消去、LU分解都用实例
这个进行比较。
知识理论
解线性方程组的方法大致分为直接法和迭代法。
直接法是指假设计算过程中不产生舍入误差，经过有限次的运算可求的方程组精确解的方法。
方程（2-1）
记为：
AX=b；
一、高斯（Gauss）消元法
（1）.Gauss消元法是最基本的一种方法。
先逐次消去变量，将方程组化成同解的上三角方程组。
消元过程：
先逐次消去变量，将方程组化成同解的上三角方程组
基本思想
回代过程：
按方程的相反顺序求解三角形方程组，得到原方程组的解。
（2）Gauss消去法的求解思路为：
若先让第一个方程组保持不变，利用它消去其余方程组中的,使之变成一个关于变元的n-1阶方程组。
按照
（1）中的思路继续运算得到更为低阶的方程组。
经过n-1步的消元后，得到一个三角方程。
利用求解公式回代得到线性方程组的解。
1
消元过程：
第一次消元，
（i=1,2,3……，n）.将（2-1）中第i个方程减去第一个方程乘以
（i=1,2,3……，n），完成第一次消元，
（2-2）
其中：
；
：
简记为：
其中：
按上述方法完成n-1次消元后得到。
同解的三角方程组：
简记为：
2回代过程：
按逆序逐步回代得到方程的解。
（3）算法：
（5）Matlab程序代码：
function[RA,RB,n,X]=gaus（A,b）
B=[A,b];
n=length（b）;
RA=rank（A）;
RB=rank（B）;
zhicha=RB-RA;
if（zhicha~=0）
disp（'因为RA≠RB，所以此方程组无解'）;
return;
end
ifRA==RB
ifRA==n
disp（'因为RA=RB=n,所以此方程有唯一解'）;
X=zeros（n,1）;
forp=1:
n-1
fork=p+1:
n
m=B（k,p）/B（p,p）;
B（k,p:
n+1）=B（k,p:
n+1）-m*B（p,p:
n+1）;
end
end
b=B（1:
n,n+1）;
A=B（1:
n,1:
n）;
X（n）=b（n）/A（n,n）;
forq=n-1:
-1:
1
X（q）=（b（q）-sum（A（q,q+1:
n）*X（q+1:
n）））/A（q,q）;
end
else
disp（'因为RA=RBend
end
运行检测：
>>A=[111;12-33;-183-1];
b=[6;15;-15];
[RA,RB,n,X]=gaus（A,b）
因为RA=RB=n,所以此方程有唯一解
>>RA=
3
RB=
3
n=
3
X=
1.0000
2.0000
3.0000
二、主元素法
1.列主元素法消元
（1）基本思想：
在每次消元前，在要消去未知数的系数中找到绝对值最大的系数作为主元，通过对换行将其换到对角线上，然后进行消元.
（2）消元过程：
与Gauss很类似，每次对对角线换成最大的值，后面过程与Gauss基本相同。
如此经过n-1步,（2-1）的增广矩阵[A,b]被化为上三角矩阵；
回代过程：
同Gauss算法一样回代求解。
（3）算法：
det<-1
对于k=1,2,3,…n-1;
|aik,k|=max|aik|（k≤i≤n）
（i）如果aik=0,则停止计算（det（A）=0）
（ii）如果aik=k,则zhuan（i）
换行akj<-->aik,j（j=k,k+1,…n）
bk<-->bik
后边就是Gauss消元
（4）.Matlab程序
functionX=liezhu（A,b）
B=[A,b];
n=length（b）;
RA=rank（A）;
RB=rank（B）;
ifRA==RB
ifRA==n
disp（'因为系数矩阵与增广矩阵的秩相同且为n，此方程有唯一解'）
X=zeros（n,1）;
forp=1:
n-1
[y,j]=max（abs（B（p:
n,p）））;
C=B（p,:
）;
B（p,:
）=B（j+p-1,:
）;
B（j+p-1,:
）=C;
end
forp=1:
n-1
fork=p+1:
n
m=B（k,p）/B（p,p）;
B（k,p:
n+1）=B（k,p:
n+1）-m*B（p,p:
n+1）;
end
end
b=B（1:
n,n+1）;
A=B（1:
n,1:
n）;
X（n）=b（n）/A（n,n）;
forq=n-1:
-1:
1
X（q）=（b（q）-sum（A（q,q+1:
n）*X（q+1:
n）））/A（q,q）;
end
else
disp（'因为系数矩阵与增广矩阵的秩相同且小于n，方程组有无穷解'）
end
else
disp（'因为系数矩阵与增广矩阵的秩不相同，所以此方程组无解'）
end
end
测试结果：
>>clearall
clc
A=[111;12-33;-183-1];
b=[615-15]';
liezhu（A,b）
>>因为系数矩阵与增广矩阵的秩相同且为n，此方程有唯一解
ans=
1.0000
2.0000
3.0000
例2测试结果
>>clearall
clc
A=[0.501.13.1;2.04.50.36;5.00.966.5];
b=[6.00.0200.96]';
liezhu（A,b）
>>ans=
-2.6000
1.0000
2.0000
三、直接三角分解法
1.高斯消元的矩阵表示
第一次消元
经过K次消元后得到：
消元过程：
则
（1）LU分解：
，
因为A的各元素已知，可求出L和U的各元素。
（2）看A的第k行，有
，再看A的第k列，有
。
这样就可以计算计算出U的第k行和L的第k列的元素。
（3）从k=1到k=n交替计算U和L，就能逐步计算出U和L的全部元素。
（4）分两步求解LUx=b，先解方程组Ly=b，因为L是下三角矩阵，求解只要逐步向前回代。
第二步是解方程组Ux=y，因为U是上三角矩阵，用逐次向后回代的方法求出x。
3.Matlab程序
function[L,U,Y,X]=Doolittle（A,b）
n=length（b）;
L=zeros（n,n）;
U=zeros（n,n）;
%先将A进行LU分解
fori=1:
n
%把第一列第一行求出
forj=i:
n
ifi==1
A（i,j）=A（i,j）;
ifj>i
A（j,i）=A（j,i）/A（i,i）;
end
end
end
%下面计算U的第i行值
forj=i:
n
ifi~=1
m=0;
fork=1:
i-1
m=m+A（i,k）*A（k,j）;
end
A（i,j）=A（i,j）-m;
end
end
%下面计算L的第i列值
forj=i:
n
ifi~=1&&j>i
l=0;
fork=1:
i-1
l=l+A（j,k）*A（k,i）;
end
A（j,i）=（A（j,i）-l）/A（i,i）;
end
end
end
%将LU的值分别加入
fori=1:
n
forj=i:
n
U（i,j）=A（i,j）;
end
forj=1:
i
ifi==j
L（i,j）=1;
else
L（i,j）=A（i,j）;
end
end
end
%LY=B，UX=Y，进行求解
Ab=[L,b]
Y=zeros（n,1）;
Y
（1）=Ab（1,n+1）/Ab（1,1）;
fori=2:
n
forj=1:
i-1
Y（i）=Y（i）+Ab（i,j）*Y（j）;
end
Y（i）=（Ab（i,n+1）-Y（i））/Ab（i,i）;
end
Ab=[U,Y];
X=zeros（n,1）;
X（n）=Ab（n,n+1）/Ab（n,n）;
fori=n-1:
-1:
1
forj=i+1:
n
X（i）=X（i）+Ab（i,j）*X（j）;
end
X（i）=（Ab（i,n+1）-X（i））/Ab（i,i）;
end
end
运行检测：
>>A=[111;12-33;-183-1];
b=[6;15;-15];
[L,U,Y,X]=Doolittle（A,b）
运行结果：
>>Ab=
1.0000006.0000
12.00001.0000015.0000
-18.0000-1.40001.0000-15.0000
L=
1.000000
12.00001.00000
-18.0000-1.40001.0000
U=
1.00001.00001.0000
0-15.0000-9.0000
004.4000
Y=
6.0000
-57.0000
13.2000
X=
1.0000
2.0000
3.0000
四、解三角方程组的追赶法
1.设系数矩阵为三对角矩阵
则方程组Ax=f称为三对角方程组。
2.设矩阵A非奇异，A有Crout分解A=LU，其中L为下三角矩阵，U为单位上三角矩阵，记
先依次求出L，U中的元素后，令Ux=y，先求解下三角方程组Ly=f得出y，再求解上三角方程组Ux=y。
求解三对角方程组的2追赶法将矩阵三角分解的计算与求解两个三角方程组的计算放在一起，使算法吏为紧凑。
其计算公式为：
3.算法：
1.输入a=（a2,…，an）,b=（b1,…，bn）,c=（c1,…，cn-1）,d=（d1,…，dn）,n
2.对i=2,3，…，n
ai/bi-1=>ai（li）
bi-ci-1ai=>bi（ui）
di-aidi-1=>di（yi）
3.dn/bn=>dn（xn）
4.对i=n-1,n-2…，1
（di-cidi+1）/bi=>di（xi）
5.输出d=x,停机
4.Matlab程序：
clc;
clear;
%如果需要验证其他矩阵的解，只需更改A、B和n的值即可
A=[2,1,0,0;1,3,1,0;0,1,4,-1;0,0,2,3]
B=[3;5;4;5]
n=4;
%获得A中不为0的数存到相应数组中
fori=1:
n
forj=1:
n
ifi==j
b（i）=A（i,j）;
end
ifi==j-1
c（j-1）=A（i,j）;
end
ifj==i-1
a（i）=A（i,j）;
end
end
end
%将变换后的值存入相应矩阵中
u
（1）=b
（1）;
fori=2:
n;
l（i）=a（i）/u（i-1）;
u（i）=b（i）-l（i）*c（i-1）;
end
fori=1:
n
forj=1:
n
ifi==j
L（i,j）=1;
U（i,j）=u（i）;
end
ifi==j-1
U（i,j）=c（i）;
end
ifj==i-1
L（i,j）=l（i）;
end
end
end
%LY=B，UX=Y，进行求解
disp（'增广矩阵'）;
AB=[L,B]
B=zeros（n,1）;
B
（1）=AB（1,n+1）/AB（1,1）;
fori=2:
n
forj=1:
i-1
B（i）=B（i）+AB（i,j）*B（j）;
end
B（i）=（AB（i,n+1）-B（i））/AB（i,i）;
end
disp（'解方程LY=B得'）;
Y=B
AB=[U,B];
B=zeros（n,1）;
B（n）=AB（n,n+1）/AB（n,n）;
fori=n-1:
-1:
1
forj=i+1:
n
B（i）=B（i）+AB（i,j）*B（j）;
end
B（i）=（AB（i,n+1）-B（i））/AB（i,i）;
end
disp（'方程的解为：
'）;
X=B
运行结果：
A=
2100
1310
014-1
0023
B=
3
5
4
5
增广矩阵
AB=
1.00000003.0000
0.50001.0000005.0000
00.40001.000004.0000
000.55561.00005.0000
解方程LY=B得
Y=
3.0000
3.5000
2.6000
3.5556
方程的解为：
X=
1.0000
1.0000
1.0000
1.0000
五、算法比较：
Gauss消去法是针对一般的线性方程组，与线性代数中的初等变换解线性方程组方法类似。
高斯消元法的算法复杂度是O（n3）；这就是说，如果系数矩阵的是n×n，那么高斯消元法所需要的计算量大约与n3成比例。
高斯消元法可用在任何域中。
高斯消元法对于一些矩阵来说是稳定的。
对于普遍的矩阵来说，高斯消元法通常也是稳定的。
LU分解法比较简便迅速，当解多个系数矩阵为A的线性方程组时，LU分解法就显得特别优越，只要对系数矩阵做一次LU分解，以后只要解三角形方程即可。
也可以更具系数矩阵的形状来设计算法。
追赶法只是针对系数矩阵为三对角阵的方程组，因此是一种特殊的方程组。
此方法效率较高，不过不适用于一般的线性方程组。
升级会员 