Python 實作對角線有理函數插值演算法

有理函數插值是數值分析中常用的技術，特別適合處理具有極點或快速變化的函數。對角線有理函數插值則是將分子與分母多項式的次數設為相同，能在較少節點下獲得較高精度。本文介紹如何使用 Python 實作此演算法，並提供完整程式碼供參考。

演算法原理

給定數據點 (xData[i], yData[i])，對角線有理函數插值尋找形式為 R(x) = P(x)/Q(x) 的有理函數，其中 P 與 Q 為同次多項式，且滿足 R(xData[i]) = yData[i]。演算法採用遞推方式計算插值函數在給定 x 處的值，避免顯式構造多項式係數，數值穩定性較好。

Python 實作

以下函數 rational(xData, yData, x) 回傳插值結果。程式碼使用 NumPy 陣列操作，核心邏輯為雙層迴圈更新中間值 r 與 rOld。

from numpy import zeros

def rational(xData, yData, x):
    """
    計算對角線有理函數插值 p(x)
    參數:
        xData: 節點 x 座標列表或陣列
        yData: 節點 y 座標列表或陣列
        x:     欲插值的點
    回傳:
        插值結果 p(x)
    """
    m = len(xData)
    r = yData.copy()
    rOld = zeros(m)

    for k in range(m - 1):
        for i in range(m - k - 1):
            # 若 x 極接近節點，直接回傳對應函數值避免除零
            if abs(x - xData[i + k + 1]) < 1.0e-9:
                return yData[i + k + 1]
            else:
                c1 = r[i + 1] - r[i]
                c2 = r[i + 1] - rOld[i + 1]
                c3 = (x - xData[i]) / (x - xData[i + k + 1])
                r[i] = r[i + 1] + c1 / (c3 * (1.0 - c1 / c2) - 1.0)
                rOld[i + 1] = r[i + 1]
    return r[0]

使用範例

假設我們有以下數據點：

xData = [0.0, 0.5, 1.0, 1.5, 2.0]
yData = [1.0, 0.8, 0.5, 0.3, 0.2]

欲在 x = 0.75 處插值：

result = rational(xData, yData, 0.75)
print(f"插值結果: {result}")

注意事項

• 輸入的 xData 與 yData 長度必須相同，且 xData 應單調遞增或遞減。
• 演算法在節點極接近時會直接回傳節點值，避免數值不穩定。
• 若插值點落在節點區間外，結果為外推，精度可能下降。

此實作簡潔高效，適合嵌入科學計算流程中。若需處理大規模數據，可考慮向量化優化或使用 SciPy 的 interpolate 模組。