Aが対称行列のときにAx=bを解く、反復法の最も基本的なやつです。
Aの状態数が悪い場合は、
不完全コレスキー分解等を用いて前処理を行う方法があります。
from numpy import *
import copy
# 共役勾配法によりAx=bを解く(Aは正値対称行列)
def ConjugateGradientMethod(amat, b, eps, nloop):
x = matrix(zeros((amat.shape[1],1)))
r = b - amat*x
p = copy.copy(r)
for k in range(nloop):
if linalg.norm(r) < eps:
break
alpha = (r.T*r)/(p.T*amat*p)
x += alpha[0,0]*p
rold = copy.copy(r)
r -= alpha[0,0]*amat*p
beta = (r.T*r)/(rold.T*rold)
p = r + beta[0,0]*p
return x
A = matrix([[1.0, 0.5],
[0.5, 1.0]])
b = matrix([[0.5],[0.0]])
x = ConjugateGradientMethod(A, b, 1.0e-5, 2)
print x
最終更新:2011年04月30日 02:41