NLA-5.4

2377 days 전, jhlee2chn 작성

Algorithm 5.4.1  (Preconditioned 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 # Preconditioned Conjugate Gradient x = 10*random_vector(RDF, n) tol = 1e-08 r = b-dot(A, x);k = 0;bn = linalg.norm(b) k_max = 30 M = diag(10.0*ones(n)) while sqrt(dot(r,r)) > tol*bn and k < k_max: z = linalg.solve(M, r) k = k + 1 rho = dot(r, z) if k == 1: p = z else: bet=rho/rho_old p = z + bet*p w = dot(A,p) alp = rho/dot(p, w) x = x + alp*p;r = r-alp*w;rho_old=rho 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.00000001  0.99999998  1.          1.          1.00000001  1.     
1.
  1.00000001  1.          1.          1.00000001  1.          1.
  0.99999998  1.          1.          1.          1.          1.00000001
  0.99999998  1.00000001  1.          1.          1.          1.
  1.00000001  1.          1.00000001  1.00000001  1.          1.        
1.
  0.99999999  1.          1.00000001  1.          1.00000001  0.99999999
  1.00000001  0.99999999  1.          1.00000001  1.          1.00000001
  1.          1.00000002  1.          1.          1.00000001  1.
  1.00000001  1.          1.          1.          1.          0.99999999
  1.00000001  0.99999999  1.          1.00000001  1.          0.99999999
  1.          0.99999999  1.          1.          1.          0.99999999
  1.00000001  1.00000001  0.99999999  1.          1.          1.
  1.00000001  1.00000001  1.          1.          1.          1.        
1.
  0.99999999  1.          1.          1.          1.          1.00000001
  1.00000002  0.99999999  1.00000001  0.99999998  1.          1.        
1.
  1.00000001  0.99999998  0.99999999  0.99999999  1.         
0.99999999]
||r||= 1.91648814166e-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.00000001  0.99999998  1.          1.          1.00000001  1.          1.
  1.00000001  1.          1.          1.00000001  1.          1.
  0.99999998  1.          1.          1.          1.          1.00000001
  0.99999998  1.00000001  1.          1.          1.          1.
  1.00000001  1.          1.00000001  1.00000001  1.          1.          1.
  0.99999999  1.          1.00000001  1.          1.00000001  0.99999999
  1.00000001  0.99999999  1.          1.00000001  1.          1.00000001
  1.          1.00000002  1.          1.          1.00000001  1.
  1.00000001  1.          1.          1.          1.          0.99999999
  1.00000001  0.99999999  1.          1.00000001  1.          0.99999999
  1.          0.99999999  1.          1.          1.          0.99999999
  1.00000001  1.00000001  0.99999999  1.          1.          1.
  1.00000001  1.00000001  1.          1.          1.          1.          1.
  0.99999999  1.          1.          1.          1.          1.00000001
  1.00000002  0.99999999  1.00000001  0.99999998  1.          1.          1.
  1.00000001  0.99999998  0.99999999  0.99999999  1.          0.99999999]
||r||= 1.91648814166e-07
9