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

### 17074216作业8

4.3EXERCISE

2. Apply classical Gram–Schmidt orthogonalization to fifind the full QR factorization of the

following matrices:

(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

4.Apply modifified Gram–Schmidt orthogonalization to fifind the full QR factorization of the

matrices in Exercise 2

y1=A1=[ -4 ; 2 ; 4] ,  r11=||y1||2=(16+4+16)^1/2=6

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

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

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

r12=q1A2=-3 , r22=||y2||2=9

Adding athird A3=[1 ; 0 ; 0] leads to

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

y3=y3(1)- [-2/3 ; 2/3 ; -1/3](-2/3)=[1/9 ; 2/9 ; 2/9]

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

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 reflflectors to fifind the full QR factorization of the matrices in Exercise 2

x=[2;-2;1], w=[||x||2;0;0]=[3;0;0], v=w-x=[3;0;0]-[2;-2;1]=[1;2;-1],

p=(v*v')/(v'*v)=1/6*[1 2 -1;2 4 -2;-1 -2 1],

H1=[1 0 0;0 1 0;0 0 1]-2/6*[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],

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

x1=[0;-3], w1=[3;0], v1=w1-x1=[3;0]-[0;-3]=[3;3],

p1=(v1*v1')/(v1'*v1)=1/18*[9 9;9 9]=1/2*[1 1;1 1],

H2=I-2*p1=[1 0;0 1]-[1 1;1 1]=[0 -1;-1 0],

H2*x1=[0 -1;-1 0]*[0;-3]=[3;0]=w2,

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]=[3 6;0 3;0 0]=R,

[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 2/3;-2/3 -2/3 1/3;1/3 -2/3 -2/3]*[3 6;0 3;0 0]=Q*R.

(b):

x=[-4;-2;4], w=[||x||2;0;0]=[6;0;0], v=w-x=[6;0;0]-[-4;-2;4]=[10;2;-4],

p=(v*v')/(v'*v)=1/120*[100 20 -40;20 4 -8;-40 -8 16]=1/30*[25 5 -10;5 1 -2;-10 -2 4],

H1=[1 0 0;0 1 0;0 0 1]-2/30*[25 5 -10;5 1 -2;-10 -2 4]=[-2/3 -1/3 2/3;-1/3 14/15 2/15;2/3 2/15 11/15],

H1*A=[-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/5;0 -27/5],

x1=[36/5;-27/5], w1=[9;0], v1=w1-x1=[9;0]-[36/5;-27/5]=[9/5;27/5],

p1=(v1*v1')/(v1'*v1)=1/2*[1/5 3/5;3/5 9/5],

H2=I-2*p1=[1 0;0 1]-[1/5 3/5;3/5 9/5]=[4/5 -3/5;-3/5 -4/5],

H2*x1=[4/5 -3/5;-3/5 -4/5]*[36/5;-27/5]=[9;0]=w1,

H2*H1*A=[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,

[-4 -4;-2 7;4 -5]=H1*H2*R=[-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]=Q*R.

13. Let P be the matrix defifined in (4.29). Show (a) P

2 = P (b) P is symmetric (c) P v = v.

(a)
P^2=v*v’*v*v’/(v’*v*v’*v)
Due to the result of v’*v is a real number, hence, v*v’*v*v’= v’*v*v*v’
So, we can prove P^2=P.

(b)
(v*v’)’=(v’)’*v’=v*v’
S0, P is symmetric.

(c)
Due to the result of v’*v is a real number, hence, v*v’*v= v’*v*v
So, P*v=v.

function [Q R]=Householder(A)
[m,n]=size(A);H1=eye(m);B=A;
for i=1:n
H=eye(m);
x=A(i:m,i);
y=norm(x,2);
B=eye(m);B(i,i)=y;
w=B(i:m,i);
v=w-x;
P=(v*v')/(v'*v);
H(i:m,i:m)=eye(m+1-i)-2*P;
A=H*A;
H1=H1*H;

end
Q=H1;
R=A;

4.3computer problem
1. Write a Matlab program that implements classical Gram–Schmidt to fifind 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 QR=CGS(A)

[m n]=size(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

(a)

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

>> QR=CGS(A)

R =

5.0000    0.6000

0    0.8000

Q =

0.8000   -0.6000

0.6000    0.8000

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

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

>> QR=CGS(A)

R =

1.4142    2.1213

0    0.7071

Q =

0.7071    0.7071

0.7071   -0.7071

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

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

>> QR=CGS(A)

R =

3.0000    1.0000

0    1.4142

Q =

0.6667    0.2357

0.3333   -0.9428

0.6667    0.2357

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

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

>>  QR=CGS(A)

R =

5    10     5

0     2    -2

0     0     5

Q =

0.8000         0   -0.6000

0    1.0000         0

0.6000         0    0.8000

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

>> [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 fifind 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]

a）The least squares solutions is x=[1.6154 1.6615].

Code： A=[1 1;2 1;1 2;0 3];

>> y=[3 5 5 5];

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

>> c=Q'*y'

c =

-7.3485

5.4688

0.2158

0.2139

>> x=R(1:4,1:2)c(1:4)

x =

1.6154

1.6615

Result： n=norm(y'-A*x,2)

n =

0.3038

4.4computer 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 fifitwithin 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.

(a)

function matrixA(n)

S1=sparse(1:n,1:n,1:n);

E=sparse(2:n,1:n-1,(1/2)*ones(1,n-1),n,n);

F=sparse(3:n,1:n-2,(1/2)*ones(1,n-2),n,n);

A=S1+E+E'+F+F';

spy(A)

title('The nonzero structure spy(A)','color','b')

end

>> matrixA(1000)

(b)function matrix4(n)

S1=sparse(1:n,1:n,1:n);

E=sparse(2:n,1:n-1,(1/2)*ones(1,n-1),n,n);

F=sparse(3:n,1:n-2,(1/2)*ones(1,n-2),n,n);

A=S1+E+E'+F+F';

xe=ones(n,1);

b=A*xe;

[x0,f10,rr0,it0,rv0]=pcg(A,b,1e-8,100);

[itr rnorm]=cg_it(A,b,1e-8,100);

[j y]=jacobi(A,b,100);

semilogy(0:it0,rnorm,'bo-',0:it0,rv0/norm(b),'rs-',1:j,y,'kd--',...'LineWidth',1);

xlabel('Step number'));

ylabel('Relative residual');

legend('The Conjugate Gradient Method','The Jacobi preconditioner',...'Gauss-Seidel Method');

title('Efficiency of Preconditioned Conjugate Gradient Method'))

end

>> matrix4(100)

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.

The results showed that the GMRES without preconditions were slow to converge, and the use of the Jacobi precondition was improved greatly, and the GMRES with gauss-seidel preconditions were rapidly convergent.

Code:

(a)

function [A]=Spy(n)

a=zeros(n);

for i=1:n/2

t=2*i;

a1=find(a(i,t));

a1=1/2;

a2=find(a(t,i));

a2=1/2;

end

a1=(ones(1,n-1)*1/2);

a2=(1:n);

A=diag(a1,1)+diag(a1,-1)+diag(a2)+a;

spy(A)

test:

clear

clc

n=1000;

[A]=Spy(n)

(b)

function [xk]=GMRES(A,b)

r=(M.^(-1))*(b-A*x0);

q(:,1)=r/norm(r);

for k=1:m

w=(M.^(-1))*A*q(:,k);

for j=1:k

h(j,k)=w'*q(:,j);

w=w-h(j,k)*q(:,j);

end

h(k+1,k)=norm(w);

q(:,k+1)=w/h(k+1,k);

ck1=[norm(r) zeros(1,length(r)-1)];

ck=norm(h(:,1:ck)-ck1');

xk=q(:,1:ck)+x0;

end