数值线性代数第二版徐树方高立张平文上机习题第一章实验报告.pdf

《数值线性代数第二版徐树方高立张平文上机习题第一章实验报告.pdf》由会员分享，可在线阅读，更多相关《数值线性代数第二版徐树方高立张平文上机习题第一章实验报告.pdf（17页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、数数值值线线性性代代数数第第二二版版徐徐树树方方高高立立张张平平文文上上机机习习题题第第一一章章实实验验报报告告-CAL-FENGHAI.-(YICAI)-Company One1上机习题1.先用你所熟悉的的计算机语言将不选主元和列主元 Gauss 消去法编写成通用的子程序；然后用你编写的程序求解 84 阶方程组；最后将你的计算结果与方程的精确解进行比较，并就此谈谈你对 Gauss 消去法的看法。Sol：（1）先用 matlab 将不选主元和列主元 Gauss 消去法编写成通用的子程序，得到L,U,P：不选主元 Gauss 消去法：L,U GaussLA(A)得到L,U满足A LU列主元 Ga

2、uss 消去法：L,U,P GaussCol(A)得到L,U,P满足PA LU（2）用前代法解Ly bor Pb，得y用回代法解Ux y，得x求解程序为x GaussA,b,L,U,P（P可缺省，缺省时默认为单位矩阵）（3）计算脚本为 ex1_1代码代码%算法算法 1.1.31.1.3（计算三角分解：（计算三角分解：GaussGauss 消去法）消去法）functionfunction L,U=GaussLA(A)L,U=GaussLA(A)n=length(A);n=length(A);forfor k=1:n-1 k=1:n-1 A(k+1:n,k)=A(k+1:n,k)/A(k,k);A

3、(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);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);endendU=triu(A);U=triu(A);L=tril(A);L=tril(A);L=L-diag(diag(L)+diag(ones(1,n);L=L-diag(diag(L)+diag(ones(1,n);endend%算法算法 1.2.2(1.2.2(计算列主元三角分解：列主元计算列主元三角分解：列主元 GaussGa

4、uss 消去法）消去法）functionfunction L,U,P=GaussCol(A)L,U,P=GaussCol(A)n=length(A);n=length(A);forfor k=1:n-1 k=1:n-1s,t=max(abs(A(k:n,k);s,t=max(abs(A(k:n,k);p=t+k-1;p=t+k-1;temp=A(k,1:n);temp=A(k,1:n);A(k,1:n)=A(p,1:n);A(k,1:n)=A(p,1:n);A(p,1:n)=temp;A(p,1:n)=temp;2u(k)=p;u(k)=p;ifif A(k,k)=0 A(k,k)=0 A(k

5、+1:n,k)=A(k+1:n,k)/A(k,k);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);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);elseelsebreakbreak;endendendendL=tril(A);U=triu(A);L=L-diag(diag(L)+diag(ones(1,n);L=tril(A);U=triu(A);L=L-diag(diag(L)+diag(ones(1,n);P=eye(

6、n);P=eye(n);forfor i=1:n-1 i=1:n-1 temp=P(i,:);temp=P(i,:);P(i,:)=P(u(i),:);P(i,:)=P(u(i),:);P(u(i),:)=temp;P(u(i),:)=temp;endendendend%高斯消去法解线性方程组高斯消去法解线性方程组functionfunction x=Gauss(A,b,L,U,P)x=Gauss(A,b,L,U,P)ifif nargin5 nargin5 P=eye(length(A);P=eye(length(A);endendn=length(A);n=length(A);b=P*b;

7、b=P*b;forfor j=1:n-1 j=1:n-1 b(j)=b(j)/L(j,j);b(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);endendb(n)=b(n)/L(n,n);b(n)=b(n)/L(n,n);y=b;y=b;forfor j=n:-1:2 j=n:-1:2 y(j)=y(j)/U(j,j);y(j)=y(j)/U(j,j);y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);y(1:j-1)=y(1:j-1)-y(j)*U(1:j

8、-1,j);endendy(1)=y(1)/U(1,1);y(1)=y(1)/U(1,1);x=y;x=y;endendex1_1ex1_1clc;clear;clc;clear;%第一题第一题A=6*eye(84)+diag(8*ones(1,83),-1)+diag(ones(1,83),1);A=6*eye(84)+diag(8*ones(1,83),-1)+diag(ones(1,83),1);3b=7;15*ones(82,1);14;b=7;15*ones(82,1);14;%不选主元不选主元 GaussGauss 消去法消去法L,U=GaussLA(A);L,U=GaussLA(

9、A);x1_1=Gauss(A,b,L,U);x1_1=Gauss(A,b,L,U);%列主元列主元 GaussGauss 消去法消去法L,U,P=GaussCol(A);L,U,P=GaussCol(A);x1_2=Gauss(A,b,L,U,P);x1_2=Gauss(A,b,L,U,P);%解的比较解的比较subplot(1,3,1);plot(1:84,x1_1,subplot(1,3,1);plot(1:84,x1_1,o-o-);title();title(GaussGauss););subplot(1,3,2);plot(1:84,x1_2,subplot(1,3,2);plot

10、(1:84,x1_2,.-.-);title();title(PGaussPGauss););subplot(1,3,3);plot(1:84,ones(1,84),subplot(1,3,3);plot(1:84,ones(1,84),*-*-);title();title(精确解精确解););结果为（其中 Gauss表示不选主元的 Gauss消去法，PGauss表示列主元 Gauss消去法，精确解为1,1184）：6x 108Gauss1PGauss21.8精 确 解411.611.41.2110.80.60.410.2210-21-4-605010010501000050100由图，显然

11、列主元消去法与精确解更为接近，不选主元的 Gauss消去法误差比列主元消去法大，且不如列主元消去法稳定。Gauss消去法重点在于 A的分解过程，无论 A 如何分解，后面两步的运算过程不变。2.先用你所熟悉的的计算机语言将平方根法和改进的平方根法编写成通用的子程序；然后用你编写的程序求解对称正定方程组 Ax=b。4Sol：（1）先用 matlab将平方根法和改进的平方根法编写成通用的子程序，得到 L，(D)：平方根法：L=Cholesky(A)改进的平方根法：L,D=LDLt(A)（2）求解得Ly b求解得LTx y or DLTx y求解程序为 x=Gauss(A,b,L,U,P)（U LTo

12、r U DLT，P此时缺省，缺省时默认为单位矩阵）（3）计算脚本为 ex1_2代码%算法算法 1.3.11.3.1（计算（计算 CholeskyCholesky 分解：平方根法）分解：平方根法）functionfunction L=Cholesky(A)L=Cholesky(A)n=length(A);n=length(A);forfor k=1:n k=1:n A(k,k)=sqrt(A(k,k);A(k,k)=sqrt(A(k,k);A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k)=A(k+1:n,k)/A(k,k);forfor j=k+1:n j=k+1:

13、n A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);endendendendL=tril(A);L=tril(A);endend%计算计算 LDLLDL分解：改进的平方根法分解：改进的平方根法functionfunction L,D=LDLt(A)L,D=LDLt(A)n=length(A);n=length(A);forfor j=1:n j=1:nforfor i=1:n i=1:n v(i,1)=A(j,i)*A(i,i);v(i,1)=A(j,i)*A(i,i);endend A(j,j)=A(j,

14、j)-A(j,1:j-1)*v(1:j-1,1);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);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);5endendL=tril(A);L=tril(A);D=diag(diag(A);D=diag(diag(A);L=L-diag(diag(L)+diag(ones(1,n);L=L-diag(diag(L)+diag(ones(1,n);endend

15、%高斯消去法解线性方程组高斯消去法解线性方程组functionfunction x=Gauss(A,b,L,U,P)x=Gauss(A,b,L,U,P)ifif nargin5 nargin5 P=eye(length(A);P=eye(length(A);endendn=length(A);n=length(A);b=P*b;b=P*b;forfor j=1:n-1 j=1:n-1 b(j)=b(j)/L(j,j);b(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);ende

16、ndb(n)=b(n)/L(n,n);b(n)=b(n)/L(n,n);y=b;y=b;forfor j=n:-1:2 j=n:-1:2 y(j)=y(j)/U(j,j);y(j)=y(j)/U(j,j);y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);endendy(1)=y(1)/U(1,1);y(1)=y(1)/U(1,1);x=y;x=y;endendex1_2ex1_2%第二题第二题%第一问第一问A=10*eye(100)+diag(ones(1,99),-1)+diag(ones(1,99),1

17、);A=10*eye(100)+diag(ones(1,99),-1)+diag(ones(1,99),1);b=round(100*rand(100,1);b=round(100*rand(100,1);%平方根法平方根法L=Cholesky(A);L=Cholesky(A);x1_2_1_1=Gauss(A,b,L,L);x1_2_1_1=Gauss(A,b,L,L);%改进的平方根法改进的平方根法L,D=LDLt(A);L,D=LDLt(A);x1_2_1_2=Gauss(A,b,L,D*L);x1_2_1_2=Gauss(A,b,L,D*L);%第二问第二问A=hilb(40);A=h

18、ilb(40);b=sum(A);b=sum(A);b=b;b=b;6%平方根法平方根法L=Cholesky(A);L=Cholesky(A);x1_2_2_1=Gauss(A,b,L,L);x1_2_2_1=Gauss(A,b,L,L);%改进的平方根法改进的平方根法L,D=LDLt(A);L,D=LDLt(A);x1_2_2_2=Gauss(A,b,L,D*L);x1_2_2_2=Gauss(A,b,L,D*L);结果分别为x1_2_1_1=x1_2_1_1=7.2586 7.2586 8.4143 8.4143 -0.4013 -0.4013 8.5984 8.5984 5.4177 5

19、.4177 0.2249 0.2249 2.3336 2.3336 4.4389 4.4389 8.2772 8.2772 8.7890 8.7890 -0.1667 -0.1667 8.8784 8.8784 8.3824 8.3824 3.2978 3.2978 7.6401 7.6401 0.3014 0.3014 3.3457 3.3457 8.2418 8.2418 6.2368 6.2368 8.3906 8.3906 5.8575 5.8575 -0.9656 -0.9656 7.7981 7.7981 7.9842 7.9842 5.3601 5.3601 6.4143 6.4

20、143 6.4966 6.4966 2.6201 2.6201 6.3024 6.3024 0.3563 0.3563 7.1342 7.1342 -0.6985 -0.6985 2.8506 2.85067 0.1927 0.1927 0.2227 0.2227 7.5801 7.5801 5.9762 5.9762 1.6583 1.6583 9.4409 9.4409 -1.0677 -1.0677 4.2362 4.2362 2.7060 2.7060 6.7037 6.7037 7.2570 7.2570 0.7265 0.7265 4.4778 4.4778 3.4959 3.49

21、59 5.5637 5.5637 5.8675 5.8675 6.7614 6.7614 1.5180 1.5180 6.0582 6.0582 5.9002 5.9002 0.9397 0.9397 0.7031 0.7031 4.0294 4.0294 9.0032 9.0032 1.9382 1.9382 5.6150 5.6150 0.9120 0.9120 7.2652 7.2652 1.4360 1.4360 4.3749 4.3749 5.8146 5.8146 7.4791 7.4791 8.3942 8.3942 4.5789 4.5789 0.8169 0.8169 1.2

22、523 1.2523 1.6603 1.6603 8.1448 8.1448 0.8915 0.8915 7.9401 7.9401 0.7075 0.7075 8.9849 8.9849 2.4437 2.4437 1.5777 1.57778 1.7790 1.7790 5.6319 5.6319 3.9018 3.9018 2.3506 2.3506 7.5925 7.5925 4.7245 4.7245 4.1627 4.1627 8.6483 8.6483 1.3543 1.3543 6.8087 6.8087 6.5589 6.5589 2.6027 2.6027 5.4140 5

23、.4140 0.2577 0.2577 0.0090 0.0090 4.6522 4.6522 6.4685 6.4685 8.6626 8.6626 -0.0948 -0.0948 5.2856 5.2856 4.2385 4.2385 -0.6706 -0.67063.46713.4671x1_2_1_2=x1_2_1_2=7.2586 7.2586 8.4143 8.4143 -0.4013 -0.4013 8.5984 8.5984 5.4177 5.4177 0.2249 0.2249 2.3336 2.3336 4.4389 4.4389 8.2772 8.2772 8.7890

24、8.7890 -0.1667 -0.1667 8.8784 8.8784 8.3824 8.3824 3.2978 3.2978 7.6401 7.6401 0.3014 0.3014 3.3457 3.3457 8.2418 8.2418 6.2368 6.23689 8.3906 8.3906 5.8575 5.8575 -0.9656 -0.9656 7.7981 7.7981 7.9842 7.9842 5.3601 5.3601 6.4143 6.4143 6.4966 6.4966 2.6201 2.6201 6.3024 6.3024 0.3563 0.3563 7.1342 7

25、.1342 -0.6985 -0.6985 2.8506 2.8506 0.1927 0.1927 0.2227 0.2227 7.5801 7.5801 5.9762 5.9762 1.6583 1.6583 9.4409 9.4409 -1.0677 -1.0677 4.2362 4.2362 2.7060 2.7060 6.7037 6.7037 7.2570 7.2570 0.7265 0.7265 4.4778 4.4778 3.4959 3.4959 5.5637 5.5637 5.8675 5.8675 6.7614 6.7614 1.5180 1.5180 6.0582 6.0

26、582 5.9002 5.9002 0.9397 0.9397 0.7031 0.7031 4.0294 4.0294 9.0032 9.0032 1.9382 1.9382 5.6150 5.6150 0.9120 0.9120 7.2652 7.2652 1.4360 1.4360 4.3749 4.374910 5.8146 5.8146 7.4791 7.4791 8.3942 8.3942 4.5789 4.5789 0.8169 0.8169 1.2523 1.2523 1.6603 1.6603 8.1448 8.1448 0.8915 0.8915 7.9401 7.9401

27、0.7075 0.7075 8.9849 8.9849 2.4437 2.4437 1.5777 1.5777 1.7790 1.7790 5.6319 5.6319 3.9018 3.9018 2.3506 2.3506 7.5925 7.5925 4.7245 4.7245 4.1627 4.1627 8.6483 8.6483 1.3543 1.3543 6.8087 6.8087 6.5589 6.5589 2.6027 2.6027 5.4140 5.4140 0.2577 0.2577 0.0090 0.0090 4.6522 4.6522 6.4685 6.4685 8.6626

28、 8.6626 -0.0948 -0.0948 5.2856 5.2856 4.2385 4.2385 -0.6706 -0.67063.46713.4671x1_2_2_1=x1_2_2_1=1.0e+07*1.0e+07*0.0000 0.0000 -0.0000 -0.0000 0.0001 0.000111 -0.0004 -0.0004 -0.0014 -0.0014 0.0424 0.0424 -0.2980 -0.2980 1.1419 1.1419 -2.7335 -2.7335 4.2539 4.2539 -4.3018 -4.3018 2.7733 2.7733 -1.19

29、89 -1.1989 0.5406 0.5406 -0.3688 -0.3688 0.3285 0.3285 -0.4438 -0.4438 0.4621 0.4621 -0.2513 -0.2513 0.0565 0.0565 0.0000 0.0000 -0.0051 -0.0051 0.0071 0.0071 -0.0027 -0.0027 -0.0031 -0.0031 0.0036 0.0036 -0.0019 -0.0019 0.0009 0.0009 0.0002 0.0002 -0.0002 -0.0002 -0.0006 -0.0006 0.0004 0.0004 0.000

30、1 0.0001 -0.0002 -0.0002 0.0001 0.0001 0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 x1_2_2_2=x1_2_2_2=1.0000 1.0000 1.0000 1.0000 0.9998 0.9998 1.0011 1.0011 1.0064 1.006412 0.8681 0.8681 1.8034 1.8034 -1.5693 -1.5693 5.5763 5.5763 -2.5315 -2.5315 -1.7693 -1.7693 10.48

31、83 10.4883 -6.2807 -6.2807 0.5882 0.5882 -4.7157 -4.7157 22.8299 22.8299 -19.9134 -19.9134 8.7032 8.7032 10.3265 10.3265 -25.2140 -25.2140 10.0282 10.0282 12.3882 12.3882 -1.9425 -1.9425 14.1891 14.1891 -12.0552 -12.0552 -0.5803 -0.5803 -12.4791 -12.4791 8.5652 8.5652 9.8724 9.8724 -10.5502 -10.5502

32、 16.3871 16.3871 -5.8132 -5.8132 13.4216 13.4216 11.1767 11.1767 -64.3154 -64.3154 46.3837 46.3837 12.6957 12.6957 -21.7556 -21.7556 12.1204 12.1204 -1.9342 -1.93423.用第 1 题的程序求解第 2题的两个方程组并比较所有的计算结果，然后评价各个方法的优劣。Sol：Gauss表示不选主元的 Gauss消去法，PGauss 表示列主元 Gauss消去法。13计算脚本为：%第三题第三题%第一问第一问A=10*eye(100)+diag(o

33、nes(1,99),-1)+diag(ones(1,99),1);A=10*eye(100)+diag(ones(1,99),-1)+diag(ones(1,99),1);b=round(100*rand(100,1);b=round(100*rand(100,1);%不选主元不选主元 GaussGauss 消去法消去法L,U=GaussLA(A);L,U=GaussLA(A);x1_3_1_1=Gauss(A,b,L,U);x1_3_1_1=Gauss(A,b,L,U);%列主元列主元 GaussGauss 消去法消去法L,U,P=GaussCol(A);L,U,P=GaussCol(A);

34、x1_3_1_2=Gauss(A,b,L,U,P);x1_3_1_2=Gauss(A,b,L,U,P);%第二问第二问A=hilb(40);A=hilb(40);b=sum(A);b=sum(A);b=b;b=b;%不选主元不选主元 GaussGauss 消去法消去法L,U=GaussLA(A);L,U=GaussLA(A);x1_3_2_1=Gauss(A,b,L,U);x1_3_2_1=Gauss(A,b,L,U);%列主元列主元 GaussGauss 消去法消去法L,U,P=GaussCol(A);L,U,P=GaussCol(A);x1_3_2_2=Gauss(A,b,L,U,P);x

35、1_3_2_2=Gauss(A,b,L,U,P);ex1_2;ex1_2;y1=1:100;y2=1:40;y1=1:100;y2=1:40;subplot(4,2,1);plot(y1,x1_2_1_1);title(subplot(4,2,1);plot(y1,x1_2_1_1);title(平方根法平方根法 11););subplot(4,2,2);plot(y1,x1_2_1_2);title(subplot(4,2,2);plot(y1,x1_2_1_2);title(改进的平方根法改进的平方根法 11););subplot(4,2,3);plot(y1,x1_3_1_1);titl

36、e(subplot(4,2,3);plot(y1,x1_3_1_1);title(Gauss1Gauss1););subplot(4,2,4);plot(y1,x1_3_1_2);title(subplot(4,2,4);plot(y1,x1_3_1_2);title(PGauss1PGauss1););subplot(4,2,5);plot(y2,x1_2_2_1);title(subplot(4,2,5);plot(y2,x1_2_2_1);title(平方根法平方根法 22););subplot(4,2,6);plot(y2,x1_2_2_2);title(subplot(4,2,6);

37、plot(y2,x1_2_2_2);title(改进的平方根法改进的平方根法 22););subplot(4,2,7);plot(y2,x1_3_2_1);title(subplot(4,2,7);plot(y2,x1_3_2_1);title(Gauss2Gauss2););14subplot(4,2,8);plot(y2,x1_3_2_2);title(subplot(4,2,8);plot(y2,x1_3_2_2);title(PGauss2PGauss2););平 方 根 法 11050-501020304050Gauss11050-50 x 107改 进 的 平 方 根 法 1105

38、060708090100-501020304050PGauss11050607080901001020304050平 方 根 法 260708090100-50102030405060708090100改 进 的 平 方 根 法 250050-50-505101520Gauss22000500025303540-10005101520PGauss225303540-2000510152025303540-5000510152025303540平方根法和改进的平法根法计算量更小，计算过程稳定，但使用范围窄；不选主元和列主元的 Gauss消去法计算量较大，但适用范围广。例题 1.3.2考虑对称正定

39、线性方程组 Ax=b，其中向量 b是随机生成的，其元素是服从区间0,1上均匀分布的随机数，矩阵ALLT,这里 L是随机生成的一个下三角矩阵，其元素是服从区间1,2上均匀分布的随机数。对 n=10，20，.,500 分别应用 Gauss消去法、列主元 Gauss消去法和 Cholesky分解法求解该方程组，画出它们所用的 CPU 时间，其中“Gauss”表示 Gauss消去法、“PGauss”表示列主元 Gauss消去法,“Cholesky”表示 Cholesky分解法。Sol：经试验知，对应课本上图 1.1所示的 Cholesky 分解法应为改进后的 Cholesky分解法即LDLT分解。15

40、此处所用的 CPU 时间利用 cputime测量。计算脚本为 eg1_3_2clc;clear;clc;clear;forfor i=1:50;i=1:50;n=i*10;n=i*10;b=rand(n,1);b=rand(n,1);L=tril(unifrnd(1,2,n,n);L=tril(unifrnd(1,2,n,n);A=L*L;A=L*L;t1(i)=cputime;t1(i)=cputime;L1,U1=GaussLA(A);L1,U1=GaussLA(A);x1=Gauss(A,b,L1,U1);x1=Gauss(A,b,L1,U1);t1(i)=cputime-t1(i);t

41、1(i)=cputime-t1(i);t2(i)=cputime;t2(i)=cputime;L2,U2,P2=GaussCol(A);L2,U2,P2=GaussCol(A);x2=Gauss(A,b,L2,U2,P2);x2=Gauss(A,b,L2,U2,P2);t2(i)=cputime-t2(i);t2(i)=cputime-t2(i);t3(i)=cputime;t3(i)=cputime;L3=LDLt(A);L3=LDLt(A);x3=Gauss(A,b,L3,L3);x3=Gauss(A,b,L3,L3);t3(i)=cputime-t3(i);t3(i)=cputime-t3(i);endendN=10:10:500;N=10:10:500;plot(N,t1,plot(N,t1,o-o-,N,t2,N,t2,.-.-,N,t3,N,t3,*-*-););legend(legend(GaussGauss,PGaussPGauss,CholeskyCholesky););结果为161.5GaussPGaussCholesky10.5005010015020025030035040045050017