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

### 17074188 赵文轩 作业8

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

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

y1=[-4;-2;4]  ||y1||2=(16+4+16)^1/2=6  q1=y1/||y1||2=[-2/3;-1/3;2/3]

y2=A2-q1*q1'*A2=[-4;7;-5]-[-2/3;-1/3;2/3]*(-3)=[-6;6;-3] ||y2||2=(36+36+9)^1/2=9 q2=y2/||y2||2=[-2/3;2/3;-1/3]

r11=||y1||2=6  r22=||y2||2=9   r12=q1'*A2=-3

A=[-2/3 -2/3;-1/3 2/3;2/3 -1/3]*[6 -3;0 9]=Q*R

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

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

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

q1=[-4 ;  -2 ;  4]/6= [-2/3 ; -1/3 ; 2/3]

y2=A2=[-4;7;-5]

r12= q1^T*y2=[-2/3   -1/3   2/3]* [-4;7;-5]=-3

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

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

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

y3=A3=[1;0;0]

r13= q1^T*y3=[-2/3   -1/3   2/3] [1;0;0]=-2/3

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

r23= q2^T*y3=[-2/3  2/3  -1/3] [5/9;-2/9;4/9]=-2/3

y3=y3- r23*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/2)=1/3

q3=y3/ r33=[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 reflectors to find the full QR factorization of the matrices in Exercise2.

Sol:

(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 =

0.8333    0.1667   -0.3333

0.1667    0.0333   -0.0667

-0.3333   -0.0667    0.1333

H1=I-2P=

-0.6667   -0.3333    0.6667

-0.3333    0.9333    0.1333

0.6667    0.1333    0.7333

H1*A=

6.0000   -3.0000

0   7.2000

0   -5.4000

x2=[7.2,-5.4]

w2=[||x||2;0]=[9,0]

v2=w2-x2=[1.8,5.4]

P2=v2*v2’/(v2’*v2)

=

0.1    0.3

0.3    0.9

H2=I-2*P2

=

0.8000   -0.6000

-0.6000   -0.8000

R=H2*H1*A

=

6.0000   -3.0000

-0.0000    9.0000

0.0000   -9.0000

Q=H2*H1

=

-0.6667   -0.8000    0.6000

-0.3333    0.6400   -0.4800

0.6667   -0.4800    0.3600

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

(a)Sol:

From the define of the P,we know P=v*v’/(v’*v)

Then P^2=[v*v’/(v’*v)]^2=v*v’*v*v’/(v’*v*v’*v)

It is easy to find that k=v’*v,k is a real number.

P^2=k*v*v’/(k*v’*v)=v*v’/(v’*v)=P

(b)Sol:

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

So P is symmetric

(c)Sol:

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

4.3 Computer Problems

1Write a MATLAB program that implements classical Gram-Schmidt to find the reduced QR factorization.Check

your work by comparing factorizations  of 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,m,n)

q=zeros(size(A));

r=zeros(n);

y=zeros(m,1);

for j=1:1:n

y=A(:,j);

for i=1: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

Q=q;

R=r;

Results:

Exercise 1(a)

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

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

Q =

0.8000   -0.6000

0.6000    0.8000

R =

5.0000    0.6000

0    0.8000

Check:

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

Q =

-0.8000   -0.6000

-0.6000    0.8000

R =

-5.0000   -0.6000

0    0.8000

1(b)

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

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

Q =

0.7071    0.7071

0.7071   -0.7071

R =

1.4142    2.1213

0    0.7071

Check:

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

Q =

-0.7071   -0.7071

-0.7071    0.7071

R =

-1.4142   -2.1213

0   -0.7071

5Use the MATLAB QR factorization to find the least squares solution and 2-norm error of the following inconsistent system;

(a)[1 1;2 1;1 2;0 3]*[x1;x2]=[3;5;5;5;5]

(a)

Code:

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

>> y=[ 3; 5; 5; 5];

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

>> b=Q'*y;

>> c=R(1:2,1:2)b(1:2)

c =

1.6154

1.6615

>>  e=norm(b(3:4),2)

e =

0.3038

P229 4.4

Solve Ax=b for the following A and b=[1,0,0]’ ,usingGMRES with x0=[0,0,0]’.Report all approximation xk up to and including the correct solution.

(c)[0 0 1; 0 0;0 1 0]

Solution

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

>> b=[1;0;0];

>> X=gmres(A,b)

gmres 在解的 迭代 3 处收敛，并且相对残差为 0

X =

0

0

1

Computer 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 fit within the matrix.

(a) Print the nonzero structure spy(A).

(b) Let xe be the vector of n ones.Set b=A*xe,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):

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

Results:

>> A=P(1000);

>> spy(A)

Problem 2:Let n=1000.Start with the n*n matrix A from Computer Problem 1,and add the nonzero entries A(i,2*i)=A(2*i,i)=1/2 for

1<=i<=n/2.Carry out steps (a) and (b) as in that problem.

(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

Results:

>> A=P(1000);

>> spy(A)