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

时间：2020-10-29 16:32:06　来源：勤学考试网本文已影响人

- 1 -

第三章上机习题

用你所熟悉的的计算机语言编制利用QR分解求解线性方程组和线性最小二乘问题的通用子程序，并用你编制的子程序完成下面的计算任务：

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

2+bt+cy=at3.2范数下最小的意义下拟合表2）求一个二次多项式，使得在残向量的2（中的数据；表 3.2

t i10000501005010

-1

-0.75

-0.5

0.25

0.5

0.75

y i

0.8125

0.75

1.3125

1.75

2.3125

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

y?x?ax?axax 110211121a,a,?,a分别表示税、浴室数目、占地面积、车库数目、房屋数目、居室数目、房中，1121y代表房屋价格。现根据表3.3和表3.4给出的龄、建筑类型、户型及壁炉数目，28组数据，求出模型中参数的最小二乘结果。

（表3.3和表3.4见课本P99-100）

解分析：

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

TTvvww?I?H?I?2v、。H由于等价于计算相应的变换，则计算一个HouseholderT?)/(ve||x||,v?2v?x? 。其中12 在实际计算中，22)?x-(xn20xv 时，令；为避免出现两个相近的数出现的情形，当 11||xx?||21T2)/(vvv?2vv?/vv ，相应的，规格化为变为为便于储存，将1||/||xx 为防止溢出现象，用代替? 分解：）（2QRAHnA?HHm,?，则有转化为上三角矩阵变换逐步将利用Householder1n?m?nn1- 1 -

RQ?HH?HR(1:n,:)QA?。

　，，其中n210~(H)mj?)):m,jx?A((jj?1:n 对应的，依次计算在实际计算中，从，若j)1?m?vjv(2:d()j)?1:m,jA(j，迭代结束，将，，储存到即对应的储存到jjjjI01j?HQ?HH?HQ?HH?Hmn?Q~）时（后再次计算，有，12-12n1njH0jAx?b或最小二乘问题的步骤为（3）求解线性方程组

A的QR分解；计算 iT)n(:,1:Q?Qbc?Q 计算，其中 ii111c?Rx 利用回代法求解上三角方程组 iii1（4）对第一章第一个线性方程组，由于R的结果最后一行为零，故使用前代法时不计x。

　最后一行，而用运行结果计算84

运算matlab程序为

[v,belta]=house(x) 变换1 计算Householderfunction [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

- 2 -

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

[m,n]=size(A);

Q=eye(m);

for j=1:n

if j<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);

end

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

for j=1:n

if j<m

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

x=qiandaifa(L,b) 3 解下三角形方程组的前代法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);

- 3 -

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

%第二题第二问

- 4 -

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

%列主元Gauss消去法

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

x3_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('改进的平方根法');

5 求解二次多项式 ex3_2

clear;clc;

t=[-1 -0.75 -0.5 0 0.25 0.5 0.75];

y=[1 0.8125 0.75 1 1.3125 1.75 2.3125];

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)

6 求解房产估价的线性模型 ex3_3

clear;clc;

A=xlsread('E:\temporary\专业课\数值代数\cha3_3_4.xls','A2:L29');

y=xlsread('E:\temporary\专业课\数值代数\cha3_3_4.xls','M2:M29');

[Q,R]=QRfenjie(A);

- 5 -

Q1=Q(:,1:12);

c=Q1'*y;

x=huidaifa(R,c);

x=x'

计算结果为

（1）第一章上机习题中的三个线性方程组结果对比图依次为

PGaussGauss解QR分x 10161.8

1.6141.411.221100.81-20.610.4-410.21

QR分解GaussPGauss平方根法改进的平方根

8888866666444442222200000-2-2-2-2-2050100505050500100010001000100

- 6 -

7根法改进的平GaussPGaussQR分解平方方根法x 105602500200400

42000150403003150010020200210005010100500000-200-1-50-500-2-40-100-100-1000-3-60-200-150-1500-4-80-5-2000-200-3000002040204020200402004040 分解，不选主元QR以第二个线性方程组为例，比较各方法的运行速度。依次为消去法，平方根法，改进的平方根法。Gauss消去法，列主元Gauss的Elapsed time is 0.034588 seconds.

Elapsed time is 0.006237 seconds.

Elapsed time is 0.009689 seconds.

Elapsed time is 0.030862 seconds.

Elapsed time is 0.007622 seconds.

）二次多项式的系数为（2x =

1.0000

3）房产估价的线性模型的系数为（x =

Columns 1 through 6

2.0775 0.7189 9.6802 0.1535 13.6796 1.9868

Columns 7 through 12

-0.9582 -0.4840 -0.0736 1.0187 1.4435 2.9028

结果分析