SOR法

ガウス・ザイデル法の効率を上げたSOR法をやってみた。
・http://www.akita-nct.ac.jp/yamamoto/lecture/2004/5E/linear_equations/relaxation/html/node6.html
・SOR法 - Wikipedia

import math

def sor(matrix, dim, omega):
	count = 0
	check = 1.0e-10
	oldx = []
	for i in range(dim):
		oldx.append(1.0)   #解の初期値を設定する。今回は1にしてみた。
	newx = oldx[:]

	while True:
		error = []

		if count == 10000:
			print('Error.')
			return -1

		for i in range(dim):
			newx[i] = matrix[i][-1]
			for j in range(dim):
				if j != i:
					if j < i:
						newx[i] -= matrix[i][j] * newx[j]
					else:
						newx[i] -= matrix[i][j] * oldx[j]
			newx[i] /= matrix[i][i]
			newx[i] = (1 - omega) * oldx[i] + omega * newx[i]    #ガウス・ザイデル法と違うところ

		for i in range(dim):
			error.append(math.fabs(oldx[i] - newx[i]))

		count += 1

		if max(error) < check:
			print(newx)
			return count

		oldx = newx[:]

def output(filename, ans):
	f = open(filename, 'w')
	for i in ans:
		f.write('%f %d\n' % (i[0], i[1]))
	f.close()

def main():
	dim = 0
	matrix = []
	num = []
	ans = []
	omega = 0.1
	for line in open('matrix.txt'):
		item = line.split(' ')
		for i in item:
			num.append(float(i))
		matrix.append(num)
		num = []
		dim += 1
	print('Dimensions: %d' % dim)
	while omega <= 3.0:
		count = sor(matrix, dim, omega)
		ans.append([omega, count])
		omega += 0.1
	output('sor.txt', ans)

if __name__ == '__main__':
	main()

matrix.txtの中身。

1 -1 1 5
1 2 0 1
2 0 3 9

実行結果。

Dimensions: 3
[3.0000000009420646, -1.0000000011183832, 0.9999999985088222]
[3.0000000004246195, -1.0000000004868626, 0.9999999993508506]
[3.000000000269135, -1.0000000002967036, 0.9999999996043953]
[3.0000000002026237, -1.0000000002135532, 0.9999999997152625]
[3.0000000001049036, -1.0000000001049036, 0.9999999998601288]
[3.000000000072527, -1.0000000000681017, 0.9999999999091977]
[3.0000000000627702, -1.0000000000545017, 0.9999999999273311]
[3.000000000059581, -1.0000000000466576, 0.9999999999377899]
[3.000000000022313, -1.0000000000150144, 0.9999999999799807]
[3.0000000000042535, -1.0000000000021267, 0.9999999999971644]
[3.000000000002366, -1.0000000000015759, 0.9999999999981658]
[3.000000000015998, -1.0000000000105382, 0.9999999999894444]
[2.999999999989247, -1.0000000000048093, 0.9999999999960981]
[3.000000000000676, -0.9999999999831402, 1.0000000000194054]
[2.999999999990875, -1.000000000021526, 0.9999999999761493]
[3.0000000000242766, -1.0000000000216558, 0.9999999999840582]
[2.9999999999953153, -1.000000000002716, 1.0000000000095044]
[2.9999999999238147, -0.9999999999194701, 1.0000000000700444]
[3.0000000000412554, -1.0000000000403544, 0.9999999999770918]
Error.
Error.
Error.
Error.
Error.
Error.
Error.
Error.
Error.
Error.

sor.txtをgnuplotで棒グラフにして書いたものがこちら。収束しなかったものについては-1とした。

今回の場合、加速パラメータが1.1のとき反復回数が最も少ないことが分かる。
また、加速パラメータが2以上のときは発散していることが分かる。