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

- 1 -

第三章上机习题

用你所熟悉的的计算机语言编制利用QR 分解求解线性方程组和线性最小二乘问题的

通用子程序，并用你编制的子程序完成下面的计算任务： （1）求解第一章上机习题中的三个线性方程组，并将所得的计算结果与前面的结果相比较，说明各方法的优劣；

（2）求一个二次多项式+bt+c y=at 2

，使得在残向量的2范数下最小的意义下拟合表3.2中的数据；

（3）在房产估价的线性模型

111122110x a x a x a x y ++++=

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

（1）计算一个Householder 变换H ：

由于T

T

vv I ww I H β-=-=2，则计算一个Householder 变换H 等价于计算相应的v 、β。其中)/(2,||||12v v e x x v T

=-=β。 在实际计算中，

为避免出现两个相近的数出现的情形，当01>x 时，令2

12221||||)

(-x x x x v n +++= ；

为便于储存，将v 规格化为1/v v v =，相应的，β变为)/(221v v v T

=β 为防止溢出现象，用∞||||/x x 代替 （2）QR 分解：

利用Householder 变换逐步将n m A n m ≥⨯,转化为上三角矩阵A H H H n n 11 -=Λ，则有

⎥⎦

⎤

⎢⎣⎡=0R Q A ，其中n H H H Q 21=，:),:1(n R Λ=。

在实际计算中，从n j :1=，若m j <，依次计算)),:((j m j A x =对应的)1()1()~

(+-⨯+-k m k m j H

即对应的j v ，j β，将)1:2(+-j m v j 储存到),:1(j m j A +，j β储存到)(j d ，迭代结束

后再次计算Q ，有⎥⎥⎦

⎤

⎢⎢

⎣⎡=-~001

j j j H I H ，n H H H Q 21=（m n =时1-21n H H H Q =） （3）求解线性方程组b Ax =或最小二乘问题的步骤为

i 计算A 的QR 分解；

ii 计算b Q c T

11=，其中):1(:,1n Q Q = iii 利用回代法求解上三角方程组1c Rx =

（4）对第一章第一个线性方程组，由于R 的结果最后一行为零，故使用前代法时不计最后一行，而用运行结果计算84x 。

运算matlab 程序为

1 计算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); if sigma==0 belta=0; else

alpha=sqrt(x(1)^2+sigma); if x(1)<=0

v(1)=x(1)-alpha; else

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

belta=2*v(1)^2/(sigma+v(1)^2); v=v/v(1,1); end end

2计算A的QR分解[Q,R]=QRfenjie(A)

function [Q,R]=QRfenjie(A)

[m,n]=size(A);

Q=eye(m);

for j=1:n

if j

[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);

end

end

R=triu(A(1:n,:));

for j=1:n

if j

H=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;

end

end

end

3 解下三角形方程组的前代法x=qiandaifa(L,b)

function x=qiandaifa(L,b)

n=length(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);

x=b;

end

4 求解第一章上机习题中的三个线性方程组ex3_1

clear;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);

%列主元Gauss消去法

[L,U,P]=GaussCol(A);

x1_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;

%不选主元Gauss消去法

tic;[L,U]=GaussLA(A);

x2_1=Gauss(A,b,L,U);toc;

%列主元Gauss消去法

tic;[L,U,P]=GaussCol(A);

x2_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('改进的平方根法'); %第二题第二问