# -*- coding: UTF-8 -*-
from math import *
import numpy as np
############################################
# 行列の固有値（フレーム法＋ベアストウ法）
# n : 次数
# ct : 最大繰り返し回数
# eps : 収束判定条件
# p0, q0 : x2+px+qにおけるp,qの初期値
# a : 係数（最高次から与え，値は変化する）
# b, c : 作業域（(n+1)次の配列）
# r : 結果
# A : 行列
# H1, H2 : 作業域（nxnの行列）
# return : =0 : 正常
# =1 : 収束せず
# coded by Y.Suganuma
############################################
def Frame(n, ct, eps, p0, q0, a, b, c, r, A, H1, H2) :
a[0] = 1.0
# a1の計算
a[1] = 0.0
for i1 in range(0, n) :
a[1] -= A[i1][i1]
# a2の計算
for i1 in range(0, n) :
for i2 in range(0, n) :
H1[i1][i2] = A[i1][i2]
H1[i1][i1] += a[1]
a[2] = 0.0
for i1 in range(0, n) :
for i2 in range(0, n) :
a[2] -= A[i1][i2] * H1[i2][i1]
a[2] *= 0.5
# a3からanの計算
for i1 in range(3, n+1) :
for i2 in range(0, n) :
for i3 in range(0, n) :
H2[i2][i3] = 0.0
for i4 in range(0, n) :
H2[i2][i3] += A[i2][i4] * H1[i4][i3]
H2[i2][i2] += a[i1-1]
a[i1] = 0.0
for i2 in range(0, n) :
for i3 in range(0, n) :
a[i1] -= A[i2][i3] * H2[i3][i2]
a[i1] /= i1
for i2 in range(0, n) :
for i3 in range(0, n) :
H1[i2][i3] = H2[i2][i3]
# ベアストウ法
ind = Bairstow(n, ct, eps, p0, q0, a, b, c, r, 0)
return ind
############################################
# 実係数代数方程式の解（ベアストウ法）
# n : 次数
# ct : 最大繰り返し回数
# eps : 収束判定条件
# p0, q0 : x2+px+qにおけるp,qの初期値
# a : 係数（最高次から与え，値は変化する）
# b,c : 作業域（(n+1)次の配列）
# r : 結果
# k : 結果の位置
# return : =0 : 正常
# =1 : 収束せず
# coded by Y.Suganuma
############################################
def Bairstow(n, ct, eps, p0, q0, a, b, c, r, k) :
p1 = p0
p2 = 0.0
q1 = q0
q2 = 0.0
ind = 0
count = 0
# １次の場合
if n == 1 :
if abs(a[0]) < eps :
ind = 1
else :
r[k] = complex(-a[1] / a[0], 0)
# ２次の場合
elif n == 2 :
# １次式
if abs(a[0]) < eps :
if abs(a[1]) < eps :
ind = 1
else :
r[k] = complex(-a[2] / a[1], 0)
# ２次式
else :
D = a[1] * a[1] - 4.0 * a[0] * a[2]
if D < 0.0 : # 虚数
D = sqrt(-D)
a[0] *= 2.0
r[k] = complex(-a[1] / a[0], D / a[0])
r[k+1] = complex(-a[1] / a[0], -D / a[0])
else : # 実数
D = sqrt(D)
a[0] = 1.0 / (2.0 * a[0])
r[k] = complex(a[0] * (-a[1] + D), 0)
r[k+1] = complex(a[0] * (-a[1] - D), 0)
# ３次以上の場合
else :
# 因数分解
ind = 1
while ind > 0 and count <= ct :
for i1 in range(0, n+1) :
if i1 == 0 :
b[i1] = a[i1]
elif i1 == 1 :
b[i1] = a[i1] - p1 * b[i1-1]
else :
b[i1] = a[i1] - p1 * b[i1-1] - q1 * b[i1-2]
for i1 in range(0, n+1) :
if i1 == 0 :
c[i1] = b[i1]
elif i1 == 1 :
c[i1] = b[i1] - p1 * c[i1-1]
else :
c[i1] = b[i1] - p1 * c[i1-1] - q1 * c[i1-2]
D = c[n-2] * c[n-2] - c[n-3] * (c[n-1] - b[n-1])
if fabs(D) < eps :
return ind
else :
dp = (b[n-1] * c[n-2] - b[n] * c[n-3]) / D
dq = (b[n] * c[n-2] - b[n-1] * (c[n-1] - b[n-1])) / D
p2 = p1 + dp
q2 = q1 + dq
if abs(dp) < eps and fabs(dq) < eps :
ind = 0
else :
count += 1
p1 = p2
q1 = q2
if ind == 0 :
# ２次方程式を解く
D = p2 * p2 - 4.0 * q2
if D < 0.0 : # 虚数
D = sqrt(-D)
r[k] = complex(-0.5 * p2, 0.5 * D)
r[k+1] = complex(-0.5 * p2, -0.5 * D)
else : # 実数
D = sqrt(D)
r[k] = complex(0.5 * (-p2 + D), 0)
r[k+1] = complex(0.5 * (-p2 - D), 0)
# 残りの方程式を解く
n -= 2
for i1 in range(0, n+1) :
a[i1] = b[i1]
ind = Bairstow(n, ct, eps, p0, q0, a, b, c, r, k+2)
return ind
----------------------------------
# -*- coding: UTF-8 -*-
import numpy as np
from math import *
from function import Bairstow, Frame
############################################
# 行列の固有値（フレーム法＋ベアストウ法） */
# coded by Y.Suganuma */
############################################
# データの設定
ct = 1000
eps = 1.0e-10
p0 = 0.0
q0 = 0.0
n = 3
a = np.empty(n+1, np.float)
b = np.empty(n+1, np.float)
c = np.empty(n+1, np.float)
r = np.empty(n, np.complex)
A = np.array([[7, 2, 1],[2, 1, -4],[1, -4, 2]], np.float)
H1 = np.empty((n, n), np.float)
H2 = np.empty((n, n), np.float)
# 計算
ind = Frame(n, ct, eps, p0, q0, a, b, c, r, A, H1, H2)
# 出力
if ind > 0 :
print("収束しませんでした！")
else :
for i1 in range(0, n) :
print(" " + str(r[i1]))