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

130213031+杨丽丽+作业九

作者:130213031   发布时间:2015-06-09 23:03:40   浏览次数:174

1.设x=(1,0,4,6,3,4)'   ,求一个Househoder变换和一个正数alpha ,使得Hx=(1,alpha,4,6,0,0)'

过程:

%Househoder变换

function H=house(x)
n=length(x);
elta=norm(x,inf);
x=x/elta;
sigma=x(2:n)'*x(2:n);
v=zeros(n,1);
v(2:n)=x(2:n);
if sigma==0
    beta=0;
else
    alpha=sqrt(x(1)*x(1)+sigma);
    if x(1)<=0
        v(1)=x(1)-alpha;
    else
        v(1)=-sigma/(x(1)+alpha);
    end
   beta=2*v(1)*v(1)/(sigma+v(1)*v(1));
   
v=v/v(1);
end
I=eye(n);
H=I-beta*v*v';  

%test_house

clear all
x=[1;0;4;6;3;4];
n=size(x,1);
H=house(x);
H*x;
x1=[0;3;4];
H1=house(x1);
h=eye(n);
h([2,5,6],[2,5,6])=H1;
h*x
 

结果:

ans =

    1.0000
    5.0000
    4.0000
    6.0000
         0
   -0.0000

2.假定x和y 是R^n中的两个单位向量,给出一种使用givens 变换的方法,计算一个正交阵Q,使得Qx=y

过程:

function G=givens(a,b)
if b==0
    c=1;s=0;
else
    if abs(b)>abs(a)
        t=a/b;s=1/sqrt(1+t*t);c=s*t;
    else
        t=b/a;c=1/sqrt(1+t*t);s=c*t;
    end
end

G=[c,s;-s,c];
 

%test_givens

n=6;
x1=rand(n,1);
x1=x1/norm(x1)
y1=rand(n,1);
y1=y1/norm(y1)
G=eye(n,n);
g=eye(n,n);
for i=2:n  
I=eye(n,n);
G1=givens(x1(1),x1(i))
I([1,i],[1,i])=G1;
g=I*g;
x1=I*x1;
end

for j=2:n  
I=eye(n,n);
G2=givens(y1(1),y1(j))
I([1,j],[1,j])=G2;
G=I*G;
y1=I*y1;
end

Q=inv(g)*G
x1
y1

结果:

Q =

    0.8965    0.3754   -0.1493   -0.0578    0.1680   -0.0400
   -0.4100    0.9042   -0.0757   -0.0293    0.0852   -0.0203
    0.0876    0.1061    0.9787   -0.0483    0.1405   -0.0334
    0.0434    0.0525    0.0359    0.9967    0.0249   -0.0059
   -0.1338   -0.1620   -0.1107   -0.0060    0.9712   -0.0176
    0.0284    0.0344    0.0235    0.0013    0.0304    0.9983

x1 =

    0.5691
    0.2886
    0.4758
    0.0844
    0.2508
    0.5445


y1 =

    0.4191
    0.5076
    0.3469
    0.0189
    0.4492
    0.4941

 

3.例3.3.1

过程:

clear all
format short
b=[3;2;-1];
x=[1;sqrt(3)];
% r=zeros(3,1);
% y=zeros(3,1);
 for k=[10^5,10^7,10^9]
A=[k/sqrt(2)+1/2,-k/sqrt(6)+sqrt(3)/2;-k/sqrt(2),k/sqrt(6);k/sqrt(2)-1/2,-k/sqrt(6)-sqrt(3)/2];


% cholesky分解求解
B=A'*A;
c=A'*b;
c=cholesky(B,c);
r=norm(A*c-b)
y=norm(x-c)


% guasss_pp
B=A'*A;
c=A'*b;
c=gausslu_pp(B,c);
m=norm(A*c-b)
n=norm(c-x)

end

结果:

r =

    2.4495


y =

  9.5375e-007


m =

    2.4495


n =

  2.3849e-007


r =

    2.4495


y =

    0.0020


m =

    2.4495


n =

    0.0020


r =

    3.6973


y =

    1.9583


m =

   NaN


n =

  Inf

 

%QR分解

过程:

    function [Q,R]=QR(A)
 m=size(A,1);
n=size(A,2);
for j=1:n
    I=eye(m-j+1);
    if j<m
        [v,beta]=house1(A(j:m,j));
        A(j:m,j:n)=(I-beta*v*v')*A(j:m,j:n);
        d(j)=beta;
        A(j+1:m,j)=v(2:m-j+1);
    end
end
  R=triu(A,0);
  V1=tril(A,-1);
  q=eye(n);
   for k=2:n
       I=eye(m-k+1);
       I1=eye(k-1);
       vk=ones(n-k+1,1);
       vk(2:(n-k+1))=A(k:n,k);
       HK=I-d(k)*vk*vk';
       H=diag(I1,HK);
       q=q*H;
   end
   Q=zeros(size(R,1),size(R,2));
   Q([1:size(q,1)],[1:size(q,2)])=q;
   Q
   R

结果:暂无

 

 







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

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

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