matlab上机作业报告(计算初等反射阵-用Householder变换法对矩阵A作正交分解-连续函数最佳平方逼近等)
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
-1.0000
>> H*x
ans =
4.0000
-3.7417
0.0000
0
二、设A为n阶矩阵，编写用Householder变换法对矩阵A作正交分解的程序。
1．程序功能：
给定n阶矩阵A，通过本程序用Householder变换法对矩阵A作正交分解，得出A＝QR
2．基本原理：
任一实列满秩的m×n矩阵A，可以分解成两个矩阵的乘积，即A＝QR，其中Q是具有法正交列向量的m×n矩阵，R是非奇异的n阶上三角阵。
Error: M=0
H =
1 0 0 0
0 1 00
0 0 1 0
0 0 0 1
情形2：
K值溢出：
>> x=[1,2,3,4]';
>> H=holderk(x,5)
Error: k>n
情形3：
K值为1：
>> x=[2,3,4,5]';
>> H=holderk(x,1)
H =
-0.2722 -0.4082 -0.5443 -0.6804
p=s*u(1);%计算p
if n==1 u=0; end%若x是1×1维向量，则u=0
(2)
function H=holderk(x,k)
%HOLDERK给定向量x≠0,数k，计算初等反射阵Hk，使HkX=Y，其中Y的第k+1项到最后项全为零；
%程序功能：函数holderk给定向量x≠0，数k，计算初等反射阵Hk，使HkX=Y，%程序功能：函数holder2给定向量x≠0,计算Householder初等变换阵的p,u；
情形4：
(1)K值为3：
>> x=[4,3,2,1]';
>> H=holderk(x,3)
H =
1.0000 0 00
0 1.0000 0 0
0 0 -0.8944 -0.4472
0 0 -0.4472 0.8944
检验:
>> det(H)
ans =
-1
>> H*x
ans =
4.0000
3.0000
[p,u]=holder2(x(k:n));%得到计算Householde初等变换阵的系数ρ、向量U；
H(k:n,k:n)=eye(n-k+1)-p\u*u'; %计算H(k:n,k:n)=I-p\u*u'；
5．使用示例:
情形1:
X为零向量
>> x=[0,0,0,0]';
>> H=holderk(x,1)
%输入：n维向量x；
%输出：[p,u]。p是Householder初等变换阵的系数ρ，
% u是Householder初等变换阵的向量U。
n=length(x);%得到n维向量x的维数；
p=1;u=0;%初始化p,u；
M=max(abs(x));%得到向量x的无穷范数，即x中绝对值最大的一项的绝对值；
-0.4082 0.8690 -0.1747 -0.2184
-0.5443 -0.1747 0.7671 -0.2911
-0.6804 -0.2184 -0.2911 0.6361
检验:
>> det(H)
ans =
-1.0000
>> H*x
ans =
-7.3485
0.0000
0.0000
0.0000
-2.2361
0
(2)K值为2:
>> x=[4,3,2,1]';
>> H=holderk(x,2)
H =
1.0000 0 0 0
0 -0.8018 -0.5345 -0.2673
0 -0.5345 0.8414 -0.0793
0 -0.2673 -0.0793 0.9604
>> det(H)
ans =
(7)
(8)按要求输出,结束
3．变量说明：
x－输入的n维向量；
n－n维向量x的维数；
M－M是向量x的无穷范数，即x中绝对值最大的一项的绝对值；
p－Householder初等变换阵的系数ρ；
u－Householder初等变换阵的向量U
s－向量x的二范数；
x－输入的n维向量；
n－n维向量x的维数；
p－Householder初等变换阵的系数ρ；
k，jj，ii－循环变量；
t1－计算上三角阵R的系数tj；
t2－计算正交矩阵Q的系数ti；
4．程序代码：
function [Q,A]=qrhh(A)
%QRHH用Householder变换法对n阶矩阵A作正交分解A=QR；
%程序功能：函数qrhh用Householder变换法对矩阵A作正交分解A=QR；
if M==0%如果x=0，提示出错,程序终止;
disp('Error:M=0');
return;
else
x=x/M;%规范化
end;
s=norm(x);%求x的二范数
if x(1)<0%首项为负，s值要变号
s=-s;
end
u=x;%除首项外，其余各项x,u相同
u(1)=s+x(1);%计算u的首项
一、给定向量x≠0，计算初等反射阵Hk。
1．程序功能：
给定向量x≠0，计算初等反射阵Hk。
2．基本原理：
若 的分量不全为零，则由
确定的镜面反射阵H使得 ；当 时，由
有
算法：
(1)输入x，若x为零向量，则报错
(2)将x规范化，
如果M＝0，则报错同时转出停机
否则
(3)计算 ，如果 ，则
(4)
(5)计算
(6)
算法：Biblioteka Baidu
(1)输入n阶矩阵A
(2)对 ，求Househoulder初等反射阵的 。
(3)计算上三角阵R，仍然存储在A
(4)计算正交阵Q
(5)按要求输出,结束
3．变量说明：
A－输入的n阶矩阵，同时用于存储上三角阵R；
n－矩（方）阵A的阶数；
Q－Q是具有法正交列向量的n阶矩阵；
p,u－向量A（k:n,k），对应初等反射阵的ρ,u
%输入：n维向量x,数k；
%输出：H。H是Householder初等变换阵，H*x=y，使得y的第k+1项到最后项全为零；
%引用函数：holder2；
n=length(x); %得到n维向量x的维数；
if k>n%如果k值溢出，报错；
disp('Error: k>n');
end
H=eye(n); %初始化H，并使H(1:k,1:k)=I；
%输入：n阶矩阵A；
%输出：[Q,A]。Q是具有法正交列向量的n阶矩阵，
% A（即R）是非奇异的n阶上三角阵，仍用输入的矩阵A存储。
u－Householder初等变换阵的向量U
k－数k，H*x=y，使得y的第k+1项到最后项全为零；
4．程序代码：
(1)
function [p,u]=holder2(x)
%HOLDER2给定向量x≠0,计算Householder初等变换阵的p,u
%程序功能：函数holder2给定向量x≠0,计算Householder初等变换阵的p,u；
>> H*x
ans =
4.0000
-3.7417
0.0000
0
二、设A为n阶矩阵，编写用Householder变换法对矩阵A作正交分解的程序。
1．程序功能：
给定n阶矩阵A，通过本程序用Householder变换法对矩阵A作正交分解，得出A＝QR
2．基本原理：
任一实列满秩的m×n矩阵A，可以分解成两个矩阵的乘积，即A＝QR，其中Q是具有法正交列向量的m×n矩阵，R是非奇异的n阶上三角阵。
Error: M=0
H =
1 0 0 0
0 1 00
0 0 1 0
0 0 0 1
情形2：
K值溢出：
>> x=[1,2,3,4]';
>> H=holderk(x,5)
Error: k>n
情形3：
K值为1：
>> x=[2,3,4,5]';
>> H=holderk(x,1)
H =
-0.2722 -0.4082 -0.5443 -0.6804
p=s*u(1);%计算p
if n==1 u=0; end%若x是1×1维向量，则u=0
(2)
function H=holderk(x,k)
%HOLDERK给定向量x≠0,数k，计算初等反射阵Hk，使HkX=Y，其中Y的第k+1项到最后项全为零；
%程序功能：函数holderk给定向量x≠0，数k，计算初等反射阵Hk，使HkX=Y，%程序功能：函数holder2给定向量x≠0,计算Householder初等变换阵的p,u；
情形4：
(1)K值为3：
>> x=[4,3,2,1]';
>> H=holderk(x,3)
H =
1.0000 0 00
0 1.0000 0 0
0 0 -0.8944 -0.4472
0 0 -0.4472 0.8944
检验:
>> det(H)
ans =
-1
>> H*x
ans =
4.0000
3.0000
[p,u]=holder2(x(k:n));%得到计算Householde初等变换阵的系数ρ、向量U；
H(k:n,k:n)=eye(n-k+1)-p\u*u'; %计算H(k:n,k:n)=I-p\u*u'；
5．使用示例:
情形1:
X为零向量
>> x=[0,0,0,0]';
>> H=holderk(x,1)
%输入：n维向量x；
%输出：[p,u]。p是Householder初等变换阵的系数ρ，
% u是Householder初等变换阵的向量U。
n=length(x);%得到n维向量x的维数；
p=1;u=0;%初始化p,u；
M=max(abs(x));%得到向量x的无穷范数，即x中绝对值最大的一项的绝对值；
-0.4082 0.8690 -0.1747 -0.2184
-0.5443 -0.1747 0.7671 -0.2911
-0.6804 -0.2184 -0.2911 0.6361
检验:
>> det(H)
ans =
-1.0000
>> H*x
ans =
-7.3485
0.0000
0.0000
0.0000
-2.2361
0
(2)K值为2:
>> x=[4,3,2,1]';
>> H=holderk(x,2)
H =
1.0000 0 0 0
0 -0.8018 -0.5345 -0.2673
0 -0.5345 0.8414 -0.0793
0 -0.2673 -0.0793 0.9604
>> det(H)
ans =
(7)
(8)按要求输出,结束
3．变量说明：
x－输入的n维向量；
n－n维向量x的维数；
M－M是向量x的无穷范数，即x中绝对值最大的一项的绝对值；
p－Householder初等变换阵的系数ρ；
u－Householder初等变换阵的向量U
s－向量x的二范数；
x－输入的n维向量；
n－n维向量x的维数；
p－Householder初等变换阵的系数ρ；
k，jj，ii－循环变量；
t1－计算上三角阵R的系数tj；
t2－计算正交矩阵Q的系数ti；
4．程序代码：
function [Q,A]=qrhh(A)
%QRHH用Householder变换法对n阶矩阵A作正交分解A=QR；
%程序功能：函数qrhh用Householder变换法对矩阵A作正交分解A=QR；
if M==0%如果x=0，提示出错,程序终止;
disp('Error:M=0');
return;
else
x=x/M;%规范化
end;
s=norm(x);%求x的二范数
if x(1)<0%首项为负，s值要变号
s=-s;
end
u=x;%除首项外，其余各项x,u相同
u(1)=s+x(1);%计算u的首项
一、给定向量x≠0，计算初等反射阵Hk。
1．程序功能：
给定向量x≠0，计算初等反射阵Hk。
2．基本原理：
若 的分量不全为零，则由
确定的镜面反射阵H使得 ；当 时，由
有
算法：
(1)输入x，若x为零向量，则报错
(2)将x规范化，
如果M＝0，则报错同时转出停机
否则
(3)计算 ，如果 ，则
(4)
(5)计算
(6)
算法：Biblioteka Baidu
(1)输入n阶矩阵A
(2)对 ，求Househoulder初等反射阵的 。
(3)计算上三角阵R，仍然存储在A
(4)计算正交阵Q
(5)按要求输出,结束
3．变量说明：
A－输入的n阶矩阵，同时用于存储上三角阵R；
n－矩（方）阵A的阶数；
Q－Q是具有法正交列向量的n阶矩阵；
p,u－向量A（k:n,k），对应初等反射阵的ρ,u
%输入：n维向量x,数k；
%输出：H。H是Householder初等变换阵，H*x=y，使得y的第k+1项到最后项全为零；
%引用函数：holder2；
n=length(x); %得到n维向量x的维数；
if k>n%如果k值溢出，报错；
disp('Error: k>n');
end
H=eye(n); %初始化H，并使H(1:k,1:k)=I；
%输入：n阶矩阵A；
%输出：[Q,A]。Q是具有法正交列向量的n阶矩阵，
% A（即R）是非奇异的n阶上三角阵，仍用输入的矩阵A存储。
u－Householder初等变换阵的向量U
k－数k，H*x=y，使得y的第k+1项到最后项全为零；
4．程序代码：
(1)
function [p,u]=holder2(x)
%HOLDER2给定向量x≠0,计算Householder初等变换阵的p,u
%程序功能：函数holder2给定向量x≠0,计算Householder初等变换阵的p,u；