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

(完整word版)数值线性代数第二版徐树方高立张平文上机习题第一章实验报告上机习题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);for k=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);endU=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);for k=1:n-1s,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;if A(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;endendL=tril(A);U=triu(A);L=L-diag(diag(L)+diag(ones(1,n);P=eye(n);for i=1:n-1 temp=P(i,:); P(i,:)=P(u(i),:); P(u(i),:)=temp;endend%高斯消去法解线性方程组function x=Gauss(A,b,L,U,P)if nargin5 P=eye(length(A);endn=length(A);b=P*b;for j=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);endb(n)=b(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);endy(1)=y(1)/U(1,1);x=y;endex1_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);%解的比较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分解：平方根法）function L=Cholesky(A)n=length(A);for k=1:n A(k,k)=sqrt(A(k,k); A(k+1:n,k)=A(k+1:n,k)/A(k,k); for j=k+1:n A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k); endendL=tril(A);end%计算LDL分解：改进的平方根法function L,D=LDLt(A)n=length(A);for j=1:n for i=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); endL=tril(A);D=diag(diag(A);L=L-diag(diag(L)+diag(ones(1,n);end%高斯消去法解线性方程组function x=Gauss(A,b,L,U,P)if nargin5 P=eye(length(A);endn=length(A);b=P*b;for j=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);endb(n)=b(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);endy(1)=y(1)/U(1,1);x=y;endex1_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.67063.4671x1_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.67063.4671x1_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.0000x1_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.93423. 用第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_2clc;clear; for i=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);endN=10:10:500;plot(N,t1,o-,N,t2,.-,N,t3,*-);legend(Gauss,PGauss,Cholesky);结果为