﻿ 17074233 李蒙蒙 作业8 - 计算数学达人 - 专，学者，数值代数，微分方程数值解

# 4.3 Exercises

2. Apply classical Gram–Schmidt orthogonalization to find the full QR factorization of the following matrices:

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

Sol:

y1=A1=[-4;-2;4].Then r11 = ||y1||2 =6,

q1=y1/||y1||2=[-2/3;-1/3;2/3]

y2=A2-q1q1’ A2=[-6;6;-3]

q2=y2/||y2||2=[-2/3;2/3;-1/3]

r12=q1’A2=-3 and r22=||y2||2=9

A=[-4,-4;-2,7;4,-5]=[-2/3,-2/3;-1/3,2/3;2/3,-1/3][6,-3;0,9]=QR

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

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

Sol:

y1=A1=[-4;-2;4]

r11=||y1||2=((-4)^2+(-2)^2+4^2)=6

q1=y1/r11=[-2/3;-1/3;2/3]

y2=A2-q1*(q1’*A2)

r12=q1’*A2=[-2/3,-1/3,2/3]*[-4;7;-5]=-3

y2=A2+3*q1=[-4;7;-5]+[-2;-1;2]=[-6;6;-3]

r22=||y2||2=((-6)^2+6^2+(-3)^2)=9

q2=y2/||r22=[-2/3;2/3;-1/3]

Assuming A3=[1;0;0] is linearly independent of A1 and A2, then

y31=A3-q1*(q1’*A3)

r13= r13=q1’*A3=[-2/3,-1/3,2/3]*[1;0;0]=-2/3

y31=A3+2/3*q1=[1;0;0]+2/3*[-2/3;-1/3;2/3]=[5/9;-2/9;4/9]

y3=y31-q2*(q2’*y31)

r23=q2’*y31=[-2/3,2/3,-1/3]*[5/9;-2/9;4/9]=-2/3

y3=y31+2/3*q2=[5/9;-2/9;4/9]+2/3*[-2/3;2/3;-1/3]=[1/9;2/9;2/9]

r33=||y3||2=((1/9)^2+(2/9)^2+(2/9)^2)=1/3

q3=y3/r33=[1/3;2/3;2/3]

The complete QR decomposition of matrix A obtained by modified Gram–Schmidt orthogonal method is as follows:

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

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

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

Sol:

x1=[-4;-2;4]  w1=[||x||2;0;0]=[6;0;0]    v1=w1-x1=[10;2;-4]

H1=I-2v1v1’/v1’v1=[-2/3,-1/3,2/3;-1/3,14/15,2/15;2/3,2/15,11/15]

H1A=[6,-3;0,36/5;0,-27/5]

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

H2^=I-2v2v2’/v2’v2=[4/5,-3/5;-3/5,-4/5]

H2H1A=[6,-3;0,9;0,0]

A=[-4,-4;-2,7;4,-5]=H1H2R=[-2/3,-2/3,-1/3;-1/3,2/3,-2/3;2/3,-1/3,-2/3][6,-3;0,9;0,0]

13. Let P be the matrix defined in (4.29). Show (a) P 2 = P (b) P is symmetric

(c) Pv = v.

Sol

(a):P^2=P*P=((v*v’)/(v’*v))* ((v*v’)/(v’*v))=(v*v’)/(v’*v)=P

(b):P’=((v*v’)/(v’*v))’=(v*v’)/(v’*v)=P

(c):Pv=((v*v’)*v)/(v’*v)=v

# 4.3 Computer Problems

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.

(a)[4,0;3,1]   (b) [1,2;1,1]  (c) [2,1;1,-1;2,1]  (d) [4,8,1;0,2,-2;3,6,7]

Sol:

Code:

function [q,r]=QR(A)

[m,n]=size(A);

for j=1:n

y=A(1:m,j);

for i=1:j-1

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

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

end

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

q(1:m,j)=y/r(j,j)

end

(a)

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

>> [q,r]=QR(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

(b)

>> A=[1,2;1,1];

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

q =    0.7071    0.7071

0.7071   -0.7071

r = 1.4142    2.1213

0    0.7071

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

Q =   -0.7071   -0.7071

-0.7071    0.7071

R =  -1.4142   -2.1213

0   -0.7071

(c)

>> A=[2,1;1,-1;2,1];

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

q =    0.6667    0.2357

0.3333   -0.9428

0.6667    0.2357

r =    3.0000    1.0000

0    1.4142

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

Q =   -0.6667    0.2357

-0.3333   -0.9428

-0.6667    0.2357

R =   -3.0000   -1.0000

0    1.4142

(d)

>> A=[4,8,1;0,2,-2;3,6,7];

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

q =

0.8000         0   -0.6000

0    1.0000         0

0.6000         0    0.8000

r =

5    10     5

0     2    -2

0     0     5

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

Q =

-0.8000         0   -0.6000

0    1.0000         0

-0.6000         0    0.8000

R =

-5   -10    -5

0     2    -2

0     0     5

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]

Sol:

Code:

function [q,r,x]=QR_LS(A,b)

[m,n]=size(A);

for j=1:n

y=A(1:m,j);

for i=1:j-1

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

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

end

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

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

end

r1=r(1:n,1:n);

d=q(1:m,1:n)'*b;

x=r1d;

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

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

>> [q,r,x]=QR_LS(A,b)

q =

0.4082    0.0506

0.8165   -0.2025

0.4082    0.3545

0    0.9115

r =

2.4495    2.0412

0    3.2914

x =

1.6154

1.6615

# 4.4 Exercises

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

Report all approximations x k up to and including the correct solution.

(c)

0  0  1

1  0  0

0  1  0

Sol:

x0=[0,0,0]

r=b-Ar0=[1;0;0]-[0,0,1;1,0,0;0,1,0]*[0;0;0]=[1;0;0]

q1=r/||r||2=[1;0;0]

y1=A*q1=[0,0,1;1,0,0;0,1,0]*[1;0;0]=[0;1;0]

h11=q1T*y=[1,0,0]*[0;1;0]=0

y2=y1-h11*q1=[0;1;0]

h21=||y2||2=1

q2=[0;1;0]

minimize||Hc1-[||r||2,0]||2,get c1=0

x1=Q1*c1+x0=[0,0,0]

y3=A*q2=[0,0,1;1,0,0;0,1,0]*[0;1;0]=[0;0;1]

h12=q1T*y3=0

y4=y3-h12*q1=[0;0;1]

h22=q2T*y4=0

y5=[0;0;1]

h32=||y5||2=1

q3=y5/h32=[0;0;1]

minimize||Hc2-[||r||2,0]||2,get c2=[0;0]

x2=Q2*c2+x0=[0,0,0]

y6=A*q3=[0,0,1;1,0,0;0,1,0]*[0;0;1]=[1;0;0]

h13=q1T*y6=[1,0,0]*[1;0;0]=1

y7=y6-h13*q1=[0;0;0]

h23=q2T*y7=0

y8=y7-h23*q2=[0;0;0]

h33=q3T*y8=0

y9=y8-h33*q3=[0;0;0]

h43=||y9||2=0

minimize||Hc3-[||r||2,0]||2,get c3=[0;0;1]

x3=Q3*c3+x0=[0;0;1]

3. Let A =

1  0   a13

0  1   a23

0  0    1

Prove that for any x 0 and b, GMRES converges to the exact

solution after two steps.

Pf:

Suppose x0 and b are arbitrary vector.

r=b-A*x0

q1=r/||r||2

y1=A*q1

h11=q1T*y1

y2=y1-h11*q1=A*q1-h11*q1

h21=||y2||2

q2=r2/h21

y3=A*q2=A*(A*q1-h11*q1)/||y2||2

h12=q1T*y3

y4=y3-h12*q1= A*(A*q1-h11*q1)/||y2||2-h12*q1

h22=q2T*y4

y5=y4-h22*q2

h32=||y5||2

so title prove that h32=0,or that y5 is the zero vector,which is

A*(A*q1-h11*q1)/||y2||2-h12*q1=h22*q2

# 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 fitwithin the matrix. (a) Print the nonzero structure spy(A). (b) Let x e be the vector of n ones.Set b = Ax e , 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.

function [A]=Spy(n)

for i=1:1000

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

spy(A)

function x=CG(A,b)

n=length(b);

x0=zeros(n,1);

r0=b-A*x0;

d0=r0;

for k=0:(n-1)

a=(r0'*r0)/(d0'*A*d0);

x=x0+a*d0;

rk=r0-a*A*d0;

if rk==0;

break;

end

b0=(rk'*rk)/(r0'*r0);

dk=rk+b0*d0;

x0=x;

r0=rk;

d0=dk;

end

function [x,k]=Jacobi(A,b,nmax,tol)

n=length(b);           % find n

d=diag(A);             % extract diagonal of a

r=A-diag(d);           % r is the remainder

x=zeros(n,1);          % initialize vector x

k=0;

for j=1:nmax

x=(b-r*x)./d;

k=k+1;

end

>>x=CG(A,A*ones(n,1))

X=ones(1000,1)

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.

function [A]=Spy(n)

for i=1:1000

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;

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

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

end

spy(A)