跳轉到

4.3 多項式求根#

概念說明#

多項式 p(x) = aₙxⁿ + … + a₁x + a₀ 是最特殊的一類函數：它由有限個係數完全決定，求值、微分、積分都有封閉公式，而且代數基本定理保證它恰有 n 個（計入重數的）複數根。這些結構讓多項式求根可以比一般函數求根更「系統化」：

  1. 利用合成除法（秦九韶/Horner 法）：以 O(n) 時間求 p(x₀) 的值，順便得到商式，可用於「降階」——找到一個根後除掉對應因式，繼續找下一個根。
  2. 實根的定位與隔離：斯圖姆定理（Sturm's theorem）或 Descartes 符號法則可以確定某區間內實根的個數，把根一個個「關進」互不重疊的區間，再交給牛頓法或二分法精化。
  3. 全部根一次解：Durand-Kerner 法（又稱 Weierstrass 法）同時迭代 n 個近似根，每一步是簡單的代數公式；或直接呼叫 numpy.roots，其內部使用伴隨矩陣的特徵值分解。

與一般求根相比，多項式求根要特別小心兩件事：重根（收斂變慢且數值敏感）與病態性（係數的微小擾動可能讓根大幅移動，Wilkinson 多項式是經典例子）。

演算法步驟#

以「Horner 求值 + 降階 + 牛頓法」找多項式全部實根為例：

  1. 用 Horner 法實作 p(x) 與 p′(x) 的求值。
  2. 用符號法則估計正實根個數；對實係數多項式，實根都落在 [−R, R]，其中 R = 1 + max|aᵢ/aₙ|（柯西界）。
  3. 在 [−R, R] 上掃描取樣找變號區間，每個變號區間用牛頓法（或二分法）精化出一個根。
  4. 找到根 r 後做降階：p ← p/(x − r)（用合成除法），回到第 3 步直到次數為 1。

Python 程式碼範例#

import numpy as np

def horner(coeffs, x):
    """coeffs = [a_n, ..., a_1, a_0]，回傳 p(x)"""
    y = 0.0
    for c in coeffs:
        y = y * x + c
    return y

def deflate(coeffs, r):
    """以合成除法除以 (x - r)，回傳新係數"""
    out, carry = [], 0.0
    for c in coeffs[:-1]:
        carry = c + carry * r
        out.append(carry)
    return out

def poly_real_roots(coeffs, tol=1e-10):
    roots = []
    cs = list(map(float, coeffs))
    while len(cs) > 2:
        # 在 [-R, R] 掃描找變號區間
        R = 1 + max(abs(c) for c in cs[:-1]) / abs(cs[0])
        xs = np.linspace(-R, R, 2001)
        vals = [horner(cs, x) for x in xs]
        r = None
        for i in range(len(xs) - 1):
            if vals[i] == 0:
                r = xs[i]; break
            if vals[i] * vals[i+1] < 0:
                # 牛頓法精化（導數係數由冪法則生成）
                dc = [cs[j]*(len(cs)-1-j) for j in range(len(cs)-1)]
                x = 0.5*(xs[i] + xs[i+1])
                for _ in range(60):
                    fx, dfx = horner(cs, x), horner(dc, x)
                    if abs(dfx) < 1e-14: break
                    x -= fx/dfx
                r = x; break
        if r is None:
            break
        roots.append(r)
        cs = deflate(cs, r)
    return sorted(roots)

cs = [1, -5, 5, -1]          # p(x) = x^3 - 5x^2 + 5x - 1
print("自寫法:", poly_real_roots(cs))
print("numpy :", np.roots(cs))
print("sympy :", __import__("sympy").solve(__import__("sympy").
      Symbol("x")**3 - 5*__import__("sympy").Symbol("x")**2
      + 5*__import__("sympy").Symbol("x") - 1))

三種方法都應給出根 {1, (5±√13)/2}。注意 numpy.roots 會連複數根一起回傳，而符號解（sympy）才是數學上的精確答案——這正是第 1 章「精確 vs 近似」光譜的具體展現。

複雜度與誤差分析#

Horner 求值需要 n 次乘法與 n 次加法，時間 O(n)，這是漸近最優的。掃描加牛頓的總成本約為 O(n²)（每個根一次降階，降階本身 O(n)，根最多 n 個）。Durand-Kerner 每步 O(n²)，通常數十步收斂。誤差方面，多項式求根的條件數可能非常大：對 n 次多項式，係數擾動 δ 導致的根擾動可放大到 O(δ · 10ⁿ)（Wilkinson 現象），因此高次多項式不要自己冪基底展開，改用切比雪夫基底或對根直接參數化。降階法還有誤差累積問題：每除一次，後面根的精度就下降一級，所以精密計算應改用同時迭代法或在最後對原多項式做一步牛頓「拋光」。

練習#

  1. 用 Horner 法分別在 x = 1.5 與 x = 10 處求 p(x) = x⁵ − 3x⁴ + 2x − 1 的值，並與 np.polyval 比對；順便統計你的程式做了幾次乘法。
  2. 對 Wilkinson 多項式的 n=8 版本（根為 1, 2, …, 8）：先用 np.poly 由根重建係數，再用 np.roots 求回來，觀察根的誤差有多大，並解釋病因。
  3. 實作 Durand-Kerner 法：x_i ← x_i − p(x_i) / ∏_{j≠i} (x_i − x_j)，對 p(x) = x⁴ − 1 測試，驗證四個根 {±1, ±i} 都被找出來。