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

17074186作业5

6. Find the Cholesky factorization A=RT R of each matrix.

(a)  [ 4 2 0 ,2 2 3, 0 3 10]

`function [R,Rt]=Cholesky(A)`

n=size(A);

for k=1:n

if A(k,k)<0;

end

R(k,k)=sqrt(A(k,k));

u=A(k+1:n,k)/R(k,k);

ut=A(k,k+1:n)/R(k,k);

R(k,k+1:n)=ut;

A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-u*ut;

end

Rt=R';

>> A=[4 -2 0;-2 2 -3;0 -3 10]

A =

4    -2     0

-2     2    -3

0    -3    10

>> [R,Rt]=Cholesky(A)

R =

2    -1     0

0     1    -3

0     0     1

Rt =

2     0     0

-1     1     0

0    -3     1

12. Prove that a principal submatrix of a symmetric positive-denite matrix is symmetric positive-denite. (Hint: Consider an appropriate X and use Property 2.)

13. Solve the problems by carrying out the Conjugate Gradient Method by hand. (a) [1 2 ,2 5 ] [u, v ]= [1, 1 ]

Solution:x0=[0;0];

r0=d0=[1;1];

a0=([1;1]'*[1;1])/([1;1]'*[1 2;2 5]*[1;1])=1/5;

x1= [1/5;1/5];

r1=[2/5;-2/5];

β0=(r1'*r1)/(r0'*r0)= 4/25;

d1=r1+β0*d0= [14/25;-6/25];

a1=5

x2= [3;-1];

r2=[0;0]

So r2=b-A*x2=[0;0], x2=[3;-1]

14. Solve the problems by carrying out the Conjugate Gradient Method by hand. (a) [1 1, 1 2 ] [u v] = [0 1]

solution:x0=[0;0];

r0= [0;1];

α0=([0;1]'*[0;1])/([0;1]'*[1 -1;-1 2]*[0;1])=1/2;

x1= [0;0]+1/2*[0;1]=[0;1/2];

r1=[0;1]-1/2*[1 -1;-1 2]*[0;1]=[1/2;0];

β0= ([1/2;0]'*[1/2;0])/1=1/4;

d1= [1/2;0]+1/4*[0;1]=[1/2;1/4];

α1=(1/4)/([1/2;1/4]'*[1 -1;-1 2]*[1/2;1/4])=2;

x2 =[0;1/2]+2*[1/2;1/4]=[1;1];

r2 =[1/2;0]-2*[1 -1;-1 2]*[1/2;1/4]=[0;0]

Since r2=b-A*x2=[0;0], x2=[1;1]

2. Use a Matlab version of conjugate gradient to solve the following problems: (a) [ 1 1 0 ,1 2 1 ,0 1 2] [ u, v, w] = [ 0 ,2, 3 ]  

(b)[ 1 1 0 ,1 2 1, 0 1 5 ][ u v w ] =[3 ,3 ,4]

function xc=CG(x0,A,b,n)

x1=x0;

r1=b-A*x1;

d1=r1;

for k=1:n

if r1==0

xc=x1;

end

a1=(r1'*r1)/(d1'*A*d1);

x2=x1+a1*d1;

r2=r1-a1*A*d1;

b1=(r2'*r2)/(r1'*r1);

d2=r2+b1*d1;

r1=r2;

x1=x2;

d1=d2;

end

xc=x1;

>> x0=[0;0;0];

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

b=[0;2;3];

xc=CG(x0,A,b,5)

xc =

1.0000

1.0000

1.0000

 

>> x0=[0;0;0];

A=[1 -1 0;-1 2 1;0 1 5];

b=[3;-3;4]

xc=CG(x0,A,b,5)

b =

3

-3

4

xc =

2.0000

-1.0000

1.0000

4. Solve the sparse problem of (2.45) by the Conjugate Gradient Method for (a) n=6 (b) n=12.

function [A,b]=sparsesetup(n)

e=ones(n,1);n2=n/2;

A=spdiags([-e 3*e -e],-1:1,n,n);%entries of a

c=spdiags([e/2],0,n,n);c=fliplr(c);A=A+c;

A(n2+1,n2)=-1;A(n2,n2+1)=-1;%fix up 2 entries

b=zeros(n,1);

b(1)=2.5;b(n)=2.5;b(2:n-1)=1.5;b(n2:n2+1)=1;

function xc=CG(x0,A,b,n)

x1=x0;

r1=b-A*x1;

d1=r1;

for k=1:n

if r1==0

xc=x1;

end

a1=(r1'*r1)/(d1'*A*d1);

x2=x1+a1*d1;

r2=r1-a1*A*d1;

b1=(r2'*r2)/(r1'*r1);

d2=r2+b1*d1;

r1=r2;

x1=x2;

d1=d2;

end

xc=x1;

>> n=6

n =

6

>> [A,b]=sparsesetup(n)

A =

(1,1)       3.0000

(2,1)      -1.0000

(6,1)       0.5000

(1,2)      -1.0000

(2,2)       3.0000

(3,2)      -1.0000

(5,2)       0.5000

(2,3)      -1.0000

(3,3)       3.0000

(4,3)      -1.0000

(3,4)      -1.0000

(4,4)       3.0000

(5,4)      -1.0000

(2,5)       0.5000

(4,5)      -1.0000

(5,5)       3.0000

(6,5)      -1.0000

(1,6)       0.5000

(5,6)      -1.0000

(6,6)       3.0000

b =

2.5000

1.5000

1.0000

1.0000

1.5000

2.5000

>> x0=[0;0;0;0;0;0]

x0 =

0

0

0

0

0

0

>> xc=CG(x0,A,b,n)

xc =

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

>> n=12

n =

12

>> [A,b]=sparsesetup(n)

A =

(1,1)       3.0000

(2,1)      -1.0000

(12,1)       0.5000

(1,2)      -1.0000

(2,2)       3.0000

(3,2)      -1.0000

(11,2)       0.5000

(2,3)      -1.0000

(3,3)       3.0000

(4,3)      -1.0000

(10,3)       0.5000

(3,4)      -1.0000

(4,4)       3.0000

(5,4)      -1.0000

(9,4)       0.5000

(4,5)      -1.0000

(5,5)       3.0000

(6,5)      -1.0000

(8,5)       0.5000

(5,6)      -1.0000

(6,6)       3.0000

(7,6)      -1.0000

(6,7)      -1.0000

(7,7)       3.0000

(8,7)      -1.0000

(5,8)       0.5000

(7,8)      -1.0000

(8,8)       3.0000

(9,8)      -1.0000

(4,9)       0.5000

(8,9)      -1.0000

(9,9)       3.0000

(10,9)      -1.0000

(3,10)      0.5000

(9,10)     -1.0000

(10,10)      3.0000

(11,10)     -1.0000

(2,11)      0.5000

(10,11)     -1.0000

(11,11)     3.0000

(12,11)     -1.0000

(1,12)      0.5000

(11,12)     -1.0000

(12,12)      3.0000

b =

2.5000

1.5000

1.5000

1.5000

1.5000

1.0000

1.0000

1.5000

1.5000

1.5000

1.5000

2.5000

>> x0=[0;0;0;0;0;0;0;0;0;0;0;0]

x0 =

0

0

0

0

0

0

0

0

0

0

0

0

>> xc=CG(x0,A,b,n)

xc =

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

Copyright ©2020    计算数学达人 All Right Reserved.
 技术支持：自助建站 | 领地网站建设 |短信接口 |燕窝 版权所有 © 2005-2020 lingw.net.粤ICP备16125321号 -5