数值线性代数第二版徐树方高立张平文上机知识题第一章实验报告.docx
《数值线性代数第二版徐树方高立张平文上机知识题第一章实验报告.docx》由会员分享，可在线阅读，更多相关《数值线性代数第二版徐树方高立张平文上机知识题第一章实验报告.docx（20页珍藏版）》请在冰点文库上搜索。
数值线性代数第二版徐树方高立张平文上机知识题第一章实验报告
上机习题
1.先用你所熟悉的的计算机语言将不选主元和列主元Gauss消去法编写成通用的子程序；然后用你编写的程序求解84阶方程组；最后将你的计算结果与方程的精确解进行比较，并就此谈谈你对Gauss消去法的看法。
Sol：
（1）先用matlab将不选主元和列主元Gauss消去法编写成通用的子程序，得到
：
不选主元Gauss消去法：
得到
满足
列主元Gauss消去法：
得到
满足
（2）用前代法解
，得
用回代法解
，得
求解程序为
（
可缺省，缺省时默认为单位矩阵）
（3）计算脚本为ex1_1
代码
%算法1.1.3（计算三角分解：
Gauss消去法）
function[L,U]=GaussLA（A）
n=length（A）;
fork=1:
n-1
A（k+1:
n,k）=A（k+1:
n,k）/A（k,k）;
A（k+1:
n,k+1:
n）=A（k+1:
n,k+1:
n）-A（k+1:
n,k）*A（k,k+1:
n）;
end
U=triu（A）;
L=tril（A）;
L=L-diag（diag（L））+diag（ones（1,n））;
end
%算法1.2.2（计算列主元三角分解：
列主元Gauss消去法）
function[L,U,P]=GaussCol（A）
n=length（A）;
fork=1:
n-1
[s,t]=max（abs（A（k:
n,k）））;
p=t+k-1;
temp=A（k,1:
n）;
A（k,1:
n）=A（p,1:
n）;
A（p,1:
n）=temp;
u（k）=p;
ifA（k,k）~=0
A（k+1:
n,k）=A（k+1:
n,k）/A（k,k）;
A（k+1:
n,k+1:
n）=A（k+1:
n,k+1:
n）-A（k+1:
n,k）*A（k,k+1:
n）;
else
break;
end
end
L=tril（A）;U=triu（A）;L=L-diag（diag（L））+diag（ones（1,n））;
P=eye（n）;
fori=1:
n-1
temp=P（i,:
）;
P（i,:
）=P（u（i）,:
）;
P（u（i）,:
）=temp;
end
end
%高斯消去法解线性方程组
functionx=Gauss（A,b,L,U,P）
ifnargin<5
P=eye（length（A））;
end
n=length（A）;
b=P*b;
forj=1:
n-1
b（j）=b（j）/L（j,j）;
b（j+1:
n）=b（j+1:
n）-b（j）*L（j+1:
n,j）;
end
b（n）=b（n）/L（n,n）;
y=b;
forj=n:
-1:
2
y（j）=y（j）/U（j,j）;
y（1:
j-1）=y（1:
j-1）-y（j）*U（1:
j-1,j）;
end
y
（1）=y
（1）/U（1,1）;
x=y;
end
ex1_1
clc;clear;
%第一题
A=6*eye（84）+diag（8*ones（1,83）,-1）+diag（ones（1,83）,1）;
b=[7;15*ones（82,1）;14];
%不选主元Gauss消去法
[L,U]=GaussLA（A）;
x1_1=Gauss（A,b,L,U）;
%列主元Gauss消去法
[L,U,P]=GaussCol（A）;
x1_2=Gauss（A,b,L,U,P）;
%解的比较
subplot（1,3,1）;plot（1:
84,x1_1,'o-'）;title（'Gauss'）;
subplot（1,3,2）;plot（1:
84,x1_2,'.-'）;title（'PGauss'）;
subplot（1,3,3）;plot（1:
84,ones（1,84）,'*-'）;title（'精确解'）;
结果为（其中Gauss表示不选主元的Gauss消去法，PGauss表示列主元Gauss消去法，精确解为
）：
由图，显然列主元消去法与精确解更为接近，不选主元的Gauss消去法误差比列主元消去法大，且不如列主元消去法稳定。
Gauss消去法重点在于A的分解过程，无论A如何分解，后面两步的运算过程不变。
2.先用你所熟悉的的计算机语言将平方根法和改进的平方根法编写成通用的子程序；然后用你编写的程序求解对称正定方程组Ax=b。
Sol：
（1）先用matlab将平方根法和改进的平方根法编写成通用的子程序，得到L，（D）：
平方根法：
L=Cholesky（A）
改进的平方根法：
[L,D]=LDLt（A）
（2）求解得
求解得
求解程序为x=Gauss（A,b,L,U,P）（
，
此时缺省，缺省时默认为单位矩阵）
（3）计算脚本为ex1_2
代码
%算法1.3.1（计算Cholesky分解：
平方根法）
functionL=Cholesky（A）
n=length（A）;
fork=1:
n
A（k,k）=sqrt（A（k,k））;
A（k+1:
n,k）=A（k+1:
n,k）/A（k,k）;
forj=k+1:
n
A（j:
n,j）=A（j:
n,j）-A（j:
n,k）*A（j,k）;
end
end
L=tril（A）;
end
%计算LDL‘分解：
改进的平方根法
function[L,D]=LDLt（A）
n=length（A）;
forj=1:
n
fori=1:
n
v（i,1）=A（j,i）*A（i,i）;
end
A（j,j）=A（j,j）-A（j,1:
j-1）*v（1:
j-1,1）;
A（j+1:
n,j）=（A（j+1:
n,j）-A（j+1:
n,1:
j-1）*v（1:
j-1,1））/A（j,j）;
end
L=tril（A）;
D=diag（diag（A））;
L=L-diag（diag（L））+diag（ones（1,n））;
end
%高斯消去法解线性方程组
functionx=Gauss（A,b,L,U,P）
ifnargin<5
P=eye（length（A））;
end
n=length（A）;
b=P*b;
forj=1:
n-1
b（j）=b（j）/L（j,j）;
b（j+1:
n）=b（j+1:
n）-b（j）*L（j+1:
n,j）;
end
b（n）=b（n）/L（n,n）;
y=b;
forj=n:
-1:
2
y（j）=y（j）/U（j,j）;
y（1:
j-1）=y（1:
j-1）-y（j）*U（1:
j-1,j）;
end
y
（1）=y
（1）/U（1,1）;
x=y;
end
ex1_2
%第二题
%第一问
A=10*eye（100）+diag（ones（1,99）,-1）+diag（ones（1,99）,1）;
b=round（100*rand（100,1））;
%平方根法
L=Cholesky（A）;
x1_2_1_1=Gauss（A,b,L,L'）;
%改进的平方根法
[L,D]=LDLt（A）;
x1_2_1_2=Gauss（A,b,L,D*L'）;
%第二问
A=hilb（40）;
b=sum（A）;
b=b';
%平方根法
L=Cholesky（A）;
x1_2_2_1=Gauss（A,b,L,L'）;
%改进的平方根法
[L,D]=LDLt（A）;
x1_2_2_2=Gauss（A,b,L,D*L'）;
结果分别为
x1_2_1_1=
7.2586
8.4143
-0.4013
8.5984
5.4177
0.2249
2.3336
4.4389
8.2772
8.7890
-0.1667
8.8784
8.3824
3.2978
7.6401
0.3014
3.3457
8.2418
6.2368
8.3906
5.8575
-0.9656
7.7981
7.9842
5.3601
6.4143
6.4966
2.6201
6.3024
0.3563
7.1342
-0.6985
2.8506
0.1927
0.2227
7.5801
5.9762
1.6583
9.4409
-1.0677
4.2362
2.7060
6.7037
7.2570
0.7265
4.4778
3.4959
5.5637
5.8675
6.7614
1.5180
6.0582
5.9002
0.9397
0.7031
4.0294
9.0032
1.9382
5.6150
0.9120
7.2652
1.4360
4.3749
5.8146
7.4791
8.3942
4.5789
0.8169
1.2523
1.6603
8.1448
0.8915
7.9401
0.7075
8.9849
2.4437
1.5777
1.7790
5.6319
3.9018
2.3506
7.5925
4.7245
4.1627
8.6483
1.3543
6.8087
6.5589
2.6027
5.4140
0.2577
0.0090
4.6522
6.4685
8.6626
-0.0948
5.2856
4.2385
-0.6706
3.4671
x1_2_1_2=
7.2586
8.4143
-0.4013
8.5984
5.4177
0.2249
2.3336
4.4389
8.2772
8.7890
-0.1667
8.8784
8.3824
3.2978
7.6401
0.3014
3.3457
8.2418
6.2368
8.3906
5.8575
-0.9656
7.7981
7.9842
5.3601
6.4143
6.4966
2.6201
6.3024
0.3563
7.1342
-0.6985
2.8506
0.1927
0.2227
7.5801
5.9762
1.6583
9.4409
-1.0677
4.2362
2.7060
6.7037
7.2570
0.7265
4.4778
3.4959
5.5637
5.8675
6.7614
1.5180
6.0582
5.9002
0.9397
0.7031
4.0294
9.0032
1.9382
5.6150
0.9120
7.2652
1.4360
4.3749
5.8146
7.4791
8.3942
4.5789
0.8169
1.2523
1.6603
8.1448
0.8915
7.9401
0.7075
8.9849
2.4437
1.5777
1.7790
5.6319
3.9018
2.3506
7.5925
4.7245
4.1627
8.6483
1.3543
6.8087
6.5589
2.6027
5.4140
0.2577
0.0090
4.6522
6.4685
8.6626
-0.0948
5.2856
4.2385
-0.6706
3.4671
x1_2_2_1=
1.0e+07*
0.0000
-0.0000
0.0001
-0.0004
-0.0014
0.0424
-0.2980
1.1419
-2.7335
4.2539
-4.3018
2.7733
-1.1989
0.5406
-0.3688
0.3285
-0.4438
0.4621
-0.2513
0.0565
0.0000
-0.0051
0.0071
-0.0027
-0.0031
0.0036
-0.0019
0.0009
0.0002
-0.0002
-0.0006
0.0004
0.0001
-0.0002
0.0001
0.0000
-0.0000
0.0000
-0.0000
-0.0000
x1_2_2_2=
1.0000
1.0000
0.9998
1.0011
1.0064
0.8681
1.8034
-1.5693
5.5763
-2.5315
-1.7693
10.4883
-6.2807
0.5882
-4.7157
22.8299
-19.9134
8.7032
10.3265
-25.2140
10.0282
12.3882
-1.9425
14.1891
-12.0552
-0.5803
-12.4791
8.5652
9.8724
-10.5502
16.3871
-5.8132
13.4216
11.1767
-64.3154
46.3837
12.6957
-21.7556
12.1204
-1.9342
3.用第1题的程序求解第2题的两个方程组并比较所有的计算结果，然后评价各个方法的优劣。
Sol：
Gauss表示不选主元的Gauss消去法，PGauss表示列主元Gauss消去法。
计算脚本为：
%第三题
%第一问
A=10*eye（100）+diag（ones（1,99）,-1）+diag（ones（1,99）,1）;
b=round（100*rand（100,1））;
%不选主元Gauss消去法
[L,U]=GaussLA（A）;
x1_3_1_1=Gauss（A,b,L,U）;
%列主元Gauss消去法
[L,U,P]=GaussCol（A）;
x1_3_1_2=Gauss（A,b,L,U,P）;
%第二问
A=hilb（40）;
b=sum（A）;
b=b';
%不选主元Gauss消去法
[L,U]=GaussLA（A）;
x1_3_2_1=Gauss（A,b,L,U）;
%列主元Gauss消去法
[L,U,P]=GaussCol（A）;
x1_3_2_2=Gauss（A,b,L,U,P）;
ex1_2;
y1=1:
100;y2=1:
40;
subplot（4,2,1）;plot（y1,x1_2_1_1）;title（'平方根法1'）;
subplot（4,2,2）;plot（y1,x1_2_1_2）;title（'改进的平方根法1'）;
subplot（4,2,3）;plot（y1,x1_3_1_1）;title（'Gauss1'）;
subplot（4,2,4）;plot（y1,x1_3_1_2）;title（'PGauss1'）;
subplot（4,2,5）;plot（y2,x1_2_2_1）;title（'平方根法2'）;
subplot（4,2,6）;plot（y2,x1_2_2_2）;title（'改进的平方根法2'）;
subplot（4,2,7）;plot（y2,x1_3_2_1）;title（'Gauss2'）;
subplot（4,2,8）;plot（y2,x1_3_2_2）;title（'PGauss2'）;
平方根法和改进的平法根法计算量更小，计算过程稳定，但使用范围窄；
不选主元和列主元的Gauss消去法计算量较大，但适用范围广。
例题1.3.2
考虑对称正定线性方程组Ax=b，其中向量b是随机生成的，其元素是服从区间[0,1]上均匀分布的随机数，矩阵
这里L是随机生成的一个下三角矩阵，其元素是服从区间[1,2]上均匀分布的随机数。
对n=10，20，...,500分别应用Gauss消去法、列主元Gauss消去法和Cholesky分解法求解该方程组，画出它们所用的CPU时间，其中“Gauss”表示Gauss消去法、“PGauss”表示列主元Gauss消去法,“Cholesky”表示Cholesky分解法。
Sol：
经试验知，对应课本上图1.1所示的Cholesky分解法应为改进后的Cholesky分解法即
分解。
此处所用的CPU时间利用cputime测量。
计算脚本为eg1_3_2
clc;clear;
fori=1:
50;
n=i*10;
b=rand（n,1）;
L=tril（unifrnd（1,2,n,n））;
A=L*L';
t1（i）=cputime;
[L1,U1]=GaussLA（A）;
x1=Gauss（A,b,L1,U1）;
t1（i）=cputime-t1（i）;
t2（i）=cputime;
[L2,U2,P2]=GaussCol（A）;
x2=Gauss（A,b,L2,U2,P2）;
t2（i）=cputime-t2（i）;
t3（i）=cputime;
L3=LDLt（A）;
x3=Gauss（A,b,L3,L3'）;
t3（i）=cputime-t3（i）;
end
N=10:
10:
500;
plot（N,t1,'o-',N,t2,'.-',N,t3,'*-'）;
legend（'Gauss','PGauss','Cholesky'）;
结果为
升级会员 