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

### 数值分析作业八

Chap 4.3

Exercises

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

following matrices:

Solution:

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

y1=A1=[-4;-2;4]   A2=[-4;7;-5]

r11=||y1||2=(16+4+16)^(1/2)=6

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

r12=q1’*A2=-3

y2=A2-q1(q1’*A2)=[-6;6;-3]

r22=y2/||y2||2=9

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

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

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

So,

Q=[-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 modified Gram–Schmidt orthogonalization to find the full QR factorization of the matrices in Exercise 2.

Solution:

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

y1=A1=[2;-2;1]   A2=[3;-6;0]

r11=||y1||2=(4+4+1)^(1/2)=3

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

y2=A2-q1(q1’*A2)=[-1;-2;-1]

r12=q1’*A2=6

r22=y2/||y2||2=3

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

y3=A3-q1(q1’*A3)-q2(q2’*A3)=[4/9;2/9;-4/9]

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

So,

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

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

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

y1=A1=[-4;-2;4]   A2=[-4;7;-5]

r11=||y1||2=(16+4+16)^(1/2)=6

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

r12=q1’*A2=-3

y2=A2-q1(q1’*A2)=[-6;6;-3]

r22=y2/||y2||2=9

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

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

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

So,

Q=[-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]

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

Solution:

(a)    x=[2,-2,1],so w=[||x||2,0,0]=[3,0,0],let v=w-x=[1,2,-1],

so P=v’v/vv’=[1/6,1/3,-1/3;1/3,2/3,-1/3;-1/6,-1/3,1/6]

H1=I-2P=[2/3,-2/3,1/3;-2/3,-1/3,2/3;1/3,2/3,2/3]

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

x1=[0,-3],so w1=[||x1||2]=[3,0],let v1=w1-x1=[3,3]

so P1= v1’v1/ v1v1’=[1/2,1/2;1/2,1/2]

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

Then,H2H1A=[3,6;0,3;0,0]

So,

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

(b)   x=[-4,-2,4],so w=[||x||2,0,0]=[6,0,0],let v=w-x=[10,2,-4],

so P=v’v/vv’=[1/6,1/6,-1/3;1/6,1/6,-1/3;-1/3,-1/3,2/3]

H1=I-2P= [-2/3,-1/3,2/3;-1/3,14/15,2/15;2/3,2/15,11/15]

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

x1=[ 36/5;-27/5],so w1=[||x1||2]=[9,0],let v1=w1-x1=[ 9/5;27/5]

so P1= v1’v1/ v1v1’=[4/5,-3/5;-3/5,-4/5]

H2=[ 1,0,0;0,4/5,-3/5;0,-3/5,-4/5]

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

So,

Q=[-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]

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

Proof

(a)    v=w-x,norm(w,2)=norm(x,2)

P=v*v’/v’*v,P2= (v*v’/v’*v)*( v*v’/v’*v)= v*v’/v’*v=P

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

So,P=P’;

Then,we can prove P is symmetric.

(c)    Pv=( v*v’/v’*v)*v=[v*(v’*v)]/ v’*v=v

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.

Code:

GS.m

function [q] = GS(A)

q(:,1)=A(:,1)/norm(A(:,1)) %提取矩阵的第一列

[Ahang,Alie]=size(A)  %矩阵的行和列

for j=1:Alie

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

MGS.m

function [q] = MGS(A)

q(:,1)=A(:,1)/norm(A(:,1)) %提取矩阵的第一列

[Ahang,Alie]=size(A)  %矩阵的行和列

for j=1:Alie

y=A(:,j); %临时变量 存储后面的求和

for i=1:j-1

r(i,j)=q(:,i)'*y;

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

end

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

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

end

5. Use the Matlab QR factorization to find the least squares solutions and 2-norm error of the following inconsistent systems:

Solution:

Code:

function [q] = MGS(A)

q(:,1)=A(:,1)/norm(A(:,1)) %提取矩阵的第一列

[Ahang,Alie]=size(A)  %矩阵的行和列

for j=1:Alie

y=A(:,j); %临时变量 存储后面的求和

for i=1:j-1

r(i,j)=q(:,i)'*y;

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

end

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

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

end

r =

2.4495    2.0412

0    3.2914

q =

0.4082    0.0506

0.8165   -0.2025

0.4082    0.3545

0    0.9115

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

Solution:

[h11,h12,h13;h21,h22,h23;0,h32,h33;0,0,0,h43]*[c1;c2;c3]=[1;0;0;0]

h11=1,

x1=[0.0332,0.8505,0.9668]

x2=[0.0672,0.8779,0.9696]

x3=[0,0,1]

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 ﬁ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:

(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

>> A=P(1000);

>> spy(A)

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

[n,n]=size(A);

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

Solution

(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