﻿ 130213031+杨丽丽+作业九 - 计算数学达人 - 专，学者，数值代数，微分方程数值解

### 130213031+杨丽丽+作业九

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