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

第三章上机习题用你所熟悉的的计算机语言编制利用QR分解求解线性方程组和线性最小二乘问题的通用子程序，并用你编制的子程序完成下面的计算任务：求解第一章上机习题中的三个线性方程组说明各方法的优劣；y=at2+bt+c2数据；表

yx0

ax1

ax2

ax1111aa1

,,a11

分别表示税、浴室数目、占地面积、车库数目、房屋数目、居室数目、房龄、建筑类型、户型及壁炉数目，y代表房屋价格。现根据表和表给出的28组数据，求出模型中参数的最小二乘结果。（表和表见课本P99-100）解分析：计算一个HouseholderH：HI2wwT

IvvT

，则计算一个Householder变换H等价于计算相应的、v。其中vx||x||2在实际计算中，

e2/(vTv。1

2n-(x2x2)2n为避免出现两个相近的数出现的情形，当x1

0时，令v1

x||x|| ；1 2为便于储存，将v规格化为vv/v1

变为2v2/(vTv)1t-10iy11i为防止溢出现象，用xt-10iy11iQR利用Householder变换逐步将Amn

,mn转化为上三角矩阵HHn

n1

HA，则有1R0AQ ，其中QHHH0 1 2

R(1n。~j1njmxAjm,j对应的(H

)j(mk1)(mk1)即对应的vj

vj

(2mjAj1m,jj

储存到d(j)，迭代结束后再次计算QH

Ij1

0~，QHHH

（nm时QHHH ）0Hj 0H j

1 2

1 2 n-1AxbiAQR计算

QTb，其中

Q(:,1:n)1 1 1利用回代法求解上三角方程组Rxc1R最后一行，而用运行结果计算x 。84运算matlab程序为Householder[v,belta]=house(x)function[v,belta]=house(x)n=length(x);x=x/norm(x,inf);sigma=x(2:n)'*x(2:n);v=zeros(n,1);v(2:n,1)=x(2:n);ifsigma==0belta=0;elsealpha=sqrt(x(1)^2+sigma);ifx(1)<=0v(1)=x(1)-alpha;elseend

v(1)=-sigma/(x(1)+alpha);end

end

belta=2*v(1)^2/(sigma+v(1)^2);v=v/v(1,1);AQR[Q,R]=QRfenjie(A)function[Q,R]=QRfenjie(A)[m,n]=size(A);Q=eye(m);forj=1:nifj<m[v,belta]=house(A(j:m,j));H=eye(m-j+1)-belta*v*v';A(j:m,j:n)=H*A(j:m,j:n);d(j)=belta;A(j+1:m,j)=v(2:m-j+1);endendR=triu(A(1:n,:));forj=1:nifj<mH=eye(m);temp=[1;A(j+1:m,j)];H(j:m,j:m)=H(j:m,j:m)-d(j)*temp*temp';Q=Q*H;endendend解下三角形方程组的前代法x=qiandaifa(L,b)functionx=qiandaifa(L,b)n=length(b);forj=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);x=b;endex3_1clear;clc;%第一题A=6*eye(84)+diag(8*ones(1,83),-1)+diag(ones(1,83),1);b=[7;15*ones(82,1);14];n=length(A);%QR分解[Q,R]=QRfenjie(A);c=Q'*b;x1=huidaifa(R(1:n-1,1:n-1),c(1:n-1));x1(n)=c(n)-R(n,1:n-1)*x1;Gauss[L,U]=GaussLA(A);x1_1=Gauss(A,b,L,U);Gaussx1_2=Gauss(A,b,L,U,P);%解的比较figure(1);subplot(1,3,1);plot(1:n,x1);title('QR);subplot(1,3,2);plot(1:84,x1_1);title('Gauss');subplot(1,3,3);plot(1:84,x1_2);title('PGauss');%第二题第一问A=10*eye(100)+diag(ones(1,99),-1)+diag(ones(1,99),1);b=round(100*rand(100,1));n=length(A);%QR分解tic;[Q,R]=QRfenjie(A);c=Q'*b;x2=huidaifa(R,c);toc;GaussGaussx2_2=Gauss(A,b,L,U,P);toc;%平方根法tic;L=Cholesky(A);x2_3=Gauss(A,b,L,L');toc;%改进的平方根法tic;[L,D]=LDLt(A);x2_4=Gauss(A,b,L,D*L');toc;%解的比较figure(2);subplot(1,5,1);plot(1:n,x2);title('QR);subplot(1,5,2);plot(1:n,x2_1);title('Gauss');subplot(1,5,3);plot(1:n,x2_2);title('PGauss');subplot(1,5,4);plot(1:n,x2_3);title('平方根法');subplot(1,5,5);plot(1:n,x2_4);title('改进的平方根法');%第二题第二问A=hilb(40);b=sum(A);b=b';n=length(A);[Q,R]=QRfenjie(A);c=Q'*b;x3=huidaifa(R,c);Gauss[L,U]=GaussLA(A);x3_1=Gauss(A,b,L,U);Gaussx3_2=Gauss(A,b,L,U,P);%平方根法L=Cholesky(A);x3_3=Gauss(A,b,L,L');%改进的平方根法[L,D]=LDLt(A);x3_4=Gauss(A,b,L,D*L');%解的比较figure(3);subplot(1,5,1);plot(1:n,x3);title('QR);subplot(1,5,2);plot(1:n,x3_1);title('Gauss');subplot(1,5,3);plot(1:n,x3_2);title('PGauss');subplot(1,5,4);plot(1:n,x3_3);title('平方根法');subplot(1,5,5);plot(1:n,x3_4);title('改进的平方根法');ex3_2clear;clc;t=[-1 0 ];y=[1 1 ];A=ones(7,3);A(:,1)=t'.^2;A(:,2)=t';[Q,R]=QRfenjie(A);Q1=Q(:,1:3);c=Q1'*y';x=huidaifa(R,c)ex3_3clear;clc;A=xlsread('E:\temporary\专业课\数值代数\','A2:L29');y=xlsread('E:\temporary\专业课\数值代数\','M2:M29');[Q,R]=QRfenjie(A);Q1=Q(:,1:12);c=Q1'*y;x=huidaifa(R,c);x=x'计算结果为第一章上机习题中的三个线性方程组结果对