日: 2024年3月7日

ラゲール法による多項式の求解

概要

(長田 2024) の第4章代数方程式:1ラゲール法をもとに実装する。以下のサイトも参考にした。

  1. ラゲール法による求根アルゴリズム～Pythonコード付き～ – LASCODE

方針

根 1, 2, 3, 4, 5 を持つ5次の多項式を解くことにする. 

import numpy as np

p5deg = np.polynomial.polynomial.Polynomial.fromroots([1, 2, 3, 4, 5])
print(p5deg)

-120.0 + 274.0·x - 225.0·x² + 85.0·x³ - 15.0·x⁴ + 1.0·x⁵

実装

def laguerre(f, x0, delta=0.01, max_iter=100):
    ans = []

    f = f.copy()
    while 1 < f.degree():
        n = f.degree()
        df = f.deriv()
        ddf = df.deriv()

        x = x0
        iter = 0
        while not (np.abs(f(x)) < delta):
            numerator = n * f(x)
            denom_plus = df(x) + np.sqrt((n-1)**2 * df(x)**2 - n*(n-1)*f(x)*ddf(x))
            denom_minus = df(x) - np.sqrt((n-1)**2 * df(x)**2 - n*(n-1)*f(x)*ddf(x))

            denominator = denom_plus if np.abs(denom_minus) < np.abs(denom_plus) else denom_minus
            x -= numerator / denominator

            if np.abs(f(x)) < delta:
                break
            if max_iter < iter:
                break
            iter += 1

        ans.append(x)
        f = f // np.polynomial.polynomial.Polynomial.fromroots([x])

    if f.degree() == 1:
        ans.append(-f.coef[0])
    return ans

検証

5つの根の近似値が得られた。 

roots = laguerre(p5deg, x0=0, delta=0.00001)
print(roots)

[1.0000000000000002, 1.9999999999999907, 3.0000000000000107, 4.000000000000017, 4.999999999999981]

参考文献

長田 直樹. 2024. 数値解析 非線形方程式と数値積分. 単行本. 現代数学社. https://lead.to/amazon/jp/?op=bt&la=ja&key=4768706266.