実対称行列の固有値・固有ベクトル（ヤコビ法）

# -*- coding: UTF-8 -*-
from math import *
import numpy as np

############################################
# 実対称行列の固有値・固有ベクトル（ヤコビ法）
#      n : 次数
#      ct : 最大繰り返し回数
#      eps : 収束判定条件
#      A : 対象とする行列
#      A1, A2 : 作業域（nxnの行列），A1の対角要素が固有値
#      X1, X2 : 作業域（nxnの行列），X1の各列が固有ベクトル
#      return : =0 : 正常
#               =1 : 収束せず
#      coded by Y.Suganuma
############################################

def Jacobi(n, ct, eps, A, A1, A2, X1, X2) :
			# 初期設定
	k   = 0
	ind = 1
	p   = 0
	q   = 0
	for i1 in range(0, n) :
		for i2 in range(0, n) :
			A1[i1][i2] = A[i1][i2]
			X1[i1][i2] = 0.0
		X1[i1][i1] = 1.0
			# 計算
	while ind > 0 and k < ct :
				# 最大要素の探索
		max = 0.0
		for i1 in range(0, n) :
			for i2 in range(0, n) :
				if i2 != i1 :
					if abs(A1[i1][i2]) > max :
						max = fabs(A1[i1][i2])
						p   = i1
						q   = i2
				# 収束判定
					# 収束した
		if max < eps :
			ind = 0
					# 収束しない
		else :
						# 準備
			s  = -A1[p][q]
			t  = 0.5 * (A1[p][p] - A1[q][q])
			v  = fabs(t) / sqrt(s * s + t * t)
			sn = sqrt(0.5 * (1.0 - v))
			if s*t < 0.0 :
				sn = -sn
			cs = sqrt(1.0 - sn * sn)
						# Akの計算
			for i1 in range(0, n) :
				if i1 == p :
					for i2 in range(0, n) :
						if i2 == p :
							A2[p][p] = A1[p][p] * cs * cs + A1[q][q] * sn * sn - 2.0 * A1[p][q] * sn * cs
						elif i2 == q :
							A2[p][q] = 0.0
						else :
							A2[p][i2] = A1[p][i2] * cs - A1[q][i2] * sn
				elif i1 == q :
					for i2 in range(0, n) :
						if i2 == q :
							A2[q][q] = A1[p][p] * sn * sn + A1[q][q] * cs * cs + 2.0 * A1[p][q] * sn * cs
						elif i2 == p :
							A2[q][p] = 0.0
						else :
							A2[q][i2] = A1[q][i2] * cs + A1[p][i2] * sn
				else :
					for i2 in range(0, n) :
						if i2 == p :
							A2[i1][p] = A1[i1][p] * cs - A1[i1][q] * sn
						elif i2 == q :
							A2[i1][q] = A1[i1][q] * cs + A1[i1][p] * sn
						else :
							A2[i1][i2] = A1[i1][i2]
						# Xkの計算
			for i1 in range(0, n) :
				for i2 in range(0, n) :
					if i2 == p :
						X2[i1][p] = X1[i1][p] * cs - X1[i1][q] * sn
					elif i2 == q :
						X2[i1][q] = X1[i1][q] * cs + X1[i1][p] * sn
					else :
						X2[i1][i2] = X1[i1][i2]
						# 次のステップへ
			k += 1
			for i1 in range(0, n) :
				for i2 in range(0, n) :
					A1[i1][i2] = A2[i1][i2]
					X1[i1][i2] = X2[i1][i2]

	return ind

----------------------------------

# -*- coding: UTF-8 -*-
import numpy as np
from math import *
from function import Jacobi

############################################
# 実対称行列の固有値・固有ベクトル（ヤコビ法）
#      coded by Y.Suganuma
############################################
			# データの設定
ct  = 1000
eps = 1.0e-10
n   = 3
A   = np.array([[1, 0, -1],[0, 1, -1],[-1, -1, 2]], np.float)
A1  = np.empty((n, n), np.float)
A2  = np.empty((n, n), np.float)
X1  = np.empty((n, n), np.float)
X2  = np.empty((n, n), np.float)
			# 計算
ind = Jacobi(n, ct, eps, A, A1, A2, X1, X2)
			# 出力
if ind > 0 :
	print("収束しませんでした！")
else :
	print("固有値")
	print(" ", A1[0][0], A1[1][1], A1[2][2])
	print("固有ベクトル（各列が固有ベクトル）")
	for i1 in range(0, n) :
		for i2 in range(0, n) :
			print(" " + str(X1[i1][i2]), end="")
		print()