#*******************************************/
# 行列の固有値（フレーム法＋ベアストウ法） */
# coded by Y.Suganuma */
#*******************************************/
#************************************************/
# 行列の固有値（フレーム法＋ベアストウ法） */
# n : 次数 */
# ct : 最大繰り返し回数 */
# eps : 収束判定条件 */
# p0, q0 : x2+px+qにおけるp,qの初期値 */
# a : 係数（最高次から与え，値は変化する） */
# b, c : 作業域（(n+1)次の配列） */
# rl, im : 結果の実部と虚部 */
# aa : 行列 */
# h1, h2 : 作業域（nxnの行列） */
# return : =0 : 正常 */
# =1 : 収束せず */
# coded by Y.Suganuma */
#************************************************/
def Frame(n, ct, eps, p0, q0, a, b, c, rl, im, aa, h1, h2)
a[0] = 1.0
# a1の計算
a[1] = 0.0
for i1 in 0 ... n
a[1] -= aa[i1][i1]
end
# a2の計算
for i1 in 0 ... n
for i2 in 0 ... n
h1[i1][i2] = aa[i1][i2]
end
h1[i1][i1] += a[1]
end
a[2] = 0.0
for i1 in 0 ... n
for i2 in 0 ... n
a[2] -= aa[i1][i2] * h1[i2][i1]
end
end
a[2] *= 0.5
# a3からanの計算
for i1 in 3 ... n+1
for i2 in 0 ... n
for i3 in 0 ... n
h2[i2][i3] = 0.0
for i4 in 0 ... n
h2[i2][i3] += aa[i2][i4] * h1[i4][i3]
end
end
h2[i2][i2] += a[i1-1]
end
a[i1] = 0.0
for i2 in 0 ... n
for i3 in 0 ... n
a[i1] -= aa[i2][i3] * h2[i3][i2]
end
end
a[i1] /= i1
for i2 in 0 ... n
for i3 in 0 ... n
h1[i2][i3] = h2[i2][i3]
end
end
end
# ベアストウ法
ind = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, 0)
return ind
end
#************************************************/
# 実係数代数方程式の解（ベアストウ法） */
# n : 次数 */
# ct : 最大繰り返し回数 */
# eps : 収束判定条件 */
# p0, q0 : x2+px+qにおけるp,qの初期値 */
# a : 係数（最高次から与え，値は変化する） */
# b,c : 作業域（(n+1)次の配列） */
# rl, im : 結果の実部と虚部 */
# k : 結果の位置 */
# return : =0 : 正常 */
# =1 : 収束せず */
# coded by Y.Suganuma */
#************************************************/
def Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, k)
# 初期設定
p1 = p0
p2 = 0.0
q1 = q0
q2 = 0.0
ind = 0
count = 0
#
# １次の場合
#
if n == 1
if a[0].abs() < eps
ind = 1
else
rl[k] = -a[1] / a[0]
im[k] = 0.0
end
#
# ２次の場合
#
elsif n == 2
# １次式
if a[0].abs() < eps
if a[1].abs() < eps
ind = 1
else
rl[k] = -a[2] / a[1]
im[k] = 0.0
end
# ２次式
else
d = a[1] * a[1] - 4.0 * a[0] * a[2]
if d < 0.0 # 虚数
d = Math.sqrt(-d)
a[0] *= 2.0
rl[k] = -a[1] / a[0]
rl[k+1] = -a[1] / a[0]
im[k] = d / a[0]
im[k+1] = -im[0]
else # 実数
d = Math.sqrt(d)
a[0] = 1.0 / (2.0 * a[0])
rl[k] = a[0] * (-a[1] + d)
rl[k+1] = a[0] * (-a[1] - d)
im[k] = 0.0
im[k+1] = 0.0
end
end
# ３次以上の場合
else
# 因数分解
ind = 1
while ind > 0 && count <= ct
for i1 in 0 ... n+1
if i1 == 0
b[i1] = a[i1]
elsif i1 == 1
b[i1] = a[i1] - p1 * b[i1-1]
else
b[i1] = a[i1] - p1 * b[i1-1] - q1 * b[i1-2]
end
end
for i1 in 0 ... n+1
if i1 == 0
c[i1] = b[i1]
elsif i1 == 1
c[i1] = b[i1] - p1 * c[i1-1]
else
c[i1] = b[i1] - p1 * c[i1-1] - q1 * c[i1-2]
end
end
d = c[n-2] * c[n-2] - c[n-3] * (c[n-1] - b[n-1])
if d.abs() < 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 dp.abs() < eps && dq.abs() < eps
ind = 0
else
count += 1
p1 = p2
q1 = q2
end
end
end
if ind == 0
# ２次方程式を解く
d = p2 * p2 - 4.0 * q2
if d < 0.0 # 虚数
d = Math.sqrt(-d)
rl[k] = -0.5 * p2
rl[k+1] = -0.5 * p2
im[k] = 0.5 * d
im[k+1] = -im[k]
else # 実数
d = Math.sqrt(d)
rl[k] = 0.5 * (-p2 + d)
rl[k+1] = 0.5 * (-p2 - d)
im[k] = 0.0
im[k+1] = 0.0
end
# 残りの方程式を解く
n -= 2
for i1 in 0 ... n+1
a[i1] = b[i1]
end
ind = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, k+2)
end
end
return ind
end
# データの設定
ct = 1000
eps = 1.0e-10
p0 = 0.0
q0 = 0.0
n = 3
a = Array.new(n+1)
b = Array.new(n+1)
c = Array.new(n+1)
rl = Array.new(n)
im = Array.new(n)
aa = Array.new(n)
h1 = Array.new(n)
h2 = Array.new(n)
for i1 in 0 ... n
aa[i1] = Array.new(n)
h1[i1] = Array.new(n)
h2[i1] = Array.new(n)
end
aa[0][0] = 7.0
aa[0][1] = 2.0
aa[0][2] = 1.0
aa[1][0] = 2.0
aa[1][1] = 1.0
aa[1][2] = -4.0
aa[2][0] = 1.0
aa[2][1] = -4.0
aa[2][2] = 2.0
# 計算
ind = Frame(n, ct, eps, p0, q0, a, b, c, rl, im, aa, h1, h2)
# 出力
if ind > 0
printf("収束しませんでした！\n")
else
for i1 in 0 ... n
printf(" %f i %f\n", rl[i1], im[i1])
end
end