数值分析第07次作业.docx
《数值分析第07次作业.docx》由会员分享，可在线阅读，更多相关《数值分析第07次作业.docx（23页珍藏版）》请在冰豆网上搜索。
数值分析第07次作业
实验七直接解线性方程组的数值实验
实验1：
题目：
用列主元Gauss消元法求解以下方程组：
编程及结果：
编程：
fprintf（'增广矩阵A'）
A=[0.7290.810.90.6867;1110.8338;1.3311.211.11.000]
tempo=A（3,:
）;A（3,:
）=A（1,:
）;A（1,:
）=tempo;
fprintf（'第一次选主元后的矩阵'）
A
A（2,:
）=A（2,:
）-A（1,:
）*A（2,1）/A（1,1）;
A（3,:
）=A（3,:
）-A（1,:
）*A（3,1）/A（1,1）;
fprintf（'第一次消元后的矩阵'）
A
tempo=A（3,:
）;A（3,:
）=A（2,:
）;A（2,:
）=tempo;
fprintf（'第二次选主元后的矩阵'）
A
A（3,:
）=A（3,:
）-A（2,:
）*A（3,2）/A（2,2）;
fprintf（'第二次消元后的矩阵'）
A
%回代求解方程组
x（3）=A（3,4）/A（3,3）;
x
（2）=（A（2,4）-A（2,3）*x（3））/A（2,2）;
x
（1）=（A（1,4）-A（1,3）*x（3）-A（1,2）*x
（2））/A（1,1）;
x
结果：
增广矩阵A
A=
Columns1through2
0.72900.8100
1.00001.0000
1.33101.2100
Columns3through4
0.90000.6867
1.00000.8338
1.10001.0000
第一次选主元后的矩阵
A=
Columns1through2
1.33101.2100
1.00001.0000
0.72900.8100
Columns3through4
1.10001.0000
1.00000.8338
0.90000.6867
第一次消元后的矩阵
A=
Columns1through2
1.33101.2100
00.0909
00.1473
Columns3through4
1.10001.0000
0.17360.0825
0.29750.1390
第二次选主元后的矩阵
A=
Columns1through2
1.33101.2100
00.1473
00.0909
Columns3through4
1.10001.0000
0.29750.1390
0.17360.0825
第二次消元后的矩阵
A=
Columns1through2
1.33101.2100
00.1473
00
Columns3through4
1.10001.0000
0.29750.1390
-0.0101-0.0033
x=
Columns1through2
0.22450.2814
Column3
0.3279
实验2：
题目：
写用追赶法解三对角方程组的程序，并解下列方程组：
编程及结果：
编程：
A=[2,-1,0,0;-1,2,-1,0;0,-1,2,-1;0,0,-1,2];
b=[5,-12,11,-1]';
[l,u,p]=lu（A）;%对矩阵A进行LU分解,l单位下三角，u上三角，p是置换矩阵
%l,u
y
（1）=b
（1）;
y
（2）=b
（2）-l（2,1）*y
（1）;
y（3）=b（3）-l（3,1）*y
（1）-l（3,2）*y
（2）;
y（4）=b（4）-l（4,1）*y
（1）-l（4,2）*y
（2）-l（4,3）*y（3）;
x（4）=y（4）/u（4,4）;
x（3）=（y（3）-u（3,4）*x（4））/u（3,3）;
x
（2）=（y
（2）-u（2,4）*x（4）-u（2,3）*x（3））/u（2,2）;
x
（1）=（y
（1）-u（1,4）*x（4）-u（1,3）*x（3）-u（1,2）*x
（2））/u（1,1）;
x
结果：
x=
1.0000-3.00005.00002.0000
实验3：
题目：
用LU分解及列主元高斯消去法解线性方程组
输出
中系数
分解的矩阵
和
，解向量
及
；列主元法的行交换次序，解向量
及
；比较两种方法所得的结果。
编程及结果：
1）LU分解
clc
clear
A=[10-701
-32.09999962
5-15-1
2102];
b=[85.90000151];
[l,u,p]=lu（A）;%对矩阵A进行LU分解,l单位下三角，u上三角，p是置换矩阵
l,u
y
（1）=b
（1）;
y
（2）=b
（2）-l（2,1）*y
（1）;
y（3）=b（3）-l（3,1）*y
（1）-l（3,2）*y
（2）;
y（4）=b（4）-l（4,1）*y
（1）-l（4,2）*y
（2）-l（4,3）*y（3）;
x（4）=y（4）/u（4,4）;
x（3）=（y（3）-u（3,4）*x（4））/u（3,3）;
x
（2）=（y
（2）-u（2,4）*x（4）-u（2,3）*x（3））/u（2,2）;
x
（1）=（y
（1）-u（1,4）*x（4）-u（1,3）*x（3）-u（1,2）*x
（2））/u（1,1）;
x
t=det（A）
结果：
l=
1.0000000
0.50001.000000
-0.3000-0.00001.00000
0.20000.9600-0.80001.0000
u=
10.0000-7.000001.0000
02.50005.0000-1.5000
006.00002.3000
0005.0800
x=0.1949-0.76610.96950.6882
t=
-762.0001
2）列主元高斯消去法
clear
clc
A=[10-701;-32.09999962;5-15-1;2102];
b=[85.90000151];
n=length（b）;
X=zeros（n,1）;
c=zeros（1,n）;
d1=0;
t=det（A）
fori=1:
n-1
max=abs（A（i,i））;
m=i;
forj=i+1:
n
ifmaxmax=abs（A（j,i））;
m=j;
end
end
if（m~=i）
fork=i:
n
c（k）=A（i,k）;
A（i,k）=A（m,k）;
A（m,k）=c（k）;
end
d1=b（i）;
b（i）=b（m）;
b（m）=d1;
end
fork=i+1:
n
forj=i+1:
n
A（k,j）=A（k,j）-A（i,j）*A（k,i）/A（i,i）;
end
b（k）=b（k）-b（i）*A（k,i）/A（i,i）;
A（k,i）=0;
end
A
end
X（n）=b（n）/A（n,n）;
fori=n-1:
-1:
1
sum=0;
forj=i+1:
n
sum=sum+A（i,j）*X（j）;
end
X（i）=（b（i）-sum）/A（i,i）;
end
X
>>
t=
762.0001
A=
10.0000-7.000001.0000
0-0.00006.00002.3000
02.50005.0000-1.5000
02.400001.8000
A=
10.0000-7.000001.0000
02.50005.0000-1.5000
006.00002.3000
00-4.80003.2400
A=
10.0000-7.000001.0000
02.50005.0000-1.5000
006.00002.3000
0005.0800
X=
0.0000
-1.0000
1.0000
1.0000
t=
762.0001
实验4：
题目：
用列主元高斯消去法解线性方程组
：
（1）
（2）
分别输出
，解向量
，
（1）中
的条件数，分析比较
（1）,
（2）的计算结果。
编程及结果：
（1）
编写程序：
clear
clc
A=[3.016.031.99;1.274.16-1.23;0.987-4.819.34]
b=[111]
n=length（b）;
X=zeros（n,1）;
c=zeros（1,n）;
d1=0;
fprintf（'输出det（A）'）
det（A）
fprintf（'输出条件数：
1-条件数，2-条件数，无穷-条件数'）
[cond（A,1）,cond（A）,cond（A,inf）]
fori=1:
n-1
max=abs（A（i,i））;
m=i;
forj=i+1:
n
ifmaxmax=abs（A（j,i））;
m=j;
end
end
if（m~=i）
fork=i:
n
c（k）=A（i,k）;
A（i,k）=A（m,k）;
A（m,k）=c（k）;
end
d1=b（i）;
b（i）=b（m）;
b（m）=d1;
end
fork=i+1:
n
forj=i+1:
n
A（k,j）=A（k,j）-A（i,j）*A（k,i）/A（i,i）;
end
b（k）=b（k）-b（i）*A（k,i）/A（i,i）;
A（k,i）=0;
end
end
X（n）=b（n）/A（n,n）;
fori=n-1:
-1:
1
sum=0;
forj=i+1:
n
sum=sum+A（i,j）*X（j）;
end
X（i）=（b（i）-sum）/A（i,i）;
end
X
运行结果：
A=
3.01006.03001.9900
1.27004.1600-1.2300
0.9870-4.81009.3400
b=
111
输出det（A）
ans=
-0.0305
输出条件数：
1-条件数，2-条件数，无穷-条件数
ans=
1.0e+004*
5.52283.06975.6751
X=
1.0e+003*
1.5926
-0.6319
-0.4936
（2）
编写程序：
clear
clc
A=[3.006.031.99;1.274.16-1.23;0.990-4819.34]
b=[111]
n=length（b）;
X=zeros（n,1）;
c=zeros（1,n）;
d1=0;
fprintf（'输出det（A）'）
det（A）
fprintf（'输出条件数：
1-条件数，2-条件数，无穷-条件数'）
[cond（A,1）,cond（A）,cond（A,inf）]
fori=1:
n-1
max=abs（A（i,i））;
m=i;
forj=i+1:
n
ifmaxmax=abs（A（j,i））;
m=j;
end
end
if（m~=i）
fork=i:
n
c（k）=A（i,k）;
A（i,k）=A（m,k）;
A（m,k）=c（k）;
end
d1=b（i）;
b（i）=b（m）;
b（m）=d1;
end
fork=i+1:
n
forj=i+1:
n
A（k,j）=A（k,j）-A（i,j）*A（k,i）/A（i,i）;
end
b（k）=b（k）-b（i）*A（k,i）/A（i,i）;
A（k,i）=0;
end
end
X（n）=b（n）/A（n,n）;
fori=n-1:
-1:
1
sum=0;
forj=i+1:
n
sum=sum+A（i,j）*X（j）;
end
X（i）=（b（i）-sum）/A（i,i）;
end
X
运行结果：
A=
3.00006.03001.9900
1.27004.1600-1.2300
0.9900-481.00009.3400
b=
111
输出det（A）
ans=
-2.9610e+003
输出条件数：
1-条件数，2-条件数，无穷-条件数
ans=
412.8103289.3865343.2779
X=
0.5343
-0.0065
-0.2833
分别输出
，解向量
，
（1）中
的条件数，分析比较
（1）,
（2）的计算结果。
提高题:
题目：
用高斯消元法、列主元的高斯消元法和追赶法分别求解三对角线性方程组：
其中精确解为
编程及结果：
1）高斯消元法
编写程序：
clear
clc
A=[4-100000000
-14-10000000
0-14-1000000
00-14-100000
000-14-10000
0000-14-1000
00000-14-100
000000-14-10
0000000-14-1
00000000-14];
b=[75-1326-1214-45-5];
n=length（b）;
x=b;
fori=1:
n-1
forj=i+1:
n
mi=A（j,i）/A（i,i）;
b（j）=b（j）-mi*b（i）;
fork=i:
n
A（j,k）=A（j,k）-mi*A（i,k）;
end
[Ab'];
end
end
x（n）=b（n）/A（n,n）;
fori=n-1:
-1:
1
s=0;
forj=i+1:
n
s=s+A（i,j）*x（j）;
end
x（i）=（b（i）-s）/A（i,i）;
end
x
运行结果：
x=
Columns1through6
2.00001.0000-3.00000.00001.0000-2.0000
Columns7through10
3.0000-0.00001.0000-1.0000
2）列主元的高斯消元法
编写程序：
clear
clc
A=[4-100000000
-14-10000000
0-14-1000000
00-14-100000
000-14-10000
0000-14-1000
00000-14-100
000000-14-10
0000000-14-1
00000000-14];
b=[75-1326-1214-45-5];
n=length（b）;
X=zeros（n,1）;
c=zeros（1,n）;
d1=0;
fori=1:
n-1
max=abs（A（i,i））;
m=i;
forj=i+1:
n
ifmaxmax=abs（A（j,i））;
m=j;
end
end
if（m~=i）
fork=i:
n
c（k）=A（i,k）;
A（i,k）=A（m,k）;
A（m,k）=c（k）;
end
d1=b（i）;
b（i）=b（m）;
b（m）=d1;
end
fork=i+1:
n
forj=i+1:
n
A（k,j）=A（k,j）-A（i,j）*A（k,i）/A（i,i）;
end
b（k）=b（k）-b（i）*A（k,i）/A（i,i）;
A（k,i）=0;
end
end
X（n）=b（n）/A（n,n）;
fori=n-1:
-1:
1
sum=0;
forj=i+1:
n
sum=sum+A（i,j）*X（j）;
end
X（i）=（b（i）-sum）/A（i,i）;
end
X
运行结果：
>>
X=
2.0000
1.0000
-3.0000
0.0000
1.0000
-2.0000
3.0000
-0.0000
1.0000
-1.0000
3）追赶法
编写程序：
clear
clc
A=[4-100000000
-14-10000000
0-14-1000000
00-14-100000
000-14-10000
0000-14-1000
00000-14-100
000000-14-10
0000000-14-1
00000000-14];
d=[75-1326-1214-45-5];%定义三对角矩阵A的各组成单元。
方程为Ax=d
[l,u,p]=lu（A）;
a=[-1-1-1-1-1-1-1-1-1-1];%a为-1对角线元素（2~n），
b=[4444444444];%b为A的对角线元素（1~n）
c=[-1-1-1-1-1-1-1-1-1-1];%c为+1对角线元素（1~n-1）。
n=length（b）;
u0=0;y0=0;a
（1）=-1;
L
（1）=b
（1）-a
（1）*u0;
y
（1）=（d
（1）-y0*a
（1））/L
（1）;
u
（1）=c
（1）/L
（1）;
fori=2:
（n-1）
L（i）=b（i）-a（i）*u（i-1）;
y（i）=（d（i）-y（i-1）*a（i））/L（i）;
u（i）=c（i）/L（i）;
end
L（n）=b（n）-a（n）*u（n-1）;
y（n）=（d（n）-y（n-1）*a（n））/L（n）;
x（n）=y（n）;
fori=（n-1）:
-1:
1
x（i）=y（i）-u（i）*x（i+1）;
end
x
运行结果：
>>
x=
Columns1through6
2.00001.0000-3.00000.00001.0000-2.0000
Columns7through10
3.0000-0.00001.0000-1.0000
升级会员 