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

第三章上机习题

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

用子程序，并用你编制的子程序完成下面的计算任务：

（1）求解第一章上机习题中的三个线性方程组，并将所得的计算结果与前面的结果相比较，说

明各方法的优劣；

（2）求一个二次多项式y=at 2+bt+c ，使得在残向量的 2 范数下最小的意义下拟合表中的

数据；

表

t i -1 0

y i 1 1

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

y x0 a1x1 a2 x2 a11x11

中， a1 ,a2 ,, a11分别表示税、浴室数目、占地面积、车库数目、房屋数目、居室数目、房

龄、建筑类型、户型及壁炉数目，y 代表房屋价格。现根据表和表给出的28 组数据，求出模型中参数的最小二乘结果。

（表和表见课本P99-100 ）

解分析：

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

由于H I 2ww T Ivv T，则计算一个Householder 变换 H 等价于计算相应的、 v 。

其中

v x ||

x

||

2

e1

, 2 /( T )

v v 。

在实际计算中，

为避免出现两个相近的数出现的情形，当x1 0 时，令v1 - ( x22 x n2 ) ；

x1 || x ||2

为便于储存，将v 规格化为 v v / v1，相应的，变为2v2 /(v T v)

1

为防止溢出现象，用x / || x || 代替

（2） QR分

解：

利用 Householder 变换逐步将 A m n , m n 转化为上三角矩阵H n H n 1 H 1 A ，则有

R

A Q，其中Q H1H 2 H n， R (1: n,:) 。

~

在实际计算中，从j 1: n ，若j m ，依次计算x A(( j : m, j )) 对应的( H j)( m k 1) ( m k 1) 即对应的 v j，j，将 v j (2 : m j 1) 储存到 A( j 1: m, j) ，j储存到 d ( j) ，迭代结束

后再次计算 Q ，有 H j I

j 1 0

H n（ n m 时 Q H 1H 2

~ ， Q H1H 2 H n-1 ）0 H j

（3）求解线性方程组Ax b 或最小二乘问题的步骤为

i计算 A 的QR分解；

ii计算 c1Q1T b ，其中 Q1Q (:,1: n)

iii利用回代法求解上三角方程组 Rx c1

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

运算 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);