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

130213027王者 上机作业九

作者:王者   发布时间:2015-06-09 22:08:08   浏览次数:122

P97 T3

程序代码:

function [v,B]=house1(x)

n=length(x);

N=norm(x,inf);

x=x/N;

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

v=zeros(n,1);

v(1)=1;

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

if p==0

    B=0;

else

    a=sqrt(x(1)*x(1)+p);

    if x(1)<=0

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

    else

    v(1)=-p/(x(1)+a);

    end

    B=2*v(1)*v(1)/(p+v(1)*v(1));

    v=v/v(1);

end

 上机运行结果:

clear all

format short

x=[0 3 4]';

[v,B]=house1(x)

n=length(x);

I=eye(n);

w=B*x'*v;

Hx=x-v*w'

H=I-B*v*v'

A=eye(6);

A([2,5,6],[2,5,6])=H

%A(2,2)=H(1,1);

%A(2,5)=H(1,2);

%A(2,6)=H(1,3);

%A(5,2)=H(2,1);

%A(5,5)=H(2,2);

%A(5,6)=H(2,3);

%A(6,2)=H(3,1);

%A(6,5)=H(3,2);

%A(6,6)=H(3,3)

x1=[1 0 4 6 3 4]';

A*x1

运行结果:

v =

    1.0000

   -0.6000

   -0.8000

B =

     1

Hx =

     5

     0

     0

H =

         0    0.6000    0.8000

    0.6000    0.6400   -0.4800

    0.8000   -0.4800    0.3600

A =

    1.0000         0         0         0         0         0

         0         0         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.0000

    5.0000

    4.0000

    6.0000

         0

   -0.0000

P98 T6

程序代码:

%Givens变换对单位向量x规范化

function [G]=givens(x)

n=length(x);

for i=1:n

    for j=1:n

        if i==j

            G(i,j)=1;

        else  G(i,j)=0;

    end

    end

end

        for i=n-1:-1:1

            if x(i+1,1)==0

                continue;

            else if abs(x(i+1,1))>abs(x(i,1))

        t=x(i,1)/x(i+1,1);

        s=1/sqrt(1+t*t);

        c=s*t;

          

    else

        t=x(i+1,1)/x(i,1);

        c=1/sqrt(1+t*t);

        s=c*t;

    end

            end    

    x(i,1)=c*x(i,1)+s*x(i+1,1);

    x(i+1,1)=0;

    for j=1:n

        a=G(i,j);b=G(i+1,j);

        G(i,j)=c*a+s*b;

        G(i+1,j)=c*b-s*a;

    end  

        end

end

2function [cs]=givens2(x,y)

n=length(x);

for i=1;n

    %cs*x=(y,a)'

a=sqrt(x(i,1)*x(i,1)+x(i+1,1)*x(i+1,1)-y(i,1)*y(i,1));

if a==0

    cs(i,1)=1;

    cs(i+1,1)=0;

else

    cs(i,1)=(x(i,1)*y(i,1)+a*x(i+1,1))/(norm(x,2)*norm(x,2));

    cs(i+1,1)=(x(i+1,1)*y(i,1)-a*x(i,1))/(norm(x,2)*norm(x,2));

    x(i,1)=y(i,1);

    x(i,1)=a;

end

end

end

3运行结果:

%function [G]=givens3(x,y)

n=6;

G=zeros(n,n);

A=rand(n,1);

B=rand(n,1);

%x,y为单位向量

x=A/norm(A);

%E=1

E=norm(x,2);

y=B/norm(B);

[G]=givens(x);

n=length(x);

%验证x规范,即G*x=1,0。。。)’

e=G*x;

t=0;

for i=1:n-1

    [cs]=givens2(x,y);

    for j=1:n

     c=G(i,j);

     s=G(i+1,j);

     G(i,j)=cs(1,1)*c+cs(2,1)*s;

     G(i+1,j)=cs(1,1)*s-cs(2,1)*c;

     t=t+y(i,1)*y(i,1);

    end

end

     G

     d=G*x;

     err=abs(y-d)

     if t==1

         break;

     end

上机运行结果:

G =

 

   -0.2484    0.7033    0.1195    0.0759    0.2662    0.1416

   -0.6027   -0.2850    0.0546    0.0347    0.1216    0.0647

    0.1564   -0.2116   -0.0400    0.1512    0.5301    0.2820

   -0.0406    0.0549   -0.6125   -0.2324    0.1233    0.0656

    0.0105   -0.0142    0.1589   -0.5868   -0.0100    0.2886

   -0.0035    0.0048   -0.0532    0.1966   -0.4012    0.6639

 

 

err =

 

    0.4740

    0.4093

    0.5319

    0.5784

    0.0288

    0.4540

 







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

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

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