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

上机习题 1.先用你所熟悉的的计算机语言将不选主元和列主元Gauss消去法编写成通用的子程序；然后用你编写的程序求解84阶方程组；最后将你的计算结果与方程的精确解进行比较，并就此谈谈你对Gauss消去法的看法。 Sol： （1）先用matlab将不选主元和列主元Gauss消去法编写成通用的子程序，得到： 不选主元Gauss消去法：得到满足 列主元Gauss消去法：得到满足 用前代法解，得 用回代法解，得 求解程序为（可缺省，缺省时默认为单位矩阵） （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); end U=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-1 [s,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; end end L=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; end end %高斯消去法解线性方程组 function x=Gauss(A,b,L,U,P) if nargin5 P=eye(length(A)); end n=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); end b(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); end y(1)=y(1)/U(1,1); x=y; end ex1_1 clc;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如何分解，后面两步的运算过程不变。 先用你所熟悉的的计算机语言将平方根法和改进的平方根法编写成通用的子程序；然后用你编写的程序求解对称正定方程组Ax=b。 Sol： （1）先用matlab将平方根法和改进的平方根法编写成通用的子程序，得到L，(D)： 平方根法：L=Cholesky(A) 改进的平方根法：[L,D]=LDLt(A) 求解得 求解得 求解程序为x=Gauss(A,b,L,U,P)（，此时缺省，缺省时默认为单位矩阵） 计算脚本为ex1_2 代码 %算法1.3.1（计算Cholesky分解：平方根