NLA-5.3

2027 days 전, jhlee2chn 작성

Algorithm 5.3.1  (A Practical Form of CG-Method for solving LSE with p.d. matrix)

from numpy import * n = 100 A = diag(10.0*ones(n))+diag(ones(n-1), k=1)+diag(ones(n-1), k=-1) # tridiagonal matrix 생성 b = 12.0*ones(n) b[0] = 11.0;b[n-1] = 11.0 print "A=" print A print "b=" print b # Practical Conjugate Gradient x = 10*random_vector(RDF, n) tol = 1e-08 r = b-dot(A, x);k = 0;bn = linalg.norm(b) rho = dot(r,r) k_max = 30 while sqrt(rho) > tol*bn and k < k_max: k = k + 1 if k == 1: p = r else: bet=rho/rho_old p = r + bet*p w = dot(A,p) alp = rho/dot(p, w) x = x + alp*p;r = r-alp*w;rho_old=rho;rho=dot(r,r) print "x=", x print "||r||=", linalg.norm(r, inf) print k 
       
A=
[[ 10.   1.   0. ...,   0.   0.   0.]
 [  1.  10.   1. ...,   0.   0.   0.]
 [  0.   1.  10. ...,   0.   0.   0.]
 ..., 
 [  0.   0.   0. ...,  10.   1.   0.]
 [  0.   0.   0. ...,   1.  10.   1.]
 [  0.   0.   0. ...,   0.   1.  10.]]
b=
[ 11.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12. 
12.
  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12. 
12.
  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12. 
12.
  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12. 
12.
  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12. 
12.
  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12. 
12.
  12.  12.  12.  12.  12.  12.  12.  12.  12.  11.]
x= [ 1.          1.          1.          0.99999999  1.          1.
  0.99999999  0.99999999  1.          1.          1.          1.
  0.99999999  1.00000002  1.          0.99999999  1.          0.99999999
  1.          1.00000001  1.          0.99999999  1.          1.        
1.
  1.          1.00000001  1.          1.          1.          1.
  1.00000001  1.00000002  1.          1.          0.99999998  0.99999999
  1.00000002  1.          1.          0.99999999  1.          1.00000001
  1.          1.00000002  0.99999999  1.          1.          1.00000002
  0.99999999  1.00000002  1.00000001  1.00000001  0.99999999  0.99999998
  1.          1.00000001  1.          1.00000001  1.          1.00000001
  1.          1.00000001  1.00000001  0.99999999  1.          1.00000002
  0.99999999  1.00000001  1.00000001  1.          1.00000001  0.99999998
  0.99999999  1.00000001  0.99999999  1.          0.99999999  1.00000002
  0.99999999  0.99999999  1.00000003  1.          1.          1.00000001
  1.          1.00000001  1.00000002  0.99999999  1.00000002  0.99999999
  1.          1.          0.99999998  1.          0.99999998  1.        
1.
  1.          1.        ]
||r||= 2.60275350004e-07
9
A=
[[ 10.   1.   0. ...,   0.   0.   0.]
 [  1.  10.   1. ...,   0.   0.   0.]
 [  0.   1.  10. ...,   0.   0.   0.]
 ..., 
 [  0.   0.   0. ...,  10.   1.   0.]
 [  0.   0.   0. ...,   1.  10.   1.]
 [  0.   0.   0. ...,   0.   1.  10.]]
b=
[ 11.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.
  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.
  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.
  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.
  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.
  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.  12.
  12.  12.  12.  12.  12.  12.  12.  12.  12.  11.]
x= [ 1.          1.          1.          0.99999999  1.          1.
  0.99999999  0.99999999  1.          1.          1.          1.
  0.99999999  1.00000002  1.          0.99999999  1.          0.99999999
  1.          1.00000001  1.          0.99999999  1.          1.          1.
  1.          1.00000001  1.          1.          1.          1.
  1.00000001  1.00000002  1.          1.          0.99999998  0.99999999
  1.00000002  1.          1.          0.99999999  1.          1.00000001
  1.          1.00000002  0.99999999  1.          1.          1.00000002
  0.99999999  1.00000002  1.00000001  1.00000001  0.99999999  0.99999998
  1.          1.00000001  1.          1.00000001  1.          1.00000001
  1.          1.00000001  1.00000001  0.99999999  1.          1.00000002
  0.99999999  1.00000001  1.00000001  1.          1.00000001  0.99999998
  0.99999999  1.00000001  0.99999999  1.          0.99999999  1.00000002
  0.99999999  0.99999999  1.00000003  1.          1.          1.00000001
  1.          1.00000001  1.00000002  0.99999999  1.00000002  0.99999999
  1.          1.          0.99999998  1.          0.99999998  1.          1.
  1.          1.        ]
||r||= 2.60275350004e-07
9