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

1、上机习题1先用你所熟悉的的计算机语言将不选主元和列主元Gauss消去法编写成通用的子程序；然后用你编写的程序求解84阶方程组；最后将你的计算结果与方程的精确解进行比较，并就此谈谈你对Gauss消去法的看法。Sol :(1) 先用matlab 将不选主元和列主元 Gauss消去法编写成通用的子程序，得到L,U,P：不选主元Gauss消去法：L,U GaussLA(A)得到L,U满足A LU列主元 Gauss 消去法：L,U,P GaussCol(A)得到 L,U ,P 满足 PA LU(2) 用前代法解Ly b or Pb，得y用回代法解Ux y，得x求解程序为x Gauss A,b, L,U,

2、P ( P可缺省，缺省时默认为单位矩阵)(3) 计算脚本为ex1_1代码%算法(计算三角分解：Gauss消去法)fun ctio nL,U=GaussLA(A)n=len gth(A);for k=1:n-1A(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( on es(1, n);end%算法122(计算列主元三角分解：列主元 Gauss消去法)fun ctio nL,U,P=

3、GaussCol(A)n=len gth(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)=0A(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);elsebreak ;endendL=tril(A);U=triu(A);L=L-diag(diag(L)+diag(o nes(1, n);P=eye(

4、 n);for i=1:n-1temp=P(i,:);P(i,:)=P(u(i),:);P(u(i),:)=temp;endend%高斯消去法解线性方程组fun ctio nx=Gauss(A,b,L,U,P)if nargin<5P=eye(le ngth(A);endn=len gth(A);b=P*b;for j=1:n-1b(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:2y(j)=y(j)/u(j,j)； y(1:j-1)=y(1:j-1)-y(

5、j)*U(1:j-1,j);endy(1)=y(1)/U(1,1);x=y;end ex1_1clc;clear;%第一题A=6*eye(84)+diag(8*o nes(1,83),-1)+diag(o nes(1,83),1);b=7;15*o nes(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);%解的比较'o-' );title( 'Gauss');subplot(1,3,1);p

6、lot(1:84,x1_1.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消去法，精确解为1,11 84）：由图，显然列主元消去法与精确解更为接近，不选主元的Gauss消去法误差比列主元消去法大，且不如列主元消去法稳定。Gauss消去法重点在于 A的分解过程，无论 A如何

7、分解，后面两步的运算过程不变。2先用你所熟悉的的计算机语言将平方根法和改进的平方根法编写成通用的子程序；然后用你编写的程序求解对称正定方程组Ax=b。Sol :（1 ）先用matlab将平方根法和改进的平方根法编写成通用的子程序，得到L，（D）:平方根法：L=Cholesky(A)改进的平方根法：L,D=LDLt(A)(2 )求解得Ly b求解得 LTxy or DLTxy求解程序为x=Gauss(A,b,L,U,P) ( ULT or U DLT , P此时缺省，缺省时默认为单位矩阵)(3)计算脚本为ex1_2代码%算法(计算Cholesky 分解：平方根法)fun ctio nL=Chol

8、esky(A)n=len gth(A);for k=1:nA(k,k)=sqrt(A(k,k);A(k+1: n,k)=A(k+1: n,k)/A(k,k);for j=k+1:nA(j: n,j)=A(j: n,j)-A(j: n,k)*A(j,k);endendL=tril(A);end%计算LDL '分解：改进的平方根法fun ctio nL,D=LDLt(A)n=len gth(A);for j=1:nfor i=1:nv(i,1)=A(j,i)*A(i,i);endA(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1,1);A(j+1: n,j)=(A(j+1:

9、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( on es(1, n);end%高斯消去法解线性方程组fun ctio nx=Gauss(A,b,L,U,P)if nargin<5P=eye(le ngth(A);endn=len gth(A);b=P*b;for j=1:n-1b(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:-

10、1:2y(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;end ex1_2%第二题%第一问A=10*eye(100)+diag(o nes(1,99),-1)+diag(o nes(1,99),1);b=ro un d(100*ra nd(100,1);%平方根法L=Cholesky(A);x1_2_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

11、);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.25868.4143-0.40138.59845.41770.22492.33364.43898.27728.7890-0.16678.87848.38243.29787.64010.30143.34578.24186.23688.39065.8575-0.96567.79817.98425.36016.41436.49662.62016.302

12、40.35637.1342-0.69852.85060.19270.22277.58015.97621.65839.4409-1.06774.23622.70606.70377.25700.72654.47783.49595.56375.86756.76141.51806.05825.90020.93970.70314.02949.00321.93825.61500.91207.26521.43604.37495.81467.47918.39424.57890.81691.25231.66038.14480.89157.94010.70758.98492.44371.57771.77905.6

13、3193.90182.35067.59254.72454.16278.64831.35436.80876.55892.60275.41400.25770.00904.65226.46858.6626-0.09485.28564.2385-0.67063.4671x1_2_1_2 =7.25868.4143-0.40138.59845.41770.22492.33364.43898.27728.7890-0.16678.87848.38243.29787.64010.30143.34578.24186.23688.39065.8575-0.96567.79817.98425.36016.4143

14、6.49662.62016.30240.35637.1342-0.69852.85060.19270.22277.58015.97621.65839.4409-1.06774.23622.70606.70377.25700.72654.47783.49595.56375.86756.76141.51806.05825.90020.93970.70314.02949.00321.93825.61500.91207.26521.43604.37495.81467.47918.39424.57890.81691.25231.66038.14480.89157.94010.70758.98492.44

15、371.57771.77905.63193.90182.35067.59254.72454.16278.64831.35436.80876.55892.60275.41400.25770.00904.65226.46858.6626-0.09485.28564.2385-0.67063.4671x1_2_2_1 =1.0e+07 *0.0000-0.00000.0001-0.0004-0.00140.0424-0.29801.1419-2.73354.2539-4.30182.7733-1.19890.5406-0.36880.3285-0.44380.4621-0.25130.05650.0

16、000 -0.00510.0071-0.0027-0.00310.0036-0.00190.00090.0002-0.0002-0.00060.00040.0001-0.00020.00010.0000 -0.00000.0000-0.0000-0.0000x1_2_2_2 =1.00001.00000.99981.00111.00640.86811.8034-1.56935.5763-2.5315-1.769310.4883-6.28070.5882-4.715722.8299-19.91348.703210.3265-25.214010.028212.3882-1.942514.1891-

17、12.0552-0.5803-12.47918.56529.8724-10.550216.3871-5.813213.421611.1767-64.315446.383712.6957-21.755612.1204-1.9342然后评价各个方法的优3.用第1题的程序求解第2题的两个方程组并比较所有的计算结果, 劣。Sol :Gauss 表示不选主元的 Gauss消去法，PGauss 表示列主元 Gauss 消去法。计算脚本为：%第三题%第一问A=10*eye(100)+diag(o nes(1,99),-1)+diag(o nes(1,99),1);b=ro un d(100*ra nd(10

18、0,1);%不选主元 Gauss消去法L,U=GaussLA(A);x1_3_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(

19、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);t

20、itle('改进的平方根法2');subplot(4,2,7);plot(y2,x1_3_2_1);title('Gauss2');平方根法1改进的平方根法110Gaussl010 2030405060708090100010203010000-57X 10平方根法2051015202530354050改进的平方根法2510152025303540Gauss2-50-1000PGauss2'PGauss2');subplot(4,2,8);plot(y2,x1_3_2_2);title(平方根法和改进的平法根法计算量更小，计算过程稳定，但使用范围窄;不选主元和列主元的 Gauss消去法计算量较大，但适用范围广。例题考虑对称正定线性方程组Ax=b，其中向量b是随机生成的，其元素是服从区间0,1上均匀分布的随机数，矩阵 A LLT,这里L是随机生成的一个下三角矩阵，其元素是服从区间1,2上均匀分布的随机数。对n=10 ，20， .,500分别应用 Gauss消去法、列主元 Gauss消去法和 Cholesky 分解法求解该方程组，画出它们所用的 CPU时间，其中