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

17074221 黄剑锋 作业八 4.3 4.4

作者:17074221   发布时间:2019-06-09 21:38:44   浏览次数:108

 4.3

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

   -0.3333    0.9397

   -0.6667   -0.2921

r =

 

    6.0000    3.6667

         0    8.7496

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

q =

 

   -0.6667   -0.1778

   -0.3333    0.9397

   -0.6667   -0.2921

r =

 

    6.0000    3.6667

         0    8.7496

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.

function [q,r]=cgs(A)

[n,m]=size(A);

for j=1:m

    y=A(1:n,j);

    for i=1:(j-1)

        r(i,j)=(q(1:n,i))'*A(1:n,j);

        y=y-r(i,j)*q(1:n,i)

    end

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

    q(1:n,j)=y/r(j,j);

    A(1:n,j)

end

>> A=[4 0;3 1]

 

A =

 

     4     0

     3     1

 

>> [q,r]=cgs(A)

 

ans =

 

     4

     3

 

 

y =

 

  -0.480000000000000

   0.640000000000000

 

 

ans =

 

     0

     1

 

 

q =

 

   0.800000000000000  -0.600000000000000

   0.600000000000000   0.800000000000000

 

 

r =

 

   5.000000000000000   0.600000000000000

                   0   0.800000000000000

 

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

 

Q =

 

  -0.800000000000000  -0.600000000000000

  -0.600000000000000   0.800000000000000

 

 

R =

 

  -5.000000000000000  -0.600000000000000

                   0   0.800000000000000

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)

[r,c]=size(A);

R=A;

Q=eye(r);

for j=1:c

    x=R(j:r,j);

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

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

    v=w-x;

    H=eye(r);

    H(j:r,j:r)=eye(r-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.408248290463863   0.050636968354183   0.671535958632810   0.616286205453556

   0.816496580927726  -0.202547873416733  -0.537189441417618   0.061116252642976

   0.408248290463863   0.354458778479283   0.402842924202425  -0.738518710739509

                   0   0.911465430375300  -0.313344121873347   0.266544987794162

 

 

R =

 

   2.449489742783178   2.041241452319315

  -0.000000000000000   3.291402943021916

  -0.000000000000000  -0.000000000000000

  -0.000000000000000   0.000000000000000

>> x=Q'*b

 

x =

 

   7.348469228349534

   5.468792582251799

  -0.223845319544270

  -0.205428735151186

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

 

err =

 

  13.1420734810042304.4

 

4.4

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.

 

 

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=A*xe,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);

 

Problem 2:

Let n=1000.Start with the n*n matrix A from Computer Problem 1,and add the nonzero entries A(i,2*i)=A(2*i,i)=1/2 for

1<=i<=n/2.Carry out steps (a) and (b) as in that problem.

 

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

for i=1:n/2

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

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

end

 

Results:

(b):

Code and Results:

the same as problem 1







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

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

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