ヤコビ法（連立方程式）

ヤコビ法をやってみた。固有値を求める方ではなく連立方程式を解く方の。
次のサイトを参考にした。
・http://next1.cc.it-hiroshima.ac.jp/MULTIMEDIA/numeanal2/node13.html

import math

def jacobi(a, b, oldx):
	newx = [0.0,0.0,0.0,0.0]
	d = len(oldx)
	count = 0
	check = 1.0e-20
	while True:
		error = []
		if count == 10000:    #収束しなかった場合はここでループを抜ける
			print('Error.')
			break
		for j in range(d):
			newx[j] = b[j]
			for k in range(d):
				if k != j:
					newx[j] -= a[j][k] * oldx[k]
			newx[j] = newx[j] / a[j][j]
		for j in range(d):
			error.append(math.fabs(oldx[j] - newx[j]))    #それぞれの誤差をリストにまとめる
			oldx[j] = newx[j]
		count += 1
		if max(error) < check:    #リスト内で一番大きいのが設定値より低くなったら計算終了
			print('Count: %d' % count)
			return

def main():
	a = [[9.0,2.0,1.0,1.0],[2.0,8.0,-2.0,1.0],[-1.0,-2.0,7.0,-2.0],[1.0,-1.0,-2.0,6.0]]    #係数行列
	b = [20,16,8,17]    #各式の解
	ans = [0.0, 0.0, 0.0, 0.0]    #結果を入れるリスト
	jacobi(a, b, ans)
	print ans

if __name__ == '__main__':
	main()

実行結果はこちら。

Count: 48
[1.0, 2.0, 3.0, 4.0]