ベアストウ法

ベアストウ法を勉強したのでPythonで書いてみた。
ベアストウ法を使うと、代数方程式にしか使えないけど複素解も含むすべての解を求めることができる。無理やり2次式で因数分解して積の形で表すって感じになる。詳しくはググってください。

コードを書く際、次のサイトを参考にした。
http://www.geocities.jp/supermisosan/bairstow.html
http://www.astro.phys.s.chiba-u.ac.jp/~miyaji/class/keisan-buturi/hb-kaisetsu.html

import math
import sys

def answer(b, c, ans):
    d = b * b - 4 * c    #判別式 以下dの値で場合分け
    if d > 0:
        re1 = (-1 * b + math.sqrt(d)) / 2.0    #2次方程式の解の公式を使う
        re2 = (-1 * b - math.sqrt(d)) / 2.0
        ans.append(str(re1))
        ans.append(str(re2))
    elif d == 0:
        re1 = -1 * b / 2.0
        ans.append(str(re1))
    else:
        re1 = re2 = -1 * b / 2.0
        im1 = math.sqrt(-1 * d) / 2.0
        im2 = -1 * im1
        ans.append(str(re1) + ' + ' + str(im1) + 'i')
        ans.append(str(re1) + ' - ' + str(im1) + 'i')

def bairstow(a, n):
    m = 0    #計算回数を数えるための変数
    b = range(n + 1)
    c = range(n + 1)
    u = v = 1.0    #2次因子x^2-u*x-vのu,vの初期値
    chk = 0.00000001    #誤差の範囲
    while True:
        m += 1
        if m > 1000000:    #指定回数を超えても収束しなかったら終了
            print ('Error!')
            sys.exit(1)
        b[0] = a[0]
        b[1] = a[1] - u * b[0]
        for i in range(2, n + 1):
            b[i] = a[i] - u * b[i - 1] - v * b[i - 2]
        c[0] = b[0]
        c[1] = b[1] - u * c[0]
        for i in range(2, n + 1):
            c[i] = b[i] - u * c[i - 1] - v * c[i - 2]
        h = c[n - 2] * c[n - 2] - c[n - 3] * (c[n - 1] - b[n - 1])
        du = (b[n - 1] * c[n - 2] - b[n] * c[n - 3]) / h
        dv = (b[n] * c[n - 2] - b[n - 1] * (c[n - 1] - b[n - 1])) / h
        u += du
        v += dv
        if du < chk and du > -1 * chk and  dv < chk and dv > -1 * chk:
            break
    for i in range(n - 1):
        a[i] = b[i]
    return u, v

def main():
    #(x-5)(x-4)(x-3)(x-2)(x-1)=0 を解く
    a = [1.0,-15.0,85.0,-225.0,274.0,-120.0]
    n = 5    #解く方程式の次数
    k = []
    p = a[0]
    q = a[1]
    while True:
        if n >= 3:
            p, q = bairstow(a, n)
            answer(p, q, k)
            n -= 2
        if n == 0:
            break
        if n == 1:
            res = -a[1] / a[0]
            k.append(str(res))
            break
        if n == 2:
            answer(p, q, k)
            break
    print k

if __name__ == '__main__':
    main()

実行結果はこちら。

['2.0', '1.0', '4.0', '3.0', '5.00000000044']

u,vを求める際、2変数のニュートン法を用いている。


なんかもう殆どコピーした感じになってて無理やり感が漂うけどまあ……