跳轉到

7.1 拉格朗日插值#

概念說明#

給定 n+1 個資料點 (x₀, y₀), …, (xₙ, yₙ)（xᵢ 互異），是否存在唯一的多項式 p，次數 ≤ n，通過所有點？答案是肯定的，而且有兩種經典建構法：牛頓的差商法（適合動態加點）與拉格朗日法（公式最優雅）。

拉格朗日的核心想法是基底函數：先造出 n+1 個「開關函數」

Lᵢ(x) = ∏_{j≠i} (x − xⱼ)/(xᵢ − xⱼ)

每個 Lᵢ 在 xᵢ 處等於 1、在其他所有節點 xⱼ（j≠i）處等於 0。於是 p(x) = Σ yᵢLᵢ(x) 自動通過每一點——在 xₖ 處，只有 Lₖ 不為零，恰好貢獻 yₖ。唯一性則由代數事實保證：次數 ≤ n 的多項式若有 n+1 個相異根，必為零多項式。

但高次插值有個惡名昭彰的陷阱：Runge 現象。對 f(x) = 1/(1+25x²) 在 [−1, 1] 上取等距節點做高次插值，次數越高，曲線在兩端震盪得越兇（n = 20 時震幅可達真值的數十倍）。罪魁不是拉格朗日公式本身，而是節點分布：解藥是把節點集中在端點附近，即切比雪夫節點 xᵢ = cos((2i+1)π/(2n+2))。採用切比雪夫節點後，誤差一致有界且隨 n 指數下降，插值從「危險品」變成「精密工具」。

演算法步驟#

  1. 輸入節點 xᵢ 與函數值 yᵢ。
  2. 對每個 i 建構 Lᵢ（分母先算好，避免重複乘除）。
  3. 對查詢點 x：計算 Σ yᵢ·Lᵢ(x)。實作上用重心法（barycentric）更穩更快： p(x) = [Σ wᵢyᵢ/(x−xᵢ)] / [Σ wᵢ/(x−xᵢ)]，其中 wᵢ = 1/∏_{j≠i}(xᵢ−xⱼ)。

Python 程式碼範例#

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

def lagrange_eval(xs, ys, x):
    """樸素拉格朗日求值，O(n^2)"""
    n = len(xs)
    total = 0.0
    for i in range(n):
        term = ys[i]
        for j in range(n):
            if j != i:
                term *= (x - xs[j]) / (xs[i] - xs[j])
        total += term
    return total

def barycentric_weights(xs):
    n = len(xs)
    w = np.ones(n)
    for j in range(n):
        diff = xs[j] - np.delete(xs, j)
        w[j] = 1.0 / np.prod(diff)
    return w

def barycentric_eval(xs, ys, w, x):
    diff = x - xs
    hit = np.where(np.abs(diff) < 1e-14)[0]
    if len(hit):                       # x 恰為節點
        return ys[hit[0]]
    return np.sum(w * ys / diff) / np.sum(w / diff)

f = lambda x: 1 / (1 + 25 * x**2)      # Runge 函數

for n in [8, 16]:
    xe = np.linspace(-1, 1, n + 1)                    # 等距節點
    i = np.arange(n + 1)
    xc = np.cos((2 * i + 1) * np.pi / (2 * n + 2))    # 切比雪夫節點
    xx = np.linspace(-1, 1, 400)
    for name, xs in [("等距", xe), ("切比雪夫", xc)]:
        ys = f(xs)
        w = barycentric_weights(np.sort(xs))
        xs_s = np.sort(xs)
        approx = np.array([barycentric_eval(xs_s, f(xs_s), w, t) for t in xx])
        print(f"n={n:2d} {name}節點 最大誤差 = {np.abs(approx - f(xx)).max():.3e}")

# 與 sympy 驗證唯一性：三點插值
import sympy as sp
X = sp.symbols("x")
p = sp.interpolate([(0, 1), (1, 3), (2, 2)], X)
print("通過 (0,1),(1,3),(2,2) 的多項式:", sp.expand(p))

誤差分析#

插值誤差有精確公式：

f(x) − p(x) = f⁽ⁿ⁺¹⁾(ξ)/(n+1)! · ∏(x − xᵢ)

誤差由兩個因子相乘：函數的高階導數（多光滑）與節點多項式 ∏(x−xᵢ)（節點怎麼分布）。等距節點時，節點多項式在區間兩端的值爆炸性成長（Lebesgue 常數 ~2ⁿ），這就是 Runge 現象的代數根源；切比雪夫節點把節點多項式的最大值壓到 2^(−n)，使誤差 ≈ O(2⁻ⁿ/(n+1)!)，呈譜級收斂。重心公式的優勢在此浮現：換節點不用重算基底，權重 wᵢ 只跟節點有關；對切比雪夫節點甚至有封閉權重公式。計算複雜度：重心法求一點 O(n)、建權重 O(n²)；樸素拉格朗日每點 O(n²)。牛頓差商法另有一個實用性質：新增一個資料點只需加一項，不必從頭重算。

練習#

  1. 對 Runge 函數畫出 n = 4, 8, 16, 24 的等距插值曲線（matplotlib），目視 Runge 現象；同圖畫切比雪夫版本對照。
  2. 實作牛頓差商插值（建差商表 O(n²)，求值用巢狀乘法 O(n)），新增一個節點後驗證「舊結果不變、只需加一項」。
  3. 用插值證明平方根表：對 f(x)=√x 在 x = 1, 1.21, 1.44 三點做二次插值，估計 √1.3，與真值比較，並用誤差公式解釋精度為何約 10⁻³。