﻿ 17074223-陈涛涛-作业八：4.3 and4.4 - 计算数学达人 - 专，学者，数值代数，微分方程数值解

### 17074223-陈涛涛-作业八：4.3 and4.4

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 ﬁt 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题程序同上