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

17074216作业8

作者:17074216丁龙   发布时间:2019-06-09 21:40:52   浏览次数:94

4.3EXERCISE

2. Apply classical Gram–Schmidt orthogonalization to fifind the full QR factorization of the

following matrices:

 

(b)A=[-4 -4;-2 7;4 -5]
ANSWER:

We need to find Householder reflectors.
Let x=[-4;-2;4]and w=[||x||2;0;0]=[6;0;0]
Let v=w-x=[10;2;-4]
And P=v*v/(v*v)
P =
    5/6    1/6   -1/3
    1/6    1/15   -1/15
    -1/3   -1/15    2/15

H1=I-2P=
        -2/3   -1/3     2/3
        -1/3    14/15   2/15
         2/3    2/15    11/15

H1*A=
6    -3
0 108/15
0    -81/15

x=             w=             v=
108/15          9              27/15
-81/15           0              81/15

P=
   1/10   3/10
   3/10   9/10

H2=
4/5    -3/5
-3/5    -4/5

Q=H1*H2=
            -2/3    -2/3   -1/3
            -1/3     2/3   -2/3
             2/3    -1/3   -2/3
R=
6   -3
0 9
0 0

4.Apply modifified Gram–Schmidt orthogonalization to fifind 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 reflflectors to fifind the full QR factorization of the matrices in Exercise 2

Answer:

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

p=(v*v')/(v'*v)=1/6*[1 2 -1;2 4 -2;-1 -2 1],

H1=[1 0 0;0 1 0;0 0 1]-2/6*[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],

H1*A=[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],

x1=[0;-3], w1=[3;0], v1=w1-x1=[3;0]-[0;-3]=[3;3],

p1=(v1*v1')/(v1'*v1)=1/18*[9 9;9 9]=1/2*[1 1;1 1],

H2=I-2*p1=[1 0;0 1]-[1 1;1 1]=[0 -1;-1 0],

H2*x1=[0 -1;-1 0]*[0;-3]=[3;0]=w2,

H2*H1*A=[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,

[2 3;-2 -6;1 0]=H1*H2*R=[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 2/3;-2/3 -2/3 1/3;1/3 -2/3 -2/3]*[3 6;0 3;0 0]=Q*R.

 

(b):

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

p=(v*v')/(v'*v)=1/120*[100 20 -40;20 4 -8;-40 -8 16]=1/30*[25 5 -10;5 1 -2;-10 -2 4],

H1=[1 0 0;0 1 0;0 0 1]-2/30*[25 5 -10;5 1 -2;-10 -2 4]=[-2/3 -1/3 2/3;-1/3 14/15 2/15;2/3 2/15 11/15],

H1*A=[-2/3 -1/3 2/3;-1/3 14/15 2/15;2/3 2/15 11/15]*[-4 -4;-2 7;4 -5]=[6 -3;0 36/5;0 -27/5],

x1=[36/5;-27/5], w1=[9;0], v1=w1-x1=[9;0]-[36/5;-27/5]=[9/5;27/5],

p1=(v1*v1')/(v1'*v1)=1/2*[1/5 3/5;3/5 9/5],

H2=I-2*p1=[1 0;0 1]-[1/5 3/5;3/5 9/5]=[4/5 -3/5;-3/5 -4/5],

H2*x1=[4/5 -3/5;-3/5 -4/5]*[36/5;-27/5]=[9;0]=w1,

H2*H1*A=[1 0 0;0 4/5 -3/5;0 -3/5 -4/5]*[-2/3 -1/3 2/3;-1/3 14/15 2/15;2/3 2/15 11/15]*[-4 -4;-2 7;4 -5]=[6 -3;0 9;0 0]=R,

[-4 -4;-2 7;4 -5]=H1*H2*R=[-2/3 -1/3 2/3;-1/3 14/15 2/15;2/3 2/15 11/15]*[1 0 0;0 4/5 -3/5;0 -3/5 -4/5]*[6 -3;0 9;0 0]

                 =[-2/3 -2/3 -1/3;-1/3 2/3 -2/3;2/3 -1/3 -2/3]*[6 -3;0 9;0 0]=Q*R.

 

 

 



13. Let P be the matrix defifined in (4.29). Show (a) P

2 = P (b) P is symmetric (c) P v = v.

Answer:

(a)
P^2=v*v’*v*v’/(v’*v*v’*v)
Due to the result of v’*v is a real number, hence, v*v’*v*v’= v’*v*v*v’
So, we can prove P^2=P.

(b)
(v*v’)’=(v’)’*v’=v*v’
S0, P is symmetric.

(c)
Due to the result of v’*v is a real number, hence, v*v’*v= v’*v*v
So, P*v=v.



function [Q R]=Householder(A)
[m,n]=size(A);H1=eye(m);B=A;
for i=1:n
   H=eye(m);
   x=A(i:m,i); 
   y=norm(x,2);
   B=eye(m);B(i,i)=y;
   w=B(i:m,i);
   v=w-x;
   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;




4.3computer problem
1. Write a Matlab program that implements classical Gram–Schmidt to fifind the reduced QR

factorization. Check your work by comparing factorizations 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


Answer:

function QR=CGS(A)

 

[m n]=size(A);

 

for j=1:n

 

        y=A(:,j)

 

    for i=1:j-1

 

            R(i,j)=Q(:,i)'*A(:,j)

 

            y=y-R(i,j)*Q(:,i)

 

    end

 

    R(j,j)=norm(y,2)

 

    Q(:,j)=y/R(j,j)

 

end


(a)

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

>> QR=CGS(A)

R =

    5.0000    0.6000

         0    0.8000

Q =

    0.8000   -0.6000

        0.6000    0.8000

 

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

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

 

Q =

 

   -0.8000   -0.6000

   -0.6000    0.8000

 

 

R =

 

   -5.0000   -0.6000

         0    0.8000
(b)

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

>> QR=CGS(A)

R =

    1.4142    2.1213

         0    0.7071

Q =

    0.7071    0.7071

0.7071   -0.7071

 

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

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

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];

>> QR=CGS(A)

R =

    3.0000    1.0000

         0    1.4142

Q =

    0.6667    0.2357

    0.3333   -0.9428

0.6667    0.2357

 

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

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

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];

>>  QR=CGS(A)

R =

     5    10     5

     0     2    -2

     0     0     5

Q =

    0.8000         0   -0.6000

         0    1.0000         0

0.6000         0    0.8000

 

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

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

Q =

   -0.8000         0   -0.6000

         0    1.0000         0

   -0.6000         0    0.8000

R =

    -5   -10    -5

     0     2    -2

     0     0     5

5. Use the Matlab QR factorization to fifind the least squares solutions and 2-norm error of the

following inconsistent systems:

Answer:

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

a)The least squares solutions is x=[1.6154 1.6615].

    Code: A=[1 1;2 1;1 2;0 3];

        >> y=[3 5 5 5];

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

        >> c=Q'*y'

 

c =

 

   -7.3485

    5.4688

    0.2158

    0.2139

 

     >> x=R(1:4,1:2)c(1:4)

 

x =

 

    1.6154

    1.6615

Result: n=norm(y'-A*x,2)

 

n =

 

    0.3038

 

4.4computer problem

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) = A(i + 2, i) = 1/2 for all i that fifitwithin the matrix. (a) Print the nonzero structure spy(A). (b) Let xe be the vector of n 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)

function matrixA(n)

S1=sparse(1:n,1:n,1:n);

E=sparse(2:n,1:n-1,(1/2)*ones(1,n-1),n,n);

F=sparse(3:n,1:n-2,(1/2)*ones(1,n-2),n,n);

A=S1+E+E'+F+F';

spy(A)

title('The nonzero structure spy(A)','color','b')

end 

>> matrixA(1000)

 

(b)function matrix4(n)

S1=sparse(1:n,1:n,1:n);

E=sparse(2:n,1:n-1,(1/2)*ones(1,n-1),n,n);

F=sparse(3:n,1:n-2,(1/2)*ones(1,n-2),n,n);

A=S1+E+E'+F+F';

xe=ones(n,1);

b=A*xe;

[x0,f10,rr0,it0,rv0]=pcg(A,b,1e-8,100);

[itr rnorm]=cg_it(A,b,1e-8,100);

[j y]=jacobi(A,b,100);

semilogy(0:it0,rnorm,'bo-',0:it0,rv0/norm(b),'rs-',1:j,y,'kd--',...'LineWidth',1);

xlabel('Step number'));

ylabel('Relative residual');

legend('The Conjugate Gradient Method','The Jacobi preconditioner',...'Gauss-Seidel Method');

    title('Efficiency of Preconditioned Conjugate Gradient Method'))

end

>> matrix4(100)

 

 

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.

 

answer:

The results showed that the GMRES without preconditions were slow to converge, and the use of the Jacobi precondition was improved greatly, and the GMRES with gauss-seidel preconditions were rapidly convergent.

 

Code:

(a)

function [A]=Spy(n)

a=zeros(n);

for i=1:n/2

    t=2*i;

    a1=find(a(i,t)); 

    a1=1/2; 

    a2=find(a(t,i)); 

    a2=1/2;

end

a1=(ones(1,n-1)*1/2);

a2=(1:n);

A=diag(a1,1)+diag(a1,-1)+diag(a2)+a;

spy(A)

 

test:

clear

clc

n=1000;

[A]=Spy(n)

 

(b)

function [xk]=GMRES(A,b)

r=(M.^(-1))*(b-A*x0);

q(:,1)=r/norm(r);

for k=1:m

    w=(M.^(-1))*A*q(:,k);

    for j=1:k

        h(j,k)=w'*q(:,j);

        w=w-h(j,k)*q(:,j);

    end

    h(k+1,k)=norm(w);

    q(:,k+1)=w/h(k+1,k);

    ck1=[norm(r) zeros(1,length(r)-1)];

    ck=norm(h(:,1:ck)-ck1');

    xk=q(:,1:ck)+x0;

end

 







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

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

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