﻿ 17074198 作业八 - 计算数学达人 - 专，学者，数值代数，微分方程数值解

### 17074198 作业八

4.3Exercises

2. Problem:Apply classical GramSchmidt orthogonalization to find the full QR factorization of the following matrices: (b)[-4,-4;-2,7;4,-5]

Solution:

(b)Set y1=A1=[-4;-2;4].Then r11=||y1||2=sqrt((-4)^2+(-2)^2+4^2)=6,and the first unit vector is q1=y1/||y1||2=[-2/3;-1/3;2/3].

To find the second unit vector,set

y2=A2-q1*q1*A2=[-4;7;-5]-[-2/3;-1/3;2/3]*(-3)=[-6;6;-3] and

q2=y2/||y2||2=(1/9)*[-6;6;-3]=[-2/3;2/3;-1/3].

r12=q1*A2=-3,r22=||y2||2=9.

y3=A3-q1*q1*A3-q2*q2*A3=[1;0;0]-[-2/3;-1/3;2/3]*(-2/3)-[-2/3;2/3;-1/3]*(-2/3)

=[1/9;2/9;2/9] and

q3=y3/||y3||2=[1/9;2/9;2/9]*3=[1/3;2/3;2/3].

Putting the parts together,we obtain the full QR factorization

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

[-1/3 2/3 2/3][0 9]

[2/3 -1/3 2/3][0 0]

4. Problem:Apply modified GramSchmidt orthogonalization to find the full QR factorization of the matrices in Exercise 2.

Solution:

(a) [2,3;-2,-6;1,0]

Set y1=A1=[2;-2;1].Then r11=||y1||2=sqrt(2^2+(-2)^2+1^2)=3,and the first unit vector is q1=y1/||y1||2=[2/3;-2/3;1/3].

To find the second unit vector,set

y2=A2-q1*q1*A2=[3;-6;0]-[2/3;-2/3;1/3]*6=[-1;-2;-2] and

q2=y2/||y2||2=(1/3)*[-1;-2;-2]=[-1/3;-2/3;-2/3].

r12=q1*A2=6,r22=||y2||2=3.

y3=A3-q1*q1*A3=[1;0;0]-[2/3;-2/3;1/3]*(2/3)=[5/9;4/9;-2/9],

y3=y3-q2*q2y3=[5/9;4/9;-2/9]-[-1/3;-2/3;-2/3]*(-1/3)=[3/9;2/9;-4/9]

and q3=y3/||y3||2=9/(29)*[3/9;2/9;-4/9]=[3/29;2/29;-4/29].

Putting the parts together,we obtain the full QR factorization

A=[2,3;-2,-6;1,0]=[2/3 -1/3 3/29][3 6]=QR

[-2/3 -2/3 2/29][0 3]

[1/3 -2/3 -4/29][0 0]

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

Set y1=A1=[-4;-2;4].Then r11=||y1||2=sqrt((-4)^2+(-2)^2+4^2)=6,and the first unit vector is q1=y1/||y1||2=[-2/3;-1/3;2/3].

To find the second unit vector,set

y2=A2-q1*q1*A2=[-4;7;-5]-[-2/3;-1/3;2/3]*(-3)=[-6;6;-3] and

q2=y2/||y2||2=(1/9)*[-6;6;-3]=[-2/3;2/3;-1/3].

r12=q1*A2=-3,r22=||y2||2=9.

y3=A3-q1*q1*A3=[1;0;0]-[-2/3;-1/3;2/3]*(-2/3)=[5/9;-2/9;4/9]

y3=y3-q2*q2y3=[5/9;-2/9;4/9]-[-2/3;2/3;-1/3]*(-2/3)=[1/9;2/9;2/9]

and q3=y3/||y3||2=3*[1/9;2/9;2/9]=[1/3;2/3;2/3].

Putting the parts together,we obtain the full QR factorization

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

[-1/3 2/3 2/3][0 9]

[2/3 -1/3 2/3][0 0]

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

Solution:

(a)We need to find a Householder reflector that moves the first column x = [2,-2,1] to

the vector w = [||x||2,0,0].Set v = w x = [3,0,0] [2,-2,1] = [1,2,-1]. Referring

to Theorem 4.4, we have

H1=I-2v*v/v*v=[1,0,0;0,1,0;0,0,1]-2[1,2,-1][1;2;-1]/[1;2;-1][1,2,-1]

=[1,0,0;0,1,0;0,0,1]-(1/3)*[1,2,-1;2,4,-2;-1,-2,1]

=[2/3,-2/3,1/3;-2/3,-1/3,2/3;1/3,2/3,2/3]

and H1*A=[2/3,-2/3,1/3;-2/3,-1/3,2/3;1/3,2/3,2/3][2,3;-2,-6;1,0]

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

The remaining step is to move the vector xˆ = [0,-3] to wˆ = [3,0]. Calculating H^2 from Theorem 4.4 yields

H^2=[1,0;0,1]-2[3,3][3;3]/[3;3][3,3]

=[1,0;0,1]-(1/9)*[9,9;9,9]

=[0,-1;-1,0]

H^2*x^=[0,-1;-1,0][0,-3]=[3;0]

H2*H1*A=[1,0,0;0,0,-1;0,-1,0][2/3,-2/3,1/3;-2/3,-1/3,2/3;1/3,2/3,2/3][2,3;-2,-6;1,0]

=[2/3,-2/3,1/3;-1/3,-2/3,-2/3;2/3,1/3,-2/3][2,3;-2,-6;1,0]

=[3,6;0,3;0,0]=R

Multiplying both sides on the left by H1^−1* H2 ^-1= H1*H2 yields the QR factorization:

[2,3;-2,-6;1,0]=H1*H2*R

=[2/3,-2/3,1/3;-2/3,-1/3,2/3;1/3,2/3,2/3][1,0,0;0,0,-1;0,-1,0][3,6;0,3;0,0]

=[2/3 -1/3 3/29][3 6]=QR

[-2/3 -2/3 2/29][0 3]

[1/3 -2/3 -4/29][0 0]

(b)We need to find a Householder reflector that moves the first column x=[-4;-2;4] to the vector w=[||x||2,0,0].Set v=w-x=[6;0;0]-[-4;-2;4]=[10;2;-4].Referring to Theorem 4.4,we have

H1=[1 0 0;0 1 0;0 0 1]-(2/120)*[100 20 -40;20 4 -8;-40 -8 16]=[ 1 0 0;0 1 0;0 0 1]-[5/3 1/3 -2/3;1/3 1/15 -2/15;-2/3 -2/15 4/15]=[-2/3 -1/3 2/3;-1/3 14/15 2/15;2/3 2/15 11/15]

And H1A=[-2/3 -1/3 2/3;-1/3 14/15 2/15;2/3 2/15 11/15][ -4,-4;-2,7;4,-5]=[6 -3; 0 36/15;0 -81/15]

The remaining step is to move the vector x^=[36/5;-27/5] to w^=[ 9;0]. Calculating H2^ from Theorem4.4 yields

V=w-x=[9;0]-[36/5;-27/5] =[9/5;27/5]

H2=I-2*v*v/v*v=[1 0;0 1]-(2*25/810)*[81/25 243/25;243/25 729/25]=[4/5 -3/5;-3/5 -4/5]

[4/5 -3/5;-3/5 -4/5][36/5;-27/5]=[9;0]

H2H1A=[1 0 0; 0 4/5 -3/5; 0 -3/5 -4/5][ -2/3 -1/3 2/3;-1/3 14/15 2/15;2/3 2/15 11/15][ -4,-4;-2,7;4,-5]=[6 -3;0 9;0 0]=R

Multiplying both sides on the left by H1^(-1)H2^(-1)=H1*H2 yields the QR factorization:

[-4,-4;-2,7;4,-5]=H1H2R=[-2/3 -1/3 2/3;-1/3 14/15 2/15;2/3 2/15 11/15][ 1 0 0; 0 4/5 -3/5; 0 -3/5 -4/5 ][ 6 -3;0 9;0 0]=[-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

13. Problem:Let P be the matrix defined in (4.29). Show (a) P^2 = P (b) P is symmetric (c) P v = v.

Solution:

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

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

(c) P*v=(v*v/v*v)*v=v

4.3Computer Problems

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

Code:function [Q,R]=CGS(A,n)

Q=A;

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

Result:

(a)

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

>> n=2;

>> [Q,R]=CGS(A,n)

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];

>> n=2

n =

2

>> [Q,R]=CGS(A,n)

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];

>> n=2;

>> [Q,R]=CGS(A,n)

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];

>> n=3;

>> [Q,R]=CGS(A,n)

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.0000   -0.6000

0   -1.0000   -0.0000

-0.6000   -0.0000    0.8000

R =

-5.0000  -10.0000   -5.0000

0   -2.0000    2.0000

0         0    5.0000

5.Problem: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]

Code:

function [Q,R,x]=CGS(A,n,b)

Q=A;

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

R1=R(1:n,:);

d=Q'*b;

d1=d(1:n,:);

x=R1d1;

Result:>> [Q,R,x]=CGS(A,n,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 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]

X=x0+xadd=x0+Qkck=[0;0;0]+[q1 q2 q3][ck1;ck2;ck3]= [0;0;0]+[q1 q2 q3 q4][h11 h12 h13;h21 h22 h23; 0 h32 h33;0 0 h43]

H11=

Ax=b

xk=x0+Qk*c=

||r||2=1

AQk=Qk+1*HK

[h11 h12 h13;h12 h22 h23;0 h32 h24;0 0 h25][c1 c2 c3]=[||r||2;0;0]=[1;0;0]

H11=q1’Aq1

3, 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 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.

A=sparse(1000,1000)

A =

全零稀疏矩阵: 1000×1000

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

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.

code:

A=sparse(1000,1000);

>>   for i=1:500

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

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

end

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