现在时间是:
当前位置:首 页 >> 数值分析(英)>> 教学区>> 文章列表

17074190作业8

作者:17074190   发布时间:2019-06-09 21:00:18   浏览次数:104

 4.3 Exercises

问题:2.Apply classical Gram-Schmidt orthogonalization to find the full QR factorization of the following matrices:  (b)[-4 -4;-2 7;4 -5]

Answer:

y1=A1=[ -4 ; 2 ; 4] , r11=||y1||2=(16+4+16)^1/2=6

q1=y1/||y1||2=(1/6)y1=[-2/3 ; -1/3 ; 2/3 ]

y2=A2-q1(q1A2)=[-4 ; 7 ; -5]+[-2/3 ; -1/3 ; 2/3]3=[-6 ; 6 ; -3]

q2=y2/||y2||2=(1/9)[-6 ; 6 ; -3]=[-2/3 ; 2/3 ; -1/3]

r12=q1A2=-3 , r22=||y2||2=9

Adding athird A3=[1 ; 0 ; 0] leads to

y3=A3-q1(q1A3)-q2(q2A3)

 =[1 ; 0 ; 0]+[-2/3 ; -1/3 ; 2/3](2/3)+[-2/3 ; 2/3 ; -1/3](2/3)=[1/9 ; 2/9 ; 2/9]

q3=y3/||y3||2=3[1/9 ; 2/9 ; 2/9]=[1/3 ; 2/3 ; 2/3]

A=[-4 -4 ; -2 7 ; 4 -5]

=[-2/3 -2/3 1/3 ; -1/3 2/3 2/3 ; 2/3 -1/3 2/3][6 -3 ; 0 9 ; 0 0]=QR.

 

 

问题:4.Apply modified Gram-Schmidt orthogonalization to find the full QR factorization of the matrices in Exercise 2.

Answer:

y1=A1=[ -4 ; 2 ; 4] , r11=||y1||2=(16+4+16)^1/2=6

q1=y1/||y1||2=(1/6)y1=[-2/3 ; -1/3 ; 2/3 ]

y2=A2-q1(q1A2)=[-4 ; 7 ; -5]+[-2/3 ; -1/3 ; 2/3]3=[-6 ; 6 ; -3]

q2=y2/||y2||2=(1/9)[-6 ; 6 ; -3]=[-2/3 ; 2/3 ; -1/3]

r12=q1A2=-3 , r22=||y2||2=9

Adding athird A3=[1 ; 0 ; 0] leads to

y3(1)=[1 ; 0 ; 0]- [-2/3 ; -1/3 ; 2/3 ]q1A3=[5/9 ; -2/9 ; 4/9]

y3=y3(1)- [-2/3 ; 2/3 ; -1/3](-2/3)=[1/9 ; 2/9 ; 2/9]

q3=y3/||y3||2=3[1/9 ; 2/9 ; 2/9]=[1/3 ; 2/3 ; 2/3]

A=[-4 -4 ; -2 7 ; 4 -5]

=[-2/3 -2/3 1/3 ; -1/3 2/3 2/3 ; 2/3 -1/3 2/3][6 -3 ; 0 9 ; 0 0]=QR.

 

 

问题:6. Apply Householder reflectors orthogonalization to find the full QR factorization of the matrices in Exercise 2.

(a)[2 3;-2 -6;1 0]

Answer:

x=[2;-2;1],w=[3;0;0],v=w-x=[1;2;-1]

H1=I-2vv’/v’v

  =[1 0 0; 0 1 0; 0 0 1]-1/3[1 2 -1; 2 4 -2; -1 -2 1]

  =[2/3 -2/3 1/3; -2/3 -1/3 2/3; 1/3 2/3 2/3]

H1A=[2/3 -2/3 1/3; -2/3 -1/3 2/3; 1/3 2/3 2/3][2 3; -2 -6; 1 0]

   =[3 6; 0 0; 0 -3]

X~=[0;-3],w~=[3;0],v~=w~-x~=[3;3]

H2~=I-2P=[1 0;0 1]-1/9[9 9;9 9]

        =[0 -1;-1 0]

H2H1A=[1 0 0; 0 0 -1; 0 -1 0][2/3 -2/3 1/3; -2/3 -1/3 2/3; 1/3 2/3 2/3][2 3; -2 -6; 1 0]

     =[3 6; 0 3; 0 0]

     =R

A=H1H2R

=[2/3 -2/3 1/3; -2/3 -1/3 2/3; 1/3 2/3 2/3][1 0 0; 0 0 -1; 0 -1 0][3 6; 0 3; 0 0]

=[2/3 -1/3 1/3; -2/3 -2/3 1/3; 1/3 -2/3 -2/3][3 6; 0 3; 0 0]

=QR

(b)[-4 -4;-2 7;4 -5]

Answer:

x=[-4;-2;4],w=[6;0;0],v=w-x=[10;2;-4]

H1=I-2P=[1 0 0; 0 1 0; 0 0 1]-1/60[100 20 -40; 20 4 -8; -40 -8 16]

       =[-2/3 -1/3 2/3; -1/3 14/15 2/15; 2/3 2/15 11/15]

H1A=[6 -3; 0 36/5; 0 -27/5]

x~=[36/5;-27/5],w~=[9;0],v~=w~-x~=[9/5;27/5]

H2~=I-2P=[1 0;0 1]-5/81[81/25 243/25;243/25 729/25]

        =[1 0;0 1]-[1/5 3/5;3/5 9/5]

        =[4/5 -3/5; -3/5 -4/5]

H2H1A=[1 0 0;0 4/5 -3/5;0 -3/5 -4/5][6 -3;0 36/5;0 -27/5]

      =[6 -3;0 9;0 0]

      =R

A= H1H2R

 =[-2/3 -4/3 -1/3;-1/3 2/3 -2/3;2/3 -1/3 -2/3][6 -3;0 9;0 0]

 =QR

 

 

问题:13.Let P be the matrix defined in(4.29).Show (a) P^2=P (b) P is symmetric (c)Pv=v.

Answer:

(a)P=vv’/v’v

  P^2=(vv’/v’v)(vv’/v’v)

     =vv’vv’/(v’v)^2=v(v’v)v’/(v’v)^2=vv’/v’v=P

(b)P’=(vv’/v’v)’=(1/v’v)(vv’)’ =(1/v’v)(vv’)=P

So, P is symmetric.

(c)Pv=(vv’/v’v)v=v(v’v)/v’v=v.

 

 

4.3 Computer Problems

1.Write a Matlab program that implements classical GramSchmidt to find the reduced QR factorization.Check your work by comparing factorization of the matrices in Exercise 1 with the Matlab qr(A,0) command or equivalent.The factorization is unique up to signs of the entries of Q and R.

Code:

function GS(A)

m=size(A,1);         

n=size(A,2);        

for j=1:n

    for i=1:m

        B(i)=A(i,j);

    end

    y=B';

    for i=1:j-1

        r(i,j)=q*B';

        y=y-r(i,j)*q';

    end

    r(j,j)=norm(y,2);

    for i=1:m

        q(i)=y(i)/r(j,j);

        Q(i,j)=q(i);

    end

end

disp ('Q='); disp(Q);

disp ('R='); disp(r);

 

Results:

(a)

>> A=[4 0;3 1];

>> GS(A)

Q=

    0.8000   -0.6000

    0.6000    0.8000

 

R=

    5.0000    0.6000

         0    0.8000

 

(b)

>> A=[1 2;1 1];

>> GS(A)

Q=

    0.7071    0.7071

    0.7071   -0.7071

 

R=

    1.4142    2.1213

         0    0.7071

 

(c)

>> A=[2 1;1 -1;2 1];

>> GS(A)

Q=

    0.6667    0.2357

    0.3333   -0.9428

    0.6667    0.2357

 

R=

    3.0000    1.0000

         0    1.4142

 

(d)

>> A=[4 8 1;0 2 -2;3 6 -7];

>> GS(A)

Q=

    0.8000         0    0.1361

         0    1.0000    0.2722

    0.6000         0   -0.9526

 

R=

    5.0000   10.0000   -2.0000

         0    2.0000   -2.0000

         0         0    7.3485

 

 

Problem 5.

Use the MATLAB QR factorization to find the least squares solutions and 2-norm error of the

following inconsistent systems:

(a)[1 1;2 1;1 2;0 3]*[x1;x2]=[3;5;5;5].

Answer:

>> A=[1 1;2 1;1 2;0 3];

>> y=[3;5;5;5];

>> [Q,R]=qr(A);

>> b=Q'*y

 

b =

 

   -7.3485

    5.4688

    0.2158

    0.2139

 

>> c=Rb

 

c =

 

    1.6154

    1.6615

 

>> e=sqrt(0.2158.^2+0.2139.^2)

 

e =

 

    0.3038

 

 

4.4 Exercises

问题:3.Let A=[1 0 a13;0 1 a23;0 0 1].Prove that for any x0 and b,GMRES converges to the exact solution after two steps.

 

 

4.4 Computer Problems

1.Let A be the n*n matrix with n=1000 and entries A(i,i)=i,A(i,i+1)=A(i+1,i)=1/2,A(i,i+2)=1/2 for all i that fit within the matrix.(a) Print  the zeros structure spy(A).(b) Let xe be the vector of ones.Set  b=Axe,and apply the conjugate Gradient Method,without preconditioner,with the Jacobi preconditioner,and with the Gauss-Seidel preconditioner.Compare errors of the  three runs in a plot versus step number.

Answer:

(a):
Code:

function A=P(n)
for i=1:n-2
    A(i,i)=i;
    A(i,i+1)=1/2;
    A(i+1,i)=1/2;
    A(i,i+2)=1/2;
    A(i+2,i)=1/2;
end

Results:

>> A=P(1000);
>> spy(A)

(b):
Code

1

function A=P(n)
for i=1:n-2
    A(i,i)=i;
    A(i,i+1)=1/2;
    A(i+1,i)=1/2;
    A(i,i+2)=1/2;
    A(i+2,i)=1/2;
end

2
function xc=CG(x0,A,n)
m=size(A,1);
x=ones(m,1);
b=A*x;
x1=x0;
r1=b-A*x1;
d1=r1;
for k=1:n
    if r1==0
        xc=x1;
    end
    a1=(r1'*r1)/(d1'*A*d1);
    x2=x1+a1*d1;
    r2=r1-a1*A*d1;
    b1=(r2'*r2)/(r1'*r1);
    d2=r2+b1*d1;
    r1=r2;
    x1=x2;
    d1=d2;
end
xc=x1;

3
function xc=JCG(x0,A,n)
m=size(A,1);
M=zeros(m,m);
for i=1:m
    M(i,i)=A(i,i);
end
x=ones(m,1);
b=A*x;
x1=x0;
r1=b-A*x1;
z1=inv(M)*r1;
d1=z1;
for k=1:n
    if r1==0
        xc=x1;
    end
    a1=(r1'*z1)/(d1'*A*d1);
    x2=x1+a1*d1;
    r2=r1-a1*A*d1;
    z2=inv(M)*r2;
    b1=(r2'*z2)/(r1'*z1);
    d2=z2+b1*d1;
    r1=r2;
    z1=z2;
    x1=x2;
    d1=d2;
end
xc=x1;

4
function xc=GCG(x0,A,n)
m=size(A,1);
M=zeros(m,m);
L=M;
U=M;
D=M;
for i=1:m
   for j=1:m
       if i>j
           L(i,j)=A(i,j);
       end
       if i<j
           U(i,j)=A(i,j);
       end
       if i==j
           D(i,j)=A(i,j);
       end
   end
end
M=(D+L)*D.^(-1)*(D+U);          
x=ones(m,1);
b=A*x;
x1=x0;
r1=b-A*x1;
z1=M.^(-1)*r1;
d1=z1;
for k=1:n
    if r1==0
        xc=x1;
    end
    a1=(r1'*z1)/(d1'*A*d1);
    x2=x1+a1*d1;
    r2=r1-a1*A*d1;
    z2=M.^(-1)*r2;
    b1=(r2'*z2)/(r1'*z1);
    d2=z2+b1*d1;
    r1=r2;
    z1=z2;
    x1=x2;
    d1=d2;
end
xc=x1;

Results:

>> A=P(1000);
>> x0=zeros(1000,1);
>> xc=CG(x0,A,50);

>> A=P(1000);
>> x0=zeros(1000,1);
>> xc=JCG(x0,A,30);

>> A=P(1000);
>> x0=zeros(1000,1);
>> xc=GCG(x0,A,30);

 

 

2.Let n=1000.Start with the n*n matrix A from computer Problem 1,and add the nonzero entries A(i,2i)=A(2i,i)=1/2 for 1<=i<=n/2.Carry out steps (a)and(b) as in that problem.

 

function [Q R]=Householder(A)

 [m,n]=size(A);H1=eye(m);

 for i=1:n
   H=eye(m);
   b=A(i:m,i); 
   y=norm(b,2);
   B=eye(m);B(i,i)=y;
   w=B(i:m,i);
   v=w-b;
   P=(v*v')/(v'*v);
   H(i:m,i:m)=eye(m+1-i)-2*P;
   A=H*A;
   H1=H1*H;
   
end
Q=H1;
R=A;
 (3)Result:
>> A=[4 2 3 0;-2 3 -1 1;1 3 -4 2;1 0 1 -1;3 1 3 -2];
>> b=[10;0;2;0;5];
>> [Q,R]=Householder(A)
 Q =
    0.7184    0.2115    0.2259    0.6086   -0.1332
   -0.3592    0.7685    0.5228   -0.0516   -0.0666
    0.1796    0.5993   -0.7688   -0.1320    0.0133
    0.1796   -0.0564    0.0525   -0.4071   -0.8922
0.5388    0.0494    0.2862   -0.6662    0.4261

 

 

 


R =

    5.5678    1.4368    3.5921   -1.2572
   -0.0000    4.5755   -2.4393    1.9247
   -0.0000   -0.0000    4.1408   -1.6395
   -0.0000   -0.0000   -0.0000    1.4237
    0.0000    0.0000    0.0000   -0.0000
>>x=R(Q'*b)
x =
    1.2739
    0.6885
    1.2124
    1.7497
>> error=norm(b-A*x)
error =
0.8256
So, 2-norm error of the inconsistent system is 0.8256.

 







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

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

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