最大固有値と固有ベクトル（べき乗法）

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

############################################
# 最大固有値と固有ベクトル（べき乗法）
#      i : 何番目の固有値かを示す
#      n : 次数
#      m : 丸め誤差の修正回数
#      ct : 最大繰り返し回数
#      eps : 収束判定条件
#      A : 対象とする行列
#      B : nxnの行列，最初は，単位行列
#      C : 作業域，nxnの行列
#      r : 固有値
#      v : 各行が固有ベクトル（nxn）
#      v0 : 固有ベクトルの初期値
#      v1,v2 : 作業域（n次元ベクトル）
#      return : 求まった固有値の数
#      coded by Y.Suganuma
############################################

def power(i, n, m, ct, eps, A, B, C, r, v, v0, v1, v2) :
			# 初期設定
	ind = i
	k   = 0
	l1  = 0.0
	for i1 in range(0, n) :
		l1     += v0[i1] * v0[i1]
		v1[i1]  = v0[i1]
	l1 = sqrt(l1)
			# 繰り返し計算
	while k < ct :
				# 丸め誤差の修正
		if k%m == 0 :
			l2 = 0.0
			for i1 in range(0, n) :
				v2[i1] = 0.0
				for i2 in range(0, n) :
					v2[i1] += B[i1][i2] * v1[i2]
				l2 += v2[i1] * v2[i1]
			l2 = sqrt(l2)
			for i1 in range(0, n) :
				v1[i1] = v2[i1] / l2
				# 次の近似
		l2 = 0.0
		for i1 in range(0, n) :
			v2[i1] = 0.0
			for i2 in range(0, n) :
				v2[i1] += A[i1][i2] * v1[i2]
			l2 += v2[i1] * v2[i1]
		l2 = sqrt(l2)
		for i1 in range(0, n) :
			v2[i1] /= l2
				# 収束判定
					# 収束した場合
		if abs((l2-l1)/l1) < eps :
			k1 = -1
			for i1 in range(0, n) :
				if abs(v2[i1]) > 0.001 :
					k1 = i1
					if v2[k1]*v1[k1] < 0.0 :
						l2 = -l2
					break
			k    = ct
			r[i] = l2
			for i1 in range(0, n) :
				v[i][i1] = v2[i1]
			if i == n-1 :
				ind = i + 1
			else :
				for i1 in range(0, n) :
					for i2 in range(0, n) :
						C[i1][i2] = 0.0
						for i3 in range(0, n) :
							if i1 == i3 :
								x = A[i1][i3] - l2
							else :
								x = A[i1][i3]
							C[i1][i2] += x * B[i3][i2]
				for i1 in range(0, n) :
					for i2 in range(0, n) :
						B[i1][i2] = C[i1][i2]
				for i1 in range(0, n) :
					v1[i1] = 0.0
					for i2 in range(0, n) :
						v1[i1] += B[i1][i2] * v0[i2]
				for i1 in range(0, n) :
					v0[i1] = v1[i1]
				ind = power(i+1, n, m, ct, eps, A, B, C, r, v, v0, v1, v2)
					# 収束しない場合
		else :
			for i1 in range(0, n) :
				v1[i1] = v2[i1]
			l1 = l2
			k += 1

	return ind

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

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

############################################
# 最大固有値と固有ベクトル（べき乗法）
#      coded by Y.Suganuma
############################################
			# データの設定
ct  = 200
eps = 1.0e-10
n   = 3
m   = 15
A   = np.array([[11, 6, -2],[-2, 18, 1],[-12, 24, 13]], np.float)
B   = np.array([[1, 0, 0],[0, 1, 0],[0, 0, 1]], np.float)
C   = np.empty((n, n), np.float)
v   = np.empty((n, n), np.float)
r   = np.empty(n, np.float)
v0  = np.array([1, 1, 1], np.float)
v1  = np.empty(n, np.float)
v2  = np.empty(n, np.float)
			# 計算
ind = power(0, n, m, ct, eps, A, B, C, r, v, v0, v1, v2)
			# 出力
if ind == 0 :
	print("固有値が求まりませんでした！")
else :
	for i1 in range(0, ind) :
		print("固有値 " + str(r[i1]), end="")
		print(" 固有ベクトル ", end="")
		for i2 in range(0, n) :
			print(" " + str(v[i1][i2]), end="")
		print()