最新Gauss列主元素消去法实验.docx
《最新Gauss列主元素消去法实验.docx》由会员分享，可在线阅读，更多相关《最新Gauss列主元素消去法实验.docx（32页珍藏版）》请在冰豆网上搜索。
最新Gauss列主元素消去法实验
Lab06．Gauss列主元素消去法实验
【实验目的和要求】
1．使学生深入理解并掌握Gauss消去法和Gauss列主元素消去法步骤；
2．通过对Gauss消去法和Gauss列主元素消去法的程序设计，以提高学生程序设计的能力；
3．对具体问题，分别用Gauss消去法和Gauss列主元素消去法求解。
通过对结果的分析比较，使学生感受Gauss列主元素消去法优点。
【实验内容】
1．根据Matlab语言特点，描述Gauss消去法和Gauss列主元素消去法步骤。
2．编写用不选主元的直接三角分解法解线性方程组Ax=b的M文件。
要求输出Ax=b中矩阵A及向量b，A=LU分解的L与U，detA及解向量x。
3．编写用Gauss列主元素消去法解线性方程组Ax=b的M文件。
要求输出Ax=b中矩阵A及向量b、PA=LU分解的L与U、detA及解向量x，交换顺序。
4．给定方程组
（1）
（2）
先用编写的程序计算，再将
（1）中的系数3.01改为3.00，0.987改为0.990；将
（2）中的系数2.099999改为2.1，5.900001改为9.5,再用Gauss列主元素消去法解，并将两次计算的结果进行比较。
【实验仪器与软件】
1．CPU主频在1GHz以上，内存在128Mb以上的PC；
2．Matlab6.0及以上版本。
实验讲评：
实验成绩：
评阅教师：
200年月日
Lab06．Gauss列主元素消去法实验
第一题：
1、算法描述：
Ⅰ、Gauss消去法由书上定理5可知设Ax=b，其中A∈R
^（n*n）
（1）如果
，则可通过高斯消去法将Ax=b约化为等价的
角形线性方程组，且计算公式为：
1消元计算（k=1,2，….，n-1）
2回带公式
（2）如果A为非奇异矩阵，则可通过高斯消去法将方程组Ax=b约化方程组为上三角矩阵
以上消元和回代过程总的乘除法次数为
，加减法次数为
以上过程就叫高斯消去法。
Ⅱ、Gauss列主元素消去法设Ax=b.本算法用A的具有行交换的列主元素消去法，消元结果冲掉A，乘法
冲掉
，计算解x冲掉常数项b，行列式存放在det中。
1、det
1
2、对于k=1,2，…，n-1
（1）按列主元
（2）如果
，则停止（det（A）=0）
（3）如果
=k则转（4）
换行：
（4）消元计算对于i=k+1，…，n
①
=
/
②对于j=k+1，…，n
③
（5）
3、如果
，则计算停止（det（A）=0）
4、回代求解
（1）
（2）对于i=n-1，…，2,1
5、
第二题：
编写用不选主元的直接三角分解法解线性方程组Ax=b的M文件。
要求输出Ax=b中矩阵A及向量b，A=LU分解的L与U，detA及解向量x。
编写M文件如下：
function[x,L,U]=Doolittle（A,b）
N=size（A）;
n=N
（1）;
L=eye（n,n）;
U=zeros（n,n）;
U（1,1:
n）=A（1,1:
n）;
L（1:
n,1）=A（1:
n,1）/U（1,1）;
fork=2:
n
fori=k:
n
U（k,i）=A（k,i）-L（k,1:
（k-1））*U（1:
（k-1）,i）
end
forj=（k+1）:
n
L（j,k）=（A（j,k）-L（j,1:
（k-1））*U（1:
（k-1）,k））/U（k,k）
end
end
y=SolveDownTriangle（L,b）;
x=SolveUpTriangle（U,y）;
functiony=SolveDownTriangle（L,b）
N=size（L）;
n=N
（1）;
fori=1:
n
if（i>1）
s=L（i,1:
（i-1））*y（1:
（i-1）,1）;
else
s=0;
end
y（i,1）=（b（i）-s）/L（i,i）;
end
functionx=SolveUpTriangle（U,y）
N=size（U）;
n=N
（1）;
fori=n:
-1:
1
if（is=U（i,（i+1）:
n）*x（（1+i）:
n,1）;
else
s=0;
end
x（i,1）=（y（i）-s）/U（i,i）;
end
第三题：
编写用Gauss列主元素消去法解线性方程组Ax=b的M文件。
要求输出Ax=b中矩阵A及向量b、PA=LU分解的L与U、detA及解向量x，交换顺序。
编写M文件如下：
function[x,det,index]=gauss（A,b）
[n,m]=size（A）;
nb=length（b）;
ifn~=m
error（'TherowsandcolumnsofnatrixAmustbeequal!
'）;
return;
end
ifm~=nb
error（'ThecolumnsofAmustbeequalthelengthb!
'）;
return;
end
index=1;det=1;x=zeros（n,1）;
fork=1:
n-1
maxa=0;
fori=k:
n
ifabs（A（i,k））>maxa
maxa=abs（A（i,k））;r=i;
end
end
ifmaxa<1e-10
index=0;return;
end
ifr>k
forj=k:
n
z=A（k,j）;A（k,j）=A（r,j）;A（r,j）=z;
end
z=b（k）;b（k）=b（r）;b（r）=z;det=-det;
end
fori=k+1:
n
m=A（i,k）/A（k,k）;
forj=k+1:
n
A（i,j）=A（i,j）-m*A（k,j）;
end
b（i）=b（i）-m*b（k）;
end
det=det*A（k,k）;
end
det=det*A（n,n）;
ifabs（A（n,n））<1e-10
index=0;return;
end
fork=n:
-1:
1
forj=k+1:
n
b（k）=b（k）-A（k,j）*x（j）;
end
x（k）=b（k）/A（k,k）;
end
第四题：
Ⅰ、运用Gauss列主元素消去法解上述程序计算矩阵：
⑴、计算第一个矩阵
1、计算源程序：
计算程序如下：
function[x,det,index]=gauss（A,b）
A=[3.01,6.03,1.99;1.27,4.16,-1.23;0.987,-4.81,9.34];b=[1;1;1];
[n,m]=size（A）;
nb=length（b）;
ifn~=m
error（'TherowsandcolumnsofnatrixAmustbeequal!
'）;
return;
end
ifm~=nb
error（'ThecolumnsofAmustbeequalthelengthb!
'）;
return;
end
index=1;det=1;x=zeros（n,1）;
fork=1:
n-1
maxa=0;
fori=k:
n
ifabs（A（i,k））>maxa
maxa=abs（A（i,k））;r=i;
end
end
ifmaxa<1e-10
index=0;return;
end
ifr>k
forj=k:
n
z=A（k,j）;A（k,j）=A（r,j）;A（r,j）=z;
end
z=b（k）;b（k）=b（r）;b（r）=z;det=-det;
end
fori=k+1:
n
m=A（i,k）/A（k,k）;
forj=k+1:
n
A（i,j）=A（i,j）-m*A（k,j）;
end
b（i）=b（i）-m*b（k）;
end
det=det*A（k,k）;
end
det=det*A（n,n）;
ifabs（A（n,n））<1e-10
index=0;return;
end
fork=n:
-1:
1
forj=k+1:
n
b（k）=b（k）-A（k,j）*x（j）;
end
x（k）=b（k）/A（k,k）;
end
运行结果如下：
ans=
1.0e+003*
1.5926
-0.6319
-0.4936
>>
②、计算修改后的程序：
计算程序如下：
function[x,det,index]=gauss（A,b）
A=[3.00,6.03,1.99;1.27,4.16,-1.23;0.990,-4.81,9.34];b=[1;1;1];
[n,m]=size（A）;
nb=length（b）;
ifn~=m
error（'TherowsandcolumnsofnatrixAmustbeequal!
'）;
return;
end
ifm~=nb
error（'ThecolumnsofAmustbeequalthelengthb!
'）;
return;
end
index=1;det=1;x=zeros（n,1）;
fork=1:
n-1
maxa=0;
fori=k:
n
ifabs（A（i,k））>maxa
maxa=abs（A（i,k））;r=i;
end
end
ifmaxa<1e-10
index=0;return;
end
ifr>k
forj=k:
n
z=A（k,j）;A（k,j）=A（r,j）;A（r,j）=z;
end
z=b（k）;b（k）=b（r）;b（r）=z;det=-det;
end
fori=k+1:
n
m=A（i,k）/A（k,k）;
forj=k+1:
n
A（i,j）=A（i,j）-m*A（k,j）;
end
b（i）=b（i）-m*b（k）;
end
det=det*A（k,k）;
end
det=det*A（n,n）;
ifabs（A（n,n））<1e-10
index=0;return;
end
fork=n:
-1:
1
forj=k+1:
n
b（k）=b（k）-A（k,j）*x（j）;
end
x（k）=b（k）/A（k,k）;
end
计算结果如下：
ans=
119.5273
-47.1426
-36.8403
>>
⑵：
计算第二个矩阵：
①：
计算源程序：
计算程序如下：
function[x,det,index]=gauss（A,b）
A=[10,-7,0,1;-3,2.09999,6,2;5,-1,5,-1;2,1,0,2];b=[8;5.900001;5;1];
[n,m]=size（A）;
nb=length（b）;
ifn~=m
error（'TherowsandcolumnsofnatrixAmustbeequal!
'）;
return;
end
ifm~=nb
error（'ThecolumnsofAmustbeequalthelengthb!
'）;
return;
end
index=1;det=1;x=zeros（n,1）;
fork=1:
n-1
maxa=0;
fori=k:
n
ifabs（A（i,k））>maxa
maxa=abs（A（i,k））;r=i;
end
end
ifmaxa<1e-10
index=0;return;
end
ifr>k
forj=k:
n
z=A（k,j）;A（k,j）=A（r,j）;A（r,j）=z;
end
z=b（k）;b（k）=b（r）;b（r）=z;det=-det;
end
fori=k+1:
n
m=A（i,k）/A（k,k）;
forj=k+1:
n
A（i,j）=A（i,j）-m*A（k,j）;
end
b（i）=b（i）-m*b（k）;
end
det=det*A（k,k）;
end
det=det*A（n,n）;
ifabs（A（n,n））<1e-10
index=0;return;
end
fork=n:
-1:
1
forj=k+1:
n
b（k）=b（k）-A（k,j）*x（j）;
end
x（k）=b（k）/A（k,k）;
end
计算结果如下：
ans=
0.0000
-1.0000
1.0000
1.0000
>>
②、计算修改后的程序：
计算程序如下：
function[x,det,index]=gauss（A,b）A=[10,-7,0,1;-3,2.1,6,2;5,-1,5,-1;2,1,0,2];b=[8;9.5;5;1];
[n,m]=size（A）;
nb=length（b）;
ifn~=m
error（'TherowsandcolumnsofnatrixAmustbeequal!
'）;
return;
end
ifm~=nb
error（'ThecolumnsofAmustbeequalthelengthb!
'）;
return;
end
index=1;det=1;x=zeros（n,1）;
fork=1:
n-1
maxa=0;
fori=k:
n
ifabs（A（i,k））>maxa
maxa=abs（A（i,k））;r=i;
end
end
ifmaxa<1e-10
index=0;return;
end
ifr>k
forj=k:
n
z=A（k,j）;A（k,j）=A（r,j）;A（r,j）=z;
end
z=b（k）;b（k）=b（r）;b（r）=z;det=-det;
end
fori=k+1:
n
m=A（i,k）/A（k,k）;
forj=k+1:
n
A（i,j）=A（i,j）-m*A（k,j）;
end
b（i）=b（i）-m*b（k）;
end
det=det*A（k,k）;
end
det=det*A（n,n）;
ifabs（A（n,n））<1e-10
index=0;return;
end
fork=n:
-1:
1
forj=k+1:
n
b（k）=b（k）-A（k,j）*x（j）;
end
x（k）=b（k）/A（k,k）;
end
计算结果如下：
ans=
-0.3543
-1.4252
1.3827
1.5669
>>
Ⅱ、运用不选主元的直接三角分解法解上述程：
1、计算第一个矩阵：
1、计算源程序：
计算程序如下：
function[x,L,U]=Doolittle（A,b）
A=[3.01,6.03,1.99;1.27,4.16,-1.23;0.987,-4.81,9.34];b=[1,1,1];
N=size（A）;
n=N
（1）;
L=eye（n,n）;
U=zeros（n,n）;
U（1,1:
n）=A（1,1:
n）;
L（1:
n,1）=A（1:
n,1）/U（1,1）;
fork=2:
n
fori=k:
n
U（k,i）=A（k,i）-L（k,1:
（k-1））*U（1:
（k-1）,i）
end
forj=（k+1）:
n
L（j,k）=（A（j,k）-L（j,1:
（k-1））*U（1:
（k-1）,k））/U（k,k）
end
end
y=SolveDownTriangle（L,b）;
x=SolveUpTriangle（U,y）;
functiony=SolveDownTriangle（L,b）
N=size（L）;
n=N
（1）;
fori=1:
n
if（i>1）
s=L（i,1:
（i-1））*y（1:
（i-1）,1）;
else
s=0;
end
y（i,1）=（b（i）-s）/L（i,i）;
end
functionx=SolveUpTriangle（U,y）
N=size（U）;
n=N
（1）;
fori=n:
-1:
1
if（is=U（i,（i+1）:
n）*x（（1+i）:
n,1）;
else
s=0;
end
x（i,1）=（y（i）-s）/U（i,i）;
end
计算结果如下：
U=
3.01006.03001.9900
01.61580
000
U=
3.01006.03001.9900
01.6158-2.0696
000
L=
1.000000
0.42191.00000
0.3279-4.20061.0000
U=
3.01006.03001.9900
01.6158-2.0696
00-0.0063
ans=
1.0e+003*
1.5926
-0.6319
-0.4936
>>
2、计算修改后的程序：
计算程序如下：
function[x,L,U]=Doolittle（A,b）
A=[3.00,6.03,1.99;1.27,4.16,-1.23;0.990,-4.81,9.34];b=[1,1,1];
N=size（A）;
n=N
（1）;
L=eye（n,n）;
U=zeros（n,n）;
U（1,1:
n）=A（1,1:
n）;
L（1:
n,1）=A（1:
n,1）/U（1,1）;
fork=2:
n
fori=k:
n
U（k,i）=A（k,i）-L（k,1:
（k-1））*U（1:
（k-1）,i）
end
forj=（k+1）:
n
L（j,k）=（A（j,k）-L（j,1:
（k-1））*U（1:
（k-1）,k））/U（k,k）
end
end
y=SolveDownTriangle（L,b）;
x=SolveUpTriangle（U,y）;
functiony=SolveDownTriangle（L,b）
N=size（L）;
n=N
（1）;
fori=1:
n
if（i>1）
s=L（i,1:
（i-1））*y（1:
（i-1）,1）;
else
s=0;
end
y（i,1）=（b（i）-s）/L（i,i）;
end
functionx=SolveUpTriangle（U,y）
N=size（U）;
n=N
（1）;
fori=n:
-1:
1
if（is=U（i,（i+1）:
n）*x（（1+i）:
n,1）;
else
s=0;
end
x（i,1）=（y（i）-s）/U（i,i）;
end
计算结果：
U=
3.00006.03001.9900
01.60730
000
U=
3.00006.03001.9900
01.6073-2.0724
000
L=
1.000000
0.42331.00000
0.3300-4.23061.0000
U=
3.00006.03001.9900
01.6073-2.0724
00-0.0844
ans=
119.5273
-47.1426
-36.8403
>>
⑵、计算第二个矩阵的程序：
①、计算源程序：
计算程序如下：
function[x,L,U]=Doolittle（A,b）
A=[10,-7,0,1;-3,2.09999,6,2;5,-1,5,-1;2,1,0,2];b=[8,5.900001,5,1];
N=size（A）;
n=N
（1）;
L=eye（n,n）;
U=zeros（n,n）;
U（1,1:
n）=A（1,1:
n）;
L（1:
n,1）=A（1:
n,1）/U（1,1）;
fork=2:
n
fori=k:
n
U（k,i）=A（k,i）-L（k,1:
（k-1））*U（1:
（k-1）,i）
end
forj=（k+1）:
n
L（j,k）=（A（j,k）-L（j,1:
（k-1））*U（1:
（k-1）,k））/U（k,k）
end
end
y=SolveDownTriangle（L,b）;
x=SolveUpTriangle（U,y）;
functiony=SolveDownTriangle（L,b）
N=size（L）;
n=N
（1）;
fori=1:
n
if（i>1）
s=L（i,1:
（i-1））*y（1:
（i-1）,1）;
else
s=0;
end
y（i,1）=（b（i）-s）/L（i,i）;
end
functionx=SolveUpTriangle（U,y）
N=size（U）;
n=N
（1）;
fori=n:
-1:
1
if（is=U（i,（i+1）:
n）*x（（1+i）:
n,1）;
else
s=0;
end
x（i,1）=（y（i）-s）/U（i,i）;
end
计算结果如下：
U=
10.0000-7.000001.0000
0-0.000000
0000
0000
U=
10.0000-7.000001.0000
0-0.00006.00000
0000
0000
U=
10.0000-7.000001.0000
0-0.00006.00002.3000
0000
0000
L=
1.0e+005*
0.0000000
-0.00000.000000
0.0000-2.50000.00000
0.0000000.0000
L=
1.0e+005*
0.0000000
-0.00000.000000
0.0000-2.50000.00000
0.0000-2.400000.0000
U=
1.0e+006*
0.0000-0.000000.0000
0-0.00000.00000.0000
001.50000
0000
U=
1.0e+006*
0.0000-0.000000.0000
0
升级会员 