ベアストウ法
ベアストウ法を勉強したので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変数のニュートン法を用いている。
なんかもう殆どコピーした感じになってて無理やり感が漂うけどまあ……