﻿ 17074221 黄剑锋 作业八 4.3 4.4 - 计算数学达人 - 专，学者，数值代数，微分方程数值解

### 17074221 黄剑锋 作业八 4.3 4.4

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.

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

(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