数值线性代数第二版徐树方高立张平文上机习题第三章实验报告.docx
《数值线性代数第二版徐树方高立张平文上机习题第三章实验报告.docx》由会员分享,可在线阅读,更多相关《数值线性代数第二版徐树方高立张平文上机习题第三章实验报告.docx(11页珍藏版)》请在冰豆网上搜索。
数值线性代数第二版徐树方高立张平文上机习题第三章实验报告
第三章上机习题
用你所熟悉的的计算机语言编制利用QR分解求解线性方程组和线性最小二乘问题的通用子程序,并用你编制的子程序完成下面的计算任务:
(1)求解第一章上机习题中的三个线性方程组,并将所得的计算结果与前面的结果相比较,说明各方法的优劣;
(2)求一个二次多项式
,使得在残向量的2范数下最小的意义下拟合表3.2中的数据;
表3.2
ti
-1
-0.75
-0.5
0
0.25
0.5
0.75
yi
1
0.8125
0.75
1
1.3125
1.75
2.3125
(3)在房产估价的线性模型
中,
分别表示税、浴室数目、占地面积、车库数目、房屋数目、居室数目、房龄、建筑类型、户型及壁炉数目,
代表房屋价格。
现根据表3.3和表3.4给出的28组数据,求出模型中参数的最小二乘结果。
(表3.3和表3.4见课本P99-100)
解分析:
(1)计算一个Householder变换H:
由于
,则计算一个Householder变换H等价于计算相应的
。
其中
。
在实际计算中,
为避免出现两个相近的数出现的情形,当
时,令
;
为便于储存,将
规格化为
,相应的,
变为
为防止溢出现象,用
代替
(2)QR分解:
利用Householder变换逐步将
转化为上三角矩阵
,则有
,其中
,
。
在实际计算中,从
,若
,依次计算
对应的
即对应的
,
,将
储存到
,
储存到
,迭代结束后再次计算
,有
,
(
时
)
(3)求解线性方程组
或最小二乘问题的步骤为
i计算
的QR分解;
ii计算
,其中
iii利用回代法求解上三角方程组
(4)对第一章第一个线性方程组,由于R的结果最后一行为零,故使用前代法时不计最后一行,而用运行结果计算
。
运算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);
ifsigma==0
belta=0;
else
alpha=sqrt(x
(1)^2+sigma);
ifx
(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计算
的QR分解[Q,R]=QRfenjie(A)
function[Q,R]=QRfenjie(A)
[m,n]=size(A);
Q=eye(m);
forj=1:
n
ifj[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,:
));
forj=1:
n
ifjH=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)
functionx=qiandaifa(L,b)
n=length(b);
forj=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('改进的平方根法');
%第二题第二问
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.500.250.50.75];
y=[10.81250.7511.31251.752.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);
Q1=Q(:
1:
12);
c=Q1'*y;
x=huidaifa(R,c);
x=x'
计算结果为
(1)第一章上机习题中的三个线性方程组结果对比图依次为
以第二个线性方程组为例,比较各方法的运行速度。
依次为QR分解,不选主元的Gauss消去法,列主元Gauss消去法,平方根法,改进的平方根法。
Elapsedtimeis0.034588seconds.
Elapsedtimeis0.006237seconds.
Elapsedtimeis0.009689seconds.
Elapsedtimeis0.030862seconds.
Elapsedtimeis0.007622seconds.
(2)二次多项式的系数为
x=
1.0000
1.0000
1.0000
(3)房产估价的线性模型的系数为
x=
Columns1through6
2.07750.71899.68020.153513.67961.9868
Columns7through12
-0.9582-0.4840-0.07361.01871.44352.9028
结果分析
对第一章上机习题中的第二个线性方程组利用五种求解方法求解所需时间可知,不选主元的Gauss消去法,列主元Gauss消去法,改进的平方根法较快,所需时间大致在一个数量级,QR分解,平方根法,所需时间较慢,所需时间在一个数量级上。