數(shù)值線性代數(shù)第二版徐樹(shù)方高立張平文上機(jī)習(xí)題第三章實(shí)驗(yàn)報(bào)告

第三章上機(jī)習(xí)題用你所熟悉的的計(jì)算機(jī)語(yǔ)言編制利用QR分解求解線性方程組和線性最小二乘問(wèn)題的通用子程序，并用你編制的子程序完成下面的計(jì)算任務(wù)：求解第一章上機(jī)習(xí)題中的三個(gè)線性方程組說(shuō)明各方法的優(yōu)劣；y=at2+bt+c2數(shù)據(jù)；表

yx0

ax1

ax2

ax1111aa1

,,a11

分別表示稅、浴室數(shù)目、占地面積、車(chē)庫(kù)數(shù)目、房屋數(shù)目、居室數(shù)目、房齡、建筑類(lèi)型、戶型及壁爐數(shù)目，y代表房屋價(jià)格?，F(xiàn)根據(jù)表和表給出的28組數(shù)據(jù)，求出模型中參數(shù)的最小二乘結(jié)果。（表和表見(jiàn)課本P99-100）解分析：計(jì)算一個(gè)HouseholderH：HI2wwT

IvvT

，則計(jì)算一個(gè)Householder變換H等價(jià)于計(jì)算相應(yīng)的、v。其中vx||x||2在實(shí)際計(jì)算中，

e2/(vTv。1

2n-(x2x2)2n為避免出現(xiàn)兩個(gè)相近的數(shù)出現(xiàn)的情形，當(dāng)x1

0時(shí)，令v1

x||x|| ；1 2為便于儲(chǔ)存，將v規(guī)格化為vv/v1

變?yōu)?v2/(vTv)1t-10iy11i為防止溢出現(xiàn)象，用xt-10iy11iQR利用Householder變換逐步將Amn

,mn轉(zhuǎn)化為上三角矩陣HHn

n1

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

R(1n。~j1njmxAjm,j對(duì)應(yīng)的(H

)j(mk1)(mk1)即對(duì)應(yīng)的vj

vj

(2mjAj1m,jj

儲(chǔ)存到d(j)，迭代結(jié)束后再次計(jì)算QH

Ij1

0~，QHHH

（nm時(shí)QHHH ）0Hj 0H j

1 2

1 2 n-1AxbiAQR計(jì)算

QTb，其中

Q(:,1:n)1 1 1利用回代法求解上三角方程組Rxc1R最后一行，而用運(yùn)行結(jié)果計(jì)算x 。84運(yùn)算matlab程序?yàn)镠ouseholder[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');%第二題第一問(wèn)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;%改進(jìn)的平方根法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('改進(jìn)的平方根法');%第二題第二問(wèn)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');%改進(jìn)的平方根法[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('改進(jìn)的平方根法');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\專(zhuān)業(yè)課\數(shù)值代數(shù)\','A2:L29');y=xlsread('E:\temporary\專(zhuān)業(yè)課\數(shù)值代數(shù)\','M2:M29');[Q,R]=QRfenjie(A);Q1=Q(:,1:12);c=Q1'*y;x=huidaifa(R,c);x=x'計(jì)算結(jié)果為第一章上機(jī)習(xí)題中的三個(gè)線性方程組結(jié)果對(duì)