现在时间是:
当前位置:首 页 >> 数值分析(英)>> 教学区>> 文章列表

数值分析作业八

作者:孙思姣   发布时间:2019-06-13 20:09:24   浏览次数:73

 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]

adding a third vector A3=[1;0;0]

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]

adding a third vector A3=[1;0;0]

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]

adding a third vector A3=[1;0;0]

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

Answer:

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

Answer

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

 







上一篇:没有了    下一篇:没有了

Copyright ©2020    计算数学达人 All Right Reserved.

技术支持:自助建站 | 领地网站建设 |短信接口 |燕窝 版权所有 © 2005-2020 lingw.net.粤ICP备16125321号 -5