LU分解

そういえばLU分解やってなかった事に気づいたのでやってみた。つぎのサイトを参考にした。
・http://next1.cc.it-hiroshima.ac.jp/MULTIMEDIA/numeanal1/node32.html
・LU分解 - [物理のかぎしっぽ]

import copy

def lu(matrix, dim):
	up = copy.deepcopy(matrix)    #元のmatrixを壊さないためにコピーする

	low = []    #下三角行列を作る(最初は対角成分のみ1で他は0)
	for i in range(dim):
		tmp = []
		for j in range(dim):
			if i == j:
				tmp.append(1.0)
			else:
				tmp.append(0.0)
		low.append(tmp)

	for i in range(dim):    #ここからLU分解
		p = up[i][i]
		for j in range(i+1, dim):
			q = up[j][i]
			low[j][i] = q / p
			for k in range(dim):
				up[j][k] -= up[i][k] * q / p

	return up, low

def main():
	dim = 0
	matrix = []
	num = []
	for line in open('matrix.txt'):
		item = line.split(' ')
		for i in item:
			num.append(float(i))
		matrix.append(num)
		num = []
		dim += 1
	u, l = lu(matrix, dim)
	for i in u:
		print i
　　　　print('\n')
	for i in l:
		print i

if __name__ == '__main__':
    main()

例としてこんな感じのmatrix.txtを用意した。

1 1 0 3
2 1 -1 1
3 -1 -1 2
-1 2 3 -1

実行結果はこちら。

[1.0, 1.0, 0.0, 3.0]
[0.0, -1.0, -1.0, -5.0]
[0.0, 0.0, 3.0, 13.0]
[0.0, 0.0, 0.0, -13.0]


[1.0, 0.0, 0.0, 0.0]
[2.0, 1.0, 0.0, 0.0]
[3.0, 4.0, 1.0, 0.0]
[-1.0, -3.0, 0.0, 1.0]

上半分が上三角行列、下半分が下三角行列。