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

17074223-陈涛涛-作业八:4.3 and4.4

作者:17074223   发布时间:2019-06-09 21:24:38   浏览次数:103

4.3Exercises

2.Apply classical GramSchmidt orthogonalization to nd the full QR factorization of the following matrices:

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

Q =

 

   -0.6667   -0.6667

   -0.3333    0.6667

0.6667   -0.3333

R =

 

    6.0000   -3.0000

         0    9.0000

 

4.    Apply modied GramSchmidt orthogonalization to nd the full QR factorization of the matrices in Exercise 2.

Q =

 

   -0.6667   -0.6667

   -0.3333    0.6667

    0.6667   -0.3333

 

R =

 

    6.0000   -3.0000

         0    9.0000

 

6.Apply Householder reectors to nd the full QR factorization of the matrices in Exercise 2.

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

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

13. Let P be the matrix dened in (4.29). Show (a) P 2 = P (b) P is symmetric (c) Pv = v.

a)P=v*v'/v'*v

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

So P is symmetric

(c)Pv=(v*v'/v'*v)*v=1/v'v*(v'*v*v)=(v'*v/v'*v)*v=v

Computer problem

1.    Write a Matlab program that implements classical GramSchmidt to nd 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.

Code:

function [Q,R]=CGram(A)

n=size(A,2);

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=[4 0;3 1]

>> [Q,R]=CGram(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

 

5.Use the Matlab QR factorization to nd 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]

function [Q,R]=Householder(A)

[m,n]=size(A);

R=A;

Q=eye(m);

for j=1:n

    x=R(j:m,j);

    w=zeros(m-j+1,1);

    w(1,1)=norm(x);

    v=w-x;

    H=eye(m);

    H(j:m,j:m)=eye(m-j+1)-2*(v*v')/(v'*v);

    R=H*R;

    Q=Q*H';

end

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

>> b=[3;5;5;5];

>> [Q,R]=Householder(A)

Q =

 

    0.4082    0.0506    0.6715    0.6163

    0.8165   -0.2025   -0.5372    0.0611

    0.4082    0.3545    0.4028   -0.7385

         0    0.9115   -0.3133    0.2665

 

 

R =

 

    2.4495    2.0412

   -0.0000    3.2914

   -0.0000   -0.0000

   -0.0000    0.0000

>> err=norm(b-A.*x)

 

err =

 

  13.142073481004230

4.4EXERCISES

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]

(c) x1 = [0.0332, 0.8505, 0.9668],x2 = [0.0672, 0.8479, 0.9696],x3 = [0, 0, 1]

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

Solution:

Code:

(a) A=zeros(1000);

for i=1:1000

    A(i,i)=i;

end

for i=1:999

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

end

for i=1:998

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

end

(b) function x=cg (A,x0)

n=size(A,1);

x=ones(n,1);

b=A*x;

r=b-A*x0;

d=r;

for k=1:n

    if r==0

        x=x0;

    end

 a=(r'*r)/(d'*A*d);

 x=x+a*d;

 R=r-a*A*d;

 B=(R'*R)/(r'*r);

 d=R+B*d;

 r=R;

end

 

function x=jcg (A,x0)

n=size(A,1);

x=ones(n,1);

b=A*x;

r=b-A*x0;

m=diag(diag(A));

d=mr;

z=mr;

for k=0:n-1

    if r==0,stop,end

 a=(r'*z)/(d'*A*d);

 x=x+a*d;

 R=r-a*A*d;

 Z=mR;

 B=(R'*Z)/(r'*z);

 d=Z+B*d;

  z=mR;

 r=R;

end

 

function x=gcg (A,x0)

n=size(A,1);

x=ones(n,1);

b=A*x;

r=b-A*x0;

D=diag(diag(A));

L=tril(A,-1);

U=tril(A,1);

m=(D+L)/D(D+U);

d=mr;

z=mr;

for k=0:n-1

    if r==0,stop,end

 a=(r'*z)/(d'*A*d);

 x=x+a*d;

 R=r-a*A*d;

 Z=mR;

 B=(R'*Z)/(r'*z);

 d=Z+B*d;

  z=mR;

 r=R;

end

solution:

(a)    >> spy(A)

(b)    >>x0=zeros(n,1)

>> x=cg (A,x0)

 

>>x0=zeros(n,1)

>> x=jcg (A,x0)

 

>>x0=zeros(n,1)

>> x=gcg (A,x0)

 

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 1i n/2. Carry out steps (a) and (b) as in that problem.

(a) A=zeros(1000);

for i=1:1000

    A(i,i)=i;

end

for i=1:999

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

end

for i=1:998       

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

end

for i=1:500

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

End

(b)b题程序同上

 

 

 







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

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

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