Python 使用 Laguerre 方法求多項式所有根

求解多項式根是科學計算常見問題。Laguerre 方法是一種強大的迭代演算法，能同時找到多項式的所有實根與複根。本文介紹演算法原理，並提供完整 Python 實作，包含多項式求值、去多項式化與根的提取流程。

Laguerre 方法原理

對多項式 P(x) = a[0] + a[1]x + ... + a[n]x^n，Laguerre 迭代公式為：

x_{k+1} = x_k - n / (G ± sqrt((n-1)(nH - G^2)))

其中 G = P'(x)/P(x)，H = G^2 - P''(x)/P(x)。選擇分母較大的符號以增強穩定性。迭代直到 |dx| < tol 或 |P(x)| < tol。

完整 Python 實作

程式碼包含三個核心函數：evalPoly（多項式求值與導數）、deflPoly（去多項式化）、polyRoots（主流程）。

from numpy import zeros, complex
from cmath import sqrt
from random import random

def evalPoly(a, x):
    """回傳 P(x), P'(x), P''(x)"""
    n = len(a) - 1
    p = a[n]
    dp = 0.0 + 0.0j
    ddp = 0.0 + 0.0j
    for i in range(n-1, -1, -1):
        ddp = ddp * x + 2 * dp
        dp = dp * x + p
        p = p * x + a[i]
    return p, dp, ddp

def deflPoly(a, root):
    """以 (x - root) 去多項式化，回傳降次係數陣列"""
    n = len(a) - 1
    b = [(0.0 + 0.0j)] * n
    b[n-1] = a[n]
    for i in range(n-2, -1, -1):
        b[i] = a[i+1] + root * b[i+1]
    return b

def polyRoots(a, tol=1.0e-12):
    """回傳多項式所有根的陣列"""
    n = len(a) - 1
    roots = zeros(n, dtype=complex)
    for i in range(n):
        # Laguerre 迭代尋找一根
        x = random()  # 隨機初始值
        for _ in range(30):
            p, dp, ddp = evalPoly(a, x)
            if abs(p) < tol:
                break
            g = dp / p
            h = g*g - ddp / p
            f = sqrt((n - 1) * (n * h - g*g))
            if abs(g + f) > abs(g - f):
                dx = n / (g + f)
            else:
                dx = n / (g - f)
            x -= dx
            if abs(dx) < tol:
                break
        # 若虛部極小視為實根
        if abs(x.imag) < tol:
            x = x.real
        roots[i] = x
        # 去多項式化
        a = deflPoly(a, x)
    return roots

使用範例

求解 x^3 - 6x^2 + 11x - 6 = 0（根為 1, 2, 3）：

coeff = [-6, 11, -6, 1]  # 係數由低次到高次
roots = polyRoots(coeff)
print("根:", roots)  # 輸出: [1. 2. 3.]

注意事項

• 係數陣列 a 索引 0 為常數項，依次遞增。
• 隨機初始值可能導致不同執行順序，但結果相同。
• 高次多項式（>10）可能面臨數值不穩定，建議使用成熟函式庫（如 NumPy roots）。

此實作展示經典數值演算法的 Python 轉譯，適合教學與中小規模問題。