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

1、数值线性代数第二版徐树方高立张平文上机习的题目第一章实验报告材料上机习题1.先用你所熟悉的的计算机语言将不选主元和列主元Gauss消去法编写成通用的子程序；然后用你编写的程序求解84阶方程组；最后将你的计算结果与方程的准确解进展比拟，并就此谈谈你对Gauss消去法的看法。Sol：1先用matlab将不选主元和列主元Gauss消去法编写成通用的子程序，得到：不选主元Gauss消去法：得到满足列主元Gauss消去法：得到满足（2）用前代法解，得用回代法解，得求解程序为可缺省，缺省时默认为单位矩阵3计算脚本为ex1_1代码%算法1.1.3计算三角分解：Gauss消去法function L,U=Gau

2、ssLA(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

3、)=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);elsebreak;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,

4、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,8

5、3),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

6、消去法，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)

7、，此时缺省，缺省时默认为单位矩阵（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:nfor i=1:n v(i,1)=A(j,i)*A(i,i

8、);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

9、)*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

10、,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 =x1_2_1_2 =x1_2_2_1 = 1.0e+07 *x1_2_2_2 =3.用第1题的程序求解第2题的两个方程组并比拟所有的计算结果，然后评价各个方法的优劣。Sol：Gauss表示不选主元的Gauss消去法，PGauss表示列主元Gauss消去法。计算脚本为：%第三题%第一问A=10*eye(100)+diag

11、(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,

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

13、(改良的平方根法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);平方根法和改良的平法根法计算量更小，计算过程稳定，但使用X围窄；不选主元和列主元的Gauss消去法计算量较大，但适用X围广。考虑对称正定线性方程组Ax=b，其中向量b是随机生成的，其元素是服从区间0,1上均匀分布的随机数，矩阵,这里L是随机生成的一个下三角矩阵，其元素是服从区间1,2上均匀分布的随机数。对n=10，20，.,500分别应用Gauss消去法、列主元Gauss消去法和Cho

14、lesky分解法求解该方程组，画出它们所用的CPU时间，其中“Gauss表示Gauss消去法、“PGauss表示列主元Gauss消去法,“Cholesky表示Cholesky分解法。Sol：分解。此处所用的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);结果为