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

# 作业八

### EXERCISE 4.3(QR Factorization)

2. Apply classical Gram-Schmidt orthogonalization to find the full QR factorization of the matrix A=[-4 -4;-2 7;4 -5];

Solution

Denote A1=[-4;-2;4] and A2=[-4;7;-5];adding the third vector A3=[3;0;0];

Then we have y1=A1 and q1=y1/norm(y1,2)=[-2/3;-1/3;2/3];

r(1,1)=norm(y1,2)=6,r12=q1’*A2=-3;

then we can get that y2=A2-q1*r12=[-6;6;3]and q2=[-2/3;2/3;-1/3]; and r22=norm(y2,2)=9;

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

Above all we can get that the QR factorization of A is that A=Q*R=[q1,q2,q3]*R=[-2/3,-2/3,1/3;-1/3,2/3,2/3;2/3,-1/3,2/3]*[6,-3;0,9;0,0];

4. Apply the modified Gram-Schmidt orthogonalization to find the full QR factorization of the matrices in EX2.(A=[-4 -4;-2 7;4 -5])

Solution:

Denote A1=[-4;-2;4] and A2=[-4;7;-5];adding the third vector A3=[3;0;0];

Then we have y1=A1 and q1=y1/norm(y1,2)=[-2/3;-1/3;2/3];

r(1,1)=norm(y1,2)=6,r12=q1’*A2=-3;

then we can get that y2=A2-q1*r12=[-6;6;3]and q2=[-2/3;2/3;-1/3]; and r22=norm(y2,2)=9;

y31=A3-q1*(q1’*A3)=[5/3;-2/3;4/3] and then y3=y31-q2*(q2’*y31)=[1/3;0;0]

Above all we can get that the QR factorization of A is that A=Q*R=[q1,q2,q3]*R=[-2/3,-1/3,2/3;-1/3,2/3,0;2/3,-1/3,0]* [6,-3;0,9;0,0];

6.Apply Householder reflectors to find the full QR factorization of the matrix in EX2.( A=[-4 -4;-2 7;4 -5])

Solution:

Denote x1 the first column of A that x1=[-4;-2;4], and then we can let w1=[norm(x,2);0;0]=[6;0;0].

Let v1=w1-x1=[10;2;-4]. We have that H1=[1,0,0;0,1,0;0,0,1]-2(v*v’)/v’*v

=[-2/3,-1/3,2/3;-1/3,14/15,2/15;2/3,2/15,11/15].

Then we have H1*A=[6,-3;0,36/5;0,-27/5]. Denote x2=[36/5;-27/5]; let w2=[norm(x2,x);0;0] =[9;0]. Then, v2=w2-x2=[9/5;27/5]. Note that h2=H2(2:3,2:3), where H2(1,1)=1, H2(2:3,1)=zeros(2,1) and H2(1,2:3)=zeros(1,2).

h2=[1,0;0,1]-2*(v2*v2’)/(v2’*v2)=[4/5,-3/5;-3/5,-4/5]. That is H2=[1,0,0;0,4/5,-3/5;0,-3/5,-4/5].

Then H2*H1*A=[6,-3;0,9;0,0]; and Q=H1*H2=[-2/3,-2/3,-1/3;-2/3,2/3,-2/3;2/3,-1/3,-2/3], R=[6,-3;0,9;0,0]. That is A=Q*R=[-2/3,-2/3,-1/3;-2/3,2/3,-2/3;2/3,-1/3,-2/3]* [6,-3;0,9;0,0].

13. Let P be the matrix that P=v*v’/(v’*v) where v is a vector of size n*1. Show (a) P^2=P;(b) P is a symmetric; (c)P*v=v.

Proof: (a):we can know that v’*v is a real number;

P*P=[(v*v’)/(v’*v)]* [(v*v’)/(v’*v)]=(v*v’*v*v’)/(v’*v*v’*v)= (v*(v’*v)*v’)/(v’*v*v’*v)=v*v’/(v’*v)

=P; that is P^2=P.

(b) P’=[(v*v’)/(v’*v)]’=(v*v’)’/(v’*v)=(v*v’)/(v’*v)=P, so P is symmetric.

(c) P*v=[(v*v’)/(v’*v)]*v=(v*v’*v)/(v’*v)= v*(v’*v)/(v’*v)=v; that is P*v=v.

### COMPUTER PROBLEMS 4.2

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

Solution

The code of classical QR factorization:

%classical Gram-Schmidt to the reduced QR factorization

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

The checking code:

%Exercise 4.3.1

A1=[4 0;3 1];A2=[1 2;1 1];A3=[2 1;1 -1;2 1];A4=[4 8 1;0 2 -2;3 6 7];

[q11,r11]=qr(A1,0);[q12,r12]=reducedQR(A1);

[q21,r21]=qr(A2,0);[q22,r22]=reducedQR(A2);

[q31,r31]=qr(A3,0);[q32,r32]=reducedQR(A3);

[q41,r41]=qr(A4,0);[q42,r42]=reducedQR(A4);

The result of the checking is:

q11 =[-0.8000,-0.6000;-0.6000,0.8000]   r11 =[-5.0000,-0.6000;0,0.8000]

q12 =[0.8000,-0.6000;0.6000,0.8000]    r12 =[5.0000,0.6000;0,0.8000]

q21 =[-0.7071,-0.7071; -0.7071,0.7071]   r21 =[-1.4142,-2.1213;0,-0.7071]

q22 =[0.7071,0.7071;0.7071,-0.7071]    r22 =[1.4142,2.1213;0,0.7071]

q31 =[-0.6667,0.2357; -0.3333,-0.9428; -0.6667,0.2357]   r31 =[-3.0000,-1.0000;0,1.4142]

q32 =[0.6667,0.2357;0.3333,-0.9428; 0.6667, 0.2357]     r32 =[3.0000,1.0000;0,1.4142;]

q41 =[-0.8,0,-0.6;0,-1, 0; -0.6,0,0.8]  r41 =[-5,-10,-5; 0,-2,2;0,0,5]

q42 =[0.8,0,-0.6;0,1,0;0.6,0,0.8]   r42 =[5,10,5; 0,2,-2;0,0,5]

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

Solution:

The code of the problem is:

A1=[1 1;2 1;1 2;0 3];b1=[3;5;5;5];

[Q,R]=qr(A1);k=2;

n=size(A1,1);

b=inv(Q)*b1;

for i=k:-1:1

for j=i+1:k

b(i)=b(i)-x(j)*R(i,j);

end

x(i)=b(i)/R(i,i);

end

The least squares solution of the equation is: x =[1.6154,1.6615]

And the 2-norm error is norm((A1*x-b1),2)= 0.3038.

### EXERCISE4.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.(A=[0 0 1;1 0 0;0 1 0])

Solution:

First, r=b-A*x0=[1;0;0] and q1=r/norm(r,2)=[1;0;0]; then we have A*q1=[0;1;0] and h11=q1’ *A*q1=0; y1=A*q1-h11*q1=[0;1;0], h21=norm(y1)=1 and q2=[0;1;0]. x1=[0;0;0].

Then A*q2=[0;0;1], h12=q1’*A*q2=0 and y21=A*q2-0*q1=[0;0;1]; h22=q2’*y21=0, y2=y21- h22*q2=[0;0;1],h32=1 and q3=[0;0;1]. x2=[0;1;0].

Then A*q3=[1;0;0],h13=1, y31=A*q3-h13*q1=[0,0,0]’;h23=0,y3=[0;0;0],and h33=0. x3=[0;0;1]

The steps stop and the correct solution is 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.

Solution:

### COMPUTER PROBLEMS 4.4

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:

The code of the question is:

%CP4.4.1

n=1000;

A=diag(1:n)+diag(1:(n-1),1)+diag(1/2*ones(1,n-1),-1)+diag(1/2*ones(1,(n-2)),2)+diag(ones(1,n-2),-2);//生成矩阵

spy(A);//显示矩阵分布

xe=ones(n,1);x0=zeros(n,1);

b=A*xe;

x1=CG(A,b,x0,100);

x2=Jacobi(A,b,x0,100);

x3=Gauss_Seidel(A,b,x0,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 1i n/2. Carry out steps (a) and (b) as in that problem.

Solution:

%CO4-4-2

n=1000;

A=diag(1:n)+diag(1:(n-1),1)+diag(1/2*ones(1,n-1),-1)+diag(1/2*ones(1,(n-2)),2)+diag(ones(1,n-2),-2);

for i=1:(n/2)

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

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

end

spy(A);

xe=ones(n,1);x0=zeros(n,1);

b=A*xe;

x1=CG(A,b,x0,100);

x2=Jacobi(A,b,x0,100);

x3=Gauss_Seidel(A,b,x0,100);