现在时间是:
当前位置:首 页 >> 矩阵计算>> 教学区>> 文章列表

130213034+张燕+作业9(Householder变换,Givens变换,QR分解)

作者:130213034   发布时间:2015-06-09 21:57:03   浏览次数:443

上机作业9Householder变换,Givens变换,QR分解

 题目:

.P97   3

.P98    6

.P96    3.3.1

 

 

程序1Householder变换

house.m

function [v,beta]=house(x)

%Householder变换

%x的除了第一行外,全化为0

n=length(x);

elta=norm(x,inf);

x=x/elta;

sigma=x(2:n)'*x(2:n);

v=zeros(n,1);    %注意将v变成列向量!!!!

v(1)=1;

v(2:n)=x(2:n);

if sigma==0

    beta=0;

else

    alpha=sqrt(v(1)^2+sigma);

    if x(1)<=0

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

    else

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

    end

    beta=2*v(1)^2/(sigma+v(1)^2)

    v=v/v(1);

end

 

 

程序2

test_house.m

%test House

clear all

 format long

 %x的除了第一行外,全化为0

x=[1;0;4;6;3;4];

n=length(x);

I=eye(n);

[v,beta]=house(x)

H=I-beta*v*v';

H*x

 

 

运行结果:

v =

 

   1.000000000000000

                   0

  -0.510740824224823

  -0.766111236337235

  -0.383055618168617

  -0.510740824224823

 

 

beta =

 

   0.886772296585541

 

 

ans =

 

   8.831760866327848

                   0

   0.000000000000000

   0.000000000000000

                   0

   0.000000000000000

 

程序3Givens变换

givens.m

function [c,s]=givens(a,b)

%Givens变换

%用数a将数b化为0

if b==0

    c=1;

    s=0;

else

    if abs(b)>abs(a)

        tau=a/b;

        s=1/sqrt(1+tau^2);

        c=s*tau;

    else

        tau=b/a;

        c=1/sqrt(1+tau^2);

        s=c*tau;

    end

end

 

 

程序4

test_givens.m

%test givens 

%x的第二行的0,将x5行的3化为0

x=[1,0,4,6,3,4]';

[c1,s1]=givens(0,3);

I=eye(6);

G1=I;

G1(2,2)=c1;

G1(5,5)=c1;

G1(2,5)=s1;

G1(5,2)=-G1(2,5);

B=G1*x

 

运行结果:

B =

 

     1

     3

     4

     6

     0

     4    

 

程序4QR分解

HQR.m

%calculate QR分解:Householder方法

function [v,beta,R]=HQR(A)

n=size(A,2);

m=size(A,1);

I=eye(m);

for j=1:n

    if j<m

        [v,beta]=house(A(j:m,j));

        A(j:m,j:n)=(I(m-j+1)-beta*v*v')*A(j:m,j:n);

        d(j)=beta;

        A(j+1:m,j)=v(2:m-j+1);

    end

end

A;

R=triu(A,0);

for i=1:n

    v(1,i)=1%vi的第1个数 ,?列向量

    v(2:m-i+1,i)=A(i+1:m,i)  %vi的第j个数 ,?列向量

end

beta=d;

 

 

 

程序5:

%test HQR

function [x]=HQRjie(A,b)

[v,beta,R]=HQR(A)

[m,n]=size(A)

H=cell(n,1);

HH=cell(n,1)

for k=1:n

   H{k}=eye(n-k+2)-beta(k)*v(1:m-k+1,k)*(v(1:m-k+1,k))';

end

 

HH{1}=H{1};

H;

for i=2:n

    HH{i}=blkdiag(eye(i-1),H{i});

end

 

B=HH{1};

HH;

for j=2:n

    B=HH{j}*B;

end

%创建Q'  

QQ=B;

%检验

x=[1;sqrt(3)];

err=QQ*A*x-QQ*b

 

bb=QQ*b

  

   

%此时方程Ax=b,变为Q'Ax=Q'b,Rx=bb

 

%下用回代法求解上三角方程组

 

for j=n:-1:2

    bb(j)=bb(j)/R(j,j);

    bb(1:j-1)=bb(1:j-1)-bb(j)*R(1:j-1,j);

end

bb(1)=bb(1)/R(1,1);

x=bb;

 

 

程序6   test_HQRjie

k=1e5;

U=[1/sqrt(3),1/sqrt(2);-1/sqrt(3),0;1/sqrt(3),-1/sqrt(2)]

V=[sqrt(3)/2,1/2;-1/2,sqrt(3)/2]

S=[sqrt(2)*k,0;0,sqrt(2)]

A=U*S*V'

b=[3;2;-1];

 

[x]=HQRjie(A,b)

 

运行结果:

x =

 

  1.0e-004 *

 

    0.3518

    0.3728

         0

 (结果不对...)

 

 

 

第一题P97 3

1.1程序:

%P97 3,Householder

x=[1,0,4,6,3,4]';

n=length(x);

I=eye(n);

%H*x=[1,alpha,4,6,0,0]';

' alpha的值为'

alpha=sqrt(norm(x)^2-(1^2+4^2+6^2+0^2+0^2))

Hx=[1,alpha,4,6,0,0]'

v=x-Hx;

w=v/norm(v);

'所求Householder变换为'

H=I-2*w*w'

 

%checkout the result

HH=H*x;

HH-Hx

 

1.2运行结果:

ans =

 

alpha的值为

 

 

alpha =

 

    5.0000

 

 

Hx =

 

    1.0000

    5.0000

    4.0000

    6.0000

         0

         0

 

 

ans =

 

所求Householder变换为

 

 

H =

 

    1.0000         0         0         0         0         0

         0   -0.0000         0         0    0.6000    0.8000

         0         0    1.0000         0         0         0

         0         0         0    1.0000         0         0

         0    0.6000         0         0    0.6400   -0.4800

         0    0.8000         0         0   -0.4800    0.3600

 

 

ans =

 

  1.0e-014 *

 

         0

   -0.1776

         0

         0

    0.1110

    0.1332

 

 

 

 







上一篇:没有了    下一篇:没有了

Copyright ©2019    计算数学达人 All Right Reserved.

技术支持:自助建站 | 领地网站建设 |短信接口 |燕窝 版权所有 © 2005-2019 lingw.net.粤ICP备16125321号 -5