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

17074197 数值分析第八次作业

作者:   发布时间:2019-06-09 20:57:22   浏览次数:98

 P224

Exercise 

Problem 2.

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

 

Solution :

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(q1’A2)=[-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=q1’A2=-3 , r22=||y2||2=9

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

y3=A3-q1(q1’A3)-q2(q2’A3)

   =[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.

 

 

 

Problem 4.

Apply modified Gram–Schmidt orthogonalization to find the full QR factorization of the matrices in Exercise 2. 

 

Solution :

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(q1’A2)=[-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=q1’A2=-3 , r22=||y2||2=9

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

y3(1)=[1 ; 0 ; 0]- [-2/3 ; -1/3 ; 2/3 ]q1’A3=[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.

 

Problem 6.

Apply Householder reflectors to find the full QR factorization of the matrices in Exercise 2. 

 

Solution :

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

x1=[-4 -2 4]';

w1=[norm(x1) 0 0]'=[6 0 0]';

v1=w1-x1=[10 2 -4];

P1=(v1*v1')/(v1'*v1)=[5/6 1/6 -1/3;1/6 1/30 -1/15;-1/3 -1/15 2/15]; 

H1=I-2*P1=[1 0 0;0 1 0;0 0 1]-2*P1=[-2/3 -1/3 2/3;-1/3 14/15 2/15;2/3 2/15 11/15];     

H1*A=[6 -3;0 36/5;0 -27/5];

x2=[36/5 -27/5]';

w2=[norm(x2) 0]'=[9 0]';

v2=w2-x2=[9/5 27/5]';

P2=(v2*v2')/(v2'*v2)=[1/10 3/10;3/10 9/10];  

H2'=I-2*P2=[1 0;0 1]-2*P2=[4/5 -3/5;-3/5 -4/5];

H2=[1 0 0;0 4/5 -3/5;0 -3/5 -4/5];

R=H2*H1*A=[6 -3;0 9;0 0];

Q=H1*H2=[-2/3 -2/3  -1/3;-1/3 2/3  -2/3;2/3 -1/3 -2/3]. 

 

  

Problem 13.

Let P be the matrix defined in (4.29). Show (a) P2 = P (b) P is symmetric (c) Pv = v. 

 

Solution :

(a)

P.^2=[(v*v')/(v'*v)]*[(v*v')/(v'*v)]=(v*v'*v*v')/(v'*v).^2=(v*v')/(v'*v)=P.

 

(b)

P'=(v*v')'/(v'*v)=((v')'*v')/(v'*v)=(v*v')/(v'*v)=P.

 

(c)

Pv=((v*v')/(v'*v))*v=(v*v'*v)/(v'*v)=v.

 

Computer problem

Problem 1.

Write a Matlab program that implements classical Gram–Schmidt to find 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. 

 

Solution :

function CGS(A)

n=size(A,2);         

m=size(A,1);        

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

 

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

>> CGS(A)

Q=

    0.8000   -0.6000

    0.6000    0.8000

 

R=

    5.0000    0.6000

         0    0.8000

 

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

 

Q =

 

   -0.8000   -0.6000

   -0.6000    0.8000

 

 

R =

 

   -5.0000   -0.6000

         0    0.8000

 

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]
(a)

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

 

P229

Exercise

Problem 1.

Solve Ax = b for the following A and b = [1,0,0]T , using GMRES with x0 = [0,0,0]T . Report all

approximations xk up to and including the correct solution. 

(c)[0 0 1;1 0 0;0 1 0]

 

 

Problem 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. 

 

Computer problem

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 fit within 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. 

 

(a)

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

>> A=P(1000);

>> spy(A)

 

(b)

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

 

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;

 

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;

 

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

 

 

Problem 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. 

(a)

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

for i=1:n/2

    A(i,2*i)=1/2;

    A(2*i,i)=1/2;

end

 

(b) 

The same as Problem 1.

 

 

 







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

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

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