数值线性代数第二版徐树方高立张平文上机习题第一章实验报告材料Word文档格式.docx

1、s,t=max（abs（A（k:n,k）;p=t+k-1;temp=A（k,1:A（k,1:n）=A（p,1:A（p,1:n）=temp;u（k）=p;if A（k,k）=0else break;P=eye（n）;for i=1: temp=P（i,:）; P（i,:）=P（u（i）,: P（u（i）,:）=temp;%高斯消去法解线性方程组function x=Gauss（A,b,L,U,P）if nargin5 P=eye（length（A）;b=P*b;for j=1: b（j）=b（j）/L（j,j）; b（j+1:n）=b（j+1:n）-b（j）*L（j+1:n,j）;b（n）=b（

2、n）/L（n,n）;y=b;for j=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）;y（1）=y（1）/U（1,1）;x=y;ex1_1clc;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）;%解的比较su

3、bplot（1,3,1）;plot（1:84,x1_1,o-title（Gausssubplot（1,3,2）;84,x1_2,.-PGausssubplot（1,3,3）;84,ones（1,84）,*-精确解结果为（其中Gauss表示不选主元的Gauss消去法，PGauss表示列主元Gauss消去法，精确解为）：由图，显然列主元消去法与精确解更为接近，不选主元的Gauss消去法误差比列主元消去法大，且不如列主元消去法稳定。Gauss消去法重点在于A的分解过程，无论A如何分解，后面两步的运算过程不变。2.先用你所熟悉的的计算机语言将平方根法和改进的平方根法编写成通用的子程序；然后用你编写的程

4、序求解对称正定方程组Ax=b。（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分解：平方根法）function L=Cholesky（A）n A（k,k）=sqrt（A（k,k）; for j=k+1: A（j:n,j）=A（j:n,j）-A（j:n,k）*A（j,k）; end%计算LDL分解：改进的平方根法funct

5、ion L,D=LDLt（A） for i=1: v（i,1）=A（j,i）*A（i,i）; 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,1）/A（j,j）;D=diag（diag（A）;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_

6、2_1_2=Gauss（A,b,L,D*L%第二问A=hilb（40）;b=sum（A）;b=b;x1_2_2_1=Gauss（A,b,L,Lx1_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.49

7、66 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

8、.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.67063.4671x1_2_1_2 =x1_2_2_1 = 1.0e+07 * 0.0000 -0.0000 0.0001 -0.0004 -0.0014 0.0424 -0.2980

9、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.0051 0.0071 -0.0027 -0.0031 0.0036 -0.0019 0.0009 0.0002 -0.0002 -0.0006 0.0004x1_2_2_2 = 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.829

10、9 -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.93423.用第1题的程序求解第2题的两个方程组并比较所有的计算结果，然后评价各个方法的优劣。Gauss表示不选主元的Gauss消去法，PGauss表示列主元Gauss消去法。计算脚本为：%第三题x1_3_1_1=

11、Gauss（A,b,L,U）;x1_3_1_2=Gauss（A,b,L,U,P）;x1_3_2_1=Gauss（A,b,L,U）;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）;平方根法1subplot（4,2,2）;plot（y1,x1_2_1_2）;改进的平方根法1subplot（4,2,3）;plot（y1,x1_3_1_1）;Gauss1subplot（4,2,4）;plot（y1,x1_3_1_2）;PGauss1subplot（4,2,5）;plot（y2,x1_2_

12、2_1）;平方根法2subplot（4,2,6）;plot（y2,x1_2_2_2）;改进的平方根法2subplot（4,2,7）;plot（y2,x1_3_2_1）;Gauss2subplot（4,2,8）;plot（y2,x1_3_2_2）;PGauss2平方根法和改进的平法根法计算量更小，计算过程稳定，但使用围窄；不选主元和列主元的Gauss消去法计算量较大，但适用围广。例题1.3.2考虑对称正定线性方程组Ax=b，其中向量b是随机生成的，其元素是服从区间0,1上均匀分布的随机数，矩阵,这里L是随机生成的一个下三角矩阵，其元素是服从区间1,2上均匀分布的随机数。对n=10，20，.,50

13、0分别应用Gauss消去法、列主元Gauss消去法和Cholesky分解法求解该方程组，画出它们所用的CPU时间，其中“Gauss”表示Gauss消去法、“PGauss”表示列主元Gauss消去法,“Cholesky”表示Cholesky分解法。经试验知，对应课本上图1.1所示的Cholesky分解法应为改进后的Cholesky分解法即分解。此处所用的CPU时间利用cputime测量。计算脚本为eg1_3_250; 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）;N=10:10:500;plot（N,t1,N,t2,N,t3,legend（,Cholesky结果为